feat(mapper): finished porting over mapper

src/include/serif/mapping/coefficients.hpp

@@ -0,0 +1,34 @@
+#pragma once
+
+#include <functional>
+
+#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<double(const mfem::Vector& x)>;
+
+        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;
+    };
+}

src/include/serif/mapping/compactification/compactification.hpp

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <string_view>
+#include <memory>
 
 #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<const ExteriorDomainMap> clone() const = 0;
     };
 }

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<const ExteriorDomainMap> clone() const override;
 
         [[nodiscard]] double get_reference_stellar_radius() const noexcept;
         [[nodiscard]] double get_reference_infinity_radius() const noexcept;

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 {

src/include/serif/mapping/data.hpp

@@ -6,7 +6,7 @@
 #include "mfem/mesh/mesh.hpp"
 
 #include <cstdint>
-#include <unordered_map>
+#include <map>
 #include <tuple>
 
 #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<std::tuple<int, double, double, double>, DMFixedMeshGeometry> points;
+        std::map<std::tuple<int, double, double, double>, 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<int, double, double, double>. 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};
     };

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<compactification::ExteriorDomainMap> 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<const compactification::ExteriorDomainMap> exterior_map
+        );
+
+        std::unique_ptr<DomainMapper> 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
         );
 

src/include/serif/mapping/evaluator.hpp

@@ -0,0 +1,93 @@
+#pragma once
+#include <memory>
+
+#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<int> m_displacement_dofs;
+        mfem::Array<int> m_compactification_dofs;
+        mfem::Vector m_element_displacement;
+        mfem::Vector m_element_compactification;
+        std::unique_ptr<ElementDisplacementData> m_displacement_data;
+        std::unique_ptr<ElementCompactificationData> m_compactification_data;
+        int m_cached_element_id{-1};
+    };
+}

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
+    );
+}

src/include/serif/mapping/volume_cache.hpp

@@ -0,0 +1,28 @@
+#pragma once
+
+#include "mfem/linalg/densemat.hpp"
+
+#include "serif/mapping/context.hpp"
+#include <vector>
+
+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<double> m_data;
+        int m_point_count{0};
+        int m_dimension{0};
+        int m_point_stride{0};
+    };
+
+}

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 <typename Scalar>
+        requires requires (Scalar x) {{std::isfinite(x)} -> std::convertible_to<bool>; }
+    bool is_finite(const Scalar scalar) {
+        return std::isfinite(scalar);
+    }
 }

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

src/include/serif/utils/error/errors.hpp

@@ -38,7 +38,7 @@ namespace serif::utils::errors {
         throw ExceptionT(error_msg);
     }
 
-    template <SERiFErrorCode code>
+    template <SERiFErrorCode code, misc::concepts::IsException ExceptionT = std::logic_error>
     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<std::logic_error, code>(message, args...);
+                serif_error<ExceptionT, code>(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 <SERiFErrorCode code, misc::concepts::IsException ExceptionT = std::logic_error>
+    void serif_ensure(const bool condition, const std::string_view message, auto... args) {
+        if (!condition) {
+            serif_error<ExceptionT, code>(message, args...);
+        }
+    }
+
+    template <SERiFErrorCode code, types::ct_error_string message, bool condition>
+    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<message, code>);
+        }
+    }
 }

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<SERiFErrorCode::INVALID_MAPPING>(m_mapping.evaluate_point(transformation, ip, context) == MappingStatus::VALID, "Mapping evaluation failed during PhysicalPositionCoefficient evaluation.");
+        return m_f(context.physical_position);
+    }
+}

src/lib/serif/mapping/compactification/kelvin.cpp

@@ -203,6 +203,10 @@ namespace serif::mapping::compactification {
         return MappingStatus::VALID;
     }
 
+    std::unique_ptr<const ExteriorDomainMap> KelvinCompactification::clone() const {
+        return std::make_unique<const KelvinCompactification>(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;

src/lib/serif/mapping/data.cpp

@@ -88,7 +88,7 @@ namespace serif::mapping {
         }
 
         if (element.GetDerivType() != mfem::FiniteElement::GRAD) {
-            serif_error<std::invalid_argument, SERiFErrorCode::INVALID_COMPACTIFICATION_DERIVATIVE_VALUE>("Compactification coordinate finite element must provide a gradient");
+            serif_error<std::invalid_argument, SERiFErrorCode::INVALID_COMPACTIFICATION_DERIVATIVE_TYPE>("Compactification coordinate finite element must provide a gradient");
         }
 
         if (element.GetDof() <= 0) {

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 <array>
+#include <tuple>
+
+// 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<const compactification::ExteriorDomainMap> exterior_map
+    ) :
+    m_geometryCache(std::make_unique<DMGeometryCache>()),
+    m_options(options),
+    m_exterior_map(std::move(exterior_map)) {
+        if (m_options.dimension <= 0) {
+            serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DIMENSION>("Domain Mapper requires a positive, non-zero, dimension. Received dimension: {}", m_options.dimension);
+        }
+
+        if (m_options.vacuum_element_attribute_ID <= 0) {
+            serif_error<std::invalid_argument, SERiFErrorCode::INVALID_VACUUM_ELEMENT_ATTRIBUTE_ID>("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<std::invalid_argument, SERiFErrorCode::MISSING_EXTERIOR_MAP>("Domain Mapper requires a valid exterior map. One was not provided or a nullptr was provided.");
+        }
+    }
+
+    std::unique_ptr<DomainMapper> DomainMapper::clone() const {
+        if (!m_exterior_map) {
+            serif_error<std::logic_error, SERiFErrorCode::MISSING_EXTERIOR_MAP>("Domain Mapper clone failed because the exterior map is missing or a nullptr.");
+        }
+
+        auto result = std::make_unique<DomainMapper>(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<std::invalid_argument, SERiFErrorCode::INVALID_DIMENSION>("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<std::invalid_argument, SERiFErrorCode::INVALID_DIMENSION>("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<std::invalid_argument, SERiFErrorCode::INVALID_ELEMENT_GEOMETRY>("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<std::invalid_argument, SERiFErrorCode::INVALID_DIMENSION>("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<std::invalid_argument, SERiFErrorCode::INVALID_DIMENSION>("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<std::invalid_argument, SERiFErrorCode::INVALID_DIMENSION>("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<std::invalid_argument, SERiFErrorCode::INVALID_ELEMENT_GEOMETRY>("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<std::invalid_argument, SERiFErrorCode::INVALID_COMPACTIFICATION_RANK>("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<std::invalid_argument, SERiFErrorCode::INVALID_COMPACTIFICATION_TYPE>("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<std::invalid_argument, SERiFErrorCode::INVALID_COMPACTIFICATION_DERIVATIVE_TYPE>("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<std::invalid_argument, SERiFErrorCode::INVALID_DOF_COUNT>("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<std::invalid_argument, SERiFErrorCode::INVALID_DOF_COUNT>("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<std::invalid_argument, SERiFErrorCode::INVALID_DIMENSION>("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<std::invalid_argument, SERiFErrorCode::COMPACTIFICATION_MISMATCH>("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<SERiFErrorCode::MISSING_FACE_ELEMENT>(transformations.Elem1 != nullptr, "Face element transformation for ELEMENT_1 is null.");
+            return *transformations.Elem1;
+        }
+
+        serif_assert<SERiFErrorCode::MISSING_FACE_ELEMENT>(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();
 
     }
 }

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<SERiFErrorCode::MISSING_DISPLACEMENT_FINITE_ELEMENT_SPACE>(m_displacement_space != nullptr, "Displacement grid function is not associated with a finite-element space.");
+        serif_assert<SERiFErrorCode::MISSING_DISPLACEMENT_FINITE_ELEMENT_SPACE>(m_compactification_space != nullptr, "Compactification-coordinate grid function is not associated with a finite-element space.");
+        serif_assert<SERiFErrorCode::MISMATCHED_MESHING_MAPS>(m_displacement_space->GetMesh() == m_compactification_space->GetMesh(), "Displacement and compactification-coordinate grid functions are not associated with the same mesh.");
+        serif_assert<SERiFErrorCode::INVALID_DIMENSION>(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<SERiFErrorCode::INVALID_COMPACTIFICATION_RANK>(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<SERiFErrorCode::INVALID_DIMENSION>(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<SERiFErrorCode::CACHE_FE_BINDING_ERROR>(m_displacement.FESpace() == m_displacement_space, "Displacement grid function is no longer associated with the cached finite-element space.");
+        serif_assert<SERiFErrorCode::CACHE_FE_BINDING_ERROR>(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<SERiFErrorCode::INVALID_ELEMENT_ID_ERROR>(element_id >= 0, "Element ID must be non-negative. Instead {} was provided.", element_id);
+        serif_ensure<SERiFErrorCode::INVALID_ELEMENT_ID_ERROR>(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<ElementDisplacementData>(element_displacement_data_from_element_vdofs(displacement_element, m_element_displacement));
+        m_compactification_data = std::make_unique<ElementCompactificationData>(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<SERiFErrorCode::ELEMENT_NULLPTR_ERROR>(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<SERiFErrorCode::INVALID_DOMAIN_ERROR>(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<SERiFErrorCode::INVALID_DOMAIN_ERROR>(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<SERiFErrorCode::INVALID_DOMAIN_ERROR>(evaluate_point(transformation, integration_point, context) == MappingStatus::VALID, "Failed to evaluate the mapping at the given point.");
+        physical_position = context.physical_position;
+    }
+}

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<SERiFErrorCode::INVALID_DIMENSION>(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<SERiFErrorCode::INVALID_DIMENSION>(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<SERiFErrorCode::INVALID_DIMENSION>(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<SERiFErrorCode::INVALID_DIMENSION>(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<SERiFErrorCode::INVALID_DIMENSION>(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<SERiFErrorCode::INVALID_DIMENSION>(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<SERiFErrorCode::JACOBIAN_SHAPE_ERROR>((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<SERiFErrorCode::INVERTED_JACOBIAN>(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<SERiFErrorCode::JACOBIAN_SHAPE_ERROR>((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<SERiFErrorCode::INVERTED_JACOBIAN>(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<SERiFErrorCode::INVALID_DIMENSION>(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<SERiFErrorCode::INVALID_DIMENSION>(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<SERiFErrorCode::INVALID_DIMENSION>(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<SERiFErrorCode::INVALID_DIMENSION>(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<SERiFErrorCode::INVERTED_JACOBIAN>(determinate > 0.0, "Mapping jacobian determinant is non-positive ({}). The mapping may be inverted or degenerate.", determinate);
+        serif_assert<SERiFErrorCode::NON_FINITE_DETERMINATE>(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);
+            }
+        }
+    }
+}

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<SERiFErrorCode::MATRIX_SHAPE_ERROR>(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<SERiFErrorCode::INVALID_DIMENSION>(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<std::size_t>(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<SERiFErrorCode::OUT_OF_RANGE_ERROR>(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<std::size_t>(point) * m_point_stride;
+    }
+
+    void VolumeMappingCache::store(int point, VolumeMappingContext &context) {
+        validate_point(point);
+        double *data = m_data.data() + static_cast<std::size_t>(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;
+    }
+}

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(

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() {