#include <algorithm>
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <limits>
#include <memory>
#include <mfem.hpp>
#include <stdexcept>
#include <stroid/stroid.h>

import mean_field;
import test_helpers;
using namespace mean_field;
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>(
            mapping::compactification::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.GetElement(),
                  make_constant_compactification_dofs(displacement.GetElement())
              ),
              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 &GetElement() const {
            return *displacement_space.GetFE(0);
        }

        [[nodiscard]] mfem::Vector MakeZeroElementDofs() const {
            mfem::Vector element_dofs(GetElement().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_determinant, WithinAbs(expected.mapping_determinant, 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 &GetElement() 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::byVDIM
            );
            direction_dofs = make_function_element_dofs(
                element, transformation, evaluate_quadratic_direction, mfem::Ordering::byVDIM
            );
            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 evaluate_reference_hdiv_field(const mfem::Vector &position) {
        const double x = position(0);
        const double y = position(1);
        const double z = position(2);
        return make_vector(x * x + 0.1 * y, y * y - 0.2 * z, z * z + 0.3 * x);
    }

    double evaluate_reference_hdiv_divergence(const mfem::Vector &position) {
        return 2.0 * (position(0) + position(1) + position(2));
    }

    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 Preserves MFEM Ordering",
    tags::unit &tags::mapping
) {
    SingleElementFixture fixture;

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

    mfem::Vector by_vdim_dofs(dof_count * dimension);
    mfem::Vector by_nodes_dofs(dof_count * dimension);

    for (int i = 0; i < dof_count; ++i) {
        for (int component = 0; component < dimension; ++component) {
            const double value                       = 100.0 * component + i + 1.0;
            by_vdim_dofs(component + i * dimension)  = value;
            by_nodes_dofs(i + component * dof_count) = value;
        }
    }

    const mapping::ElementDisplacementData by_vdim_data(element, by_vdim_dofs, mfem::Ordering::byVDIM);
    const mapping::ElementDisplacementData by_nodes_data(element, by_nodes_dofs, mfem::Ordering::byNODES);

    REQUIRE(&by_vdim_data.GetElement() == &element);
    REQUIRE(&by_nodes_data.GetElement() == &element);
    REQUIRE(by_vdim_data.GetDimension() == dimension);
    REQUIRE(by_nodes_data.GetDimension() == dimension);
    REQUIRE(by_vdim_data.GetDofCount() == dof_count);
    REQUIRE(by_nodes_data.GetDofCount() == dof_count);
    REQUIRE(by_vdim_data.GetOrdering() == mfem::Ordering::byVDIM);
    REQUIRE(by_nodes_data.GetOrdering() == mfem::Ordering::byNODES);

    check_matrix(by_vdim_data.GetDofMatrix(), by_nodes_data.GetDofMatrix(), 0.0);

    for (int i = 0; i < dof_count; ++i) {
        for (int component = 0; component < dimension; ++component) {
            CHECK_THAT(by_vdim_data.GetDofMatrix()(i, component), WithinAbs(100.0 * component + i + 1.0, 0.0));
        }
    }
}

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

    const mfem::FiniteElement &element = fixture.GetElement();
    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, mfem::Ordering::byVDIM), std::invalid_argument
    );
    CHECK_THROWS_AS(
        mapping::ElementDisplacementData(element, incomplete_dofs, mfem::Ordering::byVDIM), std::invalid_argument
    );
}

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

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

    workspace.SetDimension(2);
    REQUIRE(workspace.GetDimension() == 2);

    workspace.SetDimension(dimension);
    REQUIRE(workspace.GetDimension() == dimension);

    CHECK_THROWS_AS(workspace.SetDimension(0), std::invalid_argument);
    CHECK_THROWS_AS(workspace.SetDimension(-1), std::invalid_argument);
    CHECK_THROWS_AS(mapping::DomainMapper::Workspace(0), std::invalid_argument);
}

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

    REQUIRE(mapper.GetDimension() == dimension);
    REQUIRE(mapper.GetExteriorMap().GetName() == "KelvinCompactification");

    const utils::DomainMapperOptions invalid_dimension{.dimension = 0, .vacuum_element_attribute = 3};
    const utils::DomainMapperOptions invalid_attribute{.dimension = dimension, .vacuum_element_attribute = 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.GetElement();
    const mfem::Vector zero_dofs       = fixture.MakeZeroElementDofs();
    const mapping::ElementDisplacementData displacement(element, zero_dofs, mfem::Ordering::byVDIM);
    const ElementMappingDataOwner element_data(displacement);

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

    mfem::ElementTransformation *transformation = fixture.mesh.GetElementTransformation(0);
    REQUIRE_FALSE(mapper.IsCompactifiedElement(*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.EvaluatePoint(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_determinant, WithinAbs(1.0, tolerance));
    }
}

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

    const mfem::FiniteElement &element = fixture.GetElement();
    const mfem::Vector zero_dofs       = fixture.MakeZeroElementDofs();
    const mapping::ElementDisplacementData displacement(element, zero_dofs, mfem::Ordering::byVDIM);
    const ElementMappingDataOwner element_data(displacement);

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute = 3}, make_kelvin_compactification()
    );
    mapping::DomainMapper::Workspace 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.EvaluateVolume(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_inv, expected_inverse);
        CHECK_THAT(context.quadrature.detJ, WithinAbs(1.0, tolerance));
        CHECK_THAT(context.quadrature.weight, WithinAbs(expected_weight, tolerance));
        CHECK(context.quadrature.weight > 0.0);
    }
}

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

    const mfem::FiniteElement &element = fixture.GetElement();
    const mfem::Vector zero_dofs       = fixture.MakeZeroElementDofs();
    const mapping::ElementDisplacementData displacement(element, zero_dofs, mfem::Ordering::byVDIM);
    const ElementMappingDataOwner element_data(displacement);

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute = 3}, make_kelvin_compactification()
    );
    mapping::DomainMapper::Workspace 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.EvaluateFace(
                    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.GetElement();
    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);
    const mfem::Vector displacement_offset        = make_vector(0.07, -0.04, 0.03);
    const mfem::Vector element_dofs               = make_affine_element_dofs(
        element, *transformation, displacement_gradient, displacement_offset, mfem::Ordering::byVDIM
    );
    const mapping::ElementDisplacementData displacement(element, element_dofs, mfem::Ordering::byVDIM);
    const ElementMappingDataOwner element_data(displacement);

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute = 3}, make_kelvin_compactification()
    );
    mapping::DomainMapper::Workspace 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);

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

        REQUIRE(
            mapper.EvaluatePoint(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_determinant, WithinAbs(deformation_determinant, tolerance));
    }
}

TEST_CASE(
    "Stateless Domain Mapper Produces Equivalent Results For Both MFEM "
    "Orderings",
    tags::unit &tags::mapping
) {
    SingleElementFixture fixture;

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

    const mfem::DenseMatrix displacement_gradient = make_affine_displacement_gradient();
    const mfem::Vector displacement_offset        = make_vector(0.07, -0.04, 0.03);

    const mfem::Vector by_vdim_dofs               = make_affine_element_dofs(
        element, *transformation, displacement_gradient, displacement_offset, mfem::Ordering::byVDIM
    );
    const mfem::Vector by_nodes_dofs = make_affine_element_dofs(
        element, *transformation, displacement_gradient, displacement_offset, mfem::Ordering::byNODES
    );

    const mapping::ElementDisplacementData by_vdim_displacement(element, by_vdim_dofs, mfem::Ordering::byVDIM);
    const mapping::ElementDisplacementData by_nodes_displacement(element, by_nodes_dofs, mfem::Ordering::byNODES);
    const ElementMappingDataOwner by_vdim_data(by_vdim_displacement);
    const ElementMappingDataOwner by_nodes_data(by_nodes_displacement);

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

    mapping::VolumeMappingContext by_vdim_context;
    mapping::VolumeMappingContext by_nodes_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);

        REQUIRE(
            mapper.EvaluateVolume(by_vdim_data.Get(), *transformation, integration_point, workspace, by_vdim_context) ==
            mapping::MappingStatus::valid
        );
        REQUIRE(
            mapper.EvaluateVolume(
                by_nodes_data.Get(), *transformation, integration_point, workspace, by_nodes_context
            ) == mapping::MappingStatus::valid
        );

        CAPTURE(q);
        check_point_context(by_vdim_context.mapping, by_nodes_context.mapping);
        check_matrix(by_vdim_context.quadrature.J_inv, by_nodes_context.quadrature.J_inv);
        CHECK_THAT(by_vdim_context.quadrature.detJ, WithinAbs(by_nodes_context.quadrature.detJ, tolerance));
        CHECK_THAT(by_vdim_context.quadrature.weight, WithinAbs(by_nodes_context.quadrature.weight, tolerance));
    }
}

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

    const mfem::FiniteElement &element            = fixture.GetElement();
    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);
    const mfem::Vector displacement_offset        = make_vector(0.07, -0.04, 0.03);
    const mfem::Vector element_dofs               = make_affine_element_dofs(
        element, *transformation, displacement_gradient, displacement_offset, mfem::Ordering::byVDIM
    );
    const mapping::ElementDisplacementData displacement(element, element_dofs, mfem::Ordering::byVDIM);
    const ElementMappingDataOwner element_data(displacement);

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

    const double mapping_determinant = mapping_jacobian.Det();
    REQUIRE(mapping_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);
        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_determinant;

        REQUIRE(
            mapper.EvaluateVolume(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_inv, expected_inverse);
        CHECK_THAT(context.quadrature.detJ, WithinAbs(mapping_determinant, tolerance));
        CHECK_THAT(context.quadrature.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.GetElement();
    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);
    const mfem::Vector displacement_offset              = make_vector(0.07, -0.04, 0.03);
    const mfem::Vector element_dofs                     = make_affine_element_dofs(
        element, *element_transformation, displacement_gradient, displacement_offset, mfem::Ordering::byVDIM
    );
    const mapping::ElementDisplacementData displacement(element, element_dofs, mfem::Ordering::byVDIM);
    const ElementMappingDataOwner element_data(displacement);

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

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute = 3}, make_kelvin_compactification()
    );
    mapping::DomainMapper::Workspace 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_determinant;

            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.EvaluateFace(
                    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.GetElement();
    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;

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

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

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

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute = 3}, make_kelvin_compactification()
    );
    mapping::DomainMapper::Workspace 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.EvaluatePoint(element_data_a.Get(), *transformation, integration_point, workspace, first_a) ==
        mapping::MappingStatus::valid
    );
    REQUIRE(
        mapper.EvaluatePoint(element_data_b.Get(), *transformation, integration_point, workspace, result_b) ==
        mapping::MappingStatus::valid
    );
    REQUIRE(
        mapper.EvaluatePoint(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.GetElement();
    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 = 3}, make_kelvin_compactification()
    );
    mapping::DomainMapper::Workspace workspace(dimension);
    mapping::MappingPointContext context;

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

    CHECK(
        mapper.EvaluatePoint(non_finite_data.Get(), *transformation, integration_point, workspace, context) ==
        mapping::MappingStatus::non_finite_input
    );

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

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

    CHECK(
        mapper.EvaluatePoint(singular_data.Get(), *transformation, integration_point, workspace, context) ==
        mapping::MappingStatus::non_positive_determinant
    );

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

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

    CHECK(
        mapper.EvaluatePoint(inverted_data.Get(), *transformation, integration_point, workspace, context) ==
        mapping::MappingStatus::non_positive_determinant
    );

    mapping::DomainMapper::Workspace wrong_workspace(2);
    CHECK_THROWS_AS(
        mapper.EvaluatePoint(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, mfem::Ordering::byVDIM
    );
    const ElementMappingDataOwner two_dimensional_data(two_dimensional_displacement);

    CHECK_THROWS_AS(
        mapper.EvaluatePoint(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.GetElement();
    mfem::ElementTransformation *transformation = fixture.mesh.GetElementTransformation(0);

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

    mapping::DomainMapper mapper(
        {.dimension = dimension, .vacuum_element_attribute = 3}, make_kelvin_compactification()
    );
    mapping::DomainMapper::Workspace 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.EvaluatePoint(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_determinant, 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.GetElement();
    mfem::ElementTransformation *transformation = fixture.mesh.GetElementTransformation(0);
    const QuadraticMappingData mapping_data(element, *transformation);

    const mapping::ElementDisplacementData base_displacement(element, mapping_data.base_dofs, mfem::Ordering::byVDIM);
    const mapping::ElementDisplacementData direction(element, mapping_data.direction_dofs, mfem::Ordering::byVDIM);
    const mapping::ElementDisplacementData plus_displacement(element, mapping_data.plus_dofs, mfem::Ordering::byVDIM);
    const mapping::ElementDisplacementData minus_displacement(element, mapping_data.minus_dofs, mfem::Ordering::byVDIM);
    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 = 3}, make_kelvin_compactification()
    );
    mapping::DomainMapper::Workspace 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.EvaluatePoint(base_data.Get(), *transformation, integration_point, workspace, base_context) ==
            mapping::MappingStatus::valid
        );
        REQUIRE(
            mapper.EvaluatePointVariation(
                base_data.Get(), direction, *transformation, integration_point, base_context, workspace, variation
            ) == mapping::MappingStatus::valid
        );
        REQUIRE(
            mapper.EvaluatePoint(plus_data.Get(), *transformation, integration_point, workspace, plus_context) ==
            mapping::MappingStatus::valid
        );
        REQUIRE(
            mapper.EvaluatePoint(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_determinant - minus_context.mapping_determinant) / (2.0 * difference_step);
        CHECK_THAT(
            determinant_finite_difference, WithinAbs(variation.mapping_determinant_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.GetElement();
    mfem::ElementTransformation *transformation = fixture.mesh.GetElementTransformation(0);
    const QuadraticMappingData mapping_data(element, *transformation);

    const mapping::ElementDisplacementData base_displacement(element, mapping_data.base_dofs, mfem::Ordering::byVDIM);
    const mapping::ElementDisplacementData direction(element, mapping_data.direction_dofs, mfem::Ordering::byVDIM);
    const mapping::ElementDisplacementData plus_displacement(element, mapping_data.plus_dofs, mfem::Ordering::byVDIM);
    const mapping::ElementDisplacementData minus_displacement(element, mapping_data.minus_dofs, mfem::Ordering::byVDIM);
    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 = 3}, make_kelvin_compactification()
    );
    mapping::DomainMapper::Workspace 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.EvaluateVolume(base_data.Get(), *transformation, integration_point, workspace, base_context) ==
            mapping::MappingStatus::valid
        );
        REQUIRE(
            mapper.EvaluateVolumeVariation(
                base_data.Get(), direction, *transformation, integration_point, base_context, workspace, variation
            ) == mapping::MappingStatus::valid
        );
        REQUIRE(
            mapper.EvaluateVolume(plus_data.Get(), *transformation, integration_point, workspace, plus_context) ==
            mapping::MappingStatus::valid
        );
        REQUIRE(
            mapper.EvaluateVolume(minus_data.Get(), *transformation, integration_point, workspace, minus_context) ==
            mapping::MappingStatus::valid
        );

        CAPTURE(q);
        check_matrix_central_difference(
            plus_context.quadrature.J_inv, minus_context.quadrature.J_inv, variation.inverse_element_jacobian_variation,
            difference_step, linearization_tolerance
        );

        const double determinant_finite_difference =
            (plus_context.quadrature.detJ - minus_context.quadrature.detJ) / (2.0 * difference_step);
        const double weight_finite_difference =
            (plus_context.quadrature.weight - minus_context.quadrature.weight) / (2.0 * difference_step);

        CHECK_THAT(
            determinant_finite_difference,
            WithinAbs(variation.mapping.mapping_determinant_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.GetElement();
    mfem::ElementTransformation *element_transformation = fixture.mesh.GetElementTransformation(0);
    const QuadraticMappingData mapping_data(element, *element_transformation);

    const mapping::ElementDisplacementData base_displacement(element, mapping_data.base_dofs, mfem::Ordering::byVDIM);
    const mapping::ElementDisplacementData direction(element, mapping_data.direction_dofs, mfem::Ordering::byVDIM);
    const mapping::ElementDisplacementData plus_displacement(element, mapping_data.plus_dofs, mfem::Ordering::byVDIM);
    const mapping::ElementDisplacementData minus_displacement(element, mapping_data.minus_dofs, mfem::Ordering::byVDIM);
    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 = 3}, make_kelvin_compactification()
    );
    mapping::DomainMapper::Workspace 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.EvaluateFace(
                    base_data.Get(), *transformation, mapping::FaceElementSide::element_1, integration_point, workspace,
                    base_context
                ) == mapping::MappingStatus::valid
            );
            REQUIRE(
                mapper.EvaluateFaceVariation(
                    base_data.Get(), direction, *transformation, mapping::FaceElementSide::element_1, integration_point,
                    base_context, workspace, variation
                ) == mapping::MappingStatus::valid
            );
            REQUIRE(
                mapper.EvaluateFace(
                    plus_data.Get(), *transformation, mapping::FaceElementSide::element_1, integration_point, workspace,
                    plus_context
                ) == mapping::MappingStatus::valid
            );
            REQUIRE(
                mapper.EvaluateFace(
                    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.GetElement();
    mfem::ElementTransformation *transformation     = fixture.mesh.GetElementTransformation(0);
    const mfem::IntegrationPoint &integration_point = mfem::Geometries.GetCenter(transformation->GetGeometryType());
    const QuadraticMappingData mapping_data(element, *transformation);

    const mapping::ElementDisplacementData base_displacement(element, mapping_data.base_dofs, mfem::Ordering::byVDIM);
    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, mfem::Ordering::byVDIM
    );

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

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

    REQUIRE(
        mapper.EvaluatePoint(base_data.Get(), *transformation, integration_point, workspace, base_context) ==
        mapping::MappingStatus::valid
    );
    CHECK(
        mapper.EvaluatePointVariation(
            base_data.Get(), non_finite_direction, *transformation, integration_point, base_context, workspace,
            variation
        ) == mapping::MappingStatus::non_finite_input
    );

    SingleElementFixture linear_fixture;
    const mfem::FiniteElement &linear_element = linear_fixture.GetElement();
    const mfem::Vector linear_direction_dofs  = linear_fixture.MakeZeroElementDofs();
    const mapping::ElementDisplacementData incompatible_direction(
        linear_element, linear_direction_dofs, mfem::Ordering::byVDIM
    );

    CHECK_THROWS_AS(
        mapper.EvaluatePointVariation(
            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_determinant    = context.mapping_jacobian.Det();
    REQUIRE(context.mapping_determinant > 0.0);

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

    const mfem::Vector reference_flux          = make_vector(0.7, -0.4, 1.1);
    const mfem::Vector reference_test_flux     = make_vector(-0.2, 0.9, 0.5);
    const mfem::Vector reference_gradient      = make_vector(0.3, -0.6, 0.8);
    const mfem::Vector reference_test_gradient = make_vector(-0.7, 0.2, 0.4);

    mfem::Vector physical_flux;
    mfem::Vector recovered_flux;
    mfem::Vector physical_test_flux;
    mfem::Vector physical_gradient;
    mfem::Vector recovered_gradient;
    mfem::Vector physical_test_gradient;

    mapping::MapHDivFluxToPhysical(context, reference_flux, physical_flux);
    mapping::MapPhysicalFluxToHDivReference(context, physical_flux, recovered_flux);
    mapping::MapHDivFluxToPhysical(context, reference_test_flux, physical_test_flux);

    mapping::MapReferenceGradientToPhysical(context, reference_gradient, physical_gradient);
    mapping::MapPhysicalGradientToReference(context, physical_gradient, recovered_gradient);
    mapping::MapReferenceGradientToPhysical(context, reference_test_gradient, physical_test_gradient);

    check_vector(recovered_flux, reference_flux, transform_tolerance);
    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::MapReferenceVectorGradientToPhysical(context, reference_vector_gradient, physical_vector_gradient);
    mapping::MapPhysicalVectorGradientToReference(context, physical_vector_gradient, recovered_vector_gradient);
    check_matrix(recovered_vector_gradient, reference_vector_gradient, transform_tolerance);

    mfem::DenseMatrix hdiv_mass_tensor;
    mfem::DenseMatrix diffusion_tensor;
    mapping::ComputeHDivMassTensor(context, hdiv_mass_tensor);
    mapping::ComputeScalarDiffusionTensor(context, diffusion_tensor);

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

    const double physical_hdiv_inner_product  = context.mapping_determinant * (physical_flux * physical_test_flux);
    const double reference_hdiv_inner_product = reference_flux * mass_action;
    const double physical_gradient_inner_product =
        context.mapping_determinant * (physical_gradient * physical_test_gradient);
    const double reference_gradient_inner_product = reference_gradient * diffusion_action;

    CHECK_THAT(physical_hdiv_inner_product, WithinAbs(reference_hdiv_inner_product, transform_tolerance));
    CHECK_THAT(physical_gradient_inner_product, WithinAbs(reference_gradient_inner_product, transform_tolerance));

    const double reference_divergence = 0.73;
    const double physical_divergence  = mapping::MapHDivDivergenceToPhysical(context, reference_divergence);
    CHECK_THAT(physical_divergence, WithinAbs(reference_divergence / context.mapping_determinant, transform_tolerance));
}

TEST_CASE(
    "Mapped Hdiv Flux Preserves Physical Face Flux",
    tags::unit &tags::mapping
) {
    constexpr double flux_tolerance = 2.0e-11;

    SingleElementFixture fixture;

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

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

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

    const mfem::Vector reference_flux = make_vector(0.7, -0.4, 1.1);

    for (int boundary_element = 0; boundary_element < fixture.mesh.GetNBE(); ++boundary_element) {
        mfem::FaceElementTransformations *transformation = fixture.mesh.GetBdrFaceTransformations(boundary_element);
        REQUIRE(transformation != 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);
            mapping::FaceMappingContext context;
            mfem::Vector physical_flux;

            REQUIRE(
                mapper.EvaluateFace(
                    element_data.Get(), *transformation, mapping::FaceElementSide::element_1, integration_point,
                    workspace, context
                ) == mapping::MappingStatus::valid
            );

            mapping::MapHDivFluxToPhysical(context.mapping, reference_flux, physical_flux);

            const double reference_integrated_flux =
                (reference_flux * context.reference_normal) * context.reference_surface_weight;
            const double physical_integrated_flux =
                (physical_flux * context.quadrature.normal) * context.physical_surface_weight;

            CAPTURE(boundary_element, q);
            CHECK_THAT(physical_integrated_flux, WithinAbs(reference_integrated_flux, flux_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>(
            mapping::compactification::options::KelvinCompactificationOptions{
                .r_star_ref = r_star, .r_inf_ref = r_infinity
            }
        );

    mapping::DomainMapper mapper({.dimension = dimension, .vacuum_element_attribute = 3}, std::move(exterior_map));
    mapping::DomainMapper::Workspace 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::ElementDisplacementDataFromElementVDofs(*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.EvaluateVolume(element_data.Get(), *transformation, integration_point, workspace, context) ==
                mapping::MappingStatus::valid
            );

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

            REQUIRE(context.mapping.mapping_determinant > 0.0);
            REQUIRE(context.quadrature.weight > 0.0);

            mfem::Vector reference_flux = make_vector(
                0.4 + context.mapping.reference_position(0), -0.3 + 0.5 * context.mapping.reference_position(1),
                0.7 - 0.2 * context.mapping.reference_position(2)
            );

            mfem::Vector physical_flux;
            mfem::Vector recovered_flux;
            mapping::MapHDivFluxToPhysical(context.mapping, reference_flux, physical_flux);
            mapping::MapPhysicalFluxToHDivReference(context.mapping, physical_flux, recovered_flux);
            check_vector(recovered_flux, reference_flux, transform_tolerance);

            const mfem::Vector reference_gradient = make_vector(0.3, -0.5, 0.8);
            mfem::Vector physical_gradient;
            mfem::Vector recovered_gradient;
            mapping::MapReferenceGradientToPhysical(context.mapping, reference_gradient, physical_gradient);
            mapping::MapPhysicalGradientToReference(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          = coordinate,
                    .compactification_coordinate_gradient = compactification_gradient
                };

                REQUIRE(mapper.GetExteriorMap().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.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);

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

    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, mfem::Ordering::byVDIM);
    const mapping::ElementDisplacementData direction(element, direction_dofs, mfem::Ordering::byVDIM);
    const mapping::ElementDisplacementData plus_displacement(element, plus_dofs, mfem::Ordering::byVDIM);
    const mapping::ElementDisplacementData minus_displacement(element, minus_dofs, mfem::Ordering::byVDIM);

    const 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 = 3}, make_kelvin_compactification()
    );
    mapping::DomainMapper::Workspace 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.EvaluateVolume(base_data.Get(), *transformation, integration_point, workspace, base_context) ==
            mapping::MappingStatus::valid
        );
        REQUIRE(
            mapper.EvaluateVolumeVariation(
                base_data.Get(), direction, *transformation, integration_point, base_context, workspace, variation
            ) == mapping::MappingStatus::valid
        );
        REQUIRE(
            mapper.EvaluateVolume(plus_data.Get(), *transformation, integration_point, workspace, plus_context) ==
            mapping::MappingStatus::valid
        );
        REQUIRE(
            mapper.EvaluateVolume(minus_data.Get(), *transformation, integration_point, workspace, minus_context) ==
            mapping::MappingStatus::valid
        );

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

        const 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
        );
        const mfem::DenseMatrix inverse_element_jacobian_difference =
            matrix_central_difference(plus_context.quadrature.J_inv, minus_context.quadrature.J_inv);

        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)
        );
        CHECK_THAT(
            relative_matrix_difference(
                inverse_element_jacobian_difference, variation.inverse_element_jacobian_variation
            ),
            Catch::Matchers::WithinAbs(0.0, relative_tolerance)
        );

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

        check_centered_difference(
            variation.mapping.mapping_determinant_variation, determinant_difference,
            plus_context.mapping.mapping_determinant, minus_context.mapping.mapping_determinant, kelvin_difference_step,
            relative_tolerance
        );
        check_centered_difference(
            variation.weight_variation, weight_difference, plus_context.quadrature.weight,
            minus_context.quadrature.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.EvaluateFace(
                    base_data.Get(), *face_transformation, mapping::FaceElementSide::element_1, integration_point,
                    workspace, base_context
                ) == mapping::MappingStatus::valid
            );
            REQUIRE(
                mapper.EvaluateFaceVariation(
                    base_data.Get(), direction, *face_transformation, mapping::FaceElementSide::element_1,
                    integration_point, base_context, workspace, variation
                ) == mapping::MappingStatus::valid
            );
            REQUIRE(
                mapper.EvaluateFace(
                    plus_data.Get(), *face_transformation, mapping::FaceElementSide::element_1, integration_point,
                    workspace, plus_context
                ) == mapping::MappingStatus::valid
            );
            REQUIRE(
                mapper.EvaluateFace(
                    minus_data.Get(), *face_transformation, mapping::FaceElementSide::element_1, integration_point,
                    workspace, minus_context
                ) == mapping::MappingStatus::valid
            );

            CAPTURE(boundary_element, q);

            const 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>(
            mapping::compactification::options::KelvinCompactificationOptions{
                .r_star_ref = r_star, .r_inf_ref = r_infinity
            }
        );

    mapping::DomainMapper mapper({.dimension = dimension, .vacuum_element_attribute = 3}, std::move(exterior_map));
    mapping::DomainMapper::Workspace 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::ElementDisplacementDataFromElementVDofs(*displacement_space.GetFE(element_1), displacement_dofs_1);
        const mapping::ElementDisplacementData displacement_2 =
            mapping::ElementDisplacementDataFromElementVDofs(*displacement_space.GetFE(element_2), displacement_dofs_2);
        const mapping::ElementDisplacementData direction_1 =
            mapping::ElementDisplacementDataFromElementVDofs(*displacement_space.GetFE(element_1), direction_dofs_1);
        const mapping::ElementDisplacementData direction_2 =
            mapping::ElementDisplacementDataFromElementVDofs(*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.EvaluateFace(
                element_data_1.Get(), *transformation, mapping::FaceElementSide::element_1, integration_point,
                workspace, context_1
            );
            auto status_2 = mapper.EvaluateFace(
                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.EvaluateFaceVariation(
                    element_data_1.Get(), direction_1, *transformation, mapping::FaceElementSide::element_1,
                    integration_point, context_1, workspace, variation_1
                ) == mapping::MappingStatus::valid
            );
            REQUIRE(
                mapper.EvaluateFaceVariation(
                    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);

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

            mapping::MapPhysicalFluxToHDivReference(context_1.mapping, physical_flux, reference_flux_1);
            mapping::MapPhysicalFluxToHDivReference(context_2.mapping, physical_flux, reference_flux_2);

            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 Hdiv Field Satisfies The Divergence Theorem",
    tags::unit &tags::mapping
) {
    constexpr double divergence_tolerance = 2.0e-9;

    QuadraticElementFixture fixture;

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

    const mfem::Vector displacement_dofs                = make_function_element_dofs(
        element, *element_transformation, evaluate_quadratic_displacement, mfem::Ordering::byVDIM
    );
    const mapping::ElementDisplacementData displacement(element, displacement_dofs, mfem::Ordering::byVDIM);
    const ElementMappingDataOwner element_data(displacement);

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

    double reference_divergence_integral     = 0.0;
    double physical_divergence_integral      = 0.0;
    double physical_boundary_flux            = 0.0;

    const mfem::IntegrationRule &volume_rule = mfem::IntRules.Get(element_transformation->GetGeometryType(), 10);

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

        REQUIRE(
            mapper.EvaluateVolume(element_data.Get(), *element_transformation, integration_point, workspace, context) ==
            mapping::MappingStatus::valid
        );

        element_transformation->SetIntPoint(&integration_point);

        const double reference_divergence = evaluate_reference_hdiv_divergence(context.mapping.reference_position);
        const double physical_divergence  = mapping::MapHDivDivergenceToPhysical(context.mapping, reference_divergence);
        const double reference_weight     = integration_point.weight * element_transformation->Weight();

        reference_divergence_integral += reference_divergence * reference_weight;
        physical_divergence_integral += physical_divergence * context.quadrature.weight;
    }

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

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

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

            REQUIRE(
                mapper.EvaluateFace(
                    element_data.Get(), *transformation, mapping::FaceElementSide::element_1, integration_point,
                    workspace, context
                ) == mapping::MappingStatus::valid
            );

            const mfem::Vector reference_flux = evaluate_reference_hdiv_field(context.mapping.reference_position);
            mfem::Vector physical_flux;
            mapping::MapHDivFluxToPhysical(context.mapping, reference_flux, physical_flux);

            physical_boundary_flux += (physical_flux * context.quadrature.normal) * context.physical_surface_weight;
        }
    }

    constexpr double analytic_reference_integral = 216.0;

    INFO("Analytic reference divergence integral = " << analytic_reference_integral);
    INFO("Computed reference divergence integral = " << reference_divergence_integral);
    INFO("Computed physical divergence integral = " << physical_divergence_integral);
    INFO("Computed physical boundary flux = " << physical_boundary_flux);

    CHECK_THAT(reference_divergence_integral, WithinAbs(analytic_reference_integral, divergence_tolerance));
    CHECK_THAT(physical_divergence_integral, WithinAbs(analytic_reference_integral, divergence_tolerance));
    CHECK_THAT(physical_boundary_flux, WithinAbs(analytic_reference_integral, divergence_tolerance));
}

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_determinant    = context.mapping_jacobian.Det();
    REQUIRE(context.mapping_determinant > 0.0);

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

    const mfem::Vector reference_field      = make_vector(0.7, -0.4, 1.1);
    const 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;

    const 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::MapHCurlFieldToPhysical(context, reference_field, physical_field);
    mapping::MapPhysicalFieldToHCurlReference(context, physical_field, recovered_field);
    mapping::MapHCurlFieldToPhysical(context, reference_test_field, physical_test_field);
    mapping::MapHCurlCurlToPhysical(context, reference_curl, physical_curl);
    mapping::MapPhysicalCurlToHCurlReference(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);

    const 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::ComputeHCurlMassTensor(context, mass_tensor);
    mapping::ComputeHCurlCurlTensor(context, curl_tensor);

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

    const double physical_mass_inner_product  = context.mapping_determinant * (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));

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

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

    const double physical_curl_inner_product  = context.mapping_determinant * (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 = 3}, make_kelvin_compactification()
        );

        mapping::DomainMapper::Workspace 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::ElementDisplacementDataFromElementVDofs(displacement_element, element_displacement);

            REQUIRE(displacement_data.GetOrdering() == mfem::Ordering::byNODES);

            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.GetDofMatrix().MultTranspose(shape, computed_value);

                mfem::MultAtB(displacement_data.GetDofMatrix(), 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.EvaluatePoint(element_data.Get(), *transformation, integration_point, workspace, context);

                REQUIRE(status == mean_field::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_determinant, Catch::Matchers::WithinAbs(expected_determinant, mapping_tolerance)
                );

                REQUIRE(context.mapping_determinant > 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::mapping_evaluator_unit
) {
    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 = 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.EvaluatePoint(*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);

    const 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.EvaluatePoint(*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);

    mfem::Vector physical_position;
    evaluator.GetPhysicalPoint(*transformation, integration_point, physical_position);
    check_vector(physical_position, expected_position);
}

TEST_CASE(
    "Grid Function Mapping Evaluator Invalidates A Cached Compactification Element",
    tags::mapping_evaluator_unit
) {
    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 = 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.EvaluatePoint(*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.InvalidateCache();
    evaluator.InvalidateCache();

    mapping::MappingPointContext refreshed_context;
    REQUIRE(
        evaluator.EvaluatePoint(*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::mapping_evaluator_unit
) {
    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;

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

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

    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.EvaluatePoint(*transformation, integration_point, context) == mapping::MappingStatus::valid);

    displacement.SetSpace(&alternate_displacement_space);

    CHECK_THROWS_AS(evaluator.Refresh(), std::invalid_argument);
    CHECK_THROWS_AS(evaluator.EvaluatePoint(*transformation, integration_point, context), std::invalid_argument);
}