module;
#include <cstdint>
#include <expected>
#include <memory>
#include <mfem.hpp>
#include <stdexcept>
#include <vector>

export module mean_field:operators.prepared_hdiv_mass;
export import :fem;
export import :field.mfem;
export import :mapping.domain_mapper;
import :fem.reference_tables;

export namespace mean_field::operators {
    enum class HDivMassPreparationRejectionReason : std::uint8_t { invalid_mapping, non_finite_arithmetic };

    struct HDivMassPreparationRejection final {
        HDivMassPreparationRejectionReason reason{HDivMassPreparationRejectionReason::invalid_mapping};
        mapping::MappingStatus mappingStatus{mapping::MappingStatus::valid};
    };

    using HDivMassPreparationResult = std::expected<void, HDivMassPreparationRejection>;

    [[noreturn]] inline void throwHDivMassPreparationRejection(const HDivMassPreparationRejection &rejection) {
        if (rejection.reason == HDivMassPreparationRejectionReason::non_finite_arithmetic) {
            throw std::domain_error("Prepared H(div) mass data contained non-finite arithmetic.");
        }
        throw std::domain_error("Prepared H(div) mass data could not map the candidate geometry.");
    }

    class PreparedMappedHDivMassOperator final : public mfem::Operator {
    public:
        PreparedMappedHDivMassOperator(
            const fem::FEM &f,
            const mapping::DomainMapper &domain_mapper
        );

        void Prepare(const mfem::Vector &displacement);
        void PreparePrimal(const mfem::Vector &displacement);
        [[nodiscard]] HDivMassPreparationResult TryPrepare(const mfem::Vector &displacement);
        [[nodiscard]] HDivMassPreparationResult TryPreparePrimal(const mfem::Vector &displacement);
        void Mult(
            const mfem::Vector &gravity_gradient,
            mfem::Vector &action
        ) const override;
        void MultDisplacementVariationTrue(
            const mfem::Vector &gravityGradientTrue,
            const mfem::Vector &displacementVariationTrue,
            mfem::Vector &actionVariationTrue
        ) const;
        void AssembleDiagonal(mfem::Vector &diagonal) const override;
        void AssembleTrueDiagonal(mfem::Vector &diagonal) const;

        [[nodiscard]] bool IsPrepared() const noexcept;
        [[nodiscard]] bool HasVariationData() const noexcept;
        [[nodiscard]] std::uint64_t GetPreparationCount() const noexcept;

        [[nodiscard]] const field::FieldDofMap &GetFluxMap() const noexcept;
        [[nodiscard]] const field::FieldDofMap &GetDisplacementMap() const noexcept;

        template <typename Visitor> void VisitMappedGeometryRules(Visitor &&visitor) const {
            for (const ElementVariationData &data : m_variationElements) {
                visitor(data.elementId, *data.integrationRule);
            }
        }

    private:
        enum class PreparationMode : std::uint8_t { primal, linearization };

        struct ElementVariationData {
            int elementId{-1};
            mfem::Array<int> gravityGradientDofs;
            mfem::Array<int> displacementDofs;
            mfem::Array<int> compactificationDofs;
            mfem::DofTransformation *gravityGradientDofTransformation{nullptr};
            mfem::DofTransformation *displacementDofTransformation{nullptr};
            mfem::Vector baseDisplacement;
            mfem::Vector compactification;
            const mfem::IntegrationRule *integrationRule{nullptr};
            std::shared_ptr<const fem::VectorReferenceTable> gravityReferenceTable;
            // Fixed computational-mesh Piola factor, separate from J_map.
            mfem::DenseMatrix meshPiolaJacobians;
            mfem::Vector referenceWeights;
            mfem::DenseMatrix frozenMappingData;
        };

        [[nodiscard]] mapping::MappingStatus PrepareVariationData();
        [[nodiscard]] HDivMassPreparationResult TryPrepareImpl(
            const mfem::Vector &displacement,
            PreparationMode mode
        );

        const fem::FEM &m_fem;
        const mapping::DomainMapper &m_domain_mapper;

        field::FieldDofMap m_flux_map;
        field::FieldDofMap m_displacement_map;

        mfem::Array<int> m_stellar_marker;
        mfem::Array<int> m_vacuum_marker;

        std::unique_ptr<mfem::MatrixCoefficient> m_stellar_mass_coefficient;
        std::unique_ptr<mfem::MatrixCoefficient> m_vacuum_mass_coefficient;
        std::unique_ptr<mfem::ParBilinearForm> m_stellar_mass_form;
        std::unique_ptr<mfem::ParBilinearForm> m_vacuum_mass_form;
        mutable mfem::Vector m_flux_true;
        mutable mfem::Vector m_action_true;
        mutable mfem::Vector m_domain_action_true;
        mutable mfem::Vector m_flux_local;
        mutable mfem::Vector m_action_local;
        mutable mfem::Vector m_domain_action_local;
        mfem::Vector m_displacement_true;
        std::vector<ElementVariationData> m_variationElements;

        mutable mapping::DomainMapper::Workspace m_variationWorkspace;
        mutable mapping::VolumeMappingContext m_baseMappingContext;
        mutable mapping::VolumeMappingVariation m_mappingVariation;
        mutable mfem::Vector m_gravityGradientLocal;
        mutable mfem::Vector m_displacementVariationLocal;
        mutable mfem::Vector m_localVariationAction;
        mutable mfem::Vector m_elementGravityGradient;
        mutable mfem::Vector m_elementDisplacementVariation;
        mutable mfem::Vector m_elementVariationAction;
        mutable mfem::Vector m_gravityGradientValue;
        mutable mfem::Vector m_gravityReferenceCellValue;
        mutable mfem::Vector m_referenceCellDual;
        mutable mfem::Vector m_massTensorVariationAction;
        mutable mfem::DenseMatrix m_meshPiolaJacobian;
        mutable mfem::DenseMatrix m_gravityGradientShape;
        mutable mfem::DenseMatrix m_massTensorVariation;
        std::uint64_t m_preparation_count{0};
        bool m_is_prepared{false};
        bool m_has_variation_data{false};
        bool m_single_rank{true};
    };
} // namespace mean_field::operators