MeanField/libmeanfield/interface/fem.cppm

module;

#include <memory>
#include <string>

#include <mfem.hpp>
#include <stroid/stroid.h>

export module mean_field:fem;

export import :boundary.contexts;
export import :mapping.domain_mapper;
export import :utils.misc;
export import :utils.user;
export import :quadrature.mfem;
export import :field.mfem;
export import :fem.reference_tables;

export namespace mean_field::fem {
    using GravityField      = field::Field<field::Gravity>;
    using DisplacementField = field::Field<field::Displacement>;
    using DensityField      = field::Field<field::Density>;
    using EnthalpyField     = field::Field<field::Enthalpy>;

    struct FEM {
        // =====================================================================
        // Mesh
        // =====================================================================

        stroid::StroidMesh smesh;
        std::unique_ptr<mfem::ParMesh> mesh;
        std::unique_ptr<mfem::ParMesh> logicalReferenceMesh;

        // =====================================================================
        // Compile-time field descriptors
        // =====================================================================

        GravityField gravityField;
        DisplacementField displacementField;
        DensityField densityField;
        EnthalpyField enthalpyField;

        // =====================================================================
        // Gravity field
        //
        // Collection members are declared before their corresponding spaces so
        // that the spaces are destroyed first.
        // =====================================================================

        std::unique_ptr<mfem::FiniteElementCollection> gravityPotentialFec;

        std::unique_ptr<mfem::ParFiniteElementSpace> gravityPotentialFes;

        std::unique_ptr<mfem::FiniteElementCollection> gravityFluxFec;

        std::unique_ptr<mfem::ParFiniteElementSpace> gravityFluxFes;

        // =====================================================================
        // Displacement field
        // =====================================================================

        std::unique_ptr<mfem::FiniteElementCollection> displacementFec;

        std::unique_ptr<mfem::ParFiniteElementSpace> displacementFes;

        std::unique_ptr<mfem::ParGridFunction> displacement;

        /*
         * Scalar companion of the vector displacement space. Only its
         * StellarSurface true DOFs become surface-deformation coordinates.
         * Sharing displacementFec guarantees identical scalar basis
         * functions without duplicating the finite-element collection.
         */
        std::unique_ptr<mfem::ParFiniteElementSpace> surfaceDeformationFes;

        // =====================================================================
        // Density field
        // =====================================================================

        std::unique_ptr<mfem::FiniteElementCollection> densityFec;

        std::unique_ptr<mfem::ParFiniteElementSpace> densityFes;

        // =====================================================================
        // Specific-enthalpy field
        // =====================================================================

        std::unique_ptr<mfem::FiniteElementCollection> enthalpyFec;

        std::unique_ptr<mfem::ParFiniteElementSpace> enthalpyFes;

        // =====================================================================
        // Compactification coordinate
        // =====================================================================

        std::unique_ptr<mfem::H1_FECollection> compactificationFec;

        std::unique_ptr<mfem::ParFiniteElementSpace> compactificationFes;

        std::unique_ptr<mfem::ParGridFunction> compactificationCoordinate;

        // =====================================================================
        // Domain mapping
        // =====================================================================

        std::unique_ptr<mapping::DomainMapper> domainMapperStateless;

        // =====================================================================
        // Global diagnostics
        // =====================================================================

        mfem::Vector com;
        mfem::DenseMatrix Q;

        // =====================================================================
        // Boundary context
        // =====================================================================

        boundary::BoundaryContext boundaryContext;

        std::unique_ptr<quadrature::RuleFactory> quadratureFactory;

        // =====================================================================
        // Validation
        // =====================================================================

        [[nodiscard]] bool okay() const {
            return mesh != nullptr && logicalReferenceMesh != nullptr &&

                   gravityPotentialFec != nullptr && gravityPotentialFes != nullptr && gravityFluxFec != nullptr &&
                   gravityFluxFes != nullptr &&

                   displacementFec != nullptr && displacementFes != nullptr && displacement != nullptr &&
                   surfaceDeformationFes != nullptr &&

                   densityFec != nullptr && densityFes != nullptr &&

                   enthalpyFec != nullptr && enthalpyFes != nullptr &&

                   compactificationFec != nullptr && compactificationFes != nullptr &&
                   compactificationCoordinate != nullptr &&

                   domainMapperStateless != nullptr && quadratureFactory != nullptr;
        }

        [[nodiscard]] bool has_mapping() const {
            return domainMapperStateless != nullptr && displacement != nullptr && compactificationCoordinate != nullptr;
        }

        [[nodiscard]] const ReferenceTableCache &GetReferenceTables() const {
            return *m_reference_tables;
        }

    private:
        std::unique_ptr<ReferenceTableCache> m_reference_tables{std::make_unique<ReferenceTableCache>()};
    };

    FEM setup_fem(
        const std::string &filename,
        const utils::Args &args,
        int extraRefine = 0
    );
} // namespace mean_field::fem