#pragma once

// Include after `import mean_field;`.  This is deliberately an experiment-only
// observer: it uses the production estimator and mapper without changing either.
#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <limits>
#include <span>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>

#include <mfem.hpp>
#include <mpi.h>

namespace experiment {
    namespace geometry_detail {
        inline double Frobenius(const mfem::DenseMatrix &matrix) {
            double squared = 0.0;
            for (int row = 0; row < matrix.Height(); ++row) {
                for (int column = 0; column < matrix.Width(); ++column) {
                    squared += matrix(row, column) * matrix(row, column);
                }
            }
            return std::sqrt(squared);
        }

        inline double Difference(const mfem::DenseMatrix &left, const mfem::DenseMatrix &right) {
            double squared = 0.0;
            for (int row = 0; row < left.Height(); ++row) {
                for (int column = 0; column < left.Width(); ++column) {
                    const double difference = left(row, column) - right(row, column);
                    squared += difference * difference;
                }
            }
            return std::sqrt(squared);
        }

        inline std::array<double, 4> DeterminantCoefficients(
            const mfem::DenseMatrix &base, const mfem::DenseMatrix &direction
        ) {
            std::array<double, 4> coefficients{};
            const int dimension = base.Height();
            mfem::DenseMatrix selected(dimension);
            for (int mask = 0; mask < (1 << dimension); ++mask) {
                int degree = 0;
                for (int column = 0; column < dimension; ++column) {
                    const bool useDirection = (mask & (1 << column)) != 0;
                    degree += useDirection ? 1 : 0;
                    for (int row = 0; row < dimension; ++row) {
                        selected(row, column) = useDirection ? direction(row, column) : base(row, column);
                    }
                }
                coefficients[static_cast<std::size_t>(degree)] += selected.Det();
            }
            return coefficients;
        }

        inline double Polynomial(const std::array<double, 4> &coefficients, const double alpha) {
            return ((coefficients[3] * alpha + coefficients[2]) * alpha + coefficients[1]) * alpha + coefficients[0];
        }

        inline std::ofstream OpenCsv(const std::filesystem::path &path) {
            std::ofstream stream(path);
            if (!stream) {
                throw std::runtime_error("Cannot open geometry diagnostic output: " + path.string());
            }
            stream << std::setprecision(17);
            return stream;
        }

        inline void Coordinates(std::ostream &stream, const mfem::Vector &position) {
            for (int component = 0; component < 3; ++component) {
                stream << ',' << (component < position.Size() ? position(component) : 0.0);
            }
        }

        struct ElementReport final {
            int element{-1};
            int globalRule{-1};
            mean_field::deformation::LargestSafeNewtonStepSizeEstimate estimate{};
        };

        // Reference-mesh probes deliberately avoid DomainMapper and inverse
        // Jacobians: an exact vertex may be singular even though interior
        // quadrature points remain valid.
        inline void InspectReferenceCorners(
            mfem::ParMesh &mesh, const std::filesystem::path &outputDirectory, const std::string &label
        ) {
            auto csv = OpenCsv(outputDirectory / (label + "_reference_corner_probes.csv"));
            csv << "element,attribute,vertex,inward_fraction,xi,eta,zeta,reference_x,reference_y,reference_z,"
                   "reference_radius,reference_det,reference_sigma_min,reference_sigma_mid,reference_sigma_max\n";
            mfem::Vector position;
            for (const int element : {0, 9, 18, 27, 36, 45, 54, 63}) {
                if (element >= mesh.GetNE()) {
                    continue;
                }
                auto *transformation = mesh.GetElementTransformation(element);
                const auto *vertices = mfem::Geometries.GetVertices(transformation->GetGeometryType());
                const auto &center = mfem::Geometries.GetCenter(transformation->GetGeometryType());
                for (int vertex = 0; vertex < vertices->GetNPoints(); ++vertex) {
                    const auto &corner = vertices->IntPoint(vertex);
                    for (const double epsilon : {0.0, 1.0e-5, 1.0e-4, 0.001, 0.01, 0.05, 0.1}) {
                        mfem::IntegrationPoint point;
                        point.Set3(
                            (1.0 - epsilon) * corner.x + epsilon * center.x,
                            (1.0 - epsilon) * corner.y + epsilon * center.y,
                            (1.0 - epsilon) * corner.z + epsilon * center.z
                        );
                        transformation->SetIntPoint(&point);
                        transformation->Transform(point, position);
                        const auto &jacobian = transformation->Jacobian();
                        csv << element << ',' << transformation->Attribute << ',' << vertex << ',' << epsilon << ','
                            << point.x << ',' << point.y << ',' << point.z;
                        Coordinates(csv, position);
                        csv << ',' << position.Norml2() << ',' << jacobian.Det() << ',' << jacobian.CalcSingularvalue(2)
                            << ',' << jacobian.CalcSingularvalue(1) << ',' << jacobian.CalcSingularvalue(0) << '\n';
                    }
                }
            }
        }
    } // namespace geometry_detail

    struct GeometryReport final {
        double boundaryStep{1.0};
        double safeStep{1.0};
        int limitingElement{-1};
        int limitedElements{0};
        int elementsWithinOnePartPerMillion{0};
        int elementsWithinOnePercent{0};
        double maximumRelativeMappingError{0.0};
        double maximumRelativeDeterminantError{0.0};
        double maximumReportedPointGradient{0.0};
        int invalidDirectSamples{0};
        std::uint64_t sampleCount{0};
    };

    // Sandbox-specific inspection of the negative-corner core element.  The
    // continuous comparison describes a logical r_L^power radial extension before
    // nodal interpolation, with logical stellar-surface radius equal to one.
    // The measured direction is always the supplied production FE field.
    inline void InspectCoreDiagonal(
        const mean_field::fem::FEM &finiteElements,
        const mfem::Vector &volumeDirection,
        const double cornerSurfaceAmplitude,
        const std::filesystem::path &outputDirectory,
        const std::string &label,
        const double radialPower = 2.0
    ) {
        using namespace geometry_detail;
        auto &space = *finiteElements.displacementFes;
        auto &mesh = *finiteElements.mesh;
        auto &logicalMesh = *finiteElements.logicalReferenceMesh;
        int ranks = 0;
        MPI_Comm_size(space.GetComm(), &ranks);
        if (ranks != 1 || mesh.GetNE() < 1 || logicalMesh.GetNE() != mesh.GetNE() || mesh.GetAttribute(0) != 1 ||
            volumeDirection.Size() != space.GetTrueVSize() || !std::isfinite(cornerSurfaceAmplitude) ||
            !std::isfinite(radialPower) || radialPower <= 0.0) {
            throw std::invalid_argument("Core-diagonal probe requires a single-rank compatible core displacement.");
        }
        auto *physicalTransformation = mesh.GetElementTransformation(0);
        auto *logicalTransformation = logicalMesh.GetElementTransformation(0);
        if (physicalTransformation->GetGeometryType() != mfem::Geometry::CUBE) {
            std::cout << "CoreDiagonal[" << label << "]: skipped unrecognized non-hex core element\n";
            return;
        }
        mfem::Vector physicalPosition(3), logicalPosition(3);
        for (const double parameter : {0.0, 0.5, 1.0}) {
            mfem::IntegrationPoint point;
            point.Set3(parameter, parameter, parameter);
            logicalTransformation->Transform(point, logicalPosition);
            physicalTransformation->Transform(point, physicalPosition);
            bool recognized = logicalPosition(0) < 0.0 && physicalPosition(0) < 0.0;
            for (int component = 1; component < 3; ++component) {
                recognized = recognized && std::abs(logicalPosition(component) - logicalPosition(0)) < 1.0e-11 &&
                             std::abs(physicalPosition(component) - physicalPosition(0)) < 1.0e-11;
            }
            // Element 0 in sandbox.smesh covers logical diagonal [-1/4,-1/8].
            const double expectedLogical = -0.25 + 0.125 * parameter;
            recognized = recognized && std::abs(logicalPosition(0) - expectedLogical) < 1.0e-11;
            if (!recognized) {
                std::cout << "CoreDiagonal[" << label << "]: skipped unrecognized sandbox core diagonal\n";
                return;
            }
        }

        mfem::ParGridFunction direction(&space);
        direction.SetFromTrueDofs(volumeDirection);
        mfem::Array<int> dofs;
        auto *dofTransformation = space.GetElementVDofs(0, dofs);
        mfem::Vector elementDirection;
        direction.GetSubVector(dofs, elementDirection);
        if (dofTransformation != nullptr) {
            dofTransformation->InvTransformPrimal(elementDirection);
        }
        const auto &finiteElement = *space.GetFE(0);
        const mean_field::mapping::ElementDisplacementData directionData(finiteElement, elementDirection, space.GetOrdering());
        mfem::Vector shape(finiteElement.GetDof()), value(3), directionAlongDiagonal(3), physicalAlongDiagonal(3);
        mfem::Vector logicalAlongDiagonal(3), unitDiagonal(3), radial(3), radialDerivative(3);
        unitDiagonal = 1.0;
        mfem::DenseMatrix derivativeShape(finiteElement.GetDof(), 3), directionGradientHat(3);
        auto csv = OpenCsv(outputDirectory / (label + "_core_diagonal.csv"));
        csv << "element,s,logical_x,logical_y,logical_z,reference_x,reference_y,reference_z,logical_radius,physical_radius,"
               "corner_surface_amplitude,drL_ds,dr_ds,actual_radial_displacement,desired_radial_displacement,"
               "actual_du_radial_ds,desired_du_radial_ds,actual_radial_gradient,desired_radial_gradient,"
               "actual_du_x_ds,actual_du_y_ds,actual_du_z_ds,reference_det,reference_sigma_min,reference_sigma_max,"
               "radial_power,relative_det_coefficient_0,relative_det_coefficient_1,relative_det_coefficient_2,relative_det_coefficient_3\n";
        for (const double parameter : {0.0, 0.005, 0.010885670926971493, 0.02, 0.05, 0.1, 0.2,
                                      0.276393202250021, 0.5, 0.723606797749979, 0.9, 1.0}) {
            mfem::IntegrationPoint point;
            point.Set3(parameter, parameter, parameter);
            logicalTransformation->SetIntPoint(&point);
            physicalTransformation->SetIntPoint(&point);
            logicalTransformation->Transform(point, logicalPosition);
            physicalTransformation->Transform(point, physicalPosition);
            const auto &referenceJacobian = physicalTransformation->Jacobian();
            referenceJacobian.Mult(unitDiagonal, physicalAlongDiagonal);
            logicalTransformation->Jacobian().Mult(unitDiagonal, logicalAlongDiagonal);
            finiteElement.CalcShape(point, shape);
            finiteElement.CalcDShape(point, derivativeShape);
            directionData.GetDofMatrix().MultTranspose(shape, value);
            mfem::MultAtB(directionData.GetDofMatrix(), derivativeShape, directionGradientHat);
            directionGradientHat.Mult(unitDiagonal, directionAlongDiagonal);

            const double physicalRadius = physicalPosition.Norml2();
            const double logicalRadius = std::max({std::abs(logicalPosition(0)), std::abs(logicalPosition(1)), std::abs(logicalPosition(2))});
            const double logicalRadiusDerivative = -(logicalAlongDiagonal(0) + logicalAlongDiagonal(1) + logicalAlongDiagonal(2)) / 3.0;
            const double nan = std::numeric_limits<double>::quiet_NaN();
            double physicalRadiusDerivative = nan, actualRadial = nan, actualRadialDerivative = nan;
            if (physicalRadius > 0.0) {
                radial = physicalPosition;
                radial /= physicalRadius;
                physicalRadiusDerivative = radial * physicalAlongDiagonal;
                radialDerivative = physicalAlongDiagonal;
                radialDerivative.Add(-physicalRadiusDerivative, radial);
                radialDerivative /= physicalRadius;
                actualRadial = radial * value;
                actualRadialDerivative = radial * directionAlongDiagonal + radialDerivative * value;
            }
            const double desiredRadial = cornerSurfaceAmplitude * std::pow(logicalRadius, radialPower);
            const double desiredRadialDerivative = radialPower * cornerSurfaceAmplitude *
                                                   std::pow(logicalRadius, radialPower - 1.0) * logicalRadiusDerivative;
            const double actualGradient = std::isfinite(physicalRadiusDerivative) && physicalRadiusDerivative != 0.0
                ? actualRadialDerivative / physicalRadiusDerivative : nan;
            const double desiredGradient = std::isfinite(physicalRadiusDerivative) && physicalRadiusDerivative != 0.0
                ? desiredRadialDerivative / physicalRadiusDerivative : nan;
            csv << "0," << parameter;
            Coordinates(csv, logicalPosition);
            Coordinates(csv, physicalPosition);
            csv << ',' << logicalRadius << ',' << physicalRadius << ',' << cornerSurfaceAmplitude << ','
                << logicalRadiusDerivative << ',' << physicalRadiusDerivative << ',' << actualRadial << ',' << desiredRadial
                << ',' << actualRadialDerivative << ',' << desiredRadialDerivative << ',' << actualGradient << ',' << desiredGradient;
            Coordinates(csv, directionAlongDiagonal);
            const double referenceDeterminant = referenceJacobian.Det();
            csv << ',' << referenceDeterminant << ',' << referenceJacobian.CalcSingularvalue(2) << ','
                << referenceJacobian.CalcSingularvalue(0) << ',' << radialPower;
            // Direct total-element determinant, normalized by undeformed
            // reference volume. This remains evaluable at a vertex without
            // constructing a potentially ill-conditioned inverse Jacobian.
            const auto determinant = DeterminantCoefficients(referenceJacobian, directionGradientHat);
            for (const double coefficient : determinant) {
                csv << ',' << (referenceDeterminant != 0.0 ? coefficient / referenceDeterminant : nan);
            }
            csv << '\n';
        }
    }

    // Per-element calls to the collective production estimator are intentional.
    // Restrict this diagnostic to one rank: element counts differ on MPI ranks,
    // so independently iterating them would mismatch estimator collectives.
    inline GeometryReport InspectGeometry(
        const mean_field::mapping::DomainMapper &mapper,
        const mfem::ParFiniteElementSpace &displacementSpace,
        const mfem::ParGridFunction &compactification,
        const mfem::Vector &acceptedVolumeDisplacement,
        const mfem::Vector &volumeDirection,
        const std::span<const mean_field::deformation::NewtonStepGeometryRule> productionRules,
        const std::filesystem::path &outputDirectory,
        const std::string &label,
        const std::size_t detailedElementCount = 12
    ) {
        using namespace mean_field;
        using namespace geometry_detail;
        int ranks = 0;
        MPI_Comm_size(displacementSpace.GetComm(), &ranks);
        if (ranks != 1) {
            throw std::invalid_argument("Per-element geometry experiment requires exactly one MPI rank.");
        }
        if (mapper.GetDimension() != 3) {
            throw std::invalid_argument("Geometry experiment currently requires three spatial dimensions.");
        }
        std::filesystem::create_directories(outputDirectory);
        auto *mesh = displacementSpace.GetParMesh();
        if (label.find("uniform") != std::string::npos) {
            InspectReferenceCorners(*mesh, outputDirectory, label);
        }
        std::vector<std::vector<deformation::NewtonStepGeometryRule>> elementRules(mesh->GetNE());
        std::vector<std::vector<int>> ruleIndices(mesh->GetNE());
        for (std::size_t index = 0; index < productionRules.size(); ++index) {
            const auto &rule = productionRules[index];
            elementRules.at(rule.element).push_back(rule);
            ruleIndices.at(rule.element).push_back(static_cast<int>(index));
        }

        GeometryReport report;
        std::vector<ElementReport> elements;
        elements.reserve(mesh->GetNE());
        for (int element = 0; element < mesh->GetNE(); ++element) {
            if (elementRules[element].empty()) {
                continue;
            }
            const auto estimate = deformation::estimate_largest_safe_newton_step_size(
                mapper, displacementSpace, compactification, acceptedVolumeDisplacement, volumeDirection,
                elementRules[element]
            );
            const int globalRule = estimate.limitingRule < 0 ? -1 : ruleIndices[element].at(estimate.limitingRule);
            elements.push_back({element, globalRule, estimate});
            report.sampleCount += estimate.sampledQuadraturePointCount;
            if (estimate.limitedByGeometry) {
                ++report.limitedElements;
                if (report.limitingElement < 0 || estimate.boundaryStepSize < report.boundaryStep) {
                    report.boundaryStep = estimate.boundaryStepSize;
                    report.safeStep = estimate.stepSize;
                    report.limitingElement = element;
                }
            }
        }
        std::stable_sort(elements.begin(), elements.end(), [](const auto &left, const auto &right) {
            return left.estimate.boundaryStepSize < right.estimate.boundaryStepSize;
        });
        for (const auto &element : elements) {
            if (element.estimate.limitedByGeometry) {
                report.elementsWithinOnePartPerMillion +=
                    element.estimate.boundaryStepSize <= report.boundaryStep * (1.0 + 1.0e-6);
                report.elementsWithinOnePercent +=
                    element.estimate.boundaryStepSize <= report.boundaryStep * 1.01;
            }
        }

        auto elementCsv = OpenCsv(outputDirectory / (label + "_geometry_elements.csv"));
        elementCsv << "rank,element,attribute,compactified,limited,boundary_step,safe_step,rule_index,point_index,"
                      "rule_points,xi,eta,zeta,reference_x,reference_y,reference_z,physical_x,physical_y,physical_z,"
                      "physical_radius,min_det_accepted_samples,min_det_full_step_samples,limiter_det_accepted,"
                      "limiter_sigma_min_accepted,limiter_sigma_max_accepted,limiter_direction_gradient_frobenius,"
                      "limiter_direction_radial,limiter_direction_tangential,"
                      "limiter_radial_gradient,limiter_tangential_gradient_trace,"
                      "reference_element_det,reference_element_sigma_min,reference_element_sigma_max,"
                      "total_physical_element_det,total_physical_element_sigma_min,total_physical_element_sigma_max,"
                      "det_coefficient_0,det_coefficient_1,det_coefficient_2,det_coefficient_3\n";
        auto sampleCsv = OpenCsv(outputDirectory / (label + "_geometry_mapping_checks.csv"));
        sampleCsv << "element,rule_index,point_index,alpha,boundary_fraction,status,polynomial_det,direct_det,"
                     "relative_det_error,relative_mapping_matrix_error,"
                     "mapped_sigma_min,mapped_sigma_max,displacement_det,displacement_sigma_min\n";
        auto matrixCsv = OpenCsv(outputDirectory / (label + "_geometry_limiting_matrices.csv"));
        matrixCsv << "element,matrix,row,column,value\n";

        mfem::Vector acceptedLocal(displacementSpace.GetVSize());
        mfem::Vector directionLocal(displacementSpace.GetVSize());
        const auto *prolongation = displacementSpace.GetProlongationMatrix();
        if (prolongation != nullptr) {
            prolongation->Mult(acceptedVolumeDisplacement, acceptedLocal);
            prolongation->Mult(volumeDirection, directionLocal);
        } else {
            acceptedLocal = acceptedVolumeDisplacement;
            directionLocal = volumeDirection;
        }
        const auto *compactificationSpace = compactification.ParFESpace();
        mfem::Array<int> displacementDofs, compactificationDofs;
        mfem::Vector acceptedElement, directionElement, compactificationElement, trialElement;
        mapping::DomainMapper::Workspace workspace(mapper.GetDimension());
        mapping::MappingPointContext base, direct;
        mapping::MappingPointVariation variation;

        for (std::size_t sortedIndex = 0; sortedIndex < elements.size(); ++sortedIndex) {
            const auto &entry = elements[sortedIndex];
            const int element = entry.element;
            auto *transformation = mesh->GetElementTransformation(element);
            const auto *displacementFe = displacementSpace.GetFE(element);
            const auto *compactificationFe = compactificationSpace->GetFE(element);
            auto *displacementTransform = displacementSpace.GetElementVDofs(element, displacementDofs);
            auto *compactificationTransform = compactificationSpace->GetElementDofs(element, compactificationDofs);
            acceptedLocal.GetSubVector(displacementDofs, acceptedElement);
            directionLocal.GetSubVector(displacementDofs, directionElement);
            compactification.GetSubVector(compactificationDofs, compactificationElement);
            if (displacementTransform != nullptr) {
                displacementTransform->InvTransformPrimal(acceptedElement);
                displacementTransform->InvTransformPrimal(directionElement);
            }
            if (compactificationTransform != nullptr) {
                compactificationTransform->InvTransformPrimal(compactificationElement);
            }
            const mapping::ElementDisplacementData baseData(
                *displacementFe, acceptedElement, displacementSpace.GetOrdering()
            );
            const mapping::ElementDisplacementData directionData(
                *displacementFe, directionElement, displacementSpace.GetOrdering()
            );
            const mapping::ElementCompactificationData compactificationData(*compactificationFe, compactificationElement);
            const mapping::ElementMappingData mappingData{baseData, compactificationData};
            const auto &point = entry.globalRule >= 0
                ? productionRules[entry.globalRule].integrationRule->IntPoint(entry.estimate.limitingQuadraturePoint)
                : mfem::Geometries.GetCenter(transformation->GetGeometryType());
            if (mapper.EvaluatePoint(mappingData, *transformation, point, workspace, base) != mapping::MappingStatus::valid ||
                mapper.EvaluatePointVariation(mappingData, directionData, *transformation, point, base, workspace, variation) !=
                    mapping::MappingStatus::valid) {
                throw std::runtime_error("Geometry diagnostic could not evaluate element " + std::to_string(element));
            }
            const auto coefficients = DeterminantCoefficients(base.mapping_jacobian, variation.mapping_jacobian_variation);
            const mfem::DenseMatrix referenceJacobian(transformation->Jacobian());
            mfem::DenseMatrix totalPhysicalJacobian(3);
            mfem::Mult(base.mapping_jacobian, referenceJacobian, totalPhysicalJacobian);
            const double radius = base.physical_position.Norml2();
            double radialDirection = 0.0;
            if (radius > 0.0) {
                radialDirection = (base.physical_position * variation.physical_position_variation) / radius;
            }
            const double directionNorm = variation.physical_position_variation.Norml2();
            const double tangentialDirection = std::sqrt(std::max(0.0, directionNorm * directionNorm - radialDirection * radialDirection));
            mfem::DenseMatrix physicalGradient(3);
            mfem::Mult(variation.mapping_jacobian_variation, base.inverse_mapping_jacobian, physicalGradient);
            double radialGradient = 0.0;
            if (radius > 0.0) {
                for (int row = 0; row < 3; ++row) {
                    for (int column = 0; column < 3; ++column) {
                        radialGradient += base.physical_position(row) * physicalGradient(row, column) *
                                          base.physical_position(column) / (radius * radius);
                    }
                }
            }
            const double tangentialTrace = physicalGradient(0, 0) + physicalGradient(1, 1) + physicalGradient(2, 2) - radialGradient;
            const double gradientNorm = Frobenius(variation.displacement_jacobian_variation);
            report.maximumReportedPointGradient = std::max(report.maximumReportedPointGradient, gradientNorm);
            elementCsv << "0," << element << ',' << transformation->Attribute << ',' << base.compactified << ','
                       << entry.estimate.limitedByGeometry << ',' << entry.estimate.boundaryStepSize << ','
                       << entry.estimate.stepSize << ',' << entry.globalRule << ',' << entry.estimate.limitingQuadraturePoint << ','
                       << (entry.globalRule >= 0 ? productionRules[entry.globalRule].integrationRule->GetNPoints() : 0) << ','
                       << point.x << ',' << point.y << ',' << point.z;
            Coordinates(elementCsv, base.reference_position);
            Coordinates(elementCsv, base.physical_position);
            elementCsv << ',' << radius << ',' << entry.estimate.minimumDeterminantAtAcceptedState << ','
                       << entry.estimate.minimumDeterminantAtMaximumStepSize << ',' << base.mapping_determinant << ','
                       << base.mapping_jacobian.CalcSingularvalue(2) << ',' << base.mapping_jacobian.CalcSingularvalue(0) << ','
                       << gradientNorm << ',' << radialDirection << ',' << tangentialDirection << ',' << radialGradient << ','
                       << tangentialTrace << ',' << referenceJacobian.Det() << ',' << referenceJacobian.CalcSingularvalue(2)
                       << ',' << referenceJacobian.CalcSingularvalue(0) << ',' << totalPhysicalJacobian.Det() << ','
                       << totalPhysicalJacobian.CalcSingularvalue(2) << ',' << totalPhysicalJacobian.CalcSingularvalue(0);
            for (const double coefficient : coefficients) {
                elementCsv << ',' << coefficient;
            }
            elementCsv << '\n';

            if (sortedIndex >= detailedElementCount) {
                continue;
            }
            const std::array<std::pair<const char *, const mfem::DenseMatrix *>, 6> matrices{{
                {"base_mapping", &base.mapping_jacobian},
                {"direction_mapping", &variation.mapping_jacobian_variation},
                {"direction_displacement", &variation.displacement_jacobian_variation},
                {"physical_direction_gradient", &physicalGradient},
                {"reference_element", &referenceJacobian},
                {"total_physical_element", &totalPhysicalJacobian}
            }};
            for (const auto &[name, matrix] : matrices) {
                for (int row = 0; row < 3; ++row) {
                    for (int column = 0; column < 3; ++column) {
                        matrixCsv << element << ',' << name << ',' << row << ',' << column << ',' << (*matrix)(row, column) << '\n';
                    }
                }
            }
            for (const double fraction : {0.0, 0.25, 0.5, 0.9, 0.99, 0.999}) {
                const double alpha = fraction * entry.estimate.boundaryStepSize;
                trialElement = acceptedElement;
                trialElement.Add(alpha, directionElement);
                const mapping::ElementDisplacementData trialData(*displacementFe, trialElement, displacementSpace.GetOrdering());
                const mapping::ElementMappingData trialMappingData{trialData, compactificationData};
                const auto status = mapper.EvaluatePoint(trialMappingData, *transformation, point, workspace, direct);
                mfem::DenseMatrix affine(base.mapping_jacobian);
                affine.Add(alpha, variation.mapping_jacobian_variation);
                const double predictedDeterminant = Polynomial(coefficients, alpha);
                const double nan = std::numeric_limits<double>::quiet_NaN();
                double relativeMatrixError = nan, relativeDeterminantError = nan;
                double minimumSingular = nan, maximumSingular = nan, displacementDeterminant = nan, displacementMinimumSingular = nan;
                if (status == mapping::MappingStatus::valid) {
                    relativeMatrixError = Difference(direct.mapping_jacobian, affine) / std::max(Frobenius(affine), 1.0e-300);
                    relativeDeterminantError = std::abs(direct.mapping_determinant - predictedDeterminant) /
                                              std::max(std::abs(base.mapping_determinant), 1.0e-300);
                    minimumSingular = direct.mapping_jacobian.CalcSingularvalue(2);
                    maximumSingular = direct.mapping_jacobian.CalcSingularvalue(0);
                    // DomainMapper's displacement_jacobian already includes I.
                    const mfem::DenseMatrix &displacementMapping = direct.displacement_jacobian;
                    displacementDeterminant = displacementMapping.Det();
                    displacementMinimumSingular = displacementMapping.CalcSingularvalue(2);
                    report.maximumRelativeMappingError = std::max(report.maximumRelativeMappingError, relativeMatrixError);
                    report.maximumRelativeDeterminantError = std::max(report.maximumRelativeDeterminantError, relativeDeterminantError);
                } else {
                    ++report.invalidDirectSamples;
                }
                sampleCsv << element << ',' << entry.globalRule << ',' << entry.estimate.limitingQuadraturePoint << ','
                          << alpha << ',' << fraction << ',' << static_cast<int>(status) << ',' << predictedDeterminant << ','
                          << (status == mapping::MappingStatus::valid ? direct.mapping_determinant : nan) << ','
                          << relativeDeterminantError << ',' << relativeMatrixError << ',' << minimumSingular << ','
                          << maximumSingular << ',' << displacementDeterminant << ',' << displacementMinimumSingular << '\n';
            }
        }
        std::cout << "Geometry[" << label << "]: boundary=" << report.boundaryStep << " safe=" << report.safeStep
                  << " limiter=" << report.limitingElement << " limited_elements=" << report.limitedElements
                  << " ties_1ppm=" << report.elementsWithinOnePartPerMillion << " ties_1pct=" << report.elementsWithinOnePercent
                  << " samples=" << report.sampleCount << " max_mapping_error=" << report.maximumRelativeMappingError
                  << " max_det_error=" << report.maximumRelativeDeterminantError
                  << " invalid_direct_samples=" << report.invalidDirectSamples << std::endl;
        return report;
    }
} // namespace experiment