perf(allocations): reduced overall allocations by 95%, increaseed jacobian applicatin by 2x

This commit uses global pre allocated work space to dramatically reduce memory usage and allocation time

libmeanfield/impl/deformation/safe_newton_step.cpp

@@ -0,0 +1,650 @@
+module;
+
+#include <algorithm>
+#include <array>
+#include <bit>
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <limits>
+#include <mfem.hpp>
+#include <mpi.h>
+#include <numeric>
+#include <stdexcept>
+
+module mean_field;
+
+namespace {
+    using Coefficients = std::array<double, 4>;
+
+    enum class InputFailure : int {
+        none,
+        invalid_options,
+        incompatible_space,
+        invalid_vector,
+        invalid_rule,
+        non_affine_exterior_map
+    };
+
+    enum class EvaluationFailure : int {
+        none,
+        invalid_accepted_mapping,
+        accepted_determinant_below_floor,
+        invalid_mapping_variation,
+        non_finite_polynomial
+    };
+
+    [[nodiscard]] bool vector_is_finite(const mfem::Vector &vector) noexcept {
+        for (int index = 0; index < vector.Size(); ++index) {
+            if (!std::isfinite(vector(index))) {
+                return false;
+            }
+        }
+        return true;
+    }
+
+    [[nodiscard]] bool matrix_is_finite(const mfem::DenseMatrix &matrix) noexcept {
+        for (int row = 0; row < matrix.Height(); ++row) {
+            for (int column = 0; column < matrix.Width(); ++column) {
+                if (!std::isfinite(matrix(row, column))) {
+                    return false;
+                }
+            }
+        }
+        return true;
+    }
+
+    [[nodiscard]] bool
+    options_are_valid(const mean_field::deformation::LargestSafeNewtonStepSizeOptions &options) noexcept {
+        return std::isfinite(options.maximumStepSize) && options.maximumStepSize > 0.0 &&
+               std::isfinite(options.determinantFloor) && options.determinantFloor >= 0.0 &&
+               std::isfinite(options.fractionToBoundarySafety) && options.fractionToBoundarySafety > 0.0 &&
+               options.fractionToBoundarySafety < 1.0;
+    }
+
+    void require_mpi_success(
+        const int status,
+        const char *operation
+    ) {
+        if (status != MPI_SUCCESS) {
+            throw std::runtime_error(operation);
+        }
+    }
+
+    [[nodiscard]] int collective_maximum(
+        const int localValue,
+        const MPI_Comm communicator,
+        const char *operation
+    ) {
+        int globalValue = 0;
+        require_mpi_success(MPI_Allreduce(&localValue, &globalValue, 1, MPI_INT, MPI_MAX, communicator), operation);
+        return globalValue;
+    }
+
+    [[nodiscard]] double selected_entry(
+        const mfem::DenseMatrix &base,
+        const mfem::DenseMatrix &direction,
+        const unsigned int directionColumnMask,
+        const int row,
+        const int column,
+        const double maximumStepSize
+    ) noexcept {
+        if ((directionColumnMask & (1U << static_cast<unsigned int>(column))) != 0U) {
+            return maximumStepSize * direction(row, column);
+        }
+        return base(row, column);
+    }
+
+    [[nodiscard]] double selected_column_determinant(
+        const mfem::DenseMatrix &base,
+        const mfem::DenseMatrix &direction,
+        const unsigned int directionColumnMask,
+        const int dimension,
+        const double maximumStepSize
+    ) noexcept {
+        const auto entry = [&](const int row, const int column) {
+            return selected_entry(base, direction, directionColumnMask, row, column, maximumStepSize);
+        };
+
+        if (dimension == 1) {
+            return entry(0, 0);
+        }
+        if (dimension == 2) {
+            return entry(0, 0) * entry(1, 1) - entry(0, 1) * entry(1, 0);
+        }
+
+        return entry(0, 0) * (entry(1, 1) * entry(2, 2) - entry(1, 2) * entry(2, 1)) -
+               entry(0, 1) * (entry(1, 0) * entry(2, 2) - entry(1, 2) * entry(2, 0)) +
+               entry(0, 2) * (entry(1, 0) * entry(2, 1) - entry(1, 1) * entry(2, 0));
+    }
+
+    [[nodiscard]] Coefficients determinant_polynomial(
+        const mfem::DenseMatrix &base,
+        const mfem::DenseMatrix &direction,
+        const int dimension,
+        const double maximumStepSize,
+        const double determinantFloor
+    ) noexcept {
+        Coefficients coefficients{};
+        const unsigned int termCount = 1U << static_cast<unsigned int>(dimension);
+        for (unsigned int mask = 0; mask < termCount; ++mask) {
+            const int degree = std::popcount(mask);
+            coefficients[static_cast<std::size_t>(degree)] +=
+                selected_column_determinant(base, direction, mask, dimension, maximumStepSize);
+        }
+        coefficients[0] -= determinantFloor;
+        return coefficients;
+    }
+
+    [[nodiscard]] double evaluate_polynomial(
+        const Coefficients &coefficients,
+        const double parameter
+    ) noexcept {
+        return std::fma(
+            parameter, std::fma(parameter, std::fma(parameter, coefficients[3], coefficients[2]), coefficients[1]),
+            coefficients[0]
+        );
+    }
+
+    [[nodiscard]] int polynomial_degree(const Coefficients &coefficients) noexcept {
+        double scale = 0.0;
+        for (const double coefficient : coefficients) {
+            scale = std::max(scale, std::abs(coefficient));
+        }
+        const double tolerance = 64.0 * std::numeric_limits<double>::epsilon() * scale;
+        for (int degree = 3; degree > 0; --degree) {
+            if (std::abs(coefficients[static_cast<std::size_t>(degree)]) > tolerance) {
+                return degree;
+            }
+        }
+        return 0;
+    }
+
+    void append_unit_interval_root(
+        std::array<
+            double,
+            2> &roots,
+        int &rootCount,
+        const double root
+    ) noexcept {
+        if (!std::isfinite(root) || root <= 0.0 || root >= 1.0) {
+            return;
+        }
+        if (rootCount > 0 && std::abs(root - roots[0]) <= 64.0 * std::numeric_limits<double>::epsilon()) {
+            return;
+        }
+        roots[static_cast<std::size_t>(rootCount)] = root;
+        ++rootCount;
+    }
+
+    [[nodiscard]] int derivative_critical_points(
+        const Coefficients &coefficients,
+        const int degree,
+        std::array<
+            double,
+            2> &criticalPoints
+    ) noexcept {
+        int count = 0;
+        if (degree == 2) {
+            append_unit_interval_root(criticalPoints, count, -coefficients[1] / (2.0 * coefficients[2]));
+        } else if (degree == 3) {
+            const double quadratic             = 3.0 * coefficients[3];
+            const double linear                = 2.0 * coefficients[2];
+            const double constant              = coefficients[1];
+            const double discriminant          = std::fma(linear, linear, -4.0 * quadratic * constant);
+            const double discriminantScale     = linear * linear + std::abs(4.0 * quadratic * constant);
+            const double discriminantTolerance = 64.0 * std::numeric_limits<double>::epsilon() * discriminantScale;
+
+            if (discriminant >= -discriminantTolerance) {
+                const double squareRoot = std::sqrt(std::max(0.0, discriminant));
+                if (squareRoot == 0.0) {
+                    append_unit_interval_root(criticalPoints, count, -linear / (2.0 * quadratic));
+                } else {
+                    const double q = -0.5 * (linear + std::copysign(squareRoot, linear));
+                    append_unit_interval_root(criticalPoints, count, q / quadratic);
+                    append_unit_interval_root(criticalPoints, count, constant / q);
+                }
+            }
+        }
+
+        std::sort(criticalPoints.begin(), criticalPoints.begin() + count);
+        return count;
+    }
+
+    [[nodiscard]] double bisect_first_nonpositive_value(
+        const Coefficients &coefficients,
+        double lower,
+        double upper
+    ) noexcept {
+        for (int iteration = 0; iteration < 80; ++iteration) {
+            const double middle = std::midpoint(lower, upper);
+            if (evaluate_polynomial(coefficients, middle) > 0.0) {
+                lower = middle;
+            } else {
+                upper = middle;
+            }
+        }
+        return upper;
+    }
+
+    [[nodiscard]] double first_boundary_parameter(const Coefficients &coefficients) noexcept {
+        const int degree = polynomial_degree(coefficients);
+        if (degree == 0) {
+            return std::numeric_limits<double>::infinity();
+        }
+
+        double coefficientScale = 0.0;
+        for (const double coefficient : coefficients) {
+            coefficientScale += std::abs(coefficient);
+        }
+        const double valueTolerance = 128.0 * std::numeric_limits<double>::epsilon() * coefficientScale;
+
+        std::array<double, 2> criticalPoints{};
+        const int criticalPointCount = derivative_critical_points(coefficients, degree, criticalPoints);
+        std::array<double, 4> intervalEnds{};
+        intervalEnds[0] = 0.0;
+        for (int index = 0; index < criticalPointCount; ++index) {
+            intervalEnds[static_cast<std::size_t>(index + 1)] = criticalPoints[static_cast<std::size_t>(index)];
+        }
+        intervalEnds[static_cast<std::size_t>(criticalPointCount + 1)] = 1.0;
+
+        for (int interval = 0; interval <= criticalPointCount; ++interval) {
+            const double lower      = intervalEnds[static_cast<std::size_t>(interval)];
+            const double upper      = intervalEnds[static_cast<std::size_t>(interval + 1)];
+            const double upperValue = evaluate_polynomial(coefficients, upper);
+            if (upperValue <= 0.0) {
+                return bisect_first_nonpositive_value(coefficients, lower, upper);
+            }
+            if (upperValue <= valueTolerance) {
+                // A repeated root only touches zero.  Floating-point evaluation
+                // at the derivative root may land a few ulps above it.
+                return upper;
+            }
+        }
+        return std::numeric_limits<double>::infinity();
+    }
+
+    void true_to_local(
+        const mfem::ParFiniteElementSpace &finiteElementSpace,
+        const mfem::Vector &trueVector,
+        mfem::Vector &localVector
+    ) {
+        localVector.SetSize(finiteElementSpace.GetVSize());
+        const mfem::Operator *prolongation = finiteElementSpace.GetProlongationMatrix();
+        if (prolongation != nullptr) {
+            prolongation->Mult(trueVector, localVector);
+        } else {
+            localVector = trueVector;
+        }
+    }
+} // namespace
+
+namespace mean_field::deformation {
+    void LargestSafeNewtonStepSizeOptions::Validate() const {
+        if (!std::isfinite(maximumStepSize) || maximumStepSize <= 0.0) {
+            throw std::invalid_argument("A geometry preflight requires a finite, positive maximum step size.");
+        }
+        if (!std::isfinite(determinantFloor) || determinantFloor < 0.0) {
+            throw std::invalid_argument("A geometry preflight requires a finite, non-negative determinant floor.");
+        }
+        if (!std::isfinite(fractionToBoundarySafety) || fractionToBoundarySafety <= 0.0 ||
+            fractionToBoundarySafety >= 1.0) {
+            throw std::invalid_argument(
+                "A geometry preflight requires a finite fraction-to-boundary safety factor strictly between zero "
+                "and one."
+            );
+        }
+    }
+
+    LargestSafeNewtonStepSizeEstimate estimate_largest_safe_newton_step_size(
+        const mapping::DomainMapper &domainMapper,
+        const mfem::ParFiniteElementSpace &displacementSpace,
+        const mfem::ParGridFunction &compactificationCoordinate,
+        const mfem::Vector &acceptedVolumeDisplacement,
+        const mfem::Vector &volumeNewtonDirection,
+        const std::span<const NewtonStepGeometryRule> geometryRules,
+        const LargestSafeNewtonStepSizeOptions &options
+    ) {
+        const MPI_Comm communicator = displacementSpace.GetComm();
+        if (communicator == MPI_COMM_NULL) {
+            throw std::invalid_argument("A geometry preflight requires a valid displacement communicator.");
+        }
+
+        const mfem::FiniteElementSpace *compactificationSpace = compactificationCoordinate.FESpace();
+        mfem::Mesh *mesh                                      = displacementSpace.GetMesh();
+        InputFailure localInputFailure                        = InputFailure::none;
+        const auto recordInputFailure                         = [&](const InputFailure failure) {
+            localInputFailure =
+                static_cast<InputFailure>(std::max(static_cast<int>(localInputFailure), static_cast<int>(failure)));
+        };
+
+        if (!options_are_valid(options)) {
+            recordInputFailure(InputFailure::invalid_options);
+        }
+        const int dimension                 = domainMapper.GetDimension();
+        const mfem::Ordering::Type ordering = displacementSpace.GetOrdering();
+        if (mesh == nullptr || compactificationSpace == nullptr || compactificationSpace->GetMesh() != mesh ||
+            dimension < 1 || dimension > 3 || (mesh != nullptr && mesh->SpaceDimension() != dimension) ||
+            displacementSpace.GetVDim() != dimension ||
+            (compactificationSpace != nullptr && compactificationSpace->GetVDim() != 1) ||
+            (compactificationSpace != nullptr &&
+             compactificationCoordinate.Size() != compactificationSpace->GetVSize()) ||
+            (ordering != mfem::Ordering::byNODES && ordering != mfem::Ordering::byVDIM)) {
+            recordInputFailure(InputFailure::incompatible_space);
+        }
+        if (acceptedVolumeDisplacement.Size() != displacementSpace.GetTrueVSize() ||
+            volumeNewtonDirection.Size() != displacementSpace.GetTrueVSize() ||
+            !vector_is_finite(acceptedVolumeDisplacement) || !vector_is_finite(volumeNewtonDirection)) {
+            recordInputFailure(InputFailure::invalid_vector);
+        }
+        if (geometryRules.size() > static_cast<std::size_t>(std::numeric_limits<int>::max())) {
+            recordInputFailure(InputFailure::invalid_rule);
+        }
+
+        std::uint64_t localPointCount = 0;
+        if (mesh != nullptr) {
+            for (const NewtonStepGeometryRule &entry : geometryRules) {
+                if (entry.element < 0 || entry.element >= mesh->GetNE() || entry.integrationRule == nullptr ||
+                    entry.integrationRule->GetNPoints() <= 0) {
+                    recordInputFailure(InputFailure::invalid_rule);
+                    continue;
+                }
+                localPointCount += static_cast<std::uint64_t>(entry.integrationRule->GetNPoints());
+                mfem::ElementTransformation *transformation    = mesh->GetElementTransformation(entry.element);
+                const mfem::FiniteElement *displacementElement = displacementSpace.GetFE(entry.element);
+                const mfem::FiniteElement *compactificationElement =
+                    compactificationSpace != nullptr ? compactificationSpace->GetFE(entry.element) : nullptr;
+                if (transformation == nullptr || displacementElement == nullptr || compactificationElement == nullptr ||
+                    transformation->GetSpaceDim() != dimension || displacementElement->GetDim() != dimension ||
+                    compactificationElement->GetDim() != dimension ||
+                    displacementElement->GetGeomType() != compactificationElement->GetGeomType() ||
+                    displacementElement->GetRangeType() != mfem::FiniteElement::SCALAR ||
+                    displacementElement->GetMapType() != mfem::FiniteElement::VALUE ||
+                    displacementElement->GetDerivType() != mfem::FiniteElement::GRAD ||
+                    compactificationElement->GetRangeType() != mfem::FiniteElement::SCALAR ||
+                    compactificationElement->GetMapType() != mfem::FiniteElement::VALUE ||
+                    compactificationElement->GetDerivType() != mfem::FiniteElement::GRAD) {
+                    recordInputFailure(InputFailure::invalid_rule);
+                } else if (
+                    domainMapper.IsCompactifiedElement(*transformation) &&
+                    !domainMapper.GetExteriorMap().IsAffineInDisplacement()
+                ) {
+                    recordInputFailure(InputFailure::non_affine_exterior_map);
+                }
+            }
+        }
+
+        const int globalInputFailure = collective_maximum(
+            static_cast<int>(localInputFailure), communicator,
+            "The geometry preflight could not validate its distributed inputs."
+        );
+        if (globalInputFailure != static_cast<int>(InputFailure::none)) {
+            switch (static_cast<InputFailure>(globalInputFailure)) {
+            case InputFailure::invalid_options:
+                throw std::invalid_argument("The geometry preflight options are invalid on at least one rank.");
+            case InputFailure::incompatible_space:
+                throw std::invalid_argument(
+                    "The geometry preflight requires compatible displacement and compactification spaces in one to "
+                    "three dimensions."
+                );
+            case InputFailure::invalid_vector:
+                throw std::invalid_argument(
+                    "The geometry preflight received an incompatible or non-finite true-DOF displacement vector."
+                );
+            case InputFailure::invalid_rule:
+                throw std::invalid_argument("The geometry preflight received an invalid local quadrature rule.");
+            case InputFailure::non_affine_exterior_map:
+                throw std::invalid_argument(
+                    "The geometry preflight requires compactified mappings that are affine in displacement."
+                );
+            case InputFailure::none:
+                break;
+            }
+        }
+
+        const std::array<double, 3> localOptions{
+            options.maximumStepSize, options.determinantFloor, options.fractionToBoundarySafety
+        };
+        std::array<double, 3> minimumOptions{};
+        std::array<double, 3> maximumOptions{};
+        require_mpi_success(
+            MPI_Allreduce(
+                localOptions.data(), minimumOptions.data(), static_cast<int>(localOptions.size()), MPI_DOUBLE, MPI_MIN,
+                communicator
+            ),
+            "The geometry preflight could not compare its distributed options."
+        );
+        require_mpi_success(
+            MPI_Allreduce(
+                localOptions.data(), maximumOptions.data(), static_cast<int>(localOptions.size()), MPI_DOUBLE, MPI_MAX,
+                communicator
+            ),
+            "The geometry preflight could not compare its distributed options."
+        );
+        if (minimumOptions != maximumOptions) {
+            throw std::invalid_argument("The geometry preflight requires identical options on every rank.");
+        }
+
+        std::uint64_t globalPointCount = 0;
+        require_mpi_success(
+            MPI_Allreduce(&localPointCount, &globalPointCount, 1, MPI_UINT64_T, MPI_SUM, communicator),
+            "The geometry preflight could not count its distributed samples."
+        );
+        if (globalPointCount == 0) {
+            throw std::invalid_argument("The geometry preflight requires at least one quadrature point globally.");
+        }
+
+        mfem::Vector acceptedLocal;
+        mfem::Vector directionLocal;
+        true_to_local(displacementSpace, acceptedVolumeDisplacement, acceptedLocal);
+        true_to_local(displacementSpace, volumeNewtonDirection, directionLocal);
+
+        mapping::DomainMapper::Workspace workspace(dimension);
+        mapping::MappingPointContext mappingContext;
+        mapping::MappingPointVariation mappingVariation;
+        mfem::Array<int> displacementDofs;
+        mfem::Array<int> compactificationDofs;
+        mfem::Vector elementAcceptedDisplacement;
+        mfem::Vector elementDirection;
+        mfem::Vector elementCompactification;
+
+        double localBoundaryStep      = std::numeric_limits<double>::infinity();
+        double localMinimumAtAccepted = std::numeric_limits<double>::infinity();
+        double localMinimumAtMaximum  = std::numeric_limits<double>::infinity();
+        Coefficients localLimitingCoefficients{};
+        int localLimitingElement                 = -1;
+        int localLimitingRule                    = -1;
+        int localLimitingPoint                   = -1;
+        EvaluationFailure localEvaluationFailure = EvaluationFailure::none;
+        const auto recordEvaluationFailure       = [&](const EvaluationFailure failure) {
+            localEvaluationFailure = static_cast<EvaluationFailure>(
+                std::max(static_cast<int>(localEvaluationFailure), static_cast<int>(failure))
+            );
+        };
+
+        for (std::size_t ruleIndex = 0; ruleIndex < geometryRules.size(); ++ruleIndex) {
+            const NewtonStepGeometryRule &entry         = geometryRules[ruleIndex];
+            mfem::ElementTransformation *transformation = mesh->GetElementTransformation(entry.element);
+            mfem::DofTransformation *displacementDofTransformation =
+                displacementSpace.GetElementVDofs(entry.element, displacementDofs);
+            mfem::DofTransformation *compactificationDofTransformation =
+                compactificationSpace->GetElementDofs(entry.element, compactificationDofs);
+
+            acceptedLocal.GetSubVector(displacementDofs, elementAcceptedDisplacement);
+            directionLocal.GetSubVector(displacementDofs, elementDirection);
+            compactificationCoordinate.GetSubVector(compactificationDofs, elementCompactification);
+            if (displacementDofTransformation != nullptr) {
+                displacementDofTransformation->InvTransformPrimal(elementAcceptedDisplacement);
+                displacementDofTransformation->InvTransformPrimal(elementDirection);
+            }
+            if (compactificationDofTransformation != nullptr) {
+                compactificationDofTransformation->InvTransformPrimal(elementCompactification);
+            }
+
+            const mfem::FiniteElement &displacementElement     = *displacementSpace.GetFE(entry.element);
+            const mfem::FiniteElement &compactificationElement = *compactificationSpace->GetFE(entry.element);
+            const mapping::ElementDisplacementData acceptedData(
+                displacementElement, elementAcceptedDisplacement, displacementSpace.GetOrdering()
+            );
+            const mapping::ElementDisplacementData directionData(
+                displacementElement, elementDirection, displacementSpace.GetOrdering()
+            );
+            const mapping::ElementCompactificationData compactificationData(
+                compactificationElement, elementCompactification
+            );
+            const mapping::ElementMappingData elementData{
+                .displacement = acceptedData, .compactification = compactificationData
+            };
+
+            for (int point = 0; point < entry.integrationRule->GetNPoints(); ++point) {
+                const mfem::IntegrationPoint &integrationPoint = entry.integrationRule->IntPoint(point);
+                const mapping::MappingStatus mappingStatus     = domainMapper.EvaluatePoint(
+                    elementData, *transformation, integrationPoint, workspace, mappingContext
+                );
+                if (mappingStatus != mapping::MappingStatus::valid) {
+                    recordEvaluationFailure(EvaluationFailure::invalid_accepted_mapping);
+                    continue;
+                }
+                if (mappingContext.mapping_determinant <= options.determinantFloor) {
+                    recordEvaluationFailure(EvaluationFailure::accepted_determinant_below_floor);
+                    continue;
+                }
+
+                const mapping::MappingStatus variationStatus = domainMapper.EvaluatePointVariation(
+                    elementData, directionData, *transformation, integrationPoint, mappingContext, workspace,
+                    mappingVariation
+                );
+                if (variationStatus != mapping::MappingStatus::valid) {
+                    recordEvaluationFailure(EvaluationFailure::invalid_mapping_variation);
+                    continue;
+                }
+                if (!matrix_is_finite(mappingContext.mapping_jacobian) ||
+                    !matrix_is_finite(mappingVariation.mapping_jacobian_variation)) {
+                    recordEvaluationFailure(EvaluationFailure::invalid_mapping_variation);
+                    continue;
+                }
+
+                Coefficients coefficients = determinant_polynomial(
+                    mappingContext.mapping_jacobian, mappingVariation.mapping_jacobian_variation, dimension,
+                    options.maximumStepSize, options.determinantFloor
+                );
+                // Use the mapper's own determinant at the accepted state to
+                // avoid a second, slightly different round-off path.
+                coefficients[0]                   = mappingContext.mapping_determinant - options.determinantFloor;
+                const double determinantAtMaximum = evaluate_polynomial(coefficients, 1.0) + options.determinantFloor;
+                if (!std::isfinite(determinantAtMaximum)) {
+                    recordEvaluationFailure(EvaluationFailure::non_finite_polynomial);
+                    continue;
+                }
+
+                localMinimumAtAccepted         = std::min(localMinimumAtAccepted, mappingContext.mapping_determinant);
+                localMinimumAtMaximum          = std::min(localMinimumAtMaximum, determinantAtMaximum);
+
+                const double boundaryParameter = first_boundary_parameter(coefficients);
+                if (std::isfinite(boundaryParameter)) {
+                    const double boundaryStep = options.maximumStepSize * boundaryParameter;
+                    if (boundaryStep < localBoundaryStep) {
+                        localBoundaryStep         = boundaryStep;
+                        localLimitingCoefficients = coefficients;
+                        localLimitingElement      = entry.element;
+                        localLimitingRule         = static_cast<int>(ruleIndex);
+                        localLimitingPoint        = point;
+                    }
+                }
+            }
+        }
+
+        const int globalEvaluationFailure = collective_maximum(
+            static_cast<int>(localEvaluationFailure), communicator,
+            "The geometry preflight could not combine its distributed mapping status."
+        );
+        if (globalEvaluationFailure != static_cast<int>(EvaluationFailure::none)) {
+            switch (static_cast<EvaluationFailure>(globalEvaluationFailure)) {
+            case EvaluationFailure::invalid_accepted_mapping:
+                throw std::domain_error(
+                    "The geometry preflight received an accepted displacement with an invalid mapped geometry."
+                );
+            case EvaluationFailure::accepted_determinant_below_floor:
+                throw std::domain_error(
+                    "The accepted displacement does not lie strictly above the requested determinant floor."
+                );
+            case EvaluationFailure::invalid_mapping_variation:
+                throw std::domain_error("The geometry preflight could not evaluate the mapping direction.");
+            case EvaluationFailure::non_finite_polynomial:
+                throw std::domain_error("The geometry preflight produced a non-finite determinant polynomial.");
+            case EvaluationFailure::none:
+                break;
+            }
+        }
+
+        double globalMinimumAtAccepted = 0.0;
+        double globalMinimumAtMaximum  = 0.0;
+        require_mpi_success(
+            MPI_Allreduce(&localMinimumAtAccepted, &globalMinimumAtAccepted, 1, MPI_DOUBLE, MPI_MIN, communicator),
+            "The geometry preflight could not reduce its accepted-state determinant."
+        );
+        require_mpi_success(
+            MPI_Allreduce(&localMinimumAtMaximum, &globalMinimumAtMaximum, 1, MPI_DOUBLE, MPI_MIN, communicator),
+            "The geometry preflight could not reduce its maximum-step determinant."
+        );
+
+        int rank = 0;
+        require_mpi_success(MPI_Comm_rank(communicator, &rank), "The geometry preflight could not identify its rank.");
+        struct BoundaryLocation {
+            double step;
+            int rank;
+        };
+        const BoundaryLocation localLocation{.step = localBoundaryStep, .rank = rank};
+        BoundaryLocation globalLocation{};
+        require_mpi_success(
+            MPI_Allreduce(&localLocation, &globalLocation, 1, MPI_DOUBLE_INT, MPI_MINLOC, communicator),
+            "The geometry preflight could not select its limiting point."
+        );
+
+        const bool limitedByGeometry = std::isfinite(globalLocation.step);
+        std::array<int, 3> limitingLocation{-1, -1, -1};
+        Coefficients limitingCoefficients{};
+        if (limitedByGeometry) {
+            if (rank == globalLocation.rank) {
+                limitingLocation     = {localLimitingElement, localLimitingRule, localLimitingPoint};
+                limitingCoefficients = localLimitingCoefficients;
+            }
+            require_mpi_success(
+                MPI_Bcast(
+                    limitingLocation.data(), static_cast<int>(limitingLocation.size()), MPI_INT, globalLocation.rank,
+                    communicator
+                ),
+                "The geometry preflight could not broadcast its limiting location."
+            );
+            require_mpi_success(
+                MPI_Bcast(
+                    limitingCoefficients.data(), static_cast<int>(limitingCoefficients.size()), MPI_DOUBLE,
+                    globalLocation.rank, communicator
+                ),
+                "The geometry preflight could not broadcast its limiting polynomial."
+            );
+        }
+
+        const double boundaryStepSize = limitedByGeometry ? globalLocation.step : options.maximumStepSize;
+        const double stepSize =
+            limitedByGeometry ? options.fractionToBoundarySafety * boundaryStepSize : options.maximumStepSize;
+        const double limitingPointDeterminantAtStepSize =
+            limitedByGeometry ? evaluate_polynomial(limitingCoefficients, stepSize / options.maximumStepSize) +
+                                    options.determinantFloor
+                              : globalMinimumAtMaximum;
+
+        return {
+            .stepSize                            = stepSize,
+            .boundaryStepSize                    = boundaryStepSize,
+            .minimumDeterminantAtAcceptedState   = globalMinimumAtAccepted,
+            .minimumDeterminantAtMaximumStepSize = globalMinimumAtMaximum,
+            .limitingPointDeterminantAtStepSize  = limitingPointDeterminantAtStepSize,
+            .sampledQuadraturePointCount         = globalPointCount,
+            .limitedByGeometry                   = limitedByGeometry,
+            .limitingRank                        = limitedByGeometry ? globalLocation.rank : -1,
+            .limitingElement                     = limitingLocation[0],
+            .limitingRule                        = limitingLocation[1],
+            .limitingQuadraturePoint             = limitingLocation[2]
+        };
+    }
+} // namespace mean_field::deformation

libmeanfield/impl/fem/reference_tables.cpp

@@ -0,0 +1,159 @@
+module;
+
+#include <array>
+#include <cmath>
+#include <functional>
+#include <map>
+#include <memory>
+#include <mutex>
+#include <stdexcept>
+#include <utility>
+#include <vector>
+
+#include <mfem.hpp>
+
+module mean_field;
+import :fem.reference_tables;
+
+namespace mean_field::fem {
+    namespace {
+        struct ReferenceTableKey {
+            const mfem::FiniteElement *element;
+            std::vector<std::array<double, 4>> points;
+
+            bool operator<(const ReferenceTableKey &other) const {
+                if (element != other.element)
+                    return std::less<const mfem::FiniteElement *>{}(element, other.element);
+                return points < other.points;
+            }
+        };
+
+        ReferenceTableKey make_key(
+            const mfem::FiniteElement &element,
+            const mfem::IntegrationRule &rule
+        ) {
+            ReferenceTableKey key{.element = &element, .points = {}};
+            key.points.reserve(rule.GetNPoints());
+            for (int q = 0; q < rule.GetNPoints(); ++q) {
+                const auto &point = rule.IntPoint(q);
+                const std::array<double, 4> values{
+                    point.x, element.GetDim() > 1 ? point.y : 0.0, element.GetDim() > 2 ? point.z : 0.0, point.weight
+                };
+                for (const double value : values) {
+                    if (!std::isfinite(value))
+                        throw std::invalid_argument("Reference table quadrature entries must be finite.");
+                }
+                key.points.push_back(values);
+            }
+            return key;
+        }
+    } // namespace
+
+    struct ReferenceTableCache::Storage {
+        std::mutex mutex;
+        std::map<ReferenceTableKey, std::shared_ptr<const ScalarReferenceTable>> scalar_tables;
+        std::map<ReferenceTableKey, std::shared_ptr<const VectorReferenceTable>> vector_tables;
+    };
+
+    ReferenceTableCache::ReferenceTableCache() : m_storage(std::make_unique<Storage>()) {
+    }
+    ReferenceTableCache::~ReferenceTableCache() = default;
+
+    std::shared_ptr<const ScalarReferenceTable> ReferenceTableCache::GetScalarTable(
+        const mfem::FiniteElement &element,
+        const mfem::IntegrationRule &rule
+    ) const {
+        if (element.GetRangeType() != mfem::FiniteElement::SCALAR)
+            throw std::invalid_argument("A scalar reference table requires a scalar finite element.");
+        auto key = make_key(element, rule);
+        const std::lock_guard lock(m_storage->mutex);
+        if (const auto found = m_storage->scalar_tables.find(key); found != m_storage->scalar_tables.end())
+            return found->second;
+        auto table = std::shared_ptr<const ScalarReferenceTable>(new ScalarReferenceTable(element, rule));
+        m_storage->scalar_tables.emplace(std::move(key), table);
+        return table;
+    }
+
+    std::shared_ptr<const VectorReferenceTable> ReferenceTableCache::GetVectorTable(
+        const mfem::FiniteElement &element,
+        const mfem::IntegrationRule &rule
+    ) const {
+        if (element.GetRangeType() != mfem::FiniteElement::VECTOR)
+            throw std::invalid_argument("A vector reference table requires a vector finite element.");
+        auto key = make_key(element, rule);
+        const std::lock_guard lock(m_storage->mutex);
+        if (const auto found = m_storage->vector_tables.find(key); found != m_storage->vector_tables.end())
+            return found->second;
+        auto table = std::shared_ptr<const VectorReferenceTable>(new VectorReferenceTable(element, rule));
+        m_storage->vector_tables.emplace(std::move(key), table);
+        return table;
+    }
+
+    ScalarReferenceTable::ScalarReferenceTable(
+        const mfem::FiniteElement &element,
+        const mfem::IntegrationRule &rule
+    )
+        : m_values(
+              rule.GetNPoints(),
+              element.GetDof()
+          ),
+          m_dimension(element.GetDim()) {
+        mfem::Vector values(element.GetDof());
+        if (element.GetDerivType() == mfem::FiniteElement::GRAD)
+            m_gradients.resize(rule.GetNPoints());
+        for (int q = 0; q < rule.GetNPoints(); ++q) {
+            const auto &point = rule.IntPoint(q);
+            element.CalcShape(point, values);
+            for (int dof = 0; dof < element.GetDof(); ++dof)
+                m_values(q, dof) = values(dof);
+            if (!m_gradients.empty()) {
+                auto &gradient = m_gradients[q];
+                gradient.SetSize(element.GetDof(), m_dimension);
+                element.CalcDShape(point, gradient);
+            }
+        }
+    }
+
+    const mfem::DenseMatrix &ScalarReferenceTable::GetValues() const {
+        return m_values;
+    }
+    const mfem::DenseMatrix &ScalarReferenceTable::GetGradients(const int point) const {
+        return m_gradients.at(point);
+    }
+    int ScalarReferenceTable::GetPointCount() const {
+        return m_values.Height();
+    }
+    int ScalarReferenceTable::GetDofCount() const {
+        return m_values.Width();
+    }
+    int ScalarReferenceTable::GetDimension() const {
+        return m_dimension;
+    }
+
+    VectorReferenceTable::VectorReferenceTable(
+        const mfem::FiniteElement &element,
+        const mfem::IntegrationRule &rule
+    )
+        : m_dof_count(element.GetDof()),
+          m_dimension(element.GetRangeDim()) {
+        m_values.resize(rule.GetNPoints());
+        for (int q = 0; q < rule.GetNPoints(); ++q) {
+            auto &values = m_values[q];
+            values.SetSize(m_dof_count, m_dimension);
+            element.CalcVShape(rule.IntPoint(q), values);
+        }
+    }
+
+    const mfem::DenseMatrix &VectorReferenceTable::GetValues(const int point) const {
+        return m_values.at(point);
+    }
+    int VectorReferenceTable::GetPointCount() const {
+        return static_cast<int>(m_values.size());
+    }
+    int VectorReferenceTable::GetDofCount() const {
+        return m_dof_count;
+    }
+    int VectorReferenceTable::GetDimension() const {
+        return m_dimension;
+    }
+} // namespace mean_field::fem

libmeanfield/impl/mapping/compactification/kelvin.cpp

@@ -54,6 +54,7 @@ namespace mean_field::mapping::compactification {
         if (!std::isfinite(compactification_coordinate))
             return MappingStatus::non_finite_input;
 
+        // How close we will allow the code to get to compactified infinity
         const double tolerance = m_options.coordinate_tolerance;
 
         if (compactification_coordinate < -tolerance || compactification_coordinate > 1.0 + tolerance) {
@@ -68,13 +69,21 @@ namespace mean_field::mapping::compactification {
             return MappingStatus::at_compactified_infinity;
         }
 
+        // Here we need to do some transformations from the options defined on the mesh to useful computational coordinates
+        // r_inf_ref is the computational / reference radius of the infinity surface (the edge of the entire domain) and r_star_ref is the radius of the spherical
+        // stellar model inscribed within. Therefore radial extent is the computational radial distance between the stellar surface and the infinity surface.
+        // Note that this is separate from the parameterize compactification coordinate.
         const double radial_extent        = m_options.r_inf_ref - m_options.r_star_ref;
+
+        // This places us at the correct spot in computational space given the current compactification coordinate. Say you have compactification = 0.5,
+        // an r_star_ref of 2 and a radial extent of 3, this this will place you at 2 + 0.5 * 3 = 3.5 in computational space, which is half way between the stellar surface and the infinity surface.
         const double computational_radius = m_options.r_star_ref + coordinate * radial_extent;
 
         if (!std::isfinite(computational_radius) || computational_radius <= 0.0) {
             return MappingStatus::invalid_reference_radius;
         }
 
+        // Invert the exterior coordinate so it runs from 0 at the star to 1 at compactified infinity
         const double one_minus_coordinate = 1.0 - coordinate;
         const double denominator          = computational_radius * one_minus_coordinate;
 
@@ -82,6 +91,11 @@ namespace mean_field::mapping::compactification {
             return MappingStatus::non_finite_result;
         }
 
+        // The scale here is the factor which stretches the finite computational domain into the infinite physical domain. Properties we need this to have
+        // include that it should go to 1 at the stellar surface and go to infinity at the compactified infinity.
+        // Mathematically this is scale = |r|/|x| where r is the physical radius and x is the computational radius.
+        // Put another way, scale is the ratio of the target physical radius for the current compactification coordinate
+        // to the current mesh radius.
         const double scale            = m_options.r_star_ref / denominator;
         const double scale_derivative = scale * (1.0 / one_minus_coordinate - radial_extent / computational_radius);
 
@@ -119,6 +133,8 @@ namespace mean_field::mapping::compactification {
         }
 
         RadialFactors factors;
+
+        // The key radial factors we use are the scale (which stretches the finite computational domain into the infinite physical domain) and the scale derivative (which is used to compute the mapping jacobian).
         const MappingStatus factor_status = ComputeRadialFactors(input.compactification_coordinate, factors);
         if (factor_status != MappingStatus::valid)
             return factor_status;
@@ -127,12 +143,21 @@ namespace mean_field::mapping::compactification {
         result.mapping_jacobian.SetSize(dimension, dimension);
 
         for (int i = 0; i < dimension; ++i) {
+            // Note how the physical position is just the product of the displaced position and the scale factor.
             result.physical_position(i) = factors.scale * input.displaced_position(i);
 
+            // The mapping jacobian comes from trivial application of the product rule
+            // recall: r_i = scale * x_i where r is the physical position and x is the displaced position.
+            // then we can differentiate wrt. X_j holding nothing fixed. Note the capital X here, this is the mesh coordinate not the displaced position.
+            // Lets call this jacobian F
+            // F = \frac{\partial r_i}{\partial X_{j}}
+            // F then tells us how the physical position changes as we move along mesh coordinates
+            // Lets then apply this to the function we have for the kelvin compactification
+            // F = scale * \frac{\partial x_i}{\partial X_j} + x_i * \frac{\partial scale}{\partial X_j}
+            // Below you can see the displacement jacobian (\frac{\partial x_i}{\partial X_j}) and the scale derivative (\frac{\partial scale}{\partial X_j}) being applied to compute the mapping jacobian.
             for (int j = 0; j < dimension; ++j) {
                 const double scale_gradient = factors.scale_derivative * input.compactification_coordinate_gradient(j);
-                result.mapping_jacobian(i, j) =
-                    factors.scale * input.displacement_jacobian(i, j) + input.displaced_position(i) * scale_gradient;
+                result.mapping_jacobian(i, j) = factors.scale * input.displacement_jacobian(i, j) + input.displaced_position(i) * scale_gradient;
             }
         }
 

libmeanfield/impl/mapping/domain_mapper.cpp

@@ -167,6 +167,7 @@ namespace mean_field::mapping {
 
         m_field_value.SetSize(dimension);
         m_field_jacobian.SetSize(dimension, dimension);
+        m_reference_field_jacobian.SetSize(dimension, dimension);
 
         m_compactification_point.coordinate = 0.0;
         m_compactification_point.coordinate_gradient.SetSize(dimension);
@@ -538,6 +539,10 @@ namespace mean_field::mapping {
         return MappingStatus::valid;
     }
 
+    /**
+     * @brief Evaluate a displacement field dof matrix at a given integration point and compute what the displacement of that point is and what the gradient of the the displacement is with respect to the computational coordinates / reference frame.
+     * @note There is actually nothing in this function preventing some field other than displacement from being passed through here; this should maybe be tightened.
+     */
     void DomainMapper::EvaluateField(
         const ElementDisplacementData &field,
         mfem::ElementTransformation &transformation,
@@ -551,22 +556,30 @@ namespace mean_field::mapping {
         const mfem::DenseMatrix &dof_matrix = field.GetDofMatrix();
 
         workspace.m_shape.SetSize(element.GetDof());
-        workspace.m_mesh_dshape.SetSize(element.GetDof(), m_options.dimension);
 
         element.CalcShape(integration_point, workspace.m_shape);
-        if (inverse_mesh_jacobian != nullptr) {
-            workspace.m_reference_dshape.SetSize(element.GetDof(), m_options.dimension);
-            element.CalcDShape(integration_point, workspace.m_reference_dshape);
-            mfem::Mult(workspace.m_reference_dshape, *inverse_mesh_jacobian, workspace.m_mesh_dshape);
-        } else {
-            element.CalcPhysDShape(transformation, workspace.m_mesh_dshape);
-        }
-
         value.SetSize(m_options.dimension);
         dof_matrix.MultTranspose(workspace.m_shape, value);
-
         jacobian.SetSize(m_options.dimension, m_options.dimension);
-        mfem::MultAtB(dof_matrix, workspace.m_mesh_dshape, jacobian);
+
+        if (inverse_mesh_jacobian != nullptr || element.GetMapType() == mfem::FiniteElement::VALUE) {
+            workspace.m_reference_dshape.SetSize(element.GetDof(), m_options.dimension);
+            element.CalcDShape(integration_point, workspace.m_reference_dshape);
+            // Contract DOFs before applying fixed-mesh geometry. This is the
+            // same DOF^T * (Dshape * J_mesh^-1), without transforming every
+            // basis gradient. The scratch matrix must not alias the cached
+            // inverse supplied by EvaluateVolumeVariation.
+            mfem::MultAtB(dof_matrix, workspace.m_reference_dshape, workspace.m_reference_field_jacobian);
+            const mfem::DenseMatrix &inverseMeshJacobian =
+                inverse_mesh_jacobian != nullptr ? *inverse_mesh_jacobian : transformation.InverseJacobian();
+            mfem::Mult(workspace.m_reference_field_jacobian, inverseMeshJacobian, jacobian);
+        } else {
+            // Retain the original finite-element-specific physical-gradient
+            // path for mapping types without the ordinary VALUE pullback.
+            workspace.m_mesh_dshape.SetSize(element.GetDof(), m_options.dimension);
+            element.CalcPhysDShape(transformation, workspace.m_mesh_dshape);
+            mfem::MultAtB(dof_matrix, workspace.m_mesh_dshape, jacobian);
+        }
     }
 
     MappingStatus DomainMapper::EvaluatePoint(
@@ -594,6 +607,7 @@ namespace mean_field::mapping {
         context.reference_position.SetSize(m_options.dimension);
         transformation.Transform(integration_point, context.reference_position);
 
+        // Get the displacement field value and its Jacobian at the integration point. Note these are in the workspace to avoid repeated allocations.
         EvaluateField(
             element_data.displacement, transformation, integration_point, workspace, workspace.m_field_value,
             workspace.m_field_jacobian, nullptr
@@ -605,17 +619,32 @@ namespace mean_field::mapping {
         }
 
         context.displaced_position.SetSize(m_options.dimension);
+
+        // Get the position of the point in physical space by adding the displacement to the reference position. Note MFEM really dislikes raw arithmetic operators
+        // so we need to first assign the reference position then use the in place += operator.
         context.displaced_position = context.reference_position;
         context.displaced_position += workspace.m_field_value;
 
         context.displacement_jacobian.SetSize(m_options.dimension, m_options.dimension);
         context.displacement_jacobian = workspace.m_field_jacobian;
+
+        // Ensure that the diagonal of the displacement Jacobian is incremented by 1.0 to account for the identity mapping from reference to physical space.
+        // recall that r = x + d (where d is the workspace.m_field_value and x is context.reference_position) then we can differentiate this
+        // component wise to find the gradient of the displaced position wrt. the mesh coordinate (reference position). E.g as you move along
+        // the mesh coordinate how much does the physical coordinate change and in what direction. Lets call this F
+        // F = \frac{\partial r_i}{\partial x_j} where r is the displaced position and x is the mesh position.
+        // We then have F = \frac{\partial x_i}{\partial x_j} + \frac{\partial d_i}{x_j} where d is the displacement (recall r = x + d)
+        // By definition the first term is the identity matrix. The second term we get out of EvaluateField. Thus why we need to add the identity matrix here
         for (int i = 0; i < m_options.dimension; ++i)
             context.displacement_jacobian(i, i) += 1.0;
 
         context.compactified = IsCompactifiedElement(transformation);
 
+        // This branch only runs for vacuum elements
         if (context.compactified) {
+            // There are two things that we need to the mapping. First is a reference coordinate which stroid embeds into the mesh at mesh generation time, this is
+            // parameterized from 0 - 1 where 0 is the model surface and 1 is the mesh exterior (what will becomes the compactified infinity, note also we never actually evaluate at s=1; rather we define some arbitrary small tolerance to approach s=1). Lets call this s. We also need
+            // the gradient of s as we move along the mesh coordinates. All of this is stashes within workspace.m_compactification_point.
             const MappingStatus coordinate_status = EvaluateCompactificationCoordinate(
                 element_data.compactification, transformation, integration_point, workspace,
                 workspace.m_compactification_point, nullptr
@@ -632,6 +661,8 @@ namespace mean_field::mapping {
                 .compactification_coordinate_gradient = workspace.m_compactification_point.coordinate_gradient
             };
 
+            // This apply whatever the exterior map is to generate the new physical exterior coordinate and jacobian between physical and reference space.
+            // In general we have only implemented a kelvin mapping; however, in future additional mappings may be implemented.
             const MappingStatus exterior_status = m_exterior_map->Evaluate(exterior_input, workspace.m_exterior_result);
             if (exterior_status != MappingStatus::valid)
                 return exterior_status;
@@ -643,6 +674,7 @@ namespace mean_field::mapping {
             context.mapping_jacobian  = context.displacement_jacobian;
         }
 
+        // Validation work
         if (!vector_is_finite(context.physical_position) || !matrix_is_finite(context.mapping_jacobian))
             return MappingStatus::non_finite_result;
 
@@ -650,9 +682,12 @@ namespace mean_field::mapping {
         if (!std::isfinite(context.mapping_determinant))
             return MappingStatus::non_finite_result;
         if (context.mapping_determinant <= 0.0)
+            // This is the most common error we see come out of this function, specifically it is common when we try to deform the mesh too much in one step.
             return MappingStatus::non_positive_determinant;
 
         context.inverse_mapping_jacobian.SetSize(m_options.dimension, m_options.dimension);
+
+        // It can be useful to have the inverse jacobian, here we just use MFEM's build in inverse tooling.
         mfem::CalcInverse(context.mapping_jacobian, context.inverse_mapping_jacobian);
 
         if (!matrix_is_finite(context.inverse_mapping_jacobian))

libmeanfield/impl/mapping/prepared_cache.cpp

@@ -0,0 +1,128 @@
+module;
+
+#include <algorithm>
+#include <cstddef>
+#include <mfem.hpp>
+#include <stdexcept>
+
+module mean_field;
+import :mapping.prepared_cache;
+import :mapping.types;
+
+namespace mean_field::mapping {
+    namespace {
+        void pack_vector(
+            double *&destination,
+            const mfem::Vector &vector,
+            const int dimension
+        ) {
+            if (vector.Size() != dimension)
+                throw std::invalid_argument("Prepared mapping vector dimension mismatch.");
+            std::copy_n(vector.HostRead(), dimension, destination);
+            destination += dimension;
+        }
+
+        void pack_matrix(
+            double *&destination,
+            const mfem::DenseMatrix &matrix,
+            const int dimension
+        ) {
+            if (matrix.Height() != dimension || matrix.Width() != dimension)
+                throw std::invalid_argument("Prepared mapping matrix dimension mismatch.");
+            std::copy_n(matrix.HostRead(), dimension * dimension, destination);
+            destination += dimension * dimension;
+        }
+
+        void unpack_vector(
+            const double *&source,
+            mfem::Vector &vector,
+            const int dimension
+        ) {
+            vector.SetSize(dimension);
+            std::copy_n(source, dimension, vector.HostWrite());
+            source += dimension;
+        }
+
+        void unpack_matrix(
+            const double *&source,
+            mfem::DenseMatrix &matrix,
+            const int dimension
+        ) {
+            matrix.SetSize(dimension);
+            std::copy_n(source, dimension * dimension, matrix.HostWrite());
+            source += dimension * dimension;
+        }
+    } // namespace
+
+    void VolumeMappingCache::SetSize(
+        const int point_count,
+        const int dimension
+    ) {
+        if (point_count < 0 || dimension < 1 || dimension > 3)
+            throw std::invalid_argument("Prepared mapping storage requires nonnegative point count and dimension 1-3.");
+        const int stride = 3 * dimension + 4 * dimension * dimension + 4;
+        m_data.resize(static_cast<std::size_t>(point_count) * stride);
+        m_point_count  = point_count;
+        m_dimension    = dimension;
+        m_point_stride = stride;
+    }
+
+    const double *VolumeMappingCache::GetPointData(const int point) const {
+        if (point < 0 || point >= m_point_count)
+            throw std::out_of_range("Prepared mapping quadrature point is out of range.");
+        return m_data.data() + static_cast<std::size_t>(point) * m_point_stride;
+    }
+
+    void VolumeMappingCache::Store(
+        const int point,
+        const VolumeMappingContext &context
+    ) {
+        // Validate the index through the same checked accessor used by readers.
+        (void)GetPointData(point);
+        double *data = m_data.data() + static_cast<std::size_t>(point) * m_point_stride;
+        pack_vector(data, context.mapping.reference_position, m_dimension);
+        pack_vector(data, context.mapping.displaced_position, m_dimension);
+        pack_vector(data, context.mapping.physical_position, m_dimension);
+        pack_matrix(data, context.mapping.displacement_jacobian, m_dimension);
+        pack_matrix(data, context.mapping.mapping_jacobian, m_dimension);
+        pack_matrix(data, context.mapping.inverse_mapping_jacobian, m_dimension);
+        pack_matrix(data, context.quadrature.J_inv, m_dimension);
+        *data++ = context.mapping.mapping_determinant;
+        *data++ = context.mapping.compactified ? 1.0 : 0.0;
+        *data++ = context.quadrature.detJ;
+        *data   = context.quadrature.weight;
+    }
+
+    void VolumeMappingCache::Load(
+        const int point,
+        VolumeMappingContext &context
+    ) const {
+        const double *data = GetPointData(point);
+        unpack_vector(data, context.mapping.reference_position, m_dimension);
+        unpack_vector(data, context.mapping.displaced_position, m_dimension);
+        unpack_vector(data, context.mapping.physical_position, m_dimension);
+        unpack_matrix(data, context.mapping.displacement_jacobian, m_dimension);
+        unpack_matrix(data, context.mapping.mapping_jacobian, m_dimension);
+        unpack_matrix(data, context.mapping.inverse_mapping_jacobian, m_dimension);
+        unpack_matrix(data, context.quadrature.J_inv, m_dimension);
+        context.mapping.mapping_determinant = *data++;
+        context.mapping.compactified        = *data++ != 0.0;
+        context.quadrature.detJ             = *data++;
+        context.quadrature.weight           = *data;
+    }
+
+    void VolumeMappingCache::LoadInverseJacobian(
+        const int point,
+        mfem::DenseMatrix &inverse
+    ) const {
+        const double *data = GetPointData(point) + 3 * m_dimension + 3 * m_dimension * m_dimension;
+        unpack_matrix(data, inverse, m_dimension);
+    }
+
+    int VolumeMappingCache::GetPointCount() const {
+        return m_point_count;
+    }
+    int VolumeMappingCache::GetDimension() const {
+        return m_dimension;
+    }
+} // namespace mean_field::mapping

libmeanfield/impl/operators/prepared_angular_momentum.cpp

@@ -311,13 +311,12 @@ namespace mean_field::operators {
             const mfem::FiniteElement &densityElement = *m_fem.densityFes->GetFE(elementId);
             const mfem::IntegrationRule &integrationRule =
                 get_moment_of_inertia_rule(m_fem, densityElement, *transformation);
-            data.quadraturePoints.resize(integrationRule.GetNPoints());
-            for (int quadraturePoint = 0; quadraturePoint < integrationRule.GetNPoints(); ++quadraturePoint) {
-                QuadraturePointData &point = data.quadraturePoints[quadraturePoint];
-                point.integrationPoint     = integrationRule.IntPoint(quadraturePoint);
-                point.densityShape.SetSize(densityElement.GetDof());
-                densityElement.CalcShape(point.integrationPoint, point.densityShape);
-            }
+            data.integrationRule = &integrationRule;
+            data.densityBasis    = m_fem.GetReferenceTables().GetScalarTable(densityElement, integrationRule);
+            data.mappingContexts.SetSize(integrationRule.GetNPoints(), m_fem.mesh->Dimension());
+            data.density.SetSize(integrationRule.GetNPoints());
+            data.quadratureWeights.SetSize(integrationRule.GetNPoints());
+            data.cylindricalRadiusSquared.SetSize(integrationRule.GetNPoints());
         }
         int globalStellarElementCount = 0;
         MFEM_VERIFY(
@@ -344,6 +343,7 @@ namespace mean_field::operators {
             return mapping::MappingStatus::non_finite_result;
         }
         mapping::DomainMapper::Workspace workspace(m_fem.mesh->Dimension());
+        mapping::VolumeMappingContext mappingContext;
 
         for (ElementPAData &data : m_elements) {
             displacementLocal.GetSubVector(data.displacementDofs, data.baseDisplacement);
@@ -372,21 +372,24 @@ namespace mean_field::operators {
                 .displacement = displacementData, .compactification = compactificationData
             };
             mfem::ElementTransformation *transformation = m_fem.mesh->GetElementTransformation(data.elementId);
-            for (QuadraturePointData &point : data.quadraturePoints) {
+            for (int quadraturePoint = 0; quadraturePoint < data.integrationRule->GetNPoints(); ++quadraturePoint) {
                 const mapping::MappingStatus status = m_domainMapper.EvaluateVolume(
-                    mappingData, *transformation, point.integrationPoint, workspace, point.mappingContext
+                    mappingData, *transformation, data.integrationRule->IntPoint(quadraturePoint), workspace,
+                    mappingContext
                 );
                 if (status != mapping::MappingStatus::valid) {
                     return status;
                 }
-                if (point.mappingContext.mapping.compactified) {
+                if (mappingContext.mapping.compactified) {
                     return mapping::MappingStatus::at_compactified_infinity;
                 }
-                point.cylindricalRadiusSquared =
-                    CylindricalRadiusSquared(point.mappingContext.mapping.physical_position);
-                if (!std::isfinite(point.cylindricalRadiusSquared)) {
+                data.cylindricalRadiusSquared(quadraturePoint) =
+                    CylindricalRadiusSquared(mappingContext.mapping.physical_position);
+                if (!std::isfinite(data.cylindricalRadiusSquared(quadraturePoint))) {
                     return mapping::MappingStatus::non_finite_result;
                 }
+                data.mappingContexts.Store(quadraturePoint, mappingContext);
+                data.quadratureWeights(quadraturePoint) = mappingContext.quadrature.weight;
             }
         }
         return std::nullopt;
@@ -411,11 +414,9 @@ namespace mean_field::operators {
             if (!is_finite_vector(elementDensity)) {
                 return false;
             }
-            for (QuadraturePointData &point : data.quadraturePoints) {
-                point.density = elementDensity * point.densityShape;
-                if (!std::isfinite(point.density)) {
-                    return false;
-                }
+            data.densityBasis->GetValues().Mult(elementDensity, data.density);
+            if (!is_finite_vector(data.density)) {
+                return false;
             }
         }
         return true;
@@ -424,9 +425,9 @@ namespace mean_field::operators {
     std::optional<AngularMomentumPreparationRejection> PreparedAngularMomentumOperator::TryAssembleResidual() {
         double localMomentOfInertia = 0.0;
         for (const ElementPAData &data : m_elements) {
-            for (const QuadraturePointData &point : data.quadraturePoints) {
-                localMomentOfInertia +=
-                    point.density * point.cylindricalRadiusSquared * point.mappingContext.quadrature.weight;
+            for (int quadraturePoint = 0; quadraturePoint < data.density.Size(); ++quadraturePoint) {
+                localMomentOfInertia += data.density(quadraturePoint) * data.cylindricalRadiusSquared(quadraturePoint) *
+                                        data.quadratureWeights(quadraturePoint);
             }
         }
         m_momentOfInertia = GlobalSum(localMomentOfInertia);
@@ -472,15 +473,18 @@ namespace mean_field::operators {
             "Angular-momentum density action has the wrong true-vector size."
         );
         true_to_local(*m_fem.densityFes, densityVariation, m_densityVariationLocal);
+        mfem::Vector quadratureDensityVariation;
         double localAction = 0.0;
         for (const ElementPAData &data : m_elements) {
             m_densityVariationLocal.GetSubVector(data.densityDofs, m_elementDensityVariation);
             if (data.densityDofTransformation != nullptr) {
                 data.densityDofTransformation->InvTransformPrimal(m_elementDensityVariation);
             }
-            for (const QuadraturePointData &point : data.quadraturePoints) {
-                localAction += (m_elementDensityVariation * point.densityShape) * point.cylindricalRadiusSquared *
-                               point.mappingContext.quadrature.weight;
+            quadratureDensityVariation.SetSize(data.integrationRule->GetNPoints());
+            data.densityBasis->GetValues().Mult(m_elementDensityVariation, quadratureDensityVariation);
+            for (int quadraturePoint = 0; quadraturePoint < quadratureDensityVariation.Size(); ++quadraturePoint) {
+                localAction += quadratureDensityVariation(quadraturePoint) *
+                               data.cylindricalRadiusSquared(quadraturePoint) * data.quadratureWeights(quadraturePoint);
             }
         }
         return localAction;
@@ -496,6 +500,7 @@ namespace mean_field::operators {
         true_to_local(*m_fem.displacementFes, displacementVariation, m_displacementVariationLocal);
         mapping::DomainMapper::Workspace workspace(m_fem.mesh->Dimension());
         mapping::VolumeMappingVariation variation;
+        mapping::VolumeMappingContext mappingContext;
         double localAction = 0.0;
         for (const ElementPAData &data : m_elements) {
             m_displacementVariationLocal.GetSubVector(data.displacementDofs, m_elementDisplacementVariation);
@@ -515,20 +520,22 @@ namespace mean_field::operators {
                 .displacement = baseDisplacementData, .compactification = compactificationData
             };
             mfem::ElementTransformation *transformation = m_fem.mesh->GetElementTransformation(data.elementId);
-            for (const QuadraturePointData &point : data.quadraturePoints) {
+            for (int quadraturePoint = 0; quadraturePoint < data.integrationRule->GetNPoints(); ++quadraturePoint) {
+                data.mappingContexts.Load(quadraturePoint, mappingContext);
                 const mapping::MappingStatus status = m_domainMapper.EvaluateVolumeVariation(
-                    mappingData, directionData, *transformation, point.integrationPoint, point.mappingContext,
-                    workspace, variation
+                    mappingData, directionData, *transformation, data.integrationRule->IntPoint(quadraturePoint),
+                    mappingContext, workspace, variation
                 );
                 MFEM_VERIFY(
                     status == mapping::MappingStatus::valid,
                     "Mapped angular-momentum variation is invalid. Element: " << data.elementId
                 );
                 const double radiusSquaredVariation = CylindricalRadiusSquaredVariation(
-                    point.mappingContext.mapping.physical_position, variation.mapping.physical_position_variation
+                    mappingContext.mapping.physical_position, variation.mapping.physical_position_variation
                 );
-                localAction += point.density * (radiusSquaredVariation * point.mappingContext.quadrature.weight +
-                                                point.cylindricalRadiusSquared * variation.weight_variation);
+                localAction += data.density(quadraturePoint) *
+                               (radiusSquaredVariation * data.quadratureWeights(quadraturePoint) +
+                                data.cylindricalRadiusSquared(quadraturePoint) * variation.weight_variation);
             }
         }
         return localAction;

libmeanfield/impl/operators/prepared_barotropic_closure.cpp

@@ -450,8 +450,8 @@ namespace mean_field::operators {
         m_displacementMap.scatter(m_context.GetDisplacement(), m_baseDisplacementTrue);
 
         m_isPrepared = false;
-        m_elements.clear();
         m_elements.reserve(m_fem.mesh->GetNE());
+        std::size_t preparedElementCount{0};
 
         mfem::Vector baseDensityLocal;
         mfem::Vector baseEnthalpyLocal;
@@ -462,6 +462,7 @@ namespace mean_field::operators {
         true_to_local(*m_fem.displacementFes, m_baseDisplacementTrue, displacementLocal);
 
         mapping::DomainMapper::Workspace workspace(m_fem.mesh->Dimension());
+        mapping::VolumeMappingContext mappingContext;
 
         mfem::Array<int> compactificationDofs;
 
@@ -485,8 +486,10 @@ namespace mean_field::operators {
                 continue;
             }
 
-            m_elements.emplace_back();
-            ElementPAData &data            = m_elements.back();
+            if (preparedElementCount == m_elements.size()) {
+                m_elements.emplace_back();
+            }
+            ElementPAData &data            = m_elements[preparedElementCount++];
             data.elementId                 = elementId;
 
             data.densityDofTransformation  = m_fem.densityFes->GetElementDofs(elementId, data.densityDofs);
@@ -533,13 +536,15 @@ namespace mean_field::operators {
 
             const mfem::IntegrationRule &integrationRule =
                 get_eos_rule(m_fem, m_equationOfState, densityElement, enthalpyElement, *transformation);
+            data.integrationRule           = &integrationRule;
 
             const int quadraturePointCount = integrationRule.GetNPoints();
             const int densityDofCount      = densityElement.GetDof();
             const int enthalpyDofCount     = enthalpyElement.GetDof();
 
-            data.densityBasis.SetSize(quadraturePointCount, densityDofCount);
-            data.enthalpyBasis.SetSize(quadraturePointCount, enthalpyDofCount);
+            data.densityBasis              = m_fem.GetReferenceTables().GetScalarTable(densityElement, integrationRule);
+            data.enthalpyBasis     = m_fem.GetReferenceTables().GetScalarTable(enthalpyElement, integrationRule);
+            data.displacementBasis = m_fem.GetReferenceTables().GetScalarTable(displacementElement, integrationRule);
             data.inverseElementJacobians.SetSize(
                 quadraturePointCount, m_fem.mesh->Dimension() * m_fem.mesh->Dimension()
             );
@@ -555,8 +560,6 @@ namespace mean_field::operators {
 
                 transformation->SetIntPoint(&integrationPoint);
 
-                mapping::VolumeMappingContext mappingContext;
-
                 const mapping::MappingStatus mappingStatus = m_domainMapper.EvaluateVolume(
                     mappingData, *transformation, integrationPoint, workspace, mappingContext
                 );
@@ -585,8 +588,8 @@ namespace mean_field::operators {
                     }
                 }
 
-                densityElement.CalcShape(integrationPoint, densityShape);
-                enthalpyElement.CalcShape(integrationPoint, enthalpyShape);
+                data.densityBasis->GetValues().GetRow(quadraturePoint, densityShape);
+                data.enthalpyBasis->GetValues().GetRow(quadraturePoint, enthalpyShape);
 
                 if (!vector_is_finite(densityShape) || !vector_is_finite(enthalpyShape)) {
                     retain_higher_priority_rejection(
@@ -595,13 +598,6 @@ namespace mean_field::operators {
                     continue;
                 }
 
-                for (int densityDof = 0; densityDof < densityDofCount; ++densityDof) {
-                    data.densityBasis(quadraturePoint, densityDof) = densityShape(densityDof);
-                }
-                for (int enthalpyDof = 0; enthalpyDof < enthalpyDofCount; ++enthalpyDof) {
-                    data.enthalpyBasis(quadraturePoint, enthalpyDof) = enthalpyShape(enthalpyDof);
-                }
-
                 const double density          = elementBaseDensity * densityShape;
                 const double enthalpy         = elementBaseEnthalpy * enthalpyShape;
                 const double quadratureWeight = mappingContext.quadrature.weight;
@@ -665,6 +661,7 @@ namespace mean_field::operators {
                 data.weightedEnthalpyDerivative(quadraturePoint) = weightedEnthalpyDerivative;
             }
         }
+        m_elements.resize(preparedElementCount);
 
         if (auto globalRejection = synchronize_rejection(localRejection, m_fem.densityFes->GetComm());
             globalRejection.has_value()) {
@@ -689,7 +686,7 @@ namespace mean_field::operators {
 
         for (const ElementPAData &data : m_elements) {
             elementResidual.SetSize(data.densityDofs.Size());
-            data.densityBasis.MultTranspose(data.weightedResidual, elementResidual);
+            data.densityBasis->GetValues().MultTranspose(data.weightedResidual, elementResidual);
 
             if (data.densityDofTransformation != nullptr) {
                 data.densityDofTransformation->TransformDual(elementResidual);
@@ -719,7 +716,7 @@ namespace mean_field::operators {
             elementDiagonal = 0.0;
             for (int trialDof = 0; trialDof < data.densityDofs.Size(); ++trialDof) {
                 for (int quadraturePoint = 0; quadraturePoint < data.quadratureWeights.Size(); ++quadraturePoint) {
-                    const double basis = data.densityBasis(quadraturePoint, trialDof);
+                    const double basis = data.densityBasis->GetValues()(quadraturePoint, trialDof);
                     elementDiagonal(trialDof) += data.quadratureWeights(quadraturePoint) * basis * basis;
                 }
             }
@@ -842,8 +839,8 @@ namespace mean_field::operators {
             quadratureEnthalpyVariation.SetSize(data.quadratureWeights.Size());
             quadratureAction.SetSize(data.quadratureWeights.Size());
 
-            data.densityBasis.Mult(elementDensityVariation, quadratureDensityVariation);
-            data.enthalpyBasis.Mult(elementEnthalpyVariation, quadratureEnthalpyVariation);
+            data.densityBasis->GetValues().Mult(elementDensityVariation, quadratureDensityVariation);
+            data.enthalpyBasis->GetValues().Mult(elementEnthalpyVariation, quadratureEnthalpyVariation);
 
             for (int quadraturePoint = 0; quadraturePoint < quadratureAction.Size(); ++quadraturePoint) {
                 quadratureAction(quadraturePoint) =
@@ -852,7 +849,7 @@ namespace mean_field::operators {
             }
 
             elementAction.SetSize(data.densityDofs.Size());
-            data.densityBasis.MultTranspose(quadratureAction, elementAction);
+            data.densityBasis->GetValues().MultTranspose(quadratureAction, elementAction);
 
             if (data.densityDofTransformation != nullptr) {
                 data.densityDofTransformation->TransformDual(elementAction);
@@ -906,14 +903,14 @@ namespace mean_field::operators {
                 "Prepared barotropic closure inverse-Jacobian data has an incompatible size."
             );
 
-            m_referenceDShape.SetSize(displacementElement.GetDof(), dimension);
             m_referenceDisplacementJacobian.SetSize(dimension, dimension);
             m_quadratureDisplacementAction.SetSize(integrationRule.GetNPoints());
 
             for (int quadraturePoint = 0; quadraturePoint < integrationRule.GetNPoints(); ++quadraturePoint) {
-                const mfem::IntegrationPoint &integrationPoint = integrationRule.IntPoint(quadraturePoint);
-                displacementElement.CalcDShape(integrationPoint, m_referenceDShape);
-                mfem::MultAtB(directionDofs, m_referenceDShape, m_referenceDisplacementJacobian);
+                mfem::MultAtB(
+                    directionDofs, data.displacementBasis->GetGradients(quadraturePoint),
+                    m_referenceDisplacementJacobian
+                );
 
                 double logarithmicJacobianVariation{0.0};
                 for (int row = 0; row < dimension; ++row) {
@@ -933,7 +930,7 @@ namespace mean_field::operators {
             }
 
             m_elementDisplacementAction.SetSize(data.densityDofs.Size());
-            data.densityBasis.MultTranspose(m_quadratureDisplacementAction, m_elementDisplacementAction);
+            data.densityBasis->GetValues().MultTranspose(m_quadratureDisplacementAction, m_elementDisplacementAction);
 
             if (data.densityDofTransformation != nullptr) {
                 data.densityDofTransformation->TransformDual(m_elementDisplacementAction);

libmeanfield/impl/operators/prepared_gravity_displacement_force.cpp

@@ -13,6 +13,7 @@ module mean_field;
 
 import :operators.kernels.gravity_displacement_force;
 import :operators.prepared_gravity_displacement_force;
+import :fem.reference_tables;
 
 namespace {
     using DomainSchema = mean_field::utils::domain::CoreEnvelopeVacuumDomainSchema;
@@ -249,6 +250,17 @@ namespace mean_field::operators {
             const mfem::FiniteElement &displacementElement     = *m_fem.displacementFes->GetFE(elementId);
             const mfem::FiniteElement &compactificationElement = *m_fem.compactificationFes->GetFE(elementId);
             data.integrationRule                               = &get_gravity_force_rule(m_fem, *transformation);
+            data.densityReferenceTable =
+                m_fem.GetReferenceTables().GetScalarTable(densityElement, *data.integrationRule);
+            data.displacementReferenceTable =
+                m_fem.GetReferenceTables().GetScalarTable(displacementElement, *data.integrationRule);
+            if (gravityGradientElement.GetMapType() == mfem::FiniteElement::H_DIV &&
+                gravityGradientElement.GetDim() == dimension && gravityGradientElement.GetRangeDim() == dimension &&
+                transformation->GetSpaceDim() == dimension) {
+                data.gravityReferenceTable =
+                    m_fem.GetReferenceTables().GetVectorTable(gravityGradientElement, *data.integrationRule);
+                data.meshPiolaJacobians.SetSize(data.integrationRule->GetNPoints(), dimension * dimension);
+            }
 
             const mapping::ElementDisplacementData displacementData =
                 mapping::ElementDisplacementDataFromElementVDofs(displacementElement, elementBaseDisplacement);
@@ -282,13 +294,17 @@ namespace mean_field::operators {
                     "Prepared gravity force encountered compactification on a stellar element."
                 );
 
-                densityElement.CalcShape(integrationPoint, densityShape);
+                for (int dof = 0; dof < densityElement.GetDof(); ++dof) {
+                    densityShape(dof) = data.densityReferenceTable->GetValues()(quadraturePoint, dof);
+                }
                 gravityGradientElement.CalcVShape(*transformation, gravityGradientShape);
                 gravityGradientShape.MultTranspose(elementBaseGravityGradient, baseGravityReferenceValue);
                 data.baseDensityValues(quadraturePoint)      = elementBaseDensity * densityShape;
                 data.referenceWeights(quadraturePoint)       = integrationPoint.weight * transformation->Weight();
 
                 const mfem::DenseMatrix &inverseMeshJacobian = transformation->InverseJacobian();
+                const mfem::DenseMatrix &meshJacobian        = transformation->Jacobian();
+                const double inverseMeshWeight               = 1.0 / transformation->Weight();
                 for (int row = 0; row < dimension; ++row) {
                     data.baseGravityReferenceValues(quadraturePoint, row) = baseGravityReferenceValue(row);
                     for (int column = 0; column < dimension; ++column) {
@@ -296,6 +312,10 @@ namespace mean_field::operators {
                         data.mappingJacobians(quadraturePoint, entry) =
                             mappingContext.mapping.mapping_jacobian(row, column);
                         data.inverseMeshJacobians(quadraturePoint, entry) = inverseMeshJacobian(row, column);
+                        if (data.gravityReferenceTable != nullptr) {
+                            data.meshPiolaJacobians(quadraturePoint, entry) =
+                                inverseMeshWeight * meshJacobian(row, column);
+                        }
                     }
                 }
 
@@ -308,7 +328,9 @@ namespace mean_field::operators {
                     for (int column = 0; column < dimension; ++column) {
                         const int entry = row * dimension + column;
                         if (!std::isfinite(data.mappingJacobians(quadraturePoint, entry)) ||
-                            !std::isfinite(data.inverseMeshJacobians(quadraturePoint, entry))) {
+                            !std::isfinite(data.inverseMeshJacobians(quadraturePoint, entry)) ||
+                            (data.gravityReferenceTable != nullptr &&
+                             !std::isfinite(data.meshPiolaJacobians(quadraturePoint, entry)))) {
                             return std::unexpected(non_finite_rejection());
                         }
                     }
@@ -462,6 +484,10 @@ namespace mean_field::operators {
 
         for (const ElementPAData &data : m_elements) {
             MFEM_VERIFY(data.integrationRule != nullptr, "Prepared gravity force has no integration rule.");
+            MFEM_VERIFY(
+                data.densityReferenceTable != nullptr && data.displacementReferenceTable != nullptr,
+                "Prepared gravity force has no reference basis tables."
+            );
 
             m_densityVariationLocal.GetSubVector(data.densityDofs, m_elementDensityVariation);
             m_gravityGradientVariationLocal.GetSubVector(data.gravityGradientDofs, m_elementGravityGradientVariation);
@@ -490,13 +516,14 @@ namespace mean_field::operators {
             m_densityShape.SetSize(densityElement.GetDof());
             m_displacementShape.SetSize(scalarDisplacementDofCount);
             m_gravityGradientShape.SetSize(gravityGradientElement.GetDof(), dimension);
-            m_referenceDisplacementDShape.SetSize(scalarDisplacementDofCount, dimension);
             m_referenceDisplacementJacobian.SetSize(dimension, dimension);
             m_displacementJacobianVariation.SetSize(dimension, dimension);
             m_mappingJacobian.SetSize(dimension, dimension);
             m_inverseMeshJacobian.SetSize(dimension, dimension);
+            m_meshPiolaJacobian.SetSize(dimension, dimension);
             m_baseGravityReferenceValue.SetSize(dimension);
             m_gravityVariationReferenceValue.SetSize(dimension);
+            m_gravityVariationReferenceCellValue.SetSize(dimension);
             m_mappedBaseGravity.SetSize(dimension);
             m_mappedGravityVariation.SetSize(dimension);
             m_mappedGeometryVariation.SetSize(dimension);
@@ -506,16 +533,28 @@ namespace mean_field::operators {
 
             for (int quadraturePoint = 0; quadraturePoint < data.integrationRule->GetNPoints(); ++quadraturePoint) {
                 const mfem::IntegrationPoint &integrationPoint = data.integrationRule->IntPoint(quadraturePoint);
-                densityElement.CalcShape(integrationPoint, m_densityShape);
-                displacementElement.CalcShape(integrationPoint, m_displacementShape);
-                displacementElement.CalcDShape(integrationPoint, m_referenceDisplacementDShape);
-                mfem::MultAtB(directionDofs, m_referenceDisplacementDShape, m_referenceDisplacementJacobian);
+                const mfem::DenseMatrix &referenceDisplacementDShape =
+                    data.displacementReferenceTable->GetGradients(quadraturePoint);
+                mfem::MultAtB(directionDofs, referenceDisplacementDShape, m_referenceDisplacementJacobian);
+                const mfem::DenseMatrix &densityValues      = data.densityReferenceTable->GetValues();
+                const mfem::DenseMatrix &displacementValues = data.displacementReferenceTable->GetValues();
+                for (int dof = 0; dof < densityElement.GetDof(); ++dof) {
+                    m_densityShape(dof) = densityValues(quadraturePoint, dof);
+                }
+                for (int dof = 0; dof < scalarDisplacementDofCount; ++dof) {
+                    m_displacementShape(dof) = displacementValues(quadraturePoint, dof);
+                }
 
-                transformation->SetIntPoint(&integrationPoint);
-                gravityGradientElement.CalcVShape(*transformation, m_gravityGradientShape);
-                m_gravityGradientShape.MultTranspose(
-                    m_elementGravityGradientVariation, m_gravityVariationReferenceValue
-                );
+                if (data.gravityReferenceTable != nullptr) {
+                    data.gravityReferenceTable->GetValues(quadraturePoint)
+                        .MultTranspose(m_elementGravityGradientVariation, m_gravityVariationReferenceCellValue);
+                } else {
+                    transformation->SetIntPoint(&integrationPoint);
+                    gravityGradientElement.CalcVShape(*transformation, m_gravityGradientShape);
+                    m_gravityGradientShape.MultTranspose(
+                        m_elementGravityGradientVariation, m_gravityVariationReferenceValue
+                    );
+                }
 
                 for (int row = 0; row < dimension; ++row) {
                     m_baseGravityReferenceValue(row) = data.baseGravityReferenceValues(quadraturePoint, row);
@@ -523,8 +562,14 @@ namespace mean_field::operators {
                         const int entry                    = row * dimension + column;
                         m_mappingJacobian(row, column)     = data.mappingJacobians(quadraturePoint, entry);
                         m_inverseMeshJacobian(row, column) = data.inverseMeshJacobians(quadraturePoint, entry);
+                        if (data.gravityReferenceTable != nullptr) {
+                            m_meshPiolaJacobian(row, column) = data.meshPiolaJacobians(quadraturePoint, entry);
+                        }
                     }
                 }
+                if (data.gravityReferenceTable != nullptr) {
+                    m_meshPiolaJacobian.Mult(m_gravityVariationReferenceCellValue, m_gravityVariationReferenceValue);
+                }
                 mfem::Mult(m_referenceDisplacementJacobian, m_inverseMeshJacobian, m_displacementJacobianVariation);
                 m_mappingJacobian.Mult(m_baseGravityReferenceValue, m_mappedBaseGravity);
                 m_mappingJacobian.Mult(m_gravityVariationReferenceValue, m_mappedGravityVariation);

libmeanfield/impl/operators/prepared_gravity_source.cpp

@@ -182,19 +182,17 @@ namespace {
                 .displacement = *m_displacement_data, .compactification = *m_compactification_data
             };
 
-            mean_field::mapping::VolumeMappingContext mapping_context;
-
             const mean_field::mapping::MappingStatus status = m_domain_mapper.EvaluateVolume(
-                mapping_data, transformation, integration_point, m_workspace, mapping_context
+                mapping_data, transformation, integration_point, m_workspace, m_mapping_context
             );
 
             if (status != mean_field::mapping::MappingStatus::valid) {
                 m_mappingFailure = status;
                 return 0.0;
             }
-            const double mapping_determinant = mapping_context.mapping.mapping_determinant;
+            const double mapping_determinant = m_mapping_context.mapping.mapping_determinant;
 
-            m_inverse_element_jacobian       = mapping_context.quadrature.J_inv;
+            m_inverse_element_jacobian       = m_mapping_context.quadrature.J_inv;
 
             const double value               = 4.0 * std::numbers::pi * mean_field::utils::G * mapping_determinant;
             if (!std::isfinite(value)) {
@@ -269,6 +267,7 @@ namespace {
         std::unique_ptr<mean_field::mapping::ElementCompactificationData> m_compactification_data;
 
         mean_field::mapping::DomainMapper::Workspace m_workspace;
+        mean_field::mapping::VolumeMappingContext m_mapping_context;
         mfem::DenseMatrix m_inverse_element_jacobian;
         int m_cached_element_id{-1};
         mean_field::mapping::MappingStatus m_mappingFailure{mean_field::mapping::MappingStatus::valid};
@@ -394,8 +393,8 @@ namespace mean_field::operators {
         m_has_variation_data = false;
         m_displacement_true.SetSize(m_displacement_map.full_size());
         m_displacement_map.scatter(displacement, m_displacement_true);
-        m_elements.clear();
         m_elements.reserve(m_fem.mesh->GetNE());
+        std::size_t prepared_element_count{0};
 
         FrozenMappedGravitySourceCoefficient source_coefficient(m_fem, m_domain_mapper, m_displacement_true);
         bool localNonFiniteQuadrature = false;
@@ -407,8 +406,10 @@ namespace mean_field::operators {
                 continue;
             }
 
-            m_elements.emplace_back();
-            ElementPAData &data             = m_elements.back();
+            if (prepared_element_count == m_elements.size()) {
+                m_elements.emplace_back();
+            }
+            ElementPAData &data             = m_elements[prepared_element_count++];
 
             data.element_id                 = element_id;
 
@@ -438,37 +439,59 @@ namespace mean_field::operators {
 
             const int potential_dof_count    = potential_element.GetDof();
 
-            data.density_basis.SetSize(quadrature_point_count, density_dof_count);
-
-            data.potential_basis.SetSize(quadrature_point_count, potential_dof_count);
+            if (density_element.GetMapType() == mfem::FiniteElement::VALUE) {
+                data.density_reference = m_fem.GetReferenceTables().GetScalarTable(density_element, integration_rule);
+                data.density_basis.SetSize(0, 0);
+            } else {
+                data.density_reference.reset();
+                data.density_basis.SetSize(quadrature_point_count, density_dof_count);
+            }
+            if (potential_element.GetMapType() == mfem::FiniteElement::VALUE) {
+                data.potential_reference =
+                    m_fem.GetReferenceTables().GetScalarTable(potential_element, integration_rule);
+                data.potential_basis.SetSize(0, 0);
+            } else {
+                data.potential_reference.reset();
+                data.potential_basis.SetSize(quadrature_point_count, potential_dof_count);
+            }
 
             const int dimension = m_fem.mesh->Dimension();
             if (mode == PreparationMode::linearization) {
                 data.inverse_element_jacobians.SetSize(quadrature_point_count, dimension * dimension);
+                data.displacement_reference = m_fem.GetReferenceTables().GetScalarTable(
+                    *m_fem.displacementFes->GetFE(element_id), integration_rule
+                );
             }
 
             data.quadrature_data.SetSize(quadrature_point_count);
 
-            mfem::Vector density_shape(density_dof_count);
-            mfem::Vector potential_shape(potential_dof_count);
+            mfem::Vector density_shape;
+            mfem::Vector potential_shape;
+            if (!data.density_reference) {
+                density_shape.SetSize(density_dof_count);
+            }
+            if (!data.potential_reference) {
+                potential_shape.SetSize(potential_dof_count);
+            }
 
             for (int quadrature_point = 0; quadrature_point < quadrature_point_count; ++quadrature_point) {
                 const mfem::IntegrationPoint &integration_point = integration_rule.IntPoint(quadrature_point);
 
                 transformation.SetIntPoint(&integration_point);
 
-                // CalcPhysShape matches the scalar mixed-mass discretization,
-                // including the finite-element map type.
-                density_element.CalcPhysShape(transformation, density_shape);
-
-                potential_element.CalcPhysShape(transformation, potential_shape);
-
-                for (int i = 0; i < density_dof_count; ++i) {
-                    data.density_basis(quadrature_point, i) = density_shape(i);
+                // VALUE maps use the shared reference basis. Preserve the
+                // physical-shape evaluation for every other scalar map type.
+                if (!data.density_reference) {
+                    density_element.CalcPhysShape(transformation, density_shape);
+                    for (int i = 0; i < density_dof_count; ++i) {
+                        data.density_basis(quadrature_point, i) = density_shape(i);
+                    }
                 }
-
-                for (int i = 0; i < potential_dof_count; ++i) {
-                    data.potential_basis(quadrature_point, i) = potential_shape(i);
+                if (!data.potential_reference) {
+                    potential_element.CalcPhysShape(transformation, potential_shape);
+                    for (int i = 0; i < potential_dof_count; ++i) {
+                        data.potential_basis(quadrature_point, i) = potential_shape(i);
+                    }
                 }
 
                 const double coefficient_value = source_coefficient.Eval(transformation, integration_point);
@@ -505,6 +528,7 @@ namespace mean_field::operators {
                 break;
             }
         }
+        m_elements.resize(prepared_element_count);
 
         const bool localNonFiniteArithmetic = source_coefficient.HasNonFiniteArithmetic() || localNonFiniteQuadrature;
         auto preparationResult              = synchronize_preparation_failure(
@@ -555,7 +579,7 @@ namespace mean_field::operators {
             m_quadrature_action.SetSize(data.quadrature_data.Size());
 
             // B_density * x_e
-            data.density_basis.Mult(m_element_input, m_quadrature_action);
+            data.GetDensityBasis().Mult(m_element_input, m_quadrature_action);
 
             // D * B_density * x_e
             for (int q = 0; q < m_quadrature_action.Size(); ++q) {
@@ -565,7 +589,7 @@ namespace mean_field::operators {
             m_element_action.SetSize(data.potential_dofs.Size());
 
             // B_potential^T * D * B_density * x_e
-            data.potential_basis.MultTranspose(m_quadrature_action, m_element_action);
+            data.GetPotentialBasis().MultTranspose(m_quadrature_action, m_element_action);
 
             if (data.potential_dof_transformation != nullptr) {
                 data.potential_dof_transformation->TransformDual(m_element_action);
@@ -641,15 +665,15 @@ namespace mean_field::operators {
                 "Prepared gravity source inverse-Jacobian data has an incompatible size."
             );
 
-            m_reference_displacement_dshape.SetSize(displacement_element.GetDof(), dimension);
             m_reference_displacement_jacobian.SetSize(dimension, dimension);
             m_quadrature_variation_action.SetSize(data.integration_rule->GetNPoints());
-            data.density_basis.Mult(m_element_density, m_quadrature_variation_action);
+            data.GetDensityBasis().Mult(m_element_density, m_quadrature_variation_action);
 
             for (int quadrature_point = 0; quadrature_point < data.integration_rule->GetNPoints(); ++quadrature_point) {
-                const mfem::IntegrationPoint &integration_point = data.integration_rule->IntPoint(quadrature_point);
-                displacement_element.CalcDShape(integration_point, m_reference_displacement_dshape);
-                mfem::MultAtB(direction_dofs, m_reference_displacement_dshape, m_reference_displacement_jacobian);
+                mfem::MultAtB(
+                    direction_dofs, data.displacement_reference->GetGradients(quadrature_point),
+                    m_reference_displacement_jacobian
+                );
 
                 double logarithmic_jacobian_variation{0.0};
                 for (int row = 0; row < dimension; ++row) {
@@ -669,7 +693,7 @@ namespace mean_field::operators {
             }
 
             m_element_variation_action.SetSize(data.potential_dofs.Size());
-            data.potential_basis.MultTranspose(m_quadrature_variation_action, m_element_variation_action);
+            data.GetPotentialBasis().MultTranspose(m_quadrature_variation_action, m_element_variation_action);
 
             if (data.potential_dof_transformation != nullptr) {
                 data.potential_dof_transformation->TransformDual(m_element_variation_action);
@@ -714,7 +738,7 @@ namespace mean_field::operators {
 
             m_quadrature_action.SetSize(data.quadrature_data.Size());
 
-            data.potential_basis.Mult(m_element_input, m_quadrature_action);
+            data.GetPotentialBasis().Mult(m_element_input, m_quadrature_action);
 
             for (int q = 0; q < m_quadrature_action.Size(); ++q) {
                 m_quadrature_action(q) *= data.quadrature_data(q);
@@ -722,7 +746,7 @@ namespace mean_field::operators {
 
             m_element_action.SetSize(data.density_dofs.Size());
 
-            data.density_basis.MultTranspose(m_quadrature_action, m_element_action);
+            data.GetDensityBasis().MultTranspose(m_quadrature_action, m_element_action);
 
             if (data.density_dof_transformation != nullptr) {
                 data.density_dof_transformation->TransformDual(m_element_action);

libmeanfield/impl/operators/prepared_hdiv_mass.cpp

@@ -13,6 +13,8 @@ module;
 #include <mpi.h>
 
 module mean_field;
+
+import :fem.reference_tables;
 import :operators.prepared_hdiv_mass;
 
 namespace {
@@ -551,6 +553,15 @@ namespace mean_field::operators {
             );
 
             data.integrationRule = &get_hdiv_mass_rule(m_fem, m_domain_mapper, gravityGradientElement, *transformation);
+            const int dimension  = m_domain_mapper.GetDimension();
+            if (gravityGradientElement.GetMapType() == mfem::FiniteElement::H_DIV &&
+                gravityGradientElement.GetDim() == dimension && gravityGradientElement.GetRangeDim() == dimension &&
+                transformation->GetSpaceDim() == dimension) {
+                data.gravityReferenceTable =
+                    m_fem.GetReferenceTables().GetVectorTable(gravityGradientElement, *data.integrationRule);
+                data.meshPiolaJacobians.SetSize(data.integrationRule->GetNPoints(), dimension * dimension);
+                data.referenceWeights.SetSize(data.integrationRule->GetNPoints());
+            }
             data.frozenMappingData.SetSize(
                 data.integrationRule->GetNPoints(), frozen_mapping_width(m_domain_mapper.GetDimension())
             );
@@ -573,6 +584,25 @@ namespace mean_field::operators {
                     return status;
                 }
                 freeze_mapping_context(mappingContext, quadraturePoint, data.frozenMappingData);
+                if (data.gravityReferenceTable != nullptr) {
+                    // CalcVShape_RT = reference_shape * J_mesh^T / Weight.
+                    // Cache only this small factor, never the mapped basis.
+                    const double meshWeight                = transformation->Weight();
+                    const mfem::DenseMatrix &meshJacobian  = transformation->Jacobian();
+                    const double inverseMeshWeight         = 1.0 / meshWeight;
+                    data.referenceWeights(quadraturePoint) = integrationPoint.weight * meshWeight;
+                    for (int row = 0; row < dimension; ++row) {
+                        for (int column = 0; column < dimension; ++column) {
+                            const double entry = inverseMeshWeight * meshJacobian(row, column);
+                            if (!std::isfinite(entry))
+                                return mapping::MappingStatus::non_finite_result;
+                            data.meshPiolaJacobians(quadraturePoint, row * dimension + column) = entry;
+                        }
+                    }
+                    if (!std::isfinite(data.referenceWeights(quadraturePoint))) {
+                        return mapping::MappingStatus::non_finite_result;
+                    }
+                }
             }
         }
         return mapping::MappingStatus::valid;
@@ -830,7 +860,10 @@ namespace mean_field::operators {
             m_elementVariationAction.SetSize(gravityGradientElement.GetDof());
             m_elementVariationAction = 0.0;
             m_gravityGradientValue.SetSize(dimension);
+            m_gravityReferenceCellValue.SetSize(dimension);
+            m_referenceCellDual.SetSize(dimension);
             m_massTensorVariationAction.SetSize(dimension);
+            m_meshPiolaJacobian.SetSize(dimension, dimension);
             m_gravityGradientShape.SetSize(gravityGradientElement.GetDof(), dimension);
             m_massTensorVariation.SetSize(dimension, dimension);
 
@@ -853,12 +886,33 @@ namespace mean_field::operators {
                     m_baseMappingContext.mapping, m_mappingVariation.mapping, m_massTensorVariation
                 );
 
-                transformation->SetIntPoint(&integrationPoint);
-                gravityGradientElement.CalcVShape(*transformation, m_gravityGradientShape);
-                m_gravityGradientShape.MultTranspose(m_elementGravityGradient, m_gravityGradientValue);
-                m_massTensorVariation.Mult(m_gravityGradientValue, m_massTensorVariationAction);
-                const double referenceWeight = integrationPoint.weight * transformation->Weight();
-                m_gravityGradientShape.AddMult(m_massTensorVariationAction, m_elementVariationAction, referenceWeight);
+                if (data.gravityReferenceTable != nullptr) {
+                    const mfem::DenseMatrix &referenceShape = data.gravityReferenceTable->GetValues(quadraturePoint);
+                    referenceShape.MultTranspose(m_elementGravityGradient, m_gravityReferenceCellValue);
+                    for (int row = 0; row < dimension; ++row) {
+                        for (int column = 0; column < dimension; ++column) {
+                            m_meshPiolaJacobian(row, column) =
+                                data.meshPiolaJacobians(quadraturePoint, row * dimension + column);
+                        }
+                    }
+                    m_meshPiolaJacobian.Mult(m_gravityReferenceCellValue, m_gravityGradientValue);
+                    m_massTensorVariation.Mult(m_gravityGradientValue, m_massTensorVariationAction);
+                    // Move the test-side Piola transform onto the three-vector
+                    // dual before applying the reference basis transpose.
+                    m_meshPiolaJacobian.MultTranspose(m_massTensorVariationAction, m_referenceCellDual);
+                    referenceShape.AddMult(
+                        m_referenceCellDual, m_elementVariationAction, data.referenceWeights(quadraturePoint)
+                    );
+                } else {
+                    transformation->SetIntPoint(&integrationPoint);
+                    gravityGradientElement.CalcVShape(*transformation, m_gravityGradientShape);
+                    m_gravityGradientShape.MultTranspose(m_elementGravityGradient, m_gravityGradientValue);
+                    m_massTensorVariation.Mult(m_gravityGradientValue, m_massTensorVariationAction);
+                    const double referenceWeight = integrationPoint.weight * transformation->Weight();
+                    m_gravityGradientShape.AddMult(
+                        m_massTensorVariationAction, m_elementVariationAction, referenceWeight
+                    );
+                }
             }
 
             if (data.gravityGradientDofTransformation != nullptr) {

libmeanfield/impl/operators/prepared_hydrostatic_equilibrium.cpp

@@ -488,9 +488,6 @@ namespace mean_field::operators {
         m_elements.clear();
         m_elements.reserve(m_fem.mesh->GetNE());
 
-        mfem::Vector enthalpyShape;
-        mfem::Vector gravityPotentialShape;
-
         for (int elementId = 0; elementId < m_fem.mesh->GetNE(); ++elementId) {
             mfem::ElementTransformation *transformation = m_fem.mesh->GetElementTransformation(elementId);
 
@@ -528,34 +525,10 @@ namespace mean_field::operators {
             data.integrationRule =
                 &get_hydrostatic_rule(m_fem, enthalpyElement, gravityPotentialElement, *transformation);
 
-            const int quadraturePointCount     = data.integrationRule->GetNPoints();
-
-            const int enthalpyDofCount         = enthalpyElement.GetDof();
-
-            const int gravityPotentialDofCount = gravityPotentialElement.GetDof();
-
-            data.enthalpyBasis.SetSize(quadraturePointCount, enthalpyDofCount);
-
-            data.gravityPotentialBasis.SetSize(quadraturePointCount, gravityPotentialDofCount);
-
-            enthalpyShape.SetSize(enthalpyDofCount);
-            gravityPotentialShape.SetSize(gravityPotentialDofCount);
-
-            for (int quadraturePoint = 0; quadraturePoint < quadraturePointCount; ++quadraturePoint) {
-                const mfem::IntegrationPoint &integrationPoint = data.integrationRule->IntPoint(quadraturePoint);
-
-                enthalpyElement.CalcShape(integrationPoint, enthalpyShape);
-
-                gravityPotentialElement.CalcShape(integrationPoint, gravityPotentialShape);
-
-                for (int dof = 0; dof < enthalpyDofCount; ++dof) {
-                    data.enthalpyBasis(quadraturePoint, dof) = enthalpyShape(dof);
-                }
-
-                for (int dof = 0; dof < gravityPotentialDofCount; ++dof) {
-                    data.gravityPotentialBasis(quadraturePoint, dof) = gravityPotentialShape(dof);
-                }
-            }
+            const fem::ReferenceTableCache &referenceTables = m_fem.GetReferenceTables();
+            data.enthalpyReferenceTable = referenceTables.GetScalarTable(enthalpyElement, *data.integrationRule);
+            data.gravityPotentialReferenceTable =
+                referenceTables.GetScalarTable(gravityPotentialElement, *data.integrationRule);
         }
     }
 
@@ -565,6 +538,7 @@ namespace mean_field::operators {
         true_to_local(*m_fem.displacementFes, m_context.GetDisplacementTrue(), displacementLocal);
 
         mapping::DomainMapper::Workspace workspace(m_fem.mesh->Dimension());
+        mapping::VolumeMappingContext mappingContext;
 
         mfem::Array<int> compactificationDofs;
 
@@ -615,16 +589,14 @@ namespace mean_field::operators {
 
             data.quadratureWeights.SetSize(quadraturePointCount);
 
-            data.baseMappingContexts.resize(quadraturePointCount);
+            data.baseMappingContexts.SetSize(quadraturePointCount, m_fem.mesh->Dimension());
 
             for (int quadraturePoint = 0; quadraturePoint < quadraturePointCount; ++quadraturePoint) {
                 const mfem::IntegrationPoint &integrationPoint = data.integrationRule->IntPoint(quadraturePoint);
 
                 transformation->SetIntPoint(&integrationPoint);
 
-                mapping::VolumeMappingContext &mappingContext = data.baseMappingContexts[quadraturePoint];
-
-                const mapping::MappingStatus mappingStatus    = m_domainMapper.EvaluateVolume(
+                const mapping::MappingStatus mappingStatus = m_domainMapper.EvaluateVolume(
                     mappingData, *transformation, integrationPoint, workspace, mappingContext
                 );
 
@@ -641,6 +613,7 @@ namespace mean_field::operators {
                     return mapping::MappingStatus::non_positive_determinant;
                 }
 
+                data.baseMappingContexts.Store(quadraturePoint, mappingContext);
                 data.quadratureWeights(quadraturePoint) = quadratureWeight;
 
                 for (int component = 0; component < m_fem.mesh->Dimension(); ++component) {
@@ -660,15 +633,17 @@ namespace mean_field::operators {
 
     bool PreparedHydrostaticEquilibriumOperator::PrepareAlgebraicJacobianBlocks() {
         for (ElementPAData &data : m_elements) {
-            const int quadraturePointCount     = data.quadratureWeights.Size();
+            const mfem::DenseMatrix &enthalpyBasis         = data.GetEnthalpyBasis();
+            const mfem::DenseMatrix &gravityPotentialBasis = data.GetGravityPotentialBasis();
+            const int quadraturePointCount                 = data.quadratureWeights.Size();
 
-            const int enthalpyDofCount         = data.enthalpyBasis.Width();
+            const int enthalpyDofCount                     = enthalpyBasis.Width();
 
-            const int gravityPotentialDofCount = data.gravityPotentialBasis.Width();
+            const int gravityPotentialDofCount             = gravityPotentialBasis.Width();
 
             MFEM_VERIFY(
-                data.enthalpyBasis.Height() == quadraturePointCount &&
-                    data.gravityPotentialBasis.Height() == quadraturePointCount,
+                enthalpyBasis.Height() == quadraturePointCount &&
+                    gravityPotentialBasis.Height() == quadraturePointCount,
                 "Prepared hydrostatic algebraic Jacobian has "
                 "inconsistent quadrature data."
             );
@@ -687,18 +662,18 @@ namespace mean_field::operators {
                 const double quadratureWeight = data.quadratureWeights(quadraturePoint);
 
                 for (int testDof = 0; testDof < enthalpyDofCount; ++testDof) {
-                    const double weightedTestBasis = quadratureWeight * data.enthalpyBasis(quadraturePoint, testDof);
+                    const double weightedTestBasis = quadratureWeight * enthalpyBasis(quadraturePoint, testDof);
 
                     data.bernoulliConstantJacobian(testDof) -= weightedTestBasis;
 
                     for (int trialDof = 0; trialDof < enthalpyDofCount; ++trialDof) {
                         data.enthalpyJacobian(testDof, trialDof) +=
-                            weightedTestBasis * data.enthalpyBasis(quadraturePoint, trialDof);
+                            weightedTestBasis * enthalpyBasis(quadraturePoint, trialDof);
                     }
 
                     for (int trialDof = 0; trialDof < gravityPotentialDofCount; ++trialDof) {
                         data.gravityPotentialJacobian(testDof, trialDof) +=
-                            weightedTestBasis * data.gravityPotentialBasis(quadraturePoint, trialDof);
+                            weightedTestBasis * gravityPotentialBasis(quadraturePoint, trialDof);
                     }
                 }
             }
@@ -795,9 +770,9 @@ namespace mean_field::operators {
 
             quadratureGravityPotential.SetSize(quadraturePointCount);
 
-            data.enthalpyBasis.Mult(elementEnthalpy, quadratureEnthalpy);
+            data.GetEnthalpyBasis().Mult(elementEnthalpy, quadratureEnthalpy);
 
-            data.gravityPotentialBasis.Mult(elementGravityPotential, quadratureGravityPotential);
+            data.GetGravityPotentialBasis().Mult(elementGravityPotential, quadratureGravityPotential);
 
             MFEM_VERIFY(
                 data.rotationPotential.Size() == quadraturePointCount, "Prepared hydrostatic base state has stale "
@@ -836,7 +811,8 @@ namespace mean_field::operators {
 
             MFEM_VERIFY(
                 data.baseDisplacementData.has_value() && data.compactificationData.has_value() &&
-                    static_cast<int>(data.baseMappingContexts.size()) == quadraturePointCount &&
+                    data.baseMappingContexts.GetPointCount() == quadraturePointCount &&
+                    data.baseMappingContexts.GetDimension() == dimension &&
                     data.rotationGradient.Height() == quadraturePointCount &&
                     data.rotationGradient.Width() == dimension &&
                     data.hydrostaticImbalance.Size() == quadraturePointCount,
@@ -855,7 +831,7 @@ namespace mean_field::operators {
         for (const ElementPAData &data : m_elements) {
             elementResidual.SetSize(data.enthalpyDofs.Size());
 
-            data.enthalpyBasis.MultTranspose(data.weightedResidual, elementResidual);
+            data.GetEnthalpyBasis().MultTranspose(data.weightedResidual, elementResidual);
 
             if (data.enthalpyDofTransformation != nullptr) {
                 data.enthalpyDofTransformation->TransformDual(elementResidual);
@@ -1072,7 +1048,7 @@ namespace mean_field::operators {
                                                      data.rotationPotential(quadraturePoint);
             }
             elementAction.SetSize(data.enthalpyDofs.Size());
-            data.enthalpyBasis.MultTranspose(weightedVariation, elementAction);
+            data.GetEnthalpyBasis().MultTranspose(weightedVariation, elementAction);
             if (data.enthalpyDofTransformation != nullptr) {
                 data.enthalpyDofTransformation->TransformDual(elementAction);
             }
@@ -1204,6 +1180,7 @@ namespace mean_field::operators {
         mfem::Vector elementDisplacementVariation;
         mfem::Vector weightedQuadratureVariation;
         mfem::Vector elementAction;
+        mapping::VolumeMappingContext mappingContext;
         mapping::VolumeMappingVariation variation;
 
         for (const ElementPAData &data : m_elements) {
@@ -1239,7 +1216,7 @@ namespace mean_field::operators {
             const int quadraturePointCount = data.integrationRule->GetNPoints();
 
             MFEM_VERIFY(
-                static_cast<int>(data.baseMappingContexts.size()) == quadraturePointCount &&
+                data.baseMappingContexts.GetPointCount() == quadraturePointCount &&
                     data.quadratureWeights.Size() == quadraturePointCount &&
                     data.hydrostaticImbalance.Size() == quadraturePointCount &&
                     data.rotationGradient.Height() == quadraturePointCount &&
@@ -1252,10 +1229,10 @@ namespace mean_field::operators {
 
             for (int quadraturePoint = 0; quadraturePoint < quadraturePointCount; ++quadraturePoint) {
                 const mfem::IntegrationPoint &integrationPoint = data.integrationRule->IntPoint(quadraturePoint);
+                data.baseMappingContexts.Load(quadraturePoint, mappingContext);
 
-                const mapping::MappingStatus mappingStatus     = m_domainMapper.EvaluateVolumeVariation(
-                    mappingData, directionData, *transformation, integrationPoint,
-                    data.baseMappingContexts[quadraturePoint], workspace, variation
+                const mapping::MappingStatus mappingStatus = m_domainMapper.EvaluateVolumeVariation(
+                    mappingData, directionData, *transformation, integrationPoint, mappingContext, workspace, variation
                 );
 
                 MFEM_VERIFY(
@@ -1288,7 +1265,7 @@ namespace mean_field::operators {
 
             elementAction.SetSize(data.enthalpyDofs.Size());
 
-            data.enthalpyBasis.MultTranspose(weightedQuadratureVariation, elementAction);
+            data.GetEnthalpyBasis().MultTranspose(weightedQuadratureVariation, elementAction);
 
             if (data.enthalpyDofTransformation != nullptr) {
                 data.enthalpyDofTransformation->TransformDual(elementAction);

libmeanfield/impl/operators/prepared_mass_normalization.cpp

@@ -420,17 +420,12 @@ namespace mean_field::operators {
 
             const mfem::IntegrationRule &integrationRule =
                 get_mass_normalization_rule(m_fem, densityElement, *transformation);
+            data.integrationRule = &integrationRule;
 
-            data.quadraturePoints.resize(integrationRule.GetNPoints());
-
-            for (int quadraturePoint = 0; quadraturePoint < integrationRule.GetNPoints(); ++quadraturePoint) {
-                QuadraturePointData &point = data.quadraturePoints[quadraturePoint];
-
-                point.integrationPoint     = integrationRule.IntPoint(quadraturePoint);
-
-                point.densityShape.SetSize(densityElement.GetDof());
-                densityElement.CalcShape(point.integrationPoint, point.densityShape);
-            }
+            data.densityBasis    = m_fem.GetReferenceTables().GetScalarTable(densityElement, integrationRule);
+            data.mappingContexts.SetSize(integrationRule.GetNPoints(), m_fem.mesh->Dimension());
+            data.density.SetSize(integrationRule.GetNPoints());
+            data.quadratureWeights.SetSize(integrationRule.GetNPoints());
         }
 
         int globalStellarElementCount = 0;
@@ -459,6 +454,7 @@ namespace mean_field::operators {
         true_to_local(*m_fem.displacementFes, displacement, displacementLocal);
 
         mapping::DomainMapper::Workspace workspace(m_fem.mesh->Dimension());
+        mapping::VolumeMappingContext mappingContext;
         std::optional<MassNormalizationPreparationRejection> rejection;
 
         for (ElementPAData &data : m_elements) {
@@ -491,9 +487,10 @@ namespace mean_field::operators {
 
             mfem::ElementTransformation *transformation = m_fem.mesh->GetElementTransformation(data.elementId);
 
-            for (QuadraturePointData &point : data.quadraturePoints) {
+            for (int quadraturePoint = 0; quadraturePoint < data.integrationRule->GetNPoints(); ++quadraturePoint) {
                 const mapping::MappingStatus status = m_domainMapper.EvaluateVolume(
-                    mappingData, *transformation, point.integrationPoint, workspace, point.mappingContext
+                    mappingData, *transformation, data.integrationRule->IntPoint(quadraturePoint), workspace,
+                    mappingContext
                 );
 
                 MFEM_VERIFY(
@@ -505,6 +502,9 @@ namespace mean_field::operators {
                         rejection, {.reason        = MassNormalizationPreparationRejectionReason::mapping_failure,
                                     .mappingStatus = status}
                     );
+                } else {
+                    data.mappingContexts.Store(quadraturePoint, mappingContext);
+                    data.quadratureWeights(quadraturePoint) = mappingContext.quadrature.weight;
                 }
             }
         }
@@ -537,9 +537,9 @@ namespace mean_field::operators {
                 data.densityDofTransformation->InvTransformPrimal(elementDensity);
             }
 
-            for (QuadraturePointData &point : data.quadraturePoints) {
-                point.density = elementDensity * point.densityShape;
-                if (!std::isfinite(point.density)) {
+            data.densityBasis->GetValues().Mult(elementDensity, data.density);
+            for (int quadraturePoint = 0; quadraturePoint < data.density.Size(); ++quadraturePoint) {
+                if (!std::isfinite(data.density(quadraturePoint))) {
                     retain_higher_priority_rejection(
                         rejection,
                         {.reason = MassNormalizationPreparationRejectionReason::non_finite_density_interpolation}
@@ -556,8 +556,8 @@ namespace mean_field::operators {
         std::optional<MassNormalizationPreparationRejection> localRejection;
 
         for (const ElementPAData &data : m_elements) {
-            for (const QuadraturePointData &point : data.quadraturePoints) {
-                const double contribution = point.density * point.mappingContext.quadrature.weight;
+            for (int quadraturePoint = 0; quadraturePoint < data.density.Size(); ++quadraturePoint) {
+                const double contribution = data.density(quadraturePoint) * data.quadratureWeights(quadraturePoint);
                 if (!std::isfinite(contribution) || !std::isfinite(localMass + contribution)) {
                     localRejection = {.reason = MassNormalizationPreparationRejectionReason::non_finite_assembled_mass};
                     continue;
@@ -615,6 +615,7 @@ namespace mean_field::operators {
         true_to_local(*m_fem.densityFes, densityVariation, densityVariationLocal);
 
         mfem::Vector elementDensityVariation;
+        mfem::Vector quadratureDensityVariation;
         double localAction = 0.0;
 
         for (const ElementPAData &data : m_elements) {
@@ -624,8 +625,10 @@ namespace mean_field::operators {
                 data.densityDofTransformation->InvTransformPrimal(elementDensityVariation);
             }
 
-            for (const QuadraturePointData &point : data.quadraturePoints) {
-                localAction += (elementDensityVariation * point.densityShape) * point.mappingContext.quadrature.weight;
+            quadratureDensityVariation.SetSize(data.integrationRule->GetNPoints());
+            data.densityBasis->GetValues().Mult(elementDensityVariation, quadratureDensityVariation);
+            for (int quadraturePoint = 0; quadraturePoint < quadratureDensityVariation.Size(); ++quadraturePoint) {
+                localAction += quadratureDensityVariation(quadraturePoint) * data.quadratureWeights(quadraturePoint);
             }
         }
 
@@ -650,6 +653,7 @@ namespace mean_field::operators {
 
         mapping::DomainMapper::Workspace workspace(m_fem.mesh->Dimension());
         mapping::VolumeMappingVariation variation;
+        mapping::VolumeMappingContext mappingContext;
 
         mfem::Vector elementDisplacementVariation;
         double localAction = 0.0;
@@ -681,10 +685,11 @@ namespace mean_field::operators {
 
             mfem::ElementTransformation *transformation = m_fem.mesh->GetElementTransformation(data.elementId);
 
-            for (const QuadraturePointData &point : data.quadraturePoints) {
+            for (int quadraturePoint = 0; quadraturePoint < data.integrationRule->GetNPoints(); ++quadraturePoint) {
+                data.mappingContexts.Load(quadraturePoint, mappingContext);
                 const mapping::MappingStatus status = m_domainMapper.EvaluateVolumeVariation(
-                    mappingData, directionData, *transformation, point.integrationPoint, point.mappingContext,
-                    workspace, variation
+                    mappingData, directionData, *transformation, data.integrationRule->IntPoint(quadraturePoint),
+                    mappingContext, workspace, variation
                 );
 
                 MFEM_VERIFY(
@@ -694,7 +699,7 @@ namespace mean_field::operators {
                                                                  << ", status: " << static_cast<int>(status)
                 );
 
-                localAction += point.density * variation.weight_variation;
+                localAction += data.density(quadraturePoint) * variation.weight_variation;
             }
         }
 
@@ -789,14 +794,13 @@ namespace mean_field::operators {
         localDual = 0.0;
 
         mfem::Vector elementDual;
+        mfem::Vector weightedDual;
 
         for (const ElementPAData &data : m_elements) {
             elementDual.SetSize(data.densityDofs.Size());
-            elementDual = 0.0;
-
-            for (const QuadraturePointData &point : data.quadraturePoints) {
-                elementDual.Add(residualDual * point.mappingContext.quadrature.weight, point.densityShape);
-            }
+            weightedDual = data.quadratureWeights;
+            weightedDual *= residualDual;
+            data.densityBasis->GetValues().MultTranspose(weightedDual, elementDual);
 
             if (data.densityDofTransformation != nullptr) {
                 data.densityDofTransformation->TransformDual(elementDual);
@@ -821,6 +825,7 @@ namespace mean_field::operators {
 
         mapping::DomainMapper::Workspace workspace(m_fem.mesh->Dimension());
         mapping::VolumeMappingVariation variation;
+        mapping::VolumeMappingContext mappingContext;
         mfem::Vector elementDirection;
         mfem::Vector elementDual;
 
@@ -852,10 +857,11 @@ namespace mean_field::operators {
 
                 double elementDofAction = 0.0;
 
-                for (const QuadraturePointData &point : data.quadraturePoints) {
+                for (int quadraturePoint = 0; quadraturePoint < data.integrationRule->GetNPoints(); ++quadraturePoint) {
+                    data.mappingContexts.Load(quadraturePoint, mappingContext);
                     const mapping::MappingStatus status = m_domainMapper.EvaluateVolumeVariation(
-                        mappingData, directionData, *transformation, point.integrationPoint, point.mappingContext,
-                        workspace, variation
+                        mappingData, directionData, *transformation, data.integrationRule->IntPoint(quadraturePoint),
+                        mappingContext, workspace, variation
                     );
 
                     MFEM_VERIFY(
@@ -864,7 +870,7 @@ namespace mean_field::operators {
                             << data.elementId << ", status: " << static_cast<int>(status)
                     );
 
-                    elementDofAction += point.density * variation.weight_variation;
+                    elementDofAction += data.density(quadraturePoint) * variation.weight_variation;
                 }
 
                 elementDual(elementDof) = residualDual * elementDofAction;

libmeanfield/impl/operators/prepared_pressure_force.cpp

@@ -492,9 +492,6 @@ namespace mean_field::operators {
 
         m_elements.reserve(m_fem.mesh->GetNE());
 
-        mfem::Vector enthalpyShape;
-        mfem::DenseMatrix displacementDShape;
-
         for (int elementId = 0; elementId < m_fem.mesh->GetNE(); ++elementId) {
             mfem::ElementTransformation *transformation = m_fem.mesh->GetElementTransformation(elementId);
 
@@ -559,29 +556,11 @@ namespace mean_field::operators {
                 "an unexpected vector DOF count."
             );
 
-            data.enthalpyBasis.SetSize(quadraturePointCount, enthalpyDofCount);
-
-            data.referenceTestGradients.resize(quadraturePointCount);
-
+            const fem::ReferenceTableCache &referenceTables = m_fem.GetReferenceTables();
+            data.enthalpyReferenceTable = referenceTables.GetScalarTable(enthalpyElement, *data.integrationRule);
+            data.displacementReferenceTable =
+                referenceTables.GetScalarTable(displacementElement, *data.integrationRule);
             data.physicalTestGradients.resize(quadraturePointCount);
-
-            enthalpyShape.SetSize(enthalpyDofCount);
-
-            displacementDShape.SetSize(scalarDisplacementDofCount, dimension);
-
-            for (int quadraturePoint = 0; quadraturePoint < quadraturePointCount; ++quadraturePoint) {
-                const mfem::IntegrationPoint &integrationPoint = data.integrationRule->IntPoint(quadraturePoint);
-
-                enthalpyElement.CalcShape(integrationPoint, enthalpyShape);
-
-                displacementElement.CalcDShape(integrationPoint, displacementDShape);
-
-                for (int enthalpyDof = 0; enthalpyDof < enthalpyDofCount; ++enthalpyDof) {
-                    data.enthalpyBasis(quadraturePoint, enthalpyDof) = enthalpyShape(enthalpyDof);
-                }
-
-                data.referenceTestGradients[quadraturePoint] = displacementDShape;
-            }
         }
     }
 
@@ -591,6 +570,7 @@ namespace mean_field::operators {
         true_to_local(*m_fem.displacementFes, m_baseDisplacementTrue, displacementLocal);
 
         mapping::DomainMapper::Workspace workspace(m_fem.mesh->Dimension());
+        mapping::VolumeMappingContext mappingContext;
 
         mfem::Vector elementDisplacement;
         mfem::Vector elementCompactification;
@@ -634,14 +614,15 @@ namespace mean_field::operators {
             const int quadraturePointCount = data.integrationRule->GetNPoints();
 
             MFEM_VERIFY(
-                static_cast<int>(data.referenceTestGradients.size()) == quadraturePointCount,
+                data.displacementReferenceTable != nullptr &&
+                    data.displacementReferenceTable->GetPointCount() == quadraturePointCount,
                 "Prepared pressure-force geometry has inconsistent "
                 "static gradient data."
             );
 
             data.quadratureWeights.SetSize(quadraturePointCount);
 
-            data.baseMappingContexts.resize(quadraturePointCount);
+            data.baseMappingContexts.SetSize(quadraturePointCount, m_fem.mesh->Dimension());
 
             data.physicalTestGradients.resize(quadraturePointCount);
 
@@ -650,9 +631,7 @@ namespace mean_field::operators {
 
                 transformation->SetIntPoint(&integrationPoint);
 
-                mapping::VolumeMappingContext &mappingContext = data.baseMappingContexts[quadraturePoint];
-
-                const mapping::MappingStatus mappingStatus    = m_domainMapper.EvaluateVolume(
+                const mapping::MappingStatus mappingStatus = m_domainMapper.EvaluateVolume(
                     mappingData, *transformation, integrationPoint, workspace, mappingContext
                 );
 
@@ -669,11 +648,13 @@ namespace mean_field::operators {
                     return mapping_rejection(mapping::MappingStatus::non_positive_determinant);
                 }
 
-                data.quadratureWeights(quadraturePoint)        = quadratureWeight;
+                data.baseMappingContexts.Store(quadraturePoint, mappingContext);
+                data.quadratureWeights(quadraturePoint) = quadratureWeight;
 
-                const mfem::DenseMatrix &referenceTestGradient = data.referenceTestGradients[quadraturePoint];
+                const mfem::DenseMatrix &referenceTestGradient =
+                    data.displacementReferenceTable->GetGradients(quadraturePoint);
 
-                mfem::DenseMatrix &physicalTestGradient        = data.physicalTestGradients[quadraturePoint];
+                mfem::DenseMatrix &physicalTestGradient = data.physicalTestGradients[quadraturePoint];
 
                 MFEM_VERIFY(
                     referenceTestGradient.Width() == mappingContext.quadrature.J_inv.Height() &&
@@ -714,9 +695,11 @@ namespace mean_field::operators {
                 data.enthalpyDofTransformation->InvTransformPrimal(elementEnthalpy);
             }
 
-            const int quadraturePointCount                 = data.enthalpyBasis.Height();
+            const mfem::DenseMatrix &enthalpyBasis         = data.GetEnthalpyBasis();
 
-            const int enthalpyDofCount                     = data.enthalpyBasis.Width();
+            const int quadraturePointCount                 = enthalpyBasis.Height();
+
+            const int enthalpyDofCount                     = enthalpyBasis.Width();
 
             const mfem::FiniteElement &displacementElement = *m_fem.displacementFes->GetFE(data.elementId);
 
@@ -734,7 +717,7 @@ namespace mean_field::operators {
 
             quadratureEnthalpy.SetSize(quadraturePointCount);
 
-            data.enthalpyBasis.Mult(elementEnthalpy, quadratureEnthalpy);
+            enthalpyBasis.Mult(elementEnthalpy, quadratureEnthalpy);
 
             data.pressure.SetSize(quadraturePointCount);
 
@@ -814,8 +797,8 @@ namespace mean_field::operators {
                         data.elementResidual(vectorDof) -= residualContribution;
 
                         for (int enthalpyDof = 0; enthalpyDof < enthalpyDofCount; ++enthalpyDof) {
-                            const double jacobianContribution = pressureDerivative * weightedTestGradient *
-                                                                data.enthalpyBasis(quadraturePoint, enthalpyDof);
+                            const double jacobianContribution =
+                                pressureDerivative * weightedTestGradient * enthalpyBasis(quadraturePoint, enthalpyDof);
                             if (!std::isfinite(jacobianContribution)) {
                                 materialFailure = non_finite_rejection();
                                 continue;
@@ -841,8 +824,10 @@ namespace mean_field::operators {
 
             MFEM_VERIFY(
                 data.baseDisplacementData.has_value() && data.compactificationData.has_value() &&
-                    static_cast<int>(data.baseMappingContexts.size()) == quadraturePointCount &&
-                    static_cast<int>(data.referenceTestGradients.size()) == quadraturePointCount &&
+                    data.baseMappingContexts.GetPointCount() == quadraturePointCount &&
+                    data.baseMappingContexts.GetDimension() == dimension &&
+                    data.displacementReferenceTable != nullptr &&
+                    data.displacementReferenceTable->GetPointCount() == quadraturePointCount &&
                     static_cast<int>(data.physicalTestGradients.size()) == quadraturePointCount &&
                     data.quadratureWeights.Size() == quadraturePointCount &&
                     data.pressure.Size() == quadraturePointCount,
@@ -852,7 +837,7 @@ namespace mean_field::operators {
 
             for (int quadraturePoint = 0; quadraturePoint < quadraturePointCount; ++quadraturePoint) {
                 MFEM_VERIFY(
-                    data.referenceTestGradients[quadraturePoint].Width() == dimension &&
+                    data.displacementReferenceTable->GetGradients(quadraturePoint).Width() == dimension &&
                         data.physicalTestGradients[quadraturePoint].Width() == dimension,
                     "Prepared pressure-force displacement Jacobian has "
                     "a gradient with the wrong dimension."
@@ -978,6 +963,7 @@ namespace mean_field::operators {
         mfem::Vector elementAction;
 
         mfem::DenseMatrix referenceDisplacementJacobian;
+        mfem::DenseMatrix inverseElementJacobian;
         mfem::DenseMatrix inverseElementJacobianVariation;
         mfem::DenseMatrix matrixTemporary;
         mfem::DenseMatrix physicalTestGradientVariation;
@@ -1016,7 +1002,7 @@ namespace mean_field::operators {
             );
 
             MFEM_VERIFY(
-                static_cast<int>(data.baseMappingContexts.size()) == quadraturePointCount &&
+                data.baseMappingContexts.GetPointCount() == quadraturePointCount &&
                     data.pressure.Size() == quadraturePointCount,
                 "Prepared pressure-force displacement Jacobian has "
                 "stale quadrature data."
@@ -1032,12 +1018,12 @@ namespace mean_field::operators {
             physicalTestGradientVariation.SetSize(scalarDisplacementDofCount, dimension);
 
             for (int quadraturePoint = 0; quadraturePoint < quadraturePointCount; ++quadraturePoint) {
-                const mfem::DenseMatrix &referenceTestGradient = data.referenceTestGradients[quadraturePoint];
+                const mfem::DenseMatrix &referenceTestGradient =
+                    data.displacementReferenceTable->GetGradients(quadraturePoint);
 
                 mfem::MultAtB(directionDofs, referenceTestGradient, referenceDisplacementJacobian);
 
-                const mfem::DenseMatrix &inverseElementJacobian =
-                    data.baseMappingContexts[quadraturePoint].quadrature.J_inv;
+                data.baseMappingContexts.LoadInverseJacobian(quadraturePoint, inverseElementJacobian);
                 mfem::Mult(inverseElementJacobian, referenceDisplacementJacobian, matrixTemporary);
 
                 double logarithmicJacobianVariation{0.0};

libmeanfield/impl/operators/prepared_stellar_equilibrium.cpp

@@ -1068,6 +1068,29 @@ namespace mean_field::operators {
         return m_generatedVolumeDisplacement;
     }
 
+    void PreparedStellarEquilibriumOperator::BuildVolumeDisplacementDirection(
+        const mfem::Vector &surfaceDeformationDirection,
+        mfem::Vector &volumeDisplacementDirection
+    ) const {
+        VerifyPrepared();
+        MFEM_VERIFY(
+            surfaceDeformationDirection.Size() == m_domainDeformation.parameterCount(),
+            "PreparedStellarEquilibriumOperator received a surface-deformation direction with the wrong size."
+        );
+        validate_finite_vector(
+            surfaceDeformationDirection,
+            "PreparedStellarEquilibriumOperator received a non-finite surface-deformation direction."
+        );
+        MFEM_VERIFY(
+            volumeDisplacementDirection.Size() == m_domainDeformation.volumeDisplacementSize(),
+            "PreparedStellarEquilibriumOperator received a volume-displacement workspace with the wrong size."
+        );
+
+        m_domainDeformation.applyJacobian(
+            m_surfaceDeformationParameters, surfaceDeformationDirection, volumeDisplacementDirection
+        );
+    }
+
     const mfem::Vector &PreparedStellarEquilibriumOperator::GetFullMechanicalResidual() const {
         VerifyPrepared();
         return m_fullMechanicalResidual;