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

#include <catch2/catch_test_macros.hpp>

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

import mean_field;
import test_helpers;

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

    using Schema          = domain::CoreEnvelopeVacuumDomainSchema;

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

    [[nodiscard]] auto makePreparedDomainDeformation(mean_field::fem::FEM &fem) {
        const field::ScalarBoundaryDofMap surfaceDofMap =
            field::make_stellar_surface_scalar_dof_map<Schema>(*fem.surfaceDeformationFes);
        const deformation::SurfaceDeformationCompilationContext surfaceContext{
            *fem.surfaceDeformationFes, surfaceDofMap
        };
        deformation::PreparedNodalRadialSurface surface = deformation::compileSurfaceDeformationPrescription(
            deformation::NodalRadialSurface{referenceCenter(fem.mesh->SpaceDimension())}, surfaceContext
        );

        const deformation::RadialDeformationExtensionCompilationContext extensionContext =
            deformation::makeRadialDeformationExtensionCompilationContext<Schema>(
                *fem.surfaceDeformationFes, *fem.displacementFes, *fem.logicalReferenceMesh
            );
        deformation::PreparedPowerLawRadialInteriorExtension interior =
            deformation::compileInteriorDeformationExtension(
                deformation::PowerLawRadialInteriorExtension{}, extensionContext
            );
        deformation::PreparedFixedInfinityRadialVacuumExtension vacuum = deformation::compileVacuumDeformationExtension(
            deformation::FixedInfinityRadialVacuumExtension{}, extensionContext
        );

        return deformation::composePreparedDomainDeformation(
            std::move(surface), std::move(interior), std::move(vacuum), *fem.surfaceDeformationFes,
            *fem.displacementFes, *fem.logicalReferenceMesh
        );
    }

    [[nodiscard]] int volumeVectorDof(
        const int scalarTrueDof,
        const int component,
        const int scalarTrueDofCount
    ) {
        return scalarTrueDof + component * scalarTrueDofCount;
    }

    [[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 globalInnerProduct(
        const mfem::Vector &first,
        const mfem::Vector &second,
        MPI_Comm communicator
    ) {
        REQUIRE(first.Size() == second.Size());
        const double localValue = first * second;
        double globalValue      = 0.0;
        MPI_Allreduce(&localValue, &globalValue, 1, MPI_DOUBLE, MPI_SUM, communicator);
        return globalValue;
    }

    class AnalyticNonlinearSurface final {
    public:
        [[nodiscard]] deformation::SurfaceDeformationDescriptor descriptor() const noexcept {
            return {
                .name                               = "AnalyticNonlinearSurface",
                .spatialDimension                   = 3,
                .motionKind                         = deformation::SurfaceMotionKind::Radial,
                .linearOnReferenceGeometry          = false,
                .requiresStarShapedReferenceSurface = false,
                .hasExactDerivativeTranspose        = true,
                .hasExactPullbackDerivative         = true,
                .translationTreatment               = deformation::GeometricGaugeTreatment::Retained,
                .orientationTreatment               = deformation::GeometricGaugeTreatment::Retained
            };
        }

        [[nodiscard]] int parameterCount() const noexcept {
            return 2;
        }

        [[nodiscard]] int surfaceDisplacementSize() const noexcept {
            return 2;
        }

        void buildSurfaceDisplacement(
            const mfem::Vector &parameters,
            mfem::Vector &surfaceDisplacement
        ) const {
            surfaceDisplacement(0) = parameters(0) * parameters(0) + parameters(1);
            surfaceDisplacement(1) = parameters(0) * parameters(1);
        }

        void applyJacobian(
            const mfem::Vector &parameters,
            const mfem::Vector &parameterDirection,
            mfem::Vector &surfaceDisplacementDirection
        ) const {
            surfaceDisplacementDirection(0) = 2.0 * parameters(0) * parameterDirection(0) + parameterDirection(1);
            surfaceDisplacementDirection(1) =
                parameters(1) * parameterDirection(0) + parameters(0) * parameterDirection(1);
        }

        void applyJacobianTranspose(
            const mfem::Vector &parameters,
            const mfem::Vector &surfaceDisplacementDual,
            mfem::Vector &parameterDual
        ) const {
            parameterDual(0) =
                2.0 * parameters(0) * surfaceDisplacementDual(0) + parameters(1) * surfaceDisplacementDual(1);
            parameterDual(1) = surfaceDisplacementDual(0) + parameters(0) * surfaceDisplacementDual(1);
        }

        void applyPullbackDerivative(
            const mfem::Vector &,
            const mfem::Vector &parameterDirection,
            const mfem::Vector &surfaceDisplacementDual,
            mfem::Vector &parameterDualAction
        ) const {
            parameterDualAction(0) = 2.0 * parameterDirection(0) * surfaceDisplacementDual(0) +
                                     parameterDirection(1) * surfaceDisplacementDual(1);
            parameterDualAction(1) = parameterDirection(0) * surfaceDisplacementDual(1);
        }
    };

    class AnalyticNonlinearExtensionKernel {
    public:
        AnalyticNonlinearExtensionKernel(
            const int scalarTrueDofCount,
            const double coefficientScale
        )
            : m_scalarTrueDofCount(scalarTrueDofCount),
              m_coefficientScale(coefficientScale) {
        }

        [[nodiscard]] int scalarTrueDofCount() const noexcept {
            return m_scalarTrueDofCount;
        }

        [[nodiscard]] int volumeDisplacementSize() const noexcept {
            return 3 * m_scalarTrueDofCount;
        }

        [[nodiscard]] int sharedScalarDof() const noexcept {
            return m_scalarTrueDofCount / 2;
        }

        void build(
            const mfem::Vector &surfaceDisplacement,
            mfem::Vector &volumeDisplacement,
            const bool stellar
        ) const {
            volumeDisplacement = 0.0;
            for (int vectorDof = 0; vectorDof < volumeDisplacementSize(); ++vectorDof) {
                if (!hasSupport(vectorDof % m_scalarTrueDofCount, stellar)) {
                    continue;
                }
                double linearFirst  = 0.0;
                double linearSecond = 0.0;
                double bilinear     = 0.0;
                coefficients(vectorDof, linearFirst, linearSecond, bilinear);
                volumeDisplacement(vectorDof) = linearFirst * surfaceDisplacement(0) +
                                                linearSecond * surfaceDisplacement(1) +
                                                bilinear * surfaceDisplacement(0) * surfaceDisplacement(1);
            }
        }

        void applyJacobian(
            const mfem::Vector &surfaceDisplacement,
            const mfem::Vector &surfaceDirection,
            mfem::Vector &volumeDirection,
            const bool stellar
        ) const {
            volumeDirection = 0.0;
            for (int vectorDof = 0; vectorDof < volumeDisplacementSize(); ++vectorDof) {
                if (!hasSupport(vectorDof % m_scalarTrueDofCount, stellar)) {
                    continue;
                }
                double linearFirst  = 0.0;
                double linearSecond = 0.0;
                double bilinear     = 0.0;
                coefficients(vectorDof, linearFirst, linearSecond, bilinear);
                volumeDirection(vectorDof) = (linearFirst + bilinear * surfaceDisplacement(1)) * surfaceDirection(0) +
                                             (linearSecond + bilinear * surfaceDisplacement(0)) * surfaceDirection(1);
            }
        }

        void applyJacobianTranspose(
            const mfem::Vector &surfaceDisplacement,
            const mfem::Vector &volumeDual,
            mfem::Vector &surfaceDual,
            const bool stellar
        ) const {
            surfaceDual = 0.0;
            for (int vectorDof = 0; vectorDof < volumeDisplacementSize(); ++vectorDof) {
                if (!hasSupport(vectorDof % m_scalarTrueDofCount, stellar)) {
                    continue;
                }
                double linearFirst  = 0.0;
                double linearSecond = 0.0;
                double bilinear     = 0.0;
                coefficients(vectorDof, linearFirst, linearSecond, bilinear);
                surfaceDual(0) += (linearFirst + bilinear * surfaceDisplacement(1)) * volumeDual(vectorDof);
                surfaceDual(1) += (linearSecond + bilinear * surfaceDisplacement(0)) * volumeDual(vectorDof);
            }
        }

        void applyPullbackDerivative(
            const mfem::Vector &surfaceDirection,
            const mfem::Vector &volumeDual,
            mfem::Vector &surfaceDualAction,
            const bool stellar
        ) const {
            surfaceDualAction = 0.0;
            for (int vectorDof = 0; vectorDof < volumeDisplacementSize(); ++vectorDof) {
                if (!hasSupport(vectorDof % m_scalarTrueDofCount, stellar)) {
                    continue;
                }
                double linearFirst  = 0.0;
                double linearSecond = 0.0;
                double bilinear     = 0.0;
                coefficients(vectorDof, linearFirst, linearSecond, bilinear);
                surfaceDualAction(0) += bilinear * surfaceDirection(1) * volumeDual(vectorDof);
                surfaceDualAction(1) += bilinear * surfaceDirection(0) * volumeDual(vectorDof);
            }
        }

        [[nodiscard]] bool hasSupport(
            const int scalarTrueDof,
            const bool stellar
        ) const noexcept {
            return stellar ? scalarTrueDof <= sharedScalarDof() : scalarTrueDof >= sharedScalarDof();
        }

    private:
        void coefficients(
            const int vectorDof,
            double &linearFirst,
            double &linearSecond,
            double &bilinear
        ) const noexcept {
            linearFirst  = m_coefficientScale * (0.01 + 0.001 * static_cast<double>(vectorDof % 7));
            linearSecond = m_coefficientScale * (-0.02 + 0.002 * static_cast<double>(vectorDof % 5));
            bilinear     = m_coefficientScale * 0.0005 * static_cast<double>(1 + vectorDof % 3);
        }

        int m_scalarTrueDofCount;
        double m_coefficientScale;
    };

    class AnalyticNonlinearInteriorExtension final {
    public:
        explicit AnalyticNonlinearInteriorExtension(
            const int scalarTrueDofCount,
            const bool supportEnabled = true
        )
            : m_kernel(
                  scalarTrueDofCount,
                  1.0
              ),
              m_supportEnabled(supportEnabled) {
        }

        [[nodiscard]] deformation::InteriorDeformationExtensionDescriptor descriptor() const noexcept {
            return {
                .name                        = "AnalyticNonlinearInteriorExtension",
                .spatialDimension            = 3,
                .linearOnReferenceGeometry   = false,
                .requiresRadialFoliation     = false,
                .requiresAuxiliarySolve      = false,
                .hasExactDerivativeTranspose = true,
                .hasExactPullbackDerivative  = true,
                .centerBehavior              = deformation::InteriorCenterBehavior::FixedAtReferenceCenter
            };
        }

        [[nodiscard]] int surfaceDisplacementSize() const noexcept {
            return 2;
        }
        [[nodiscard]] int interiorDisplacementSize() const noexcept {
            return m_kernel.volumeDisplacementSize();
        }
        [[nodiscard]] int scalarTrueDofCount() const noexcept {
            return m_kernel.scalarTrueDofCount();
        }
        [[nodiscard]] bool hasStellarSupport(const int scalarTrueDof) const {
            return m_supportEnabled && m_kernel.hasSupport(scalarTrueDof, true);
        }
        void buildInteriorDisplacement(
            const mfem::Vector &surface,
            mfem::Vector &volume
        ) const {
            m_kernel.build(surface, volume, true);
        }
        void applyJacobian(
            const mfem::Vector &surface,
            const mfem::Vector &direction,
            mfem::Vector &volume
        ) const {
            m_kernel.applyJacobian(surface, direction, volume, true);
        }
        void applyJacobianTranspose(
            const mfem::Vector &surface,
            const mfem::Vector &volume,
            mfem::Vector &dual
        ) const {
            m_kernel.applyJacobianTranspose(surface, volume, dual, true);
        }
        void applyPullbackDerivative(
            const mfem::Vector &,
            const mfem::Vector &direction,
            const mfem::Vector &volume,
            mfem::Vector &dual
        ) const {
            m_kernel.applyPullbackDerivative(direction, volume, dual, true);
        }

    private:
        AnalyticNonlinearExtensionKernel m_kernel;
        bool m_supportEnabled;
    };

    class AnalyticNonlinearVacuumExtension final {
    public:
        explicit AnalyticNonlinearVacuumExtension(const int scalarTrueDofCount)
            : m_kernel(
                  scalarTrueDofCount,
                  -0.7
              ) {
        }

        [[nodiscard]] deformation::VacuumDeformationExtensionDescriptor descriptor() const noexcept {
            return {
                .name                        = "AnalyticNonlinearVacuumExtension",
                .spatialDimension            = 3,
                .linearOnReferenceGeometry   = false,
                .requiresRadialFoliation     = false,
                .requiresAuxiliarySolve      = false,
                .hasExactDerivativeTranspose = true,
                .hasExactPullbackDerivative  = true,
                .outerBoundaryBehavior       = deformation::VacuumOuterBoundaryBehavior::FixedAtReferenceInfinity
            };
        }

        [[nodiscard]] int surfaceDisplacementSize() const noexcept {
            return 2;
        }
        [[nodiscard]] int vacuumDisplacementSize() const noexcept {
            return m_kernel.volumeDisplacementSize();
        }
        [[nodiscard]] int scalarTrueDofCount() const noexcept {
            return m_kernel.scalarTrueDofCount();
        }
        [[nodiscard]] bool hasVacuumSupport(const int scalarTrueDof) const {
            return m_kernel.hasSupport(scalarTrueDof, false);
        }
        void buildVacuumDisplacement(
            const mfem::Vector &surface,
            mfem::Vector &volume
        ) const {
            m_kernel.build(surface, volume, false);
        }
        void applyJacobian(
            const mfem::Vector &surface,
            const mfem::Vector &direction,
            mfem::Vector &volume
        ) const {
            m_kernel.applyJacobian(surface, direction, volume, false);
        }
        void applyJacobianTranspose(
            const mfem::Vector &surface,
            const mfem::Vector &volume,
            mfem::Vector &dual
        ) const {
            m_kernel.applyJacobianTranspose(surface, volume, dual, false);
        }
        void applyPullbackDerivative(
            const mfem::Vector &,
            const mfem::Vector &direction,
            const mfem::Vector &volume,
            mfem::Vector &dual
        ) const {
            m_kernel.applyPullbackDerivative(direction, volume, dual, false);
        }

    private:
        AnalyticNonlinearExtensionKernel m_kernel;
    };
} // namespace domain_deformation_test_utils

TEST_CASE(
    "Prepared Domain Deformation Composes Surface Interior And Vacuum Maps With Explicit Ownership",
    tags::domain_deformation_composition
) {
    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());

    auto prepared = domain_deformation_test_utils::makePreparedDomainDeformation(fem);
    STATIC_CHECK(deformation::PreparedDomainDeformationOperator<decltype(prepared)>);

    const deformation::DomainDeformationDescriptor descriptor = prepared.descriptor();
    REQUIRE(descriptor.isValid());
    CHECK(descriptor.linearOnReferenceGeometry);
    CHECK_FALSE(descriptor.requiresAuxiliarySolve);
    CHECK(descriptor.supportsExactNewtonLinearization());
    CHECK(prepared.parameterCount() == prepared.surfaceDeformationPrescription().parameterCount());
    CHECK(prepared.surfaceDisplacementSize() == prepared.surfaceDeformationPrescription().surfaceDisplacementSize());
    CHECK(prepared.volumeDisplacementSize() == fem.displacementFes->GetTrueVSize());
    CHECK(prepared.matchesCurrentDiscretization());

    const deformation::DomainDeformationDiscretizationDependencies &dependencies =
        prepared.discretizationDependencies();
    CHECK(dependencies.physicalMeshIdentity == fem.mesh.get());
    CHECK(dependencies.logicalReferenceMeshIdentity == fem.logicalReferenceMesh.get());
    CHECK(dependencies.surfaceScalarSpaceIdentity == fem.surfaceDeformationFes.get());
    CHECK(dependencies.volumeDisplacementSpaceIdentity == fem.displacementFes.get());
    CHECK(dependencies.isCurrent());

    const deformation::DomainDeformationCompositionReport &composition = prepared.compositionReport();
    CHECK(composition.scalarTrueDofCount == prepared.scalarTrueDofCount());
    CHECK(composition.assignedScalarDofCount() == prepared.scalarTrueDofCount());
    CHECK(composition.stellarInteriorOwnedScalarDofCount > 0);
    CHECK(composition.vacuumOwnedScalarDofCount > 0);
    CHECK(composition.sharedSurfaceScalarDofCount > 0);

    int countedStellarOwners = 0;
    int countedVacuumOwners  = 0;
    int countedSharedDofs    = 0;
    for (int scalarTrueDof = 0; scalarTrueDof < prepared.scalarTrueDofCount(); ++scalarTrueDof) {
        if (prepared.volumeOwner(scalarTrueDof) == deformation::VolumeDeformationOwner::StellarInterior) {
            ++countedStellarOwners;
        } else {
            ++countedVacuumOwners;
        }
        countedSharedDofs += prepared.isSharedSurfaceDof(scalarTrueDof) ? 1 : 0;
    }
    CHECK(countedStellarOwners == composition.stellarInteriorOwnedScalarDofCount);
    CHECK(countedVacuumOwners == composition.vacuumOwnedScalarDofCount);
    CHECK(countedSharedDofs == composition.sharedSurfaceScalarDofCount);

    mfem::Vector zeroParameters(prepared.parameterCount());
    zeroParameters = 0.0;
    mfem::Vector composedVolume(prepared.volumeDisplacementSize());
    prepared.buildVolumeDisplacement(zeroParameters, composedVolume);
    CHECK(composedVolume.Norml2() == 0.0);

    mfem::Vector parameters(prepared.parameterCount());
    for (int parameter = 0; parameter < parameters.Size(); ++parameter) {
        const double index    = static_cast<double>(parameter + 1);
        parameters(parameter) = 0.012 * std::sin(0.19 * index) - 0.004 * std::cos(0.31 * index);
    }
    prepared.buildVolumeDisplacement(parameters, composedVolume);

    mfem::Vector surfaceDisplacement(prepared.surfaceDisplacementSize());
    mfem::Vector interiorVolume(prepared.volumeDisplacementSize());
    mfem::Vector vacuumVolume(prepared.volumeDisplacementSize());
    prepared.surfaceDeformationPrescription().buildSurfaceDisplacement(parameters, surfaceDisplacement);
    prepared.stellarInteriorExtension().buildInteriorDisplacement(surfaceDisplacement, interiorVolume);
    prepared.vacuumExtension().buildVacuumDisplacement(surfaceDisplacement, vacuumVolume);

    constexpr double tolerance = 2.0e-12;
    for (int scalarTrueDof = 0; scalarTrueDof < prepared.scalarTrueDofCount(); ++scalarTrueDof) {
        for (int component = 0; component < prepared.spatialDimension(); ++component) {
            const int volumeDof =
                domain_deformation_test_utils::volumeVectorDof(scalarTrueDof, component, prepared.scalarTrueDofCount());
            const double expected =
                prepared.volumeOwner(scalarTrueDof) == deformation::VolumeDeformationOwner::StellarInterior
                    ? interiorVolume(volumeDof)
                    : vacuumVolume(volumeDof);
            CHECK(std::abs(composedVolume(volumeDof) - expected) <= tolerance);
            if (prepared.isSharedSurfaceDof(scalarTrueDof)) {
                CHECK(std::abs(interiorVolume(volumeDof) - vacuumVolume(volumeDof)) <= tolerance);
            }
        }
    }

    const deformation::PreparedDomainDeformationActionStatistics &statistics = prepared.actionStatistics();
    CHECK(statistics.volumeBuildApplications == 2);
    CHECK(statistics.jacobianApplications == 0);
    CHECK(statistics.jacobianTransposeApplications == 0);
    CHECK(statistics.pullbackDerivativeApplications == 0);
    CHECK(statistics.geometryInspections == 0);

    mfem::Vector wrongParameters(prepared.parameterCount() + 1);
    mfem::Vector wrongVolume(prepared.volumeDisplacementSize() + 1);
    CHECK_THROWS_AS(prepared.buildVolumeDisplacement(wrongParameters, composedVolume), std::invalid_argument);
    CHECK_THROWS_AS(prepared.buildVolumeDisplacement(parameters, wrongVolume), std::invalid_argument);
    CHECK_THROWS_AS(prepared.volumeOwner(-1), std::out_of_range);
}

TEST_CASE(
    "Prepared Domain Deformation Jacobian Matches Centered Difference And Its Pullback Preserves Virtual Work",
    tags::domain_deformation_linearization
) {
    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());

    auto prepared = domain_deformation_test_utils::makePreparedDomainDeformation(fem);
    mfem::Vector parameters(prepared.parameterCount());
    mfem::Vector direction(prepared.parameterCount());
    for (int parameter = 0; parameter < parameters.Size(); ++parameter) {
        const double index    = static_cast<double>(parameter + 1);
        parameters(parameter) = 0.008 * std::sin(0.13 * index);
        direction(parameter)  = std::cos(0.17 * index) - 0.25 * std::sin(0.29 * 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 plusVolume(prepared.volumeDisplacementSize());
    mfem::Vector minusVolume(prepared.volumeDisplacementSize());
    mfem::Vector jacobianAction(prepared.volumeDisplacementSize());
    prepared.buildVolumeDisplacement(plusParameters, plusVolume);
    prepared.buildVolumeDisplacement(minusParameters, minusVolume);
    prepared.applyJacobian(parameters, direction, jacobianAction);

    mfem::Vector centeredDifference(plusVolume);
    centeredDifference -= minusVolume;
    centeredDifference /= 2.0 * step;
    CHECK(domain_deformation_test_utils::relativeError(jacobianAction, centeredDifference) < 3.0e-10);

    mfem::Vector volumeDual(prepared.volumeDisplacementSize());
    for (int dof = 0; dof < volumeDual.Size(); ++dof) {
        const double index = static_cast<double>(dof + 1);
        volumeDual(dof)    = std::sin(0.07 * index) + 0.4 * std::cos(0.11 * index);
    }
    mfem::Vector parameterDual(prepared.parameterCount());
    prepared.applyJacobianTranspose(parameters, volumeDual, parameterDual);

    const double volumeWork =
        domain_deformation_test_utils::globalInnerProduct(jacobianAction, volumeDual, fem.mesh->GetComm());
    const double parameterWork =
        domain_deformation_test_utils::globalInnerProduct(direction, parameterDual, fem.mesh->GetComm());
    const double workScale = std::max({1.0, std::abs(volumeWork), std::abs(parameterWork)});
    CHECK(std::abs(volumeWork - parameterWork) <= 8.0e-13 * workScale);

    mfem::Vector pullbackAction(prepared.parameterCount());
    prepared.applyPullbackDerivative(parameters, direction, volumeDual, pullbackAction);
    CHECK(pullbackAction.Norml2() == 0.0);

    mfem::Vector plusTranspose(prepared.parameterCount());
    mfem::Vector minusTranspose(prepared.parameterCount());
    prepared.applyJacobianTranspose(plusParameters, volumeDual, plusTranspose);
    prepared.applyJacobianTranspose(minusParameters, volumeDual, minusTranspose);
    mfem::Vector transposeCenteredDifference(plusTranspose);
    transposeCenteredDifference -= minusTranspose;
    transposeCenteredDifference /= 2.0 * step;
    CHECK(transposeCenteredDifference.Norml2() <= 1.0e-12 * std::max(1.0, parameterDual.Norml2()));
}

TEST_CASE(
    "Nonlinear Domain Deformation Pullback Matches The Directional Derivative Of Its Complete Transpose",
    tags::domain_deformation_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 int scalarTrueDofCount = fem.displacementFes->GetTrueVSize() / fem.mesh->SpaceDimension();
    CHECK_THROWS_AS(
        (deformation::composePreparedDomainDeformation(
            domain_deformation_test_utils::AnalyticNonlinearSurface{},
            domain_deformation_test_utils::AnalyticNonlinearInteriorExtension{scalarTrueDofCount, false},
            domain_deformation_test_utils::AnalyticNonlinearVacuumExtension{scalarTrueDofCount},
            *fem.surfaceDeformationFes, *fem.displacementFes, *fem.logicalReferenceMesh
        )),
        std::invalid_argument
    );
    CHECK_THROWS_AS(
        (deformation::composePreparedDomainDeformation(
            domain_deformation_test_utils::AnalyticNonlinearSurface{},
            domain_deformation_test_utils::AnalyticNonlinearInteriorExtension{scalarTrueDofCount},
            domain_deformation_test_utils::AnalyticNonlinearVacuumExtension{scalarTrueDofCount - 1},
            *fem.surfaceDeformationFes, *fem.displacementFes, *fem.logicalReferenceMesh
        )),
        std::invalid_argument
    );
    auto prepared = deformation::composePreparedDomainDeformation(
        domain_deformation_test_utils::AnalyticNonlinearSurface{},
        domain_deformation_test_utils::AnalyticNonlinearInteriorExtension{scalarTrueDofCount},
        domain_deformation_test_utils::AnalyticNonlinearVacuumExtension{scalarTrueDofCount}, *fem.surfaceDeformationFes,
        *fem.displacementFes, *fem.logicalReferenceMesh
    );
    REQUIRE_FALSE(prepared.descriptor().linearOnReferenceGeometry);

    mfem::Vector parameters(2);
    parameters(0) = 0.37;
    parameters(1) = -0.21;
    mfem::Vector parameterDirection(2);
    parameterDirection(0) = -0.42;
    parameterDirection(1) = 0.63;

    mfem::Vector volumeDual(prepared.volumeDisplacementSize());
    for (int vectorDof = 0; vectorDof < volumeDual.Size(); ++vectorDof) {
        const double index    = static_cast<double>(vectorDof + 1);
        volumeDual(vectorDof) = std::sin(0.013 * index) - 0.3 * std::cos(0.021 * index);
    }

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

    mfem::Vector plusVolume(prepared.volumeDisplacementSize());
    mfem::Vector minusVolume(prepared.volumeDisplacementSize());
    mfem::Vector jacobianAction(prepared.volumeDisplacementSize());
    prepared.buildVolumeDisplacement(plusParameters, plusVolume);
    prepared.buildVolumeDisplacement(minusParameters, minusVolume);
    prepared.applyJacobian(parameters, parameterDirection, jacobianAction);
    mfem::Vector centeredJacobian(plusVolume);
    centeredJacobian -= minusVolume;
    centeredJacobian /= 2.0 * step;
    CHECK(domain_deformation_test_utils::relativeError(jacobianAction, centeredJacobian) < 2.0e-10);

    mfem::Vector parameterDual(2);
    prepared.applyJacobianTranspose(parameters, volumeDual, parameterDual);
    const double volumeWork =
        domain_deformation_test_utils::globalInnerProduct(jacobianAction, volumeDual, fem.mesh->GetComm());
    const double parameterWork =
        domain_deformation_test_utils::globalInnerProduct(parameterDirection, parameterDual, fem.mesh->GetComm());
    CHECK(
        std::abs(volumeWork - parameterWork) <= 2.0e-12 * std::max({1.0, std::abs(volumeWork), std::abs(parameterWork)})
    );

    mfem::Vector pullbackAction(2);
    prepared.applyPullbackDerivative(parameters, parameterDirection, volumeDual, pullbackAction);
    mfem::Vector plusTranspose(2);
    mfem::Vector minusTranspose(2);
    prepared.applyJacobianTranspose(plusParameters, volumeDual, plusTranspose);
    prepared.applyJacobianTranspose(minusParameters, volumeDual, minusTranspose);
    mfem::Vector centeredPullback(plusTranspose);
    centeredPullback -= minusTranspose;
    centeredPullback /= 2.0 * step;
    CHECK(domain_deformation_test_utils::relativeError(pullbackAction, centeredPullback) < 3.0e-9);

    const int sharedScalarDof = scalarTrueDofCount / 2;
    REQUIRE(prepared.isSharedSurfaceDof(sharedScalarDof));
    CHECK(prepared.volumeOwner(sharedScalarDof) == deformation::VolumeDeformationOwner::StellarInterior);

    mfem::Vector surfaceDisplacement(2);
    mfem::Vector interiorVolume(prepared.volumeDisplacementSize());
    mfem::Vector vacuumVolume(prepared.volumeDisplacementSize());
    mfem::Vector composedVolume(prepared.volumeDisplacementSize());
    prepared.surfaceDeformationPrescription().buildSurfaceDisplacement(parameters, surfaceDisplacement);
    prepared.stellarInteriorExtension().buildInteriorDisplacement(surfaceDisplacement, interiorVolume);
    prepared.vacuumExtension().buildVacuumDisplacement(surfaceDisplacement, vacuumVolume);
    prepared.buildVolumeDisplacement(parameters, composedVolume);
    for (int component = 0; component < prepared.spatialDimension(); ++component) {
        const int vectorDof = sharedScalarDof + component * scalarTrueDofCount;
        CHECK(composedVolume(vectorDof) == interiorVolume(vectorDof));
        CHECK(interiorVolume(vectorDof) != vacuumVolume(vectorDof));
    }
}

TEST_CASE(
    "Prepared Domain Deformation Accepts Orientation Preserving Shapes And Rejects Folded Volume Maps",
    tags::domain_deformation_geometry
) {
    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());

    auto prepared = domain_deformation_test_utils::makePreparedDomainDeformation(fem);
    mfem::Vector volumeDisplacement(prepared.volumeDisplacementSize());

    mfem::Vector zeroParameters(prepared.parameterCount());
    zeroParameters = 0.0;
    const mean_field::deformation::DomainDeformationGeometryReport referenceReport =
        prepared.buildValidatedVolumeDisplacement(zeroParameters, volumeDisplacement, 0.99);
    CHECK(std::abs(referenceReport.minimumJacobianDeterminant - 1.0) <= 64.0 * std::numeric_limits<double>::epsilon());
    CHECK_THROWS_AS(
        prepared.buildValidatedVolumeDisplacement(zeroParameters, volumeDisplacement, 1.01), std::domain_error
    );
    CHECK_THROWS_AS(
        prepared.buildValidatedVolumeDisplacement(zeroParameters, volumeDisplacement, -0.01), std::invalid_argument
    );
    CHECK_THROWS_AS(
        prepared.buildValidatedVolumeDisplacement(
            zeroParameters, volumeDisplacement, std::numeric_limits<double>::quiet_NaN()
        ),
        std::invalid_argument
    );

    mfem::Vector boundedParameters(prepared.parameterCount());
    boundedParameters = 0.02;
    const mean_field::deformation::DomainDeformationGeometryReport boundedReport =
        prepared.buildValidatedVolumeDisplacement(boundedParameters, volumeDisplacement);
    CHECK(boundedReport.isOrientationPreserving());
    CHECK(boundedReport.minimumJacobianDeterminant > 0.0);

    mfem::Vector foldingParameters(prepared.parameterCount());
    foldingParameters = -2.0;
    prepared.buildVolumeDisplacement(foldingParameters, volumeDisplacement);
    const mean_field::deformation::DomainDeformationGeometryReport foldingReport =
        prepared.inspectMappedGeometry(volumeDisplacement);
    CHECK_FALSE(foldingReport.isOrientationPreserving());
    CHECK_THROWS_AS(
        prepared.buildValidatedVolumeDisplacement(foldingParameters, volumeDisplacement), std::domain_error
    );

    mfem::Vector nonfiniteParameters(zeroParameters);
    nonfiniteParameters(0) = std::numeric_limits<double>::quiet_NaN();
    CHECK_THROWS_AS(
        prepared.buildValidatedVolumeDisplacement(nonfiniteParameters, volumeDisplacement), std::domain_error
    );

    CHECK(prepared.actionStatistics().geometryInspections == 6);
}

TEST_CASE(
    "Prepared Domain Deformation Rejects Actions After Its Discretization Becomes Stale",
    tags::domain_deformation_composition
) {
    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());

    auto prepared = domain_deformation_test_utils::makePreparedDomainDeformation(fem);
    mfem::Vector parameters(prepared.parameterCount());
    mfem::Vector parameterDirection(prepared.parameterCount());
    mfem::Vector parameterDual(prepared.parameterCount());
    mfem::Vector volume(prepared.volumeDisplacementSize());
    parameters         = 0.0;
    parameterDirection = 0.0;
    volume             = 0.0;

    REQUIRE(prepared.matchesCurrentDiscretization());
    fem.mesh->UniformRefinement();
    REQUIRE_FALSE(prepared.matchesCurrentDiscretization());
    CHECK_THROWS_AS(prepared.buildVolumeDisplacement(parameters, volume), std::logic_error);
    CHECK_THROWS_AS(prepared.applyJacobian(parameters, parameterDirection, volume), std::logic_error);
    CHECK_THROWS_AS(prepared.applyJacobianTranspose(parameters, volume, parameterDual), std::logic_error);
    CHECK_THROWS_AS(
        prepared.applyPullbackDerivative(parameters, parameterDirection, volume, parameterDual), std::logic_error
    );
    CHECK_THROWS_AS(prepared.inspectMappedGeometry(volume), std::logic_error);
}