MeanField_SERiF/src/include/serif/model/equilibrium_plan.hpp

#pragma once

#include "serif/model/equilibrium_equations.hpp"
#include "serif/model/plan.hpp"
#include "serif/model/detail/surface_plan.hpp"
#include "serif/operators/concepts.hpp"
#include "serif/equations/surface_potential.hpp"
#include "serif/equations/primal_gravity.hpp"

namespace serif::model
{
    namespace detail
    {
        template <typename Shape>
        struct ShapePlan;

        template <>
        struct ShapePlan<NoShapeForce>
        {
            using Coefficients = utils::dependencies::KeyList<>;
            using Values = utils::dependencies::KeyList<>;
            using Rows = utils::dependencies::KeyList<>;
        };

        /** @brief The isobar level couples only existing scalar potential, rotation and mass coordinates. */
        template <>
        struct ShapePlan<equations::SurfacePotentialBalance>
        {
            using Coefficients = typename OperandTerms<equations::SurfacePotentialBalance::Operands>::Type;
            using Values = utils::dependencies::KeyList<discretization::blocks::surface_deformation::coefficients>;
            using Rows = utils::dependencies::KeyList<discretization::blocks::surface_deformation::shape_equilibrium>;
        };

    }

    /**
     * @brief Collect the explicitly selected physical coefficients and supported constraint coordinates.
     * @note This plan verifies representation, not mathematical closure. Equal units never merge distinct
     * physical identities, and prescribed coefficients remain available as graph inputs without becoming unknowns.
     */
    template <typename EOS, EquilibriumEquationDescription Description, utils::dependencies::IsKeyList Prescribed = utils::dependencies::KeyList<>,
              utils::dependencies::IsKeyList AdditionalFields = utils::dependencies::KeyList<>>
    struct EquilibriumPlan final
    {
        static constexpr bool has_shape = !std::same_as<typename Description::ShapeForce, NoShapeForce>;
        static constexpr bool has_surface = !std::same_as<typename Description::SurfaceCondition, NoSurfaceCondition>;
        static constexpr bool has_primal_gravity = []<typename... Equations>(utils::types::TypeList<Equations...>) {
            return (std::same_as<Equations, equations::PrimalPoissonBalance> || ...);
        }(typename Description::InteriorEquations::Equations{});
        using Interior = detail::InteriorPlan<typename Description::InteriorEquations>;
        using Shape = detail::ShapePlan<typename Description::ShapeForce>;
        using ResidualEquations = typename Interior::EquationsList;
        using AvailableFields = utils::dependencies::CombinedKeyList<typename Interior::Coefficients, typename Shape::Coefficients, AdditionalFields>;
        using SolvedFields = typename detail::RemoveKeys<AvailableFields, Prescribed>::Type;
        using Residuals = utils::dependencies::CombinedKeyList<typename Interior::Rows, typename Shape::Rows>;
        using Surface = detail::SurfacePlan<typename Description::SurfaceCondition, EOS, AvailableFields, Residuals, SolvedFields, typename Interior::EquationsList>;
        using Coefficients = utils::dependencies::CombinedKeyList<typename Interior::Coefficients, typename Shape::Coefficients, typename Surface::Coefficients>;
        static_assert(!has_shape || !utils::dependencies::ContainsKey<discretization::blocks::surface_deformation::coefficients, Prescribed>,
            "Prescribe geometry by selecting NoShapeForce and supplying a displacement field. An active shape-force equation uses independent solved surface coordinates.");
        using Values = utils::dependencies::CombinedKeyList<typename detail::RemoveKeys<Coefficients, Prescribed>::Type, typename Shape::Values>;
        // Geometry is supplied separately even when no equation solves for the shape.
        using PhysicalInputs = utils::dependencies::CombinedKeyList<Coefficients, Values, utils::dependencies::KeyList<discretization::blocks::surface_deformation::coefficients>>;
        using InactiveRows = std::conditional_t<has_shape, utils::dependencies::KeyList<>, utils::dependencies::KeyList<discretization::blocks::surface_deformation::shape_equilibrium>>;
        using Layout = discretization::blocks::SystemLayout<Values, Residuals>;
    };
}