#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>

#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <limits>
#include <memory>

#include <boost/math/quadrature/gauss_kronrod.hpp>
#include <mfem.hpp>
#include <mpi.h>

import mean_field;
import test_helpers;

using namespace mean_field;

namespace {
    double global_vector_norm(
        const mfem::Vector &vector,
        MPI_Comm communicator
    ) {
        const double local_norm_squared = vector * vector;
        double global_norm_squared      = 0.0;
        MPI_Allreduce(&local_norm_squared, &global_norm_squared, 1, MPI_DOUBLE, MPI_SUM, communicator);
        return std::sqrt(global_norm_squared);
    }

    double global_vector_dot(
        const mfem::Vector &lhs,
        const mfem::Vector &rhs,
        MPI_Comm communicator
    ) {
        const double local_dot = lhs * rhs;
        double global_dot      = 0.0;
        MPI_Allreduce(&local_dot, &global_dot, 1, MPI_DOUBLE, MPI_SUM, communicator);
        return global_dot;
    }

    double global_relative_vector_error(
        const mfem::Vector &computed,
        const mfem::Vector &reference,
        MPI_Comm communicator
    ) {
        mfem::Vector difference(computed);
        difference -= reference;
        return global_vector_norm(difference, communicator) /
               std::max(global_vector_norm(reference, communicator), std::numeric_limits<double>::epsilon());
    }

    struct GravitationalEnergies {
        double binding;
        double virial;
    };

    struct HomogeneousEllipsoidAnalytic {
        double coefficient_x;
        double coefficient_y;
        double coefficient_z;
        double energy_kernel;
    };

    class HomogeneousEllipsoidHDivCoefficient : public mfem::VectorCoefficient {
    public:
        HomogeneousEllipsoidHDivCoefficient(
            const mapping::DomainMapper &domain_mapping,
            const double density,
            const HomogeneousEllipsoidAnalytic &analytic
        )
            : VectorCoefficient(3),
              domain_mapping(domain_mapping),
              density(density),
              coefficient_x(analytic.coefficient_x),
              coefficient_y(analytic.coefficient_y),
              coefficient_z(analytic.coefficient_z) {
        }

        void Eval(
            mfem::Vector &value,
            mfem::ElementTransformation &transformation,
            const mfem::IntegrationPoint &integration_point
        ) override {
            transformation.SetIntPoint(&integration_point);

            mfem::Vector x_physical(3);
            mfem::Vector field_physical(3);
            mfem::DenseMatrix map_jacobian(3, 3);
            mfem::DenseMatrix inverse_map_jacobian(3, 3);

            domain_mapping.GetPhysicalPoint(transformation, integration_point, x_physical);

            field_physical(0) = 2.0 * M_PI * utils::G * density * coefficient_x * x_physical(0);
            field_physical(1) = 2.0 * M_PI * utils::G * density * coefficient_y * x_physical(1);
            field_physical(2) = 2.0 * M_PI * utils::G * density * coefficient_z * x_physical(2);

            domain_mapping.ComputeJacobian(transformation, map_jacobian);
            const double map_determinant = map_jacobian.Det();
            MFEM_VERIFY(map_determinant > 0.0, "Domain mapping has a non-positive Jacobian determinant.");

            mfem::CalcInverse(map_jacobian, inverse_map_jacobian);
            inverse_map_jacobian.Mult(field_physical, value);
            value *= map_determinant;
        }

    private:
        const mapping::DomainMapper &domain_mapping;
        double density;
        double coefficient_x;
        double coefficient_y;
        double coefficient_z;
    };

    template <typename GravitySolutionType>
    GravitationalEnergies compute_gravitational_energies(
        fem::FEM &f,
        const mfem::GridFunction &rho,
        const GravitySolutionType &gravity_solution,
        const int quadrature_order
    ) {
        const int dim                = f.mesh->Dimension();
        double local_bind_integral   = 0.0;
        double local_virial_integral = 0.0;

        mfem::Vector x_physical(dim);
        mfem::Vector grad_phi_element(dim);
        mfem::Vector grad_phi_physical(dim);
        mfem::DenseMatrix map_jacobian(dim, dim);

        for (int elem_id = 0; elem_id < f.mesh->GetNE(); ++elem_id) {
            if (f.mesh->GetAttribute(elem_id) == 3) {
                continue;
            }

            mfem::ElementTransformation *transformation = f.mesh->GetElementTransformation(elem_id);
            const mfem::IntegrationRule &integration_rule =
                mfem::IntRules.Get(transformation->GetGeometryType(), quadrature_order);

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

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

                if (f.has_mapping()) {
                    const double map_determinant = f.mapping->ComputeDetJ(*transformation, integration_point);
                    MFEM_VERIFY(
                        map_determinant > 0.0, "Domain mapping has a non-positive Jacobian "
                                               "determinant."
                    );

                    weight *= map_determinant;
                    f.mapping->GetPhysicalPoint(*transformation, integration_point, x_physical);
                    gravity_solution.gradPhi.GetVectorValue(elem_id, integration_point, grad_phi_element);
                    f.mapping->ComputeJacobian(*transformation, map_jacobian);
                    map_jacobian.Mult(grad_phi_element, grad_phi_physical);
                    grad_phi_physical /= map_determinant;
                } else {
                    transformation->Transform(integration_point, x_physical);
                    gravity_solution.gradPhi.GetVectorValue(elem_id, integration_point, grad_phi_physical);
                }

                const double rho_value     = rho.GetValue(elem_id, integration_point);
                const double phi_value     = gravity_solution.phi.GetValue(elem_id, integration_point);
                double radius_dot_gradient = 0.0;

                for (int d = 0; d < dim; ++d) {
                    radius_dot_gradient += (x_physical(d) - f.com(d)) * grad_phi_physical(d);
                }

                local_bind_integral += rho_value * phi_value * weight;
                local_virial_integral += rho_value * radius_dot_gradient * weight;
            }
        }

        const double local_w_bind = 0.5 * local_bind_integral;
        const double local_w_vir  = -local_virial_integral;
        double global_w_bind      = 0.0;
        double global_w_vir       = 0.0;
        MPI_Comm communicator     = f.densityFes->GetComm();

        MPI_Allreduce(&local_w_bind, &global_w_bind, 1, MPI_DOUBLE, MPI_SUM, communicator);
        MPI_Allreduce(&local_w_vir, &global_w_vir, 1, MPI_DOUBLE, MPI_SUM, communicator);

        return {.binding = global_w_bind, .virial = global_w_vir};
    }

    void zero_vacuum_density(
        const fem::FEM &f,
        mfem::GridFunction &rho
    ) {
        for (int i = 0; i < f.vacuumDensityTdofs.Size(); ++i) {
            rho(f.vacuumDensityTdofs[i]) = 0.0;
        }
    }

    int get_gravity_quadrature_order(const fem::FEM &f) {
        return 2 * std::max(f.gravityPotentialFes->GetMaxElementOrder(), f.gravityFluxFes->GetMaxElementOrder()) + 8;
    }

    double compute_ellipsoid_coefficient(
        const double normalized_axis_x,
        const double normalized_axis_y,
        const double normalized_axis_z,
        const double target_axis_squared
    ) {
        auto integrand = [=](const double t) {
            if (t <= 0.0 || t >= 1.0) {
                return 0.0;
            }

            const double one_minus_t   = 1.0 - t;
            const double s             = t / one_minus_t;
            const double s_squared     = s * s;
            const double ds_squared_dt = 2.0 * s / (one_minus_t * one_minus_t);
            const double delta         = std::sqrt(
                (normalized_axis_x * normalized_axis_x + s_squared) *
                (normalized_axis_y * normalized_axis_y + s_squared) *
                (normalized_axis_z * normalized_axis_z + s_squared)
            );

            return normalized_axis_x * normalized_axis_y * normalized_axis_z * ds_squared_dt /
                   ((target_axis_squared + s_squared) * delta);
        };

        double integration_error = 0.0;
        return boost::math::quadrature::gauss_kronrod<double, 61>::integrate(
            integrand, 0.0, 1.0, 15, 1.0e-13, &integration_error
        );
    }

    double compute_ellipsoid_energy_kernel(
        const double normalized_axis_x,
        const double normalized_axis_y,
        const double normalized_axis_z,
        const double length_scale
    ) {
        auto integrand = [=](const double t) {
            if (t <= 0.0) {
                return 0.0;
            }
            if (t >= 1.0) {
                return 2.0;
            }

            const double one_minus_t   = 1.0 - t;
            const double s             = t / one_minus_t;
            const double s_squared     = s * s;
            const double ds_squared_dt = 2.0 * s / (one_minus_t * one_minus_t);
            const double delta         = std::sqrt(
                (normalized_axis_x * normalized_axis_x + s_squared) *
                (normalized_axis_y * normalized_axis_y + s_squared) *
                (normalized_axis_z * normalized_axis_z + s_squared)
            );

            return ds_squared_dt / delta;
        };

        double integration_error            = 0.0;
        const double dimensionless_integral = boost::math::quadrature::gauss_kronrod<double, 61>::integrate(
            integrand, 0.0, 1.0, 15, 1.0e-13, &integration_error
        );

        return dimensionless_integral / length_scale;
    }

    HomogeneousEllipsoidAnalytic compute_homogeneous_ellipsoid_analytic(
        const double semi_axis_x,
        const double semi_axis_y,
        const double semi_axis_z
    ) {
        const double length_scale      = std::cbrt(semi_axis_x * semi_axis_y * semi_axis_z);
        const double normalized_axis_x = semi_axis_x / length_scale;
        const double normalized_axis_y = semi_axis_y / length_scale;
        const double normalized_axis_z = semi_axis_z / length_scale;
        const double coefficient_x     = compute_ellipsoid_coefficient(
            normalized_axis_x, normalized_axis_y, normalized_axis_z, normalized_axis_x * normalized_axis_x
        );
        const double coefficient_y = compute_ellipsoid_coefficient(
            normalized_axis_x, normalized_axis_y, normalized_axis_z, normalized_axis_y * normalized_axis_y
        );
        const double coefficient_z = compute_ellipsoid_coefficient(
            normalized_axis_x, normalized_axis_y, normalized_axis_z, normalized_axis_z * normalized_axis_z
        );
        const double energy_kernel =
            compute_ellipsoid_energy_kernel(normalized_axis_x, normalized_axis_y, normalized_axis_z, length_scale);

        return {
            .coefficient_x = coefficient_x,
            .coefficient_y = coefficient_y,
            .coefficient_z = coefficient_z,
            .energy_kernel = energy_kernel
        };
    }

    struct GravitySolutionComparison {
        double relative_gradient_difference;
        double relative_potential_difference;
    };

    template <
        typename LeftGravitySolution,
        typename RightGravitySolution>
    GravitySolutionComparison compare_gravity_solutions(
        fem::FEM &f,
        const LeftGravitySolution &left_solution,
        const RightGravitySolution &right_solution,
        const int quadrature_order
    ) {
        const int dimension                       = f.mesh->Dimension();

        double local_gradient_difference_squared  = 0.0;
        double local_left_gradient_norm_squared   = 0.0;
        double local_right_gradient_norm_squared  = 0.0;
        double local_potential_difference_squared = 0.0;
        double local_left_potential_norm_squared  = 0.0;
        double local_right_potential_norm_squared = 0.0;

        mfem::Vector left_gradient_element(dimension);
        mfem::Vector right_gradient_element(dimension);
        mfem::Vector left_gradient_physical(dimension);
        mfem::Vector right_gradient_physical(dimension);
        mfem::Vector gradient_difference(dimension);
        mfem::DenseMatrix map_jacobian(dimension, dimension);

        for (int element_id = 0; element_id < f.mesh->GetNE(); ++element_id) {
            if (f.mesh->GetAttribute(element_id) == 3) {
                continue;
            }

            mfem::ElementTransformation *transformation = f.mesh->GetElementTransformation(element_id);
            const mfem::IntegrationRule &integration_rule =
                mfem::IntRules.Get(transformation->GetGeometryType(), quadrature_order);

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

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

                left_solution.gradPhi.GetVectorValue(element_id, integration_point, left_gradient_element);
                right_solution.gradPhi.GetVectorValue(element_id, integration_point, right_gradient_element);

                if (f.has_mapping()) {
                    const double map_determinant = f.mapping->ComputeDetJ(*transformation, integration_point);
                    MFEM_VERIFY(
                        map_determinant > 0.0, "Domain mapping has a non-positive Jacobian "
                                               "determinant."
                    );

                    weight *= map_determinant;
                    f.mapping->ComputeJacobian(*transformation, map_jacobian);

                    map_jacobian.Mult(left_gradient_element, left_gradient_physical);
                    map_jacobian.Mult(right_gradient_element, right_gradient_physical);

                    left_gradient_physical /= map_determinant;
                    right_gradient_physical /= map_determinant;
                } else {
                    left_gradient_physical  = left_gradient_element;
                    right_gradient_physical = right_gradient_element;
                }

                gradient_difference = left_gradient_physical;
                gradient_difference -= right_gradient_physical;

                local_gradient_difference_squared += (gradient_difference * gradient_difference) * weight;
                local_left_gradient_norm_squared += (left_gradient_physical * left_gradient_physical) * weight;
                local_right_gradient_norm_squared += (right_gradient_physical * right_gradient_physical) * weight;

                const double left_potential       = left_solution.phi.GetValue(element_id, integration_point);
                const double right_potential      = right_solution.phi.GetValue(element_id, integration_point);
                const double potential_difference = left_potential - right_potential;

                local_potential_difference_squared += potential_difference * potential_difference * weight;
                local_left_potential_norm_squared += left_potential * left_potential * weight;
                local_right_potential_norm_squared += right_potential * right_potential * weight;
            }
        }

        const std::array<double, 6> local_values{local_gradient_difference_squared, local_left_gradient_norm_squared,
                                                 local_right_gradient_norm_squared, local_potential_difference_squared,
                                                 local_left_potential_norm_squared, local_right_potential_norm_squared};

        std::array<double, 6> global_values{};

        MPI_Allreduce(
            local_values.data(), global_values.data(), static_cast<int>(local_values.size()), MPI_DOUBLE, MPI_SUM,
            f.densityFes->GetComm()
        );

        const double gradient_scale_squared  = 0.5 * (global_values[1] + global_values[2]);
        const double potential_scale_squared = 0.5 * (global_values[4] + global_values[5]);

        MFEM_VERIFY(gradient_scale_squared > 0.0, "Cannot compare gravity solutions with zero gradient norm.");
        MFEM_VERIFY(potential_scale_squared > 0.0, "Cannot compare gravity solutions with zero potential norm.");

        return {
            .relative_gradient_difference  = std::sqrt(global_values[0] / gradient_scale_squared),
            .relative_potential_difference = std::sqrt(global_values[3] / potential_scale_squared)
        };
    }

    struct GravityResidualMetrics {
        double relative_total;
        double relative_gradient;
        double relative_poisson;
    };

    double gravity_test_global_norm(
        const mfem::Vector &vector,
        MPI_Comm communicator
    ) {
        const double local_norm_squared = vector * vector;
        double global_norm_squared      = 0.0;
        MPI_Allreduce(&local_norm_squared, &global_norm_squared, 1, MPI_DOUBLE, MPI_SUM, communicator);
        return std::sqrt(global_norm_squared);
    }

    void gravity_test_get_true_dofs(
        const mfem::ParFiniteElementSpace &finite_element_space,
        const mfem::GridFunction &grid_function,
        mfem::Vector &true_dofs
    ) {
        MFEM_VERIFY(grid_function.Size() == finite_element_space.GetVSize(), "Grid function has the wrong local size.");

        true_dofs.SetSize(finite_element_space.GetTrueVSize());

        const mfem::Operator *restriction = finite_element_space.GetRestrictionMatrix();

        if (restriction != nullptr) {
            restriction->Mult(grid_function, true_dofs);
        } else {
            MFEM_VERIFY(grid_function.Size() == true_dofs.Size(), "Local and true sizes do not match.");
            true_dofs = grid_function;
        }
    }

    GravityResidualMetrics measure_gravity_residual(
        const mfem::Vector &residual,
        const mfem::Array<int> &offsets,
        const double right_hand_side_norm,
        MPI_Comm communicator
    ) {
        MFEM_VERIFY(offsets.Size() == 3, "Gravity residual must contain two blocks.");
        MFEM_VERIFY(residual.Size() == offsets.Last(), "Gravity residual has the wrong size.");
        MFEM_VERIFY(right_hand_side_norm > 0.0, "Gravity right-hand side must be nonzero.");

        mfem::Vector gradient_residual(offsets[1] - offsets[0]);
        mfem::Vector poisson_residual(offsets[2] - offsets[1]);

        for (int i = 0; i < gradient_residual.Size(); ++i) {
            gradient_residual(i) = residual(offsets[0] + i);
        }

        for (int i = 0; i < poisson_residual.Size(); ++i) {
            poisson_residual(i) = residual(offsets[1] + i);
        }

        return {
            .relative_total    = gravity_test_global_norm(residual, communicator) / right_hand_side_norm,
            .relative_gradient = gravity_test_global_norm(gradient_residual, communicator) / right_hand_side_norm,
            .relative_poisson  = gravity_test_global_norm(poisson_residual, communicator) / right_hand_side_norm
        };
    }

    enum class ExteriorMonopoleMapping { legacy, stateless };

    struct ExteriorMonopoleShellMetrics {
        long long quadrature_points{0};
        double minimum_radius{std::numeric_limits<double>::infinity()};
        double maximum_radius{0.0};
        double potential_rms_error{0.0};
        double radial_field_rms_error{0.0};
        double tangential_field_rms{0.0};
    };

    struct ExteriorMonopoleShellAccumulator {
        long long quadrature_points{0};
        double minimum_radius{std::numeric_limits<double>::infinity()};
        double maximum_radius{0.0};
        double reference_weight{0.0};
        double potential_error_squared{0.0};
        double radial_field_error_squared{0.0};
        double tangential_field_squared{0.0};
    };

    constexpr std::array<double, 6> exterior_shell_boundaries{0.0, 0.25, 0.50, 0.75, 0.90, 1.0};

    int get_exterior_shell(const double compactification_coordinate) {
        REQUIRE(std::isfinite(compactification_coordinate));
        REQUIRE(compactification_coordinate >= -1.0e-12);
        REQUIRE(compactification_coordinate <= 1.0 + 1.0e-12);

        const double coordinate = std::clamp(compactification_coordinate, 0.0, std::nextafter(1.0, 0.0));

        for (int shell = 0; shell < static_cast<int>(exterior_shell_boundaries.size()) - 1; ++shell) {
            if (coordinate < exterior_shell_boundaries[shell + 1]) {
                return shell;
            }
        }

        return static_cast<int>(exterior_shell_boundaries.size()) - 2;
    }

    std::array<
        ExteriorMonopoleShellMetrics,
        5>
    measure_exterior_monopole_shells(
        fem::FEM &f,
        const physics::GravitySolution &solution,
        const mfem::GridFunction &displacement,
        const ExteriorMonopoleMapping mapping_path,
        const double mass
    ) {
        REQUIRE(f.mesh != nullptr);
        REQUIRE(f.gravityFluxFes != nullptr);
        REQUIRE(f.displacementFes != nullptr);
        REQUIRE(f.compactificationFes != nullptr);
        REQUIRE(f.compactificationCoordinate != nullptr);
        REQUIRE(f.mapping != nullptr);
        REQUIRE(f.domainMapperStateless != nullptr);

        constexpr int shell_count = static_cast<int>(exterior_shell_boundaries.size()) - 1;

        std::array<ExteriorMonopoleShellAccumulator, shell_count> local_shells{};

        mapping::DomainMapperStateless::Workspace workspace(f.mesh->Dimension());

        const int vacuum_attribute = f.domainMapperStateless->GetVacuumElementAttribute();

        const int quadrature_order = get_gravity_quadrature_order(f);

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

            REQUIRE(transformation != nullptr);

            if (transformation->Attribute != vacuum_attribute) {
                continue;
            }

            const mfem::FiniteElement &displacement_element     = *f.displacementFes->GetFE(element_id);

            const mfem::FiniteElement &compactification_element = *f.compactificationFes->GetFE(element_id);

            mfem::Array<int> displacement_dofs;
            mfem::Array<int> compactification_dofs;

            mfem::DofTransformation *displacement_dof_transformation =
                f.displacementFes->GetElementVDofs(element_id, displacement_dofs);

            mfem::DofTransformation *compactification_dof_transformation =
                f.compactificationFes->GetElementDofs(element_id, compactification_dofs);

            mfem::Vector element_displacement;
            mfem::Vector element_compactification;

            displacement.GetSubVector(displacement_dofs, element_displacement);

            f.compactificationCoordinate->GetSubVector(compactification_dofs, element_compactification);

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

            if (compactification_dof_transformation != nullptr) {
                compactification_dof_transformation->InvTransformPrimal(element_compactification);
            }

            const mapping::ElementDisplacementData displacement_data(
                displacement_element, element_displacement, f.displacementFes->GetOrdering()
            );

            const mapping::ElementCompactificationData compactification_data(
                compactification_element, element_compactification
            );

            const mapping::ElementMappingData mapping_data{
                .displacement = displacement_data, .compactification = compactification_data
            };

            mfem::Vector compactification_shape(compactification_element.GetDof());

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

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

                transformation->SetIntPoint(&integration_point);

                compactification_element.CalcShape(integration_point, compactification_shape);

                const double compactification_coordinate = element_compactification * compactification_shape;

                const int shell                          = get_exterior_shell(compactification_coordinate);

                mfem::Vector reference_field(3);
                mfem::Vector physical_field(3);
                mfem::Vector physical_position(3);

                solution.gradPhi.GetVectorValue(element_id, integration_point, reference_field);

                if (mapping_path == ExteriorMonopoleMapping::stateless) {
                    mapping::VolumeMappingContext mapping_context;

                    const mapping::MappingStatus status = f.domainMapperStateless->EvaluateVolume(
                        mapping_data, *transformation, integration_point, workspace, mapping_context
                    );

                    CAPTURE(element_id, q, compactification_coordinate, static_cast<int>(status));

                    REQUIRE(status == mean_field::mapping::MappingStatus::valid);

                    physical_position = mapping_context.mapping.physical_position;

                    mean_field::mapping::MapHDivFluxToPhysical(
                        mapping_context.mapping, reference_field, physical_field
                    );
                } else {
                    /*
                     * The legacy path intentionally does not use the exterior
                     * coordinate to construct its physical mapping.
                     */
                    f.mapping->GetPhysicalPoint(*transformation, integration_point, physical_position);

                    mfem::DenseMatrix mapping_jacobian(3);
                    f.mapping->ComputeJacobian(*transformation, mapping_jacobian);

                    const double mapping_determinant = mapping_jacobian.Det();

                    REQUIRE(std::isfinite(mapping_determinant));
                    REQUIRE(mapping_determinant > 0.0);

                    mapping_jacobian.Mult(reference_field, physical_field);

                    physical_field /= mapping_determinant;
                }

                const double radius = physical_position.Norml2();

                CAPTURE(element_id, q, shell, compactification_coordinate, radius);

                REQUIRE(std::isfinite(radius));
                REQUIRE(radius > 0.0);

                mfem::Vector radial_unit_vector(physical_position);
                radial_unit_vector /= radius;

                const double numerical_radial_field = physical_field * radial_unit_vector;

                mfem::Vector tangential_field(physical_field);
                tangential_field.Add(-numerical_radial_field, radial_unit_vector);

                const double numerical_potential     = solution.phi.GetValue(element_id, integration_point);

                /*
                 * For an exterior monopole:
                 *
                 *     phi       = -GM/r
                 *     grad(phi) =  GM r_hat/r^2
                 *
                 * These scaled quantities should therefore be one, one, and
                 * zero respectively. They remain well-conditioned as r -> inf.
                 */
                const double scaled_potential        = -radius * numerical_potential / (utils::G * mass);

                const double scaled_radial_field     = radius * radius * numerical_radial_field / (utils::G * mass);

                const double scaled_tangential_field = radius * radius * tangential_field.Norml2() / (utils::G * mass);

                REQUIRE(std::isfinite(scaled_potential));
                REQUIRE(std::isfinite(scaled_radial_field));
                REQUIRE(std::isfinite(scaled_tangential_field));

                /*
                 * Use the finite reference-domain measure for averaging. A
                 * physical L2 norm of phi over an infinite three-dimensional
                 * exterior domain is not finite.
                 */
                const double reference_weight                 = integration_point.weight * transformation->Weight();

                ExteriorMonopoleShellAccumulator &accumulator = local_shells[shell];

                ++accumulator.quadrature_points;

                accumulator.minimum_radius = std::min(accumulator.minimum_radius, radius);

                accumulator.maximum_radius = std::max(accumulator.maximum_radius, radius);

                accumulator.reference_weight += reference_weight;

                accumulator.potential_error_squared += reference_weight * std::pow(scaled_potential - 1.0, 2);

                accumulator.radial_field_error_squared += reference_weight * std::pow(scaled_radial_field - 1.0, 2);

                accumulator.tangential_field_squared +=
                    reference_weight * scaled_tangential_field * scaled_tangential_field;
            }
        }

        MPI_Comm communicator = f.gravityFluxFes->GetComm();

        std::array<ExteriorMonopoleShellMetrics, shell_count> metrics{};

        for (int shell = 0; shell < shell_count; ++shell) {
            long long global_points = 0;

            MPI_Allreduce(
                &local_shells[shell].quadrature_points, &global_points, 1, MPI_LONG_LONG, MPI_SUM, communicator
            );

            double local_sums[4]{
                local_shells[shell].reference_weight, local_shells[shell].potential_error_squared,
                local_shells[shell].radial_field_error_squared, local_shells[shell].tangential_field_squared
            };

            double global_sums[4]{};

            MPI_Allreduce(local_sums, global_sums, 4, MPI_DOUBLE, MPI_SUM, communicator);

            double global_minimum_radius = 0.0;
            double global_maximum_radius = 0.0;

            MPI_Allreduce(
                &local_shells[shell].minimum_radius, &global_minimum_radius, 1, MPI_DOUBLE, MPI_MIN, communicator
            );

            MPI_Allreduce(
                &local_shells[shell].maximum_radius, &global_maximum_radius, 1, MPI_DOUBLE, MPI_MAX, communicator
            );

            REQUIRE(global_points > 0);
            REQUIRE(global_sums[0] > 0.0);

            metrics[shell] = {
                .quadrature_points      = global_points,
                .minimum_radius         = global_minimum_radius,
                .maximum_radius         = global_maximum_radius,
                .potential_rms_error    = std::sqrt(global_sums[1] / global_sums[0]),
                .radial_field_rms_error = std::sqrt(global_sums[2] / global_sums[0]),
                .tangential_field_rms   = std::sqrt(global_sums[3] / global_sums[0])
            };
        }

        return metrics;
    }

    class StatelessProjectionGeometry {
    public:
        StatelessProjectionGeometry(
            const fem::FEM &f,
            const mapping::DomainMapperStateless &domain_mapper,
            const mfem::GridFunction &displacement
        )
            : m_fem(f),
              m_domain_mapper(domain_mapper),
              m_displacement(displacement),
              m_workspace(f.mesh->Dimension()) {
        }

        mapping::MappingStatus Evaluate(
            mfem::ElementTransformation &transformation,
            const mfem::IntegrationPoint &integration_point,
            mapping::MappingPointContext &context,
            const bool permit_infinity_limit
        ) {
            m_last_evaluation_used_infinity_limit = false;

            const int element_id                  = transformation.ElementNo;

            MFEM_VERIFY(element_id >= 0, "Projection coefficient received an invalid element number.");

            const mfem::FiniteElement &displacement_element     = *m_fem.displacementFes->GetFE(element_id);

            const mfem::FiniteElement &compactification_element = *m_fem.compactificationFes->GetFE(element_id);

            mfem::Array<int> displacement_dofs;
            mfem::Array<int> compactification_dofs;

            mfem::DofTransformation *displacement_dof_transformation =
                m_fem.displacementFes->GetElementVDofs(element_id, displacement_dofs);

            mfem::DofTransformation *compactification_dof_transformation =
                m_fem.compactificationFes->GetElementDofs(element_id, compactification_dofs);

            mfem::Vector element_displacement;
            mfem::Vector element_compactification;

            m_displacement.GetSubVector(displacement_dofs, element_displacement);

            m_fem.compactificationCoordinate->GetSubVector(compactification_dofs, element_compactification);

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

            if (compactification_dof_transformation != nullptr) {
                compactification_dof_transformation->InvTransformPrimal(element_compactification);
            }

            const mapping::ElementDisplacementData displacement_data(
                displacement_element, element_displacement, m_fem.displacementFes->GetOrdering()
            );

            const mapping::ElementCompactificationData compactification_data(
                compactification_element, element_compactification
            );

            mfem::Vector requested_compactification_shape(compactification_element.GetDof());

            compactification_element.CalcShape(integration_point, requested_compactification_shape);

            const double requested_compactification_coordinate =
                element_compactification * requested_compactification_shape;

            constexpr double infinity_candidate_tolerance = 1.0e-8;

            const bool requested_infinity_limit =
                std::isfinite(requested_compactification_coordinate) &&
                requested_compactification_coordinate >= 1.0 - infinity_candidate_tolerance;

            const mapping::ElementMappingData mapping_data{
                .displacement = displacement_data, .compactification = compactification_data
            };

            mapping::MappingStatus status =
                m_domain_mapper.EvaluatePoint(mapping_data, transformation, integration_point, m_workspace, context);

            if (status == mapping::MappingStatus::valid) {
                transformation.SetIntPoint(&integration_point);
                return status;
            }

            if (!permit_infinity_limit || !m_domain_mapper.IsCompactifiedElement(transformation)) {
                transformation.SetIntPoint(&integration_point);
                return status;
            }

            const bool retryable_boundary_status =
                status == mapping::MappingStatus::at_compactified_infinity ||
                status == mapping::MappingStatus::outside_reference_domain ||
                status == mapping::MappingStatus::non_finite_result ||
                (requested_infinity_limit && status == mapping::MappingStatus::non_positive_determinant);

            if (!retryable_boundary_status) {
                transformation.SetIntPoint(&integration_point);
                return status;
            }

            const mfem::IntegrationPoint &element_center = mfem::Geometries.GetCenter(transformation.GetGeometryType());

            /*
             * Use the nearest admissible point. Starting extremely close to
             * the requested point preserves the limiting RT trace, while the
             * larger fallbacks accommodate the mapper's infinity guard.
             */
            constexpr std::array<double, 9> inward_fractions{1.0e-12, 1.0e-11, 1.0e-10, 1.0e-9, 1.0e-8,
                                                             1.0e-7,  1.0e-6,  1.0e-5,  1.0e-4};

            for (const double inward_fraction : inward_fractions) {
                mfem::IntegrationPoint inward_point;

                inward_point.x = (1.0 - inward_fraction) * integration_point.x + inward_fraction * element_center.x;

                inward_point.y = (1.0 - inward_fraction) * integration_point.y + inward_fraction * element_center.y;

                inward_point.z = (1.0 - inward_fraction) * integration_point.z + inward_fraction * element_center.z;

                inward_point.weight = integration_point.weight;

                status =
                    m_domain_mapper.EvaluatePoint(mapping_data, transformation, inward_point, m_workspace, context);

                if (status == mapping::MappingStatus::valid) {
                    m_last_evaluation_used_infinity_limit = true;
                    transformation.SetIntPoint(&integration_point);
                    return status;
                }

                const bool still_retryable =
                    status == mapping::MappingStatus::at_compactified_infinity ||
                    status == mapping::MappingStatus::outside_reference_domain ||
                    status == mapping::MappingStatus::non_finite_result ||
                    (requested_infinity_limit && status == mapping::MappingStatus::non_positive_determinant);
                if (!still_retryable) {
                    break;
                }
            }

            transformation.SetIntPoint(&integration_point);
            return status;
        }

        [[nodiscard]]
        bool LastEvaluationUsedInfinityLimit() const noexcept {
            return m_last_evaluation_used_infinity_limit;
        }

    private:
        const fem::FEM &m_fem;
        const mapping::DomainMapperStateless &m_domain_mapper;
        const mfem::GridFunction &m_displacement;
        mapping::DomainMapperStateless::Workspace m_workspace;
        bool m_last_evaluation_used_infinity_limit{false};
    };
    class StatelessMonopolePotentialCoefficient final : public mfem::Coefficient {
    public:
        StatelessMonopolePotentialCoefficient(
            const fem::FEM &f,
            const mapping::DomainMapperStateless &domain_mapper,
            const mfem::GridFunction &displacement,
            const double mass,
            const double stellar_radius
        )
            : m_geometry(
                  f,
                  domain_mapper,
                  displacement
              ),
              m_vacuum_attribute(domain_mapper.GetVacuumElementAttribute()),
              m_mass(mass),
              m_stellar_radius(stellar_radius) {
        }

        double Eval(
            mfem::ElementTransformation &transformation,
            const mfem::IntegrationPoint &integration_point
        ) override {
            mapping::MappingPointContext context;

            const mapping::MappingStatus status = m_geometry.Evaluate(transformation, integration_point, context, true);

            MFEM_VERIFY(
                status == mean_field::mapping::MappingStatus::valid,
                "Stateless monopole-potential projection failed."
                    << "\nMapping status = " << static_cast<int>(status)
                    << "\nElement ID = " << transformation.ElementNo
                    << "\nElement attribute = " << transformation.Attribute << "\nIntegration point = <"
                    << integration_point.x << ", " << integration_point.y << ", " << integration_point.z << ">"
            );

            if (m_geometry.LastEvaluationUsedInfinityLimit()) {
                return 0.0;
            }

            const double radius = context.physical_position.Norml2();

            MFEM_VERIFY(
                std::isfinite(radius) && radius > 0.0, "Monopole projection encountered an invalid physical radius."
            );

            if (transformation.Attribute == m_vacuum_attribute) {
                return -utils::G * m_mass / radius;
            }

            return -utils::G * m_mass / (2.0 * m_stellar_radius * m_stellar_radius * m_stellar_radius) *
                   (3.0 * m_stellar_radius * m_stellar_radius - radius * radius);
        }

    private:
        StatelessProjectionGeometry m_geometry;
        int m_vacuum_attribute;
        double m_mass;
        double m_stellar_radius;
    };

    class StatelessMonopoleHDivCoefficient final : public mfem::VectorCoefficient {
    public:
        StatelessMonopoleHDivCoefficient(
            const fem::FEM &f,
            const mapping::DomainMapperStateless &domain_mapper,
            const mfem::GridFunction &displacement,
            const double mass,
            const double stellar_radius
        )
            : VectorCoefficient(f.mesh->Dimension()),
              m_geometry(
                  f,
                  domain_mapper,
                  displacement
              ),
              m_vacuum_attribute(domain_mapper.GetVacuumElementAttribute()),
              m_mass(mass),
              m_stellar_radius(stellar_radius) {
        }

        void Eval(
            mfem::Vector &value,
            mfem::ElementTransformation &transformation,
            const mfem::IntegrationPoint &integration_point
        ) override {
            mapping::MappingPointContext context;

            const mapping::MappingStatus status = m_geometry.Evaluate(transformation, integration_point, context, true);

            MFEM_VERIFY(
                status == mean_field::mapping::MappingStatus::valid,
                "Stateless monopole H(div) projection failed."
                    << "\nMapping status = " << static_cast<int>(status)
                    << "\nElement ID = " << transformation.ElementNo
                    << "\nElement attribute = " << transformation.Attribute << "\nIntegration point = <"
                    << integration_point.x << ", " << integration_point.y << ", " << integration_point.z << ">"
                    << "\nInfinity-limit evaluation attempted = " << m_geometry.LastEvaluationUsedInfinityLimit()
            );
            const mfem::Vector &displaced_position = context.displaced_position;

            const double computational_radius      = displaced_position.Norml2();

            MFEM_VERIFY(
                std::isfinite(computational_radius) && computational_radius > 0.0,
                "Monopole H(div) projection encountered an invalid displaced "
                "computational radius."
            );

            const double displacement_determinant = context.displacement_jacobian.Det();

            MFEM_VERIFY(
                std::isfinite(displacement_determinant) && displacement_determinant > 0.0,
                "Monopole H(div) projection encountered an invalid "
                "displacement "
                "Jacobian determinant."
            );

            mfem::DenseMatrix inverse_displacement_jacobian;

            inverse_displacement_jacobian.SetSize(
                context.displacement_jacobian.Height(), context.displacement_jacobian.Width()
            );

            mfem::CalcInverse(context.displacement_jacobian, inverse_displacement_jacobian);

            /*
             * Pull the radial field back only through the regular displacement
             * map.
             *
             * In the compactified vacuum, the Kelvin scale and its radial
             * derivative cancel exactly from the three-dimensional H(div) Piola
             * pullback of the inverse-square monopole field:
             *
             *   det(J) J^{-1} (GM x / |x|^3)
             *     = GM det(A) A^{-1} y / |y|^3.
             *
             * This is also the finite reference-space limit at compactified
             * infinity.
             */
            inverse_displacement_jacobian.Mult(displaced_position, value);

            double radial_denominator = 0.0;

            if (transformation.Attribute == m_vacuum_attribute) {
                radial_denominator = computational_radius * computational_radius * computational_radius;
            } else {
                radial_denominator = m_stellar_radius * m_stellar_radius * m_stellar_radius;
            }

            value *= mean_field::utils::G * m_mass * displacement_determinant / radial_denominator;

            for (int component = 0; component < value.Size(); ++component) {
                MFEM_VERIFY(
                    std::isfinite(value(component)), "Monopole H(div) projection produced a non-finite "
                                                     "reference flux."
                );
            }
        }

    private:
        StatelessProjectionGeometry m_geometry;
        int m_vacuum_attribute;
        double m_mass;
        double m_stellar_radius;
    };

} // namespace

TEST_CASE(
    "Uniform Potential Matches Analytic",
    tags::gravity_analytic_initialization
) {
    auto args  = test_utils::setup_args();
    fem::FEM f = fem::setup_fem(args.mesh_file, args, 0);
    f.mapping->ResetDisplacement();
    physics::update_stiffness_matrix(f);

    const double radius          = utils::RADIUS;
    const double mass            = utils::MASS;
    const double analytic_volume = (4.0 / 3.0) * M_PI * std::pow(radius, 3.0);
    const double density         = mass / analytic_volume;

    mfem::GridFunction rho_uniform(f.densityFes.get());
    rho_uniform = density;
    zero_vacuum_density(f, rho_uniform);
    analysis::conserve_mass(f, rho_uniform, mass);

    f.com                                = analysis::get_com(f, rho_uniform);
    f.Q                                  = physics::compute_quadrupole_moment_tensor(f, rho_uniform, f.com);

    const auto gravity_solution          = physics::grav_potential(f, args, rho_uniform);
    constexpr double potential_tolerance = utils::APPROX_MAX_ACCEPTABLE_POTENTIAL_ERROR_SI_BURNING;
    double local_max_abs_error           = 0.0;
    double local_max_rel_error           = 0.0;

    const int num_elements_to_test       = std::min(30, f.mesh->GetNE());
    for (int elem_id = 0; elem_id < num_elements_to_test; ++elem_id) {
        mfem::ElementTransformation *transformation     = f.mesh->GetElementTransformation(elem_id);
        const mfem::IntegrationRule &integration_rule   = mfem::IntRules.Get(transformation->GetGeometryType(), 2);
        const mfem::IntegrationPoint &integration_point = integration_rule.IntPoint(0);
        transformation->SetIntPoint(&integration_point);

        mfem::Vector x_physical;
        f.mapping->GetPhysicalPoint(*transformation, integration_point, x_physical);

        const double radial_coordinate = x_physical.Norml2();
        if (radial_coordinate < 1.0e-9) {
            continue;
        }

        const double phi_analytic   = -(utils::G * mass / (2.0 * std::pow(radius, 3.0))) *
                                      (3.0 * radius * radius - radial_coordinate * radial_coordinate);
        const double phi_fem        = gravity_solution.phi.GetValue(elem_id, integration_point);
        const double absolute_error = std::abs(phi_fem - phi_analytic);
        const double relative_error = absolute_error / std::abs(phi_analytic);

        local_max_abs_error         = std::max(local_max_abs_error, absolute_error);
        local_max_rel_error         = std::max(local_max_rel_error, relative_error);
        CHECK_THAT(relative_error, Catch::Matchers::WithinAbs(0.0, 0.1 * potential_tolerance));
    }

    double global_max_abs_error = 0.0;
    double global_max_rel_error = 0.0;
    MPI_Comm communicator       = f.densityFes->GetComm();
    MPI_Allreduce(&local_max_abs_error, &global_max_abs_error, 1, MPI_DOUBLE, MPI_MAX, communicator);
    MPI_Allreduce(&local_max_rel_error, &global_max_rel_error, 1, MPI_DOUBLE, MPI_MAX, communicator);

    const int quadrature_order = get_gravity_quadrature_order(f);
    const GravitationalEnergies energies =
        compute_gravitational_energies(f, rho_uniform, gravity_solution, quadrature_order);
    const double analytic_binding_energy = -(3.0 / 5.0) * utils::G * mass * mass / radius;
    const double relative_binding_error =
        std::abs(energies.binding - analytic_binding_energy) / std::abs(analytic_binding_energy);
    const double relative_virial_error =
        std::abs(energies.virial - analytic_binding_energy) / std::abs(analytic_binding_energy);
    const double relative_consistency_error = std::abs(energies.binding - energies.virial) / std::abs(energies.binding);

    INFO("Analytic binding energy = " << analytic_binding_energy);
    INFO("Computed binding energy = " << energies.binding);
    INFO("Computed virial energy = " << energies.virial);
    INFO("Relative virial consistency error = " << relative_consistency_error);

    constexpr double energy_tolerance      = 1.0e-5;
    constexpr double consistency_tolerance = 1.0e-6;

    CHECK_THAT(global_max_rel_error, Catch::Matchers::WithinAbs(0.0, 0.1 * potential_tolerance));
    CHECK_THAT(global_max_abs_error, Catch::Matchers::WithinAbs(0.0, 0.1 * potential_tolerance));
    CHECK_THAT(relative_binding_error, Catch::Matchers::WithinAbs(0.0, energy_tolerance));
    CHECK_THAT(relative_virial_error, Catch::Matchers::WithinAbs(0.0, energy_tolerance));
    CHECK_THAT(relative_consistency_error, Catch::Matchers::WithinAbs(0.0, consistency_tolerance));
}

TEST_CASE(
    "Parabolic Density Virial Self-Consistency",
    tags::gravity_analytic_initialization
) {
    auto args  = test_utils::setup_args();
    fem::FEM f = fem::setup_fem(args.mesh_file, args, 0);
    f.mapping->ResetDisplacement();
    physics::update_stiffness_matrix(f);

    const double radius          = utils::RADIUS;
    const double mass            = utils::MASS;
    const double central_density = (15.0 * mass) / (8.0 * M_PI * std::pow(radius, 3.0));

    auto parabolic_rho           = [central_density, radius](const mfem::Vector &x) {
        const double radial_coordinate = x.Norml2();
        return central_density * (1.0 - radial_coordinate * radial_coordinate / (radius * radius));
    };

    std::unique_ptr<mfem::Coefficient> rho_coeff;
    if (f.has_mapping()) {
        rho_coeff = std::make_unique<mapping::PhysicalPositionFunctionCoefficient>(*f.mapping, parabolic_rho);
    } else {
        rho_coeff = std::make_unique<mfem::FunctionCoefficient>(parabolic_rho);
    }

    mfem::GridFunction rho_grid(f.densityFes.get());
    rho_grid.ProjectCoefficient(*rho_coeff);
    zero_vacuum_density(f, rho_grid);
    analysis::conserve_mass(f, rho_grid, mass);

    f.com                       = analysis::get_com(f, rho_grid);
    f.Q                         = physics::compute_quadrupole_moment_tensor(f, rho_grid, f.com);

    const auto gravity_solution = physics::grav_potential(f, args, rho_grid);
    const int quadrature_order  = get_gravity_quadrature_order(f);
    const GravitationalEnergies energies =
        compute_gravitational_energies(f, rho_grid, gravity_solution, quadrature_order);
    const double analytic_binding_energy = -(5.0 / 7.0) * utils::G * mass * mass / radius;
    const double relative_binding_error =
        std::abs(energies.binding - analytic_binding_energy) / std::abs(analytic_binding_energy);
    const double relative_virial_error =
        std::abs(energies.virial - analytic_binding_energy) / std::abs(analytic_binding_energy);
    const double relative_consistency_error = std::abs(energies.binding - energies.virial) / std::abs(energies.binding);

    INFO("Analytic binding energy = " << analytic_binding_energy);
    INFO("Computed binding energy = " << energies.binding);
    INFO("Computed virial energy = " << energies.virial);
    INFO("Relative virial consistency error = " << relative_consistency_error);

    constexpr double analytic_tolerance    = 1.0e-5;
    constexpr double consistency_tolerance = 1.0e-6;

    CHECK_THAT(relative_binding_error, Catch::Matchers::WithinAbs(0.0, analytic_tolerance));
    CHECK_THAT(relative_virial_error, Catch::Matchers::WithinAbs(0.0, analytic_tolerance));
    CHECK_THAT(relative_consistency_error, Catch::Matchers::WithinAbs(0.0, consistency_tolerance));
}

TEST_CASE(
    "Rational Density Virial Self-Consistency",
    tags::gravity_consistency_initialization
) {
    auto args  = test_utils::setup_args();
    fem::FEM f = fem::setup_fem(args.mesh_file, args, 0);
    f.mapping->ResetDisplacement();
    physics::update_stiffness_matrix(f);

    const double radius            = utils::RADIUS;
    const double mass              = utils::MASS;

    // Larger values are more centrally concentrated and generally harder for a
    // polynomial to represent. A regression would be considered if this test
    // does not pass for concentrations <= 16.0.
    constexpr double concentration = 16.0;
    const double density_scale     = mass / std::pow(radius, 3.0);

    auto rational_rho              = [radius, density_scale](const mfem::Vector &x) {
        const double normalized_radius_squared = (x * x) / (radius * radius);
        if (normalized_radius_squared >= 1.0) {
            return 0.0;
        }

        const double denominator = 1.0 + concentration * normalized_radius_squared;
        return density_scale * (1.0 - normalized_radius_squared) / (denominator * denominator);
    };

    std::unique_ptr<mfem::Coefficient> rho_coeff;
    if (f.has_mapping()) {
        rho_coeff = std::make_unique<mapping::PhysicalPositionFunctionCoefficient>(*f.mapping, rational_rho);
    } else {
        rho_coeff = std::make_unique<mfem::FunctionCoefficient>(rational_rho);
    }

    mfem::GridFunction rho_grid(f.densityFes.get());
    rho_grid.ProjectCoefficient(*rho_coeff);
    zero_vacuum_density(f, rho_grid);
    analysis::conserve_mass(f, rho_grid, mass);

    f.com                       = analysis::get_com(f, rho_grid);
    f.Q                         = physics::compute_quadrupole_moment_tensor(f, rho_grid, f.com);

    const auto gravity_solution = physics::grav_potential(f, args, rho_grid);
    const int quadrature_order  = get_gravity_quadrature_order(f);
    const GravitationalEnergies energies =
        compute_gravitational_energies(f, rho_grid, gravity_solution, quadrature_order);

    REQUIRE(energies.binding < 0.0);
    REQUIRE(energies.virial < 0.0);

    const double relative_consistency_error = std::abs(energies.binding - energies.virial) / std::abs(energies.binding);
    INFO("W_bind = " << energies.binding);
    INFO("W_vir = " << energies.virial);
    INFO("Relative virial consistency error = " << relative_consistency_error);

    constexpr double virial_tolerance = 1.0e-5;
    CHECK_THAT(relative_consistency_error, Catch::Matchers::WithinAbs(0.0, virial_tolerance));
}

TEST_CASE(
    "Homogeneous Ellipsoid Analytic Gravity",
    tags::gravity_analytic_initialization
) {
    auto args                = test_utils::setup_args();
    fem::FEM f               = fem::setup_fem(args.mesh_file, args, 0);

    const double radius      = utils::RADIUS;
    const double mass        = utils::MASS;
    constexpr double x_scale = 1.15;
    constexpr double y_scale = 0.95;
    constexpr double z_scale = 1.0 / (x_scale * y_scale);
    assert(std::abs(x_scale * y_scale * z_scale - 1.0) < 1.0e-14);

    const double semi_axis_x = x_scale * radius;
    const double semi_axis_y = y_scale * radius;
    const double semi_axis_z = z_scale * radius;

    auto affine_displacement = [](const mfem::Vector &x, mfem::Vector &displacement_value) {
        displacement_value.SetSize(3);
        displacement_value(0) = (x_scale - 1.0) * x(0);
        displacement_value(1) = (y_scale - 1.0) * x(1);
        displacement_value(2) = (z_scale - 1.0) * x(2);
    };

    mfem::VectorFunctionCoefficient displacement_coeff(3, affine_displacement);
    mfem::ParGridFunction displacement(f.displacementFes.get());
    displacement.ProjectCoefficient(displacement_coeff);
    f.mapping->SetDisplacement(displacement);
    physics::update_stiffness_matrix(f);

    const double analytic_volume = (4.0 / 3.0) * M_PI * semi_axis_x * semi_axis_y * semi_axis_z;
    const double density         = mass / analytic_volume;

    mfem::GridFunction rho_grid(f.densityFes.get());
    rho_grid = density;
    zero_vacuum_density(f, rho_grid);

    const double projected_mass    = analysis::domain_integrate_grid_function(f, rho_grid, utils::DOMAINS::STELLAR);
    const double numerical_density = density * mass / projected_mass;
    analysis::conserve_mass(f, rho_grid, mass);

    f.com = analysis::get_com(f, rho_grid);
    f.Q   = physics::compute_quadrupole_moment_tensor(f, rho_grid, f.com);

    const HomogeneousEllipsoidAnalytic analytic =
        compute_homogeneous_ellipsoid_analytic(semi_axis_x, semi_axis_y, semi_axis_z);
    const double coefficient_sum = analytic.coefficient_x + analytic.coefficient_y + analytic.coefficient_z;

    INFO("A_x = " << analytic.coefficient_x);
    INFO("A_y = " << analytic.coefficient_y);
    INFO("A_z = " << analytic.coefficient_z);
    INFO("A_x + A_y + A_z = " << coefficient_sum);
    REQUIRE_THAT(coefficient_sum, Catch::Matchers::WithinAbs(2.0, 1.0e-11));

    mfem::DenseMatrix analytic_quadrupole(3, 3);
    analytic_quadrupole = 0.0;
    analytic_quadrupole(0, 0) =
        (mass / 5.0) * (2.0 * semi_axis_x * semi_axis_x - semi_axis_y * semi_axis_y - semi_axis_z * semi_axis_z);
    analytic_quadrupole(1, 1) =
        (mass / 5.0) * (2.0 * semi_axis_y * semi_axis_y - semi_axis_x * semi_axis_x - semi_axis_z * semi_axis_z);
    analytic_quadrupole(2, 2) =
        (mass / 5.0) * (2.0 * semi_axis_z * semi_axis_z - semi_axis_x * semi_axis_x - semi_axis_y * semi_axis_y);

    mfem::DenseMatrix quadrupole_difference(f.Q);
    quadrupole_difference -= analytic_quadrupole;
    const double relative_quadrupole_error = quadrupole_difference.FNorm() / analytic_quadrupole.FNorm();
    INFO("Relative quadrupole error = " << relative_quadrupole_error);

    HomogeneousEllipsoidHDivCoefficient analytic_field_coefficient(*f.mapping, numerical_density, analytic);
    mfem::ParGridFunction analytic_field_projection(f.gravityFluxFes.get());
    analytic_field_projection = 0.0;

    for (int i = 0; i < f.mesh->attributes.Size(); ++i) {
        const int attribute = f.mesh->attributes[i];

        if (attribute != 3) {
            analytic_field_projection.ProjectCoefficient(analytic_field_coefficient, attribute);
        }
    }

    const auto gravity_solution                        = physics::grav_potential(f, args, rho_grid);
    const int quadrature_order                         = get_gravity_quadrature_order(f);
    double local_field_error_squared                   = 0.0;
    double local_field_norm_squared                    = 0.0;
    double local_projection_error_squared              = 0.0;
    double local_gravity_projection_difference_squared = 0.0;

    mfem::Vector x_physical(3);
    mfem::Vector grad_phi_element(3);
    mfem::Vector grad_phi_physical(3);
    mfem::Vector grad_phi_analytic(3);
    mfem::Vector grad_phi_difference(3);
    mfem::Vector projected_field_element(3);
    mfem::Vector projected_field_physical(3);
    mfem::Vector projected_field_difference(3);
    mfem::DenseMatrix map_jacobian(3, 3);
    mfem::Vector gravity_projection_difference(3);

    for (int elem_id = 0; elem_id < f.mesh->GetNE(); ++elem_id) {
        if (f.mesh->GetAttribute(elem_id) == 3) {
            continue;
        }

        mfem::ElementTransformation *transformation = f.mesh->GetElementTransformation(elem_id);
        const mfem::IntegrationRule &integration_rule =
            mfem::IntRules.Get(transformation->GetGeometryType(), quadrature_order);

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

            const double map_determinant = f.mapping->ComputeDetJ(*transformation, integration_point);
            MFEM_VERIFY(map_determinant > 0.0, "Domain mapping has a non-positive Jacobian determinant.");

            const double weight = transformation->Weight() * integration_point.weight * map_determinant;
            f.mapping->GetPhysicalPoint(*transformation, integration_point, x_physical);
            gravity_solution.gradPhi.GetVectorValue(elem_id, integration_point, grad_phi_element);
            f.mapping->ComputeJacobian(*transformation, map_jacobian);
            map_jacobian.Mult(grad_phi_element, grad_phi_physical);
            grad_phi_physical /= map_determinant;

            analytic_field_projection.GetVectorValue(elem_id, integration_point, projected_field_element);
            map_jacobian.Mult(projected_field_element, projected_field_physical);
            projected_field_physical /= map_determinant;

            grad_phi_analytic(0) = 2.0 * M_PI * utils::G * numerical_density * analytic.coefficient_x * x_physical(0);
            grad_phi_analytic(1) = 2.0 * M_PI * utils::G * numerical_density * analytic.coefficient_y * x_physical(1);
            grad_phi_analytic(2) = 2.0 * M_PI * utils::G * numerical_density * analytic.coefficient_z * x_physical(2);

            projected_field_difference = projected_field_physical;
            projected_field_difference -= grad_phi_analytic;
            local_projection_error_squared += (projected_field_difference * projected_field_difference) * weight;

            gravity_projection_difference = grad_phi_physical;
            gravity_projection_difference -= projected_field_physical;
            local_gravity_projection_difference_squared +=
                (gravity_projection_difference * gravity_projection_difference) * weight;

            grad_phi_difference = grad_phi_physical;
            grad_phi_difference -= grad_phi_analytic;
            local_field_error_squared += (grad_phi_difference * grad_phi_difference) * weight;
            local_field_norm_squared += (grad_phi_analytic * grad_phi_analytic) * weight;
        }
    }

    double global_field_error_squared                   = 0.0;
    double global_field_norm_squared                    = 0.0;
    double global_projection_error_squared              = 0.0;
    double global_gravity_projection_difference_squared = 0.0;
    MPI_Comm communicator                               = f.densityFes->GetComm();
    MPI_Allreduce(&local_field_error_squared, &global_field_error_squared, 1, MPI_DOUBLE, MPI_SUM, communicator);
    MPI_Allreduce(&local_field_norm_squared, &global_field_norm_squared, 1, MPI_DOUBLE, MPI_SUM, communicator);
    MPI_Allreduce(
        &local_projection_error_squared, &global_projection_error_squared, 1, MPI_DOUBLE, MPI_SUM, communicator
    );
    MPI_Allreduce(
        &local_gravity_projection_difference_squared, &global_gravity_projection_difference_squared, 1, MPI_DOUBLE,
        MPI_SUM, communicator
    );

    const double relative_field_error = std::sqrt(global_field_error_squared / global_field_norm_squared);
    const GravitationalEnergies energies =
        compute_gravitational_energies(f, rho_grid, gravity_solution, quadrature_order);
    const double analytic_binding_energy = -(3.0 / 10.0) * utils::G * mass * mass * analytic.energy_kernel;
    const double relative_binding_energy_error =
        std::abs(energies.binding - analytic_binding_energy) / std::abs(analytic_binding_energy);
    const double relative_virial_energy_error =
        std::abs(energies.virial - analytic_binding_energy) / std::abs(analytic_binding_energy);
    const double relative_consistency_error = std::abs(energies.binding - energies.virial) / std::abs(energies.binding);
    const double relative_projection_error  = std::sqrt(global_projection_error_squared / global_field_norm_squared);
    const double gravity_to_projection_error_ratio = relative_projection_error > 0.0
                                                         ? relative_field_error / relative_projection_error
                                                         : std::numeric_limits<double>::infinity();
    const double relative_gravity_projection_difference =
        std::sqrt(global_gravity_projection_difference_squared / global_field_norm_squared);
    const double projection_gap_ratio = relative_gravity_projection_difference / relative_projection_error;

    INFO("Analytic binding energy = " << analytic_binding_energy);
    INFO("Computed binding energy = " << energies.binding);
    INFO("Computed virial energy = " << energies.virial);
    INFO("Relative field L2 error = " << relative_field_error);
    INFO("Relative binding energy error = " << relative_binding_energy_error);
    INFO("Relative virial energy error = " << relative_virial_energy_error);
    INFO("Relative virial consistency error = " << relative_consistency_error);
    INFO("Relative gravity-to-RT-projection difference = " << relative_gravity_projection_difference);
    INFO("Gravity-to-projection gap ratio = " << projection_gap_ratio);

    INFO("Relative RT projection L2 error = " << relative_projection_error);
    INFO("Gravity-to-projection error ratio = " << gravity_to_projection_error_ratio);
    REQUIRE(std::isfinite(relative_projection_error));

    constexpr double quadrupole_tolerance  = 2.0e-4;
    constexpr double field_tolerance       = 1.0e-5;
    constexpr double energy_tolerance      = 1.0e-5;
    constexpr double consistency_tolerance = 1.0e-5;

    CHECK_THAT(relative_quadrupole_error, Catch::Matchers::WithinAbs(0.0, quadrupole_tolerance));
    CHECK_THAT(relative_field_error, Catch::Matchers::WithinAbs(0.0, field_tolerance));
    CHECK_THAT(relative_binding_energy_error, Catch::Matchers::WithinAbs(0.0, energy_tolerance));
    CHECK_THAT(relative_virial_energy_error, Catch::Matchers::WithinAbs(0.0, energy_tolerance));
    CHECK_THAT(relative_consistency_error, Catch::Matchers::WithinAbs(0.0, consistency_tolerance));
}

TEST_CASE(
    "Deformed Rational Density Virial Self-Consistency",
    tags::gravity_consistency_initialization
) {
    auto args                = test_utils::setup_args();
    fem::FEM f               = fem::setup_fem(args.mesh_file, args, 0);

    const double radius      = utils::RADIUS;
    const double mass        = utils::MASS;
    constexpr double x_scale = 1.15;
    constexpr double y_scale = 0.95;
    constexpr double z_scale = 1.0 / (x_scale * y_scale);
    assert(std::abs(x_scale * y_scale * z_scale - 1.0) < 1.0e-14);

    const double semi_axis_x = x_scale * radius;
    const double semi_axis_y = y_scale * radius;
    const double semi_axis_z = z_scale * radius;

    auto affine_displacement = [](const mfem::Vector &x, mfem::Vector &displacement_value) {
        displacement_value.SetSize(3);
        displacement_value(0) = (x_scale - 1.0) * x(0);
        displacement_value(1) = (y_scale - 1.0) * x(1);
        displacement_value(2) = (z_scale - 1.0) * x(2);
    };

    mfem::VectorFunctionCoefficient displacement_coeff(3, affine_displacement);
    mfem::ParGridFunction displacement(f.displacementFes.get());
    displacement.ProjectCoefficient(displacement_coeff);
    f.mapping->SetDisplacement(displacement);
    physics::update_stiffness_matrix(f);

    constexpr double concentration = 16.0;
    const double density_scale     = mass / std::pow(radius, 3.0);

    auto ellipsoidal_rho           = [semi_axis_x, semi_axis_y, semi_axis_z, density_scale](const mfem::Vector &x) {
        const double ellipsoidal_radius_squared = x(0) * x(0) / (semi_axis_x * semi_axis_x) +
                                                  x(1) * x(1) / (semi_axis_y * semi_axis_y) +
                                                  x(2) * x(2) / (semi_axis_z * semi_axis_z);

        if (ellipsoidal_radius_squared >= 1.0) {
            return 0.0;
        }

        const double denominator = 1.0 + concentration * ellipsoidal_radius_squared;
        return density_scale * (1.0 - ellipsoidal_radius_squared) / (denominator * denominator);
    };

    std::unique_ptr<mfem::Coefficient> rho_coeff;
    if (f.has_mapping()) {
        rho_coeff = std::make_unique<mapping::PhysicalPositionFunctionCoefficient>(*f.mapping, ellipsoidal_rho);
    } else {
        rho_coeff = std::make_unique<mfem::FunctionCoefficient>(ellipsoidal_rho);
    }

    mfem::GridFunction rho_grid(f.densityFes.get());
    rho_grid.ProjectCoefficient(*rho_coeff);
    zero_vacuum_density(f, rho_grid);
    analysis::conserve_mass(f, rho_grid, mass);

    f.com                              = analysis::get_com(f, rho_grid);
    f.Q                                = physics::compute_quadrupole_moment_tensor(f, rho_grid, f.com);

    const double normalized_quadrupole = f.Q.FNorm() / (mass * radius * radius);
    INFO("Normalized quadrupole = " << normalized_quadrupole);
    REQUIRE(normalized_quadrupole > 1.0e-3);

    const auto gravity_solution = physics::grav_potential(f, args, rho_grid);
    const int quadrature_order  = get_gravity_quadrature_order(f);
    const GravitationalEnergies energies =
        compute_gravitational_energies(f, rho_grid, gravity_solution, quadrature_order);

    REQUIRE(energies.binding < 0.0);
    REQUIRE(energies.virial < 0.0);

    const double relative_consistency_error = std::abs(energies.binding - energies.virial) / std::abs(energies.binding);
    INFO("W_bind = " << energies.binding);
    INFO("W_vir = " << energies.virial);
    INFO("Relative virial consistency error = " << relative_consistency_error);

    constexpr double virial_tolerance = 1.0e-5;
    CHECK_THAT(relative_consistency_error, Catch::Matchers::WithinAbs(0.0, virial_tolerance));
}

TEST_CASE(
    "New Gravity Potential Matches Uniform Sphere Analytic",
    tags::gravity_analytic
) {
    auto args        = test_utils::setup_args();
    args.p.rtol      = 1.0e-13;
    args.p.max_iters = std::max(args.p.max_iters, 1000);

    fem::FEM f       = fem::setup_fem(args.mesh_file, args, 0);

    f.mapping->ResetDisplacement();
    physics::update_stiffness_matrix(f);

    mfem::ParGridFunction displacement(f.displacementFes.get());
    displacement                 = 0.0;

    const double radius          = utils::RADIUS;
    const double mass            = utils::MASS;
    const double analytic_volume = (4.0 / 3.0) * M_PI * std::pow(radius, 3.0);
    const double density         = mass / analytic_volume;

    mfem::GridFunction rho_uniform(f.densityFes.get());
    rho_uniform = density;
    zero_vacuum_density(f, rho_uniform);
    analysis::conserve_mass(f, rho_uniform, mass);

    f.com                                           = analysis::get_com(f, rho_uniform);
    f.Q                                             = physics::compute_quadrupole_moment_tensor(f, rho_uniform, f.com);

    const physics::GravitySolution gravity_solution = physics::grav_potential_new(f, args, rho_uniform, displacement);

    constexpr double potential_tolerance            = utils::APPROX_MAX_ACCEPTABLE_POTENTIAL_ERROR_SI_BURNING;
    double local_maximum_absolute_error             = 0.0;
    double local_maximum_relative_error             = 0.0;

    const int elements_to_test                      = std::min(30, f.mesh->GetNE());

    for (int element_id = 0; element_id < elements_to_test; ++element_id) {
        if (f.mesh->GetAttribute(element_id) == 3) {
            continue;
        }

        mfem::ElementTransformation *transformation     = f.mesh->GetElementTransformation(element_id);
        const mfem::IntegrationRule &integration_rule   = mfem::IntRules.Get(transformation->GetGeometryType(), 2);
        const mfem::IntegrationPoint &integration_point = integration_rule.IntPoint(0);
        transformation->SetIntPoint(&integration_point);

        mfem::Vector physical_position;
        f.mapping->GetPhysicalPoint(*transformation, integration_point, physical_position);

        const double radial_coordinate = physical_position.Norml2();

        if (radial_coordinate < 1.0e-9) {
            continue;
        }

        const double analytic_potential = -(utils::G * mass / (2.0 * std::pow(radius, 3.0))) *
                                          (3.0 * radius * radius - radial_coordinate * radial_coordinate);

        const double computed_potential = gravity_solution.phi.GetValue(element_id, integration_point);
        const double absolute_error     = std::abs(computed_potential - analytic_potential);
        const double relative_error     = absolute_error / std::abs(analytic_potential);

        local_maximum_absolute_error    = std::max(local_maximum_absolute_error, absolute_error);
        local_maximum_relative_error    = std::max(local_maximum_relative_error, relative_error);
    }

    double global_maximum_absolute_error = 0.0;
    double global_maximum_relative_error = 0.0;

    MPI_Allreduce(
        &local_maximum_absolute_error, &global_maximum_absolute_error, 1, MPI_DOUBLE, MPI_MAX, f.densityFes->GetComm()
    );

    MPI_Allreduce(
        &local_maximum_relative_error, &global_maximum_relative_error, 1, MPI_DOUBLE, MPI_MAX, f.densityFes->GetComm()
    );

    const int quadrature_order = get_gravity_quadrature_order(f);
    const GravitationalEnergies energies =
        compute_gravitational_energies(f, rho_uniform, gravity_solution, quadrature_order);
    const double analytic_binding_energy = -(3.0 / 5.0) * utils::G * mass * mass / radius;
    const double relative_binding_error =
        std::abs(energies.binding - analytic_binding_energy) / std::abs(analytic_binding_energy);
    const double relative_virial_error =
        std::abs(energies.virial - analytic_binding_energy) / std::abs(analytic_binding_energy);
    const double relative_consistency_error = std::abs(energies.binding - energies.virial) / std::abs(energies.binding);

    INFO("Global maximum absolute potential error = " << global_maximum_absolute_error);
    INFO("Global maximum relative potential error = " << global_maximum_relative_error);
    INFO("Analytic binding energy = " << analytic_binding_energy);
    INFO("New-solver binding energy = " << energies.binding);
    INFO("New-solver virial energy = " << energies.virial);
    INFO("Relative binding-energy error = " << relative_binding_error);
    INFO("Relative virial-energy error = " << relative_virial_error);
    INFO("Relative virial consistency error = " << relative_consistency_error);

    REQUIRE(energies.binding < 0.0);
    REQUIRE(energies.virial < 0.0);

    constexpr double energy_tolerance      = 1.0e-5;
    constexpr double consistency_tolerance = 1.0e-6;

    CHECK_THAT(global_maximum_relative_error, Catch::Matchers::WithinAbs(0.0, 0.1 * potential_tolerance));
    CHECK_THAT(global_maximum_absolute_error, Catch::Matchers::WithinAbs(0.0, 0.1 * potential_tolerance));
    CHECK_THAT(relative_binding_error, Catch::Matchers::WithinAbs(0.0, energy_tolerance));
    CHECK_THAT(relative_virial_error, Catch::Matchers::WithinAbs(0.0, energy_tolerance));
    CHECK_THAT(relative_consistency_error, Catch::Matchers::WithinAbs(0.0, consistency_tolerance));
}

TEST_CASE(
    "New Gravity Potential Matches Legacy Solver On Homogeneous Ellipsoid",
    tags::gravity_legacy
) {
    auto args                = test_utils::setup_args();
    args.p.rtol              = 1.0e-13;
    args.p.max_iters         = std::max(args.p.max_iters, 1000);

    fem::FEM f               = fem::setup_fem(args.mesh_file, args, 0);

    const double radius      = utils::RADIUS;
    const double mass        = utils::MASS;

    constexpr double x_scale = 1.15;
    constexpr double y_scale = 0.95;
    constexpr double z_scale = 1.0 / (x_scale * y_scale);

    REQUIRE_THAT(x_scale * y_scale * z_scale, Catch::Matchers::WithinAbs(1.0, 1.0e-14));

    const double semi_axis_x = x_scale * radius;
    const double semi_axis_y = y_scale * radius;
    const double semi_axis_z = z_scale * radius;

    auto affine_displacement = [](const mfem::Vector &position, mfem::Vector &value) {
        value.SetSize(3);
        value(0) = (x_scale - 1.0) * position(0);
        value(1) = (y_scale - 1.0) * position(1);
        value(2) = (z_scale - 1.0) * position(2);
    };

    mfem::VectorFunctionCoefficient displacement_coefficient(3, affine_displacement);
    mfem::ParGridFunction displacement(f.displacementFes.get());
    displacement.ProjectCoefficient(displacement_coefficient);

    f.mapping->SetDisplacement(displacement);
    physics::update_stiffness_matrix(f);

    REQUIRE(f.gravityContext.block_A != nullptr);
    REQUIRE(f.gravityContext.source_form != nullptr);
    REQUIRE(f.domainMapperStateless != nullptr);

    const double analytic_volume = (4.0 / 3.0) * M_PI * semi_axis_x * semi_axis_y * semi_axis_z;
    const double density         = mass / analytic_volume;

    mfem::GridFunction rho_grid(f.densityFes.get());
    rho_grid = density;
    zero_vacuum_density(f, rho_grid);

    const double projected_mass    = analysis::domain_integrate_grid_function(f, rho_grid, utils::DOMAINS::STELLAR);
    const double numerical_density = density * mass / projected_mass;
    analysis::conserve_mass(f, rho_grid, mass);

    f.com                                          = analysis::get_com(f, rho_grid);
    f.Q                                            = physics::compute_quadrupole_moment_tensor(f, rho_grid, f.com);

    const physics::GravitySolution legacy_solution = physics::grav_potential(f, args, rho_grid);
    const physics::GravitySolution new_solution    = physics::grav_potential_new(f, args, rho_grid, displacement);

    using gravity_form                             = utils::blocks::gravity_field_form;

    constexpr auto gravity_gradient_residual_block =
        utils::blocks::get_residual_block<gravity_form>(utils::blocks::gravity_field.gradient_term);

    constexpr auto gravity_poisson_residual_block =
        utils::blocks::get_residual_block<gravity_form>(utils::blocks::gravity_field.poisson_term);

    using DomainSchema                   = utils::domain::CoreEnvelopeVacuumDomainSchema;
    const field::FieldDofMap density_map = field::make_field_dof_map<field::Density, DomainSchema>(*f.densityFes);
    const field::FieldDofMap displacement_map =
        field::make_field_dof_map<field::Displacement, DomainSchema>(*f.displacementFes);
    const field::FieldDofMap gravity_flux_map =
        field::make_field_dof_map<field::Gravity, DomainSchema>(*f.gravityFluxFes);
    const field::FieldDofMap gravity_potential_map =
        field::make_field_dof_map<field::Gravity, DomainSchema>(*f.gravityPotentialFes);

    const std::array<int, gravity_form::value_block_count> value_sizes{
        density_map.reduced_size(), displacement_map.reduced_size(), gravity_flux_map.reduced_size(),
        gravity_potential_map.reduced_size()
    };

    const std::array<int, gravity_form::residual_block_count> residual_sizes{
        gravity_flux_map.reduced_size(), gravity_potential_map.reduced_size()
    };

    const utils::blocks::form_layout<gravity_form> gravity_layout(value_sizes, residual_sizes);
    const int gravity_system_size   = gravity_layout.residual_offsets().Last();
    const int gravity_gradient_size = gravity_layout.size(gravity_gradient_residual_block);
    const int gravity_poisson_size  = gravity_layout.size(gravity_poisson_residual_block);

    mfem::Vector density_true;
    mfem::Vector displacement_true;
    gravity_test_get_true_dofs(*f.densityFes, rho_grid, density_true);
    gravity_test_get_true_dofs(*f.displacementFes, displacement, displacement_true);

    REQUIRE(density_true.Size() == f.densityFes->GetTrueVSize());
    REQUIRE(displacement_true.Size() == f.displacementFes->GetTrueVSize());
    const mfem::Vector reduced_density      = density_map.gather(density_true);
    const mfem::Vector reduced_displacement = displacement_map.gather(displacement_true);

    operators::context::gravity_field::GravityFieldLinearizationContext residual_linearization_context(
        f, *f.domainMapperStateless
    );

    operators::GravityFieldJacobianOperator residual_jacobian(
        f, *f.domainMapperStateless, residual_linearization_context, gravity_layout.value_offsets(),
        gravity_layout.residual_offsets()
    );

    operators::GravityFieldOperator residual_gravity_operator(
        f, *f.domainMapperStateless, residual_linearization_context, gravity_layout.value_offsets(), residual_jacobian
    );

    operators::context::gravity_field::GravityFieldGeometryContext residual_geometry_context(
        f, *f.domainMapperStateless
    );

    operators::ReducedGravityFieldOperator new_reduced_operator(
        residual_gravity_operator, residual_geometry_context, reduced_displacement
    );

    REQUIRE(new_reduced_operator.Width() == gravity_system_size);
    REQUIRE(new_reduced_operator.Height() == gravity_system_size);

    mfem::Vector new_right_hand_side;
    new_reduced_operator.BuildRightHandSide(reduced_density, new_right_hand_side);

    REQUIRE(new_right_hand_side.Size() == gravity_system_size);

    REQUIRE(f.gravityContext.source_form->Width() == reduced_density.Size());
    REQUIRE(f.gravityContext.source_form->Height() == gravity_poisson_size);

    mfem::Vector legacy_source_action(f.gravityContext.source_form->Height());
    legacy_source_action = 0.0;

    f.gravityContext.source_form->Mult(reduced_density, legacy_source_action);

    REQUIRE(legacy_source_action.Size() == gravity_poisson_size);

    mfem::BlockVector legacy_right_hand_side(gravity_layout.residual_offsets());
    legacy_right_hand_side                                          = 0.0;
    legacy_right_hand_side.GetBlock(gravity_poisson_residual_block) = legacy_source_action;

    REQUIRE(legacy_right_hand_side.Size() == gravity_system_size);

    MPI_Comm communicator                    = f.densityFes->GetComm();

    const double new_right_hand_side_norm    = gravity_test_global_norm(new_right_hand_side, communicator);
    const double legacy_right_hand_side_norm = gravity_test_global_norm(legacy_right_hand_side, communicator);

    REQUIRE(new_right_hand_side_norm > 0.0);
    REQUIRE(legacy_right_hand_side_norm > 0.0);

    mfem::Vector right_hand_side_difference(new_right_hand_side);
    right_hand_side_difference -= legacy_right_hand_side;

    const double right_hand_side_scale =
        std::max(0.5 * (new_right_hand_side_norm + legacy_right_hand_side_norm), std::numeric_limits<double>::min());

    const double relative_right_hand_side_difference =
        gravity_test_global_norm(right_hand_side_difference, communicator) / right_hand_side_scale;

    mfem::BlockVector legacy_gravity_state(gravity_layout.residual_offsets());
    mfem::BlockVector new_gravity_state(gravity_layout.residual_offsets());

    legacy_gravity_state = 0.0;
    new_gravity_state    = 0.0;

    {
        mfem::Vector gradient_true;
        mfem::Vector potential_true;

        legacy_solution.gradPhi.GetTrueDofs(gradient_true);
        legacy_solution.phi.GetTrueDofs(potential_true);

        REQUIRE(gradient_true.Size() == gravity_gradient_size);
        REQUIRE(potential_true.Size() == gravity_poisson_size);

        legacy_gravity_state.GetBlock(gravity_gradient_residual_block) = gradient_true;
        legacy_gravity_state.GetBlock(gravity_poisson_residual_block)  = potential_true;

        new_solution.gradPhi.GetTrueDofs(gradient_true);
        new_solution.phi.GetTrueDofs(potential_true);

        REQUIRE(gradient_true.Size() == gravity_gradient_size);
        REQUIRE(potential_true.Size() == gravity_poisson_size);

        new_gravity_state.GetBlock(gravity_gradient_residual_block) = gradient_true;
        new_gravity_state.GetBlock(gravity_poisson_residual_block)  = potential_true;
    }

    REQUIRE(legacy_gravity_state.Size() == gravity_system_size);
    REQUIRE(new_gravity_state.Size() == gravity_system_size);

    REQUIRE(f.gravityContext.block_A->Width() == gravity_system_size);
    REQUIRE(f.gravityContext.block_A->Height() == gravity_system_size);

    mfem::Vector legacy_action_at_legacy_solution(gravity_system_size);
    mfem::Vector legacy_action_at_new_solution(gravity_system_size);
    mfem::Vector new_action_at_legacy_solution(gravity_system_size);
    mfem::Vector new_action_at_new_solution(gravity_system_size);

    legacy_action_at_legacy_solution = 0.0;
    legacy_action_at_new_solution    = 0.0;
    new_action_at_legacy_solution    = 0.0;
    new_action_at_new_solution       = 0.0;

    f.gravityContext.block_A->Mult(legacy_gravity_state, legacy_action_at_legacy_solution);

    f.gravityContext.block_A->Mult(new_gravity_state, legacy_action_at_new_solution);

    new_reduced_operator.Mult(legacy_gravity_state, new_action_at_legacy_solution);

    new_reduced_operator.Mult(new_gravity_state, new_action_at_new_solution);

    REQUIRE(legacy_action_at_legacy_solution.Size() == gravity_system_size);
    REQUIRE(legacy_action_at_new_solution.Size() == gravity_system_size);
    REQUIRE(new_action_at_legacy_solution.Size() == gravity_system_size);
    REQUIRE(new_action_at_new_solution.Size() == gravity_system_size);

    mfem::Vector legacy_residual_at_legacy_solution(legacy_action_at_legacy_solution);
    mfem::Vector legacy_residual_at_new_solution(legacy_action_at_new_solution);
    mfem::Vector new_residual_at_legacy_solution(new_action_at_legacy_solution);
    mfem::Vector new_residual_at_new_solution(new_action_at_new_solution);

    legacy_residual_at_legacy_solution -= legacy_right_hand_side;
    legacy_residual_at_new_solution -= legacy_right_hand_side;
    new_residual_at_legacy_solution -= new_right_hand_side;
    new_residual_at_new_solution -= new_right_hand_side;

    const GravityResidualMetrics legacy_at_legacy = measure_gravity_residual(
        legacy_residual_at_legacy_solution, gravity_layout.residual_offsets(), legacy_right_hand_side_norm, communicator
    );

    const GravityResidualMetrics legacy_at_new = measure_gravity_residual(
        legacy_residual_at_new_solution, gravity_layout.residual_offsets(), legacy_right_hand_side_norm, communicator
    );

    const GravityResidualMetrics new_at_legacy = measure_gravity_residual(
        new_residual_at_legacy_solution, gravity_layout.residual_offsets(), new_right_hand_side_norm, communicator
    );

    const GravityResidualMetrics new_at_new = measure_gravity_residual(
        new_residual_at_new_solution, gravity_layout.residual_offsets(), new_right_hand_side_norm, communicator
    );

    mfem::Vector operator_gap_at_legacy_solution(new_residual_at_legacy_solution);
    operator_gap_at_legacy_solution -= legacy_residual_at_legacy_solution;

    mfem::Vector operator_gap_at_new_solution(new_residual_at_new_solution);
    operator_gap_at_new_solution -= legacy_residual_at_new_solution;

    const GravityResidualMetrics gap_at_legacy = measure_gravity_residual(
        operator_gap_at_legacy_solution, gravity_layout.residual_offsets(), right_hand_side_scale, communicator
    );

    const GravityResidualMetrics gap_at_new = measure_gravity_residual(
        operator_gap_at_new_solution, gravity_layout.residual_offsets(), right_hand_side_scale, communicator
    );

    INFO("New/legacy right-hand-side difference = " << relative_right_hand_side_difference);
    INFO(
        "Legacy operator at legacy solution: total = " << legacy_at_legacy.relative_total
                                                       << ", gradient = " << legacy_at_legacy.relative_gradient
                                                       << ", Poisson = " << legacy_at_legacy.relative_poisson
    );
    INFO(
        "Legacy operator at new solution: total = " << legacy_at_new.relative_total
                                                    << ", gradient = " << legacy_at_new.relative_gradient
                                                    << ", Poisson = " << legacy_at_new.relative_poisson
    );
    INFO(
        "New operator at legacy solution: total = " << new_at_legacy.relative_total
                                                    << ", gradient = " << new_at_legacy.relative_gradient
                                                    << ", Poisson = " << new_at_legacy.relative_poisson
    );
    INFO(
        "New operator at new solution: total = " << new_at_new.relative_total
                                                 << ", gradient = " << new_at_new.relative_gradient
                                                 << ", Poisson = " << new_at_new.relative_poisson
    );
    INFO(
        "Operator gap at legacy solution: total = " << gap_at_legacy.relative_total
                                                    << ", gradient = " << gap_at_legacy.relative_gradient
                                                    << ", Poisson = " << gap_at_legacy.relative_poisson
    );
    INFO(
        "Operator gap at new solution: total = " << gap_at_new.relative_total
                                                 << ", gradient = " << gap_at_new.relative_gradient
                                                 << ", Poisson = " << gap_at_new.relative_poisson
    );

    REQUIRE(std::isfinite(legacy_at_legacy.relative_total));
    REQUIRE(std::isfinite(legacy_at_new.relative_total));
    REQUIRE(std::isfinite(new_at_legacy.relative_total));
    REQUIRE(std::isfinite(new_at_new.relative_total));
    REQUIRE(std::isfinite(gap_at_legacy.relative_total));
    REQUIRE(std::isfinite(gap_at_new.relative_total));

    constexpr double source_parity_tolerance     = 1.0e-12;
    constexpr double diagonal_residual_tolerance = 1.0e-8;
    constexpr double poisson_gap_tolerance       = 1.0e-11;

    CHECK_THAT(relative_right_hand_side_difference, Catch::Matchers::WithinAbs(0.0, source_parity_tolerance));

    CHECK_THAT(legacy_at_legacy.relative_total, Catch::Matchers::WithinAbs(0.0, diagonal_residual_tolerance));

    CHECK_THAT(new_at_new.relative_total, Catch::Matchers::WithinAbs(0.0, diagonal_residual_tolerance));

    CHECK_THAT(gap_at_legacy.relative_poisson, Catch::Matchers::WithinAbs(0.0, poisson_gap_tolerance));

    CHECK_THAT(gap_at_new.relative_poisson, Catch::Matchers::WithinAbs(0.0, poisson_gap_tolerance));

    const int quadrature_order = get_gravity_quadrature_order(f);

    const GravitySolutionComparison comparison =
        compare_gravity_solutions(f, legacy_solution, new_solution, quadrature_order);

    const GravitationalEnergies legacy_energies =
        compute_gravitational_energies(f, rho_grid, legacy_solution, quadrature_order);

    const GravitationalEnergies new_energies =
        compute_gravitational_energies(f, rho_grid, new_solution, quadrature_order);

    const HomogeneousEllipsoidAnalytic analytic =
        compute_homogeneous_ellipsoid_analytic(semi_axis_x, semi_axis_y, semi_axis_z);

    const double analytic_binding_energy = -(3.0 / 10.0) * utils::G * mass * mass * analytic.energy_kernel;

    const double new_binding_analytic_error =
        std::abs(new_energies.binding - analytic_binding_energy) / std::abs(analytic_binding_energy);

    const double new_virial_analytic_error =
        std::abs(new_energies.virial - analytic_binding_energy) / std::abs(analytic_binding_energy);

    const double new_virial_consistency_error =
        std::abs(new_energies.binding - new_energies.virial) / std::abs(new_energies.binding);

    const double relative_binding_difference =
        std::abs(new_energies.binding - legacy_energies.binding) /
        std::max(
            0.5 * (std::abs(new_energies.binding) + std::abs(legacy_energies.binding)),
            std::numeric_limits<double>::min()
        );

    const double relative_virial_difference =
        std::abs(new_energies.virial - legacy_energies.virial) /
        std::max(
            0.5 * (std::abs(new_energies.virial) + std::abs(legacy_energies.virial)), std::numeric_limits<double>::min()
        );

    INFO("Numerical density = " << numerical_density);
    INFO("Analytic binding energy = " << analytic_binding_energy);
    INFO("Legacy binding energy = " << legacy_energies.binding);
    INFO("New binding energy = " << new_energies.binding);
    INFO("Legacy virial energy = " << legacy_energies.virial);
    INFO("New virial energy = " << new_energies.virial);
    INFO("New binding analytic error = " << new_binding_analytic_error);
    INFO("New virial analytic error = " << new_virial_analytic_error);
    INFO("New virial consistency error = " << new_virial_consistency_error);
    INFO("New/legacy physical gradient difference = " << comparison.relative_gradient_difference);
    INFO("New/legacy physical potential difference = " << comparison.relative_potential_difference);
    INFO("New/legacy binding-energy difference = " << relative_binding_difference);
    INFO("New/legacy virial-energy difference = " << relative_virial_difference);

    REQUIRE(legacy_energies.binding < 0.0);
    REQUIRE(legacy_energies.virial < 0.0);
    REQUIRE(new_energies.binding < 0.0);
    REQUIRE(new_energies.virial < 0.0);

    REQUIRE(std::isfinite(comparison.relative_gradient_difference));
    REQUIRE(std::isfinite(comparison.relative_potential_difference));

    constexpr double gradient_parity_tolerance  = 5.0e-4;
    constexpr double potential_parity_tolerance = 1.0e-4;
    constexpr double energy_parity_tolerance    = 1.0e-5;
    constexpr double analytic_energy_tolerance  = 1.0e-5;
    constexpr double consistency_tolerance      = 1.0e-5;

    CHECK_THAT(comparison.relative_gradient_difference, Catch::Matchers::WithinAbs(0.0, gradient_parity_tolerance));

    CHECK_THAT(comparison.relative_potential_difference, Catch::Matchers::WithinAbs(0.0, potential_parity_tolerance));

    CHECK_THAT(relative_binding_difference, Catch::Matchers::WithinAbs(0.0, energy_parity_tolerance));

    CHECK_THAT(relative_virial_difference, Catch::Matchers::WithinAbs(0.0, energy_parity_tolerance));

    CHECK_THAT(new_binding_analytic_error, Catch::Matchers::WithinAbs(0.0, analytic_energy_tolerance));

    CHECK_THAT(new_virial_analytic_error, Catch::Matchers::WithinAbs(0.0, analytic_energy_tolerance));

    CHECK_THAT(new_virial_consistency_error, Catch::Matchers::WithinAbs(0.0, consistency_tolerance));
}

TEST_CASE(
    "New Gravity Potential Deformed Rational Density Virial Self-Consistency",
    tags::gravity_consistency
) {
    auto args                = test_utils::setup_args();
    args.p.rtol              = 1.0e-13;
    args.p.max_iters         = std::max(args.p.max_iters, 1000);

    fem::FEM f               = fem::setup_fem(args.mesh_file, args, 0);

    const double radius      = utils::RADIUS;
    const double mass        = utils::MASS;

    constexpr double x_scale = 1.15;
    constexpr double y_scale = 0.95;
    constexpr double z_scale = 1.0 / (x_scale * y_scale);

    REQUIRE_THAT(x_scale * y_scale * z_scale, Catch::Matchers::WithinAbs(1.0, 1.0e-14));

    const double semi_axis_x = x_scale * radius;
    const double semi_axis_y = y_scale * radius;
    const double semi_axis_z = z_scale * radius;

    auto affine_displacement = [](const mfem::Vector &position, mfem::Vector &value) {
        value.SetSize(3);
        value(0) = (x_scale - 1.0) * position(0);
        value(1) = (y_scale - 1.0) * position(1);
        value(2) = (z_scale - 1.0) * position(2);
    };

    mfem::VectorFunctionCoefficient displacement_coefficient(3, affine_displacement);

    mfem::ParGridFunction displacement(f.displacementFes.get());

    displacement.ProjectCoefficient(displacement_coefficient);

    f.mapping->SetDisplacement(displacement);
    physics::update_stiffness_matrix(f);

    constexpr double concentration = 16.0;
    const double density_scale     = mass / std::pow(radius, 3.0);

    auto ellipsoidal_density = [semi_axis_x, semi_axis_y, semi_axis_z, density_scale](const mfem::Vector &position) {
        const double ellipsoidal_radius_squared = position(0) * position(0) / (semi_axis_x * semi_axis_x) +
                                                  position(1) * position(1) / (semi_axis_y * semi_axis_y) +
                                                  position(2) * position(2) / (semi_axis_z * semi_axis_z);

        if (ellipsoidal_radius_squared >= 1.0) {
            return 0.0;
        }

        const double denominator = 1.0 + concentration * ellipsoidal_radius_squared;

        return density_scale * (1.0 - ellipsoidal_radius_squared) / (denominator * denominator);
    };

    mapping::PhysicalPositionFunctionCoefficient density_coefficient(*f.mapping, ellipsoidal_density);

    mfem::GridFunction density(f.densityFes.get());

    density.ProjectCoefficient(density_coefficient);

    zero_vacuum_density(f, density);
    analysis::conserve_mass(f, density, mass);

    f.com                              = analysis::get_com(f, density);
    f.Q                                = physics::compute_quadrupole_moment_tensor(f, density, f.com);

    const double normalized_quadrupole = f.Q.FNorm() / (mass * radius * radius);

    INFO("Normalized quadrupole = " << normalized_quadrupole);

    REQUIRE(normalized_quadrupole > 1.0e-3);

    const physics::GravitySolution gravity_solution = physics::grav_potential_new(f, args, density, displacement);

    const int base_quadrature_order                 = get_gravity_quadrature_order(f);

    const GravitationalEnergies base_energies =
        compute_gravitational_energies(f, density, gravity_solution, base_quadrature_order);

    const GravitationalEnergies medium_energies =
        compute_gravitational_energies(f, density, gravity_solution, base_quadrature_order + 4);

    const GravitationalEnergies fine_energies =
        compute_gravitational_energies(f, density, gravity_solution, base_quadrature_order + 8);

    REQUIRE(fine_energies.binding < 0.0);
    REQUIRE(fine_energies.virial < 0.0);

    const double base_consistency_error =
        std::abs(base_energies.binding - base_energies.virial) / std::abs(base_energies.binding);

    const double medium_consistency_error =
        std::abs(medium_energies.binding - medium_energies.virial) / std::abs(medium_energies.binding);

    const double fine_consistency_error =
        std::abs(fine_energies.binding - fine_energies.virial) / std::abs(fine_energies.binding);

    const double binding_quadrature_change =
        std::abs(fine_energies.binding - medium_energies.binding) / std::abs(fine_energies.binding);

    const double virial_quadrature_change =
        std::abs(fine_energies.virial - medium_energies.virial) / std::abs(fine_energies.virial);

    INFO("Base-order consistency error = " << base_consistency_error);

    INFO("Medium-order consistency error = " << medium_consistency_error);

    INFO("Fine-order consistency error = " << fine_consistency_error);

    INFO("Medium-to-fine binding-energy change = " << binding_quadrature_change);

    INFO("Medium-to-fine virial-energy change = " << virial_quadrature_change);

    constexpr double virial_tolerance                = 1.0e-5;
    constexpr double diagnostic_quadrature_tolerance = 1.0e-7;

    CHECK_THAT(fine_consistency_error, Catch::Matchers::WithinAbs(0.0, virial_tolerance));

    CHECK_THAT(binding_quadrature_change, Catch::Matchers::WithinAbs(0.0, diagnostic_quadrature_tolerance));

    CHECK_THAT(virial_quadrature_change, Catch::Matchers::WithinAbs(0.0, diagnostic_quadrature_tolerance));
}

TEST_CASE(
    "New Gravity Potential Resolves Exterior Monopole By Compactification "
    "Shell",
    tags::gravity_analytic
) {
    auto args        = test_utils::setup_args();
    args.p.rtol      = 1.0e-13;
    args.p.max_iters = std::max(args.p.max_iters, 1000);

    fem::FEM f       = fem::setup_fem(args.mesh_file, args, 0);

    REQUIRE(f.mapping != nullptr);
    REQUIRE(f.domainMapperStateless != nullptr);
    REQUIRE(f.compactificationCoordinate != nullptr);

    const double radius          = utils::RADIUS;
    const double mass            = utils::MASS;
    const double analytic_volume = (4.0 / 3.0) * M_PI * radius * radius * radius;
    const double density         = mass / analytic_volume;

    mfem::ParGridFunction displacement(f.displacementFes.get());
    displacement = 0.0;

    f.mapping->ResetDisplacement();
    physics::update_stiffness_matrix(f);

    mfem::GridFunction rho_uniform(f.densityFes.get());
    rho_uniform = density;

    zero_vacuum_density(f, rho_uniform);
    analysis::conserve_mass(f, rho_uniform, mass);

    f.com                                          = analysis::get_com(f, rho_uniform);
    f.Q                                            = physics::compute_quadrupole_moment_tensor(f, rho_uniform, f.com);

    const physics::GravitySolution legacy_solution = physics::grav_potential(f, args, rho_uniform);

    const physics::GravitySolution new_solution    = physics::grav_potential_new(f, args, rho_uniform, displacement);

    const std::array<ExteriorMonopoleShellMetrics, 5> legacy_metrics =
        measure_exterior_monopole_shells(f, legacy_solution, displacement, ExteriorMonopoleMapping::legacy, mass);

    const std::array<ExteriorMonopoleShellMetrics, 5> new_metrics =
        measure_exterior_monopole_shells(f, new_solution, displacement, ExteriorMonopoleMapping::stateless, mass);

    double maximum_new_potential_error    = 0.0;
    double maximum_new_radial_field_error = 0.0;
    double maximum_new_tangential_field   = 0.0;

    for (int shell = 0; shell < 5; ++shell) {
        DYNAMIC_SECTION(
            "Exterior coordinate in [" << exterior_shell_boundaries[shell] << ", "
                                       << exterior_shell_boundaries[shell + 1] << ")"
        ) {
            const ExteriorMonopoleShellMetrics &legacy  = legacy_metrics[shell];

            const ExteriorMonopoleShellMetrics &current = new_metrics[shell];

            maximum_new_potential_error    = std::max(maximum_new_potential_error, current.potential_rms_error);

            maximum_new_radial_field_error = std::max(maximum_new_radial_field_error, current.radial_field_rms_error);

            maximum_new_tangential_field   = std::max(maximum_new_tangential_field, current.tangential_field_rms);

            INFO("Shell = " << shell);
            INFO(
                "Exterior-coordinate interval = [" << exterior_shell_boundaries[shell] << ", "
                                                   << exterior_shell_boundaries[shell + 1] << ")"
            );

            INFO("Legacy physical-radius range = [" << legacy.minimum_radius << ", " << legacy.maximum_radius << "]");

            INFO("New physical-radius range = [" << current.minimum_radius << ", " << current.maximum_radius << "]");

            INFO("Legacy scaled-potential RMS error = " << legacy.potential_rms_error);

            INFO("New scaled-potential RMS error = " << current.potential_rms_error);

            INFO("Legacy scaled-radial-field RMS error = " << legacy.radial_field_rms_error);

            INFO("New scaled-radial-field RMS error = " << current.radial_field_rms_error);

            INFO("Legacy scaled-tangential-field RMS = " << legacy.tangential_field_rms);

            INFO("New scaled-tangential-field RMS = " << current.tangential_field_rms);

            REQUIRE(legacy.quadrature_points > 0);
            REQUIRE(current.quadrature_points > 0);

            REQUIRE(std::isfinite(current.minimum_radius));
            REQUIRE(std::isfinite(current.maximum_radius));
            REQUIRE(std::isfinite(current.potential_rms_error));
            REQUIRE(std::isfinite(current.radial_field_rms_error));
            REQUIRE(std::isfinite(current.tangential_field_rms));

            CHECK(current.minimum_radius > 0.0);
            CHECK(current.maximum_radius > current.minimum_radius);

            /*
             * These are broad regression bounds, not the virial target.
             * The measured values will determine whether a monopole lift
             * is warranted.
             */
            CHECK(current.potential_rms_error < 1.0e-2);
            CHECK(current.radial_field_rms_error < 1.0e-2);
            CHECK(current.tangential_field_rms < 1.0e-2);
        }
    }

    for (int shell = 1; shell < 5; ++shell) {
        CHECK(new_metrics[shell].minimum_radius >= new_metrics[shell - 1].minimum_radius);

        CHECK(new_metrics[shell].maximum_radius > new_metrics[shell - 1].maximum_radius);
    }

    INFO("Maximum new scaled-potential shell error = " << maximum_new_potential_error);

    INFO("Maximum new scaled-radial-field shell error = " << maximum_new_radial_field_error);

    INFO("Maximum new scaled-tangential-field shell amplitude = " << maximum_new_tangential_field);
}

TEST_CASE(
    "New Exterior Monopole Error Is Separated From Finite Element Projection "
    "Floor",
    tags::gravity_analytic_accuracy
) {
    auto args        = test_utils::setup_args();
    args.p.rtol      = 1.0e-13;
    args.p.max_iters = std::max(args.p.max_iters, 1000);

    fem::FEM f       = fem::setup_fem(args.mesh_file, args, 0);

    REQUIRE(f.mapping != nullptr);
    REQUIRE(f.domainMapperStateless != nullptr);
    REQUIRE(f.compactificationCoordinate != nullptr);

    const double stellar_radius  = utils::RADIUS;

    const double mass            = utils::MASS;

    const double analytic_volume = (4.0 / 3.0) * M_PI * stellar_radius * stellar_radius * stellar_radius;

    const double density         = mass / analytic_volume;

    mfem::ParGridFunction displacement(f.displacementFes.get());
    displacement = 0.0;

    f.mapping->ResetDisplacement();
    physics::update_stiffness_matrix(f);

    mfem::GridFunction rho_uniform(f.densityFes.get());
    rho_uniform = density;

    zero_vacuum_density(f, rho_uniform);

    analysis::conserve_mass(f, rho_uniform, mass);

    f.com = analysis::get_com(f, rho_uniform);

    f.Q   = physics::compute_quadrupole_moment_tensor(f, rho_uniform, f.com);

    const physics::GravitySolution numerical_solution = physics::grav_potential_new(f, args, rho_uniform, displacement);

    StatelessMonopolePotentialCoefficient analytic_potential_coefficient(
        f, *f.domainMapperStateless, displacement, mass, stellar_radius
    );

    StatelessMonopoleHDivCoefficient analytic_field_coefficient(
        f, *f.domainMapperStateless, displacement, mass, stellar_radius
    );

    physics::GravitySolution analytic_projection(f);
    analytic_projection.phi     = 0.0;
    analytic_projection.gradPhi = 0.0;

    analytic_projection.phi.ProjectCoefficient(analytic_potential_coefficient);

    analytic_projection.gradPhi.ProjectCoefficient(analytic_field_coefficient);

    const std::array<ExteriorMonopoleShellMetrics, 5> numerical_metrics =
        measure_exterior_monopole_shells(f, numerical_solution, displacement, ExteriorMonopoleMapping::stateless, mass);

    const std::array<ExteriorMonopoleShellMetrics, 5> projection_metrics = measure_exterior_monopole_shells(
        f, analytic_projection, displacement, ExteriorMonopoleMapping::stateless, mass
    );

    mfem::Vector numerical_gradient_true;
    mfem::Vector numerical_potential_true;
    mfem::Vector projected_gradient_true;
    mfem::Vector projected_potential_true;

    numerical_solution.gradPhi.GetTrueDofs(numerical_gradient_true);

    numerical_solution.phi.GetTrueDofs(numerical_potential_true);

    analytic_projection.gradPhi.GetTrueDofs(projected_gradient_true);

    analytic_projection.phi.GetTrueDofs(projected_potential_true);

    MPI_Comm communicator = f.gravityFluxFes->GetComm();

    const double gradient_projection_gap =
        global_relative_vector_error(numerical_gradient_true, projected_gradient_true, communicator);

    const double potential_projection_gap =
        global_relative_vector_error(numerical_potential_true, projected_potential_true, communicator);

    double maximum_numerical_potential_error  = 0.0;
    double maximum_projected_potential_error  = 0.0;
    double maximum_numerical_radial_error     = 0.0;
    double maximum_projected_radial_error     = 0.0;
    double maximum_numerical_tangential_field = 0.0;
    double maximum_projected_tangential_field = 0.0;

    std::ostringstream report;

    report << "Global numerical/projection gradient DOF gap = " << gradient_projection_gap << '\n'
           << "Global numerical/projection potential DOF gap = " << potential_projection_gap << '\n';

    for (int shell = 0; shell < 5; ++shell) {
        const ExteriorMonopoleShellMetrics &numerical = numerical_metrics[shell];

        const ExteriorMonopoleShellMetrics &projected = projection_metrics[shell];

        maximum_numerical_potential_error = std::max(maximum_numerical_potential_error, numerical.potential_rms_error);

        maximum_projected_potential_error = std::max(maximum_projected_potential_error, projected.potential_rms_error);

        maximum_numerical_radial_error    = std::max(maximum_numerical_radial_error, numerical.radial_field_rms_error);

        maximum_projected_radial_error    = std::max(maximum_projected_radial_error, projected.radial_field_rms_error);

        maximum_numerical_tangential_field =
            std::max(maximum_numerical_tangential_field, numerical.tangential_field_rms);

        maximum_projected_tangential_field =
            std::max(maximum_projected_tangential_field, projected.tangential_field_rms);

        report << "Shell " << shell << " xi=[" << exterior_shell_boundaries[shell] << ", "
               << exterior_shell_boundaries[shell + 1] << "):\n"
               << "  radius range = [" << numerical.minimum_radius << ", " << numerical.maximum_radius << "]\n"
               << "  numerical potential error = " << numerical.potential_rms_error << '\n'
               << "  projected potential error = " << projected.potential_rms_error << '\n'
               << "  numerical radial-field error = " << numerical.radial_field_rms_error << '\n'
               << "  projected radial-field error = " << projected.radial_field_rms_error << '\n'
               << "  numerical tangential field = " << numerical.tangential_field_rms << '\n'
               << "  projected tangential field = " << projected.tangential_field_rms << '\n';
    }

    INFO(report.str());

    REQUIRE(std::isfinite(gradient_projection_gap));
    REQUIRE(std::isfinite(potential_projection_gap));

    REQUIRE(maximum_numerical_potential_error > 0.0);
    REQUIRE(maximum_projected_potential_error > 0.0);
    REQUIRE(maximum_numerical_radial_error > 0.0);
    REQUIRE(maximum_projected_radial_error > 0.0);

    /*
     * Broad guards against a broken projection. These are not the final
     * physical acceptance thresholds.
     */
    CHECK(maximum_projected_potential_error < 5.0e-2);
    CHECK(maximum_projected_radial_error < 5.0e-3);
    CHECK(maximum_projected_tangential_field < 5.0e-3);

    CHECK(maximum_numerical_potential_error < 5.0e-2);
    CHECK(maximum_numerical_radial_error < 5.0e-3);
    CHECK(maximum_numerical_tangential_field < 5.0e-3);

    /*
     * These are the decisive comparisons. If either fails, the solved
     * field is farther from the direct FE representation than it is from
     * the continuum monopole, indicating an operator-consistency issue
     * rather than a simple approximation floor.
     */
    CHECK(gradient_projection_gap < maximum_numerical_radial_error);

    CHECK(potential_projection_gap < maximum_numerical_potential_error);
}

TEST_CASE(
    "New Gravity Potential Matches Analytic Interior Potential For A Deformed "
    "Homogeneous Star",
    tags::gravity_analytic_accuracy
) {
    auto args        = test_utils::setup_args();
    args.p.rtol      = 1.0e-13;
    args.p.atol      = 1.0e-14;
    args.p.max_iters = std::max(args.p.max_iters, 2000);

    fem::FEM fem     = fem::setup_fem(args.mesh_file, args, 0);

    REQUIRE(fem.mapping != nullptr);
    REQUIRE(fem.domainMapperStateless != nullptr);

    const double radius      = utils::RADIUS;
    const double mass        = utils::MASS;

    constexpr double x_scale = 1.15;
    constexpr double y_scale = 0.95;
    constexpr double z_scale = 1.0 / (x_scale * y_scale);

    const double semi_axis_x = x_scale * radius;
    const double semi_axis_y = y_scale * radius;
    const double semi_axis_z = z_scale * radius;

    REQUIRE_THAT(x_scale * y_scale * z_scale, Catch::Matchers::WithinAbs(1.0, 1.0e-14));

    auto displacement_function = [](const mfem::Vector &position, mfem::Vector &value) {
        value.SetSize(3);
        value(0) = (x_scale - 1.0) * position(0);
        value(1) = (y_scale - 1.0) * position(1);
        value(2) = (z_scale - 1.0) * position(2);
    };

    mfem::VectorFunctionCoefficient displacement_coefficient(3, displacement_function);
    mfem::ParGridFunction displacement(fem.displacementFes.get());
    displacement.ProjectCoefficient(displacement_coefficient);

    fem.mapping->SetDisplacement(displacement);
    physics::update_stiffness_matrix(fem);

    const double analytic_volume = (4.0 / 3.0) * M_PI * semi_axis_x * semi_axis_y * semi_axis_z;
    const double density_value   = mass / analytic_volume;

    mfem::GridFunction density(fem.densityFes.get());
    density = density_value;
    zero_vacuum_density(fem, density);

    const double projected_mass    = analysis::domain_integrate_grid_function(fem, density, utils::DOMAINS::STELLAR);
    const double numerical_density = density_value * mass / projected_mass;
    analysis::conserve_mass(fem, density, mass);

    fem.com = analysis::get_com(fem, density);
    fem.Q   = physics::compute_quadrupole_moment_tensor(fem, density, fem.com);

    const HomogeneousEllipsoidAnalytic analytic =
        compute_homogeneous_ellipsoid_analytic(semi_axis_x, semi_axis_y, semi_axis_z);

    auto analytic_potential = [numerical_density, semi_axis_x, semi_axis_y, semi_axis_z,
                               analytic](const mfem::Vector &position) {
        const double potential_kernel = semi_axis_x * semi_axis_y * semi_axis_z * analytic.energy_kernel -
                                        analytic.coefficient_x * position(0) * position(0) -
                                        analytic.coefficient_y * position(1) * position(1) -
                                        analytic.coefficient_z * position(2) * position(2);

        return -M_PI * utils::G * numerical_density * potential_kernel;
    };

    mapping::PhysicalPositionFunctionCoefficient analytic_potential_coefficient(*fem.mapping, analytic_potential);
    const int vacuum_attribute = fem.domainMapperStateless->GetVacuumElementAttribute();
    mfem::ParGridFunction projected_potential(fem.gravityPotentialFes.get());
    projected_potential.ProjectCoefficient(analytic_potential_coefficient);
    const physics::GravitySolution solution      = physics::grav_potential_new(fem, args, density, displacement);

    const int quadrature_order                   = get_gravity_quadrature_order(fem);
    double local_solution_error_squared          = 0.0;
    double local_projection_error_squared        = 0.0;
    double local_solution_projection_gap_squared = 0.0;
    double local_analytic_norm_squared           = 0.0;
    double local_projected_norm_squared          = 0.0;
    double local_maximum_relative_error          = 0.0;

    mfem::Vector physical_position(3);
    mfem::DenseMatrix mapping_jacobian(3);

    for (int element_id = 0; element_id < fem.mesh->GetNE(); ++element_id) {
        mfem::ElementTransformation *transformation = fem.mesh->GetElementTransformation(element_id);
        if (transformation->Attribute == vacuum_attribute) {
            continue;
        }

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

        for (int quadrature_point_id = 0; quadrature_point_id < rule.GetNPoints(); ++quadrature_point_id) {
            const mfem::IntegrationPoint &point = rule.IntPoint(quadrature_point_id);
            transformation->SetIntPoint(&point);

            fem.mapping->GetPhysicalPoint(*transformation, point, physical_position);
            fem.mapping->ComputeJacobian(*transformation, mapping_jacobian);
            const double mapping_determinant = mapping_jacobian.Det();
            MFEM_VERIFY(
                mapping_determinant > 0.0, "Deformed potential test encountered a "
                                           "non-positive mapping determinant."
            );

            const double expected_potential        = analytic_potential(physical_position);
            const double computed_potential        = solution.phi.GetValue(element_id, point);
            const double projected_potential_value = projected_potential.GetValue(element_id, point);
            const double weight                    = point.weight * transformation->Weight() * mapping_determinant;

            local_solution_error_squared +=
                weight * (computed_potential - expected_potential) * (computed_potential - expected_potential);
            local_projection_error_squared += weight * (projected_potential_value - expected_potential) *
                                              (projected_potential_value - expected_potential);
            local_solution_projection_gap_squared += weight * (computed_potential - projected_potential_value) *
                                                     (computed_potential - projected_potential_value);
            local_analytic_norm_squared += weight * expected_potential * expected_potential;
            local_projected_norm_squared += weight * projected_potential_value * projected_potential_value;
            local_maximum_relative_error = std::max(
                local_maximum_relative_error,
                std::abs(computed_potential - expected_potential) /
                    std::max(std::abs(expected_potential), std::numeric_limits<double>::epsilon())
            );
        }
    }

    const std::array<double, 5> local_values{
        local_solution_error_squared, local_projection_error_squared, local_solution_projection_gap_squared,
        local_analytic_norm_squared, local_projected_norm_squared
    };
    std::array<double, 5> global_values{};
    MPI_Allreduce(
        local_values.data(), global_values.data(), static_cast<int>(local_values.size()), MPI_DOUBLE, MPI_SUM,
        fem.densityFes->GetComm()
    );

    double maximum_relative_error = 0.0;
    MPI_Allreduce(
        &local_maximum_relative_error, &maximum_relative_error, 1, MPI_DOUBLE, MPI_MAX, fem.densityFes->GetComm()
    );

    const double solution_relative_error   = std::sqrt(global_values[0] / global_values[3]);
    const double projection_relative_error = std::sqrt(global_values[1] / global_values[3]);
    const double solution_projection_gap   = std::sqrt(global_values[2] / global_values[4]);

    INFO(
        "Ellipsoid coefficients = (" << analytic.coefficient_x << ", " << analytic.coefficient_y << ", "
                                     << analytic.coefficient_z << ")"
    );
    INFO("Analytic interior-potential L2 relative error = " << solution_relative_error);
    INFO("Analytic-potential FE projection L2 relative error = " << projection_relative_error);
    INFO("New-solver / analytic-potential projection relative gap = " << solution_projection_gap);
    INFO("Maximum interior pointwise relative potential error = " << maximum_relative_error);

    REQUIRE(std::isfinite(solution_relative_error));
    REQUIRE(std::isfinite(projection_relative_error));
    REQUIRE(std::isfinite(solution_projection_gap));
    REQUIRE(std::isfinite(maximum_relative_error));

    /*
     * The field test establishes an O(1e-3) representation floor on this
     * mesh. These are intentionally accuracy-regression guards, not claims
     * of analytic convergence of the unrefined RT/L2 representation.
     */
    CHECK(solution_relative_error < 1.0e-5);
    CHECK(projection_relative_error < 1.0e-5);
    CHECK(solution_projection_gap < 1.0e-5);
    CHECK(maximum_relative_error < 1.0e-5);
}

struct FerrersN1Analytic {
    double potential_constant;
    std::array<double, 3> first_coefficients;
    std::array<std::array<double, 3>, 3> second_coefficients;
};

class FerrersVacuumMaskedCoefficient final : public mfem::Coefficient {
public:
    FerrersVacuumMaskedCoefficient(
        Coefficient &coefficient,
        const int vacuum_attribute
    )
        : m_coefficient(coefficient),
          m_vacuum_attribute(vacuum_attribute) {
    }

    double Eval(
        mfem::ElementTransformation &transformation,
        const mfem::IntegrationPoint &integration_point
    ) override {
        if (transformation.Attribute == m_vacuum_attribute) {
            return 0.0;
        }

        return m_coefficient.Eval(transformation, integration_point);
    }

private:
    Coefficient &m_coefficient;
    int m_vacuum_attribute;
};

static double compute_ferrers_n1_coefficient(
    const std::array<
        double,
        3> &semi_axes,
    const int first_denominator_axis,
    const int second_denominator_axis
) {
    const double axis_product = semi_axes[0] * semi_axes[1] * semi_axes[2];

    const double length_scale = std::cbrt(axis_product);

    auto integrand            = [semi_axes, axis_product, length_scale, first_denominator_axis,
                                 second_denominator_axis](const double t) {
        if (t <= 0.0 || t >= 1.0) {
            return 0.0;
        }

        /*
         * Map u in [0, infinity) to t in [0, 1]:
         *
         *     u = L^2 [t / (1 - t)]^2.
         */
        const double one_minus_t = 1.0 - t;
        const double s           = t / one_minus_t;

        const double u           = length_scale * length_scale * s * s;

        const double du_dt       = length_scale * length_scale * 2.0 * s / (one_minus_t * one_minus_t);

        const double delta       = std::sqrt(
            (semi_axes[0] * semi_axes[0] + u) * (semi_axes[1] * semi_axes[1] + u) * (semi_axes[2] * semi_axes[2] + u)
        );

        double value = axis_product * du_dt / delta;

        if (first_denominator_axis >= 0) {
            value /= semi_axes[first_denominator_axis] * semi_axes[first_denominator_axis] + u;
        }

        if (second_denominator_axis >= 0) {
            value /= semi_axes[second_denominator_axis] * semi_axes[second_denominator_axis] + u;
        }

        return value;
    };

    double integration_error = 0.0;

    return boost::math::quadrature::gauss_kronrod<double, 61>::integrate(
        integrand, 0.0, 1.0, 15, 1.0e-13, &integration_error
    );
}

static FerrersN1Analytic compute_ferrers_n1_analytic(
    const double semi_axis_x,
    const double semi_axis_y,
    const double semi_axis_z
) {
    const std::array<double, 3> semi_axes{semi_axis_x, semi_axis_y, semi_axis_z};

    FerrersN1Analytic analytic{
        .potential_constant  = compute_ferrers_n1_coefficient(semi_axes, -1, -1),
        .first_coefficients  = {},
        .second_coefficients = {}
    };

    for (int axis = 0; axis < 3; ++axis) {
        analytic.first_coefficients[axis] = compute_ferrers_n1_coefficient(semi_axes, axis, -1);
    }

    for (int first_axis = 0; first_axis < 3; ++first_axis) {
        for (int second_axis = first_axis; second_axis < 3; ++second_axis) {
            const double coefficient = compute_ferrers_n1_coefficient(semi_axes, first_axis, second_axis);

            analytic.second_coefficients[first_axis][second_axis] = coefficient;

            analytic.second_coefficients[second_axis][first_axis] = coefficient;
        }
    }

    return analytic;
}

static double evaluate_ferrers_n1_potential(
    const mfem::Vector &position,
    const double central_density,
    const FerrersN1Analytic &analytic
) {
    const std::array<double, 3> coordinate_squared{
        position(0) * position(0), position(1) * position(1), position(2) * position(2)
    };

    /*
     * Expansion of
     *
     *  -pi G rho_c abc / 2
     *      integral [(1 - m^2(u))^2 / Delta(u)] du.
     */
    double potential_kernel = analytic.potential_constant;

    for (int first_axis = 0; first_axis < 3; ++first_axis) {
        potential_kernel -= 2.0 * analytic.first_coefficients[first_axis] * coordinate_squared[first_axis];

        for (int second_axis = 0; second_axis < 3; ++second_axis) {
            potential_kernel += analytic.second_coefficients[first_axis][second_axis] * coordinate_squared[first_axis] *
                                coordinate_squared[second_axis];
        }
    }

    return -0.5 * M_PI * utils::G * central_density * potential_kernel;
}

static void evaluate_ferrers_n1_gradient(
    const mfem::Vector &position,
    const double central_density,
    const FerrersN1Analytic &analytic,
    mfem::Vector &gradient
) {
    gradient.SetSize(3);

    const std::array<double, 3> coordinate_squared{
        position(0) * position(0), position(1) * position(1), position(2) * position(2)
    };

    for (int axis = 0; axis < 3; ++axis) {
        double coefficient = analytic.first_coefficients[axis];

        for (int other_axis = 0; other_axis < 3; ++other_axis) {
            coefficient -= analytic.second_coefficients[axis][other_axis] * coordinate_squared[other_axis];
        }

        /*
         * The solver stores grad(Phi), which points outward for a
         * negative gravitational potential.
         */
        gradient(axis) = 2.0 * M_PI * utils::G * central_density * position(axis) * coefficient;
    }
}

TEST_CASE(
    "New Gravity Potential Matches Analytic Ferrers Ellipsoid",
    tags::gravity_analytic_accuracy
) {
    auto args        = test_utils::setup_args();

    args.p.rtol      = 1.0e-13;
    args.p.atol      = 1.0e-14;
    args.p.max_iters = std::max(args.p.max_iters, 2000);

    fem::FEM fem     = fem::setup_fem(args.mesh_file, args, 0);

    REQUIRE(fem.mapping != nullptr);
    REQUIRE(fem.domainMapperStateless != nullptr);

    const double radius      = utils::RADIUS;

    const double mass        = utils::MASS;

    constexpr double x_scale = 1.0;
    constexpr double y_scale = 1.0;
    constexpr double z_scale = 1.0 / (x_scale * y_scale);

    const double semi_axis_x = x_scale * radius;
    const double semi_axis_y = y_scale * radius;
    const double semi_axis_z = z_scale * radius;

    REQUIRE_THAT(x_scale * y_scale * z_scale, Catch::Matchers::WithinAbs(1.0, 1.0e-14));

    auto displacement_function = [](const mfem::Vector &position, mfem::Vector &value) {
        value.SetSize(3);

        value(0) = (x_scale - 1.0) * position(0);
        value(1) = (y_scale - 1.0) * position(1);
        value(2) = (z_scale - 1.0) * position(2);
    };

    mfem::VectorFunctionCoefficient displacement_coefficient(3, displacement_function);

    mfem::ParGridFunction displacement(fem.displacementFes.get());

    displacement.ProjectCoefficient(displacement_coefficient);

    fem.mapping->SetDisplacement(displacement);

    physics::update_stiffness_matrix(fem);

    const int vacuum_attribute   = fem.domainMapperStateless->GetVacuumElementAttribute();

    /*
     * For rho = rho_c (1 - m^2), the exact mass is
     *
     *     M = 8 pi a b c rho_c / 15.
     */
    const double central_density = 15.0 * mass / (8.0 * M_PI * semi_axis_x * semi_axis_y * semi_axis_z);

    auto density_function = [central_density, semi_axis_x, semi_axis_y, semi_axis_z](const mfem::Vector &position) {
        const double ellipsoidal_radius_squared = position(0) * position(0) / (semi_axis_x * semi_axis_x) +
                                                  position(1) * position(1) / (semi_axis_y * semi_axis_y) +
                                                  position(2) * position(2) / (semi_axis_z * semi_axis_z);

        return central_density * std::max(0.0, 1.0 - ellipsoidal_radius_squared);
    };

    mapping::PhysicalPositionFunctionCoefficient physical_density_coefficient(*fem.mapping, density_function);

    FerrersVacuumMaskedCoefficient stellar_density_coefficient(physical_density_coefficient, vacuum_attribute);

    mfem::GridFunction density(fem.densityFes.get());

    density = 0.0;

    density.ProjectCoefficient(stellar_density_coefficient);

    const double projected_mass = analysis::domain_integrate_grid_function(fem, density, utils::DOMAINS::STELLAR);

    REQUIRE(std::isfinite(projected_mass));
    REQUIRE(projected_mass > 0.0);

    /*
     * Keep the projected source at exactly the requested mass. Because
     * projection and scaling are linear, this also gives the central
     * density appropriate to the represented source.
     */
    const double density_scale = mass / projected_mass;

    density *= density_scale;

    const double represented_central_density = central_density * density_scale;

    fem.com                                  = analysis::get_com(fem, density);

    fem.Q                                    = physics::compute_quadrupole_moment_tensor(fem, density, fem.com);

    const double normalized_quadrupole       = fem.Q.FNorm() / (mass * radius * radius);

    // REQUIRE(normalized_quadrupole > 1.0e-3);

    const FerrersN1Analytic analytic         = compute_ferrers_n1_analytic(semi_axis_x, semi_axis_y, semi_axis_z);

    /*
     * Independent analytic consistency checks.
     *
     * Sum(A_i) = 2 supplies the constant part of Poisson's
     * equation. The B_ij identities supply the -m^2 part.
     */
    const double first_coefficient_sum =
        analytic.first_coefficients[0] + analytic.first_coefficients[1] + analytic.first_coefficients[2];

    REQUIRE_THAT(first_coefficient_sum, Catch::Matchers::WithinAbs(2.0, 1.0e-11));

    const std::array<double, 3> semi_axes_squared{
        semi_axis_x * semi_axis_x, semi_axis_y * semi_axis_y, semi_axis_z * semi_axis_z
    };

    for (int axis = 0; axis < 3; ++axis) {
        double poisson_coefficient = 3.0 * analytic.second_coefficients[axis][axis];

        for (int other_axis = 0; other_axis < 3; ++other_axis) {
            if (other_axis != axis) {
                poisson_coefficient += analytic.second_coefficients[axis][other_axis];
            }
        }

        REQUIRE_THAT(poisson_coefficient, Catch::Matchers::WithinRel(2.0 / semi_axes_squared[axis], 1.0e-10));
    }

    const physics::GravitySolution solution = physics::grav_potential_new(fem, args, density, displacement);

    auto analytic_potential_function        = [represented_central_density, analytic](const mfem::Vector &position) {
        return evaluate_ferrers_n1_potential(position, represented_central_density, analytic);
    };

    mapping::PhysicalPositionFunctionCoefficient physical_potential_coefficient(
        *fem.mapping, analytic_potential_function
    );

    /*
     * This wrapper is required because scalar ProjectCoefficient has no
     * attribute overload. It also prevents evaluation of the quartic
     * interior formula in the compactified vacuum.
     */
    FerrersVacuumMaskedCoefficient stellar_potential_coefficient(physical_potential_coefficient, vacuum_attribute);

    mfem::ParGridFunction projected_potential(fem.gravityPotentialFes.get());

    projected_potential = 0.0;

    projected_potential.ProjectCoefficient(stellar_potential_coefficient);

    const int quadrature_order                   = get_gravity_quadrature_order(fem) + 4;

    double local_potential_error_squared         = 0.0;
    double local_projection_error_squared        = 0.0;
    double local_solution_projection_gap_squared = 0.0;
    double local_potential_norm_squared          = 0.0;
    double local_projection_norm_squared         = 0.0;
    double local_field_error_squared             = 0.0;
    double local_field_norm_squared              = 0.0;
    double local_maximum_potential_error         = 0.0;

    mfem::Vector physical_position(3);
    mfem::Vector reference_field(3);
    mfem::Vector physical_field(3);
    mfem::Vector analytic_field(3);
    mfem::Vector field_difference(3);

    mfem::DenseMatrix mapping_jacobian(3, 3);

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

        if (transformation->Attribute == vacuum_attribute) {
            continue;
        }

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

        for (int quadrature_point_id = 0; quadrature_point_id < rule.GetNPoints(); ++quadrature_point_id) {
            const mfem::IntegrationPoint &point = rule.IntPoint(quadrature_point_id);

            transformation->SetIntPoint(&point);

            fem.mapping->GetPhysicalPoint(*transformation, point, physical_position);

            fem.mapping->ComputeJacobian(*transformation, mapping_jacobian);

            const double mapping_determinant = mapping_jacobian.Det();

            MFEM_VERIFY(
                std::isfinite(mapping_determinant) && mapping_determinant > 0.0,
                "Ferrers test encountered an invalid mapping determinant."
            );

            const double expected_potential =
                evaluate_ferrers_n1_potential(physical_position, represented_central_density, analytic);

            evaluate_ferrers_n1_gradient(physical_position, represented_central_density, analytic, analytic_field);

            const double computed_potential        = solution.phi.GetValue(element_id, point);

            const double projected_potential_value = projected_potential.GetValue(element_id, point);

            solution.gradPhi.GetVectorValue(element_id, point, reference_field);

            mapping_jacobian.Mult(reference_field, physical_field);

            physical_field /= mapping_determinant;

            field_difference = physical_field;

            field_difference -= analytic_field;

            const double weight                         = point.weight * transformation->Weight() * mapping_determinant;

            const double potential_error                = computed_potential - expected_potential;

            const double projection_error               = projected_potential_value - expected_potential;

            const double solution_projection_difference = computed_potential - projected_potential_value;

            local_potential_error_squared += weight * potential_error * potential_error;

            local_projection_error_squared += weight * projection_error * projection_error;

            local_solution_projection_gap_squared +=
                weight * solution_projection_difference * solution_projection_difference;

            local_potential_norm_squared += weight * expected_potential * expected_potential;

            local_projection_norm_squared += weight * projected_potential_value * projected_potential_value;

            local_field_error_squared += weight * (field_difference * field_difference);

            local_field_norm_squared += weight * (analytic_field * analytic_field);

            local_maximum_potential_error = std::max(
                local_maximum_potential_error,
                std::abs(potential_error) /
                    std::max(std::abs(expected_potential), std::numeric_limits<double>::epsilon())
            );
        }
    }

    const std::array<double, 7> local_values{
        local_potential_error_squared, local_projection_error_squared, local_solution_projection_gap_squared,
        local_potential_norm_squared,  local_projection_norm_squared,  local_field_error_squared,
        local_field_norm_squared
    };

    std::array<double, 7> global_values{};

    MPI_Allreduce(
        local_values.data(), global_values.data(), static_cast<int>(local_values.size()), MPI_DOUBLE, MPI_SUM,
        fem.densityFes->GetComm()
    );

    double maximum_potential_error = 0.0;

    MPI_Allreduce(
        &local_maximum_potential_error, &maximum_potential_error, 1, MPI_DOUBLE, MPI_MAX, fem.densityFes->GetComm()
    );

    REQUIRE(global_values[3] > 0.0);
    REQUIRE(global_values[4] > 0.0);
    REQUIRE(global_values[6] > 0.0);

    const double potential_relative_error  = std::sqrt(global_values[0] / global_values[3]);

    const double projection_relative_error = std::sqrt(global_values[1] / global_values[3]);

    const double solution_projection_gap   = std::sqrt(global_values[2] / global_values[4]);

    const double field_relative_error      = std::sqrt(global_values[5] / global_values[6]);

    INFO("Projected mass before normalization = " << projected_mass);

    INFO("Density normalization factor = " << density_scale);

    INFO("Normalized quadrupole = " << normalized_quadrupole);

    INFO("Ferrers potential L2 relative error = " << potential_relative_error);

    INFO("Ferrers potential FE-projection relative error = " << projection_relative_error);

    INFO("New-solver / Ferrers-potential projection gap = " << solution_projection_gap);

    INFO("Ferrers field L2 relative error = " << field_relative_error);

    INFO("Maximum interior pointwise potential relative error = " << maximum_potential_error);

    REQUIRE(std::isfinite(potential_relative_error));
    REQUIRE(std::isfinite(projection_relative_error));
    REQUIRE(std::isfinite(solution_projection_gap));
    REQUIRE(std::isfinite(field_relative_error));
    REQUIRE(std::isfinite(maximum_potential_error));

    /*
     * Initial characterization guards. Unlike the homogeneous case,
     * this exact potential is quartic, so the FE representation floor
     * will generally be higher. Record the values before deciding
     * whether tighter regression thresholds are appropriate.
     */
    CHECK(potential_relative_error < 1.0e-5);
    CHECK(projection_relative_error < 1.0e-5);
    CHECK(solution_projection_gap < 1.0e-5);
    CHECK(field_relative_error < 1.0e-5);
    CHECK(maximum_potential_error < 1.0e-5);
}