From 806268ccb867a5e426c4324299cd4a849ba90cf2 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Sun, 20 Sep 2026 16:33:06 -0400 Subject: [PATCH] Progress on bringing the DomainMapper over --- src/include/serif/mapping/data.hpp | 127 ++++++++++++++ src/include/serif/mapping/domain_mapper.hpp | 161 ++++++++++++++++++ src/include/serif/utils/error/codes.hpp | 6 + .../user_options/mapping/domain_mapper.hpp | 8 + src/lib/serif/mapping/data.cpp | 153 +++++++++++++++++ src/meson.build | 1 + .../stroid/build-check/meson.build | 7 + .../stroid/build-config/mfem/meson.build | 3 + tests/sandbox/mapping_sandbox.cpp | 1 + 9 files changed, 467 insertions(+) create mode 100644 src/include/serif/mapping/data.hpp create mode 100644 src/include/serif/mapping/domain_mapper.hpp create mode 100644 src/include/serif/utils/user_options/mapping/domain_mapper.hpp create mode 100644 src/lib/serif/mapping/data.cpp create mode 100644 subprojects/packagefiles/stroid/build-check/meson.build create mode 100644 subprojects/packagefiles/stroid/build-config/mfem/meson.build diff --git a/src/include/serif/mapping/data.hpp b/src/include/serif/mapping/data.hpp new file mode 100644 index 00000000..17622fe7 --- /dev/null +++ b/src/include/serif/mapping/data.hpp @@ -0,0 +1,127 @@ +#pragma once +#include "mfem/fem/fe/fe_base.hpp" +#include "mfem/fem/fespace.hpp" +#include "mfem/linalg/vector.hpp" +#include "mfem/linalg/densemat.hpp" +#include "mfem/mesh/mesh.hpp" + +#include +#include +#include + +#include "serif/mapping/compactification/kelvin.hpp" +#include "serif/mapping/context.hpp" +#include "serif/mapping/data.hpp" + + +namespace serif::mapping { + enum class FaceElementSide : uint8_t { + ELEMENT_1, + ELEMENT_2 + }; + + class ElementDisplacementData { + public: + ElementDisplacementData(const mfem::FiniteElement& element, mfem::Vector& displacement_dofs, mfem::Ordering::Type ordering = mfem::Ordering::byNODES); + + [[nodiscard]] const mfem::FiniteElement& get_element() const noexcept; + [[nodiscard]] const mfem::DenseMatrix &get_dof_matrix() const noexcept; + [[nodiscard]] int get_dimension() const noexcept; + [[nodiscard]] int get_dof_count() const noexcept; + [[nodiscard]] mfem::Ordering::Type get_ordering() const noexcept; + + private: + const mfem::FiniteElement *m_element; + mfem::DenseMatrix m_dof_matrix; + int m_dimension; + mfem::Ordering::Type m_ordering; + }; + + // Maybe this should be moved into the compactification module + struct CompactificationPointData { + double compactification_coordinate{0.0}; + mfem::Vector coordinate_gradient; + }; + + [[nodiscard]] ElementDisplacementData element_displacement_data_from_element_vdofs( + const mfem::FiniteElement& element, + mfem::Vector &displacement_dofs + ); + + class ElementCompactificationData { + public: + ElementCompactificationData(const mfem::FiniteElement& element, const mfem::Vector& dofs); + + [[nodiscard]] const mfem::FiniteElement& get_element() const noexcept; + [[nodiscard]] const mfem::Vector &get_dofs() const noexcept; + [[nodiscard]] int get_dof_count() const noexcept; + + private: + const mfem::FiniteElement *m_element; + mfem::Vector m_dofs; + }; + + struct ElementMappingData { + const ElementDisplacementData &displacement; + const ElementCompactificationData &compactification; + }; + + class DMWorkspace { + public: + explicit DMWorkspace (int dimension = 3); + void set_dimension(int dimension); + [[nodiscard]] int get_dimension() const noexcept; + public: // This is a workspace class, these members are kept public so we do not need indirection through getters to access them + mfem::Vector m_shape; + mfem::DenseMatrix m_reference_shape; + mfem::DenseMatrix m_mesh_dshape; + mfem::DenseMatrix m_reference_field_jacobian; + mfem::Vector m_field_value; + mfem::DenseMatrix m_field_jacobian; + + mfem::Vector m_compactification_shape; + mfem::DenseMatrix m_compactification_dshape; + CompactificationPointData m_compactification_point; + + mfem::Vector m_reference_normal; + mfem::Vector m_mapped_normal; + mfem::DenseMatrix m_full_element_jacobian; + + mfem::Vector m_vector_temp; + mfem::DenseMatrix m_matrix_temp_1; + mfem::DenseMatrix m_matrix_temp_2; + + compactification::ExteriorMapResult m_exterior_result; + compactification::ExteriorMapVariation m_exterior_variation; + private: + int m_dimension; + }; + + struct DMFixedMeshGeometry { + mfem::Vector position; + mfem::DenseMatrix jacobian; + mfem::DenseMatrix inverse_jacobian; + double weight{0.0}; + }; + + struct DMGeometryCacheStatistics { + std::uint64_t hits{0}; + std::uint64_t misses{0}; + std::size_t points{0}; + }; + + struct DMReferenceDisplacementSample { + const mfem::Vector& value; + const mfem::DenseMatrix& gradient; + }; + + struct DMGeometryCache { + bool enabled{true}; + const mfem::Mesh* mesh{nullptr}; + int32_t sequence{-1}; + int32_t nodesSequence{-1}; + std::unordered_map, DMFixedMeshGeometry> points; + uint64_t hits{0}; + uint64_t misses{0}; + }; +} diff --git a/src/include/serif/mapping/domain_mapper.hpp b/src/include/serif/mapping/domain_mapper.hpp new file mode 100644 index 00000000..67e98a0b --- /dev/null +++ b/src/include/serif/mapping/domain_mapper.hpp @@ -0,0 +1,161 @@ +#pragma once + +#include +#include + +#include "mfem/linalg/vector.hpp" +#include "mfem/linalg/densemat.hpp" +#include "mfem/fem/eltrans.hpp" +#include "mfem/fem/intrules.hpp" + +#include "serif/mapping/compactification/kelvin.hpp" +#include "serif/mapping/context.hpp" +#include "serif/mapping/data.hpp" + +#include "serif/utils/user_options/mapping/domain_mapper.hpp" + + +namespace serif::mapping { + class DomainMapper { + public: + DomainMapper( + utils::user_options::DomainMapperOptions options, + std::unique_ptr exterior_map + ); + + DomainMapper(const DomainMapper&) = delete; + DomainMapper& operator=(const DomainMapper&) = delete; + ~DomainMapper(); + DomainMapper(DomainMapper&&) = delete; + DomainMapper& operator=(DomainMapper&&) = delete; + + void set_geometry_cache_enabled(bool enabled) const; + [[nodiscard]] bool is_geometry_cache_enabled() const noexcept; + void clear_geometry_cache() const; + [[nodiscard]] DMGeometryCacheStatistics get_geometry_cache_statistics() const; + + [[nodiscard]] const DMFixedMeshGeometry* get_fixed_mesh_geometry(mfem::ElementTransformation& transformation, const mfem::IntegrationPoint& ip) const; + + [[nodiscard]] MappingStatus evaluate_point( + const ElementMappingData& element_data, + mfem::ElementTransformation &transformation, + const mfem::IntegrationPoint& ip, + DMWorkspace& workspace, + MappingPointContext& context + ) const; + + [[nodiscard]] MappingStatus evaluate_volume( + const ElementMappingData& element_data, + mfem::ElementTransformation& transformation, + const mfem::IntegrationPoint& ip, + DMWorkspace& workspace, + VolumeMappingContext& context, + const DMReferenceDisplacementSample *displacement_sample = nullptr + ) const; + + [[nodiscard]] MappingStatus evaluate_face( + const ElementMappingData& element_data, + mfem::ElementTransformation& transformation, + FaceElementSide side, + const mfem::IntegrationPoint& ip, + DMWorkspace& workspace, + FaceMappingContext& context + ) const; + + [[nodiscard]] MappingStatus evaluate_point_variation( + const ElementMappingData& element_data, + const ElementDisplacementData& direction, + mfem::ElementTransformation& transformation, + const mfem::IntegrationPoint& ip, + const MappingPointContext& context, + DMWorkspace& workspace, + MappingPointVariation& variation + ); + + [[nodiscard]] MappingStatus evaluate_volume_variation( + const ElementMappingData& element_data, + const ElementDisplacementData& direction, + mfem::ElementTransformation& transformation, + const mfem::IntegrationPoint& ip, + const VolumeMappingContext& context, + DMWorkspace& workspace, + VolumeMappingVariation& variation + ); + + [[nodiscard]] MappingStatus evaluate_face_variation( + const ElementMappingData& element_data, + const ElementDisplacementData& direction, + mfem::ElementTransformation& transformation, + FaceElementSide side, + const mfem::IntegrationPoint& ip, + const FaceMappingContext& context, + DMWorkspace& workspace, + FaceMappingVariation& variation + ); + + [[nodiscard]] bool is_compactified_element(const mfem::ElementTransformation& transformation) const noexcept; + [[nodiscard]] int get_dimension() const noexcept; + [[nodiscard]] const compactification::ExteriorDomainMap& get_exterior_map() const noexcept; + private: + [[nodiscard]] MappingStatus evaluate_point_impl( + const ElementMappingData& element_data, + mfem::ElementTransformation& transformation, + const mfem::IntegrationPoint& ip, + DMWorkspace& workspace, + MappingPointContext& context, + const DMFixedMeshGeometry* geometry, + const DMReferenceDisplacementSample* displacement_sample = nullptr + ) const; + + void validate_element_data( + const ElementMappingData& element_data + ) const; + + void evaluate_field( + const ElementDisplacementData& field, + mfem::ElementTransformation& transformation, + const mfem::IntegrationPoint& ip, + DMWorkspace& workspace, + mfem::Vector& value, + mfem::DenseMatrix& jacobian, + const mfem::DenseMatrix *inverse_mesh_jacobian + ) const; + + [[nodiscard]] MappingStatus evaluate_compactification_coordinate( + const ElementCompactificationData& compactification, + mfem::ElementTransformation& transformation, + const mfem::IntegrationPoint& ip, + DMWorkspace& workspace, + CompactificationPointData& point_data, + const mfem::DenseMatrix *inverse_mesh_jacobian + ) const; + + [[nodiscard]] MappingStatus evaluate_point_variation_impl( + const ElementMappingData& element_data, + const ElementDisplacementData& direction, + mfem::ElementTransformation& transformation, + const mfem::IntegrationPoint& ip, + const MappingPointContext& context, + DMWorkspace& workspace, + MappingPointVariation& variation, + const mfem::DenseMatrix *inverse_mesh_jacobian + ) const; + + // These are static since they do not use any class members or workspace. Really these could and + // maybe should be free functions... + [[nodiscard]] static mfem::ElementTransformation& select_face_element_transformation( + mfem::FaceElementTransformations& transformations, + FaceElementSide side + ); + + [[nodiscard]] static const mfem::IntegrationPoint& select_face_element_integration_point( + mfem::FaceElementTransformations& transformations, + FaceElementSide side + ); + + private: + mutable std::unique_ptr m_geometryCache; + utils::user_options::DomainMapperOptions m_options; + std::unique_ptr m_exterior_map; + }; +} \ No newline at end of file diff --git a/src/include/serif/utils/error/codes.hpp b/src/include/serif/utils/error/codes.hpp index 73fa596e..cc9f02dc 100644 --- a/src/include/serif/utils/error/codes.hpp +++ b/src/include/serif/utils/error/codes.hpp @@ -51,6 +51,12 @@ namespace serif::utils::errors { NON_FINITE_MAPPING_RESULT, NON_POSITIVE_DETERMINATE, + INVALID_DOF_COUNT, + INVALID_DOF_ORDERING, + INVALID_COMPACTIFICATION_RANK, + INVALID_COMPACTIFICATION_TYPE, + INVALID_COMPACTIFICATION_DERIVATIVE_VALUE, + IMPROPERLY_ORDERED_REFERENCE_RADII = 450, INVALID_COORDINATE_TOLERANCE, diff --git a/src/include/serif/utils/user_options/mapping/domain_mapper.hpp b/src/include/serif/utils/user_options/mapping/domain_mapper.hpp new file mode 100644 index 00000000..36bdbd43 --- /dev/null +++ b/src/include/serif/utils/user_options/mapping/domain_mapper.hpp @@ -0,0 +1,8 @@ +#pragma once + +namespace serif::utils::user_options { + struct DomainMapperOptions { + int dimension{3}; + int vacuum_element_attribute_ID{3}; + }; +} \ No newline at end of file diff --git a/src/lib/serif/mapping/data.cpp b/src/lib/serif/mapping/data.cpp new file mode 100644 index 00000000..6de660db --- /dev/null +++ b/src/lib/serif/mapping/data.cpp @@ -0,0 +1,153 @@ +#include + +#include "mfem/linalg/vector.hpp" +#include "mfem/fem/fe/fe_base.hpp" +#include "mfem/fem/fespace.hpp" + +#include "serif/mapping/data.hpp" + +#include "serif/utils/error/errors.hpp" + +namespace serif::mapping { + using namespace serif::utils::errors; + using mfem::Ordering; + + ElementDisplacementData::ElementDisplacementData( + const mfem::FiniteElement& element, + mfem::Vector& displacement_dofs, + const Ordering::Type ordering + ) : m_element(&element), m_dimension(0), m_ordering(ordering) { + const int dof_count = element.GetDof(); + + if (dof_count <= 0) { + serif_error("ElementDisplacementData: Displacement elements must have at least one degree of freedom. Got an element with {} dofs", dof_count); + } + if (displacement_dofs.Size() <= 0 || displacement_dofs.Size() % dof_count != 0) { + serif_error("ElementDisplacementData: The displacement vector size ({}) must be a positive multiple of the element degree of freedom count ({}).", displacement_dofs.Size(), dof_count); + } + + m_dimension = displacement_dofs.Size() / dof_count; + m_dof_matrix.SetSize(dof_count, m_dimension); + + switch (ordering) { + case Ordering::byNODES: { + for (int componentID = 0; componentID < m_dimension; ++componentID) { + for (int dofID = 0; dofID < dof_count; ++dofID) { + m_dof_matrix(dofID, componentID) = displacement_dofs(dofID + componentID * dof_count); + } + } + break; + } + case Ordering::byVDIM: { + for (int dofID = 0; dofID < dof_count; ++dofID) { + for (int componentID = 0; componentID < m_dimension; ++componentID) { + m_dof_matrix(dofID, componentID) = displacement_dofs(componentID + dofID * m_dimension); + } + } + break; + } + default: + serif_error("Unsupported ordering"); + } + } + + const mfem::FiniteElement& ElementDisplacementData::get_element() const noexcept { + return *m_element; + } + + const mfem::DenseMatrix& ElementDisplacementData::get_dof_matrix() const noexcept { + return m_dof_matrix; + } + + int ElementDisplacementData::get_dimension() const noexcept { + return m_dimension; + } + + int ElementDisplacementData::get_dof_count() const noexcept { + return m_element->GetDof(); + } + + Ordering::Type ElementDisplacementData::get_ordering() const noexcept { + return m_ordering; + } + + ElementDisplacementData element_displacement_data_from_element_vdofs(const mfem::FiniteElement& element, mfem::Vector& displacement_dofs) { + return ElementDisplacementData(element, displacement_dofs, Ordering::byNODES); + } + + ElementCompactificationData::ElementCompactificationData( + const mfem::FiniteElement& element, + const mfem::Vector& dofs + ) : m_element(&element), m_dofs(dofs) { + if (element.GetRangeType() != mfem::FiniteElement::SCALAR) { + serif_error("Compactification coordinate requires a scalar finite element."); + } + + if (element.GetMapType() != mfem::FiniteElement::VALUE) { + serif_error("Compactification coordinate requires a value-mapped scalar finite element"); + } + + if (element.GetDerivType() != mfem::FiniteElement::GRAD) { + serif_error("Compactification coordinate finite element must provide a gradient"); + } + + if (element.GetDof() <= 0) { + serif_error("Compactification coordinatge finite element must have a positive non-zero number of degrees of freedom. Recived {}", element.GetDof()); + } + + if (dofs.Size() != element.GetDof()) { + serif_error("Compactification coordinate degree of freedom count ({}) does not match its finite element dof count ({})", dofs.Size(), element.GetDof()); + } + } + + const mfem::FiniteElement& ElementCompactificationData::get_element() const noexcept { + return *m_element; + } + + const mfem::Vector& ElementCompactificationData::get_dofs() const noexcept { + return m_dofs; + } + + int ElementCompactificationData::get_dof_count() const noexcept { + return m_dofs.Size(); + } + + DMWorkspace::DMWorkspace(const int dimension) { + set_dimension(dimension); + } + + void DMWorkspace::set_dimension(const int dimension) { + if (dimension <= 0) { + serif_error("Domain mapping workspace dimension must be positive and finite. Received {}", dimension); + } + + m_dimension = dimension; + + m_field_value.SetSize(dimension); + m_field_jacobian.SetSize(dimension, dimension); + m_reference_field_jacobian.SetSize(dimension, dimension); + + m_compactification_point.compactification_coordinate = 0.0; + m_compactification_point.coordinate_gradient.SetSize(dimension); + + m_reference_normal.SetSize(dimension); + m_mapped_normal.SetSize(dimension); + m_full_element_jacobian.SetSize(dimension, dimension); + + m_vector_temp.SetSize(dimension); + m_matrix_temp_1.SetSize(dimension, dimension); + m_matrix_temp_2.SetSize(dimension, dimension); + + m_exterior_result.physical_position.SetSize(dimension); + m_exterior_result.mapping_jacobian.SetSize(dimension, dimension); + + m_exterior_variation.physical_position_variation.SetSize(dimension); + m_exterior_variation.mapping_jacobian_variation.SetSize(dimension, dimension); + } + + int DMWorkspace::get_dimension() const noexcept { + return m_dimension; + } + + +} diff --git a/src/meson.build b/src/meson.build index 6d51e18e..bea4de72 100644 --- a/src/meson.build +++ b/src/meson.build @@ -6,6 +6,7 @@ serif_sources = files( 'lib/serif/discretization/quadrature/policy.cpp', 'lib/serif/eos/models/polytropic.cpp', 'lib/serif/mapping/compactification/kelvin.cpp', + 'lib/serif/mapping/data.cpp', ) serif_lib = library( diff --git a/subprojects/packagefiles/stroid/build-check/meson.build b/subprojects/packagefiles/stroid/build-check/meson.build new file mode 100644 index 00000000..f78e7b15 --- /dev/null +++ b/subprojects/packagefiles/stroid/build-check/meson.build @@ -0,0 +1,7 @@ +# A subproject must not add global arguments after the parent defines targets. +cpp = meson.get_compiler('cpp') +add_project_arguments( + cpp.get_supported_arguments('-Wno-unused-command-line-argument'), + language: 'cpp', +) +subdir('CPPC') diff --git a/subprojects/packagefiles/stroid/build-config/mfem/meson.build b/subprojects/packagefiles/stroid/build-config/mfem/meson.build new file mode 100644 index 00000000..7315b9fd --- /dev/null +++ b/subprojects/packagefiles/stroid/build-config/mfem/meson.build @@ -0,0 +1,3 @@ +# The parent overrides this dependency with its selected MFEM build. Building +# another copy here can change class layouts and cause crashes across the ABI. +mfem_dep = dependency('mfem', version: '>=4.10') diff --git a/tests/sandbox/mapping_sandbox.cpp b/tests/sandbox/mapping_sandbox.cpp index ed16a5f6..9ea7816a 100644 --- a/tests/sandbox/mapping_sandbox.cpp +++ b/tests/sandbox/mapping_sandbox.cpp @@ -5,6 +5,7 @@ #include "serif/mapping/compactification/compactification.hpp" #include "serif/mapping/status.hpp" #include "serif/mapping/compactification/kelvin.hpp" +#include "serif/mapping/domain_mapper.hpp" #include "serif/utils/user_options/mapping/compactification/kelvin.hpp" int main() {