Commit 24975a82 authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#1399 Rewrite the post-processing method that writes the dof output for each...

#1399 Rewrite the post-processing method that writes the dof output for each unknown; the bug leading to ticket #1399 is therefore resolved. I will however write a test and improve it a bit: current implementation was more a proof of concept and is wildly inefficient (some computation done for each time iteration - not an issue on my toy meshes but could be more annoying for heart models...).
parent 0cdedaf3
......@@ -72,7 +72,7 @@
</CommandLineArgument>
<CommandLineArgument
argument = "-i ${HOME}/Codes/MoReFEM/CoreLibrary/Sources/ModelInstances/Elasticity/cylinderBug.lua"
isEnabled = "NO">
isEnabled = "YES">
</CommandLineArgument>
<CommandLineArgument
argument = "-i ${HOME}/Codes/MoReFEM/CoreLibrary/Sources/ModelInstances/Elasticity/cylinderOK.lua"
......
......@@ -68,11 +68,11 @@
</CommandLineArgument>
<CommandLineArgument
argument = "-i /Volumes/Data/${USER}/MoReFEM/Results/CylinderBug/input_data.lua"
isEnabled = "NO">
isEnabled = "YES">
</CommandLineArgument>
<CommandLineArgument
argument = "-i /Volumes/Data/${USER}/MoReFEM/Results/CylinderOK/input_data.lua"
isEnabled = "YES">
isEnabled = "NO">
</CommandLineArgument>
</CommandLineArguments>
<AdditionalOptions>
......
......@@ -335,74 +335,95 @@ namespace MoReFEM::PostProcessingNS::OutputFormat
}
decltype(auto) coords_list = mesh.GetProcessorWiseCoordsList();
std::vector<std::vector<std::size_t> > match(coords_list.size());
// 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.
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;
auto coords_list_counter = 0u;
--coords_list_counter;
assert(!(!lhs_dof));
assert(!(!rhs_dof));
for (const auto& coords_ptr : coords_list)
{
++coords_list_counter;
assert(!(!coords_ptr));
const auto& coords = *coords_ptr;
decltype(auto) lhs_vertex_coords_list = lhs_dof->GetInterface().GetVertexCoordsIndexList();
decltype(auto) rhs_vertex_coords_list = rhs_dof->GetInterface().GetVertexCoordsIndexList();
auto dof_with_values_counter = 0ul;
--dof_with_values_counter;
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)!");
// Find - if any - the dofs that match these positions.
for (const auto& pair : dof_with_values)
{
++dof_with_values_counter;
const unsigned int lhs_index = lhs_vertex_coords_list.back();
const unsigned int rhs_index = rhs_vertex_coords_list.back();
decltype(auto) dof_ptr = pair.first;
assert(!(!dof_ptr));
decltype(auto) dof_coords_index_list = dof_ptr->GetInterface().GetVertexCoordsIndexList();
assert(lhs_dof->GetUnknown() == rhs_dof->GetUnknown()
&& "Filter upon unknowns should have occurred earlier");
assert(dof_coords_index_list.size() == 1ul
&& "At the moment EnsightOutput deals only "
"with P1 (from both geometric and felt point of view)!");
if (lhs_index != rhs_index)
return lhs_index < rhs_index;
const unsigned int dof_coords_index = dof_coords_index_list.back();
assert(lhs_dof->GetUnknownComponent() != rhs_dof->GetUnknownComponent());
decltype(auto) dof_coords = mesh.GetCoord(dof_coords_index);
return lhs_dof->GetUnknownComponent() < rhs_dof->GetUnknownComponent();
});
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);
}
}
}
}
unsigned int counter = 0u;
unsigned int dof_counter = 0u;
for (const auto& pair : dof_with_values)
//for (const auto& pair : dof_with_values)
for (const auto& dof_list_for_coords : match)
{
ensight_output_stream << std::setw(12) << std::scientific << std::setprecision(5)
<< pair.second;
assert(!(!pair.first));
const auto& dof = *(pair.first);
if (dof_list_for_coords.empty())
{
for (auto i = 0; i < 3; ++i)
ensight_output_stream << std::setw(12) << std::scientific << std::setprecision(5) << 0.;
if (unknown_nature == Data::UnknownNature::vectorial
&& dimension == 2u && dof.GetUnknownComponent() == 1u)
dof_counter += 3u;
}
else
{
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';
assert(dof_list_for_coords.size() == 1ul
|| dof_list_for_coords.size() == dimension);
for (const auto& dof_index : dof_list_for_coords)
{
assert(dof_index < dof_with_values.size());
const auto dof_value = dof_with_values[dof_index].second;
ensight_output_stream << std::setw(12) << std::scientific << std::setprecision(5)
<< dof_value;
}
if (unknown_nature == Data::UnknownNature::vectorial && dimension == 2u)
ensight_output_stream << std::setw(12) << std::scientific << std::setprecision(5) << 0.;
dof_counter += (unknown_nature == Data::UnknownNature::vectorial ? 3u : 1u);
}
if (++counter % 6u == 0u) // Ensight expects 6 values per line
if (dof_counter % 6u == 0u) // Ensight expects 6 values per line
ensight_output_stream << '\n';
}
if (counter % 6u != 0u)
ensight_output_stream << '\n';
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment