#include <algorithm>
#include <array>
#include <cmath>
#include <limits>
#include <stdexcept>

#include <catch2/catch_test_macros.hpp>

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

import mean_field;
import test_helpers;

namespace nodal_radial_surface_test_utils {
    namespace deformation      = mean_field::deformation;
    namespace domain           = mean_field::utils::domain;
    namespace field            = mean_field::field;

    using Schema               = domain::CoreEnvelopeVacuumDomainSchema;

    using PlanarBoundarySchema = domain::DomainSchema<
        domain::MaterialList<>,
        domain::BoundaryList<domain::BoundaryAttribute<domain::StellarSurface, 1>>,
        domain::RelationList<>>;

    [[nodiscard]] mfem::Vector referenceCenter(const int spatialDimension) {
        mfem::Vector center(spatialDimension);
        center = 0.0;
        return center;
    }

    [[nodiscard]] deformation::PreparedNodalRadialSurface makePreparedSurface(const mean_field::fem::FEM &fem) {
        const field::ScalarBoundaryDofMap surfaceDofMap =
            field::make_stellar_surface_scalar_dof_map<Schema>(*fem.surfaceDeformationFes);
        const deformation::SurfaceDeformationCompilationContext context{*fem.surfaceDeformationFes, surfaceDofMap};

        return deformation::compileSurfaceDeformationPrescription(
            deformation::NodalRadialSurface{referenceCenter(fem.mesh->SpaceDimension())}, context
        );
    }

    [[nodiscard]] mfem::Vector projectReferenceSurfacePositions(
        mfem::ParFiniteElementSpace &scalarFiniteElementSpace,
        const field::ScalarBoundaryDofMap &surfaceDofMap
    ) {
        const mfem::Mesh *mesh = scalarFiniteElementSpace.GetMesh();
        REQUIRE(mesh != nullptr);

        const int spatialDimension = mesh->SpaceDimension();
        mfem::Vector positions(spatialDimension * surfaceDofMap.local_size());
        mfem::ParGridFunction coordinateField(&scalarFiniteElementSpace);

        for (int component = 0; component < spatialDimension; ++component) {
            mfem::FunctionCoefficient coefficient([component](const mfem::Vector &position) {
                return position(component);
            });
            coordinateField.ProjectCoefficient(coefficient);

            mfem::Vector coordinateTrueDofs;
            coordinateField.GetTrueDofs(coordinateTrueDofs);
            const mfem::Vector surfaceCoordinates = surfaceDofMap.gather(coordinateTrueDofs);

            for (int surfaceDof = 0; surfaceDof < surfaceDofMap.local_size(); ++surfaceDof) {
                positions(spatialDimension * surfaceDof + component) = surfaceCoordinates(surfaceDof);
            }
        }

        return positions;
    }

    [[nodiscard]] double relativeError(
        const mfem::Vector &actual,
        const mfem::Vector &expected
    ) {
        REQUIRE(actual.Size() == expected.Size());

        mfem::Vector difference(actual);
        difference -= expected;

        return difference.Norml2() / std::max(expected.Norml2(), std::numeric_limits<double>::epsilon());
    }

    [[nodiscard]] double globalMeanReferenceRadius(const deformation::PreparedNodalRadialSurface &surface) {
        double localRadiusSum = 0.0;
        for (int parameterDof = 0; parameterDof < surface.parameterCount(); ++parameterDof) {
            localRadiusSum += surface.referenceRadius(parameterDof);
        }

        const long long localCount = surface.parameterCount();
        double globalRadiusSum     = 0.0;
        long long globalCount      = 0;
        MPI_Allreduce(&localRadiusSum, &globalRadiusSum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
        MPI_Allreduce(&localCount, &globalCount, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD);
        REQUIRE(globalCount > 0);
        return globalRadiusSum / static_cast<double>(globalCount);
    }

    struct BinaryRochePotential final {
        double primaryMass;
        double companionMass;
        double separation;

        [[nodiscard]] double operator()(
            const double x,
            const double y,
            const double z
        ) const {
            const double primaryDistance   = std::sqrt(x * x + y * y + z * z);
            const double companionOffset   = x - separation;
            const double companionDistance = std::sqrt(companionOffset * companionOffset + y * y + z * z);
            const double totalMass         = primaryMass + companionMass;
            const double centerOfMassX     = separation * companionMass / totalMass;
            const double angularSpeed2     = totalMass / (separation * separation * separation);
            const double rotationRadius2   = (x - centerOfMassX) * (x - centerOfMassX) + y * y;

            return -primaryMass / primaryDistance - companionMass / companionDistance -
                   0.5 * angularSpeed2 * rotationRadius2;
        }

        [[nodiscard]] double firstRadialIntersection(
            const double directionX,
            const double directionY,
            const double directionZ,
            const double targetPotential,
            const double referenceScale
        ) const {
            double lowerRadius = 1.0e-8 * referenceScale;
            double lowerValue =
                (*this)(lowerRadius * directionX, lowerRadius * directionY, lowerRadius * directionZ) - targetPotential;
            REQUIRE(lowerValue < 0.0);

            double upperRadius           = lowerRadius;
            double upperValue            = lowerValue;
            constexpr int bracketSamples = 512;
            for (int sample = 1; sample <= bracketSamples && upperValue <= 0.0; ++sample) {
                upperRadius = 0.45 * separation * static_cast<double>(sample) / bracketSamples;
                upperValue  = (*this)(upperRadius * directionX, upperRadius * directionY, upperRadius * directionZ) -
                              targetPotential;
            }
            REQUIRE(upperValue > 0.0);

            for (int iteration = 0; iteration < 80; ++iteration) {
                const double middleRadius = 0.5 * (lowerRadius + upperRadius);
                const double middleValue =
                    (*this)(middleRadius * directionX, middleRadius * directionY, middleRadius * directionZ) -
                    targetPotential;
                if (middleValue > 0.0) {
                    upperRadius = middleRadius;
                } else {
                    lowerRadius = middleRadius;
                }
            }
            return 0.5 * (lowerRadius + upperRadius);
        }
    };
} // namespace nodal_radial_surface_test_utils

TEST_CASE(
    "Nodal Radial Surface Satisfies The Surface Deformation Contract And Validates Its Reference Center",
    tags::nodal_radial_surface_validation
) {
    namespace deformation = mean_field::deformation;
    namespace field       = mean_field::field;

    STATIC_CHECK(deformation::SurfaceDeformationPrescription<deformation::NodalRadialSurface>);
    STATIC_CHECK(deformation::PreparedSurfaceDeformationPrescription<deformation::PreparedNodalRadialSurface>);
    STATIC_CHECK(
        deformation::SurfaceDeformationCompilable<
            deformation::NodalRadialSurface, deformation::SurfaceDeformationCompilationContext>
    );

    mfem::Vector center(3);
    center = 0.0;
    const deformation::NodalRadialSurface prescription{center};
    const deformation::SurfaceDeformationDescriptor descriptor = prescription.descriptor();

    CHECK(descriptor.name == "NodalRadialSurface");
    CHECK(descriptor.spatialDimension == 3);
    CHECK(descriptor.motionKind == deformation::SurfaceMotionKind::Radial);
    CHECK(descriptor.linearOnReferenceGeometry);
    CHECK(descriptor.requiresStarShapedReferenceSurface);
    CHECK(descriptor.supportsExactNewtonLinearization());
    CHECK(descriptor.translationTreatment == deformation::GeometricGaugeTreatment::Retained);
    CHECK(descriptor.orientationTreatment == deformation::GeometricGaugeTreatment::Retained);

    mfem::Vector emptyCenter;
    CHECK_THROWS_AS(deformation::NodalRadialSurface{emptyCenter}, std::invalid_argument);

    mfem::Vector nonfiniteCenter(3);
    nonfiniteCenter    = 0.0;
    nonfiniteCenter(1) = std::numeric_limits<double>::quiet_NaN();
    CHECK_THROWS_AS(deformation::NodalRadialSurface{nonfiniteCenter}, std::invalid_argument);

    mfem::Mesh serialMesh = mfem::Mesh::MakeCartesian2D(4, 3, mfem::Element::QUADRILATERAL, true, 2.0, 1.5);
    mfem::ParMesh mesh(MPI_COMM_WORLD, serialMesh);
    mfem::H1_FECollection finiteElementCollection(2, mesh.Dimension());
    mfem::ParFiniteElementSpace finiteElementSpace(&mesh, &finiteElementCollection);

    const field::ScalarBoundaryDofMap surfaceDofMap =
        field::make_stellar_surface_scalar_dof_map<nodal_radial_surface_test_utils::PlanarBoundarySchema>(
            finiteElementSpace
        );
    const deformation::SurfaceDeformationCompilationContext context{finiteElementSpace, surfaceDofMap};

    CHECK_THROWS_AS(deformation::compileSurfaceDeformationPrescription(prescription, context), std::invalid_argument);
}

TEST_CASE(
    "Nodal Radial Surface Resolves Oblate Rotating And Binary Roche Envelopes Between Surface Nodes",
    tags::nodal_radial_surface_analytic
) {
    mean_field::utils::Args args = test_utils::setup_args();
    mean_field::fem::FEM fem     = mean_field::fem::setup_fem(args.mesh_file, args, 0);
    REQUIRE(fem.okay());

    const mean_field::deformation::PreparedNodalRadialSurface surface =
        nodal_radial_surface_test_utils::makePreparedSurface(fem);
    REQUIRE(surface.spatialDimension() == 3);

    const double referenceScale   = nodal_radial_surface_test_utils::globalMeanReferenceRadius(surface);
    const double equatorialRadius = 1.08 * referenceScale;
    const double polarRadius      = 0.94 * referenceScale;

    mfem::Vector rotatingParameters(surface.parameterCount());
    for (int parameterDof = 0; parameterDof < surface.parameterCount(); ++parameterDof) {
        const double directionX = surface.radialDirection(parameterDof, 0);
        const double directionY = surface.radialDirection(parameterDof, 1);
        const double directionZ = surface.radialDirection(parameterDof, 2);
        const double targetRadius =
            1.0 / std::sqrt(
                      (directionX * directionX + directionY * directionY) / (equatorialRadius * equatorialRadius) +
                      directionZ * directionZ / (polarRadius * polarRadius)
                  );
        rotatingParameters(parameterDof) = targetRadius - surface.referenceRadius(parameterDof);
    }

    mfem::Vector rotatingDisplacement(surface.surfaceDisplacementSize());
    surface.buildSurfaceDisplacement(rotatingParameters, rotatingDisplacement);

    constexpr double geometricTolerance = 3.0e-12;
    for (int parameterDof = 0; parameterDof < surface.parameterCount(); ++parameterDof) {
        double movedPosition[3]{};
        for (int component = 0; component < 3; ++component) {
            const double direction   = surface.radialDirection(parameterDof, component);
            movedPosition[component] = surface.referenceRadius(parameterDof) * direction +
                                       rotatingDisplacement(surface.surfaceDisplacementDof(parameterDof, component));
        }
        const double ellipsoidLevel = (movedPosition[0] * movedPosition[0] + movedPosition[1] * movedPosition[1]) /
                                          (equatorialRadius * equatorialRadius) +
                                      movedPosition[2] * movedPosition[2] / (polarRadius * polarRadius);
        CHECK(std::abs(ellipsoidLevel - 1.0) <= geometricTolerance);
    }

    const nodal_radial_surface_test_utils::BinaryRochePotential rochePotential{
        .primaryMass = 1.0, .companionMass = 0.7, .separation = 4.0 * referenceScale
    };
    const double targetPotential = rochePotential(0.0, 0.0, referenceScale);

    mfem::Vector rocheParameters(surface.parameterCount());
    for (int parameterDof = 0; parameterDof < surface.parameterCount(); ++parameterDof) {
        const double targetRadius = rochePotential.firstRadialIntersection(
            surface.radialDirection(parameterDof, 0), surface.radialDirection(parameterDof, 1),
            surface.radialDirection(parameterDof, 2), targetPotential, referenceScale
        );
        rocheParameters(parameterDof) = targetRadius - surface.referenceRadius(parameterDof);
    }

    mfem::Vector rocheDisplacement(surface.surfaceDisplacementSize());
    surface.buildSurfaceDisplacement(rocheParameters, rocheDisplacement);

    const double potentialTolerance = 2.0e-11 * std::max(1.0, std::abs(targetPotential));
    for (int parameterDof = 0; parameterDof < surface.parameterCount(); ++parameterDof) {
        double movedPosition[3]{};
        for (int component = 0; component < 3; ++component) {
            const double direction   = surface.radialDirection(parameterDof, component);
            movedPosition[component] = surface.referenceRadius(parameterDof) * direction +
                                       rocheDisplacement(surface.surfaceDisplacementDof(parameterDof, component));
        }
        CHECK(
            std::abs(rochePotential(movedPosition[0], movedPosition[1], movedPosition[2]) - targetPotential) <=
            potentialTolerance
        );
    }

    auto domainDeformation = mean_field::deformation::compileDomainDeformation(
        mean_field::deformation::NodalRadialSurface{
            nodal_radial_surface_test_utils::referenceCenter(fem.mesh->SpaceDimension())
        },
        mean_field::deformation::PowerLawRadialInteriorExtension{},
        mean_field::deformation::FixedInfinityRadialVacuumExtension{}, fem
    );
    mfem::Vector volumeDisplacement(domainDeformation.volumeDisplacementSize());
    mfem::ParGridFunction displacementField(fem.displacementFes.get());
    mean_field::mapping::GridFunctionMappingEvaluator mappingEvaluator(
        *fem.domainMapperStateless, displacementField, *fem.compactificationCoordinate
    );

    const auto sampleRepresentationError = [&](const mfem::Vector &parameters, const auto &pointError) {
        const auto geometry = domainDeformation.buildValidatedVolumeDisplacement(parameters, volumeDisplacement);
        REQUIRE(geometry.isOrientationPreserving());
        displacementField.SetFromTrueDofs(volumeDisplacement);
        mappingEvaluator.InvalidateCache();

        double localMaximumError          = 0.0;
        double localErrorSquared          = 0.0;
        double localSurfaceArea           = 0.0;
        long long localSamples            = 0;
        const int stellarSurfaceAttribute = nodal_radial_surface_test_utils::Schema::template boundary_attribute<
            nodal_radial_surface_test_utils::domain::StellarSurface>();

        for (int boundaryElement = 0; boundaryElement < fem.mesh->GetNBE(); ++boundaryElement) {
            if (fem.mesh->GetBdrAttribute(boundaryElement) != stellarSurfaceAttribute) {
                continue;
            }

            /*
             * The stellar surface is a material interface, not an exterior
             * boundary of the complete compactified mesh. Resolve the tagged
             * boundary element to its underlying mesh face so this sampling
             * path works for both internal interfaces and true exterior
             * boundaries.
             */
            const int face = fem.mesh->GetBdrElementFaceIndex(boundaryElement);
            REQUIRE(face >= 0);
            mfem::FaceElementTransformations *transformation = fem.mesh->GetFaceElementTransformations(face);
            REQUIRE(transformation != nullptr);
            REQUIRE(transformation->Elem1 != nullptr);
            const mfem::IntegrationRule &rule =
                mfem::IntRules.Get(transformation->GetGeometryType(), 2 * mean_field::field::Displacement::vectorOrder);

            for (int point = 0; point < rule.GetNPoints(); ++point) {
                mean_field::mapping::FaceMappingContext context;
                const mfem::IntegrationPoint &integrationPoint = rule.IntPoint(point);
                REQUIRE(
                    mappingEvaluator.EvaluateFace(
                        *transformation, mean_field::mapping::FaceElementSide::element_1, integrationPoint, context
                    ) == mean_field::mapping::MappingStatus::valid
                );

                const double error = pointError(context.mapping.physical_position);
                REQUIRE(std::isfinite(error));
                REQUIRE(context.physical_surface_weight > 0.0);
                localMaximumError = std::max(localMaximumError, error);
                localErrorSquared += error * error * context.physical_surface_weight;
                localSurfaceArea += context.physical_surface_weight;
                ++localSamples;
            }
        }

        double globalMaximumError = 0.0;
        double globalErrorSquared = 0.0;
        double globalSurfaceArea  = 0.0;
        long long globalSamples   = 0;
        MPI_Allreduce(&localMaximumError, &globalMaximumError, 1, MPI_DOUBLE, MPI_MAX, fem.mesh->GetComm());
        MPI_Allreduce(&localErrorSquared, &globalErrorSquared, 1, MPI_DOUBLE, MPI_SUM, fem.mesh->GetComm());
        MPI_Allreduce(&localSurfaceArea, &globalSurfaceArea, 1, MPI_DOUBLE, MPI_SUM, fem.mesh->GetComm());
        MPI_Allreduce(&localSamples, &globalSamples, 1, MPI_LONG_LONG, MPI_SUM, fem.mesh->GetComm());
        REQUIRE(globalSamples > 0);
        REQUIRE(globalSurfaceArea > 0.0);
        return std::array<double, 2>{globalMaximumError, std::sqrt(globalErrorSquared / globalSurfaceArea)};
    };

    const auto rotatingErrors =
        sampleRepresentationError(rotatingParameters, [equatorialRadius, polarRadius](const mfem::Vector &position) {
            return std::abs(
                (position(0) * position(0) + position(1) * position(1)) / (equatorialRadius * equatorialRadius) +
                position(2) * position(2) / (polarRadius * polarRadius) - 1.0
            );
        });
    const auto rocheErrors =
        sampleRepresentationError(rocheParameters, [&rochePotential, targetPotential](const mfem::Vector &position) {
            return std::abs(rochePotential(position(0), position(1), position(2)) - targetPotential) /
                   std::max(1.0, std::abs(targetPotential));
        });

    INFO("Between-node oblate surface maximum level-set error = " << rotatingErrors[0]);
    INFO("Between-node oblate surface RMS level-set error = " << rotatingErrors[1]);
    INFO("Between-node Roche surface maximum normalized potential error = " << rocheErrors[0]);
    INFO("Between-node Roche surface RMS normalized potential error = " << rocheErrors[1]);
    CHECK(rotatingErrors[0] < 5.0e-3);
    CHECK(rotatingErrors[1] < 1.0e-3);
    CHECK(rocheErrors[0] < 2.0e-2);
    CHECK(rocheErrors[1] < 5.0e-3);

    const double companionFacingRadius =
        rochePotential.firstRadialIntersection(1.0, 0.0, 0.0, targetPotential, referenceScale);
    const double companionOpposingRadius =
        rochePotential.firstRadialIntersection(-1.0, 0.0, 0.0, targetPotential, referenceScale);
    CHECK(std::abs(companionFacingRadius - companionOpposingRadius) > 1.0e-3 * referenceScale);
}

TEST_CASE(
    "Nodal Radial Surface Produces The Requested Physical Radial Displacement At Every Surface Coordinate",
    tags::nodal_radial_surface_analytic
) {
    namespace deformation        = mean_field::deformation;

    mean_field::utils::Args args = test_utils::setup_args();
    mean_field::fem::FEM fem     = mean_field::fem::setup_fem(args.mesh_file, args, 0);
    REQUIRE(fem.okay());

    const deformation::PreparedNodalRadialSurface prepared = nodal_radial_surface_test_utils::makePreparedSurface(fem);

    REQUIRE(prepared.parameterCount() > 0);
    CHECK(prepared.spatialDimension() == fem.mesh->SpaceDimension());
    CHECK(prepared.surfaceDisplacementSize() == prepared.spatialDimension() * prepared.parameterCount());
    CHECK(
        prepared.globalSurfaceDisplacementSize() ==
        static_cast<long long>(prepared.spatialDimension()) * prepared.globalParameterCount()
    );
    CHECK(
        prepared.globalSurfaceDisplacementOffset() ==
        static_cast<long long>(prepared.spatialDimension()) * prepared.globalParameterOffset()
    );

    mfem::Vector parameters(prepared.parameterCount());
    for (int parameterDof = 0; parameterDof < parameters.Size(); ++parameterDof) {
        parameters(parameterDof) = 0.015 + 0.001 * static_cast<double>(parameterDof % 7);
    }

    mfem::Vector displacement(prepared.surfaceDisplacementSize());
    prepared.buildSurfaceDisplacement(parameters, displacement);

    const mfem::Vector referencePositions = nodal_radial_surface_test_utils::projectReferenceSurfacePositions(
        *fem.surfaceDeformationFes, prepared.surfaceDofMap()
    );

    constexpr double tolerance = 2.0e-13;

    for (int parameterDof = 0; parameterDof < prepared.parameterCount(); ++parameterDof) {
        double radiusSquared     = 0.0;
        double displacementNorm2 = 0.0;
        double radialProjection  = 0.0;

        for (int component = 0; component < prepared.spatialDimension(); ++component) {
            const int surfaceDof          = prepared.surfaceDisplacementDof(parameterDof, component);
            const double radialCoordinate = referencePositions(surfaceDof) - prepared.referenceCenter()(component);

            radiusSquared += radialCoordinate * radialCoordinate;
        }

        const double radius = std::sqrt(radiusSquared);
        REQUIRE(radius > 0.0);
        CHECK(std::abs(prepared.referenceRadius(parameterDof) - radius) <= tolerance * radius);

        for (int component = 0; component < prepared.spatialDimension(); ++component) {
            const int surfaceDof              = prepared.surfaceDisplacementDof(parameterDof, component);
            const double radialCoordinate     = referencePositions(surfaceDof) - prepared.referenceCenter()(component);
            const double expectedDirection    = radialCoordinate / radius;
            const double expectedDisplacement = parameters(parameterDof) * expectedDirection;

            CHECK(std::abs(prepared.radialDirection(parameterDof, component) - expectedDirection) <= tolerance);
            CHECK(std::abs(displacement(surfaceDof) - expectedDisplacement) <= tolerance);

            displacementNorm2 += displacement(surfaceDof) * displacement(surfaceDof);
            radialProjection += displacement(surfaceDof) * expectedDirection;
        }

        CHECK(std::abs(std::sqrt(displacementNorm2) - parameters(parameterDof)) <= tolerance);
        CHECK(std::abs(radialProjection - parameters(parameterDof)) <= tolerance);

        double movedRadiusSquared = 0.0;
        for (int component = 0; component < prepared.spatialDimension(); ++component) {
            const int surfaceDof = prepared.surfaceDisplacementDof(parameterDof, component);
            const double movedCoordinate =
                referencePositions(surfaceDof) + displacement(surfaceDof) - prepared.referenceCenter()(component);
            movedRadiusSquared += movedCoordinate * movedCoordinate;
        }
        CHECK(std::abs(std::sqrt(movedRadiusSquared) - (radius + parameters(parameterDof))) <= tolerance);
    }
}

TEST_CASE(
    "Nodal Radial Surface Jacobian Matches Centered Difference And Its Transpose Preserves Virtual Work",
    tags::nodal_radial_surface_linearization
) {
    namespace deformation        = mean_field::deformation;

    mean_field::utils::Args args = test_utils::setup_args();
    mean_field::fem::FEM fem     = mean_field::fem::setup_fem(args.mesh_file, args, 0);
    REQUIRE(fem.okay());

    const deformation::PreparedNodalRadialSurface prepared = nodal_radial_surface_test_utils::makePreparedSurface(fem);

    mfem::Vector parameters(prepared.parameterCount());
    mfem::Vector direction(prepared.parameterCount());
    for (int parameterDof = 0; parameterDof < prepared.parameterCount(); ++parameterDof) {
        const double index       = static_cast<double>(parameterDof + 1);
        parameters(parameterDof) = 0.013 * std::sin(0.37 * index);
        direction(parameterDof)  = std::cos(0.19 * index) - 0.21 * std::sin(0.43 * index);
    }

    constexpr double step = 1.0e-6;
    mfem::Vector plusParameters(parameters);
    mfem::Vector minusParameters(parameters);
    plusParameters.Add(step, direction);
    minusParameters.Add(-step, direction);

    mfem::Vector plusDisplacement(prepared.surfaceDisplacementSize());
    mfem::Vector minusDisplacement(prepared.surfaceDisplacementSize());
    mfem::Vector jacobianAction(prepared.surfaceDisplacementSize());
    prepared.buildSurfaceDisplacement(plusParameters, plusDisplacement);
    prepared.buildSurfaceDisplacement(minusParameters, minusDisplacement);
    prepared.applyJacobian(parameters, direction, jacobianAction);

    mfem::Vector centeredDifference(plusDisplacement);
    centeredDifference -= minusDisplacement;
    centeredDifference /= 2.0 * step;

    CHECK(nodal_radial_surface_test_utils::relativeError(jacobianAction, centeredDifference) < 2.0e-11);

    mfem::Vector surfaceDual(prepared.surfaceDisplacementSize());
    for (int surfaceDof = 0; surfaceDof < surfaceDual.Size(); ++surfaceDof) {
        const double index      = static_cast<double>(surfaceDof + 1);
        surfaceDual(surfaceDof) = std::sin(0.23 * index) + 0.17 * std::cos(0.31 * index);
    }

    mfem::Vector parameterDual(prepared.parameterCount());
    prepared.applyJacobianTranspose(parameters, surfaceDual, parameterDual);

    const double surfaceWork   = jacobianAction * surfaceDual;
    const double parameterWork = direction * parameterDual;
    const double workScale     = std::max({1.0, std::abs(surfaceWork), std::abs(parameterWork)});
    CHECK(std::abs(surfaceWork - parameterWork) <= 3.0e-14 * workScale);

    mfem::Vector pullbackDerivative(prepared.parameterCount());
    pullbackDerivative = 1.0;
    prepared.applyPullbackDerivative(parameters, direction, surfaceDual, pullbackDerivative);
    CHECK(pullbackDerivative.Norml2() == 0.0);

    mfem::Vector wrongParameters(prepared.parameterCount() + 1);
    mfem::Vector wrongSurfaceDisplacement(prepared.surfaceDisplacementSize() + 1);
    CHECK_THROWS_AS(prepared.buildSurfaceDisplacement(wrongParameters, plusDisplacement), std::invalid_argument);
    CHECK_THROWS_AS(prepared.buildSurfaceDisplacement(parameters, wrongSurfaceDisplacement), std::invalid_argument);
}