Commit 790b7d98 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#1399 Add a post-processing test for mesh which vertices aren't all used to...

#1399 Add a post-processing test for mesh which vertices aren't all used to bear dofs. I had to broaden slightly the interface of Ensight 6 to enable choosing the output directory.
parent 24975a82
......@@ -3294,8 +3294,6 @@
BE1D44A822560BF000911391 /* Mpi.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = Mpi.cpp; sourceTree = "<group>"; };
BE1D44AE22560CD100911391 /* Environment.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = Environment.cpp; sourceTree = "<group>"; };
BE1D44B022560CE100911391 /* Environment.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = Environment.hpp; sourceTree = "<group>"; };
BE1D44BD225B462800911391 /* cylinderBug.lua */ = {isa = PBXFileReference; lastKnownFileType = text; path = cylinderBug.lua; sourceTree = "<group>"; };
BE1D44BE225B462800911391 /* cylinderOK.lua */ = {isa = PBXFileReference; lastKnownFileType = text; path = cylinderOK.lua; sourceTree = "<group>"; };
BE1D44BF225B906600911391 /* CMakeLists.txt */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = text; path = CMakeLists.txt; sourceTree = "<group>"; };
BE1D44C0225B907000911391 /* CMakeLists.txt */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = text; path = CMakeLists.txt; sourceTree = "<group>"; };
BE1E873C1B8DBF460002EE64 /* Definitions.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = Definitions.hpp; sourceTree = "<group>"; };
......@@ -3581,6 +3579,9 @@
BE490784225D27D3000C212B /* main.cpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.cpp; path = main.cpp; sourceTree = "<group>"; };
BE490785225D27D3000C212B /* VariationalFormulationInit.cpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.cpp; path = VariationalFormulationInit.cpp; sourceTree = "<group>"; };
BE490786225D27D3000C212B /* ExpectedResults.hpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.h; path = ExpectedResults.hpp; sourceTree = "<group>"; };
BE49079B2260871B000C212B /* README.md */ = {isa = PBXFileReference; lastKnownFileType = net.daringfireball.markdown; path = README.md; sourceTree = "<group>"; };
BE49079C2260871B000C212B /* CMakeLists.txt */ = {isa = PBXFileReference; lastKnownFileType = text; path = CMakeLists.txt; sourceTree = "<group>"; };
BE49079D2260874F000C212B /* test.cpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.cpp; path = test.cpp; sourceTree = "<group>"; };
BE4940FF224250E800157863 /* FiberDensityJ1J4J6.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = FiberDensityJ1J4J6.hpp; sourceTree = "<group>"; };
BE494100224250E800157863 /* FiberDensityJ1J4J6.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = FiberDensityJ1J4J6.hxx; sourceTree = "<group>"; };
BE494101224250E800157863 /* FiberDensityJ1J4J6.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = FiberDensityJ1J4J6.cpp; sourceTree = "<group>"; };
......@@ -4410,7 +4411,6 @@
BEB0C6161EC304EF00D62905 /* README */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = text; path = README; sourceTree = "<group>"; };
BEB0C6171EC3054500D62905 /* README */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = text; path = README; sourceTree = "<group>"; };
BEB0C6181EC3054E00D62905 /* README */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = text; path = README; sourceTree = "<group>"; };
BEB1B555225C9D25001DCCC3 /* cylinderOK_ensight.lua */ = {isa = PBXFileReference; lastKnownFileType = text; path = cylinderOK_ensight.lua; sourceTree = "<group>"; };
BEB1B557225C9DD8001DCCC3 /* CMakeLists.txt */ = {isa = PBXFileReference; lastKnownFileType = text; path = CMakeLists.txt; sourceTree = "<group>"; };
BEB1B559225CB155001DCCC3 /* CMakeLists.txt */ = {isa = PBXFileReference; lastKnownFileType = text; path = CMakeLists.txt; sourceTree = "<group>"; };
BEB1B55C225CB156001DCCC3 /* README.md */ = {isa = PBXFileReference; lastKnownFileType = net.daringfireball.markdown; path = README.md; sourceTree = "<group>"; };
......@@ -6698,6 +6698,7 @@
BE5A4E1F1B677AE0006822DD /* Operators */,
BE6A983D1BC66CDC004184F2 /* Parameter */,
BEDEDE0A213002990060F00E /* ParameterInstances */,
BE4907992260851D000C212B /* PostProcessing */,
);
name = Test;
path = Sources/Test;
......@@ -7146,6 +7147,24 @@
path = Sources/Test/Operators/TestFunctions;
sourceTree = SOURCE_ROOT;
};
BE4907992260851D000C212B /* PostProcessing */ = {
isa = PBXGroup;
children = (
BE49079A22608578000C212B /* EnsightUnusedVertices */,
);
path = PostProcessing;
sourceTree = "<group>";
};
BE49079A22608578000C212B /* EnsightUnusedVertices */ = {
isa = PBXGroup;
children = (
BE49079B2260871B000C212B /* README.md */,
BE49079C2260871B000C212B /* CMakeLists.txt */,
BE49079D2260874F000C212B /* test.cpp */,
);
path = EnsightUnusedVertices;
sourceTree = "<group>";
};
BE4D0C5821A2E4F800E0D4E7 /* Instances */ = {
isa = PBXGroup;
children = (
......@@ -7998,9 +8017,6 @@
BE63B4911A31C20E003A6523 /* Elasticity */ = {
isa = PBXGroup;
children = (
BE1D44BD225B462800911391 /* cylinderBug.lua */,
BE1D44BE225B462800911391 /* cylinderOK.lua */,
BEB1B555225C9D25001DCCC3 /* cylinderOK_ensight.lua */,
BE964C7D206E8E3B00B9ED1E /* CMakeLists.txt */,
BEB0C6161EC304EF00D62905 /* README */,
BE6EA61E1EC32FFB0085E651 /* demo_input_elasticity_3d.lua */,
......@@ -42,7 +42,7 @@ target_link_libraries(MoReFEM_post_processing
target_link_libraries(${MOREFEM_TEST_TOOLS}
${ALL_LOAD_BEGIN_FLAG}
${MOREFEM_MODEL}
${MOREFEM_POST_PROCESSING}
${ALL_LOAD_END_FLAG}
${LIB_BOOST_TEST})
......
......@@ -50,7 +50,8 @@ namespace MoReFEM::PostProcessingNS::OutputFormat
const std::vector<std::string>& unknown_list,
const std::vector<unsigned int>& numbering_subset_id_list,
const Mesh& mesh,
RefinedMesh is_mesh_refined)
RefinedMesh is_mesh_refined,
const std::string& output_directory)
{
assert(unknown_list.size() == numbering_subset_id_list.size());
const auto Nunknown = unknown_list.size();
......@@ -66,7 +67,10 @@ namespace MoReFEM::PostProcessingNS::OutputFormat
throw Exception("Invalid mesh directory: '" + mesh_directory + "' doesn't exist!",
__FILE__, __LINE__);
const std::string ensight_directory = mesh_directory + "Ensight6/";
const std::string ensight_directory =
output_directory.empty()
? mesh_directory + "Ensight6/"
: output_directory;
if (FilesystemNS::Folder::DoExist(ensight_directory))
FilesystemNS::Folder::Remove(ensight_directory, __FILE__, __LINE__);
......@@ -77,7 +81,7 @@ namespace MoReFEM::PostProcessingNS::OutputFormat
numbering_subset_id_list,
mesh);
mesh.Write<MeshNS::Format::Ensight>(ensight_directory + "mesh.geo");
mesh.Write<MeshNS::Format::Ensight>(ensight_directory + "/mesh.geo");
const auto& time_iteration_list = post_processing.GetTimeIterationList();
......
......@@ -85,12 +85,17 @@ namespace MoReFEM
* \param[in] unknown_list Define here the list of unknowns for which Ensight output is sought.
* If empty, all the unknowns found in the output file are taken; if not only the selected subset is
* taken. The name specified here MUST match some of the unknowns in the unknown file.
* \param[in] output_directory Directory into which the Ensight outputs will be written. If left empty,
* a subdirectory named "Ensight6" will be created in the Mesh directory. Beware: this directory is
* removed and recreated; make sure not to provide a directory into which you have some content you
* don't want to lose.
*/
explicit Ensight6(const std::string& data_directory,
const std::vector<std::string>& unknown_list,
const std::vector<unsigned int>& numbering_subset_id_list,
const Mesh& mesh,
RefinedMesh is_mesh_refined = RefinedMesh::no);
RefinedMesh is_mesh_refined = RefinedMesh::no,
const std::string& output_directory = "");
//! Destructor.
~Ensight6() = default;
......
/*!
//
// \file
//
//
// Created by Sebastien Gilles <sebastien.gilles@inria.fr> on the Tue, 23 Dec 2014 16:13:29 +0100
// Copyright (c) Inria. All rights reserved.
//
// \ingroup PostProcessingGroup
// \addtogroup PostProcessingGroup
// \{
*/
#include <fstream>
#include "Utilities/Filesystem/File.hpp"
#include "Utilities/String/String.hpp"
#include "Geometry/Mesh/Mesh.hpp"
#include "PostProcessing/PostProcessing.hpp"
#include "PostProcessing/OutputFormat/Ensight6.hpp"
namespace MoReFEM::PostProcessingNS::OutputFormat
{
namespace // anonymous
{
void CreateCaseFile(const std::string& output_directory,
const Data::UnknownInformation::vector_const_shared_ptr& unknown_list,
const Data::TimeIteration::vector_const_unique_ptr& time_iteration_list);
void CreateUnknownFile(const Data::UnknownInformation& unknown,
const unsigned int numbering_subset_id,
const PostProcessingNS::Data::TimeIteration& time_iteration,
const std::string& data_directory,
const PostProcessing& post_processing,
RefinedMesh is_mesh_refined);
} // namespace anonymous
Ensight6::Ensight6(const std::string& data_directory,
const std::vector<std::string>& unknown_list,
const std::vector<unsigned int>& numbering_subset_id_list,
const Mesh& mesh,
RefinedMesh is_mesh_refined)
{
assert(unknown_list.size() == numbering_subset_id_list.size());
const auto Nunknown = unknown_list.size();
if (!FilesystemNS::Folder::DoExist(data_directory))
throw Exception("Invalid data directory: '" + data_directory + "' doesn't exist!",
__FILE__, __LINE__);
const std::string mesh_directory = data_directory + "/Mesh_"
+ std::to_string(mesh.GetUniqueId()) + "/";
if (!FilesystemNS::Folder::DoExist(mesh_directory))
throw Exception("Invalid mesh directory: '" + mesh_directory + "' doesn't exist!",
__FILE__, __LINE__);
const std::string ensight_directory = mesh_directory + "Ensight6/";
if (FilesystemNS::Folder::DoExist(ensight_directory))
FilesystemNS::Folder::Remove(ensight_directory, __FILE__, __LINE__);
FilesystemNS::Folder::Create(ensight_directory, __FILE__, __LINE__);
PostProcessing post_processing(data_directory,
numbering_subset_id_list,
mesh);
mesh.Write<MeshNS::Format::Ensight>(ensight_directory + "mesh.geo");
const auto& time_iteration_list = post_processing.GetTimeIterationList();
// Create the case file.
decltype(auto) complete_unknown_list = post_processing.GetExtendedUnknownList();
Data::UnknownInformation::vector_const_shared_ptr selected_unknown_list;
// Associate the more complete unknown object from the name given in input.
{
for (const auto& unknown_name : unknown_list)
{
const auto it = std::find_if(complete_unknown_list.cbegin(),
complete_unknown_list.cend(),
[&unknown_name](const auto& unknown_ptr)
{
assert(!(!unknown_ptr));
return unknown_ptr->GetName() == unknown_name;
});
if (it == complete_unknown_list.cend())
throw Exception("Unknown '" + unknown_name + "' required in the constructor was not one of "
"those used in the model.",
__FILE__, __LINE__);
selected_unknown_list.push_back(*it);
}
}
assert(selected_unknown_list.size() == Nunknown);
CreateCaseFile(ensight_directory,
selected_unknown_list,
time_iteration_list);
// Sort the unknowns by associated numbering subset.
std::multimap<unsigned int, Data::UnknownInformation::const_shared_ptr> unknown_list_per_numbering_subset;
for (auto i = 0ul; i < Nunknown; ++i)
unknown_list_per_numbering_subset.insert({numbering_subset_id_list[i], selected_unknown_list[i]});
// Iterate over time:
for (const auto& time_iteration_ptr : time_iteration_list)
{
assert(!(!time_iteration_ptr));
const auto& time_iteration = *time_iteration_ptr;
const auto range = unknown_list_per_numbering_subset.equal_range(time_iteration.GetNumberingSubsetId());
for (auto it = range.first; it != range.second; ++it)
{
const auto& unknown_list_for_numbering_subset = *it;
decltype(auto) unknown_ptr = unknown_list_for_numbering_subset.second;
assert(!(!unknown_ptr));
CreateUnknownFile(*unknown_ptr,
unknown_list_for_numbering_subset.first,
time_iteration,
ensight_directory,
post_processing,
is_mesh_refined);
}
}
}
namespace // anonymous
{
void CreateCaseFile(const std::string& output_directory,
const Data::UnknownInformation::vector_const_shared_ptr& unknown_list,
const Data::TimeIteration::vector_const_unique_ptr& time_iteration_list)
{
std::string case_file(output_directory);
case_file += "/problem.case";
std::ofstream stream;
FilesystemNS::File::Create(stream, case_file, __FILE__, __LINE__);
stream << "FORMAT" << std::endl;
stream << "type: ensight" << std::endl;
stream << "GEOMETRY" << std::endl;
stream << "model: 1 mesh.geo" << std::endl;
stream << "VARIABLE" << std::endl;
for (const auto& unknown_ptr : unknown_list)
{
assert(!(!unknown_ptr));
switch (unknown_ptr->GetNature())
{
case Data::UnknownNature::scalar:
stream << "scalar";
break;
case Data::UnknownNature::vectorial:
stream << "vector";
break;
}
const auto& name = unknown_ptr->GetName();
stream << " per node: 1 " << name << ' ' << name << ".*****.scl" << std::endl;
}
stream << "TIME" << std::endl;
stream << "time set: 1" << std::endl;
unsigned int count = 0u;
unsigned int previous_index = static_cast<unsigned int>(-1);
std::ostringstream time_values_stream;
std::ostringstream time_index_stream;
auto Ntime_iteration = 0u;
for (const auto& time_iteration_ptr : time_iteration_list)
{
assert(!(!time_iteration_ptr));
const auto& time_iteration = *time_iteration_ptr;
const auto current_iteration = time_iteration.GetIteration();
if (current_iteration == previous_index) // May happen if several unknowns for same time step.
continue;
++Ntime_iteration;
time_index_stream << current_iteration;
time_values_stream << time_iteration.GetTime();
if (++count % 5u == 0u) // lines must not exceed 79 characters; I'm truly on the very safe side here!
{
time_index_stream << std::endl;
time_values_stream << std::endl;
}
else
{
time_index_stream << ' ';
time_values_stream << ' ';
}
previous_index = current_iteration;
}
stream << "number of steps: " << Ntime_iteration << std::endl;
stream << "filename numbers: ";
stream << time_index_stream.str();
stream << std::endl;
stream << "time values: ";
stream << time_values_stream.str();
stream << std::endl;
}
void CreateUnknownFile(const Data::UnknownInformation& unknown,
const unsigned int numbering_subset_id,
const PostProcessingNS::Data::TimeIteration& time_iteration,
const std::string& output_directory,
const PostProcessing& post_processing,
RefinedMesh is_mesh_refined)
{
const auto& unknown_name = unknown.GetName();
const auto& mesh = post_processing.GetMesh();
std::ostringstream oconv;
oconv << output_directory;
oconv << '/' << unknown_name << '.' << std::setfill('0') << std::setw(5)
<< time_iteration.GetIteration() << ".scl";
std::string ensight_output_file(oconv.str());
std::ofstream ensight_output_stream;
FilesystemNS::File::Create(ensight_output_stream, ensight_output_file, __FILE__, __LINE__);
ensight_output_stream << "First line which content must be chosen (but later)" << std::endl;
const unsigned int dimension = mesh.GetDimension();
const auto unknown_nature = unknown.GetNature();
const std::size_t Nprocessor = post_processing.Nprocessor();
using DofWithValuesType = std::pair<Data::DofInformation::const_shared_ptr, double>;
std::vector<DofWithValuesType> dof_with_values;
#ifndef NDEBUG
if (is_mesh_refined == RefinedMesh::yes)
{
assert(Nprocessor == 1ul
&& "At the moment, restricted to model run sequentially.");
}
#endif // NDEBUG
decltype(auto) filename_prototype = time_iteration.GetSolutionFilename();
for (std::size_t processor = 0u; processor < Nprocessor; ++processor)
{
const auto& dof_list = post_processing.GetDofInformationList(numbering_subset_id, processor);
std::string filename = filename_prototype;
const auto check = Utilities::String::Replace("*", std::to_string(processor), filename);
assert(check == 1);
static_cast<void>(check);
const auto& solution = LoadVector(filename);
for (const auto& dof_ptr : dof_list)
{
assert(!(!dof_ptr));
const auto& dof = *dof_ptr;
if (dof.GetUnknown() != unknown_name)
continue;
if (is_mesh_refined == RefinedMesh::no)
{
// Very crude: dofs not on vertex are simply ignored.
if (dof.GetInterfaceNature() != InterfaceNS::Nature::vertex)
{
static bool first = true;
if (first)
{
std::cout << "[WARNING] As ensight handles only P1, dofs not on vertices are "
"blatantly ignored." << std::endl;
first = false;
}
continue;
}
}
const unsigned int processor_wise_index = dof.GetProcessorWiseIndex();
assert(processor_wise_index < solution.size());
dof_with_values.push_back({dof_ptr, solution[processor_wise_index]});
}
}
std::cout << "Ndof = " << dof_with_values.size() << std::endl;
// The sorting below is required to handle the parallel case; it assumes P1 or Q1.
// That's the reason refined mesh does not work with model run in parallel.
if (is_mesh_refined == RefinedMesh::no)
{
// Reminder:
// - Post-processing runs in sequential;
// - The mesh.geo file was built using that very same order.
decltype(auto) coords_list = mesh.GetProcessorWiseCoordsList();
std::vector<std::vector<std::size_t> > match(coords_list.size());
auto coords_list_counter = 0u;
--coords_list_counter;
for (const auto& coords_ptr : coords_list)
{
++coords_list_counter;
assert(!(!coords_ptr));
const auto& coords = *coords_ptr;
auto dof_with_values_counter = 0ul;
--dof_with_values_counter;
// Find - if any - the dofs that match these positions.
for (const auto& pair : dof_with_values)
{
++dof_with_values_counter;
decltype(auto) dof_ptr = pair.first;
assert(!(!dof_ptr));
decltype(auto) dof_coords_index_list = dof_ptr->GetInterface().GetVertexCoordsIndexList();
assert(dof_coords_index_list.size() == 1ul
&& "At the moment EnsightOutput deals only "
"with P1 (from both geometric and felt point of view)!");
const unsigned int dof_coords_index = dof_coords_index_list.back();
decltype(auto) dof_coords = mesh.GetCoord(dof_coords_index);
if (NumericNS::AreEqual(coords.x(), dof_coords.x())
&& NumericNS::AreEqual(coords.y(), dof_coords.y())
&& NumericNS::AreEqual(coords.z(), dof_coords.z()))
{
assert(coords_list_counter < match.size());
match[coords_list_counter].push_back(dof_with_values_counter);
}
}
}
coords_list_counter = 0u;
--coords_list_counter;
for (const auto& item : match)
{
++coords_list_counter;
std::cout << "For coords " << coords_list_counter << ": ";
Utilities::PrintContainer(item);
}
std::sort(dof_with_values.begin(), dof_with_values.end(),
[](const DofWithValuesType& lhs, const DofWithValuesType& rhs)
{
if (lhs == rhs)
return false;
const auto& lhs_dof = lhs.first;
const auto& rhs_dof = rhs.first;
assert(!(!lhs_dof));
assert(!(!rhs_dof));
decltype(auto) lhs_vertex_coords_list = lhs_dof->GetInterface().GetVertexCoordsIndexList();
decltype(auto) rhs_vertex_coords_list = rhs_dof->GetInterface().GetVertexCoordsIndexList();
assert(lhs_vertex_coords_list.size() == 1ul && "At the moment EnsightOutput deals only "
"with P1 (from both geometric and felt point of view)!");
assert(rhs_vertex_coords_list.size() == 1ul && "At the moment EnsightOutput deals only "
"with P1 (from both geometric and felt point of view)!");
const unsigned int lhs_index = lhs_vertex_coords_list.back();
const unsigned int rhs_index = rhs_vertex_coords_list.back();
assert(lhs_dof->GetUnknown() == rhs_dof->GetUnknown()
&& "Filter upon unknowns should have occurred earlier");
if (lhs_index != rhs_index)
return lhs_index < rhs_index;
assert(lhs_dof->GetUnknownComponent() != rhs_dof->GetUnknownComponent());
return lhs_dof->GetUnknownComponent() < rhs_dof->GetUnknownComponent();
});
}
unsigned int counter = 0u;
for (const auto& pair : dof_with_values)
{
ensight_output_stream << std::setw(12) << std::scientific << std::setprecision(5)
<< pair.second;
assert(!(!pair.first));
const auto& dof = *(pair.first);
if (unknown_nature == Data::UnknownNature::vectorial
&& dimension == 2u && dof.GetUnknownComponent() == 1u)
{
ensight_output_stream << std::setw(12) << std::scientific << std::setprecision(5) << 0.;
if (++counter % 6u == 0u) // Ensight expects 6 values per line
ensight_output_stream << '\n';
}
if (++counter % 6u == 0u) // Ensight expects 6 values per line
ensight_output_stream << '\n';
}
if (counter % 6u != 0u)
ensight_output_stream << '\n';
}