From 52c42ec7cd1e17dc7406debcebebc3bbec80b1bc Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Wed, 23 Sep 2026 09:29:05 -0400 Subject: [PATCH] feat(mapper): finished porting over mapper --- src/include/serif/mapping/coefficients.hpp | 34 + .../compactification/compactification.hpp | 8 +- .../serif/mapping/compactification/kelvin.hpp | 1 + src/include/serif/mapping/context.hpp | 10 +- src/include/serif/mapping/data.hpp | 6 +- src/include/serif/mapping/domain_mapper.hpp | 29 +- src/include/serif/mapping/evaluator.hpp | 93 +++ src/include/serif/mapping/transformations.hpp | 87 +++ src/include/serif/mapping/volume_cache.hpp | 28 + src/include/serif/utils/checks/finite.hpp | 8 + src/include/serif/utils/error/codes.hpp | 24 +- src/include/serif/utils/error/errors.hpp | 25 +- src/lib/serif/mapping/coefficients.cpp | 22 + .../serif/mapping/compactification/kelvin.cpp | 4 + src/lib/serif/mapping/data.cpp | 2 +- src/lib/serif/mapping/domain_mapper.cpp | 658 +++++++++++++++++- src/lib/serif/mapping/evaluator.cpp | 149 ++++ src/lib/serif/mapping/transformations.cpp | 136 ++++ src/lib/serif/mapping/volume_cache.cpp | 112 +++ src/meson.build | 7 + tests/sandbox/mapping_sandbox.cpp | 1 + 21 files changed, 1414 insertions(+), 30 deletions(-) create mode 100644 src/include/serif/mapping/coefficients.hpp create mode 100644 src/include/serif/mapping/evaluator.hpp create mode 100644 src/include/serif/mapping/transformations.hpp create mode 100644 src/include/serif/mapping/volume_cache.hpp create mode 100644 src/lib/serif/mapping/coefficients.cpp create mode 100644 src/lib/serif/mapping/evaluator.cpp create mode 100644 src/lib/serif/mapping/transformations.cpp create mode 100644 src/lib/serif/mapping/volume_cache.cpp diff --git a/src/include/serif/mapping/coefficients.hpp b/src/include/serif/mapping/coefficients.hpp new file mode 100644 index 00000000..4020496d --- /dev/null +++ b/src/include/serif/mapping/coefficients.hpp @@ -0,0 +1,34 @@ +#pragma once + +#include + +#include "mfem/fem/gridfunc.hpp" +#include "mfem/fem/coefficient.hpp" +#include "mfem/fem/eltrans.hpp" +#include "mfem/fem/intrules.hpp" +#include "mfem/linalg/vector.hpp" + +#include "serif/mapping/domain_mapper.hpp" +#include "serif/mapping/evaluator.hpp" + +namespace serif::mapping { + class PhysicalPositionFunctionCoefficient : public mfem::Coefficient { + public: + using Func = std::function; + + PhysicalPositionFunctionCoefficient( + const DomainMapper& mapper, + const mfem::GridFunction& displacement, + const mfem::GridFunction& compactification_coordinate, + Func func + ); + + double Eval( + mfem::ElementTransformation& transformation, + const mfem::IntegrationPoint& ip + ) override; + private: + Func m_f; + GridFunctionMappingEvaluator m_mapping; + }; +} diff --git a/src/include/serif/mapping/compactification/compactification.hpp b/src/include/serif/mapping/compactification/compactification.hpp index d960d6bd..47db8128 100644 --- a/src/include/serif/mapping/compactification/compactification.hpp +++ b/src/include/serif/mapping/compactification/compactification.hpp @@ -1,6 +1,7 @@ #pragma once #include +#include #include "serif/mapping/status.hpp" #include "serif/mapping/compactification/input_result.hpp" @@ -12,21 +13,18 @@ namespace serif::mapping::compactification { public: virtual ~ExteriorDomainMap() = default; - [[nodiscard]] virtual bool is_affine_in_displacement() const noexcept; - + [[nodiscard]] virtual bool is_affine_in_displacement() const noexcept = 0; [[nodiscard]] virtual MappingStatus evaluate( const ExteriorMapInput &input, ExteriorMapResult &result ) const = 0; - [[nodiscard]] virtual MappingStatus evaluate_variation( const ExteriorMapInput &input, const ExteriorMapResult &result, const ExteriorMapDirection &direction, ExteriorMapVariation &variation ) const = 0; - [[nodiscard]] virtual std::string_view get_name() const noexcept = 0; - + [[nodiscard]] virtual std::unique_ptr clone() const = 0; }; } \ No newline at end of file diff --git a/src/include/serif/mapping/compactification/kelvin.hpp b/src/include/serif/mapping/compactification/kelvin.hpp index c8a1689c..fcb69932 100644 --- a/src/include/serif/mapping/compactification/kelvin.hpp +++ b/src/include/serif/mapping/compactification/kelvin.hpp @@ -27,6 +27,7 @@ namespace serif::mapping::compactification { ) const override; [[nodiscard]] std::string_view get_name() const noexcept override; + [[nodiscard]] std::unique_ptr clone() const override; [[nodiscard]] double get_reference_stellar_radius() const noexcept; [[nodiscard]] double get_reference_infinity_radius() const noexcept; diff --git a/src/include/serif/mapping/context.hpp b/src/include/serif/mapping/context.hpp index faea57cd..087a3afc 100644 --- a/src/include/serif/mapping/context.hpp +++ b/src/include/serif/mapping/context.hpp @@ -3,11 +3,15 @@ #include "mfem/linalg/densemat.hpp" #include "mfem/linalg/vector.hpp" +// TODO: There is likely some major readability improvements that could be gained by naming these things +// less confusing names + +// TODO [Emily]: Put more comments in this file describing what the various terms are namespace serif::mapping { struct VolumeQuadratureContext { - mfem::DenseMatrix J_inv; - double detJ; - double weight; + mfem::DenseMatrix J_element_to_physical_inverse; + double mapping_determinate; + double physical_integration_weight; }; struct FaceQuadratureContext { diff --git a/src/include/serif/mapping/data.hpp b/src/include/serif/mapping/data.hpp index b9f3d3f8..29a49269 100644 --- a/src/include/serif/mapping/data.hpp +++ b/src/include/serif/mapping/data.hpp @@ -6,7 +6,7 @@ #include "mfem/mesh/mesh.hpp" #include -#include +#include #include #include "serif/mapping/compactification/kelvin.hpp" @@ -73,7 +73,7 @@ namespace serif::mapping { [[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_reference_dshape; mfem::DenseMatrix m_mesh_dshape; mfem::DenseMatrix m_reference_field_jacobian; mfem::Vector m_field_value; @@ -120,7 +120,7 @@ namespace serif::mapping { const mfem::Mesh* mesh{nullptr}; long sequence{-1}; // We explicitly use the non explicltiy sized versions here since that what MFEM uses fort these members. Alternative we could use a int64_t and verify that on the compiled system that is the same as a long long nodesSequence{-1}; - std::unordered_map, DMFixedMeshGeometry> points; + std::map, DMFixedMeshGeometry> points; // Currently this is an regular map since to use an unordered map we would need to add a std::hash specialization for tuple. Some small future optimizations may be available if we implement a fast hashing algorithm for that type. 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 index b78564af..ad9ad703 100644 --- a/src/include/serif/mapping/domain_mapper.hpp +++ b/src/include/serif/mapping/domain_mapper.hpp @@ -19,11 +19,6 @@ 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; @@ -31,6 +26,14 @@ namespace serif::mapping { DomainMapper(DomainMapper&&) = default; DomainMapper& operator=(DomainMapper&&) = default; + + DomainMapper( + utils::user_options::DomainMapperOptions options, + std::unique_ptr exterior_map + ); + + std::unique_ptr clone() const; + void set_geometry_cache_enabled(bool enabled) const; [[nodiscard]] bool is_geometry_cache_enabled() const noexcept; void clear_geometry_cache() const; @@ -57,7 +60,7 @@ namespace serif::mapping { [[nodiscard]] MappingStatus evaluate_face( const ElementMappingData& element_data, - mfem::ElementTransformation& transformation, + mfem::FaceElementTransformations& transformation, FaceElementSide side, const mfem::IntegrationPoint& ip, DMWorkspace& workspace, @@ -72,7 +75,7 @@ namespace serif::mapping { const MappingPointContext& context, DMWorkspace& workspace, MappingPointVariation& variation - ); + ) const; [[nodiscard]] MappingStatus evaluate_volume_variation( const ElementMappingData& element_data, @@ -82,18 +85,18 @@ namespace serif::mapping { const VolumeMappingContext& context, DMWorkspace& workspace, VolumeMappingVariation& variation - ); + ) const; [[nodiscard]] MappingStatus evaluate_face_variation( const ElementMappingData& element_data, const ElementDisplacementData& direction, - mfem::ElementTransformation& transformation, + mfem::FaceElementTransformations &transformation, FaceElementSide side, const mfem::IntegrationPoint& ip, const FaceMappingContext& context, DMWorkspace& workspace, FaceMappingVariation& variation - ); + ) const; [[nodiscard]] bool is_compactified_element(const mfem::ElementTransformation& transformation) const noexcept; [[nodiscard]] int get_dimension() const noexcept; @@ -113,7 +116,7 @@ namespace serif::mapping { const ElementMappingData& element_data ) const; - void evaluate_field( + void evaluate_displacement_field( const ElementDisplacementData& field, mfem::ElementTransformation& transformation, const mfem::IntegrationPoint& ip, @@ -146,12 +149,12 @@ namespace serif::mapping { // 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, + const mfem::FaceElementTransformations& transformations, FaceElementSide side ); [[nodiscard]] static const mfem::IntegrationPoint& select_face_element_integration_point( - mfem::FaceElementTransformations& transformations, + const mfem::FaceElementTransformations& transformations, FaceElementSide side ); diff --git a/src/include/serif/mapping/evaluator.hpp b/src/include/serif/mapping/evaluator.hpp new file mode 100644 index 00000000..80001d31 --- /dev/null +++ b/src/include/serif/mapping/evaluator.hpp @@ -0,0 +1,93 @@ +#pragma once +#include + +#include "mfem/fem/gridfunc.hpp" +#include "mfem/fem/eltrans.hpp" +#include "mfem/fem/intrules.hpp" +#include "mfem/linalg/vector.hpp" +#include "mfem/general/array.hpp" + +#include "serif/mapping/domain_mapper.hpp" +#include "serif/mapping/context.hpp" +#include "serif/mapping/status.hpp" + +namespace serif::mapping { + // Note this is NOT thread safe. + class GridFunctionMappingEvaluator { + public: + /* + * The evaluator references the supplied grid functions and caches copies of + * their element-local DOFs. You need to make sure you call InvalidateCache() or Refresh() after either + * grid function's values are modified. Changes to the finite-element-space sequence + * should be detected automatically. + * + * This object owns mutable workspace and cache state and *is not* thread-safe. + */ + GridFunctionMappingEvaluator( + const DomainMapper& mapper, + const mfem::GridFunction& displacement, + const mfem::GridFunction& compactification_coordinate + ); + + void invalidate_cache() noexcept; + void refresh(); + + [[nodiscard]] MappingStatus evaluate_point( + mfem::ElementTransformation &transformation, + const mfem::IntegrationPoint &integration_point, + MappingPointContext &context + ); + + [[nodiscard]] MappingStatus evaluate_volume( + mfem::ElementTransformation &transformation, + const mfem::IntegrationPoint &integration_point, + VolumeMappingContext &context + ); + + [[nodiscard]] MappingStatus evaluate_face( + mfem::FaceElementTransformations &transformation, + FaceElementSide side, + const mfem::IntegrationPoint &integration_point, + FaceMappingContext &context + ); + + [[nodiscard]] VolumeQuadratureContext get_quadrature_context( + mfem::ElementTransformation &transformation, + const mfem::IntegrationPoint &integration_point + ); + + [[nodiscard]] FaceQuadratureContext get_face_quadrature_context( + mfem::FaceElementTransformations &transformation, + const mfem::IntegrationPoint &integration_point, + FaceElementSide side = FaceElementSide::ELEMENT_1 + ); + + void get_physical_point( + mfem::ElementTransformation &transformation, + const mfem::IntegrationPoint &integration_point, + mfem::Vector &physical_position + ); + private: + void validate_field_bindings() const; + [[nodiscard]] bool invalidate_for_changed_spaces(); + void load_element(int element_id); + + private: + const DomainMapper &m_mapper; + const mfem::GridFunction &m_displacement; + const mfem::GridFunction &m_compactification_coordinate; + const mfem::FiniteElementSpace *m_displacement_space; + const mfem::FiniteElementSpace *m_compactification_space; + long m_displacement_space_sequence; + long m_compactification_space_sequence; + DMWorkspace m_workspace; + + mfem::Array m_displacement_dofs; + mfem::Array m_compactification_dofs; + mfem::Vector m_element_displacement; + mfem::Vector m_element_compactification; + std::unique_ptr m_displacement_data; + std::unique_ptr m_compactification_data; + int m_cached_element_id{-1}; + }; +} \ No newline at end of file diff --git a/src/include/serif/mapping/transformations.hpp b/src/include/serif/mapping/transformations.hpp new file mode 100644 index 00000000..02858ec6 --- /dev/null +++ b/src/include/serif/mapping/transformations.hpp @@ -0,0 +1,87 @@ +#pragma once + +#include "serif/mapping/variation.hpp" +#include "serif/mapping/context.hpp" + +namespace serif::mapping { + void map_hdiv_flux_to_physical( + const MappingPointContext& context, + const mfem::Vector& reference_flux, + mfem::Vector& physical_flux + ); + + void map_physical_flux_to_hdiv_reference( + const MappingPointContext& context, + const mfem::Vector& physical_flux, + mfem::Vector& reference_flux + ); + + void map_reference_gradient_to_physical( + const MappingPointContext& context, + const mfem::Vector& reference_gradient, + mfem::Vector& physical_gradient + ); + + void map_physical_gradient_to_reference( + const MappingPointContext& context, + const mfem::Vector& physical_gradient, + mfem::Vector& reference_gradient + ); + + void map_reference_vector_gradient_to_physical( + const MappingPointContext& context, + const mfem::DenseMatrix& reference_vector_gradient, + mfem::DenseMatrix& physical_vector_gradient + ); + + void map_physical_vector_gradient_to_reference( + const MappingPointContext& context, + const mfem::DenseMatrix& physical_vector_gradient, + mfem::DenseMatrix& reference_vector_gradient + ); + + [[nodiscard]] double map_hdiv_divergence_to_physical( + const MappingPointContext& context, + double reference_divergence + ); + + void compute_hdiv_mass_tensor( + const MappingPointContext& context, + mfem::DenseMatrix& mass_tensor + ); + + void compute_scalar_diffusion_tensor( + const MappingPointContext& context, + mfem::DenseMatrix& diffusion_tensor + ); + + void map_hcurl_field_to_physical( + const MappingPointContext& context, + const mfem::Vector& reference_field, + mfem::Vector& physical_field + ); + + void map_physical_field_to_hcurl_reference( + const MappingPointContext& context, + const mfem::Vector& physical_field, + mfem::Vector& reference_field + ); + + void map_hcurl_curl_to_physical( + const MappingPointContext& context, + const mfem::Vector &reference_curl, + mfem::Vector &physical_curl + ); + + void map_physical_curl_to_hcurl_reference( + const MappingPointContext& context, + const mfem::Vector& physical_curl, + mfem::Vector& reference_curl + ); + + void compute_hdiv_mass_tensor_variation( + const MappingPointContext& context, + const MappingPointVariation& variation, + mfem::DenseMatrix& mass_tensor_variation + ); +} \ No newline at end of file diff --git a/src/include/serif/mapping/volume_cache.hpp b/src/include/serif/mapping/volume_cache.hpp new file mode 100644 index 00000000..7c38105c --- /dev/null +++ b/src/include/serif/mapping/volume_cache.hpp @@ -0,0 +1,28 @@ +#pragma once + +#include "mfem/linalg/densemat.hpp" + +#include "serif/mapping/context.hpp" +#include + +namespace serif::mapping { + class VolumeMappingCache { + public: + void set_size(int point_count, int dimension); + void store(int point, VolumeMappingContext& context); + void load(int point, VolumeMappingContext& context) const; + void load_inverse_jacobian(int point, mfem::DenseMatrix& inverse) const; + [[nodiscard]] int get_point_count() const noexcept; + [[nodiscard]] int get_dimension() const noexcept; + + private: + void validate_point(int point) const; + [[nodiscard]] const double* get_point_data(int point) const; + + std::vector m_data; + int m_point_count{0}; + int m_dimension{0}; + int m_point_stride{0}; + }; + +} \ No newline at end of file diff --git a/src/include/serif/utils/checks/finite.hpp b/src/include/serif/utils/checks/finite.hpp index a036dde8..940b745a 100644 --- a/src/include/serif/utils/checks/finite.hpp +++ b/src/include/serif/utils/checks/finite.hpp @@ -5,6 +5,8 @@ #include "mfem/linalg/vector.hpp" #include "mfem/linalg/densemat.hpp" +#include "serif/utils/misc/concepts/numeric.hpp" + namespace serif::utils::checks { // If we have a new backend / new vector and matrix type in the future we can add new overloads inline bool is_finite(const mfem::Vector& vector) { @@ -28,4 +30,10 @@ namespace serif::utils::checks { return true; } + + template + requires requires (Scalar x) {{std::isfinite(x)} -> std::convertible_to; } + bool is_finite(const Scalar scalar) { + return std::isfinite(scalar); + } } \ 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 cc9f02dc..68cf519a 100644 --- a/src/include/serif/utils/error/codes.hpp +++ b/src/include/serif/utils/error/codes.hpp @@ -15,7 +15,7 @@ namespace serif::utils::errors { // stay consistent here. // At some point someone should go through and come up with some well defined schema for all this - enum class SERiFErrorCode : uint16_t { + enum class SERiFErrorCode : uint16_t { // 16 bit backing type, so 65535 possible error codes, which should be far far far more than enough for now UNKNOWN_ERROR = 0, // RUNTIME FAILURES (>100) @@ -30,6 +30,7 @@ namespace serif::utils::errors { // Jacobian Errors JACOBIAN_ERROR = 200, INVERTED_JACOBIAN, + NON_FINITE_DETERMINATE, // EOS Errors EOS_ERROR = 300, @@ -50,16 +51,35 @@ namespace serif::utils::errors { OUTSIDE_REFERENCE_DOMAIN, NON_FINITE_MAPPING_RESULT, NON_POSITIVE_DETERMINATE, + INVALID_MAPPING, + MISC_MAPPING_ERROR=410, INVALID_DOF_COUNT, + INVALID_ELEMENT_GEOMETRY, INVALID_DOF_ORDERING, INVALID_COMPACTIFICATION_RANK, INVALID_COMPACTIFICATION_TYPE, - INVALID_COMPACTIFICATION_DERIVATIVE_VALUE, + INVALID_COMPACTIFICATION_DERIVATIVE_TYPE, + INVALID_VACUUM_ELEMENT_ATTRIBUTE_ID, + MISSING_EXTERIOR_MAP, + COMPACTIFICATION_MISMATCH, + MISSING_FACE_ELEMENT, + JACOBIAN_SHAPE_ERROR, + MISSING_DISPLACEMENT_FINITE_ELEMENT_SPACE, + MISMATCHED_MESHING_MAPS, IMPROPERLY_ORDERED_REFERENCE_RADII = 450, INVALID_COORDINATE_TOLERANCE, + EVALUATOR_ERROR=460, + CACHE_FE_BINDING_ERROR, + INVALID_ELEMENT_ID_ERROR, + ELEMENT_NULLPTR_ERROR, + + //As yet uncategorized errors + MATRIX_SHAPE_ERROR = 500, + OUT_OF_RANGE_ERROR, + // COMPILE TIME FAILURES (>1000) // Type Invariant Failures FAILED_TYPE_INVARIANT = 1000 diff --git a/src/include/serif/utils/error/errors.hpp b/src/include/serif/utils/error/errors.hpp index 5d2f0eae..f771a616 100644 --- a/src/include/serif/utils/error/errors.hpp +++ b/src/include/serif/utils/error/errors.hpp @@ -38,7 +38,7 @@ namespace serif::utils::errors { throw ExceptionT(error_msg); } - template + template void serif_assert(const bool condition, const std::string_view message, auto... args) { // collapse to a noop if NDEBUG is defined #ifndef NDEBUG @@ -48,7 +48,7 @@ namespace serif::utils::errors { #endif if constexpr(is_debug_build) { if (!condition) { - serif_error(message, args...); + serif_error(message, args...); } } } @@ -75,4 +75,25 @@ namespace serif::utils::errors { } } } + + // The difference between serif_assert and serif_ensure is that serif_assert is a noop in release builds, while serif_ensure is always active. + // Use an ensure when you want to always check things. Use an assert when you want to check things that are true bugs + // and should not manifest in release builds. + template + void serif_ensure(const bool condition, const std::string_view message, auto... args) { + if (!condition) { + serif_error(message, args...); + } + } + + template + constexpr void serif_ensure() { + if constexpr (!condition) { + // Trigger an intentional error to display the message and code at compile time + + // Note if we move to C++26 we can replace this with a more standard and less convoluted + // formated static_assert. + [[maybe_unused]] std::size_t error_kick = sizeof(serif_error_printer); + } + } } \ No newline at end of file diff --git a/src/lib/serif/mapping/coefficients.cpp b/src/lib/serif/mapping/coefficients.cpp new file mode 100644 index 00000000..ba9f3c2b --- /dev/null +++ b/src/lib/serif/mapping/coefficients.cpp @@ -0,0 +1,22 @@ +#include "serif/mapping/coefficients.hpp" +#include "serif/utils/error/errors.hpp" + +namespace serif::mapping { + using namespace serif::utils::errors; + + PhysicalPositionFunctionCoefficient::PhysicalPositionFunctionCoefficient( + const DomainMapper &mapper, + const mfem::GridFunction &displacement, + const mfem::GridFunction &compactification_coordinate, + Func func + ) : + m_f(std::move(func)), + m_mapping(mapper, displacement, compactification_coordinate) { } + + double PhysicalPositionFunctionCoefficient::Eval(mfem::ElementTransformation &transformation, const mfem::IntegrationPoint &ip) { + transformation.SetIntPoint(&ip); + MappingPointContext context; + serif_ensure(m_mapping.evaluate_point(transformation, ip, context) == MappingStatus::VALID, "Mapping evaluation failed during PhysicalPositionCoefficient evaluation."); + return m_f(context.physical_position); + } +} diff --git a/src/lib/serif/mapping/compactification/kelvin.cpp b/src/lib/serif/mapping/compactification/kelvin.cpp index b7fc57e3..7fbf2b6b 100644 --- a/src/lib/serif/mapping/compactification/kelvin.cpp +++ b/src/lib/serif/mapping/compactification/kelvin.cpp @@ -203,6 +203,10 @@ namespace serif::mapping::compactification { return MappingStatus::VALID; } + std::unique_ptr KelvinCompactification::clone() const { + return std::make_unique(m_options); + } + MappingStatus KelvinCompactification::compute_radial_factors(double compactification_coordinate, RadialFactors& radial_factors) const { if (!std::isfinite(compactification_coordinate)) { return MappingStatus::NON_FINITE_MAPPING_INPUT; diff --git a/src/lib/serif/mapping/data.cpp b/src/lib/serif/mapping/data.cpp index 6de660db..7106a71f 100644 --- a/src/lib/serif/mapping/data.cpp +++ b/src/lib/serif/mapping/data.cpp @@ -88,7 +88,7 @@ namespace serif::mapping { } if (element.GetDerivType() != mfem::FiniteElement::GRAD) { - serif_error("Compactification coordinate finite element must provide a gradient"); + serif_error("Compactification coordinate finite element must provide a gradient"); } if (element.GetDof() <= 0) { diff --git a/src/lib/serif/mapping/domain_mapper.cpp b/src/lib/serif/mapping/domain_mapper.cpp index 4ae0d397..eec9a556 100644 --- a/src/lib/serif/mapping/domain_mapper.cpp +++ b/src/lib/serif/mapping/domain_mapper.cpp @@ -1,10 +1,51 @@ +#include "mfem/linalg/densemat.hpp" + #include "serif/mapping/domain_mapper.hpp" #include "serif/utils/error/errors.hpp" +#include "serif/utils/checks/finite.hpp" + #include +#include + +// Anonymous namespace for translation unit-local helper functions and types +namespace { + // TODO: Put something in here using variadic args to validate multiple things satisfy some predicate +} namespace serif::mapping { using namespace serif::utils::errors; + using namespace serif::utils::checks; + + DomainMapper::DomainMapper( + utils::user_options::DomainMapperOptions options, + std::unique_ptr exterior_map + ) : + m_geometryCache(std::make_unique()), + m_options(options), + m_exterior_map(std::move(exterior_map)) { + if (m_options.dimension <= 0) { + serif_error("Domain Mapper requires a positive, non-zero, dimension. Received dimension: {}", m_options.dimension); + } + + if (m_options.vacuum_element_attribute_ID <= 0) { + serif_error("Domain Mapper requires a positive, non-zero, vacuum element attribute ID. Received attribute ID: {}", m_options.vacuum_element_attribute_ID); + } + + if (!m_exterior_map) { + serif_error("Domain Mapper requires a valid exterior map. One was not provided or a nullptr was provided."); + } + } + + std::unique_ptr DomainMapper::clone() const { + if (!m_exterior_map) { + serif_error("Domain Mapper clone failed because the exterior map is missing or a nullptr."); + } + + auto result = std::make_unique(m_options, m_exterior_map->clone()); + result->set_geometry_cache_enabled(is_geometry_cache_enabled()); + return result; + } void DomainMapper::set_geometry_cache_enabled(const bool enabled) const { m_geometryCache->enabled = enabled; @@ -66,9 +107,624 @@ namespace serif::mapping { } const auto key = std::tuple{transformation.ElementNo, coordinates[0], coordinates[1], coordinates[2]}; + const auto iter = cache.points.find(key); - const auto found = cache.points.find(key); + if (iter != cache.points.end()) { + ++cache.hits; + return &iter->second; + } + DMFixedMeshGeometry geometry; + transformation.Transform(ip, geometry.position); + geometry.jacobian = transformation.Jacobian(); + geometry.inverse_jacobian = transformation.InverseJacobian(); + geometry.weight = transformation.Weight(); + ++cache.misses; + + if (!is_finite(geometry.position) || !is_finite(geometry.jacobian) || !is_finite(geometry.inverse_jacobian) || !is_finite(geometry.weight)) { + return nullptr; + } + + cache.points.emplace(key, std::move(geometry)); + return &cache.points.at(key); + } + + MappingStatus DomainMapper::evaluate_point( + const ElementMappingData &element_data, + mfem::ElementTransformation &transformation, + const mfem::IntegrationPoint &ip, + DMWorkspace &workspace, + MappingPointContext &context + ) const { + const auto* meshGeometry = get_fixed_mesh_geometry(transformation, ip); + return evaluate_point_impl(element_data, transformation, ip, workspace, context, meshGeometry); + } + + MappingStatus DomainMapper::evaluate_volume( + const ElementMappingData &element_data, + mfem::ElementTransformation &transformation, + const mfem::IntegrationPoint &ip, + DMWorkspace &workspace, + VolumeMappingContext &context, + const DMReferenceDisplacementSample *displacement_sample + ) const { + const auto *meshGeometry = get_fixed_mesh_geometry(transformation, ip); + const MappingStatus point_status = evaluate_point_impl( + element_data, + transformation, + ip, + workspace, + context.mapping, + meshGeometry, + displacement_sample + ); + + if (point_status != MappingStatus::VALID) { + return point_status; + } + + mfem::Mult(context.mapping.mapping_jacobian,meshGeometry ? meshGeometry->jacobian : transformation.Jacobian(), workspace.m_full_element_jacobian); + + context.quadrature.J_element_to_physical_inverse.SetSize(m_options.dimension, m_options.dimension); + mfem::CalcInverse(workspace.m_full_element_jacobian, context.quadrature.J_element_to_physical_inverse); + + context.quadrature.mapping_determinate = context.mapping.mapping_determinate; + context.quadrature.physical_integration_weight = ip.weight * (meshGeometry ? meshGeometry->weight : transformation.Weight()) * context.mapping.mapping_determinate; + + if (!is_finite(context.quadrature.J_element_to_physical_inverse) || !std::isfinite(context.quadrature.physical_integration_weight)) { + return MappingStatus::NON_FINITE_MAPPING_RESULT; + } + if (context.quadrature.physical_integration_weight <= 0.0) { + return MappingStatus::NON_POSITIVE_DETERMINATE; + } + + return MappingStatus::VALID; + } + + MappingStatus DomainMapper::evaluate_face( + const ElementMappingData &element_data, + mfem::FaceElementTransformations &transformation, + const FaceElementSide side, + const mfem::IntegrationPoint &ip, + DMWorkspace &workspace, + FaceMappingContext &context + ) const { + transformation.SetAllIntPoints(&ip); + mfem::ElementTransformation& element_transformation = select_face_element_transformation(transformation, side); + const mfem::IntegrationPoint& element_integration_point = select_face_element_integration_point(transformation, side); + + const MappingStatus point_status = evaluate_point(element_data, element_transformation, element_integration_point, workspace, context.mapping); + if (point_status != MappingStatus::VALID) { + return point_status; + } + + // Get the normal vector from the face in reference space + workspace.m_reference_normal.SetSize(m_options.dimension); + mfem::CalcOrtho(transformation.Jacobian(), workspace.m_reference_normal); + if (side == FaceElementSide::ELEMENT_2) { + workspace.m_reference_normal *= -1.0; // Flip the normal for the second element (e.g. the "other" side of the face) + } + + const double reference_normal_magnitude = workspace.m_reference_normal.Norml2(); + if (!is_finite(reference_normal_magnitude) || reference_normal_magnitude <= 0.0) { + return MappingStatus::NON_FINITE_MAPPING_RESULT; + } + context.reference_normal.SetSize(m_options.dimension); + context.reference_normal = workspace.m_reference_normal; + context.reference_normal /= reference_normal_magnitude; + + // Get the normal vector from the face in physical space + context.mapping.inverse_mapping_jacobian.MultTranspose(workspace.m_reference_normal, workspace.m_mapped_normal); + workspace.m_mapped_normal *= context.mapping.mapping_determinate; + + const double mapped_normal_magnitude = workspace.m_mapped_normal.Norml2(); + if (!is_finite(mapped_normal_magnitude) || mapped_normal_magnitude <= 0.0) { + return MappingStatus::NON_FINITE_MAPPING_RESULT; + } + context.quadrature.normal.SetSize(m_options.dimension); + context.quadrature.normal = workspace.m_mapped_normal; + context.quadrature.normal /= mapped_normal_magnitude; + + context.reference_surface_weight = ip.weight * reference_normal_magnitude; + context.physical_surface_weight = ip.weight * mapped_normal_magnitude; + + context.quadrature.ds = context.reference_surface_weight; + context.quadrature.v_dot_n_scale = mapped_normal_magnitude / reference_normal_magnitude; + + if (!is_finite(context.quadrature.normal) || !is_finite(context.reference_surface_weight) || !is_finite(context.physical_surface_weight) || !is_finite(context.quadrature.v_dot_n_scale)) { + return MappingStatus::NON_FINITE_MAPPING_RESULT; + } + + return MappingStatus::VALID; + } + + MappingStatus DomainMapper::evaluate_point_variation( + const ElementMappingData &element_data, + const ElementDisplacementData &direction, + mfem::ElementTransformation &transformation, + const mfem::IntegrationPoint &ip, + const MappingPointContext &context, + DMWorkspace &workspace, + MappingPointVariation &variation + ) const { + return evaluate_point_variation_impl(element_data, direction, transformation, ip, context, workspace, variation, nullptr); + } + + MappingStatus DomainMapper::evaluate_volume_variation( + const ElementMappingData &element_data, + const ElementDisplacementData &direction, + mfem::ElementTransformation &transformation, + const mfem::IntegrationPoint &ip, + const VolumeMappingContext &context, + DMWorkspace &workspace, + VolumeMappingVariation &variation + ) const { + mfem::Mult(context.quadrature.J_element_to_physical_inverse, context.mapping.mapping_jacobian, workspace.m_matrix_temp_2); + const MappingStatus point_status = evaluate_point_variation_impl(element_data, direction, transformation, ip, context.mapping, workspace, variation.mapping, &workspace.m_matrix_temp_2); + if (point_status != MappingStatus::VALID) { + return point_status; + } + + mfem::Mult(context.quadrature.J_element_to_physical_inverse, variation.mapping.mapping_jacobian_variation, workspace.m_matrix_temp_1); + + variation.inverse_element_jacobian_variation.SetSize(m_options.dimension, m_options.dimension); + mfem::Mult(workspace.m_matrix_temp_1, context.mapping.inverse_mapping_jacobian, variation.inverse_element_jacobian_variation); + + variation.inverse_element_jacobian_variation *= -1.0; // Account for the negative sign in the derivative of the inverse of a matrix + + variation.weight_variation = (context.quadrature.physical_integration_weight / context.mapping.mapping_determinate) * variation.mapping.mapping_determinate_variation; + + if (!is_finite(variation.inverse_element_jacobian_variation) || !is_finite(variation.weight_variation)) { + return MappingStatus::NON_FINITE_MAPPING_RESULT; + } + + return MappingStatus::VALID; + } + + MappingStatus DomainMapper::evaluate_face_variation( + const ElementMappingData &element_data, + const ElementDisplacementData &direction, + mfem::FaceElementTransformations &transformation, + FaceElementSide side, + const mfem::IntegrationPoint &ip, + const FaceMappingContext &context, + DMWorkspace &workspace, + FaceMappingVariation &variation + ) const { + transformation.SetAllIntPoints(&ip); + mfem::ElementTransformation &element_transformation = select_face_element_transformation(transformation, side); + const mfem::IntegrationPoint &element_integration_point = select_face_element_integration_point(transformation, side); + + const MappingStatus point_status = evaluate_point_variation(element_data, direction, element_transformation, element_integration_point, context.mapping, workspace, variation.mapping); + if (point_status != MappingStatus::VALID) { + return point_status; + } + + workspace.m_reference_normal.SetSize(m_options.dimension); + mfem::CalcOrtho(transformation.Jacobian(), workspace.m_reference_normal); + if (side == FaceElementSide::ELEMENT_2) { + workspace.m_reference_normal *= -1.0; // Flip the normal for the second element (e.g. the "other" side of the face) + } + + const double reference_normal_magnitude = workspace.m_reference_normal.Norml2(); + if (!is_finite(reference_normal_magnitude) || reference_normal_magnitude <= 0) { + return MappingStatus::NON_FINITE_MAPPING_RESULT; + } + + context.mapping.inverse_mapping_jacobian.MultTranspose(workspace.m_reference_normal, workspace.m_vector_temp); + workspace.m_mapped_normal = workspace.m_vector_temp; + workspace.m_mapped_normal *= context.mapping.mapping_determinate; + + variation.physical_normal_variation.SetSize(m_options.dimension); + variation.mapping.inverse_mapping_jacobian_variation.MultTranspose(workspace.m_reference_normal, variation.physical_normal_variation); + variation.physical_normal_variation *= context.mapping.mapping_determinate; + variation.physical_normal_variation.Add(variation.mapping.mapping_determinate_variation, workspace.m_vector_temp); + + const double mapped_normal_magnitude = workspace.m_mapped_normal.Norml2(); + if (!is_finite(mapped_normal_magnitude) || mapped_normal_magnitude <= 0) { + return MappingStatus::NON_FINITE_MAPPING_RESULT; + } + + const double mapped_normal_magnitude_variation = context.quadrature.normal * variation.physical_normal_variation; + + variation.physical_normal_variation.Add(-mapped_normal_magnitude_variation, context.quadrature.normal); + variation.physical_normal_variation /= mapped_normal_magnitude; + + variation.physical_surface_weight_variation = ip.weight * mapped_normal_magnitude_variation; + variation.normal_flux_scale_variation = mapped_normal_magnitude_variation / reference_normal_magnitude; + + if (!is_finite(variation.physical_normal_variation) || !std::isfinite(variation.physical_surface_weight_variation) || !std::isfinite(variation.normal_flux_scale_variation)) { + return MappingStatus::NON_FINITE_MAPPING_RESULT; + } + + return MappingStatus::VALID; + } + + bool DomainMapper::is_compactified_element(const mfem::ElementTransformation &transformation) const noexcept { + return transformation.Attribute == m_options.vacuum_element_attribute_ID; + } + + int DomainMapper::get_dimension() const noexcept { + return m_options.dimension; + } + + const compactification::ExteriorDomainMap& DomainMapper::get_exterior_map() const noexcept { + return *m_exterior_map; + } + + MappingStatus DomainMapper::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 + ) const { + validate_element_data(element_data); + + //TODO: Factor into its own function + //TODO: Consider if these should return or throw. + if (workspace.get_dimension() != m_options.dimension) { + serif_error("Workspace dimension does not match Domain Mapper dimension. Workspace dimension: {}, Domain Mapper dimension: {}", workspace.get_dimension(), m_options.dimension); + } + if (transformation.GetSpaceDim() != m_options.dimension) { + serif_error("Element transformation space dimension does not match Domain Mapper dimension. Element transformation space dimension: {}, Domain Mapper dimension: {}", transformation.GetSpaceDim(), m_options.dimension); + } + if (transformation.GetGeometryType() != element_data.displacement.get_element().GetGeomType()) { + serif_error("Element transformation geometry type does not match element geometry type. Element transformation geometry type: {}, Element geometry type: {}", transformation.GetGeometryType(), element_data.displacement.get_element().GetGeomType()); + } + + + context.reference_position.SetSize(m_options.dimension); + if (geometry) { + context.reference_position = geometry->position; + } else { + transformation.Transform(ip, context.reference_position); + } + + // Get the displacement field value and its Jacobian at the integration point. Note these may be in current workspace to avoid repeated allocations and Multiplies. + if (displacement_sample) { + // TODO: Generally go through and refactor these large checks into their own functions, at the very least it will be easier to read + if (element_data.displacement.get_element().GetMapType() != mfem::FiniteElement::VALUE || displacement_sample->value.Size() != m_options.dimension || displacement_sample->gradient.Height() != m_options.dimension || displacement_sample->gradient.Width() != m_options.dimension) { + return MappingStatus::INVALID_DIMENSION; + } + + workspace.m_field_value = displacement_sample->value; + workspace.m_field_jacobian.SetSize(m_options.dimension); + mfem::Mult(displacement_sample->gradient, geometry ? geometry->inverse_jacobian : transformation.InverseJacobian(), workspace.m_field_jacobian); + } else { + evaluate_displacement_field(element_data.displacement, transformation, ip, workspace, workspace.m_field_value, workspace.m_field_jacobian, geometry ? &geometry->inverse_jacobian : nullptr); + } + + if (!is_finite(context.reference_position) || !is_finite(workspace.m_field_value) || !is_finite(workspace.m_field_jacobian)) { + return MappingStatus::NON_FINITE_MAPPING_RESULT; + } + + context.displaced_position.SetSize(m_options.dimension); + + // Get the position of the point in physical space by adding the displacement to the reference position. Note MFEM really dislikes raw arithmetic operators + // so we need to first assign the reference position then use the in place += operator. + context.displaced_position = context.reference_position; + context.displaced_position += workspace.m_field_value; + + context.displacement_jacobian.SetSize(m_options.dimension, m_options.dimension); + context.displacement_jacobian = workspace.m_field_jacobian; + + // Ensure that the diagonal of the displacement Jacobian is incremented by 1.0 to account for the identity mapping from reference to physical space. + // recall that r = x + d (where d is the workspace.m_field_value and x is context.reference_position) then we can differentiate this + // component wise to find the gradient of the displaced position wrt. the mesh coordinate (reference position). E.g as you move along + // the mesh coordinate how much does the physical coordinate change and in what direction. Lets call this F + // F = \frac{\partial r_i}{\partial x_j} where r is the displaced position and x is the mesh position. + // We then have F = \frac{\partial x_i}{\partial x_j} + \frac{\partial d_i}{x_j} where d is the displacement (recall r = x + d) + // By definition the first term is the identity matrix. The second term we get out of EvaluateField. Thus why we need to add the identity matrix here + for (int i = 0; i < m_options.dimension; ++i) + context.displacement_jacobian(i, i) += 1.0; + + context.compactified = is_compactified_element(transformation); + + // This branch only runs for vacuum elements + if (context.compactified) { + // There are two things that we need to the mapping. First is a reference coordinate which stroid embeds into the mesh at mesh generation time, this is + // parameterized from 0 - 1 where 0 is the model surface and 1 is the mesh exterior (what will becomes the compactified infinity, note also we never actually evaluate at s=1; rather we define some arbitrary small tolerance to approach s=1). Lets call this s. We also need + // the gradient of s as we move along the mesh coordinates. All of this is stashes within workspace.m_compactification_point. + const MappingStatus coordinate_status = evaluate_compactification_coordinate(element_data.compactification, transformation, ip, workspace, workspace.m_compactification_point, geometry ? &geometry->inverse_jacobian : nullptr); + + if (coordinate_status != MappingStatus::VALID) { + return coordinate_status; + } + + const compactification::ExteriorMapInput exterior_input{ + .reference_position = context.reference_position, + .displaced_position = context.displaced_position, + .displacement_jacobian = context.displacement_jacobian, + .compactification_coordinate_gradient = workspace.m_compactification_point.coordinate_gradient, + .compactification_coordinate = workspace.m_compactification_point.compactification_coordinate, + }; + + // This apply whatever the exterior map is to generate the new physical exterior coordinate and jacobian between physical and reference space. + // In general we have only implemented a kelvin mapping; however, in future additional mappings may be implemented. + const MappingStatus exterior_status = m_exterior_map->evaluate(exterior_input, workspace.m_exterior_result); + if (exterior_status != MappingStatus::VALID) { + return exterior_status; + } + + context.physical_position = workspace.m_exterior_result.physical_position; + context.mapping_jacobian = workspace.m_exterior_result.mapping_jacobian; + } else { + context.physical_position = context.displaced_position; + context.mapping_jacobian = context.displacement_jacobian; + } + + // Validation work + if (!is_finite(context.physical_position) || !is_finite(context.mapping_jacobian)) { + return MappingStatus::NON_FINITE_MAPPING_RESULT; + } + + context.mapping_determinate = context.mapping_jacobian.Det(); + if (!std::isfinite(context.mapping_determinate)) { + return MappingStatus::NON_FINITE_MAPPING_RESULT; + } + if (context.mapping_determinate <= 0.0) { + // This is the most common error we see come out of this function, specifically it is common when we try to deform the mesh too much in one step. + return MappingStatus::NON_POSITIVE_DETERMINATE; + } + + context.inverse_mapping_jacobian.SetSize(m_options.dimension, m_options.dimension); + + // Finally we get the inverse jacobian. + mfem::CalcInverse(context.mapping_jacobian, context.inverse_mapping_jacobian); + + if (!is_finite(context.inverse_mapping_jacobian)) { + return MappingStatus::NON_FINITE_MAPPING_RESULT; + } + + return MappingStatus::VALID; + } + + void DomainMapper::validate_element_data(const ElementMappingData &element_data) const { + const ElementDisplacementData &displacement = element_data.displacement; + const ElementCompactificationData &compactification = element_data.compactification; + + if (displacement.get_dimension() != m_options.dimension) { + serif_error("Element displacement dimension does not match Domain Mapper dimension. Element displacement dimension: {}, Domain Mapper dimension: {}", displacement.get_dimension(), m_options.dimension); + } + + if (displacement.get_element().GetDim() != m_options.dimension) { + serif_error("Element displacement finite element dimension does not match Domain Mapper dimension. Element displacement finite element dimension: {}, Domain Mapper dimension: {}", displacement.get_element().GetDim(), m_options.dimension); + } + + if (compactification.get_element().GetDim() != m_options.dimension) { + serif_error("Element compactification finite element dimension does not match Domain Mapper dimension. Element compactification finite element dimension: {}, Domain Mapper dimension: {}", compactification.get_element().GetDim(), m_options.dimension); + } + + if (displacement.get_element().GetGeomType() != compactification.get_element().GetGeomType()) { + serif_error("Element displacement finite element geometry type does not match element compactification finite element geometry type. Element displacement finite element geometry type: {}, Element compactification finite element geometry type: {}", displacement.get_element().GetGeomType(), compactification.get_element().GetGeomType()); + } + + if (compactification.get_element().GetRangeType() != mfem::FiniteElement::SCALAR) { + serif_error("Element compactification finite element range type is not scalar. Element compactification finite element range type: {}", compactification.get_element().GetRangeType()); + } + + if (compactification.get_element().GetMapType() != mfem::FiniteElement::VALUE) { + serif_error("Element compactification finite element map type is not value. Element compactification finite element map type: {}", compactification.get_element().GetMapType()); + } + + if (compactification.get_element().GetDerivType() != mfem::FiniteElement::GRAD) { + serif_error("Element compactification finite element derivative type is not gradient. Element compactification finite element derivative type: {}", compactification.get_element().GetDerivType()); + } + + if (compactification.get_dof_count() != compactification.get_element().GetDof()) { + serif_error("Element compactification finite element DOF count does not match the number of DOFs in the finite element. Element compactification finite element DOF count: {}, Finite element DOF count: {}", compactification.get_dof_count(), compactification.get_element().GetDof()); + } + } + + // Note that there is nothing preventing this function from being used with some field other than displacement; perhaps we should tighten this? + void DomainMapper::evaluate_displacement_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 { + const mfem::FiniteElement &element = field.get_element(); + const mfem::DenseMatrix &dof_matrix = field.get_dof_matrix(); + + workspace.m_shape.SetSize(element.GetDof()); + + element.CalcShape(ip, workspace.m_shape); + value.SetSize(m_options.dimension); + dof_matrix.MultTranspose(workspace.m_shape, value); + jacobian.SetSize(m_options.dimension, m_options.dimension); + + if (inverse_mesh_jacobian != nullptr || element.GetMapType() == mfem::FiniteElement::VALUE) { + workspace.m_reference_dshape.SetSize(element.GetDof(), m_options.dimension); + element.CalcDShape(ip, workspace.m_reference_dshape); + // Contract DOFs before applying fixed-mesh geometry. + mfem::MultAtB(dof_matrix, workspace.m_reference_dshape, workspace.m_reference_field_jacobian); + const mfem::DenseMatrix &inverseMeshJacobian = inverse_mesh_jacobian != nullptr ? *inverse_mesh_jacobian : transformation.InverseJacobian(); + mfem::Mult(workspace.m_reference_field_jacobian, inverseMeshJacobian, jacobian); + } else { + // Retain the original finite-element-specific physical-gradient + // path for mapping types without the ordinary VALUE pullback. + workspace.m_mesh_dshape.SetSize(element.GetDof(), m_options.dimension); + element.CalcPhysDShape(transformation, workspace.m_mesh_dshape); + mfem::MultAtB(dof_matrix, workspace.m_mesh_dshape, jacobian); + } + } + + MappingStatus DomainMapper::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 { + const mfem::FiniteElement& element = compactification.get_element(); + const mfem::Vector& dofs = compactification.get_dofs(); + const int dof_count = element.GetDof(); + + if (workspace.get_dimension() != m_options.dimension || transformation.GetSpaceDim() != m_options.dimension || element.GetDim() != m_options.dimension) { + return MappingStatus::INVALID_DIMENSION; + } + + if (dofs.Size() != dof_count) { + return MappingStatus::INVALID_DIMENSION; + } + + for (int i = 0; i < dof_count; ++i) { + if (!is_finite(dofs(i))) { + return MappingStatus::NON_FINITE_MAPPING_INPUT; + } + } + + workspace.m_compactification_shape.SetSize(dof_count); + workspace.m_compactification_dshape.SetSize(dof_count, m_options.dimension); + + element.CalcShape(ip, workspace.m_compactification_shape); + if (inverse_mesh_jacobian != nullptr) { + workspace.m_reference_dshape.SetSize(dof_count, m_options.dimension); + element.CalcDShape(ip, workspace.m_reference_dshape); + mfem::Mult(workspace.m_reference_dshape, *inverse_mesh_jacobian, workspace.m_compactification_dshape); + } else { + element.CalcPhysDShape(transformation, workspace.m_compactification_dshape); + } + + point_data.compactification_coordinate = dofs * workspace.m_compactification_shape; + point_data.coordinate_gradient.SetSize(m_options.dimension); + workspace.m_compactification_dshape.MultTranspose(dofs, point_data.coordinate_gradient); + + if (!is_finite(point_data.compactification_coordinate)) { + return MappingStatus::NON_FINITE_MAPPING_RESULT; + } + + for (int d = 0; d < point_data.coordinate_gradient.Size(); ++d) { + if (!is_finite(point_data.coordinate_gradient(d))) { + return MappingStatus::NON_FINITE_MAPPING_RESULT; + } + } + + return MappingStatus::VALID; + } + + MappingStatus DomainMapper::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 { + validate_element_data(element_data); + const ElementMappingData direction_data { + .displacement = direction, + .compactification = element_data.compactification + }; + + validate_element_data(direction_data); + + if (element_data.displacement.get_dof_count() != direction.get_dof_count()) { + serif_error("Element displacement DOF count does not match direction DOF count. Element displacement DOF count: {}, Direction DOF count: {}", element_data.displacement.get_dof_count(), direction.get_dof_count()); + } + + if (workspace.get_dimension() != m_options.dimension) { + serif_error("Workspace dimension does not match Domain Mapper dimension. Workspace dimension: {}, Domain Mapper dimension: {}", workspace.get_dimension(), m_options.dimension); + } + + if (context.compactified != is_compactified_element(transformation)) { + serif_error("Mapping context compactification status does not match element transformation compactification status. Mapping context compactified: {}, Element transformation compactified: {}", context.compactified, is_compactified_element(transformation)); + } + + if (inverse_mesh_jacobian == nullptr) { + if (const auto* geometry = get_fixed_mesh_geometry(transformation, ip)) { + inverse_mesh_jacobian = &geometry->inverse_jacobian; + } + } + evaluate_displacement_field(direction, transformation, ip, workspace, workspace.m_field_value, workspace.m_field_jacobian, inverse_mesh_jacobian); + + if (!is_finite(workspace.m_field_jacobian) || !is_finite(workspace.m_field_value)) { + return MappingStatus::NON_FINITE_MAPPING_INPUT; + } + + variation.displacement_variation = workspace.m_field_value; + variation.displacement_jacobian_variation = workspace.m_field_jacobian; + + if (context.compactified) { + const MappingStatus coordinate_status = evaluate_compactification_coordinate(element_data.compactification, transformation, ip, workspace, workspace.m_compactification_point, inverse_mesh_jacobian); + if (coordinate_status != MappingStatus::VALID) { + return coordinate_status; + } + + const compactification::ExteriorMapInput exterior_input{ + .reference_position = context.reference_position, + .displaced_position = context.displaced_position, + .displacement_jacobian = context.displacement_jacobian, + .compactification_coordinate_gradient = workspace.m_compactification_point.coordinate_gradient, + .compactification_coordinate = workspace.m_compactification_point.compactification_coordinate, + }; + + workspace.m_exterior_result.physical_position = context.physical_position; + workspace.m_exterior_result.mapping_jacobian = context.mapping_jacobian; + + const compactification::ExteriorMapDirection exterior_direction { + .displaced_position_variation = variation.displacement_variation, + .displacement_jacobian_variation = variation.displacement_jacobian_variation + }; + + const MappingStatus exterior_status = m_exterior_map->evaluate_variation(exterior_input, workspace.m_exterior_result, exterior_direction, workspace.m_exterior_variation); + if (exterior_status != MappingStatus::VALID) { + return exterior_status; + } + + variation.physical_position_variation = workspace.m_exterior_variation.physical_position_variation; + variation.mapping_jacobian_variation = workspace.m_exterior_variation.mapping_jacobian_variation; + } else { + variation.physical_position_variation = variation.displacement_variation; + variation.mapping_jacobian_variation = variation.displacement_jacobian_variation; + } + + mfem::Mult(context.inverse_mapping_jacobian, variation.mapping_jacobian_variation, workspace.m_matrix_temp_1); + double trace = 0.0; + for (int dimID = 0; dimID < m_options.dimension; ++dimID) { + trace += workspace.m_matrix_temp_1(dimID, dimID); + } + variation.mapping_determinate_variation = context.mapping_determinate * trace; + + variation.inverse_mapping_jacobian_variation.SetSize(m_options.dimension, m_options.dimension); + mfem::Mult(workspace.m_matrix_temp_1, context.inverse_mapping_jacobian, variation.inverse_mapping_jacobian_variation); + variation.inverse_mapping_jacobian_variation *= -1.0; + + if (!is_finite(variation.physical_position_variation) || !is_finite(variation.mapping_jacobian_variation) || !std::isfinite(variation.mapping_determinate_variation) || !is_finite(variation.inverse_mapping_jacobian_variation)) { + return MappingStatus::NON_FINITE_MAPPING_RESULT; + } + + return MappingStatus::VALID; + } + + mfem::ElementTransformation &DomainMapper::select_face_element_transformation( + const mfem::FaceElementTransformations &transformations, + const FaceElementSide side + ) { + if (side == FaceElementSide::ELEMENT_1) { + serif_assert(transformations.Elem1 != nullptr, "Face element transformation for ELEMENT_1 is null."); + return *transformations.Elem1; + } + + serif_assert(transformations.Elem2 != nullptr, "Face element transformation for ELEMENT_2 is null."); + return *transformations.Elem2; + } + + const mfem::IntegrationPoint &DomainMapper::select_face_element_integration_point( + const mfem::FaceElementTransformations &transformations, + const FaceElementSide side + ) { + mfem::ElementTransformation& element_transformation = select_face_element_transformation(transformations, side); + return element_transformation.GetIntPoint(); } } diff --git a/src/lib/serif/mapping/evaluator.cpp b/src/lib/serif/mapping/evaluator.cpp new file mode 100644 index 00000000..8bd6142c --- /dev/null +++ b/src/lib/serif/mapping/evaluator.cpp @@ -0,0 +1,149 @@ +#include "serif/mapping/evaluator.hpp" + +#include "serif/utils/error/errors.hpp" + + +namespace serif::mapping { + using namespace utils::errors; + + GridFunctionMappingEvaluator::GridFunctionMappingEvaluator( + const DomainMapper &mapper, + const mfem::GridFunction &displacement, + const mfem::GridFunction &compactification_coordinate + ) : + m_mapper(mapper), + m_displacement(displacement), + m_compactification_coordinate(compactification_coordinate), + m_displacement_space(displacement.FESpace()), + m_compactification_space(compactification_coordinate.FESpace()), + m_displacement_space_sequence(m_displacement_space != nullptr ? m_displacement_space->GetSequence() : -1), + m_compactification_space_sequence(m_compactification_space != nullptr ? m_compactification_space->GetSequence() : -1), + m_workspace(mapper.get_dimension()) { + serif_assert(m_displacement_space != nullptr, "Displacement grid function is not associated with a finite-element space."); + serif_assert(m_compactification_space != nullptr, "Compactification-coordinate grid function is not associated with a finite-element space."); + serif_assert(m_displacement_space->GetMesh() == m_compactification_space->GetMesh(), "Displacement and compactification-coordinate grid functions are not associated with the same mesh."); + serif_assert(m_displacement.VectorDim() == mapper.get_dimension(), "Displacement grid function vector dimension does not match the mapper dimension. ({} vs. {})", m_displacement.VectorDim(), mapper.get_dimension()); + serif_assert(m_compactification_coordinate.VectorDim() == 1, "Compactification-coordinate grid function must be a scalar field. Instead it has vector dimension {}.", m_compactification_coordinate.VectorDim()); + serif_assert(m_displacement_space->GetMesh()->SpaceDimension() == mapper.get_dimension(), "Displacement finite-element space mesh dimension does not match the mapper dimension. ({} vs. {})", m_displacement_space->GetMesh()->SpaceDimension(), mapper.get_dimension()); + } + + void GridFunctionMappingEvaluator::invalidate_cache() noexcept { + m_displacement_data.reset(); + m_compactification_data.reset(); + m_cached_element_id = -1; + } + + void GridFunctionMappingEvaluator::validate_field_bindings() const { + serif_assert(m_displacement.FESpace() == m_displacement_space, "Displacement grid function is no longer associated with the cached finite-element space."); + serif_assert(m_compactification_coordinate.FESpace() == m_compactification_space, "Compactification-coordinate grid function is no longer associated with the cached finite-element space."); + } + + bool GridFunctionMappingEvaluator::invalidate_for_changed_spaces() { + const long displacement_sequence = m_displacement_space->GetSequence(); + const long compactification_sequence = m_compactification_space->GetSequence(); + + if (displacement_sequence == m_displacement_space_sequence && compactification_sequence == m_compactification_space_sequence) { + return false; + } + + invalidate_cache(); + + m_displacement_space_sequence = displacement_sequence; + m_compactification_space_sequence = compactification_sequence; + return true; + } + + void GridFunctionMappingEvaluator::refresh() { + validate_field_bindings(); + + if (invalidate_for_changed_spaces()) { + return; + } + + const int element_id = m_cached_element_id; + invalidate_cache(); + if (element_id >= 0) { + load_element(element_id); + } + } + + void GridFunctionMappingEvaluator::load_element(const int element_id) { + validate_field_bindings(); + (void)invalidate_for_changed_spaces(); + + if (element_id == m_cached_element_id) { + return; + } + + const mfem::FiniteElementSpace& displacement_space = *m_displacement_space; + const mfem::FiniteElementSpace& compactification_space = *m_compactification_space; + + serif_ensure(element_id >= 0, "Element ID must be non-negative. Instead {} was provided.", element_id); + serif_ensure(element_id < displacement_space.GetMesh()->GetNE(), "Element ID must be less than the number of elements in the mesh ({}). Instead {} was provided.", displacement_space.GetMesh()->GetNE(), element_id); + + mfem::DofTransformation *displacement_transformation = displacement_space.GetElementVDofs(element_id, m_displacement_dofs); + compactification_space.GetElementDofs(element_id, m_compactification_dofs); + + m_displacement.GetSubVector(m_displacement_dofs, m_element_displacement); + m_compactification_coordinate.GetSubVector(m_compactification_dofs, m_element_compactification); + + if (displacement_transformation != nullptr) { + displacement_transformation->InvTransformPrimal(m_element_displacement); + } + + const mfem::FiniteElement &displacement_element = *displacement_space.GetFE(element_id); + const mfem::FiniteElement &compactification_element = *compactification_space.GetFE(element_id); + + m_displacement_data = std::make_unique(element_displacement_data_from_element_vdofs(displacement_element, m_element_displacement)); + m_compactification_data = std::make_unique(compactification_element, m_element_compactification); + m_cached_element_id = element_id; + } + + MappingStatus GridFunctionMappingEvaluator::evaluate_point(mfem::ElementTransformation &transformation, const mfem::IntegrationPoint &integration_point, MappingPointContext &context) { + load_element(transformation.ElementNo); + const ElementMappingData data{ + .displacement = *m_displacement_data, + .compactification = *m_compactification_data + }; + return m_mapper.evaluate_point(data, transformation, integration_point, m_workspace, context); + } + + MappingStatus GridFunctionMappingEvaluator::evaluate_volume(mfem::ElementTransformation &transformation, const mfem::IntegrationPoint &integration_point, VolumeMappingContext &context) { + load_element(transformation.ElementNo); + const ElementMappingData data{ + .displacement = *m_displacement_data, + .compactification = *m_compactification_data + }; + return m_mapper.evaluate_volume(data, transformation, integration_point, m_workspace, context); + } + + MappingStatus GridFunctionMappingEvaluator::evaluate_face(mfem::FaceElementTransformations &transformations, FaceElementSide side, const mfem::IntegrationPoint &integration_point, FaceMappingContext &context) { + mfem::ElementTransformation *transformation = side == FaceElementSide::ELEMENT_1 ? transformations.Elem1 : transformations.Elem2; + serif_ensure(transformation != nullptr, "Face evaluation requested for a side that does not exist."); + + load_element(transformation->ElementNo); + const ElementMappingData data{ + .displacement = *m_displacement_data, + .compactification = *m_compactification_data + }; + return m_mapper.evaluate_face(data, transformations, side, integration_point, m_workspace, context); + } + + VolumeQuadratureContext GridFunctionMappingEvaluator::get_quadrature_context(mfem::ElementTransformation &transformation, const mfem::IntegrationPoint &integration_point) { + VolumeMappingContext context; + serif_ensure(evaluate_volume(transformation, integration_point, context) == MappingStatus::VALID, "Failed to evaluate the mapping at the given volume quadrature point."); + return context.quadrature; + } + + FaceQuadratureContext GridFunctionMappingEvaluator::get_face_quadrature_context(mfem::FaceElementTransformations &transformation, const mfem::IntegrationPoint &integration_point, FaceElementSide side) { + FaceMappingContext context; + serif_ensure(evaluate_face(transformation, side, integration_point, context) == MappingStatus::VALID, "Failed to evaluate the mapping at the given face quadrature point."); + return context.quadrature; + } + + void GridFunctionMappingEvaluator::get_physical_point(mfem::ElementTransformation &transformation, const mfem::IntegrationPoint &integration_point, mfem::Vector &physical_position) { + MappingPointContext context; + serif_ensure(evaluate_point(transformation, integration_point, context) == MappingStatus::VALID, "Failed to evaluate the mapping at the given point."); + physical_position = context.physical_position; + } +} diff --git a/src/lib/serif/mapping/transformations.cpp b/src/lib/serif/mapping/transformations.cpp new file mode 100644 index 00000000..62afa4a2 --- /dev/null +++ b/src/lib/serif/mapping/transformations.cpp @@ -0,0 +1,136 @@ +#include "serif/mapping/transformations.hpp" +#include "serif/utils/error/errors.hpp" + +#include "mfem/linalg/vector.hpp" +#include "mfem/linalg/densemat.hpp" +#include "serif/utils/checks/finite.hpp" + +namespace serif::mapping { + using namespace serif::utils::errors; + void map_hdiv_flux_to_physical(const MappingPointContext &context, const mfem::Vector &reference_flux, mfem::Vector &physical_flux) { + serif_assert(reference_flux.Size() == context.mapping_jacobian.Width(), "Reference flux vector has incompatible dimension with the mapping jacobian ({} <=> {}).", reference_flux.Size(), context.mapping_jacobian.Width()); + + physical_flux.SetSize(reference_flux.Size()); + context.mapping_jacobian.Mult(reference_flux, physical_flux); + physical_flux /= context.mapping_determinate; + } + + void map_physical_flux_to_hdiv_reference(const MappingPointContext &context, const mfem::Vector &physical_flux, mfem::Vector &reference_flux) { + serif_assert(physical_flux.Size() == context.inverse_mapping_jacobian.Height(), "Physical flux vector has incompatible dimension with the inverse mapping jacobian ({} <=> {}).", physical_flux.Size(), context.inverse_mapping_jacobian.Height()); + + reference_flux.SetSize(physical_flux.Size()); + context.inverse_mapping_jacobian.Mult(physical_flux, reference_flux); + reference_flux /= context.mapping_determinate; + } + + void map_reference_gradient_to_physical(const MappingPointContext &context, const mfem::Vector &reference_gradient, mfem::Vector &physical_gradient) { + serif_assert(reference_gradient.Size() == context.inverse_mapping_jacobian.Height(), "Reference gradient vector has incompatible dimension with the inverse mapping jacobian ({} <=> {}).", reference_gradient.Size(), context.inverse_mapping_jacobian.Height()); + + physical_gradient.SetSize(reference_gradient.Size()); + context.inverse_mapping_jacobian.Mult(reference_gradient, physical_gradient); + } + + void map_physical_gradient_to_reference(const MappingPointContext &context, const mfem::Vector &physical_gradient, mfem::Vector &reference_gradient) { + serif_assert(physical_gradient.Size() == context.mapping_jacobian.Height(), "Physical gradient vector has incompatible dimension with the mapping jacobian ({} <=> {}).", physical_gradient.Size(), context.mapping_jacobian.Height()); + + reference_gradient.SetSize(physical_gradient.Size()); + context.mapping_jacobian.MultTranspose(physical_gradient, reference_gradient); + } + + void map_reference_vector_gradient_to_physical(const MappingPointContext &context, const mfem::DenseMatrix &reference_vector_gradient, mfem::DenseMatrix &physical_vector_gradient) { + serif_assert(reference_vector_gradient.Width() == context.inverse_mapping_jacobian.Height(), "Reference vector gradient matrix has incompatible dimension with the inverse mapping jacobian ({} <=> {}).", reference_vector_gradient.Width(), context.inverse_mapping_jacobian.Height()); + + physical_vector_gradient.SetSize(reference_vector_gradient.Height(), context.inverse_mapping_jacobian.Width()); + mfem::Mult(reference_vector_gradient, context.inverse_mapping_jacobian, physical_vector_gradient); + } + + void map_physical_vector_gradient_to_reference(const MappingPointContext &context, const mfem::DenseMatrix &physical_vector_gradient, mfem::DenseMatrix &reference_vector_gradient) { + serif_assert(physical_vector_gradient.Width() == context.mapping_jacobian.Height(), "Physical vector gradient matrix has incompatible dimension with the mapping jacobian ({} <=> {}).", physical_vector_gradient.Width(), context.mapping_jacobian.Height()); + + reference_vector_gradient.SetSize(physical_vector_gradient.Height(), context.mapping_jacobian.Width()); + mfem::Mult(physical_vector_gradient, context.mapping_jacobian, reference_vector_gradient); + } + + double map_hdiv_divergence_to_physical(const MappingPointContext &context, double reference_divergence) { + return reference_divergence / context.mapping_determinate; + } + + void compute_hdiv_mass_tensor(const MappingPointContext &context, mfem::DenseMatrix &mass_tensor) { + const int dimension = context.mapping_jacobian.Height(); + serif_assert((context.mapping_jacobian.Width() == dimension) && (context.mapping_jacobian.Height() == dimension), "Mapping jacobian is not square ({} x {}). It is expected to have shape {}x{}", context.mapping_jacobian.Height(), context.mapping_jacobian.Width(), dimension, dimension); + serif_assert(context.mapping_determinate > 0.0, "Mapping jacobian determinant is non-positive ({}). The mapping may be inverted or degenerate.", context.mapping_determinate); + + mass_tensor.SetSize(dimension, dimension); + mfem::MultAtB(context.mapping_jacobian, context.mapping_jacobian, mass_tensor); + mass_tensor *= 1/context.mapping_determinate; + } + + void compute_scalar_diffusion_tensor(const MappingPointContext &context, mfem::DenseMatrix &diffusion_tensor) { + const int dimension = context.mapping_jacobian.Height(); + serif_assert((context.inverse_mapping_jacobian.Width() == dimension) && (context.inverse_mapping_jacobian.Height() == dimension), "Inverse mapping jacobian is not square ({} x {}). It is expected to have shape {}x{}", context.inverse_mapping_jacobian.Height(), context.inverse_mapping_jacobian.Width(), dimension, dimension); + serif_assert(context.mapping_determinate > 0.0, "Mapping jacobian determinant is non-positive ({}). The mapping may be inverted or degenerate.", context.mapping_determinate); + + diffusion_tensor.SetSize(dimension, dimension); + mfem::MultAtB(context.inverse_mapping_jacobian, context.inverse_mapping_jacobian, diffusion_tensor); + diffusion_tensor *= context.mapping_determinate; + } + + void map_hcurl_field_to_physical(const MappingPointContext &context, const mfem::Vector &reference_field, mfem::Vector &physical_field) { + serif_assert(reference_field.Size() == context.inverse_mapping_jacobian.Height(), "Reference field vector has incompatible dimension with the inverse mapping jacobian ({} <=> {}).", reference_field.Size(), context.inverse_mapping_jacobian.Height()); + + physical_field.SetSize(reference_field.Size()); + context.inverse_mapping_jacobian.MultTranspose(reference_field, physical_field); + } + + void map_physical_field_to_hcurl_reference(const MappingPointContext &context, const mfem::Vector &physical_field, mfem::Vector &reference_field) { + serif_assert(physical_field.Size() == context.mapping_jacobian.Height(), "Physical field vector has incompatible dimension with the mapping jacobian ({} <=> {}).", physical_field.Size(), context.mapping_jacobian.Height()); + + reference_field.SetSize(physical_field.Size()); + context.mapping_jacobian.MultTranspose(physical_field, reference_field); + } + + void map_hcurl_curl_to_physical(const MappingPointContext &context, const mfem::Vector &reference_curl, mfem::Vector &physical_curl) { + serif_assert(reference_curl.Size() == context.mapping_jacobian.Width(), "Reference curl vector has incompatible dimension with the mapping jacobian ({} <=> {}).", reference_curl.Size(), context.mapping_jacobian.Height()); + + physical_curl.SetSize(reference_curl.Size()); + context.mapping_jacobian.Mult(reference_curl, physical_curl); + physical_curl /= context.mapping_determinate; + } + + void map_physical_curl_to_hcurl_reference(const MappingPointContext &context, const mfem::Vector &physical_curl, mfem::Vector &reference_curl) { + serif_assert(physical_curl.Size() == context.inverse_mapping_jacobian.Width(), "Physical curl vector has incompatible dimension with the inverse mapping jacobian ({} <=> {}).", physical_curl.Size(), context.inverse_mapping_jacobian.Height()); + + reference_curl.SetSize(physical_curl.Size()); + context.inverse_mapping_jacobian.Mult(physical_curl, reference_curl); + reference_curl *= context.mapping_determinate; + } + + void compute_hdiv_mass_tensor_variation(const MappingPointContext &context, const MappingPointVariation &variation, mfem::DenseMatrix &mass_tensor_variation) { + const double determinate = context.mapping_determinate; + const double determinate_variation = variation.mapping_determinate_variation; + const int dimension = context.inverse_mapping_jacobian.Width(); + + mass_tensor_variation.SetSize(dimension, dimension); + + serif_assert(determinate > 0.0, "Mapping jacobian determinant is non-positive ({}). The mapping may be inverted or degenerate.", determinate); + serif_assert(utils::checks::is_finite(determinate), "Mapping jacobian determinant is non-finite ({}). The mapping may be inverted or degenerate.", determinate); + + const mfem::DenseMatrix& jacobian = context.mapping_jacobian; + const mfem::DenseMatrix& jacobian_variation = variation.mapping_jacobian_variation; + const double inverse_determinate = 1.0 / determinate; + const double determinate_scale = determinate_variation * inverse_determinate; + + for (int row = 0; row < dimension; ++row) { + for (int col = 0; col < dimension; ++col) { + double gram{0.0}; + double gram_variation{0.0}; + + for (int inner = 0; inner < dimension; ++inner) { + gram += jacobian(inner, row) * jacobian(inner, col); + gram_variation += jacobian(inner, row) * jacobian_variation(inner, col) + jacobian_variation(inner, row) * jacobian(inner, col); + } + mass_tensor_variation(row, col) = inverse_determinate * (gram_variation - determinate_scale * gram); + } + } + } +} diff --git a/src/lib/serif/mapping/volume_cache.cpp b/src/lib/serif/mapping/volume_cache.cpp new file mode 100644 index 00000000..15fca1d4 --- /dev/null +++ b/src/lib/serif/mapping/volume_cache.cpp @@ -0,0 +1,112 @@ +#include "serif/mapping/volume_cache.hpp" +#include "serif/utils/error/errors.hpp" + +#include "mfem/linalg/densemat.hpp" +#include "mfem/linalg/vector.hpp" + +namespace serif::mapping { + using namespace utils::errors; + namespace { + // Replace all of these with a fast serialization class at some point. Probably in something like a serif/serialize module + void pack_vector(double*& destination, const mfem::Vector& vector) { + const int dim = vector.Size(); + std::copy_n(vector.HostRead(), dim, destination); + destination += dim; + } + + void pack_matrix(double*& destination, const mfem::DenseMatrix& matrix) { + serif_assert(matrix.Height() == matrix.Width(), "Matrix must be square to pack"); + const int dim = matrix.Height(); + std::copy_n(matrix.HostRead(), dim * dim, destination); + destination += dim * dim; + } + + void unpack_vector(const double*& source, mfem::Vector& vector, const int dim) { + vector.SetSize(dim); + std::copy_n(source, dim, vector.HostWrite()); + source += dim; + } + + void unpack_matrix(const double*& source, mfem::DenseMatrix& matrix, const int dim) { + matrix.SetSize(dim, dim); + std::copy_n(source, dim * dim, matrix.HostWrite()); + source += dim * dim; + } + } + + void VolumeMappingCache::set_size(int point_count, int dimension) { + serif_assert(point_count < 0 || dimension < 1 || dimension > 3, "Invalid point count or dimension for volume mapping cache. Point count: {}, dimension: {}", point_count, dimension); + const int stride = 3 * dimension + 4 * dimension * dimension + 4; + m_data.resize(static_cast(point_count) * stride); + m_point_count = point_count; + m_dimension = dimension; + m_point_stride = stride; + } + + void VolumeMappingCache::validate_point(const int point) const { + serif_assert(point < 0 || point >= m_point_count, "Point index {} is out of range for volume mapping cache with {} points", point, m_point_count); + } + + const double *VolumeMappingCache::get_point_data(const int point) const { + validate_point(point); + return m_data.data() + static_cast(point) * m_point_stride; + } + + void VolumeMappingCache::store(int point, VolumeMappingContext &context) { + validate_point(point); + double *data = m_data.data() + static_cast(point) * m_point_stride; + + // What we build here is a FIFO data structure. we pack + // all the raw data from the mapping context into a single contiguous array of doubles. + // This allows us to avoid repeated allocations as the vector will grow to some appropriate size and then stay roughly there. + pack_vector(data, context.mapping.reference_position); + pack_vector(data, context.mapping.displaced_position); + pack_vector(data, context.mapping.physical_position); + + pack_matrix(data, context.mapping.displacement_jacobian); + pack_matrix(data, context.mapping.mapping_jacobian); + pack_matrix(data, context.mapping.inverse_mapping_jacobian); + pack_matrix(data, context.quadrature.J_element_to_physical_inverse); + + *data++ = context.mapping.mapping_determinate; + *data++ = context.mapping.compactified ? 1.0 : 0.0; // its possible this could cause issues, have not run into but when unpacked I worry that exact equality to some double value could be a problem. If there are wierd issues this may be a place to look, + *data++ = context.quadrature.mapping_determinate; + *data = context.quadrature.physical_integration_weight; + } + + void VolumeMappingCache::load(int point, VolumeMappingContext &context) const { + const double *data = get_point_data(point); + + // FIFO data structure so need to unpack in the same order as we packed. + // TODO: Redo this with some sort of fast serialization class so we do not need to worry about ordering + unpack_vector(data, context.mapping.reference_position, m_dimension); + unpack_vector(data, context.mapping.displaced_position, m_dimension); + unpack_vector(data, context.mapping.physical_position, m_dimension); + + unpack_matrix(data, context.mapping.displacement_jacobian, m_dimension); + unpack_matrix(data, context.mapping.mapping_jacobian, m_dimension); + unpack_matrix(data, context.mapping.inverse_mapping_jacobian, m_dimension); + unpack_matrix(data, context.quadrature.J_element_to_physical_inverse, m_dimension); + + context.mapping.mapping_determinate = *data++; + context.mapping.compactified = *data++ != 0.0; // See comment in store() + context.quadrature.mapping_determinate = *data++; + context.quadrature.physical_integration_weight = *data; + } + + void VolumeMappingCache::load_inverse_jacobian(int point, mfem::DenseMatrix &inverse) const { + // Again this is unsafe, and really should be replaced. This was a quick cache scheme that solved the direct issue I had; however, it should be replaced + // with something that has at least a modicum of safety built in (we are pretty good about this in other parts of SERiF but its about 8pm and + // Im quite tired.... + const double *data = get_point_data(point) + 3 * m_dimension + 3 * m_dimension * m_dimension; // skip to the inverse jacobian, this is just pointer math...Think back to CS-101. Do that have that in the future, are we still learning about memory like this is a PDP-11? + unpack_matrix(data, inverse, m_dimension); + } + + int VolumeMappingCache::get_point_count() const noexcept { + return m_point_count; + } + + int VolumeMappingCache::get_dimension() const noexcept { + return m_dimension; + } +} diff --git a/src/meson.build b/src/meson.build index 8a13bceb..8323c4b8 100644 --- a/src/meson.build +++ b/src/meson.build @@ -1,5 +1,8 @@ include_dir = include_directories('include') +# TODO: Really we should make a dependency graph between the various modules +# and then make sublibraries for each module that build themselves in the correct order +# and then link together at the end. serif_sources = files( 'lib/serif/discretization/domain/mesh/topology.cpp', 'lib/serif/discretization/domain/schema/validation/results.cpp', @@ -8,6 +11,10 @@ serif_sources = files( 'lib/serif/mapping/compactification/kelvin.cpp', 'lib/serif/mapping/data.cpp', 'lib/serif/mapping/domain_mapper.cpp', + 'lib/serif/mapping/transformations.cpp', + 'lib/serif/mapping/volume_cache.cpp', + 'lib/serif/mapping/evaluator.cpp', + 'lib/serif/mapping/coefficients.cpp' ) serif_lib = library( diff --git a/tests/sandbox/mapping_sandbox.cpp b/tests/sandbox/mapping_sandbox.cpp index 9ea7816a..4c64f270 100644 --- a/tests/sandbox/mapping_sandbox.cpp +++ b/tests/sandbox/mapping_sandbox.cpp @@ -6,6 +6,7 @@ #include "serif/mapping/status.hpp" #include "serif/mapping/compactification/kelvin.hpp" #include "serif/mapping/domain_mapper.hpp" +#include "serif/mapping/transformations.hpp" #include "serif/utils/user_options/mapping/compactification/kelvin.hpp" int main() {