Commit 15787a3b authored by GILLES Sebastien's avatar GILLES Sebastien
Browse files

#1520 Refactoring the GlobalCoordsQuadPt test (still in progress)

parent 0365ad3d
......@@ -4372,6 +4372,22 @@
BEA290AC240E8D9C00CC9594 /* used-but-marked-unused.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = "used-but-marked-unused.hpp"; sourceTree = "<group>"; };
BEA3248617A7E5BA00ADEB73 /* Exception.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = Exception.cpp; sourceTree = "<group>"; };
BEA3248717A7E5BA00ADEB73 /* Exception.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = Exception.hpp; sourceTree = "<group>"; };
BEA348F624A4BD4C00BF1CDE /* CMakeLists.txt */ = {isa = PBXFileReference; lastKnownFileType = text; path = CMakeLists.txt; sourceTree = "<group>"; };
BEA348F824A4BD4C00BF1CDE /* CMakeLists.txt */ = {isa = PBXFileReference; lastKnownFileType = text; path = CMakeLists.txt; sourceTree = "<group>"; };
BEA348F924A4BD4C00BF1CDE /* main_2D_quad.cpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.cpp; path = main_2D_quad.cpp; sourceTree = "<group>"; };
BEA348FA24A4BD4C00BF1CDE /* demo_3D_tetra.lua */ = {isa = PBXFileReference; lastKnownFileType = text; path = demo_3D_tetra.lua; sourceTree = "<group>"; };
BEA348FB24A4BD4C00BF1CDE /* main_3D_hexa.cpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.cpp; path = main_3D_hexa.cpp; sourceTree = "<group>"; };
BEA348FC24A4BD4C00BF1CDE /* Model.hpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.h; path = Model.hpp; sourceTree = "<group>"; };
BEA348FD24A4BD4C00BF1CDE /* InputData.hpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.h; path = InputData.hpp; sourceTree = "<group>"; };
BEA348FE24A4BD4C00BF1CDE /* Model.hxx */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.h; path = Model.hxx; sourceTree = "<group>"; };
BEA348FF24A4BD4C00BF1CDE /* main_2D_triangle.cpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.cpp; path = main_2D_triangle.cpp; sourceTree = "<group>"; };
BEA3490024A4BD4C00BF1CDE /* main_3D_tetra.cpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.cpp; path = main_3D_tetra.cpp; sourceTree = "<group>"; };
BEA3490124A4BD4C00BF1CDE /* Model.cpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.cpp; path = Model.cpp; sourceTree = "<group>"; };
BEA3490224A4BD4C00BF1CDE /* demo_2D_triangle.lua */ = {isa = PBXFileReference; lastKnownFileType = text; path = demo_2D_triangle.lua; sourceTree = "<group>"; };
BEA3490324A4BD4C00BF1CDE /* demo_2D_quad.lua */ = {isa = PBXFileReference; lastKnownFileType = text; path = demo_2D_quad.lua; sourceTree = "<group>"; };
BEA3490424A4BD4C00BF1CDE /* demo_3D_hexa.lua */ = {isa = PBXFileReference; lastKnownFileType = text; path = demo_3D_hexa.lua; sourceTree = "<group>"; };
BEA3490524A4BD4C00BF1CDE /* main_1D.cpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.cpp; path = main_1D.cpp; sourceTree = "<group>"; };
BEA3490624A4BD4C00BF1CDE /* demo_1D_edge.lua */ = {isa = PBXFileReference; lastKnownFileType = text; path = demo_1D_edge.lua; sourceTree = "<group>"; };
BEA355D617D0971500FB643B /* LuaFunction.hpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = LuaFunction.hpp; sourceTree = "<group>"; };
BEA355D717D0971500FB643B /* LuaFunction.hxx */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.h; path = LuaFunction.hxx; sourceTree = "<group>"; };
BEA4AF0F246059D200669BBB /* WritePrepartitionedData.cpp */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.cpp.cpp; path = WritePrepartitionedData.cpp; sourceTree = "<group>"; };
......@@ -7933,6 +7949,7 @@
BE5A4E201B677B23006822DD /* ConformProjector */,
BEC77E951DDDC34800F8B444 /* NonConformInterpolator */,
BE7C942A1F5FE614003D2C52 /* TestFunctions */,
BEA348F524A4BD4C00BF1CDE /* ParameterOperator */,
BEBFCE9D205919640033C4C1 /* VariationalInstances */,
);
path = Operators;
......@@ -9551,6 +9568,37 @@
path = Sources/ModelInstances;
sourceTree = "<group>";
};
BEA348F524A4BD4C00BF1CDE /* ParameterOperator */ = {
isa = PBXGroup;
children = (
BEA348F624A4BD4C00BF1CDE /* CMakeLists.txt */,
BEA348F724A4BD4C00BF1CDE /* GlobalCoordsQuadPts */,
);
path = ParameterOperator;
sourceTree = "<group>";
};
BEA348F724A4BD4C00BF1CDE /* GlobalCoordsQuadPts */ = {
isa = PBXGroup;
children = (
BEA348F824A4BD4C00BF1CDE /* CMakeLists.txt */,
BEA348FD24A4BD4C00BF1CDE /* InputData.hpp */,
BEA3490124A4BD4C00BF1CDE /* Model.cpp */,
BEA348FC24A4BD4C00BF1CDE /* Model.hpp */,
BEA348FE24A4BD4C00BF1CDE /* Model.hxx */,
BEA3490524A4BD4C00BF1CDE /* main_1D.cpp */,
BEA348F924A4BD4C00BF1CDE /* main_2D_quad.cpp */,
BEA348FF24A4BD4C00BF1CDE /* main_2D_triangle.cpp */,
BEA348FB24A4BD4C00BF1CDE /* main_3D_hexa.cpp */,
BEA3490024A4BD4C00BF1CDE /* main_3D_tetra.cpp */,
BEA348FA24A4BD4C00BF1CDE /* demo_3D_tetra.lua */,
BEA3490224A4BD4C00BF1CDE /* demo_2D_triangle.lua */,
BEA3490324A4BD4C00BF1CDE /* demo_2D_quad.lua */,
BEA3490424A4BD4C00BF1CDE /* demo_3D_hexa.lua */,
BEA3490624A4BD4C00BF1CDE /* demo_1D_edge.lua */,
);
path = GlobalCoordsQuadPts;
sourceTree = "<group>";
};
BEA4FC4118214D62002B2EA1 /* FormulationSolver */ = {
isa = PBXGroup;
children = (
......@@ -19,10 +19,15 @@ namespace MoReFEM::TestNS::GlobalCoordsQuadPt
namespace // anonymous
{
void CheckQuadratureRule(const QuadratureRulePerTopology& quadrature_rule_per_topology,
const ParameterAtQuadraturePoint<ParameterNS::Type::vector>& global_coords,
const std::string&& quadrature_order);
}
std::string_view quadrature_order);
} // namespace anonymous
Model::Model(const morefem_data_type& morefem_data)
: parent(morefem_data, create_domain_list_for_coords::yes, print_banner::no)
......@@ -40,110 +45,72 @@ namespace MoReFEM::TestNS::GlobalCoordsQuadPt
void Model::SupplInitialize()
{
quadrature_rule_per_topology_low_degree_ = std::make_unique<const QuadratureRulePerTopology>(1, 1);
// Required to enable construction of an operator after initialization step.
parent::SetClearGodOfDofTemporaryDataToFalse();
}
quadrature_rule_per_topology_medium_degree_ = std::make_unique<const QuadratureRulePerTopology>(3, 3);
quadrature_rule_per_topology_high_degree_ = std::make_unique<const QuadratureRulePerTopology>(5, 5);
void Model::CheckQuadrature(const QuadratureRulePerTopology* const quadrature_rule,
std::string_view quadrature_order) const
{
decltype(auto) god_of_dof = parent::GetGodOfDof(EnumUnderlyingType(MeshIndex::mesh));
decltype(auto) felt_space = god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::felt_space));
decltype(auto) time_manager = GetTimeManager();
namespace GPO = GlobalParameterOperatorNS;
decltype(auto) generic_unknown = UnknownManager::GetInstance(__FILE__, __LINE__)
.GetUnknown(EnumUnderlyingType(UnknownIndex::generic_unknown));
// Required to enable construction of an operator after initialization step.
parent::SetClearGodOfDofTemporaryDataToFalse();
decltype(auto) domain = DomainManager::GetInstance(__FILE__, __LINE__)
.GetDomain(EnumUnderlyingType(DomainIndex::domain), __FILE__, __LINE__);
const auto& god_of_dof = parent::GetGodOfDof(EnumUnderlyingType(MeshIndex::mesh));
const auto& felt_space = god_of_dof.GetFEltSpace(EnumUnderlyingType(FEltSpaceIndex::felt_space));
LocalVector init_vector;
const std::size_t spatial_point_size = 3;
constexpr const std::size_t spatial_point_size = 3;
init_vector.resize({spatial_point_size});
init_vector.fill(0.);
const auto& time_manager = GetTimeManager();
decltype(auto) domain = DomainManager::GetInstance(__FILE__, __LINE__)
.GetDomain(EnumUnderlyingType(DomainIndex::domain), __FILE__, __LINE__);
ParameterAtQuadraturePoint<ParameterNS::Type::vector>
global_coords_quad_pt("Global coords of quad points - Low degree",
domain,
*quadrature_rule,
init_vector,
time_manager);
namespace GPO = GlobalParameterOperatorNS;
const auto& generic_unknown_ptr = UnknownManager::GetInstance(__FILE__, __LINE__)
.GetUnknownPtr(EnumUnderlyingType(UnknownIndex::generic_unknown));
{
global_coords_quad_pt_low_degree_ =
std::make_unique<ParameterAtQuadraturePoint<ParameterNS::Type::vector>>("Global coords of quad points - Low degree",
domain,
*quadrature_rule_per_topology_low_degree_,
init_vector,
time_manager);
write_at_quad_pt_operator_low_degree_ =
std::make_unique<GPO::GlobalCoordsQuadPoints>(felt_space,
*generic_unknown_ptr,
quadrature_rule_per_topology_low_degree_.get(),
*global_coords_quad_pt_low_degree_);
write_at_quad_pt_operator_low_degree_->Update();
}
{
global_coords_quad_pt_medium_degree_ =
std::make_unique<ParameterAtQuadraturePoint<ParameterNS::Type::vector>>("Global coords of quad points - Medium degree",
domain,
*quadrature_rule_per_topology_medium_degree_,
init_vector,
time_manager);
write_at_quad_pt_operator_medium_degree_ =
std::make_unique<GPO::GlobalCoordsQuadPoints>(felt_space,
*generic_unknown_ptr,
quadrature_rule_per_topology_medium_degree_.get(),
*global_coords_quad_pt_medium_degree_);
write_at_quad_pt_operator_medium_degree_->Update();
}
GPO::GlobalCoordsQuadPoints write_at_quad_pt_operator(felt_space,
generic_unknown,
quadrature_rule,
global_coords_quad_pt);
{
global_coords_quad_pt_high_degree_ =
std::make_unique<ParameterAtQuadraturePoint<ParameterNS::Type::vector>>("Global coords of quad points - High degree",
domain,
*quadrature_rule_per_topology_high_degree_,
init_vector,
time_manager);
write_at_quad_pt_operator_high_degree_ =
std::make_unique<GPO::GlobalCoordsQuadPoints>(felt_space,
*generic_unknown_ptr,
quadrature_rule_per_topology_high_degree_.get(),
*global_coords_quad_pt_high_degree_);
write_at_quad_pt_operator_high_degree_->Update();
}
write_at_quad_pt_operator.Update();
CheckQuadratureRule(*quadrature_rule, global_coords_quad_pt, quadrature_order);
}
void Model::CheckLowOrderQuadrature() const
{
assert(!(!quadrature_rule_per_topology_low_degree_));
assert(!(!global_coords_quad_pt_low_degree_));
CheckQuadratureRule(*quadrature_rule_per_topology_low_degree_, *global_coords_quad_pt_low_degree_, "low");
auto rule = std::make_unique<const QuadratureRulePerTopology>(1, 1);
CheckQuadrature(rule.get(), "low");
}
void Model::CheckMediumOrderQuadrature() const
{
assert(!(!quadrature_rule_per_topology_medium_degree_));
assert(!(!global_coords_quad_pt_medium_degree_));
CheckQuadratureRule(*quadrature_rule_per_topology_medium_degree_, *global_coords_quad_pt_medium_degree_, "medium");
auto rule = std::make_unique<const QuadratureRulePerTopology>(3, 3);
CheckQuadrature(rule.get(), "medium");
}
void Model::CheckHighOrderQuadrature() const
{
assert(!(!quadrature_rule_per_topology_high_degree_));
assert(!(!global_coords_quad_pt_high_degree_));
CheckQuadratureRule(*quadrature_rule_per_topology_high_degree_, *global_coords_quad_pt_high_degree_, "high");
auto rule = std::make_unique<const QuadratureRulePerTopology>(5, 5);
CheckQuadrature(rule.get(), "high");
}
void Model::Forward()
......@@ -160,16 +127,18 @@ namespace MoReFEM::TestNS::GlobalCoordsQuadPt
namespace // anonymous
{
void CheckQuadratureRule(const QuadratureRulePerTopology& quadrature_rule_per_topology,
const ParameterAtQuadraturePoint<ParameterNS::Type::vector>& obtained_global_coords,
const std::string&& quadrature_order)
std::string_view quadrature_order)
{
decltype(auto) domain = DomainManager::GetInstance(__FILE__, __LINE__)
.GetDomain(EnumUnderlyingType(DomainIndex::domain), __FILE__, __LINE__);
const auto& mesh = domain.GetMesh();
const auto& ref_geom_elt_list = mesh.BagOfEltType();
decltype(auto) mesh = domain.GetMesh();
decltype(auto) ref_geom_elt_list = mesh.BagOfEltType();
std::vector<double> expected_global_coords;
for (const auto& ref_geom_elt_ptr : ref_geom_elt_list)
......@@ -180,6 +149,8 @@ namespace MoReFEM::TestNS::GlobalCoordsQuadPt
if (!domain.DoRefGeomEltMatchCriteria(ref_geom_elt))
continue;
const auto topology_identifier = ref_geom_elt.GetTopologyIdentifier();
const auto& quadrature_rule = quadrature_rule_per_topology.GetRule(ref_geom_elt.GetTopologyIdentifier());
const auto& ptr_pt_list = quadrature_rule.GetQuadraturePointList();
......@@ -190,8 +161,8 @@ namespace MoReFEM::TestNS::GlobalCoordsQuadPt
{
double expected_value;
if (ref_geom_elt.GetTopologyIdentifier() == RefGeomEltNS::TopologyNS::Type::tetrahedron
|| ref_geom_elt.GetTopologyIdentifier() == RefGeomEltNS::TopologyNS::Type::triangle)
if (topology_identifier == RefGeomEltNS::TopologyNS::Type::tetrahedron
|| topology_identifier == RefGeomEltNS::TopologyNS::Type::triangle)
{
// Integration interval is already defined on [0., 1.]^n for these topologies.
expected_value = local_coords;
......@@ -199,7 +170,7 @@ namespace MoReFEM::TestNS::GlobalCoordsQuadPt
else
{
// Integration interval is on [-1., 1.], this maps the values to the [0., 1.] interval.
expected_value = (local_coords + 1.) / 2.;
expected_value = (local_coords + 1.) * .5;
}
expected_global_coords.push_back(expected_value);
}
......@@ -213,13 +184,15 @@ namespace MoReFEM::TestNS::GlobalCoordsQuadPt
decltype(auto) environment = Utilities::Environment::GetInstance(__FILE__, __LINE__);
std::string output_dir = environment.SubstituteValues("${MOREFEM_TEST_OUTPUT_DIR}/Operators/GlobalCoordsQuadPts/");
output_dir += ref_geom_elt.GetTopologyName();
if (!Advanced::FilesystemNS::DirectoryNS::DoExist(output_dir))
Advanced::FilesystemNS::DirectoryNS::Create(output_dir, __FILE__, __LINE__);
std::string output_file(output_dir);
output_file += "/" + quadrature_order + "_order_global_coords.txt";
output_file += "/" + std::string(quadrature_order) + "_order_global_coords.txt";
obtained_global_coords.Write(output_file);
// Read from output and compare results.
......@@ -231,9 +204,9 @@ namespace MoReFEM::TestNS::GlobalCoordsQuadPt
for (int i = 0; i < 3; ++i)
getline(file_stream, line);
const auto setw_ensight = 12ul;
const auto vector_size = 3ul;
const auto ensight_precision = 1.e-6;
constexpr const auto setw_ensight = 12ul;
constexpr const auto vector_size = 3ul;
constexpr const auto ensight_precision = 1.e-6;
bool skipped_Nvalue_line = false;
std::string str_value;
......
......@@ -88,10 +88,10 @@ namespace MoReFEM::TestNS::GlobalCoordsQuadPt
/*!
* \brief Initialise the problem.
*
* This initialisation includes the resolution of the static problem.
*/
* \brief Initialise the problem.
*
* This initialisation includes the resolution of the static problem.
*/
void SupplInitialize();
......@@ -99,90 +99,66 @@ namespace MoReFEM::TestNS::GlobalCoordsQuadPt
void Forward();
/*!
* \brief Additional operations to finalize a dynamic step.
*
* Base class already update the time for next time iterations.
*/
* \brief Additional operations to finalize a dynamic step.
*
* Base class already update the time for next time iterations.
*/
void SupplFinalizeStep();
/*!
* \brief Initialise a dynamic step.
*
*/
* \brief Initialise a dynamic step.
*
*/
void SupplFinalize();
public:
/*!
* \brief Checks that the global coordinates for low order quadrature are the expected ones.
*
*/
* \brief Checks that the global coordinates for low order quadrature are the expected ones.
*
*/
void CheckLowOrderQuadrature() const;
/*!
* \brief Checks that the global coordinates for "medium" order quadrature are the expected ones.
*
*/
* \brief Checks that the global coordinates for "medium" order quadrature are the expected ones.
*
*/
void CheckMediumOrderQuadrature() const;
/*!
* \brief Checks that the global coordinates for high order quadrature are the expected ones.
*
*/
* \brief Checks that the global coordinates for high order quadrature are the expected ones.
*
*/
void CheckHighOrderQuadrature() const;
private:
/*!
* \brief Checks that the global coordinates for the quadrature given as argument are the expected ones.
*
* \param[in] quadrature_rule_per_topology \a QuadratureRule fow which the check is done for each topology.
*/
void CheckQuadrature(const QuadratureRulePerTopology* const quadrature_rule_per_topology,
std::string_view quadrature_order) const;
//! \copydoc doxygen_hide_model_SupplHasFinishedConditions_always_true
bool SupplHasFinishedConditions() const;
/*!
* \brief Part of InitializedStep() specific to Elastic model.
*
* As there are none, the body of this method is empty.
*/
* \brief Part of InitializedStep() specific to Elastic model.
*
* As there are none, the body of this method is empty.
*/
void SupplInitializeStep();
///@}
private:
//! Quadrature rule topology used for low degree integration.
QuadratureRulePerTopology::const_unique_ptr quadrature_rule_per_topology_low_degree_ = nullptr;
//! Quadrature rule topology used for "medium" degree integration.
QuadratureRulePerTopology::const_unique_ptr quadrature_rule_per_topology_medium_degree_ = nullptr;
//! Quadrature rule topology used for high degree integration.
QuadratureRulePerTopology::const_unique_ptr quadrature_rule_per_topology_high_degree_ = nullptr;
private:
//! Operator used to compute the global coordinates of the quadrature points.
GlobalParameterOperatorNS::GlobalCoordsQuadPoints::const_unique_ptr write_at_quad_pt_operator_low_degree_ = nullptr;
//! Parameter which will be updated by the operator to hold the global coordinates of the quadrature points.
ParameterAtQuadraturePoint<ParameterNS::Type::vector>::unique_ptr global_coords_quad_pt_low_degree_ = nullptr;
//! Operator used to compute the global coordinates of the quadrature points.
GlobalParameterOperatorNS::GlobalCoordsQuadPoints::const_unique_ptr write_at_quad_pt_operator_medium_degree_ = nullptr;
//! Parameter which will be updated by the operator to hold the global coordinates of the quadrature points.
ParameterAtQuadraturePoint<ParameterNS::Type::vector>::unique_ptr global_coords_quad_pt_medium_degree_ = nullptr;
//! Operator used to compute the global coordinates of the quadrature points.
GlobalParameterOperatorNS::GlobalCoordsQuadPoints::const_unique_ptr write_at_quad_pt_operator_high_degree_ = nullptr;
//! Parameter which will be updated by the operator to hold the global coordinates of the quadrature points.
ParameterAtQuadraturePoint<ParameterNS::Type::vector>::unique_ptr global_coords_quad_pt_high_degree_ = nullptr;
};
......
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