#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 {
// Concentration-16 rational profiles have a stable mixed-projection virial
// floor of approximately 1.53e-5 on the regression mesh. Increasing the
// diagnostic quadrature order changes each energy by only O(1e-13).
constexpr double rational_profile_virial_tolerance = 2.0e-5;

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_mapper,
      const mfem::GridFunction &displacement,
      const mfem::GridFunction &compactification_coordinate,
      const double density,
      const HomogeneousEllipsoidAnalytic &analytic)
      : VectorCoefficient(3),
        mapping_evaluator(domain_mapper, displacement,
                          compactification_coordinate),
        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 field_physical(3);
    mapping::MappingPointContext context;
    MFEM_VERIFY(mapping_evaluator.EvaluatePoint(
                    transformation, integration_point, context) ==
                    mapping::MappingStatus::valid,
                "Ellipsoid coefficient encountered an invalid mapping.");
    const mfem::Vector &x_physical = context.physical_position;

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

    mapping::MapPhysicalFluxToHDivReference(context, field_physical, value);
  }

private:
  mapping::GridFunctionMappingEvaluator mapping_evaluator;
  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);
  using DomainSchema = utils::domain::CoreEnvelopeVacuumDomainSchema;
  mapping::GridFunctionMappingEvaluator mapping_evaluator(
      *f.domainMapperStateless, *f.displacement,
      *f.compactificationCoordinate);

  for (int elem_id = 0; elem_id < f.mesh->GetNE(); ++elem_id) {
    if (!DomainSchema::template attribute_belongs_to<utils::domain::Stellar>(
            f.mesh->GetAttribute(elem_id))) {
      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);

      mapping::VolumeMappingContext mapping_context;
      MFEM_VERIFY(
          mapping_evaluator.EvaluateVolume(*transformation, integration_point,
                                           mapping_context) ==
              mapping::MappingStatus::valid,
          "Gravity-energy integration encountered an invalid mapping.");
      const double weight = mapping_context.quadrature.weight;
      x_physical = mapping_context.mapping.physical_position;
      gravity_solution.gradPhi.GetVectorValue(elem_id, integration_point,
                                              grad_phi_element);
      mapping::MapHDivFluxToPhysical(mapping_context.mapping,
                                     grad_phi_element, 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) {
  using DomainSchema = utils::domain::CoreEnvelopeVacuumDomainSchema;

  const field::FieldDofMap densityMap =
      field::make_field_dof_map<field::Density, DomainSchema>(*f.densityFes);

  mfem::Vector densityTrue;
  rho.GetTrueDofs(densityTrue);

  const mfem::Vector supportedDensity = densityMap.gather(densityTrue);
  densityMap.scatter(supportedDensity, densityTrue);
  rho.SetFromTrueDofs(densityTrue);
}

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 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 double mass) {
  REQUIRE(f.mesh != nullptr);
  REQUIRE(f.gravityFluxFes != nullptr);
  REQUIRE(f.displacementFes != nullptr);
  REQUIRE(f.compactificationFes != nullptr);
  REQUIRE(f.compactificationCoordinate != nullptr);
  REQUIRE(f.domainMapperStateless != 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::DomainMapper::Workspace workspace(f.mesh->Dimension());

  const int vacuum_attribute = field_dof_test_utils::vacuum_material_attribute;

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

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

      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::DomainMapper &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::DomainMapper &m_domain_mapper;
  const mfem::GridFunction &m_displacement;
  mapping::DomainMapper::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::DomainMapper &domain_mapper,
      const mfem::GridFunction &displacement, const double mass,
      const double stellar_radius)
      : m_geometry(f, domain_mapper, displacement),
        m_vacuum_attribute(field_dof_test_utils::vacuum_material_attribute),
        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::DomainMapper &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(field_dof_test_utils::vacuum_material_attribute),
        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) {
  auto args = test_utils::setup_args();
  fem::FEM f = fem::setup_fem(args.mesh_file, args, 0);
  *f.displacement = 0.0;

  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 auto gravity_solution =
      physics::solve_gravity_field(f, args, rho_uniform, displacement);
  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;
  mapping::GridFunctionMappingEvaluator mapping_evaluator(
      *f.domainMapperStateless, *f.displacement,
      *f.compactificationCoordinate);

  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;
    mapping_evaluator.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) {
  auto args = test_utils::setup_args();
  fem::FEM f = fem::setup_fem(args.mesh_file, args, 0);
  *f.displacement = 0.0;

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

  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.domainMapperStateless, *f.displacement,
        *f.compactificationCoordinate, 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::solve_gravity_field(f, args, rho_grid, displacement);
  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) {
  auto args = test_utils::setup_args();
  fem::FEM f = fem::setup_fem(args.mesh_file, args, 0);
  *f.displacement = 0.0;

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

  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.domainMapperStateless, *f.displacement,
        *f.compactificationCoordinate, 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::solve_gravity_field(f, args, rho_grid, displacement);
  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);

  CHECK_THAT(
      relative_consistency_error,
      Catch::Matchers::WithinAbs(0.0, rational_profile_virial_tolerance));
}

TEST_CASE("Homogeneous Ellipsoid Analytic Gravity", tags::gravity_analytic) {
  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.displacement = displacement;

  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.domainMapperStateless, *f.displacement,
      *f.compactificationCoordinate, 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::solve_gravity_field(f, args, rho_grid, displacement);
  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::Vector gravity_projection_difference(3);
  mapping::GridFunctionMappingEvaluator mapping_evaluator(
      *f.domainMapperStateless, *f.displacement,
      *f.compactificationCoordinate);

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

      mapping::VolumeMappingContext mapping_context;
      MFEM_VERIFY(
          mapping_evaluator.EvaluateVolume(*transformation, integration_point,
                                           mapping_context) ==
              mapping::MappingStatus::valid,
          "Ellipsoid field comparison encountered an invalid mapping.");
      const double map_determinant =
          mapping_context.mapping.mapping_determinant;
      const double weight = mapping_context.quadrature.weight;
      x_physical = mapping_context.mapping.physical_position;
      gravity_solution.gradPhi.GetVectorValue(elem_id, integration_point,
                                              grad_phi_element);
      const mfem::DenseMatrix &map_jacobian =
          mapping_context.mapping.mapping_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));

  // The RT projection itself is accurate to approximately 8.8e-4 on this
  // mesh. The solved field should remain close to that best representable
  // field, while the integrated energies have a substantially lower floor.
  constexpr double quadrupole_tolerance = 2.0e-4;
  constexpr double field_tolerance = 1.0e-3;
  constexpr double binding_energy_tolerance = 1.0e-5;
  constexpr double virial_energy_tolerance = 5.0e-5;
  constexpr double consistency_tolerance = 5.0e-5;
  constexpr double projection_gap_tolerance = 0.3;
  constexpr double projection_ratio_tolerance = 1.05;

  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, binding_energy_tolerance));
  CHECK_THAT(relative_virial_energy_error,
             Catch::Matchers::WithinAbs(0.0, virial_energy_tolerance));
  CHECK_THAT(relative_consistency_error,
             Catch::Matchers::WithinAbs(0.0, consistency_tolerance));
  CHECK(projection_gap_ratio < projection_gap_tolerance);
  CHECK(gravity_to_projection_error_ratio < projection_ratio_tolerance);
}

TEST_CASE("Deformed Rational Density Virial Self-Consistency",
          tags::gravity_consistency) {
  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.displacement = displacement;

  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.domainMapperStateless, *f.displacement,
        *f.compactificationCoordinate, 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::solve_gravity_field(f, args, rho_grid, displacement);
  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);

  CHECK_THAT(
      relative_consistency_error,
      Catch::Matchers::WithinAbs(0.0, rational_profile_virial_tolerance));
}

TEST_CASE("Gravity Field 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.displacement = 0.0;

  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::solve_gravity_field(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;
  mapping::GridFunctionMappingEvaluator mapping_evaluator(
      *f.domainMapperStateless, *f.displacement,
      *f.compactificationCoordinate);

  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;
    mapping_evaluator.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("Gravity Field 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.displacement = displacement;

  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.domainMapperStateless, *f.displacement,
      *f.compactificationCoordinate, 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::solve_gravity_field(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 diagnostic_quadrature_tolerance = 1.0e-7;

  CHECK_THAT(
      fine_consistency_error,
      Catch::Matchers::WithinAbs(0.0, rational_profile_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("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.domainMapperStateless != 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.displacement = 0.0;

  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::solve_gravity_field(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,
                                       mass);

  const std::array<ExteriorMonopoleShellMetrics, 5> projection_metrics =
      measure_exterior_monopole_shells(f, analytic_projection, displacement,
                                       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("Gravity Field 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.domainMapperStateless != 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.displacement = displacement;

  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.domainMapperStateless, *fem.displacement,
      *fem.compactificationCoordinate, analytic_potential);
  const int vacuum_attribute = field_dof_test_utils::vacuum_material_attribute;
  mfem::ParGridFunction projected_potential(fem.gravityPotentialFes.get());
  projected_potential.ProjectCoefficient(analytic_potential_coefficient);
  const physics::GravitySolution solution =
      physics::solve_gravity_field(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);
  mapping::GridFunctionMappingEvaluator mapping_evaluator(
      *fem.domainMapperStateless, *fem.displacement,
      *fem.compactificationCoordinate);

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

      mapping::MappingPointContext mapping_context;
      MFEM_VERIFY(
          mapping_evaluator.EvaluatePoint(*transformation, point,
                                          mapping_context) ==
              mapping::MappingStatus::valid,
          "Deformed potential test encountered an invalid mapping.");
      physical_position = mapping_context.physical_position;
      const double mapping_determinant =
          mapping_context.mapping_determinant;
      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));

  /*
   * On the regression mesh, the direct L2 projection floor is about 6.7e-6,
   * the mixed-solve/projection gap is about 3.3e-5, and the maximum
   * pointwise error is about 2e-4. These bounds guard those independently.
   */
  CHECK(solution_relative_error < 5.0e-5);
  CHECK(projection_relative_error < 1.0e-5);
  CHECK(solution_projection_gap < 5.0e-5);
  CHECK(maximum_relative_error < 2.5e-4);
}

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("Gravity Field 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.domainMapperStateless != 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.displacement = displacement;

  const int vacuum_attribute = field_dof_test_utils::vacuum_material_attribute;

  /*
   * 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.domainMapperStateless, *fem.displacement,
      *fem.compactificationCoordinate, 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::solve_gravity_field(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.domainMapperStateless, *fem.displacement,
      *fem.compactificationCoordinate, 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);

  mapping::GridFunctionMappingEvaluator mapping_evaluator(
      *fem.domainMapperStateless, *fem.displacement,
      *fem.compactificationCoordinate);

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

      mapping::MappingPointContext mapping_context;
      MFEM_VERIFY(
          mapping_evaluator.EvaluatePoint(*transformation, point,
                                          mapping_context) ==
              mapping::MappingStatus::valid,
          "Ferrers test encountered an invalid mapping.");
      physical_position = mapping_context.physical_position;
      const mfem::DenseMatrix &mapping_jacobian =
          mapping_context.mapping_jacobian;
      const double mapping_determinant =
          mapping_context.mapping_determinant;

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

  /*
   * On the regression mesh, the quartic analytic potential and its RT
   * gradient have representation floors of about 1.1e-4 and 8.7e-4,
   * respectively. The mixed solution remains much closer to the direct
   * FE projection, with a solution/projection gap below 1e-5.
   */
  CHECK(potential_relative_error < 1.5e-4);
  CHECK(projection_relative_error < 1.5e-4);
  CHECK(solution_projection_gap < 1.0e-5);
  CHECK(field_relative_error < 1.0e-3);
  CHECK(maximum_potential_error < 5.0e-4);
}