// Ported from MeanField tests/mapping/domain_mapper.cpp.

#include <algorithm>
#include <array>
#include <cmath>
#include <format>
#include <limits>
#include <memory>
#include <stdexcept>
#include <utility>

#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <mfem.hpp>
#include <stroid/stroid.h>

#include "serif/backend/mfem/mapping/compactification/kelvin.hpp"
#include "serif/backend/mfem/mapping/domain_mapper.hpp"
#include "serif/backend/mfem/mapping/evaluator.hpp"
#include "serif/backend/mfem/mapping/transformations.hpp"

#include "serif/tests/mapping/volume_variation.hpp"
#include "serif/tests/test_tags.hpp"

namespace mapping = serif::mapping;
namespace user_options = serif::utils::user_options;
using Catch::Matchers::WithinAbs;

namespace {
    constexpr int dimension = 3;
    constexpr double tolerance = 1.0e-12;

    std::unique_ptr<const mapping::compactification::ExteriorDomainMap> make_kelvin_compactification() {
        return std::make_unique<mapping::compactification::KelvinCompactification>(
            user_options::KelvinCompactificationOptions{.r_star_ref = 1.0, .r_inf_ref = 4.0}
        );
    }

    mfem::DenseMatrix make_identity_matrix(const int size) {
        mfem::DenseMatrix matrix(size);
        matrix = 0.0;
        for (int i = 0; i < size; ++i) {
            matrix(i, i) = 1.0;
        }
        return matrix;
    }

    mfem::Vector make_constant_compactification_dofs(const mfem::FiniteElement& element, const double value = 0.0) {
        mfem::Vector dofs(element.GetDof());
        dofs = value;
        return dofs;
    }

    template <typename Function>
    mfem::Vector make_compactification_element_dofs(
        const mfem::FiniteElement& element,
        mfem::ElementTransformation& transformation,
        Function&& function
    ) {
        const mfem::IntegrationRule& nodes = element.GetNodes();
        REQUIRE(nodes.GetNPoints() == element.GetDof());

        mfem::Vector dofs(element.GetDof());
        mfem::Vector reference_position(dimension);

        for (int i = 0; i < element.GetDof(); ++i) {
            transformation.Transform(nodes.IntPoint(i), reference_position);
            dofs(i) = function(reference_position);
        }

        return dofs;
    }

    class ElementMappingDataOwner {
    public:
        explicit ElementMappingDataOwner(const mapping::ElementDisplacementData& displacement) :
            m_compactification(
                displacement.get_element(),
                make_constant_compactification_dofs(displacement.get_element())
            ),
            m_element_data{.displacement = displacement, .compactification = m_compactification} {}

        ElementMappingDataOwner(
            const mapping::ElementDisplacementData& displacement,
            const mfem::FiniteElement& compactification_element,
            const mfem::Vector& compactification_dofs
        ) :
            m_compactification(compactification_element, compactification_dofs),
            m_element_data{.displacement = displacement, .compactification = m_compactification} {}

        ElementMappingDataOwner(const ElementMappingDataOwner&) = delete;
        ElementMappingDataOwner& operator=(const ElementMappingDataOwner&) = delete;
        ElementMappingDataOwner(ElementMappingDataOwner&&) = delete;
        ElementMappingDataOwner& operator=(ElementMappingDataOwner&&) = delete;

        [[nodiscard]] const mapping::ElementMappingData& get() const noexcept {
            return m_element_data;
        }

    private:
        mapping::ElementCompactificationData m_compactification;
        mapping::ElementMappingData m_element_data;
    };

    void check_vector(
        const mfem::Vector& actual,
        const mfem::Vector& expected,
        const double comparison_tolerance = tolerance
    ) {
        REQUIRE(actual.Size() == expected.Size());
        for (int i = 0; i < actual.Size(); ++i) {
            CHECK_THAT(actual(i), WithinAbs(expected(i), comparison_tolerance));
        }
    }

    void check_matrix(
        const mfem::DenseMatrix& actual,
        const mfem::DenseMatrix& expected,
        const double comparison_tolerance = tolerance
    ) {
        REQUIRE(actual.Height() == expected.Height());
        REQUIRE(actual.Width() == expected.Width());

        for (int i = 0; i < actual.Height(); ++i) {
            for (int j = 0; j < actual.Width(); ++j) {
                CHECK_THAT(actual(i, j), WithinAbs(expected(i, j), comparison_tolerance));
            }
        }
    }

    struct SingleElementFixture {
        mfem::Mesh mesh;
        mfem::H1_FECollection displacement_collection;
        mfem::FiniteElementSpace displacement_space;

        SingleElementFixture() :
            mesh(mfem::Mesh::MakeCartesian3D(1, 1, 1, mfem::Element::HEXAHEDRON, 2.0, 3.0, 4.0)),
            displacement_collection(1, dimension),
            displacement_space(&mesh, &displacement_collection, dimension, mfem::Ordering::byVDIM) {}

        [[nodiscard]] const mfem::FiniteElement& get_element() const {
            return *displacement_space.GetFE(0);
        }

        [[nodiscard]] mfem::Vector make_zero_element_dofs() const {
            mfem::Vector element_dofs(get_element().GetDof() * dimension);
            element_dofs = 0.0;
            return element_dofs;
        }
    };

    mfem::Vector make_vector(const double x, const double y, const double z) {
        mfem::Vector vector(3);
        vector(0) = x;
        vector(1) = y;
        vector(2) = z;
        return vector;
    }

    mfem::DenseMatrix make_affine_displacement_gradient() {
        mfem::DenseMatrix gradient(3);
        gradient(0, 0) = 0.10;
        gradient(0, 1) = 0.04;
        gradient(0, 2) = -0.02;
        gradient(1, 0) = -0.03;
        gradient(1, 1) = 0.08;
        gradient(1, 2) = 0.01;
        gradient(2, 0) = 0.02;
        gradient(2, 1) = -0.01;
        gradient(2, 2) = -0.05;
        return gradient;
    }

    mfem::Vector make_affine_element_dofs(
        const mfem::FiniteElement& element,
        mfem::ElementTransformation& transformation,
        const mfem::DenseMatrix& displacement_gradient,
        const mfem::Vector& displacement_offset,
        const mfem::Ordering::Type ordering
    ) {
        const int dof_count = element.GetDof();
        const int field_dimension = displacement_offset.Size();
        const mfem::IntegrationRule& nodes = element.GetNodes();

        REQUIRE(nodes.GetNPoints() == dof_count);

        mfem::Vector element_dofs(dof_count * field_dimension);
        mfem::Vector reference_position(field_dimension);
        mfem::Vector displacement(field_dimension);

        for (int i = 0; i < dof_count; ++i) {
            transformation.Transform(nodes.IntPoint(i), reference_position);
            displacement_gradient.Mult(reference_position, displacement);
            displacement += displacement_offset;

            for (int component = 0; component < field_dimension; ++component) {
                const int index =
                    ordering == mfem::Ordering::byNODES ? i + component * dof_count : component + i * field_dimension;
                element_dofs(index) = displacement(component);
            }
        }

        return element_dofs;
    }

    mfem::DenseMatrix make_deformation_jacobian(const mfem::DenseMatrix& displacement_gradient) {
        mfem::DenseMatrix deformation_jacobian = make_identity_matrix(displacement_gradient.Height());
        deformation_jacobian.Add(1.0, displacement_gradient);
        return deformation_jacobian;
    }

    mfem::Vector evaluate_affine_physical_position(
        const mfem::Vector& reference_position,
        const mfem::DenseMatrix& displacement_gradient,
        const mfem::Vector& displacement_offset
    ) {
        mfem::Vector physical_position(reference_position);
        mfem::Vector displacement(reference_position.Size());
        displacement_gradient.Mult(reference_position, displacement);
        physical_position += displacement;
        physical_position += displacement_offset;
        return physical_position;
    }

    void check_point_context(
        const mapping::MappingPointContext& actual,
        const mapping::MappingPointContext& expected,
        const double comparison_tolerance = tolerance
    ) {
        CHECK(actual.compactified == expected.compactified);
        check_vector(actual.reference_position, expected.reference_position, comparison_tolerance);
        check_vector(actual.displaced_position, expected.displaced_position, comparison_tolerance);
        check_vector(actual.physical_position, expected.physical_position, comparison_tolerance);
        check_matrix(actual.displacement_jacobian, expected.displacement_jacobian, comparison_tolerance);
        check_matrix(actual.mapping_jacobian, expected.mapping_jacobian, comparison_tolerance);
        check_matrix(actual.inverse_mapping_jacobian, expected.inverse_mapping_jacobian, comparison_tolerance);
        CHECK_THAT(actual.mapping_determinate, WithinAbs(expected.mapping_determinate, comparison_tolerance));
    }

    constexpr double polynomial_tolerance = 2.0e-11;
    constexpr double difference_step = 2.0e-6;

    struct QuadraticElementFixture {
        mfem::Mesh mesh;
        mfem::H1_FECollection displacement_collection;
        mfem::FiniteElementSpace displacement_space;

        QuadraticElementFixture() :
            mesh(mfem::Mesh::MakeCartesian3D(1, 1, 1, mfem::Element::HEXAHEDRON, 2.0, 3.0, 4.0)),
            displacement_collection(2, dimension),
            displacement_space(&mesh, &displacement_collection, dimension, mfem::Ordering::byVDIM) {}

        [[nodiscard]] const mfem::FiniteElement& get_element() const {
            return *displacement_space.GetFE(0);
        }
    };

    template <typename Function>
    mfem::Vector make_function_element_dofs(
        const mfem::FiniteElement& element,
        mfem::ElementTransformation& transformation,
        Function&& function,
        const mfem::Ordering::Type ordering
    ) {
        const int dof_count = element.GetDof();
        const mfem::IntegrationRule& nodes = element.GetNodes();

        REQUIRE(nodes.GetNPoints() == dof_count);

        mfem::Vector element_dofs(dof_count * dimension);
        mfem::Vector reference_position(dimension);
        mfem::Vector value(dimension);

        for (int i = 0; i < dof_count; ++i) {
            transformation.Transform(nodes.IntPoint(i), reference_position);
            function(reference_position, value);

            for (int component = 0; component < dimension; ++component) {
                const int index =
                    ordering == mfem::Ordering::byNODES ? i + component * dof_count : component + i * dimension;
                element_dofs(index) = value(component);
            }
        }

        return element_dofs;
    }

    void evaluate_quadratic_displacement(const mfem::Vector& position, mfem::Vector& displacement) {
        const double x = position(0);
        const double y = position(1);
        const double z = position(2);

        displacement.SetSize(dimension);
        displacement(0) = 0.01 + 0.010 * x * x + 0.005 * y * z;
        displacement(1) = -0.02 - 0.004 * x * y + 0.006 * z * z;
        displacement(2) = 0.015 + 0.003 * x * z - 0.002 * y * y;
    }

    mfem::DenseMatrix evaluate_quadratic_displacement_gradient(const mfem::Vector& position) {
        const double x = position(0);
        const double y = position(1);
        const double z = position(2);

        mfem::DenseMatrix gradient(dimension);
        gradient(0, 0) = 0.020 * x;
        gradient(0, 1) = 0.005 * z;
        gradient(0, 2) = 0.005 * y;
        gradient(1, 0) = -0.004 * y;
        gradient(1, 1) = -0.004 * x;
        gradient(1, 2) = 0.012 * z;
        gradient(2, 0) = 0.003 * z;
        gradient(2, 1) = -0.004 * y;
        gradient(2, 2) = 0.003 * x;
        return gradient;
    }

    void evaluate_quadratic_direction(const mfem::Vector& position, mfem::Vector& direction) {
        const double x = position(0);
        const double y = position(1);
        const double z = position(2);

        direction.SetSize(dimension);
        direction(0) = 0.020 * x - 0.010 * y * z;
        direction(1) = -0.015 * y + 0.005 * x * z;
        direction(2) = 0.010 * z + 0.004 * x * y;
    }

    mfem::DenseMatrix evaluate_quadratic_direction_gradient(const mfem::Vector& position) {
        const double x = position(0);
        const double y = position(1);
        const double z = position(2);

        mfem::DenseMatrix gradient(dimension);
        gradient(0, 0) = 0.020;
        gradient(0, 1) = -0.010 * z;
        gradient(0, 2) = -0.010 * y;
        gradient(1, 0) = 0.005 * z;
        gradient(1, 1) = -0.015;
        gradient(1, 2) = 0.005 * x;
        gradient(2, 0) = 0.004 * y;
        gradient(2, 1) = 0.004 * x;
        gradient(2, 2) = 0.010;
        return gradient;
    }

    void check_vector_central_difference(
        const mfem::Vector& plus,
        const mfem::Vector& minus,
        const mfem::Vector& expected,
        const double step,
        const double comparison_tolerance
    ) {
        REQUIRE(plus.Size() == minus.Size());
        REQUIRE(plus.Size() == expected.Size());

        for (int i = 0; i < expected.Size(); ++i) {
            const double finite_difference = (plus(i) - minus(i)) / (2.0 * step);
            CHECK_THAT(finite_difference, WithinAbs(expected(i), comparison_tolerance));
        }
    }

    void check_matrix_central_difference(
        const mfem::DenseMatrix& plus,
        const mfem::DenseMatrix& minus,
        const mfem::DenseMatrix& expected,
        const double step,
        const double comparison_tolerance
    ) {
        REQUIRE(plus.Height() == minus.Height());
        REQUIRE(plus.Width() == minus.Width());
        REQUIRE(plus.Height() == expected.Height());
        REQUIRE(plus.Width() == expected.Width());

        for (int i = 0; i < expected.Height(); ++i) {
            for (int j = 0; j < expected.Width(); ++j) {
                const double finite_difference = (plus(i, j) - minus(i, j)) / (2.0 * step);
                CHECK_THAT(finite_difference, WithinAbs(expected(i, j), comparison_tolerance));
            }
        }
    }

    struct QuadraticMappingData {
        mfem::Vector base_dofs;
        mfem::Vector direction_dofs;
        mfem::Vector plus_dofs;
        mfem::Vector minus_dofs;

        QuadraticMappingData(const mfem::FiniteElement& element, mfem::ElementTransformation& transformation) {
            base_dofs = make_function_element_dofs(
                element, transformation, evaluate_quadratic_displacement, mfem::Ordering::byNODES
            );
            direction_dofs = make_function_element_dofs(
                element, transformation, evaluate_quadratic_direction, mfem::Ordering::byNODES
            );
            plus_dofs = base_dofs;
            minus_dofs = base_dofs;
            plus_dofs.Add(difference_step, direction_dofs);
            minus_dofs.Add(-difference_step, direction_dofs);
        }
    };

    void check_scalar_relative(
        const double actual,
        const double expected,
        const double relative_tolerance,
        const double absolute_tolerance = 1.0e-11
    ) {
        CHECK_THAT(actual, WithinAbs(expected, absolute_tolerance + relative_tolerance * std::abs(expected)));
    }

    void check_vector_central_difference_relative(
        const mfem::Vector& plus,
        const mfem::Vector& minus,
        const mfem::Vector& expected,
        const double step,
        const double relative_tolerance
    ) {
        REQUIRE(plus.Size() == minus.Size());
        REQUIRE(plus.Size() == expected.Size());

        for (int i = 0; i < expected.Size(); ++i) {
            const double finite_difference = (plus(i) - minus(i)) / (2.0 * step);
            check_scalar_relative(finite_difference, expected(i), relative_tolerance);
        }
    }

    void check_matrix_central_difference_relative(
        const mfem::DenseMatrix& plus,
        const mfem::DenseMatrix& minus,
        const mfem::DenseMatrix& expected,
        const double step,
        const double relative_tolerance
    ) {
        REQUIRE(plus.Height() == minus.Height());
        REQUIRE(plus.Width() == minus.Width());
        REQUIRE(plus.Height() == expected.Height());
        REQUIRE(plus.Width() == expected.Width());

        for (int i = 0; i < expected.Height(); ++i) {
            for (int j = 0; j < expected.Width(); ++j) {
                const double finite_difference = (plus(i, j) - minus(i, j)) / (2.0 * step);
                check_scalar_relative(finite_difference, expected(i, j), relative_tolerance);
            }
        }
    }

    mfem::Vector matrix_curl(const mfem::DenseMatrix& gradient) {
        return make_vector(
            gradient(2, 1) - gradient(1, 2), gradient(0, 2) - gradient(2, 0), gradient(1, 0) - gradient(0, 1)
        );
    }

    void check_centered_difference(
        const double analytic,
        const double finite_difference,
        const double plus_value,
        const double minus_value,
        const double step,
        const double relative_tolerance = 2.0e-6,
        const double absolute_tolerance = 1.0e-11
    ) {
        const double derivative_scale = std::max(std::abs(analytic), std::abs(finite_difference));
        const double primal_scale = std::max(std::abs(plus_value), std::abs(minus_value));
        const double roundoff_tolerance = 8.0 * std::numeric_limits<double>::epsilon() * primal_scale / step;
        const double tolerance = absolute_tolerance + relative_tolerance * derivative_scale + roundoff_tolerance;

        CHECK_THAT(finite_difference, Catch::Matchers::WithinAbs(analytic, tolerance));
    }

    double relative_vector_difference(const mfem::Vector& lhs, const mfem::Vector& rhs) {
        mfem::Vector difference(lhs);
        difference -= rhs;

        const double scale = std::max({lhs.Norml2(), rhs.Norml2(), 1.0e-12});
        return difference.Norml2() / scale;
    }

    double relative_matrix_difference(const mfem::DenseMatrix& lhs, const mfem::DenseMatrix& rhs) {
        mfem::DenseMatrix difference(lhs);
        difference -= rhs;

        const double scale = std::max({lhs.FNorm(), rhs.FNorm(), 1.0e-12});
        return difference.FNorm() / scale;
    }
} // namespace

TEST_CASE(
    "Element Displacement Data Rejects Invalid Vector Sizes",
    tags::unit & tags::mapping
) {
    SingleElementFixture fixture;

    const mfem::FiniteElement& element = fixture.get_element();
    const int dof_count = element.GetDof();

    mfem::Vector empty_dofs;
    mfem::Vector incomplete_dofs(dof_count * dimension - 1);
    incomplete_dofs = 0.0;

    CHECK_THROWS_AS(
        mapping::ElementDisplacementData(element, empty_dofs), std::invalid_argument
    );
    CHECK_THROWS_AS(
        mapping::ElementDisplacementData(element, incomplete_dofs), std::invalid_argument
    );
}

TEST_CASE(
    "Domain Mapping Workspace Tracks Its Dimension",
    tags::unit & tags::mapping
) {
    mapping::DMWorkspace workspace(dimension);

    REQUIRE(workspace.get_dimension() == dimension);

    workspace.set_dimension(2);
    REQUIRE(workspace.get_dimension() == 2);

    workspace.set_dimension(dimension);
    REQUIRE(workspace.get_dimension() == dimension);

    CHECK_THROWS_AS(workspace.set_dimension(0), std::invalid_argument);
    CHECK_THROWS_AS(workspace.set_dimension(-1), std::invalid_argument);
    CHECK_THROWS_AS(mapping::DMWorkspace(0), std::invalid_argument);
}

TEST_CASE(
    "Stateless Domain Mapper Validates Its Configuration",
    tags::unit & tags::mapping
) {
    const user_options::DomainMapperOptions valid_options{.dimension = dimension, .vacuum_element_attribute_ID = 3};
    mapping::DomainMapper mapper(valid_options, make_kelvin_compactification());

    REQUIRE(mapper.get_dimension() == dimension);
    REQUIRE(mapper.get_exterior_map().get_name() == "KelvinCompactification");

    const user_options::DomainMapperOptions invalid_dimension{.dimension = 0, .vacuum_element_attribute_ID = 3};
    const user_options::DomainMapperOptions invalid_attribute{.dimension = dimension, .vacuum_element_attribute_ID = 0};

    CHECK_THROWS_AS(mapping::DomainMapper(invalid_dimension, make_kelvin_compactification()), std::invalid_argument);
    CHECK_THROWS_AS(mapping::DomainMapper(invalid_attribute, make_kelvin_compactification()), std::invalid_argument);

    std::unique_ptr<const mapping::compactification::ExteriorDomainMap> null_exterior_map;
    CHECK_THROWS_AS(mapping::DomainMapper(valid_options, std::move(null_exterior_map)), std::invalid_argument);
}

TEST_CASE(
    "Stateless Domain Mapper Preserves Identity Point Geometry",
    tags::unit & tags::mapping
) {
    SingleElementFixture fixture;

    const mfem::FiniteElement& element = fixture.get_element();
    mfem::Vector zero_dofs = fixture.make_zero_element_dofs();
    const mapping::ElementDisplacementData displacement(element, zero_dofs);
    const ElementMappingDataOwner element_data(displacement);

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );
    mapping::DMWorkspace workspace(dimension);
    mapping::MappingPointContext context;

    mfem::ElementTransformation* transformation = fixture.mesh.GetElementTransformation(0);
    REQUIRE_FALSE(mapper.is_compactified_element(*transformation));

    const mfem::IntegrationRule& integration_rule = mfem::IntRules.Get(transformation->GetGeometryType(), 4);
    const mfem::DenseMatrix identity = make_identity_matrix(dimension);

    for (int q = 0; q < integration_rule.GetNPoints(); ++q) {
        const mfem::IntegrationPoint& integration_point = integration_rule.IntPoint(q);
        mfem::Vector expected_position(dimension);
        transformation->Transform(integration_point, expected_position);

        REQUIRE(
            mapper.evaluate_point(element_data.get(), *transformation, integration_point, workspace, context) ==
            mapping::MappingStatus::VALID
        );

        CAPTURE(q);
        REQUIRE_FALSE(context.compactified);
        check_vector(context.reference_position, expected_position);
        check_vector(context.displaced_position, expected_position);
        check_vector(context.physical_position, expected_position);
        check_matrix(context.displacement_jacobian, identity);
        check_matrix(context.mapping_jacobian, identity);
        check_matrix(context.inverse_mapping_jacobian, identity);
        CHECK_THAT(context.mapping_determinate, WithinAbs(1.0, tolerance));
    }
}

TEST_CASE(
    "Stateless Domain Mapper Preserves Identity Volume Geometry",
    tags::unit & tags::mapping
) {
    SingleElementFixture fixture;

    const mfem::FiniteElement& element = fixture.get_element();
    mfem::Vector zero_dofs = fixture.make_zero_element_dofs();
    const mapping::ElementDisplacementData displacement(element, zero_dofs);
    const ElementMappingDataOwner element_data(displacement);

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );
    mapping::DMWorkspace workspace(dimension);
    mapping::VolumeMappingContext context;

    mfem::ElementTransformation* transformation = fixture.mesh.GetElementTransformation(0);
    const mfem::IntegrationRule& integration_rule = mfem::IntRules.Get(transformation->GetGeometryType(), 4);

    for (int q = 0; q < integration_rule.GetNPoints(); ++q) {
        const mfem::IntegrationPoint& integration_point = integration_rule.IntPoint(q);

        REQUIRE(
            mapper.evaluate_volume(element_data.get(), *transformation, integration_point, workspace, context) ==
            mapping::MappingStatus::VALID
        );

        transformation->SetIntPoint(&integration_point);

        mfem::DenseMatrix expected_inverse(dimension);
        mfem::CalcInverse(transformation->Jacobian(), expected_inverse);

        const double expected_weight = integration_point.weight * transformation->Weight();

        CAPTURE(q);
        REQUIRE_FALSE(context.mapping.compactified);
        check_matrix(context.quadrature.J_element_to_physical_inverse, expected_inverse);
        CHECK_THAT(context.quadrature.mapping_determinate, WithinAbs(1.0, tolerance));
        CHECK_THAT(context.quadrature.physical_integration_weight, WithinAbs(expected_weight, tolerance));
        CHECK(context.quadrature.physical_integration_weight > 0.0);
    }
}

TEST_CASE(
    "Stateless Domain Mapper Preserves Identity Face Geometry",
    tags::unit & tags::mapping
) {
    SingleElementFixture fixture;

    const mfem::FiniteElement& element = fixture.get_element();
    mfem::Vector zero_dofs = fixture.make_zero_element_dofs();
    const mapping::ElementDisplacementData displacement(element, zero_dofs);
    const ElementMappingDataOwner element_data(displacement);

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );
    mapping::DMWorkspace workspace(dimension);
    mapping::FaceMappingContext context;

    for (int boundary_element = 0; boundary_element < fixture.mesh.GetNBE(); ++boundary_element) {
        mfem::FaceElementTransformations* transformation = fixture.mesh.GetBdrFaceTransformations(boundary_element);
        REQUIRE(transformation != nullptr);
        REQUIRE(transformation->Elem1 != nullptr);

        const mfem::IntegrationRule& integration_rule = mfem::IntRules.Get(transformation->GetGeometryType(), 4);

        for (int q = 0; q < integration_rule.GetNPoints(); ++q) {
            const mfem::IntegrationPoint& integration_point = integration_rule.IntPoint(q);

            REQUIRE(
                mapper.evaluate_face(
                    element_data.get(), *transformation, mapping::FaceElementSide::ELEMENT_1, integration_point,
                    workspace, context
                ) == mapping::MappingStatus::VALID
            );

            transformation->SetAllIntPoints(&integration_point);

            mfem::Vector raw_normal(dimension);
            mfem::CalcOrtho(transformation->Jacobian(), raw_normal);

            const double raw_normal_magnitude = raw_normal.Norml2();
            mfem::Vector expected_normal(raw_normal);
            expected_normal /= raw_normal_magnitude;

            const double expected_surface_weight = integration_point.weight * raw_normal_magnitude;

            const mfem::IntegrationPoint element_integration_point = transformation->Elem1->GetIntPoint();

            mfem::Vector expected_position(dimension);
            transformation->Elem1->Transform(element_integration_point, expected_position);

            CAPTURE(boundary_element, q);
            REQUIRE_FALSE(context.mapping.compactified);
            check_vector(context.mapping.reference_position, expected_position);
            check_vector(context.mapping.displaced_position, expected_position);
            check_vector(context.mapping.physical_position, expected_position);
            check_vector(context.reference_normal, expected_normal);
            check_vector(context.quadrature.normal, expected_normal);
            CHECK_THAT(context.reference_surface_weight, WithinAbs(expected_surface_weight, tolerance));
            CHECK_THAT(context.physical_surface_weight, WithinAbs(expected_surface_weight, tolerance));
            CHECK_THAT(context.quadrature.ds, WithinAbs(expected_surface_weight, tolerance));
            CHECK_THAT(context.quadrature.v_dot_n_scale, WithinAbs(1.0, tolerance));
        }
    }
}

TEST_CASE(
    "Stateless Domain Mapper Matches Exact Affine Point Mapping",
    tags::unit & tags::mapping
) {
    SingleElementFixture fixture;

    const mfem::FiniteElement& element = fixture.get_element();
    mfem::ElementTransformation* transformation = fixture.mesh.GetElementTransformation(0);

    const mfem::DenseMatrix displacement_gradient = make_affine_displacement_gradient();
    const mfem::DenseMatrix deformation_jacobian = make_deformation_jacobian(displacement_gradient);
    mfem::Vector displacement_offset = make_vector(0.07, -0.04, 0.03);
    mfem::Vector element_dofs = make_affine_element_dofs(
        element, *transformation, displacement_gradient, displacement_offset, mfem::Ordering::byNODES
    );
    const mapping::ElementDisplacementData displacement(element, element_dofs);
    const ElementMappingDataOwner element_data(displacement);

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );
    mapping::DMWorkspace workspace(dimension);
    mapping::MappingPointContext context;

    mfem::DenseMatrix inverse_deformation_jacobian(dimension);
    mfem::CalcInverse(deformation_jacobian, inverse_deformation_jacobian);
    const double deformation_determinant = deformation_jacobian.Det();

    REQUIRE(deformation_determinant > 0.0);

    const mfem::IntegrationRule& integration_rule = mfem::IntRules.Get(transformation->GetGeometryType(), 6);

    for (int q = 0; q < integration_rule.GetNPoints(); ++q) {
        const mfem::IntegrationPoint& integration_point = integration_rule.IntPoint(q);
        mfem::Vector reference_position(dimension);
        transformation->Transform(integration_point, reference_position);

        mfem::Vector expected_position =
            evaluate_affine_physical_position(reference_position, displacement_gradient, displacement_offset);

        REQUIRE(
            mapper.evaluate_point(element_data.get(), *transformation, integration_point, workspace, context) ==
            mapping::MappingStatus::VALID
        );

        CAPTURE(q);
        REQUIRE_FALSE(context.compactified);
        check_vector(context.reference_position, reference_position);
        check_vector(context.displaced_position, expected_position);
        check_vector(context.physical_position, expected_position);
        check_matrix(context.displacement_jacobian, deformation_jacobian);
        check_matrix(context.mapping_jacobian, deformation_jacobian);
        check_matrix(context.inverse_mapping_jacobian, inverse_deformation_jacobian);
        CHECK_THAT(context.mapping_determinate, WithinAbs(deformation_determinant, tolerance));
    }
}

TEST_CASE(
    "Stateless Domain Mapper Composes Affine Volume Jacobians",
    tags::unit & tags::mapping
) {
    SingleElementFixture fixture;

    const mfem::FiniteElement& element = fixture.get_element();
    mfem::ElementTransformation* transformation = fixture.mesh.GetElementTransformation(0);

    const mfem::DenseMatrix displacement_gradient = make_affine_displacement_gradient();
    const mfem::DenseMatrix mapping_jacobian = make_deformation_jacobian(displacement_gradient);
    mfem::Vector displacement_offset = make_vector(0.07, -0.04, 0.03);
    mfem::Vector element_dofs = make_affine_element_dofs(
        element, *transformation, displacement_gradient, displacement_offset, mfem::Ordering::byNODES
    );
    const mapping::ElementDisplacementData displacement(element, element_dofs);
    const ElementMappingDataOwner element_data(displacement);

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );
    mapping::DMWorkspace workspace(dimension);
    mapping::VolumeMappingContext context;

    const double mapping_determinate = mapping_jacobian.Det();
    REQUIRE(mapping_determinate > 0.0);

    const mfem::IntegrationRule& integration_rule = mfem::IntRules.Get(transformation->GetGeometryType(), 6);

    for (int q = 0; q < integration_rule.GetNPoints(); ++q) {
        const mfem::IntegrationPoint& integration_point = integration_rule.IntPoint(q);
        transformation->SetIntPoint(&integration_point);

        mfem::DenseMatrix full_element_jacobian(dimension);
        mfem::DenseMatrix expected_inverse(dimension);
        mfem::Mult(mapping_jacobian, transformation->Jacobian(), full_element_jacobian);
        mfem::CalcInverse(full_element_jacobian, expected_inverse);

        const double expected_weight = integration_point.weight * transformation->Weight() * mapping_determinate;

        REQUIRE(
            mapper.evaluate_volume(element_data.get(), *transformation, integration_point, workspace, context) ==
            mapping::MappingStatus::VALID
        );

        CAPTURE(q);
        check_matrix(context.mapping.mapping_jacobian, mapping_jacobian);
        check_matrix(context.quadrature.J_element_to_physical_inverse, expected_inverse);
        CHECK_THAT(context.quadrature.mapping_determinate, WithinAbs(mapping_determinate, tolerance));
        CHECK_THAT(context.quadrature.physical_integration_weight, WithinAbs(expected_weight, tolerance));
    }
}

TEST_CASE(
    "Stateless Domain Mapper Applies Nanson Formula On Affine Faces",
    tags::unit & tags::mapping
) {
    SingleElementFixture fixture;

    const mfem::FiniteElement& element = fixture.get_element();
    mfem::ElementTransformation* element_transformation = fixture.mesh.GetElementTransformation(0);

    const mfem::DenseMatrix displacement_gradient = make_affine_displacement_gradient();
    const mfem::DenseMatrix mapping_jacobian = make_deformation_jacobian(displacement_gradient);
    mfem::Vector displacement_offset = make_vector(0.07, -0.04, 0.03);
    mfem::Vector element_dofs = make_affine_element_dofs(
        element, *element_transformation, displacement_gradient, displacement_offset, mfem::Ordering::byNODES
    );
    const mapping::ElementDisplacementData displacement(element, element_dofs);
    const ElementMappingDataOwner element_data(displacement);

    mfem::DenseMatrix inverse_mapping_jacobian(dimension);
    mfem::CalcInverse(mapping_jacobian, inverse_mapping_jacobian);
    const double mapping_determinate = mapping_jacobian.Det();

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );
    mapping::DMWorkspace workspace(dimension);
    mapping::FaceMappingContext context;

    for (int boundary_element = 0; boundary_element < fixture.mesh.GetNBE(); ++boundary_element) {
        mfem::FaceElementTransformations* transformation = fixture.mesh.GetBdrFaceTransformations(boundary_element);
        REQUIRE(transformation != nullptr);
        REQUIRE(transformation->Elem1 != nullptr);

        const mfem::IntegrationRule& integration_rule = mfem::IntRules.Get(transformation->GetGeometryType(), 6);

        for (int q = 0; q < integration_rule.GetNPoints(); ++q) {
            const mfem::IntegrationPoint& integration_point = integration_rule.IntPoint(q);
            transformation->SetAllIntPoints(&integration_point);

            mfem::Vector raw_normal(dimension);
            mfem::Vector mapped_normal(dimension);
            mfem::CalcOrtho(transformation->Jacobian(), raw_normal);

            inverse_mapping_jacobian.MultTranspose(raw_normal, mapped_normal);
            mapped_normal *= mapping_determinate;

            const double raw_normal_magnitude = raw_normal.Norml2();
            const double mapped_normal_magnitude = mapped_normal.Norml2();

            mfem::Vector expected_reference_normal(raw_normal);
            mfem::Vector expected_physical_normal(mapped_normal);
            expected_reference_normal /= raw_normal_magnitude;
            expected_physical_normal /= mapped_normal_magnitude;

            const double expected_reference_weight = integration_point.weight * raw_normal_magnitude;
            const double expected_physical_weight = integration_point.weight * mapped_normal_magnitude;
            const double expected_normal_scale = mapped_normal_magnitude / raw_normal_magnitude;

            REQUIRE(
                mapper.evaluate_face(
                    element_data.get(), *transformation, mapping::FaceElementSide::ELEMENT_1, integration_point,
                    workspace, context
                ) == mapping::MappingStatus::VALID
            );

            CAPTURE(boundary_element, q);
            check_matrix(context.mapping.mapping_jacobian, mapping_jacobian);
            check_vector(context.reference_normal, expected_reference_normal);
            check_vector(context.quadrature.normal, expected_physical_normal);
            CHECK_THAT(context.reference_surface_weight, WithinAbs(expected_reference_weight, tolerance));
            CHECK_THAT(context.physical_surface_weight, WithinAbs(expected_physical_weight, tolerance));
            CHECK_THAT(context.quadrature.ds, WithinAbs(expected_reference_weight, tolerance));
            CHECK_THAT(context.quadrature.v_dot_n_scale, WithinAbs(expected_normal_scale, tolerance));
        }
    }
}

TEST_CASE(
    "Stateless Domain Mapper Has No Cross State Contamination",
    tags::unit & tags::mapping
) {
    SingleElementFixture fixture;

    const mfem::FiniteElement& element = fixture.get_element();
    mfem::ElementTransformation* transformation = fixture.mesh.GetElementTransformation(0);

    const mfem::DenseMatrix gradient_a = make_affine_displacement_gradient();
    mfem::DenseMatrix gradient_b(3);
    gradient_b = 0.0;
    gradient_b(0, 0) = -0.06;
    gradient_b(0, 2) = 0.03;
    gradient_b(1, 0) = 0.02;
    gradient_b(1, 1) = 0.12;
    gradient_b(2, 1) = -0.04;
    gradient_b(2, 2) = 0.07;

    mfem::Vector offset_a = make_vector(0.07, -0.04, 0.03);
    mfem::Vector offset_b = make_vector(-0.05, 0.08, -0.02);

    mfem::Vector dofs_a =
        make_affine_element_dofs(element, *transformation, gradient_a, offset_a, mfem::Ordering::byNODES);
    mfem::Vector dofs_b =
        make_affine_element_dofs(element, *transformation, gradient_b, offset_b, mfem::Ordering::byNODES);

    const mapping::ElementDisplacementData displacement_a(element, dofs_a);
    const mapping::ElementDisplacementData displacement_b(element, dofs_b);
    const ElementMappingDataOwner element_data_a(displacement_a);
    const ElementMappingDataOwner element_data_b(displacement_b);

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );
    mapping::DMWorkspace workspace(dimension);

    const mfem::IntegrationPoint& integration_point = mfem::Geometries.GetCenter(transformation->GetGeometryType());

    mapping::MappingPointContext first_a;
    mapping::MappingPointContext result_b;
    mapping::MappingPointContext second_a;

    REQUIRE(
        mapper.evaluate_point(element_data_a.get(), *transformation, integration_point, workspace, first_a) ==
        mapping::MappingStatus::VALID
    );
    REQUIRE(
        mapper.evaluate_point(element_data_b.get(), *transformation, integration_point, workspace, result_b) ==
        mapping::MappingStatus::VALID
    );
    REQUIRE(
        mapper.evaluate_point(element_data_a.get(), *transformation, integration_point, workspace, second_a) ==
        mapping::MappingStatus::VALID
    );

    check_point_context(first_a, second_a, 0.0);

    mfem::Vector state_difference(result_b.physical_position);
    state_difference -= first_a.physical_position;
    CHECK(state_difference.Norml2() > 1.0e-3);
}

TEST_CASE(
    "Stateless Domain Mapper Reports Invalid Element States",
    tags::unit & tags::mapping
) {
    SingleElementFixture fixture;

    const mfem::FiniteElement& element = fixture.get_element();
    mfem::ElementTransformation* transformation = fixture.mesh.GetElementTransformation(0);
    const mfem::IntegrationPoint& integration_point = mfem::Geometries.GetCenter(transformation->GetGeometryType());

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );
    mapping::DMWorkspace workspace(dimension);
    mapping::MappingPointContext context;

    mfem::Vector non_finite_dofs = fixture.make_zero_element_dofs();
    non_finite_dofs(0) = std::numeric_limits<double>::quiet_NaN();
    const mapping::ElementDisplacementData non_finite_displacement(element, non_finite_dofs);
    const ElementMappingDataOwner non_finite_data(non_finite_displacement);

    CHECK(
        mapper.evaluate_point(non_finite_data.get(), *transformation, integration_point, workspace, context) ==
        mapping::MappingStatus::NON_FINITE_MAPPING_INPUT
    );

    mfem::DenseMatrix singular_gradient(dimension);
    singular_gradient = 0.0;
    for (int i = 0; i < dimension; ++i) {
        singular_gradient(i, i) = -1.0;
    }

    mfem::Vector zero_offset(dimension);
    zero_offset = 0.0;
    mfem::Vector singular_dofs =
        make_affine_element_dofs(element, *transformation, singular_gradient, zero_offset, mfem::Ordering::byNODES);
    const mapping::ElementDisplacementData singular_displacement(element, singular_dofs);
    const ElementMappingDataOwner singular_data(singular_displacement);

    CHECK(
        mapper.evaluate_point(singular_data.get(), *transformation, integration_point, workspace, context) ==
        mapping::MappingStatus::NON_POSITIVE_DETERMINATE
    );

    mfem::DenseMatrix inverted_gradient(dimension);
    inverted_gradient = 0.0;
    for (int i = 0; i < dimension; ++i) {
        inverted_gradient(i, i) = -2.0;
    }

    mfem::Vector inverted_dofs =
        make_affine_element_dofs(element, *transformation, inverted_gradient, zero_offset, mfem::Ordering::byNODES);
    const mapping::ElementDisplacementData inverted_displacement(element, inverted_dofs);
    const ElementMappingDataOwner inverted_data(inverted_displacement);

    CHECK(
        mapper.evaluate_point(inverted_data.get(), *transformation, integration_point, workspace, context) ==
        mapping::MappingStatus::NON_POSITIVE_DETERMINATE
    );

    mapping::DMWorkspace wrong_workspace(2);
    CHECK_THROWS_AS(
        mapper.evaluate_point(singular_data.get(), *transformation, integration_point, wrong_workspace, context),
        std::invalid_argument
    );

    mfem::Vector two_dimensional_dofs(element.GetDof() * 2);
    two_dimensional_dofs = 0.0;
    const mapping::ElementDisplacementData two_dimensional_displacement(
        element, two_dimensional_dofs);
    const ElementMappingDataOwner two_dimensional_data(two_dimensional_displacement);

    CHECK_THROWS_AS(
        mapper.evaluate_point(two_dimensional_data.get(), *transformation, integration_point, workspace, context),
        std::invalid_argument
    );
}

TEST_CASE(
    "Stateless Domain Mapper Matches Exact Quadratic Point Mapping",
    tags::unit & tags::mapping
) {
    QuadraticElementFixture fixture;

    const mfem::FiniteElement& element = fixture.get_element();
    mfem::ElementTransformation* transformation = fixture.mesh.GetElementTransformation(0);

    mfem::Vector element_dofs =
        make_function_element_dofs(element, *transformation, evaluate_quadratic_displacement, mfem::Ordering::byNODES);
    const mapping::ElementDisplacementData displacement(element, element_dofs);
    const ElementMappingDataOwner element_data(displacement);

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );
    mapping::DMWorkspace workspace(dimension);
    mapping::MappingPointContext context;

    const mfem::IntegrationRule& integration_rule = mfem::IntRules.Get(transformation->GetGeometryType(), 6);

    for (int q = 0; q < integration_rule.GetNPoints(); ++q) {
        const mfem::IntegrationPoint& integration_point = integration_rule.IntPoint(q);

        mfem::Vector reference_position(dimension);
        mfem::Vector expected_displacement(dimension);
        transformation->Transform(integration_point, reference_position);
        evaluate_quadratic_displacement(reference_position, expected_displacement);

        mfem::Vector expected_position(reference_position);
        expected_position += expected_displacement;

        const mfem::DenseMatrix displacement_gradient = evaluate_quadratic_displacement_gradient(reference_position);
        mfem::DenseMatrix expected_jacobian = make_identity_matrix(dimension);
        expected_jacobian.Add(1.0, displacement_gradient);

        mfem::DenseMatrix expected_inverse(dimension);
        mfem::CalcInverse(expected_jacobian, expected_inverse);
        const double expected_determinant = expected_jacobian.Det();

        REQUIRE(expected_determinant > 0.0);
        REQUIRE(
            mapper.evaluate_point(element_data.get(), *transformation, integration_point, workspace, context) ==
            mapping::MappingStatus::VALID
        );

        CAPTURE(q);
        REQUIRE_FALSE(context.compactified);
        check_vector(context.reference_position, reference_position, polynomial_tolerance);
        check_vector(context.displaced_position, expected_position, polynomial_tolerance);
        check_vector(context.physical_position, expected_position, polynomial_tolerance);
        check_matrix(context.displacement_jacobian, expected_jacobian, polynomial_tolerance);
        check_matrix(context.mapping_jacobian, expected_jacobian, polynomial_tolerance);
        check_matrix(context.inverse_mapping_jacobian, expected_inverse, polynomial_tolerance);
        CHECK_THAT(context.mapping_determinate, WithinAbs(expected_determinant, polynomial_tolerance));
    }
}

TEST_CASE(
    "Stateless Domain Mapper Point Linearization Matches Centered Differences",
    tags::unit & tags::mapping
) {
    constexpr double linearization_tolerance = 2.0e-9;

    QuadraticElementFixture fixture;

    const mfem::FiniteElement& element = fixture.get_element();
    mfem::ElementTransformation* transformation = fixture.mesh.GetElementTransformation(0);
    QuadraticMappingData mapping_data(element, *transformation);

    const mapping::ElementDisplacementData base_displacement(element, mapping_data.base_dofs);
    const mapping::ElementDisplacementData direction(element, mapping_data.direction_dofs);
    const mapping::ElementDisplacementData plus_displacement(element, mapping_data.plus_dofs);
    const mapping::ElementDisplacementData minus_displacement(element, mapping_data.minus_dofs);
    const ElementMappingDataOwner base_data(base_displacement);
    const ElementMappingDataOwner plus_data(plus_displacement);
    const ElementMappingDataOwner minus_data(minus_displacement);

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );
    mapping::DMWorkspace workspace(dimension);

    const mfem::IntegrationRule& integration_rule = mfem::IntRules.Get(transformation->GetGeometryType(), 6);

    for (int q = 0; q < integration_rule.GetNPoints(); ++q) {
        const mfem::IntegrationPoint& integration_point = integration_rule.IntPoint(q);

        mapping::MappingPointContext base_context;
        mapping::MappingPointContext plus_context;
        mapping::MappingPointContext minus_context;
        mapping::MappingPointVariation variation;

        REQUIRE(
            mapper.evaluate_point(base_data.get(), *transformation, integration_point, workspace, base_context) ==
            mapping::MappingStatus::VALID
        );
        REQUIRE(
            mapper.evaluate_point_variation(
                base_data.get(), direction, *transformation, integration_point, base_context, workspace, variation
            ) == mapping::MappingStatus::VALID
        );
        REQUIRE(
            mapper.evaluate_point(plus_data.get(), *transformation, integration_point, workspace, plus_context) ==
            mapping::MappingStatus::VALID
        );
        REQUIRE(
            mapper.evaluate_point(minus_data.get(), *transformation, integration_point, workspace, minus_context) ==
            mapping::MappingStatus::VALID
        );

        mfem::Vector reference_position(dimension);
        mfem::Vector expected_direction(dimension);
        transformation->Transform(integration_point, reference_position);
        evaluate_quadratic_direction(reference_position, expected_direction);

        const mfem::DenseMatrix expected_direction_gradient = evaluate_quadratic_direction_gradient(reference_position);

        CAPTURE(q);
        check_vector(variation.displacement_variation, expected_direction, polynomial_tolerance);
        check_vector(variation.physical_position_variation, expected_direction, polynomial_tolerance);
        check_matrix(variation.displacement_jacobian_variation, expected_direction_gradient, polynomial_tolerance);
        check_matrix(variation.mapping_jacobian_variation, expected_direction_gradient, polynomial_tolerance);

        check_vector_central_difference(
            plus_context.physical_position, minus_context.physical_position, variation.physical_position_variation,
            difference_step, linearization_tolerance
        );
        check_matrix_central_difference(
            plus_context.mapping_jacobian, minus_context.mapping_jacobian, variation.mapping_jacobian_variation,
            difference_step, linearization_tolerance
        );
        check_matrix_central_difference(
            plus_context.inverse_mapping_jacobian, minus_context.inverse_mapping_jacobian,
            variation.inverse_mapping_jacobian_variation, difference_step, linearization_tolerance
        );

        const double determinant_finite_difference =
            (plus_context.mapping_determinate - minus_context.mapping_determinate) / (2.0 * difference_step);
        CHECK_THAT(
            determinant_finite_difference, WithinAbs(variation.mapping_determinate_variation, linearization_tolerance)
        );
    }
}

TEST_CASE(
    "Stateless Domain Mapper Volume Linearization Matches Centered Differences",
    tags::unit & tags::mapping
) {
    constexpr double linearization_tolerance = 5.0e-9;

    QuadraticElementFixture fixture;

    const mfem::FiniteElement& element = fixture.get_element();
    mfem::ElementTransformation* transformation = fixture.mesh.GetElementTransformation(0);
    QuadraticMappingData mapping_data(element, *transformation);

    const mapping::ElementDisplacementData base_displacement(element, mapping_data.base_dofs);
    const mapping::ElementDisplacementData direction(element, mapping_data.direction_dofs);
    const mapping::ElementDisplacementData plus_displacement(element, mapping_data.plus_dofs);
    const mapping::ElementDisplacementData minus_displacement(element, mapping_data.minus_dofs);
    const ElementMappingDataOwner base_data(base_displacement);
    const ElementMappingDataOwner plus_data(plus_displacement);
    const ElementMappingDataOwner minus_data(minus_displacement);

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );
    mapping::DMWorkspace workspace(dimension);

    const mfem::IntegrationRule& integration_rule = mfem::IntRules.Get(transformation->GetGeometryType(), 6);

    for (int q = 0; q < integration_rule.GetNPoints(); ++q) {
        const mfem::IntegrationPoint& integration_point = integration_rule.IntPoint(q);

        mapping::VolumeMappingContext base_context;
        mapping::VolumeMappingContext plus_context;
        mapping::VolumeMappingContext minus_context;
        mapping::VolumeMappingVariation variation;

        REQUIRE(
            mapper.evaluate_volume(base_data.get(), *transformation, integration_point, workspace, base_context) ==
            mapping::MappingStatus::VALID
        );
        REQUIRE(
            serif::tests::evaluate_volume_variation(mapper,
                base_data.get(), direction, *transformation, integration_point, base_context, workspace, variation
            ) == mapping::MappingStatus::VALID
        );
        REQUIRE(
            mapper.evaluate_volume(plus_data.get(), *transformation, integration_point, workspace, plus_context) ==
            mapping::MappingStatus::VALID
        );
        REQUIRE(
            mapper.evaluate_volume(minus_data.get(), *transformation, integration_point, workspace, minus_context) ==
            mapping::MappingStatus::VALID
        );

        CAPTURE(q);
        const double determinant_finite_difference =
            (plus_context.quadrature.mapping_determinate - minus_context.quadrature.mapping_determinate) /
            (2.0 * difference_step);
        const double weight_finite_difference = (plus_context.quadrature.physical_integration_weight -
                                                    minus_context.quadrature.physical_integration_weight) /
                                                (2.0 * difference_step);

        CHECK_THAT(
            determinant_finite_difference,
            WithinAbs(variation.mapping.mapping_determinate_variation, linearization_tolerance)
        );
        CHECK_THAT(weight_finite_difference, WithinAbs(variation.weight_variation, linearization_tolerance));
    }
}

TEST_CASE(
    "Stateless Domain Mapper Face Linearization Matches Centered Differences",
    tags::unit & tags::mapping
) {
    constexpr double normal_tolerance = 2.0e-8;
    constexpr double measure_tolerance = 2.0e-8;

    QuadraticElementFixture fixture;

    const mfem::FiniteElement& element = fixture.get_element();
    mfem::ElementTransformation* element_transformation = fixture.mesh.GetElementTransformation(0);
    QuadraticMappingData mapping_data(element, *element_transformation);

    const mapping::ElementDisplacementData base_displacement(element, mapping_data.base_dofs);
    const mapping::ElementDisplacementData direction(element, mapping_data.direction_dofs);
    const mapping::ElementDisplacementData plus_displacement(element, mapping_data.plus_dofs);
    const mapping::ElementDisplacementData minus_displacement(element, mapping_data.minus_dofs);
    const ElementMappingDataOwner base_data(base_displacement);
    const ElementMappingDataOwner plus_data(plus_displacement);
    const ElementMappingDataOwner minus_data(minus_displacement);

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );
    mapping::DMWorkspace workspace(dimension);

    for (int boundary_element = 0; boundary_element < fixture.mesh.GetNBE(); ++boundary_element) {
        mfem::FaceElementTransformations* transformation = fixture.mesh.GetBdrFaceTransformations(boundary_element);
        REQUIRE(transformation != nullptr);
        REQUIRE(transformation->Elem1 != nullptr);

        const mfem::IntegrationRule& integration_rule = mfem::IntRules.Get(transformation->GetGeometryType(), 4);

        for (int q = 0; q < integration_rule.GetNPoints(); ++q) {
            const mfem::IntegrationPoint& integration_point = integration_rule.IntPoint(q);

            mapping::FaceMappingContext base_context;
            mapping::FaceMappingContext plus_context;
            mapping::FaceMappingContext minus_context;
            mapping::FaceMappingVariation variation;

            REQUIRE(
                mapper.evaluate_face(
                    base_data.get(), *transformation, mapping::FaceElementSide::ELEMENT_1, integration_point, workspace,
                    base_context
                ) == mapping::MappingStatus::VALID
            );
            REQUIRE(
                mapper.evaluate_face_variation(
                    base_data.get(), direction, *transformation, mapping::FaceElementSide::ELEMENT_1, integration_point,
                    base_context, workspace, variation
                ) == mapping::MappingStatus::VALID
            );
            REQUIRE(
                mapper.evaluate_face(
                    plus_data.get(), *transformation, mapping::FaceElementSide::ELEMENT_1, integration_point, workspace,
                    plus_context
                ) == mapping::MappingStatus::VALID
            );
            REQUIRE(
                mapper.evaluate_face(
                    minus_data.get(), *transformation, mapping::FaceElementSide::ELEMENT_1, integration_point,
                    workspace, minus_context
                ) == mapping::MappingStatus::VALID
            );

            CAPTURE(boundary_element, q);
            check_vector_central_difference(
                plus_context.quadrature.normal, minus_context.quadrature.normal, variation.physical_normal_variation,
                difference_step, normal_tolerance
            );

            const double surface_weight_finite_difference =
                (plus_context.physical_surface_weight - minus_context.physical_surface_weight) /
                (2.0 * difference_step);
            const double normal_scale_finite_difference =
                (plus_context.quadrature.v_dot_n_scale - minus_context.quadrature.v_dot_n_scale) /
                (2.0 * difference_step);

            CHECK_THAT(
                surface_weight_finite_difference,
                WithinAbs(variation.physical_surface_weight_variation, measure_tolerance)
            );
            CHECK_THAT(
                normal_scale_finite_difference, WithinAbs(variation.normal_flux_scale_variation, measure_tolerance)
            );
        }
    }
}

TEST_CASE(
    "Stateless Domain Mapper Rejects Invalid Linearization Directions",
    tags::unit & tags::mapping
) {
    QuadraticElementFixture fixture;

    const mfem::FiniteElement& element = fixture.get_element();
    mfem::ElementTransformation* transformation = fixture.mesh.GetElementTransformation(0);
    const mfem::IntegrationPoint& integration_point = mfem::Geometries.GetCenter(transformation->GetGeometryType());
    QuadraticMappingData mapping_data(element, *transformation);

    const mapping::ElementDisplacementData base_displacement(element, mapping_data.base_dofs);
    const ElementMappingDataOwner base_data(base_displacement);

    mfem::Vector non_finite_direction_dofs(mapping_data.direction_dofs);
    non_finite_direction_dofs(0) = std::numeric_limits<double>::quiet_NaN();
    const mapping::ElementDisplacementData non_finite_direction(
        element, non_finite_direction_dofs);

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );
    mapping::DMWorkspace workspace(dimension);

    mapping::MappingPointContext base_context;
    mapping::MappingPointVariation variation;

    REQUIRE(
        mapper.evaluate_point(base_data.get(), *transformation, integration_point, workspace, base_context) ==
        mapping::MappingStatus::VALID
    );
    CHECK(
        mapper.evaluate_point_variation(
            base_data.get(), non_finite_direction, *transformation, integration_point, base_context, workspace,
            variation
        ) == mapping::MappingStatus::NON_FINITE_MAPPING_INPUT
    );

    SingleElementFixture linear_fixture;
    const mfem::FiniteElement& linear_element = linear_fixture.get_element();
    mfem::Vector linear_direction_dofs = linear_fixture.make_zero_element_dofs();
    const mapping::ElementDisplacementData incompatible_direction(
        linear_element, linear_direction_dofs);

    CHECK_THROWS_AS(
        mapper.evaluate_point_variation(
            base_data.get(), incompatible_direction, *transformation, integration_point, base_context, workspace,
            variation
        ),
        std::invalid_argument
    );
}

TEST_CASE(
    "Mapping Field Transforms Satisfy Piola And Gradient Identities",
    tags::unit & tags::mapping
) {
    constexpr double transform_tolerance = 1.0e-12;

    mapping::MappingPointContext context;
    context.mapping_jacobian.SetSize(3);
    context.mapping_jacobian(0, 0) = 1.20;
    context.mapping_jacobian(0, 1) = 0.15;
    context.mapping_jacobian(0, 2) = -0.05;
    context.mapping_jacobian(1, 0) = -0.08;
    context.mapping_jacobian(1, 1) = 0.95;
    context.mapping_jacobian(1, 2) = 0.12;
    context.mapping_jacobian(2, 0) = 0.04;
    context.mapping_jacobian(2, 1) = -0.10;
    context.mapping_jacobian(2, 2) = 1.10;

    context.mapping_determinate = context.mapping_jacobian.Det();
    REQUIRE(context.mapping_determinate > 0.0);

    context.inverse_mapping_jacobian.SetSize(3);
    mfem::CalcInverse(context.mapping_jacobian, context.inverse_mapping_jacobian);

    mfem::Vector reference_gradient = make_vector(0.3, -0.6, 0.8);
    mfem::Vector reference_test_gradient = make_vector(-0.7, 0.2, 0.4);

    mfem::Vector physical_gradient;
    mfem::Vector recovered_gradient;
    mfem::Vector physical_test_gradient;

    mapping::map_reference_gradient_to_physical(context, reference_gradient, physical_gradient);
    mapping::map_physical_gradient_to_reference(context, physical_gradient, recovered_gradient);
    mapping::map_reference_gradient_to_physical(context, reference_test_gradient, physical_test_gradient);

    check_vector(recovered_gradient, reference_gradient, transform_tolerance);

    mfem::DenseMatrix reference_vector_gradient(3);
    reference_vector_gradient(0, 0) = 0.20;
    reference_vector_gradient(0, 1) = -0.10;
    reference_vector_gradient(0, 2) = 0.04;
    reference_vector_gradient(1, 0) = 0.03;
    reference_vector_gradient(1, 1) = 0.15;
    reference_vector_gradient(1, 2) = -0.08;
    reference_vector_gradient(2, 0) = -0.05;
    reference_vector_gradient(2, 1) = 0.02;
    reference_vector_gradient(2, 2) = 0.11;

    mfem::DenseMatrix physical_vector_gradient;
    mfem::DenseMatrix recovered_vector_gradient;

    mapping::map_reference_vector_gradient_to_physical(context, reference_vector_gradient, physical_vector_gradient);
    mapping::map_physical_vector_gradient_to_reference(context, physical_vector_gradient, recovered_vector_gradient);
    check_matrix(recovered_vector_gradient, reference_vector_gradient, transform_tolerance);

    mfem::DenseMatrix diffusion_tensor;
    mapping::compute_scalar_diffusion_tensor(context, diffusion_tensor);

    mfem::Vector diffusion_action(3);
    diffusion_tensor.Mult(reference_test_gradient, diffusion_action);

    const double physical_gradient_inner_product =
        context.mapping_determinate * (physical_gradient * physical_test_gradient);
    const double reference_gradient_inner_product = reference_gradient * diffusion_action;

    CHECK_THAT(physical_gradient_inner_product, WithinAbs(reference_gradient_inner_product, transform_tolerance));
}

TEST_CASE(
    "Stateless Domain Mapper Represents Strong Rotating Star Geometry",
    tags::integration & tags::mapping
) {
    constexpr double r_star = 1.0;
    constexpr double r_infinity = 4.0;
    constexpr double transform_tolerance = 2.0e-10;

    stroid::config::MeshConfig mesh_config;
    mesh_config.refinement_levels = 0;
    mesh_config.order = 2;
    mesh_config.include_external_domain = true;
    mesh_config.r_core = 0.25;
    mesh_config.r_star = r_star;
    mesh_config.r_infinity = r_infinity;
    mesh_config.flattening = 0.0;
    mesh_config.optimization_methods = stroid::config::OptimizationMethods{false, true};

    stroid::StroidMesh stroid_mesh = stroid::GenerateMesh(mesh_config);
    mfem::Mesh& mesh = *stroid_mesh.mesh;
    REQUIRE(stroid_mesh.exterior_coordinate != nullptr);
    REQUIRE(stroid_mesh.exterior_coordinate->space != nullptr);
    REQUIRE(stroid_mesh.exterior_coordinate->values != nullptr);

    mfem::FiniteElementSpace& compactification_space = *stroid_mesh.exterior_coordinate->space;
    mfem::GridFunction& compactification_coordinate = *stroid_mesh.exterior_coordinate->values;

    mfem::H1_FECollection displacement_collection(3, dimension);
    mfem::FiniteElementSpace displacement_space(&mesh, &displacement_collection, dimension, mfem::Ordering::byVDIM);
    mfem::GridFunction displacement(&displacement_space);

    auto rotating_displacement =
        [r_star, r_infinity](const mfem::Vector& reference_position, mfem::Vector& displacement_value) {
            const double x = reference_position(0);
            const double y = reference_position(1);
            const double z = reference_position(2);
            const double radius_squared = x * x + y * y + z * z;
            const double radius = std::sqrt(radius_squared);

            displacement_value.SetSize(3);
            displacement_value = 0.0;

            if (radius <= 1.0e-14) {
                return;
            }

            const double cylindrical_fraction = (x * x + y * y) / radius_squared;
            const double angular_deformation =
                0.20 * cylindrical_fraction + 0.12 * cylindrical_fraction * cylindrical_fraction;

            double radial_extension = 0.0;
            if (radius <= r_star) {
                radial_extension = radius_squared / (r_star * r_star);
            } else {
                radial_extension = std::max(0.0, (r_infinity - radius) / (r_infinity - r_star));
            }

            const double scale = radial_extension * angular_deformation;
            displacement_value(0) = scale * x;
            displacement_value(1) = scale * y;
            displacement_value(2) = scale * z;
        };

    mfem::VectorFunctionCoefficient displacement_coefficient(dimension, rotating_displacement);
    displacement.ProjectCoefficient(displacement_coefficient);

    std::unique_ptr<const mapping::compactification::ExteriorDomainMap> exterior_map =
        std::make_unique<mapping::compactification::KelvinCompactification>(
            user_options::KelvinCompactificationOptions{.r_star_ref = r_star, .r_inf_ref = r_infinity}
        );

    mapping::DomainMapper mapper({.dimension = dimension, .vacuum_element_attribute_ID = 3}, std::move(exterior_map));
    mapping::DMWorkspace workspace(dimension);

    double minimum_mapping_determinant = std::numeric_limits<double>::infinity();
    double maximum_mapping_determinant = 0.0;
    double stellar_volume = 0.0;
    double moment_x = 0.0;
    double moment_y = 0.0;
    double moment_z = 0.0;

    int stellar_elements = 0;
    int vacuum_elements = 0;

    for (int element_id = 0; element_id < mesh.GetNE(); ++element_id) {
        mfem::ElementTransformation* transformation = mesh.GetElementTransformation(element_id);
        const mfem::FiniteElement* element = displacement_space.GetFE(element_id);

        mfem::Array<int> element_vdofs;
        mfem::Vector element_dofs;
        displacement_space.GetElementVDofs(element_id, element_vdofs);
        displacement.GetSubVector(element_vdofs, element_dofs);

        const mapping::ElementDisplacementData element_displacement =
            mapping::element_displacement_data_from_element_vdofs(*element, element_dofs);

        mfem::Array<int> compactification_dof_indices;
        mfem::Vector compactification_dofs;
        compactification_space.GetElementDofs(element_id, compactification_dof_indices);
        compactification_coordinate.GetSubVector(compactification_dof_indices, compactification_dofs);

        const ElementMappingDataOwner element_data(
            element_displacement, *compactification_space.GetFE(element_id), compactification_dofs
        );
        const int quadrature_order = 2 * element->GetOrder() + 6;
        const mfem::IntegrationRule& integration_rule =
            mfem::IntRules.Get(transformation->GetGeometryType(), quadrature_order);

        if (transformation->Attribute == 3) {
            ++vacuum_elements;
        } else {
            ++stellar_elements;
        }

        for (int q = 0; q < integration_rule.GetNPoints(); ++q) {
            const mfem::IntegrationPoint& integration_point = integration_rule.IntPoint(q);
            mapping::VolumeMappingContext context;

            REQUIRE(
                mapper.evaluate_volume(element_data.get(), *transformation, integration_point, workspace, context) ==
                mapping::MappingStatus::VALID
            );

            minimum_mapping_determinant = std::min(minimum_mapping_determinant, context.mapping.mapping_determinate);
            maximum_mapping_determinant = std::max(maximum_mapping_determinant, context.mapping.mapping_determinate);

            REQUIRE(context.mapping.mapping_determinate > 0.0);
            REQUIRE(context.quadrature.physical_integration_weight > 0.0);

            mfem::Vector reference_gradient = make_vector(0.3, -0.5, 0.8);
            mfem::Vector physical_gradient;
            mfem::Vector recovered_gradient;
            mapping::map_reference_gradient_to_physical(context.mapping, reference_gradient, physical_gradient);
            mapping::map_physical_gradient_to_reference(context.mapping, physical_gradient, recovered_gradient);
            check_vector(recovered_gradient, reference_gradient, transform_tolerance);

            if (transformation->Attribute == 3) {
                transformation->SetIntPoint(&integration_point);

                mfem::Vector compactification_gradient(dimension);
                const double coordinate = compactification_coordinate.GetValue(element_id, integration_point);
                compactification_coordinate.GetGradient(*transformation, compactification_gradient);

                mapping::compactification::ExteriorMapResult direct_result;
                const mapping::compactification::ExteriorMapInput direct_input{
                    .reference_position = context.mapping.reference_position,
                    .displaced_position = context.mapping.displaced_position,
                    .displacement_jacobian = context.mapping.displacement_jacobian,
                    .compactification_coordinate_gradient = compactification_gradient,
                    .compactification_coordinate = coordinate
                };

                REQUIRE(
                    mapper.get_exterior_map().evaluate(direct_input, direct_result) == mapping::MappingStatus::VALID
                );
                check_vector(context.mapping.physical_position, direct_result.physical_position, transform_tolerance);
                check_matrix(context.mapping.mapping_jacobian, direct_result.mapping_jacobian, transform_tolerance);
            } else {
                check_vector(
                    context.mapping.physical_position, context.mapping.displaced_position, transform_tolerance
                );

                const double x = context.mapping.physical_position(0);
                const double y = context.mapping.physical_position(1);
                const double z = context.mapping.physical_position(2);
                const double weight = context.quadrature.physical_integration_weight;

                stellar_volume += weight;
                moment_x += x * x * weight;
                moment_y += y * y * weight;
                moment_z += z * z * weight;
            }
        }
    }

    REQUIRE(stellar_elements > 0);
    REQUIRE(vacuum_elements > 0);
    REQUIRE(stellar_volume > 0.0);
    REQUIRE(std::isfinite(minimum_mapping_determinant));
    REQUIRE(minimum_mapping_determinant > 0.0);

    const double moment_trace = moment_x + moment_y + moment_z;
    const double quadrupole_x = 3.0 * moment_x - moment_trace;
    const double quadrupole_y = 3.0 * moment_y - moment_trace;
    const double quadrupole_z = 3.0 * moment_z - moment_trace;
    const double normalized_quadrupole =
        std::sqrt(quadrupole_x * quadrupole_x + quadrupole_y * quadrupole_y + quadrupole_z * quadrupole_z) /
        moment_trace;
    const double axisymmetry_error = std::abs(moment_x - moment_y) / (0.5 * (moment_x + moment_y));

    INFO("Stellar volume = " << stellar_volume);
    INFO("Minimum mapping determinant = " << minimum_mapping_determinant);
    INFO("Maximum mapping determinant = " << maximum_mapping_determinant);
    INFO("Normalized geometric quadrupole = " << normalized_quadrupole);
    INFO("Axisymmetry error = " << axisymmetry_error);

    CHECK(moment_x > moment_z);
    CHECK(moment_y > moment_z);
    CHECK(normalized_quadrupole > 1.0e-2);
    CHECK(axisymmetry_error < 5.0e-2);
}

TEST_CASE(
    "Kelvin Composed Domain Mapping Linearization Matches Centered Differences",
    tags::unit & tags::mapping
) {
    constexpr double relative_tolerance = 2.0e-6;
    constexpr double kelvin_difference_step = 2.0e-4;

    mfem::Mesh mesh = mfem::Mesh::MakeCartesian3D(1, 1, 1, mfem::Element::HEXAHEDRON, 2.98, 0.02, 0.02);

    for (int vertex_id = 0; vertex_id < mesh.GetNV(); ++vertex_id) {
        double* vertex = mesh.GetVertex(vertex_id);
        vertex[0] += 1.0;
        vertex[1] -= 0.01;
        vertex[2] -= 0.01;
    }

    mesh.GetElement(0)->SetAttribute(3);
    mesh.SetAttributes();

    mfem::H1_FECollection displacement_collection(2, dimension);
    mfem::FiniteElementSpace displacement_space(&mesh, &displacement_collection, dimension, mfem::Ordering::byVDIM);

    const mfem::FiniteElement& element = *displacement_space.GetFE(0);
    mfem::ElementTransformation* transformation = mesh.GetElementTransformation(0);

    mfem::Vector base_dofs =
        make_function_element_dofs(element, *transformation, evaluate_quadratic_displacement, mfem::Ordering::byNODES);
    mfem::Vector direction_dofs =
        make_function_element_dofs(element, *transformation, evaluate_quadratic_direction, mfem::Ordering::byNODES);

    mfem::Vector plus_dofs(base_dofs);
    mfem::Vector minus_dofs(base_dofs);
    plus_dofs.Add(kelvin_difference_step, direction_dofs);
    minus_dofs.Add(-kelvin_difference_step, direction_dofs);

    const mapping::ElementDisplacementData base_displacement(element, base_dofs);
    const mapping::ElementDisplacementData direction(element, direction_dofs);
    const mapping::ElementDisplacementData plus_displacement(element, plus_dofs);
    const mapping::ElementDisplacementData minus_displacement(element, minus_dofs);

    mfem::Vector compactification_dofs =
        make_compactification_element_dofs(element, *transformation, [](const mfem::Vector& reference_position) {
            return (reference_position(0) - 1.0) / 3.0;
        });

    const ElementMappingDataOwner base_data(base_displacement, element, compactification_dofs);
    const ElementMappingDataOwner plus_data(plus_displacement, element, compactification_dofs);
    const ElementMappingDataOwner minus_data(minus_displacement, element, compactification_dofs);

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );
    mapping::DMWorkspace workspace(dimension);

    auto vector_central_difference = [](const mfem::Vector& plus_value, const mfem::Vector& minus_value) {
        mfem::Vector difference(plus_value);
        difference -= minus_value;
        difference *= 1.0 / (2.0 * kelvin_difference_step);
        return difference;
    };

    auto matrix_central_difference = [](const mfem::DenseMatrix& plus_value, const mfem::DenseMatrix& minus_value) {
        mfem::DenseMatrix difference(plus_value);
        difference -= minus_value;
        difference *= 1.0 / (2.0 * kelvin_difference_step);
        return difference;
    };

    for (const double xi : std::array{0.0, 0.25, 0.75, 0.95, 0.99}) {
        mfem::IntegrationPoint integration_point;
        integration_point.x = 3.0 * xi / 2.98;
        integration_point.y = 0.5;
        integration_point.z = 0.5;
        integration_point.weight = 0.73;

        mapping::VolumeMappingContext base_context;
        mapping::VolumeMappingContext plus_context;
        mapping::VolumeMappingContext minus_context;
        mapping::VolumeMappingVariation variation;

        REQUIRE(
            mapper.evaluate_volume(base_data.get(), *transformation, integration_point, workspace, base_context) ==
            mapping::MappingStatus::VALID
        );
        REQUIRE(
            serif::tests::evaluate_volume_variation(mapper,
                base_data.get(), direction, *transformation, integration_point, base_context, workspace, variation
            ) == mapping::MappingStatus::VALID
        );
        REQUIRE(
            mapper.evaluate_volume(plus_data.get(), *transformation, integration_point, workspace, plus_context) ==
            mapping::MappingStatus::VALID
        );
        REQUIRE(
            mapper.evaluate_volume(minus_data.get(), *transformation, integration_point, workspace, minus_context) ==
            mapping::MappingStatus::VALID
        );

        CAPTURE(xi);
        REQUIRE(base_context.mapping.compactified);

        mfem::Vector physical_position_difference =
            vector_central_difference(plus_context.mapping.physical_position, minus_context.mapping.physical_position);
        const mfem::DenseMatrix mapping_jacobian_difference =
            matrix_central_difference(plus_context.mapping.mapping_jacobian, minus_context.mapping.mapping_jacobian);
        const mfem::DenseMatrix inverse_mapping_jacobian_difference = matrix_central_difference(
            plus_context.mapping.inverse_mapping_jacobian, minus_context.mapping.inverse_mapping_jacobian
        );

        CHECK_THAT(
            relative_vector_difference(physical_position_difference, variation.mapping.physical_position_variation),
            Catch::Matchers::WithinAbs(0.0, relative_tolerance)
        );
        CHECK_THAT(
            relative_matrix_difference(mapping_jacobian_difference, variation.mapping.mapping_jacobian_variation),
            Catch::Matchers::WithinAbs(0.0, relative_tolerance)
        );
        CHECK_THAT(
            relative_matrix_difference(
                inverse_mapping_jacobian_difference, variation.mapping.inverse_mapping_jacobian_variation
            ),
            Catch::Matchers::WithinAbs(0.0, relative_tolerance)
        );

        const double determinant_difference =
            (plus_context.mapping.mapping_determinate - minus_context.mapping.mapping_determinate) /
            (2.0 * kelvin_difference_step);
        const double weight_difference = (plus_context.quadrature.physical_integration_weight -
                                             minus_context.quadrature.physical_integration_weight) /
                                         (2.0 * kelvin_difference_step);

        check_centered_difference(
            variation.mapping.mapping_determinate_variation, determinant_difference,
            plus_context.mapping.mapping_determinate, minus_context.mapping.mapping_determinate, kelvin_difference_step,
            relative_tolerance
        );
        check_centered_difference(
            variation.weight_variation, weight_difference, plus_context.quadrature.physical_integration_weight,
            minus_context.quadrature.physical_integration_weight, kelvin_difference_step, relative_tolerance
        );
    }

    for (int boundary_element = 0; boundary_element < mesh.GetNBE(); ++boundary_element) {
        mfem::FaceElementTransformations* face_transformation = mesh.GetBdrFaceTransformations(boundary_element);
        REQUIRE(face_transformation != nullptr);

        const mfem::IntegrationRule& integration_rule = mfem::IntRules.Get(face_transformation->GetGeometryType(), 4);

        for (int q = 0; q < integration_rule.GetNPoints(); ++q) {
            const mfem::IntegrationPoint& integration_point = integration_rule.IntPoint(q);

            mapping::FaceMappingContext base_context;
            mapping::FaceMappingContext plus_context;
            mapping::FaceMappingContext minus_context;
            mapping::FaceMappingVariation variation;

            REQUIRE(
                mapper.evaluate_face(
                    base_data.get(), *face_transformation, mapping::FaceElementSide::ELEMENT_1, integration_point,
                    workspace, base_context
                ) == mapping::MappingStatus::VALID
            );
            REQUIRE(
                mapper.evaluate_face_variation(
                    base_data.get(), direction, *face_transformation, mapping::FaceElementSide::ELEMENT_1,
                    integration_point, base_context, workspace, variation
                ) == mapping::MappingStatus::VALID
            );
            REQUIRE(
                mapper.evaluate_face(
                    plus_data.get(), *face_transformation, mapping::FaceElementSide::ELEMENT_1, integration_point,
                    workspace, plus_context
                ) == mapping::MappingStatus::VALID
            );
            REQUIRE(
                mapper.evaluate_face(
                    minus_data.get(), *face_transformation, mapping::FaceElementSide::ELEMENT_1, integration_point,
                    workspace, minus_context
                ) == mapping::MappingStatus::VALID
            );

            CAPTURE(boundary_element, q);

            mfem::Vector normal_difference =
                vector_central_difference(plus_context.quadrature.normal, minus_context.quadrature.normal);

            REQUIRE(normal_difference.Size() == variation.physical_normal_variation.Size());
            for (int i = 0; i < normal_difference.Size(); ++i) {
                check_centered_difference(
                    variation.physical_normal_variation(i), normal_difference(i), plus_context.quadrature.normal(i),
                    minus_context.quadrature.normal(i), kelvin_difference_step, relative_tolerance
                );
            }

            const double surface_weight_difference =
                (plus_context.physical_surface_weight - minus_context.physical_surface_weight) /
                (2.0 * kelvin_difference_step);
            const double normal_scale_difference =
                (plus_context.quadrature.v_dot_n_scale - minus_context.quadrature.v_dot_n_scale) /
                (2.0 * kelvin_difference_step);

            check_centered_difference(
                variation.physical_surface_weight_variation, surface_weight_difference,
                plus_context.physical_surface_weight, minus_context.physical_surface_weight, kelvin_difference_step,
                relative_tolerance
            );
            check_centered_difference(
                variation.normal_flux_scale_variation, normal_scale_difference, plus_context.quadrature.v_dot_n_scale,
                minus_context.quadrature.v_dot_n_scale, kelvin_difference_step, relative_tolerance
            );
        }
    }
}
TEST_CASE(
    "Stateless Domain Mapper Produces Consistent Two Sided Interface Geometry",
    tags::integration & tags::mapping
) {
    constexpr double r_star = 1.0;
    constexpr double r_infinity = 4.0;
    constexpr double interface_tolerance = 2.0e-7;

    stroid::config::MeshConfig mesh_config;
    mesh_config.refinement_levels = 0;
    mesh_config.order = 2;
    mesh_config.include_external_domain = true;
    mesh_config.r_core = 0.25;
    mesh_config.r_star = r_star;
    mesh_config.r_infinity = r_infinity;
    mesh_config.flattening = 0.0;
    mesh_config.optimization_methods = stroid::config::OptimizationMethods{false, true};

    stroid::StroidMesh stroid_mesh = stroid::GenerateMesh(mesh_config);
    mfem::Mesh& mesh = *stroid_mesh.mesh;
    REQUIRE(stroid_mesh.exterior_coordinate != nullptr);
    REQUIRE(stroid_mesh.exterior_coordinate->space != nullptr);
    REQUIRE(stroid_mesh.exterior_coordinate->values != nullptr);

    mfem::FiniteElementSpace& compactification_space = *stroid_mesh.exterior_coordinate->space;
    mfem::GridFunction& compactification_coordinate = *stroid_mesh.exterior_coordinate->values;

    mfem::H1_FECollection displacement_collection(2, dimension);
    mfem::FiniteElementSpace displacement_space(&mesh, &displacement_collection, dimension, mfem::Ordering::byVDIM);

    auto displacement_function = [](const mfem::Vector& position, mfem::Vector& value) {
        value.SetSize(3);
        value(0) = 0.08 * position(0) + 0.02 * position(1);
        value(1) = -0.03 * position(0) - 0.02 * position(1);
        value(2) = 0.04 * position(2);
    };

    auto direction_function = [](const mfem::Vector& position, mfem::Vector& value) {
        value.SetSize(3);
        value(0) = 0.02 * position(0);
        value(1) = -0.01 * position(1) + 0.005 * position(2);
        value(2) = 0.015 * position(2);
    };

    mfem::VectorFunctionCoefficient displacement_coefficient(dimension, displacement_function);
    mfem::VectorFunctionCoefficient direction_coefficient(dimension, direction_function);

    mfem::GridFunction displacement(&displacement_space);
    mfem::GridFunction direction(&displacement_space);
    displacement.ProjectCoefficient(displacement_coefficient);
    direction.ProjectCoefficient(direction_coefficient);

    std::unique_ptr<const mapping::compactification::ExteriorDomainMap> exterior_map =
        std::make_unique<mapping::compactification::KelvinCompactification>(
            user_options::KelvinCompactificationOptions{.r_star_ref = r_star, .r_inf_ref = r_infinity}
        );

    mapping::DomainMapper mapper({.dimension = dimension, .vacuum_element_attribute_ID = 3}, std::move(exterior_map));
    mapping::DMWorkspace workspace(dimension);

    int core_envelope_faces = 0;
    int stellar_vacuum_faces = 0;

    for (int face_id = 0; face_id < mesh.GetNumFaces(); ++face_id) {
        mfem::FaceElementTransformations* transformation = mesh.GetFaceElementTransformations(face_id);
        if (transformation == nullptr || transformation->Elem1 == nullptr || transformation->Elem2 == nullptr) {
            continue;
        }

        const int attribute_1 = transformation->Elem1->Attribute;
        const int attribute_2 = transformation->Elem2->Attribute;
        const bool core_envelope_interface =
            (attribute_1 == 1 && attribute_2 == 2) || (attribute_1 == 2 && attribute_2 == 1);
        const bool stellar_vacuum_interface = (attribute_1 == 3) != (attribute_2 == 3);

        if (!core_envelope_interface && !stellar_vacuum_interface) {
            continue;
        }

        if (core_envelope_interface) {
            ++core_envelope_faces;
        }
        if (stellar_vacuum_interface) {
            ++stellar_vacuum_faces;
        }

        const int element_1 = transformation->Elem1->ElementNo;
        const int element_2 = transformation->Elem2->ElementNo;

        mfem::Array<int> vdofs_1;
        mfem::Array<int> vdofs_2;
        mfem::Vector displacement_dofs_1;
        mfem::Vector displacement_dofs_2;
        mfem::Vector direction_dofs_1;
        mfem::Vector direction_dofs_2;
        mfem::Array<int> compactification_dof_indices_1;
        mfem::Array<int> compactification_dof_indices_2;
        mfem::Vector compactification_dofs_1;
        mfem::Vector compactification_dofs_2;

        displacement_space.GetElementVDofs(element_1, vdofs_1);
        displacement_space.GetElementVDofs(element_2, vdofs_2);
        displacement.GetSubVector(vdofs_1, displacement_dofs_1);
        displacement.GetSubVector(vdofs_2, displacement_dofs_2);
        direction.GetSubVector(vdofs_1, direction_dofs_1);
        direction.GetSubVector(vdofs_2, direction_dofs_2);
        compactification_space.GetElementDofs(element_1, compactification_dof_indices_1);
        compactification_space.GetElementDofs(element_2, compactification_dof_indices_2);
        compactification_coordinate.GetSubVector(compactification_dof_indices_1, compactification_dofs_1);
        compactification_coordinate.GetSubVector(compactification_dof_indices_2, compactification_dofs_2);

        const mapping::ElementDisplacementData displacement_1 = mapping::element_displacement_data_from_element_vdofs(
            *displacement_space.GetFE(element_1), displacement_dofs_1
        );
        const mapping::ElementDisplacementData displacement_2 = mapping::element_displacement_data_from_element_vdofs(
            *displacement_space.GetFE(element_2), displacement_dofs_2
        );
        const mapping::ElementDisplacementData direction_1 = mapping::element_displacement_data_from_element_vdofs(
            *displacement_space.GetFE(element_1), direction_dofs_1
        );
        const mapping::ElementDisplacementData direction_2 = mapping::element_displacement_data_from_element_vdofs(
            *displacement_space.GetFE(element_2), direction_dofs_2
        );

        const ElementMappingDataOwner element_data_1(
            displacement_1, *compactification_space.GetFE(element_1), compactification_dofs_1
        );
        const ElementMappingDataOwner element_data_2(
            displacement_2, *compactification_space.GetFE(element_2), compactification_dofs_2
        );

        const mfem::IntegrationRule& integration_rule = mfem::IntRules.Get(transformation->GetGeometryType(), 6);

        for (int q = 0; q < integration_rule.GetNPoints(); ++q) {
            const mfem::IntegrationPoint& integration_point = integration_rule.IntPoint(q);

            mapping::FaceMappingContext context_1;
            mapping::FaceMappingContext context_2;
            mapping::FaceMappingVariation variation_1;
            mapping::FaceMappingVariation variation_2;

            auto status_1 = mapper.evaluate_face(
                element_data_1.get(), *transformation, mapping::FaceElementSide::ELEMENT_1, integration_point,
                workspace, context_1
            );
            auto status_2 = mapper.evaluate_face(
                element_data_2.get(), *transformation, mapping::FaceElementSide::ELEMENT_2, integration_point,
                workspace, context_2
            );
            INFO("Element 1 ID = " << transformation->Elem1No);
            INFO("Element 2 ID = " << transformation->Elem2No);
            INFO("Element 1 attribute = " << transformation->Elem1->Attribute);
            INFO("Element 2 attribute = " << transformation->Elem2->Attribute);
            INFO("Element 2 status = " << static_cast<int>(status_2));

            mfem::Vector reference_position(dimension);
            transformation->Elem2->Transform(transformation->Elem2->GetIntPoint(), reference_position);
            INFO(
                std::format(
                    "Reference position = <{},{},{}>", reference_position(0), reference_position(1),
                    reference_position(2)
                )
            );
            INFO(std::format("Reference radius = {}", reference_position.Norml2()));
            INFO(std::format("Reference radius minus r_star = {}", reference_position.Norml2() - r_star));
            REQUIRE(status_1 == mapping::MappingStatus::VALID);
            REQUIRE(status_2 == mapping::MappingStatus::VALID);
            REQUIRE(
                mapper.evaluate_face_variation(
                    element_data_1.get(), direction_1, *transformation, mapping::FaceElementSide::ELEMENT_1,
                    integration_point, context_1, workspace, variation_1
                ) == mapping::MappingStatus::VALID
            );
            REQUIRE(
                mapper.evaluate_face_variation(
                    element_data_2.get(), direction_2, *transformation, mapping::FaceElementSide::ELEMENT_2,
                    integration_point, context_2, workspace, variation_2
                ) == mapping::MappingStatus::VALID
            );

            CAPTURE(face_id, q, attribute_1, attribute_2);
            check_vector(context_1.mapping.physical_position, context_2.mapping.physical_position, interface_tolerance);
            check_scalar_relative(
                context_1.physical_surface_weight, context_2.physical_surface_weight, interface_tolerance
            );

            mfem::Vector normal_sum(context_1.quadrature.normal);
            normal_sum += context_2.quadrature.normal;
            CHECK(normal_sum.Norml2() < interface_tolerance);

            check_vector(
                variation_1.mapping.physical_position_variation, variation_2.mapping.physical_position_variation,
                interface_tolerance
            );
            check_scalar_relative(
                variation_1.physical_surface_weight_variation, variation_2.physical_surface_weight_variation,
                interface_tolerance
            );

            mfem::Vector normal_variation_sum(variation_1.physical_normal_variation);
            normal_variation_sum += variation_2.physical_normal_variation;
            CHECK(normal_variation_sum.Norml2() < interface_tolerance);

            mfem::Vector physical_flux = make_vector(0.7, -0.4, 0.9);
            mfem::Vector reference_flux_1;
            mfem::Vector reference_flux_2;

            // Contravariant pullback det(A) A^-1 F of one physical flux from each side of the interface.
            reference_flux_1.SetSize(3);
            reference_flux_2.SetSize(3);
            context_1.mapping.inverse_mapping_jacobian.Mult(physical_flux, reference_flux_1);
            reference_flux_1 *= context_1.mapping.mapping_determinate;
            context_2.mapping.inverse_mapping_jacobian.Mult(physical_flux, reference_flux_2);
            reference_flux_2 *= context_2.mapping.mapping_determinate;

            const double flux_1 = (reference_flux_1 * context_1.reference_normal) * context_1.reference_surface_weight;
            const double flux_2 = (reference_flux_2 * context_2.reference_normal) * context_2.reference_surface_weight;

            CHECK_THAT(flux_1 + flux_2, WithinAbs(0.0, interface_tolerance));
        }
    }

    INFO("Core-envelope interface faces = " << core_envelope_faces);
    INFO("Stellar-vacuum interface faces = " << stellar_vacuum_faces);

    REQUIRE(core_envelope_faces > 0);
    REQUIRE(stellar_vacuum_faces > 0);
}

TEST_CASE(
    "Mapped Hcurl Fields Preserve Covariant Piola Identities",
    tags::unit & tags::mapping
) {
    constexpr double curl_tolerance = 2.0e-12;

    mapping::MappingPointContext context;
    context.mapping_jacobian.SetSize(3);
    context.mapping_jacobian(0, 0) = 1.20;
    context.mapping_jacobian(0, 1) = 0.15;
    context.mapping_jacobian(0, 2) = -0.05;
    context.mapping_jacobian(1, 0) = -0.08;
    context.mapping_jacobian(1, 1) = 0.95;
    context.mapping_jacobian(1, 2) = 0.12;
    context.mapping_jacobian(2, 0) = 0.04;
    context.mapping_jacobian(2, 1) = -0.10;
    context.mapping_jacobian(2, 2) = 1.10;

    context.mapping_determinate = context.mapping_jacobian.Det();
    REQUIRE(context.mapping_determinate > 0.0);

    context.inverse_mapping_jacobian.SetSize(3);
    mfem::CalcInverse(context.mapping_jacobian, context.inverse_mapping_jacobian);

    mfem::Vector reference_field = make_vector(0.7, -0.4, 1.1);
    mfem::Vector reference_test_field = make_vector(-0.2, 0.9, 0.5);

    mfem::DenseMatrix reference_gradient(3);
    reference_gradient(0, 0) = 0.20;
    reference_gradient(0, 1) = -0.10;
    reference_gradient(0, 2) = 0.04;
    reference_gradient(1, 0) = 0.03;
    reference_gradient(1, 1) = 0.15;
    reference_gradient(1, 2) = -0.08;
    reference_gradient(2, 0) = -0.05;
    reference_gradient(2, 1) = 0.02;
    reference_gradient(2, 2) = 0.11;

    mfem::Vector reference_curl = matrix_curl(reference_gradient);

    mfem::Vector physical_field;
    mfem::Vector recovered_field;
    mfem::Vector physical_test_field;
    mfem::Vector physical_curl;
    mfem::Vector recovered_curl;

    mapping::map_hcurl_field_to_physical(context, reference_field, physical_field);
    mapping::map_physical_field_to_hcurl_reference(context, physical_field, recovered_field);
    mapping::map_hcurl_field_to_physical(context, reference_test_field, physical_test_field);
    mapping::map_hcurl_curl_to_physical(context, reference_curl, physical_curl);
    mapping::map_physical_curl_to_hcurl_reference(context, physical_curl, recovered_curl);

    check_vector(recovered_field, reference_field, curl_tolerance);
    check_vector(recovered_curl, reference_curl, curl_tolerance);

    mfem::DenseMatrix temporary(3);
    mfem::DenseMatrix physical_gradient(3);
    mfem::MultAtB(context.inverse_mapping_jacobian, reference_gradient, temporary);
    mfem::Mult(temporary, context.inverse_mapping_jacobian, physical_gradient);

    mfem::Vector directly_computed_physical_curl = matrix_curl(physical_gradient);
    check_vector(directly_computed_physical_curl, physical_curl, curl_tolerance);

    mfem::DenseMatrix mass_tensor;
    mfem::DenseMatrix curl_tensor;
    mapping::compute_scalar_diffusion_tensor(context, mass_tensor);
    curl_tensor.SetSize(3);
    mfem::MultAtB(context.mapping_jacobian, context.mapping_jacobian, curl_tensor);
    curl_tensor *= 1.0 / context.mapping_determinate;

    mfem::Vector mass_action(3);
    mass_tensor.Mult(reference_test_field, mass_action);

    const double physical_mass_inner_product = context.mapping_determinate * (physical_field * physical_test_field);
    const double reference_mass_inner_product = reference_field * mass_action;

    CHECK_THAT(physical_mass_inner_product, WithinAbs(reference_mass_inner_product, curl_tolerance));

    mfem::Vector reference_test_curl = make_vector(-0.3, 0.6, 0.2);
    mfem::Vector physical_test_curl;
    mfem::Vector curl_action(3);

    mapping::map_hcurl_curl_to_physical(context, reference_test_curl, physical_test_curl);
    curl_tensor.Mult(reference_test_curl, curl_action);

    const double physical_curl_inner_product = context.mapping_determinate * (physical_curl * physical_test_curl);
    const double reference_curl_inner_product = reference_curl * curl_action;

    CHECK_THAT(physical_curl_inner_product, WithinAbs(reference_curl_inner_product, curl_tolerance));
}

TEST_CASE(
    "Element Displacement Data Matches MFEM GridFunction Evaluation",
    tags::unit & tags::mapping
) {
    constexpr int displacement_order = 2;
    constexpr int quadrature_order = 6;
    constexpr double value_tolerance = 5.0e-13;
    constexpr double gradient_tolerance = 2.0e-12;
    constexpr double mapping_tolerance = 3.0e-12;

    auto check_space_ordering = [](const mfem::Ordering::Type space_ordering) {
        mfem::Mesh mesh = mfem::Mesh::MakeCartesian3D(2, 1, 1, mfem::Element::HEXAHEDRON, 2.0, 1.5, 1.25);

        mfem::H1_FECollection displacement_collection(displacement_order, dimension);

        mfem::FiniteElementSpace displacement_space(&mesh, &displacement_collection, dimension, space_ordering);

        mfem::GridFunction displacement(&displacement_space);

        auto displacement_function = [](const mfem::Vector& position, mfem::Vector& value) {
            const double x = position(0);
            const double y = position(1);
            const double z = position(2);

            value.SetSize(dimension);

            value(0) = 0.17 + 0.11 * x - 0.07 * y + 0.03 * y * z;

            value(1) = -0.23 + 0.05 * y + 0.09 * z + 0.02 * x * z;

            value(2) = 0.31 - 0.04 * x + 0.08 * z - 0.015 * x * y;
        };

        mfem::VectorFunctionCoefficient displacement_coefficient(dimension, displacement_function);

        displacement.ProjectCoefficient(displacement_coefficient);

        mapping::DomainMapper mapper(
            {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
        );

        mapping::DMWorkspace workspace(dimension);

        REQUIRE(displacement_space.GetOrdering() == space_ordering);
        REQUIRE(displacement.VectorDim() == dimension);

        for (int element_id = 0; element_id < mesh.GetNE(); ++element_id) {
            mfem::ElementTransformation* transformation = mesh.GetElementTransformation(element_id);

            REQUIRE(transformation != nullptr);

            const mfem::FiniteElement& displacement_element = *displacement_space.GetFE(element_id);

            mfem::Array<int> element_vdofs;

            mfem::DofTransformation* dof_transformation = displacement_space.GetElementVDofs(element_id, element_vdofs);

            mfem::Vector element_displacement;

            displacement.GetSubVector(element_vdofs, element_displacement);

            if (dof_transformation != nullptr) {
                dof_transformation->InvTransformPrimal(element_displacement);
            }

            const mapping::ElementDisplacementData displacement_data =
                mapping::element_displacement_data_from_element_vdofs(displacement_element, element_displacement);


            const ElementMappingDataOwner element_data(displacement_data);

            const mfem::IntegrationRule& integration_rule =
                mfem::IntRules.Get(transformation->GetGeometryType(), quadrature_order);

            mfem::Vector shape(displacement_element.GetDof());

            mfem::DenseMatrix physical_dshape(displacement_element.GetDof(), dimension);

            mfem::Vector computed_value(dimension);
            mfem::Vector expected_value(dimension);

            mfem::DenseMatrix computed_gradient(dimension, dimension);

            mfem::DenseMatrix expected_gradient(dimension, dimension);

            for (int q = 0; q < integration_rule.GetNPoints(); ++q) {
                const mfem::IntegrationPoint& integration_point = integration_rule.IntPoint(q);

                CAPTURE(static_cast<int>(space_ordering), element_id, q);

                transformation->SetIntPoint(&integration_point);

                displacement_element.CalcShape(integration_point, shape);

                displacement_element.CalcPhysDShape(*transformation, physical_dshape);

                displacement_data.get_dof_matrix().MultTranspose(shape, computed_value);

                mfem::MultAtB(displacement_data.get_dof_matrix(), physical_dshape, computed_gradient);

                /*
                 * Use MFEM's native evaluation as the authoritative
                 * interpretation of the GridFunction.
                 */
                transformation->SetIntPoint(&integration_point);

                displacement.GetVectorValue(*transformation, integration_point, expected_value);

                transformation->SetIntPoint(&integration_point);

                displacement.GetVectorGradient(*transformation, expected_gradient);

                check_vector(computed_value, expected_value, value_tolerance);

                check_matrix(computed_gradient, expected_gradient, gradient_tolerance);

                /*
                 * Also exercise the complete stateless-mapper path.
                 * The Cartesian elements are stellar-domain elements,
                 * so the mapper should produce x + u and I + grad(u).
                 */
                mapping::MappingPointContext context;

                const mapping::MappingStatus status =
                    mapper.evaluate_point(element_data.get(), *transformation, integration_point, workspace, context);

                REQUIRE(status == serif::mapping::MappingStatus::VALID);

                REQUIRE_FALSE(context.compactified);

                mfem::Vector expected_displaced_position(context.reference_position);

                expected_displaced_position += expected_value;

                mfem::DenseMatrix expected_displacement_jacobian(expected_gradient);

                for (int d = 0; d < dimension; ++d) {
                    expected_displacement_jacobian(d, d) += 1.0;
                }

                check_vector(context.displaced_position, expected_displaced_position, mapping_tolerance);

                check_vector(context.physical_position, expected_displaced_position, mapping_tolerance);

                check_matrix(context.displacement_jacobian, expected_displacement_jacobian, mapping_tolerance);

                check_matrix(context.mapping_jacobian, expected_displacement_jacobian, mapping_tolerance);

                const double expected_determinant = expected_displacement_jacobian.Det();

                CHECK_THAT(
                    context.mapping_determinate, Catch::Matchers::WithinAbs(expected_determinant, mapping_tolerance)
                );

                REQUIRE(context.mapping_determinate > 0.0);
            }
        }
    };

    SECTION("Global finite-element-space ordering is byNODES") {
        check_space_ordering(mfem::Ordering::byNODES);
    }

    SECTION("Global finite-element-space ordering is byVDIM") {
        check_space_ordering(mfem::Ordering::byVDIM);
    }
}

TEST_CASE(
    "Grid Function Mapping Evaluator Refreshes A Cached Displacement Element",
    tags::unit & tags::mapping
) {
    mfem::Mesh mesh = mfem::Mesh::MakeCartesian3D(1, 1, 1, mfem::Element::HEXAHEDRON, 2.0, 3.0, 4.0);

    mfem::H1_FECollection displacement_collection(1, dimension);
    mfem::H1_FECollection compactification_collection(1, dimension);
    mfem::FiniteElementSpace displacement_space(&mesh, &displacement_collection, dimension, mfem::Ordering::byVDIM);
    mfem::FiniteElementSpace compactification_space(&mesh, &compactification_collection);

    mfem::GridFunction displacement(&displacement_space);
    mfem::GridFunction compactification_coordinate(&compactification_space);
    displacement = 0.0;
    compactification_coordinate = 0.0;

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );
    mapping::GridFunctionMappingEvaluator evaluator(mapper, displacement, compactification_coordinate);

    /* Refreshing before the first evaluation is a validated no-op. */
    evaluator.refresh();

    mfem::ElementTransformation* transformation = mesh.GetElementTransformation(0);
    REQUIRE(transformation != nullptr);

    mfem::IntegrationPoint integration_point;
    integration_point.Set3(0.31, 0.43, 0.57);

    mapping::MappingPointContext initial_context;
    REQUIRE(
        evaluator.evaluate_point(*transformation, integration_point, initial_context) == mapping::MappingStatus::VALID
    );

    mfem::Vector reference_position(dimension);
    transformation->Transform(integration_point, reference_position);
    check_vector(initial_context.physical_position, reference_position);

    mfem::Vector displacement_offset = make_vector(0.17, -0.09, 0.045);
    mfem::VectorFunctionCoefficient displacement_coefficient(
        dimension, [&displacement_offset](const mfem::Vector&, mfem::Vector& value) {
            value = displacement_offset;
        }
    );
    displacement.ProjectCoefficient(displacement_coefficient);

    /*
     * The same element ID is still cached. refresh must eagerly rebuild it from
     * the mutated grid function rather than retaining the old element data.
     */
    evaluator.refresh();

    mapping::MappingPointContext refreshed_context;
    REQUIRE(
        evaluator.evaluate_point(*transformation, integration_point, refreshed_context) == mapping::MappingStatus::VALID
    );

    mfem::Vector expected_position(reference_position);
    expected_position += displacement_offset;

    check_vector(refreshed_context.physical_position, expected_position);
    check_vector(refreshed_context.displaced_position, expected_position);
}

TEST_CASE(
    "Grid Function Mapping Evaluator Invalidates A Cached Compactification Element",
    tags::unit & tags::mapping
) {
    mfem::Mesh mesh = mfem::Mesh::MakeCartesian3D(1, 1, 1, mfem::Element::HEXAHEDRON, 2.0, 3.0, 4.0);
    mesh.GetElement(0)->SetAttribute(3);
    mesh.SetAttributes();

    mfem::H1_FECollection displacement_collection(1, dimension);
    mfem::H1_FECollection compactification_collection(1, dimension);
    mfem::FiniteElementSpace displacement_space(&mesh, &displacement_collection, dimension, mfem::Ordering::byVDIM);
    mfem::FiniteElementSpace compactification_space(&mesh, &compactification_collection);

    mfem::GridFunction displacement(&displacement_space);
    mfem::GridFunction compactification_coordinate(&compactification_space);
    displacement = 0.0;
    compactification_coordinate = 0.20;

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );
    mapping::GridFunctionMappingEvaluator evaluator(mapper, displacement, compactification_coordinate);

    mfem::ElementTransformation* transformation = mesh.GetElementTransformation(0);
    REQUIRE(transformation != nullptr);

    mfem::IntegrationPoint integration_point;
    integration_point.Set3(0.29, 0.37, 0.61);

    mapping::MappingPointContext initial_context;
    REQUIRE(
        evaluator.evaluate_point(*transformation, integration_point, initial_context) == mapping::MappingStatus::VALID
    );
    REQUIRE(initial_context.compactified);

    compactification_coordinate = 0.40;

    /* Lazy invalidation is idempotent and reloads on the next evaluation. */
    evaluator.invalidate_cache();
    evaluator.invalidate_cache();

    mapping::MappingPointContext refreshed_context;
    REQUIRE(
        evaluator.evaluate_point(*transformation, integration_point, refreshed_context) == mapping::MappingStatus::VALID
    );
    REQUIRE(refreshed_context.compactified);

    constexpr double coordinate = 0.40;
    constexpr double computational_radius = 1.0 + 3.0 * coordinate;
    constexpr double expected_scale = 1.0 / (computational_radius * (1.0 - coordinate));

    mfem::Vector expected_position(refreshed_context.reference_position);
    expected_position *= expected_scale;

    check_vector(refreshed_context.physical_position, expected_position);
    CHECK(relative_vector_difference(initial_context.physical_position, refreshed_context.physical_position) > 1.0e-3);
}

TEST_CASE(
    "Grid Function Mapping Evaluator Rejects Invalid Fields And Rebinding",
    tags::unit & tags::mapping
) {
    mfem::Mesh mesh = mfem::Mesh::MakeCartesian3D(1, 1, 1, mfem::Element::HEXAHEDRON, 2.0, 3.0, 4.0);
    mfem::Mesh other_mesh = mfem::Mesh::MakeCartesian3D(1, 1, 1, mfem::Element::HEXAHEDRON, 2.0, 3.0, 4.0);

    mfem::H1_FECollection collection(1, dimension);
    mfem::FiniteElementSpace displacement_space(&mesh, &collection, dimension, mfem::Ordering::byVDIM);
    mfem::FiniteElementSpace alternate_displacement_space(&mesh, &collection, dimension, mfem::Ordering::byVDIM);
    mfem::FiniteElementSpace scalar_space(&mesh, &collection);
    mfem::FiniteElementSpace other_scalar_space(&other_mesh, &collection);
    mfem::FiniteElementSpace vector_compactification_space(&mesh, &collection, dimension, mfem::Ordering::byVDIM);

    mfem::GridFunction displacement(&displacement_space);
    mfem::GridFunction scalar_displacement(&scalar_space);
    mfem::GridFunction compactification_coordinate(&scalar_space);
    mfem::GridFunction other_compactification_coordinate(&other_scalar_space);
    mfem::GridFunction vector_compactification_coordinate(&vector_compactification_space);
    mfem::GridFunction detached;
    displacement = 0.0;
    compactification_coordinate = 0.0;

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );

    CHECK_THROWS_AS(
        (mapping::GridFunctionMappingEvaluator(mapper, detached, compactification_coordinate)), std::logic_error
    );
    CHECK_THROWS_AS((mapping::GridFunctionMappingEvaluator(mapper, displacement, detached)), std::logic_error);
    CHECK_THROWS_AS(
        (mapping::GridFunctionMappingEvaluator(mapper, scalar_displacement, compactification_coordinate)),
        std::logic_error
    );
    CHECK_THROWS_AS(
        (mapping::GridFunctionMappingEvaluator(mapper, displacement, vector_compactification_coordinate)),
        std::logic_error
    );
    CHECK_THROWS_AS(
        (mapping::GridFunctionMappingEvaluator(mapper, displacement, other_compactification_coordinate)),
        std::logic_error
    );

    mapping::GridFunctionMappingEvaluator evaluator(mapper, displacement, compactification_coordinate);

    mfem::ElementTransformation* transformation = mesh.GetElementTransformation(0);
    REQUIRE(transformation != nullptr);

    mfem::IntegrationPoint integration_point;
    integration_point.Set3(0.5, 0.5, 0.5);

    mapping::MappingPointContext context;
    REQUIRE(evaluator.evaluate_point(*transformation, integration_point, context) == mapping::MappingStatus::VALID);

    displacement.SetSpace(&alternate_displacement_space);

    CHECK_THROWS_AS(evaluator.refresh(), std::logic_error);
    CHECK_THROWS_AS(evaluator.evaluate_point(*transformation, integration_point, context), std::logic_error);
}

TEST_CASE(
    "Prepared Reference Geometry Matches Direct Mapping And Current Mesh Coordinates",
    tags::unit & tags::geometry_cache
) {
    SingleElementFixture fixture;
    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );
    mapping::DMWorkspace workspace(dimension);
    const auto& point = mfem::IntRules.Get(mfem::Geometry::CUBE, 5).IntPoint(0);
    auto zero = fixture.make_zero_element_dofs();
    mapping::ElementDisplacementData displacement(fixture.get_element(), zero);
    ElementMappingDataOwner data(displacement);
    auto compare = [&](const mfem::IntegrationPoint& ip) {
        auto& transformation = *fixture.mesh.GetElementTransformation(0);
        const auto geometry = mapping::sample_fixed_mesh_geometry(transformation, ip);
        REQUIRE(geometry);
        mapping::VolumeMappingContext prepared, direct;
        REQUIRE(mapper.evaluate_volume(data.get(), transformation, ip, workspace, prepared, nullptr, &*geometry) == mapping::MappingStatus::VALID);
        REQUIRE(mapper.evaluate_volume(data.get(), transformation, ip, workspace, direct) == mapping::MappingStatus::VALID);
        check_vector(prepared.mapping.physical_position, direct.mapping.physical_position);
        check_matrix(prepared.mapping.mapping_jacobian, direct.mapping.mapping_jacobian);
        check_matrix(prepared.quadrature.J_element_to_physical_inverse, direct.quadrature.J_element_to_physical_inverse);
        CHECK_THAT(prepared.quadrature.physical_integration_weight, WithinAbs(direct.quadrature.physical_integration_weight, 1.0e-12));
        return prepared.quadrature.physical_integration_weight;
    };
    const double original_weight = compare(point);
    compare(point);
    auto changed_point = point;
    changed_point.x += 0.013;
    const double original_point_weight = compare(changed_point);
    changed_point.weight *= 0.5;
    CHECK_THAT(compare(changed_point), WithinAbs(0.5 * original_point_weight, 1.0e-12));
    fixture.mesh.Transform([](const mfem::Vector& x, mfem::Vector& y) {
        y = x;
        y(0) *= 1.2;
        y(1) += 0.17;
    });
    CHECK(compare(point) > original_weight);
}

TEST_CASE(
    "Prepared Reference Geometry Does Not Freeze Displacement Or Compactification State",
    tags::unit & tags::geometry_cache
) {
    SingleElementFixture fixture;
    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
    );
    mapping::DMWorkspace workspace(dimension);
    const auto& point = mfem::IntRules.Get(mfem::Geometry::CUBE, 3).IntPoint(0);
    auto& transformation = *fixture.mesh.GetElementTransformation(0);
    const auto geometry = mapping::sample_fixed_mesh_geometry(transformation, point);
    REQUIRE(geometry);
    for (double scale : {0.0, 0.02, -0.01}) {
        auto dofs = fixture.make_zero_element_dofs();
        dofs = scale;
        mapping::ElementDisplacementData displacement(fixture.get_element(), dofs);
        auto compactification_dofs = make_constant_compactification_dofs(fixture.get_element(), 0.1 + scale);
        mapping::ElementCompactificationData compactification(fixture.get_element(), compactification_dofs);
        mapping::ElementMappingData data{displacement, compactification};
        for (int attribute : {1, 3}) {
            transformation.Attribute = attribute;
            mapping::VolumeMappingContext prepared, direct;
            const auto prepared_status = mapper.evaluate_volume(data, transformation, point, workspace, prepared, nullptr, &*geometry);
            const auto direct_status = mapper.evaluate_volume(data, transformation, point, workspace, direct);
            REQUIRE(prepared_status == direct_status);
            if (prepared_status == mapping::MappingStatus::VALID) {
                check_vector(prepared.mapping.physical_position, direct.mapping.physical_position);
                check_matrix(prepared.mapping.mapping_jacobian, direct.mapping.mapping_jacobian);
                check_matrix(prepared.quadrature.J_element_to_physical_inverse, direct.quadrature.J_element_to_physical_inverse);
            }
        }
    }
}

TEST_CASE(
    "Reference Geometry Supports Detached Transformations And Rejects Singular Meshes",
    tags::unit & tags::geometry_cache
) {
    SingleElementFixture fixture;
    auto& transformation = *fixture.mesh.GetElementTransformation(0);
    const auto& point = mfem::IntRules.Get(mfem::Geometry::CUBE, 3).IntPoint(0);
    const auto attached = mapping::sample_fixed_mesh_geometry(transformation, point);
    REQUIRE(attached);
    transformation.mesh = nullptr;
    const auto detached = mapping::sample_fixed_mesh_geometry(transformation, point);
    REQUIRE(detached);
    check_vector(detached->position, attached->position);
    check_matrix(detached->inverse_jacobian, attached->inverse_jacobian);
    fixture.mesh.Transform([](const mfem::Vector& x, mfem::Vector& y) {
        y = x;
        y(0) = 0.0;
    });
    const auto singular = mapping::sample_fixed_mesh_geometry(*fixture.mesh.GetElementTransformation(0), point);
    REQUIRE_FALSE(singular);
    CHECK(singular.error() == mapping::MappingStatus::NON_POSITIVE_DETERMINATE);
}

TEST_CASE(
    "Supplied compactification samples preserve volume variations and reject invalid exterior inputs",
    tags::unit & tags::mapping
) {
    for (const int attribute : {1, 3}) {
        QuadraticElementFixture fixture;
        fixture.mesh.GetElement(0)->SetAttribute(attribute);
        fixture.mesh.SetAttributes();
        const mfem::FiniteElement& element = fixture.get_element();
        auto& transformation = *fixture.mesh.GetElementTransformation(0);
        QuadraticMappingData fields(element, transformation);
        const mapping::ElementDisplacementData base(element, fields.base_dofs);
        const mapping::ElementDisplacementData direction(element, fields.direction_dofs);
        const mfem::Vector compactification_dofs = make_compactification_element_dofs(
            element, transformation, [](const mfem::Vector& position) {
                return 0.1 + position[0] / 8.0;
            }
        );
        const ElementMappingDataOwner data(base, element, compactification_dofs);
        mapping::DomainMapper mapper(
            {.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification()
        );
        mapping::DMWorkspace workspace(dimension);
        mapping::CompactificationPointData sample;
        mfem::DenseMatrix inverse_mesh_jacobian(dimension);
        for (const double coordinate : {0.25, 0.5, 0.75}) {
            CAPTURE(attribute, coordinate);
            mfem::IntegrationPoint point;
            point.Set3(coordinate, 0.4, 0.6);
            point.weight = 0.73;
            mapping::VolumeMappingContext context;
            mapping::VolumeMappingVariation expected, supplied;
            REQUIRE(mapper.evaluate_volume(data.get(), transformation, point, workspace, context) == mapping::MappingStatus::VALID);
            REQUIRE(serif::tests::evaluate_volume_variation(mapper,
                data.get(), direction, transformation, point, context, workspace, expected
            ) == mapping::MappingStatus::VALID);
            mfem::Mult(context.quadrature.J_element_to_physical_inverse, context.mapping.mapping_jacobian, inverse_mesh_jacobian);
            REQUIRE(mapper.evaluate_compactification_coordinate(
                data.get().compactification, transformation, point, workspace, sample, &inverse_mesh_jacobian
            ) == mapping::MappingStatus::VALID);
            REQUIRE(serif::tests::evaluate_volume_variation(mapper,
                data.get(), direction, transformation, point, context, workspace, supplied, &sample
            ) == mapping::MappingStatus::VALID);
            check_vector(supplied.mapping.displacement_variation, expected.mapping.displacement_variation, 0.0);
            check_vector(supplied.mapping.physical_position_variation, expected.mapping.physical_position_variation, 0.0);
            check_matrix(supplied.mapping.displacement_jacobian_variation, expected.mapping.displacement_jacobian_variation, 0.0);
            check_matrix(supplied.mapping.mapping_jacobian_variation, expected.mapping.mapping_jacobian_variation, 0.0);
            check_matrix(supplied.mapping.inverse_mapping_jacobian_variation, expected.mapping.inverse_mapping_jacobian_variation, 0.0);
            CHECK(supplied.mapping.mapping_determinate_variation == expected.mapping.mapping_determinate_variation);
            CHECK(supplied.weight_variation == expected.weight_variation);

            // A caller need not establish the transformation's current integration point first.
            mfem::IntegrationPoint other_point;
            other_point.Set3(0.9, 0.2, 0.1);
            transformation.SetIntPoint(&other_point);
            mapping::CompactificationPointData independently_sampled;
            REQUIRE(mapper.evaluate_compactification_coordinate(
                data.get().compactification, transformation, point, workspace, independently_sampled
            ) == mapping::MappingStatus::VALID);
            CHECK_THAT(independently_sampled.compactification_coordinate, WithinAbs(sample.compactification_coordinate, tolerance));
            check_vector(independently_sampled.coordinate_gradient, sample.coordinate_gradient);

            const mapping::CompactificationPointData valid_sample = sample;
            sample.coordinate_gradient.SetSize(2);
            CHECK(serif::tests::evaluate_volume_variation(mapper,
                data.get(), direction, transformation, point, context, workspace, supplied, &sample
            ) == (attribute == 3 ? mapping::MappingStatus::INVALID_DIMENSION : mapping::MappingStatus::VALID));
            sample = valid_sample;
            sample.coordinate_gradient[1] = std::numeric_limits<double>::infinity();
            CHECK(serif::tests::evaluate_volume_variation(mapper,
                data.get(), direction, transformation, point, context, workspace, supplied, &sample
            ) == (attribute == 3 ? mapping::MappingStatus::NON_FINITE_MAPPING_INPUT : mapping::MappingStatus::VALID));
            sample = valid_sample;
            sample.compactification_coordinate = std::numeric_limits<double>::quiet_NaN();
            CHECK(serif::tests::evaluate_volume_variation(mapper,
                data.get(), direction, transformation, point, context, workspace, supplied, &sample
            ) == (attribute == 3 ? mapping::MappingStatus::NON_FINITE_MAPPING_INPUT : mapping::MappingStatus::VALID));
            sample.compactification_coordinate = 2.0;
            CHECK(serif::tests::evaluate_volume_variation(mapper,
                data.get(), direction, transformation, point, context, workspace, supplied, &sample
            ) == (attribute == 3 ? mapping::MappingStatus::OUTSIDE_REFERENCE_DOMAIN : mapping::MappingStatus::VALID));
            sample = valid_sample;

            mfem::DenseMatrix invalid_inverse(2);
            CHECK(mapper.evaluate_compactification_coordinate(
                data.get().compactification, transformation, point, workspace, sample, &invalid_inverse
            ) == mapping::MappingStatus::INVALID_DIMENSION);
            invalid_inverse = inverse_mesh_jacobian;
            invalid_inverse(0, 0) = std::numeric_limits<double>::quiet_NaN();
            CHECK(mapper.evaluate_compactification_coordinate(
                data.get().compactification, transformation, point, workspace, sample, &invalid_inverse
            ) == mapping::MappingStatus::NON_FINITE_MAPPING_INPUT);
        }
    }
}

TEST_CASE(
    "Same-point fixed geometry preserves direct maps and variations on curved conforming and NC meshes",
    tags::unit & tags::mapping & tags::geometry_cache
) {
    for (const bool nonconforming : {false, true}) {
        mfem::Mesh mesh = mfem::Mesh::MakeCartesian3D(2, 1, 1, mfem::Element::HEXAHEDRON, 2.0, 1.0, 1.0);
        if (nonconforming) {
            mesh.EnsureNCMesh();
            mfem::Array<int> selected{0};
            mesh.GeneralRefinement(selected);
        }
        mesh.SetCurvature(6);
        mesh.Transform([](const mfem::Vector& x, mfem::Vector& y) {
            y = x;
            y[0] += 0.03 * x[0] * x[1];
            y[1] += 0.02 * x[1] * x[2];
            y[2] += 0.01 * x[0] * x[2];
        });
        mfem::H1_FECollection collection(3, dimension);
        mfem::FiniteElementSpace space(&mesh, &collection, dimension);
        if (nonconforming) { REQUIRE(space.GetTrueVSize() < space.GetVSize()); }
        mapping::DomainMapper mapper({.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification());
        mapping::DMWorkspace direct_workspace(dimension), supplied_workspace(dimension);
        for (int element_index = 0; element_index < mesh.GetNE(); ++element_index) {
            const auto& element = *space.GetFE(element_index);
            auto& transformation = *mesh.GetElementTransformation(element_index);
            auto base_values = make_function_element_dofs(element, transformation, evaluate_quadratic_displacement, mfem::Ordering::byNODES);
            auto direction_values = make_function_element_dofs(element, transformation, evaluate_quadratic_direction, mfem::Ordering::byNODES);
            const mapping::ElementDisplacementData base(element, base_values);
            const mapping::ElementDisplacementData direction(element, direction_values);
            const auto compactification_values = make_compactification_element_dofs(element, transformation,
                [](const mfem::Vector& x) { return 0.1 + 0.03 * x[0] + 0.02 * x[1]; });
            const ElementMappingDataOwner data(base, element, compactification_values);
            for (const int attribute : {1, 3}) {
                transformation.Attribute = attribute;
                for (const double coordinate : {0.2, 0.5, 0.8}) {
                    CAPTURE(nonconforming, element_index, attribute, coordinate);
                    mfem::IntegrationPoint point, other_point;
                    point.Set3(coordinate, 0.4, 0.6);
                    other_point.Set3(0.9, 0.1, 0.3);
                    mapping::MappingPointContext direct, supplied;
                    mapping::MappingPointVariation direct_variation, supplied_variation;
                    REQUIRE(mapper.evaluate_point(data.get(), transformation, point, direct_workspace, direct) == mapping::MappingStatus::VALID);
                    REQUIRE(mapper.evaluate_point_variation(data.get(), direction, transformation, point, direct, direct_workspace, direct_variation) == mapping::MappingStatus::VALID);

                    // Reset the transformation between calls: reuse is of the owned checked sample,
                    // not a dependency on MFEM's mutable last-point Jacobian cache.
                    const auto fixed = mapping::sample_fixed_mesh_geometry(transformation, point);
                    REQUIRE(fixed);
                    transformation.SetIntPoint(&other_point);
                    REQUIRE(mapper.evaluate_point(data.get(), transformation, point, supplied_workspace, supplied, *fixed) == mapping::MappingStatus::VALID);
                    transformation.SetIntPoint(&other_point);
                    REQUIRE(mapper.evaluate_point_variation(data.get(), direction, transformation, point, supplied, supplied_workspace, supplied_variation, fixed->inverse_jacobian) == mapping::MappingStatus::VALID);
                    check_point_context(supplied, direct, 0.0);
                    check_vector(supplied_variation.displacement_variation, direct_variation.displacement_variation, 0.0);
                    check_vector(supplied_variation.physical_position_variation, direct_variation.physical_position_variation, 0.0);
                    check_matrix(supplied_variation.displacement_jacobian_variation, direct_variation.displacement_jacobian_variation, 0.0);
                    check_matrix(supplied_variation.mapping_jacobian_variation, direct_variation.mapping_jacobian_variation, 0.0);
                    check_matrix(supplied_variation.inverse_mapping_jacobian_variation, direct_variation.inverse_mapping_jacobian_variation, 0.0);
                    CHECK(supplied_variation.mapping_determinate_variation == direct_variation.mapping_determinate_variation);
                }
            }
        }
    }
}

TEST_CASE(
    "Same-point reuse rejects malformed geometry and preserves physical mapping failures",
    tags::unit & tags::mapping & tags::geometry_cache
) {
    SingleElementFixture fixture;
    auto& transformation = *fixture.mesh.GetElementTransformation(0);
    const auto& point = mfem::IntRules.Get(mfem::Geometry::CUBE, 3).IntPoint(0);
    const auto fixed = mapping::sample_fixed_mesh_geometry(transformation, point);
    REQUIRE(fixed);
    mapping::DomainMapper mapper({.dimension = dimension, .vacuum_element_attribute_ID = 3}, make_kelvin_compactification());
    mapping::DMWorkspace workspace(dimension);
    auto zero = fixture.make_zero_element_dofs();
    const mapping::ElementDisplacementData base(fixture.get_element(), zero);
    const ElementMappingDataOwner data(base);
    mapping::MappingPointContext context;
    mapping::MappingPointVariation variation;
    auto invalid = *fixed;
    invalid.position.SetSize(2);
    CHECK(mapper.evaluate_point(data.get(), transformation, point, workspace, context, invalid) == mapping::MappingStatus::INVALID_DIMENSION);
    invalid = *fixed;
    invalid.inverse_jacobian(0, 0) = std::numeric_limits<double>::infinity();
    CHECK(mapper.evaluate_point(data.get(), transformation, point, workspace, context, invalid) == mapping::MappingStatus::NON_FINITE_MAPPING_INPUT);
    invalid = *fixed;
    invalid.weight = 0.0;
    CHECK(mapper.evaluate_point(data.get(), transformation, point, workspace, context, invalid) == mapping::MappingStatus::NON_POSITIVE_DETERMINATE);
    invalid = *fixed;
    invalid.jacobian = 0.0;
    CHECK(mapper.evaluate_point(data.get(), transformation, point, workspace, context, invalid) == mapping::MappingStatus::NON_POSITIVE_DETERMINATE);
    REQUIRE(mapper.evaluate_point(data.get(), transformation, point, workspace, context) == mapping::MappingStatus::VALID);
    mfem::DenseMatrix invalid_inverse(2);
    CHECK(mapper.evaluate_point_variation(data.get(), base, transformation, point, context, workspace, variation, invalid_inverse) == mapping::MappingStatus::INVALID_DIMENSION);
    invalid_inverse = fixed->inverse_jacobian;
    invalid_inverse(0, 0) = std::numeric_limits<double>::quiet_NaN();
    CHECK(mapper.evaluate_point_variation(data.get(), base, transformation, point, context, workspace, variation, invalid_inverse) == mapping::MappingStatus::NON_FINITE_MAPPING_INPUT);

    auto inverted_values = make_affine_element_dofs(fixture.get_element(), transformation,
        [] { auto gradient = make_identity_matrix(dimension); gradient *= -2.0; return gradient; }(),
        make_vector(0.0, 0.0, 0.0), mfem::Ordering::byNODES);
    const mapping::ElementDisplacementData inverted(fixture.get_element(), inverted_values);
    const ElementMappingDataOwner inverted_data(inverted);
    CHECK(mapper.evaluate_point(inverted_data.get(), transformation, point, workspace, context) == mapping::MappingStatus::NON_POSITIVE_DETERMINATE);
    CHECK(mapper.evaluate_point(inverted_data.get(), transformation, point, workspace, context, *fixed) == mapping::MappingStatus::NON_POSITIVE_DETERMINATE);
    transformation.Attribute = 3;
    const auto invalid_compactification = make_constant_compactification_dofs(fixture.get_element(), 2.0);
    const ElementMappingDataOwner invalid_exterior(base, fixture.get_element(), invalid_compactification);
    CHECK(mapper.evaluate_point(invalid_exterior.get(), transformation, point, workspace, context) == mapping::MappingStatus::OUTSIDE_REFERENCE_DOMAIN);
    CHECK(mapper.evaluate_point(invalid_exterior.get(), transformation, point, workspace, context, *fixed) == mapping::MappingStatus::OUTSIDE_REFERENCE_DOMAIN);
}