#include <catch2/catch_test_macros.hpp>

#include <algorithm>
#include <array>
#include <concepts>
#include <cstddef>
#include <memory>
#include <mpi.h>
#include <set>
#include <stdexcept>
#include <string_view>
#include <vector>

#include <mfem.hpp>

import mean_field;
import test_helpers;

namespace field_mfem_test_utils {
    namespace field      = mean_field::field;
    namespace domain     = mean_field::utils::domain;
    namespace quadrature = mean_field::quadrature;

    using Schema         = domain::CoreEnvelopeVacuumDomainSchema;

    struct VectorL2Field {
        static constexpr std::string_view name = "test_vector_l2";

        using Support                          = field::DomainSupport<domain::All>;

        struct Vector final : field::VectorQ<field::FieldRelation::Independent, field::Disc<field::L2, 2>> { };

        using Quantities                          = field::TypeList<Vector>;

        using Constraints                         = field::TypeList<>;

        using FormList                            = field::TypeList<>;

        static constexpr bool constraintsAreValid = field::validate_constraints(Constraints{});

        static_assert(constraintsAreValid);
    };

    struct NdField {
        static constexpr std::string_view name = "test_nd";

        using Support                          = field::DomainSupport<domain::All>;

        struct Vector final : field::VectorQ<field::FieldRelation::Independent, field::Disc<field::ND, 2>> { };

        using Quantities                          = field::TypeList<Vector>;

        using Constraints                         = field::TypeList<>;

        using FormList                            = field::TypeList<>;

        static constexpr bool constraintsAreValid = field::validate_constraints(Constraints{});

        static_assert(constraintsAreValid);
    };

    using AlternateSchema = domain::DomainSchema<
        domain::MaterialList<
            domain::Material<domain::Core, 11>,
            domain::Material<domain::Envelope, 17>,
            domain::Material<domain::Vacuum, 29>>,
        domain::BoundaryList<>,
        domain::RelationList<>>;

    template <typename FieldT>
    concept CanResolveLocalSupport = requires(const mfem::FiniteElementSpace &space) {
        field::resolve_field_local_dof_support<FieldT, Schema>(space);
    };

    [[nodiscard]]
    mfem::Mesh make_two_domain_mesh(
        const int leftAttribute  = 2,
        const int rightAttribute = 3
    ) {
        mfem::Mesh mesh = mfem::Mesh::MakeCartesian2D(2, 1, mfem::Element::QUADRILATERAL, true, 2.0, 1.0);

        mesh.GetElement(0)->SetAttribute(leftAttribute);

        mesh.GetElement(1)->SetAttribute(rightAttribute);

        mesh.SetAttributes();

        return mesh;
    }

    [[nodiscard]]
    mfem::Mesh make_parallel_split_mesh() {
        constexpr int xElementCount = 4;
        constexpr int yElementCount = 2;

        mfem::Mesh mesh =
            mfem::Mesh::MakeCartesian2D(xElementCount, yElementCount, mfem::Element::QUADRILATERAL, true, 4.0, 2.0);

        for (int elementId = 0; elementId < mesh.GetNE(); ++elementId) {
            const int xIndex    = elementId % xElementCount;

            const int attribute = xIndex < 2 ? 2 : 3;

            mesh.GetElement(elementId)->SetAttribute(attribute);
        }

        mesh.SetAttributes();

        return mesh;
    }

    [[nodiscard]]
    std::vector<int> decoded_element_vdofs(
        const mfem::FiniteElementSpace &space,
        const int elementId
    ) {
        mfem::Array<int> signedVDofs;

        space.GetElementVDofs(elementId, signedVDofs);

        std::vector<int> result;

        result.reserve(static_cast<std::size_t>(signedVDofs.Size()));

        for (int index = 0; index < signedVDofs.Size(); ++index) {
            result.push_back(mfem::FiniteElementSpace::DecodeDof(signedVDofs[index]));
        }

        std::ranges::sort(result);

        result.erase(std::unique(result.begin(), result.end()), result.end());

        return result;
    }

    [[nodiscard]]
    bool contains(
        const mfem::Array<int> &values,
        const int value
    ) {
        for (int index = 0; index < values.Size(); ++index) {
            if (values[index] == value) {
                return true;
            }
        }

        return false;
    }

    [[nodiscard]]
    std::vector<int> intersection(
        const std::vector<int> &first,
        const std::vector<int> &second
    ) {
        std::vector<int> result;

        std::set_intersection(first.begin(), first.end(), second.begin(), second.end(), std::back_inserter(result));

        return result;
    }

    [[nodiscard]]
    std::vector<int> difference(
        const std::vector<int> &first,
        const std::vector<int> &second
    ) {
        std::vector<int> result;

        std::set_difference(first.begin(), first.end(), second.begin(), second.end(), std::back_inserter(result));

        return result;
    }

    [[nodiscard]]
    long long global_sum(const int localValue) {
        const long long local = static_cast<long long>(localValue);

        long long global      = 0;

        MPI_Allreduce(&local, &global, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD);

        return global;
    }
} // namespace field_mfem_test_utils

TEST_CASE(
    "Field MFEM Support Resolution Is Available Only For Domain Supported Fields",
    tags::unit &tags::field
) {
    namespace field = mean_field::field;

    STATIC_REQUIRE(field::MfemDomainField<field::Density>);

    STATIC_REQUIRE(field::MfemDomainField<field::Enthalpy>);

    STATIC_REQUIRE(field::MfemDomainField<field::Gravity>);

    STATIC_REQUIRE(field::MfemDomainField<field::Displacement>);

    STATIC_REQUIRE_FALSE(field::MfemDomainField<field::BarotropicConstant>);

    STATIC_REQUIRE(field_mfem_test_utils::CanResolveLocalSupport<field::Density>);

    STATIC_REQUIRE_FALSE(field_mfem_test_utils::CanResolveLocalSupport<field::BarotropicConstant>);

    CHECK(true);
}

TEST_CASE(
    "Field MFEM Creates The Registered Finite Element Collection Families",
    tags::unit &tags::field
) {
    namespace field             = mean_field::field;

    using DensityField          = field::Field<field::Density>;

    using GravityField          = field::Field<field::Gravity>;

    using DisplacementField     = field::Field<field::Displacement>;

    using EnthalpyField         = field::Field<field::Enthalpy>;

    auto densityCollection      = DensityField::make_fec<field::Density::Scalar>(3);

    auto potentialCollection    = GravityField::make_fec<field::Gravity::Potential>(3);

    auto fluxCollection         = GravityField::make_fec<field::Gravity::Flux>(3);

    auto displacementCollection = DisplacementField::make_fec<field::Displacement::Vector>(3);

    auto enthalpyCollection     = EnthalpyField::make_fec<field::Enthalpy::Scalar>(3);

    auto vectorL2Collection =
        field::Field<field_mfem_test_utils::VectorL2Field>::make_fec<field_mfem_test_utils::VectorL2Field::Vector>(3);

    auto ndCollection =
        field::Field<field_mfem_test_utils::NdField>::make_fec<field_mfem_test_utils::NdField::Vector>(3);

    REQUIRE(densityCollection != nullptr);

    REQUIRE(potentialCollection != nullptr);

    REQUIRE(fluxCollection != nullptr);

    REQUIRE(displacementCollection != nullptr);

    REQUIRE(enthalpyCollection != nullptr);

    REQUIRE(vectorL2Collection != nullptr);

    REQUIRE(ndCollection != nullptr);

    CHECK(dynamic_cast<mfem::L2_FECollection *>(densityCollection.get()) != nullptr);

    CHECK(dynamic_cast<mfem::L2_FECollection *>(potentialCollection.get()) != nullptr);

    CHECK(dynamic_cast<mfem::RT_FECollection *>(fluxCollection.get()) != nullptr);

    CHECK(dynamic_cast<mfem::H1_FECollection *>(displacementCollection.get()) != nullptr);

    CHECK(dynamic_cast<mfem::H1_FECollection *>(enthalpyCollection.get()) != nullptr);

    CHECK(dynamic_cast<mfem::L2_FECollection *>(vectorL2Collection.get()) != nullptr);

    CHECK(dynamic_cast<mfem::ND_FECollection *>(ndCollection.get()) != nullptr);

    CHECK_THROWS_AS((DensityField::make_fec<field::Density::Scalar>(0)), std::invalid_argument);

    CHECK_THROWS_AS((GravityField::make_fec<field::Gravity::Flux>(-1)), std::invalid_argument);
}

TEST_CASE(
    "Field MFEM Creates Parallel Spaces With Registered Dimensions Orders And Ordering",
    tags::integration &tags::field
) {
    namespace field       = mean_field::field;

    mfem::Mesh serialMesh = mfem::Mesh::MakeCartesian2D(4, 2, mfem::Element::QUADRILATERAL, true, 4.0, 2.0);

    mfem::ParMesh mesh(MPI_COMM_WORLD, serialMesh);

    auto densityFec      = field::Field<field::Density>::make_fec<field::Density::Scalar>(2);

    auto potentialFec    = field::Field<field::Gravity>::make_fec<field::Gravity::Potential>(2);

    auto fluxFec         = field::Field<field::Gravity>::make_fec<field::Gravity::Flux>(2);

    auto displacementFec = field::Field<field::Displacement>::make_fec<field::Displacement::Vector>(2);

    auto enthalpyFec     = field::Field<field::Enthalpy>::make_fec<field::Enthalpy::Scalar>(2);

    auto vectorL2Fec =
        field::Field<field_mfem_test_utils::VectorL2Field>::make_fec<field_mfem_test_utils::VectorL2Field::Vector>(2);

    auto ndFec = field::Field<field_mfem_test_utils::NdField>::make_fec<field_mfem_test_utils::NdField::Vector>(2);

    auto densitySpace   = field::Field<field::Density>::make_fespace<field::Density::Scalar>(mesh, *densityFec);

    auto potentialSpace = field::Field<field::Gravity>::make_fespace<field::Gravity::Potential>(mesh, *potentialFec);

    auto fluxSpace      = field::Field<field::Gravity>::make_fespace<field::Gravity::Flux>(mesh, *fluxFec);

    auto displacementSpace =
        field::Field<field::Displacement>::make_fespace<field::Displacement::Vector>(mesh, *displacementFec);

    auto enthalpySpace = field::Field<field::Enthalpy>::make_fespace<field::Enthalpy::Scalar>(mesh, *enthalpyFec);

    auto vectorL2Space =
        field::Field<field_mfem_test_utils::VectorL2Field>::make_fespace<field_mfem_test_utils::VectorL2Field::Vector>(
            mesh, *vectorL2Fec
        );

    auto ndSpace = field::Field<field_mfem_test_utils::NdField>::make_fespace<field_mfem_test_utils::NdField::Vector>(
        mesh, *ndFec
    );

    REQUIRE(densitySpace != nullptr);
    REQUIRE(potentialSpace != nullptr);
    REQUIRE(fluxSpace != nullptr);
    REQUIRE(displacementSpace != nullptr);
    REQUIRE(enthalpySpace != nullptr);
    REQUIRE(vectorL2Space != nullptr);
    REQUIRE(ndSpace != nullptr);

    CHECK(densitySpace->GetVDim() == 1);

    CHECK(potentialSpace->GetVDim() == 1);

    CHECK(fluxSpace->GetVDim() == 1);

    CHECK(enthalpySpace->GetVDim() == 1);

    CHECK(displacementSpace->GetVDim() == mesh.SpaceDimension());

    CHECK(vectorL2Space->GetVDim() == mesh.SpaceDimension());

    CHECK(ndSpace->GetVDim() == 1);

    CHECK(densitySpace->GetOrdering() == mfem::Ordering::byNODES);

    CHECK(potentialSpace->GetOrdering() == mfem::Ordering::byNODES);

    CHECK(fluxSpace->GetOrdering() == mfem::Ordering::byNODES);

    CHECK(enthalpySpace->GetOrdering() == mfem::Ordering::byNODES);

    /*
     * Displacement deliberately overrides the generic
     * vector-H1 rule and is part of the project's block/indexing
     * contract.
     */
    CHECK(displacementSpace->GetOrdering() == mfem::Ordering::byNODES);

    /*
     * A generic vector L2 quantity retains the ordinary backend
     * realization, demonstrating that the displacement behavior is
     * an intentional specialization rather than a global accident.
     */
    CHECK(vectorL2Space->GetOrdering() == mfem::Ordering::byVDIM);

    CHECK(ndSpace->GetOrdering() == mfem::Ordering::byNODES);

    CHECK(densitySpace->GetMaxElementOrder() == field::Density::Scalar::familyOrder);

    CHECK(potentialSpace->GetMaxElementOrder() == field::Gravity::Potential::familyOrder);

    CHECK(fluxSpace->GetMaxElementOrder() == field::Gravity::Flux::familyOrder + 1);

    CHECK(displacementSpace->GetMaxElementOrder() == field::Displacement::Vector::familyOrder);

    CHECK(enthalpySpace->GetMaxElementOrder() == field::Enthalpy::Scalar::familyOrder);
}

TEST_CASE(
    "Field MFEM Typed Queries Preserve Backend Polynomial Order Semantics",
    tags::unit &tags::field
) {
    namespace field                       = mean_field::field;
    namespace quadrature                  = mean_field::quadrature;
    namespace utils                       = mean_field::utils;

    using DensityField                    = field::Field<field::Density>;

    using GravityField                    = field::Field<field::Gravity>;

    using EnthalpyField                   = field::Field<field::Enthalpy>;

    const quadrature::Query densitySource = DensityField::make_query<field::Density::Form::ProjectionSource>(
        quadrature::QuadratureRole::projection, 3, std::array<int, 1>{4}, utils::DOMAINS::STELLAR,
        quadrature::MappingKind::general
    );

    REQUIRE(densitySource.base_order.has_value());

    /*
     * L2_2 value order 2
     * + geometry order 3
     * + dynamic coefficient order 4.
     */
    CHECK(*densitySource.base_order == 9);

    CHECK(densitySource.term == quadrature::Term::density_projection);

    CHECK(densitySource.role == quadrature::QuadratureRole::projection);

    CHECK(densitySource.domain == utils::DOMAINS::STELLAR);

    CHECK(densitySource.mapping == quadrature::MappingKind::general);

    CHECK(densitySource.geometry_weight_order == 3);

    const quadrature::Query hdivMass =
        GravityField::make_query<field::Gravity::Form::HDivMass>(quadrature::QuadratureRole::discretization, 2);

    REQUIRE(hdivMass.base_order.has_value());

    /*
     * RT_2 value order is 3, hence
     * 3 + 3 + geometry 2 = 8.
     */
    CHECK(*hdivMass.base_order == 8);

    const quadrature::Query divergence = GravityField::make_query<field::Gravity::Form::DivergenceCoupling>(
        quadrature::QuadratureRole::discretization, 2
    );

    REQUIRE(divergence.base_order.has_value());

    /*
     * div(RT_2) order 2
     * + L2_2 order 2
     * + geometry 2.
     */
    CHECK(*divergence.base_order == 6);

    const quadrature::Query pressureForce = EnthalpyField::make_query<field::Enthalpy::Form::PressureForce>(
        quadrature::QuadratureRole::discretization, 2, std::array<int, 1>{9}, utils::DOMAINS::STELLAR,
        quadrature::MappingKind::general
    );

    REQUIRE(pressureForce.base_order.has_value());

    /*
     * h value order 3
     * + grad(d) order 2
     * + geometry 2
     * + n=3 pressure extra order 9
     * = 16.
     */
    CHECK(*pressureForce.base_order == 16);

    const quadrature::Query equilibriumConstant = EnthalpyField::make_query<field::Enthalpy::Form::EquilibriumConstant>(
        quadrature::QuadratureRole::discretization, 2
    );

    REQUIRE(equilibriumConstant.base_order.has_value());

    /*
     * Global scalar C contributes zero polynomial order,
     * h contributes 3, and geometry contributes 2.
     */
    CHECK(*equilibriumConstant.base_order == 5);

    CHECK_THROWS_AS(
        (DensityField::make_query<field::Density::Form::ProjectionMass>(quadrature::QuadratureRole::projection, -1)),
        std::invalid_argument
    );

    const std::array<int, 1> negativeDynamicOrder{-1};

    CHECK_THROWS_AS(
        (EnthalpyField::make_query<field::Enthalpy::Form::PressureForce>(
            quadrature::QuadratureRole::discretization, 2, negativeDynamicOrder
        )),
        std::invalid_argument
    );
}

TEST_CASE(
    "Field MFEM Element Support Resolves Semantic Domains Through The Schema",
    tags::unit &tags::field
) {
    namespace field       = mean_field::field;

    const mfem::Mesh mesh = field_mfem_test_utils::make_two_domain_mesh();

    CHECK((field::element_is_in_field_support<field::Density, field_mfem_test_utils::Schema>(mesh, 0)));

    CHECK_FALSE((field::element_is_in_field_support<field::Density, field_mfem_test_utils::Schema>(mesh, 1)));

    CHECK((field::element_is_in_field_support<field::Enthalpy, field_mfem_test_utils::Schema>(mesh, 0)));

    CHECK_FALSE((field::element_is_in_field_support<field::Enthalpy, field_mfem_test_utils::Schema>(mesh, 1)));

    CHECK((field::element_is_in_field_support<field::Gravity, field_mfem_test_utils::Schema>(mesh, 0)));

    CHECK((field::element_is_in_field_support<field::Gravity, field_mfem_test_utils::Schema>(mesh, 1)));

    CHECK((field::element_is_in_field_support<field::Displacement, field_mfem_test_utils::Schema>(mesh, 0)));

    CHECK((field::element_is_in_field_support<field::Displacement, field_mfem_test_utils::Schema>(mesh, 1)));
}

TEST_CASE(
    "Field MFEM L2 Stellar Support Selects Exactly Stellar Element DOFs",
    tags::unit &tags::field
) {
    namespace field = mean_field::field;

    mfem::Mesh mesh = field_mfem_test_utils::make_two_domain_mesh();

    auto fec        = field::Field<field::Density>::make_fec<field::Density::Scalar>(2);

    mfem::FiniteElementSpace space(&mesh, fec.get());

    const auto support = field::resolve_field_local_dof_support<field::Density, field_mfem_test_utils::Schema>(space);

    const std::vector<int> stellarVDofs = field_mfem_test_utils::decoded_element_vdofs(space, 0);

    const std::vector<int> vacuumVDofs  = field_mfem_test_utils::decoded_element_vdofs(space, 1);

    REQUIRE_FALSE(stellarVDofs.empty());

    REQUIRE_FALSE(vacuumVDofs.empty());

    CHECK(support.activeVDofMarker.Size() == space.GetVSize());

    CHECK(support.activeVDofs.Size() + support.inactiveVDofs.Size() == space.GetVSize());

    for (const int vdof : stellarVDofs) {
        CAPTURE(vdof);

        CHECK(support.activeVDofMarker[vdof] == 1);

        CHECK(field_mfem_test_utils::contains(support.activeVDofs, vdof));

        CHECK_FALSE(field_mfem_test_utils::contains(support.inactiveVDofs, vdof));
    }

    for (const int vdof : vacuumVDofs) {
        CAPTURE(vdof);

        CHECK(support.activeVDofMarker[vdof] == 0);

        CHECK_FALSE(field_mfem_test_utils::contains(support.activeVDofs, vdof));

        CHECK(field_mfem_test_utils::contains(support.inactiveVDofs, vdof));
    }

    CHECK(support.activeVDofs.Size() == static_cast<int>(stellarVDofs.size()));

    CHECK(support.inactiveVDofs.Size() == static_cast<int>(vacuumVDofs.size()));
}

TEST_CASE(
    "Field MFEM H1 Stellar Support Keeps Shared Stellar Vacuum Trace DOFs Active",
    tags::unit &tags::field
) {
    namespace field = mean_field::field;

    mfem::Mesh mesh = field_mfem_test_utils::make_two_domain_mesh();

    auto fec        = field::Field<field::Enthalpy>::make_fec<field::Enthalpy::Scalar>(2);

    mfem::FiniteElementSpace space(&mesh, fec.get());

    const auto support = field::resolve_field_local_dof_support<field::Enthalpy, field_mfem_test_utils::Schema>(space);

    const std::vector<int> stellarVDofs    = field_mfem_test_utils::decoded_element_vdofs(space, 0);

    const std::vector<int> vacuumVDofs     = field_mfem_test_utils::decoded_element_vdofs(space, 1);

    const std::vector<int> interfaceVDofs  = field_mfem_test_utils::intersection(stellarVDofs, vacuumVDofs);

    const std::vector<int> vacuumOnlyVDofs = field_mfem_test_utils::difference(vacuumVDofs, stellarVDofs);

    REQUIRE_FALSE(interfaceVDofs.empty());

    REQUIRE_FALSE(vacuumOnlyVDofs.empty());

    for (const int vdof : stellarVDofs) {
        CAPTURE(vdof);

        CHECK(support.activeVDofMarker[vdof] == 1);
    }

    /*
     * This is the central support invariant:
     *
     * shared interface DOFs are active because they are touched
     * by a supported stellar element, even though they are also
     * touched by a vacuum element.
     */
    for (const int vdof : interfaceVDofs) {
        CAPTURE(vdof);

        CHECK(support.activeVDofMarker[vdof] == 1);

        CHECK(field_mfem_test_utils::contains(support.activeVDofs, vdof));
    }

    for (const int vdof : vacuumOnlyVDofs) {
        CAPTURE(vdof);

        CHECK(support.activeVDofMarker[vdof] == 0);

        CHECK(field_mfem_test_utils::contains(support.inactiveVDofs, vdof));
    }

    CHECK(support.activeVDofs.Size() + support.inactiveVDofs.Size() == space.GetVSize());
}

TEST_CASE(
    "Field MFEM All Domain Support Activates Every L2 And RT DOF",
    tags::unit &tags::field
) {
    namespace field   = mean_field::field;

    mfem::Mesh mesh   = field_mfem_test_utils::make_two_domain_mesh();

    auto potentialFec = field::Field<field::Gravity>::make_fec<field::Gravity::Potential>(2);

    mfem::FiniteElementSpace potentialSpace(&mesh, potentialFec.get());

    const auto potentialSupport =
        field::resolve_field_local_dof_support<field::Gravity, field_mfem_test_utils::Schema>(potentialSpace);

    CHECK(potentialSupport.activeVDofs.Size() == potentialSpace.GetVSize());

    CHECK(potentialSupport.inactiveVDofs.Size() == 0);

    for (int vdof = 0; vdof < potentialSupport.activeVDofMarker.Size(); ++vdof) {
        CHECK(potentialSupport.activeVDofMarker[vdof] == 1);
    }

    /*
     * Exercise signed/oriented MFEM element VDofs through RT as
     * well. The support resolver must DecodeDof() correctly.
     */
    auto fluxFec = field::Field<field::Gravity>::make_fec<field::Gravity::Flux>(2);

    mfem::FiniteElementSpace fluxSpace(&mesh, fluxFec.get());

    const auto fluxSupport =
        field::resolve_field_local_dof_support<field::Gravity, field_mfem_test_utils::Schema>(fluxSpace);

    CHECK(fluxSupport.activeVDofs.Size() == fluxSpace.GetVSize());

    CHECK(fluxSupport.inactiveVDofs.Size() == 0);

    for (int vdof = 0; vdof < fluxSupport.activeVDofMarker.Size(); ++vdof) {
        CHECK(fluxSupport.activeVDofMarker[vdof] == 1);
    }
}

TEST_CASE(
    "Field MFEM Support Resolution Uses Schema Material Bindings Rather Than Hard Coded IDs",
    tags::unit &tags::field
) {
    namespace field = mean_field::field;

    mfem::Mesh mesh = field_mfem_test_utils::make_two_domain_mesh(17, 29);

    auto fec        = field::Field<field::Density>::make_fec<field::Density::Scalar>(2);

    mfem::FiniteElementSpace space(&mesh, fec.get());

    const auto support =
        field::resolve_field_local_dof_support<field::Density, field_mfem_test_utils::AlternateSchema>(space);

    const std::vector<int> stellarVDofs = field_mfem_test_utils::decoded_element_vdofs(space, 0);

    const std::vector<int> vacuumVDofs  = field_mfem_test_utils::decoded_element_vdofs(space, 1);

    for (const int vdof : stellarVDofs) {
        CHECK(support.activeVDofMarker[vdof] == 1);
    }

    for (const int vdof : vacuumVDofs) {
        CHECK(support.activeVDofMarker[vdof] == 0);
    }
}

TEST_CASE(
    "Field MFEM Parallel Stellar Support Produces Consistent Local And True DOF Partitions",
    tags::integration &tags::field
) {
    namespace field       = mean_field::field;

    mfem::Mesh serialMesh = field_mfem_test_utils::make_parallel_split_mesh();

    mfem::ParMesh mesh(MPI_COMM_WORLD, serialMesh);

    auto fec   = field::Field<field::Enthalpy>::make_fec<field::Enthalpy::Scalar>(2);

    auto space = field::Field<field::Enthalpy>::make_fespace<field::Enthalpy::Scalar>(mesh, *fec);

    REQUIRE(space != nullptr);

    const auto support = field::resolve_field_dof_support<field::Enthalpy, field_mfem_test_utils::Schema>(*space);

    CHECK(support.activeVDofMarker.Size() == space->GetVSize());

    CHECK(support.activeVDofs.Size() + support.inactiveVDofs.Size() == space->GetVSize());

    CHECK(support.activeTrueDofMarker.Size() == space->GetTrueVSize());

    CHECK(support.activeTrueDofs.Size() + support.inactiveTrueDofs.Size() == space->GetTrueVSize());

    /*
     * Every local DOF touched by a supported element must be active
     * after shared-DOF synchronization.
     */
    for (int elementId = 0; elementId < mesh.GetNE(); ++elementId) {
        const int materialId = mesh.GetAttribute(elementId);

        const bool stellar =
            field_mfem_test_utils::Schema::template attribute_belongs_to<mean_field::utils::domain::Stellar>(
                materialId
            );

        if (!stellar) {
            continue;
        }

        const std::vector<int> vdofs = field_mfem_test_utils::decoded_element_vdofs(*space, elementId);

        for (const int vdof : vdofs) {
            CAPTURE(elementId, vdof);

            CHECK(support.activeVDofMarker[vdof] == 1);
        }
    }

    const long long globalActiveTrueDofs   = field_mfem_test_utils::global_sum(support.activeTrueDofs.Size());

    const long long globalInactiveTrueDofs = field_mfem_test_utils::global_sum(support.inactiveTrueDofs.Size());

    /*
     * The split mesh contains a finite stellar region and a finite
     * vacuum region with order-three H1 structure, so both categories
     * must genuinely exist globally.
     */
    CHECK(globalActiveTrueDofs > 0);

    CHECK(globalInactiveTrueDofs > 0);
}

TEST_CASE(
    "Field MFEM Parallel All Support Activates Every True Displacement DOF",
    tags::integration &tags::field
) {
    namespace field       = mean_field::field;

    mfem::Mesh serialMesh = field_mfem_test_utils::make_parallel_split_mesh();

    mfem::ParMesh mesh(MPI_COMM_WORLD, serialMesh);

    auto fec   = field::Field<field::Displacement>::make_fec<field::Displacement::Vector>(2);

    auto space = field::Field<field::Displacement>::make_fespace<field::Displacement::Vector>(mesh, *fec);

    REQUIRE(space != nullptr);

    const auto support = field::resolve_field_dof_support<field::Displacement, field_mfem_test_utils::Schema>(*space);

    CHECK(support.inactiveVDofs.Size() == 0);

    CHECK(support.activeVDofs.Size() == space->GetVSize());

    CHECK(support.inactiveTrueDofs.Size() == 0);

    CHECK(support.activeTrueDofs.Size() == space->GetTrueVSize());

    for (int index = 0; index < support.activeTrueDofMarker.Size(); ++index) {
        CHECK(support.activeTrueDofMarker[index] == 1);
    }
}