#include <concepts>
#include <string_view>
#include <type_traits>

#include <catch2/catch_test_macros.hpp>

import mean_field;

namespace {
    namespace blocks   = mean_field::utils::blocks;
    namespace eos      = mean_field::eos;
    namespace field    = mean_field::field;
    namespace material = mean_field::material;
    namespace surface  = mean_field::surface;

    struct Entropy final : eos::ThermodynamicQuantity {
        static constexpr std::string_view identifier = "mock_entropy";
    };

    struct Composition final : eos::ThermodynamicQuantity {
        static constexpr std::string_view identifier = "mock_composition";
    };

    struct Temperature final : eos::ThermodynamicQuantity {
        static constexpr std::string_view identifier = "mock_temperature";
    };

    struct EntropyField final {
        static constexpr std::string_view name = "mock_entropy";
        using PhysicalQuantity                 = Entropy;
    };

    struct CompositionField final {
        static constexpr std::string_view name = "mock_composition";
        using PhysicalQuantity                 = Composition;
    };

    struct TemperatureField final {
        static constexpr std::string_view name = "mock_temperature";
        using PhysicalQuantity                 = Temperature;
    };

    struct AlternateEntropyField final {
        static constexpr std::string_view name = "alternate_mock_entropy";
        using PhysicalQuantity                 = Entropy;
    };

    struct EntropyValue final : blocks::value_block_base { };
    struct EntropyResidual final : blocks::residual_block_base { };
    struct CompositionValue final : blocks::value_block_base { };
    struct CompositionResidual final : blocks::residual_block_base { };
    struct TemperatureValue final : blocks::value_block_base { };
    struct TemperatureResidual final : blocks::residual_block_base { };
    struct DensityValue final : blocks::value_block_base { };
    struct DensityResidual final : blocks::residual_block_base { };
    struct EnthalpyValue final : blocks::value_block_base { };
    struct EnthalpyResidual final : blocks::residual_block_base { };
    struct UnrelatedValue final : blocks::value_block_base { };
    struct UnrelatedResidual final : blocks::residual_block_base { };

    using EntropyEquation = material::ThermodynamicEquation<EntropyField, EntropyValue, EntropyResidual>;
    using CompositionEquation =
        material::ThermodynamicEquation<CompositionField, CompositionValue, CompositionResidual>;
    using TemperatureEquation =
        material::ThermodynamicEquation<TemperatureField, TemperatureValue, TemperatureResidual>;
    using DensityEquation  = material::ThermodynamicEquation<field::Density, DensityValue, DensityResidual>;
    using EnthalpyEquation = material::ThermodynamicEquation<field::Enthalpy, EnthalpyValue, EnthalpyResidual>;

    using EnthalpyFromPressureEntropyComposition =
        eos::Relation<eos::quantity::SpecificEnthalpy, eos::quantity::Pressure, Entropy, Composition>;

    class GeneralEquationOfState final {
    public:
        using Relations = eos::RelationCatalog<EnthalpyFromPressureEntropyComposition>;

        [[nodiscard]] constexpr eos::SpecificEnthalpyValue evaluate(
            EnthalpyFromPressureEntropyComposition,
            const eos::PressureValue pressure,
            const eos::QuantityValue<Entropy> entropy,
            const eos::QuantityValue<Composition> composition
        ) const noexcept {
            return eos::SpecificEnthalpyValue{
                2.0 * pressure.value() + 3.0 * entropy.value() + 5.0 * composition.value()
            };
        }

        [[nodiscard]] constexpr eos::PartialDerivative<
            eos::quantity::SpecificEnthalpy,
            Entropy>
        partialDerivative(
            EnthalpyFromPressureEntropyComposition,
            eos::WithRespectTo<Entropy>,
            eos::PressureValue,
            eos::QuantityValue<Entropy>,
            eos::QuantityValue<Composition>
        ) const noexcept {
            return eos::PartialDerivative<eos::quantity::SpecificEnthalpy, Entropy>{3.0};
        }

        [[nodiscard]] constexpr eos::PartialDerivative<
            eos::quantity::SpecificEnthalpy,
            Composition>
        partialDerivative(
            EnthalpyFromPressureEntropyComposition,
            eos::WithRespectTo<Composition>,
            eos::PressureValue,
            eos::QuantityValue<Entropy>,
            eos::QuantityValue<Composition>
        ) const noexcept {
            return eos::PartialDerivative<eos::quantity::SpecificEnthalpy, Composition>{5.0};
        }
    };

    class MissingPressureSurfaceRelationEquationOfState final {
    public:
        using Relations = eos::RelationCatalog<eos::PressureFromDensity>;

        [[nodiscard]] constexpr eos::PressureValue evaluate(
            eos::PressureFromDensity,
            const eos::DensityValue density
        ) const noexcept {
            return eos::PressureValue{density.value()};
        }
    };

    using EnthalpyFromPressureEntropy =
        eos::Relation<eos::quantity::SpecificEnthalpy, eos::quantity::Pressure, Entropy>;

    class AmbiguousEquationOfState final {
    public:
        using Relations = eos::RelationCatalog<EnthalpyFromPressureEntropyComposition, EnthalpyFromPressureEntropy>;

        [[nodiscard]] constexpr eos::SpecificEnthalpyValue evaluate(
            EnthalpyFromPressureEntropyComposition,
            const eos::PressureValue pressure,
            const eos::QuantityValue<Entropy> entropy,
            const eos::QuantityValue<Composition> composition
        ) const noexcept {
            return eos::SpecificEnthalpyValue{pressure.value() + entropy.value() + composition.value()};
        }

        [[nodiscard]] constexpr eos::SpecificEnthalpyValue evaluate(
            EnthalpyFromPressureEntropy,
            const eos::PressureValue pressure,
            const eos::QuantityValue<Entropy> entropy
        ) const noexcept {
            return eos::SpecificEnthalpyValue{pressure.value() + entropy.value()};
        }

        [[nodiscard]] constexpr eos::PartialDerivative<
            eos::quantity::SpecificEnthalpy,
            Entropy>
        partialDerivative(
            EnthalpyFromPressureEntropyComposition,
            eos::WithRespectTo<Entropy>,
            eos::PressureValue,
            eos::QuantityValue<Entropy>,
            eos::QuantityValue<Composition>
        ) const noexcept {
            return eos::PartialDerivative<eos::quantity::SpecificEnthalpy, Entropy>{1.0};
        }

        [[nodiscard]] constexpr eos::PartialDerivative<
            eos::quantity::SpecificEnthalpy,
            Composition>
        partialDerivative(
            EnthalpyFromPressureEntropyComposition,
            eos::WithRespectTo<Composition>,
            eos::PressureValue,
            eos::QuantityValue<Entropy>,
            eos::QuantityValue<Composition>
        ) const noexcept {
            return eos::PartialDerivative<eos::quantity::SpecificEnthalpy, Composition>{1.0};
        }

        [[nodiscard]] constexpr eos::PartialDerivative<
            eos::quantity::SpecificEnthalpy,
            Entropy>
        partialDerivative(
            EnthalpyFromPressureEntropy,
            eos::WithRespectTo<Entropy>,
            eos::PressureValue,
            eos::QuantityValue<Entropy>
        ) const noexcept {
            return eos::PartialDerivative<eos::quantity::SpecificEnthalpy, Entropy>{1.0};
        }
    };

    using GeneralForm = blocks::block_form<
        blocks::type_list<CompositionValue, UnrelatedValue, EntropyValue, EnthalpyValue>,
        blocks::type_list<CompositionResidual, UnrelatedResidual, EntropyResidual, EnthalpyResidual>>;

    using ReorderedAvailableEquations = material::
        ThermodynamicEquationCatalog<TemperatureEquation, EnthalpyEquation, EntropyEquation, CompositionEquation>;

    using DifferentlyReorderedAvailableEquations = material::
        ThermodynamicEquationCatalog<CompositionEquation, TemperatureEquation, EntropyEquation, EnthalpyEquation>;

    using ExpectedGeneralEquations =
        material::ThermodynamicEquationCatalog<CompositionEquation, EntropyEquation, EnthalpyEquation>;

    using GeneralCompilation =
        material::CompiledThermodynamicEquationsT<GeneralEquationOfState, GeneralForm, ReorderedAvailableEquations>;

    using GeneralCompilationFromOtherOrder = material::
        CompiledThermodynamicEquationsT<GeneralEquationOfState, GeneralForm, DifferentlyReorderedAvailableEquations>;

    using HalfPresentEquationForm = blocks::block_form<
        blocks::type_list<CompositionValue, EntropyValue, EnthalpyValue>,
        blocks::type_list<CompositionResidual, EnthalpyResidual>>;

    using NoThermodynamicEquationForm =
        blocks::block_form<blocks::type_list<UnrelatedValue>, blocks::type_list<UnrelatedResidual>>;

    using DuplicateFieldCatalog = material::ThermodynamicEquationCatalog<
        EntropyEquation,
        material::ThermodynamicEquation<EntropyField, TemperatureValue, TemperatureResidual>>;

    using DuplicateQuantityCatalog = material::ThermodynamicEquationCatalog<
        EntropyEquation,
        material::ThermodynamicEquation<AlternateEntropyField, TemperatureValue, TemperatureResidual>>;

    using DensityFromPressureEntropy = eos::Relation<eos::quantity::Density, eos::quantity::Pressure, Entropy>;

    class DensityCarrierEquationOfState final {
    public:
        using Relations = eos::RelationCatalog<DensityFromPressureEntropy>;

        [[nodiscard]] constexpr eos::DensityValue evaluate(
            DensityFromPressureEntropy,
            const eos::PressureValue pressure,
            const eos::QuantityValue<Entropy> entropy
        ) const noexcept {
            return eos::DensityValue{4.0 * pressure.value() + 2.0 * entropy.value()};
        }

        [[nodiscard]] constexpr eos::PartialDerivative<
            eos::quantity::Density,
            Entropy>
        partialDerivative(
            DensityFromPressureEntropy,
            eos::WithRespectTo<Entropy>,
            eos::PressureValue,
            eos::QuantityValue<Entropy>
        ) const noexcept {
            return eos::PartialDerivative<eos::quantity::Density, Entropy>{2.0};
        }
    };

    using DensityCarrierForm = blocks::
        block_form<blocks::type_list<EntropyValue, DensityValue>, blocks::type_list<EntropyResidual, DensityResidual>>;
    using DensityCarrierAvailableEquations =
        material::ThermodynamicEquationCatalog<DensityEquation, TemperatureEquation, EntropyEquation>;
    using DensityCarrierCompilation = material::CompiledThermodynamicEquationsT<
        DensityCarrierEquationOfState,
        DensityCarrierForm,
        DensityCarrierAvailableEquations>;

    struct DensityCarrierState final {
        double density;
        double entropy;

        [[nodiscard]] constexpr eos::DensityValue value(eos::quantity::Density) const noexcept {
            return eos::DensityValue{density};
        }

        [[nodiscard]] constexpr eos::QuantityValue<Entropy> value(Entropy) const noexcept {
            return eos::QuantityValue<Entropy>{entropy};
        }
    };
} // namespace

TEST_CASE(
    "Thermodynamic Equation Compilation Is Inferred From Fields And Canonical Problem Blocks",
    "[material][thermodynamic_equations][compilation][unit][type_contract]"
) {
    STATIC_CHECK(material::ThermodynamicField<EntropyField>);
    STATIC_CHECK(material::ThermodynamicField<field::Density>);
    STATIC_CHECK(material::ThermodynamicField<field::Enthalpy>);
    STATIC_CHECK(material::ThermodynamicEquationType<EntropyEquation>);
    STATIC_CHECK(material::ValidThermodynamicEquationCatalog<ReorderedAvailableEquations>);
    STATIC_CHECK(material::CompiledThermodynamicEquations<GeneralCompilation>);
    STATIC_CHECK(std::same_as<typename GeneralCompilation::Equations, ExpectedGeneralEquations>);
    STATIC_CHECK(std::same_as<GeneralCompilation, GeneralCompilationFromOtherOrder>);
    STATIC_CHECK(GeneralCompilation::Equations::size == 3);
    STATIC_CHECK(
        std::same_as<
            typename GeneralCompilation::StateBindings,
            surface::SurfaceStateBindings<
                surface::SurfaceStateBinding<Composition, CompositionField>,
                surface::SurfaceStateBinding<Entropy, EntropyField>,
                surface::SurfaceStateBinding<eos::quantity::SpecificEnthalpy, field::Enthalpy>>>
    );
    STATIC_CHECK(std::same_as<typename GeneralCompilation::PressureSurfaceFormulation::CarrierField, field::Enthalpy>);
}

TEST_CASE(
    "Thermodynamic Equation Compilation Rejects Incomplete Ambiguous And Missing Physics",
    "[material][thermodynamic_equations][compilation][unit][negative]"
) {
    STATIC_CHECK_FALSE(
        material::ThermodynamicEquationsCompilable<
            GeneralEquationOfState, HalfPresentEquationForm, ReorderedAvailableEquations>
    );
    STATIC_CHECK_FALSE(
        material::ThermodynamicEquationsCompilable<
            GeneralEquationOfState, NoThermodynamicEquationForm, ReorderedAvailableEquations>
    );
    STATIC_CHECK_FALSE(
        material::ThermodynamicEquationsCompilable<
            MissingPressureSurfaceRelationEquationOfState, GeneralForm, ReorderedAvailableEquations>
    );
    STATIC_CHECK_FALSE(
        material::ThermodynamicEquationsCompilable<AmbiguousEquationOfState, GeneralForm, ReorderedAvailableEquations>
    );
    STATIC_CHECK_FALSE(material::ValidThermodynamicEquationCatalog<DuplicateFieldCatalog>);
    STATIC_CHECK_FALSE(material::ValidThermodynamicEquationCatalog<DuplicateQuantityCatalog>);
}

TEST_CASE(
    "Pressure Surface Carrier Is Selected By The EOS Relation Rather Than A Preferred Field",
    "[material][thermodynamic_equations][surface][unit]"
) {
    using Formulation = DensityCarrierCompilation::PressureSurfaceFormulation;
    using Constraint  = surface::CompiledPressureSurfaceConstraintT<Formulation, DensityCarrierEquationOfState>;

    STATIC_CHECK(material::CompiledThermodynamicEquations<DensityCarrierCompilation>);
    STATIC_CHECK(std::same_as<typename Formulation::CarrierQuantity, eos::quantity::Density>);
    STATIC_CHECK(std::same_as<typename Formulation::CarrierField, field::Density>);
    STATIC_CHECK(std::same_as<typename Constraint::Relation, DensityFromPressureEntropy>);
    STATIC_CHECK(
        std::same_as<
            typename Constraint::SurfaceDependencies::StateFieldTypes, field::TypeList<field::Density, EntropyField>>
    );

    constexpr DensityCarrierEquationOfState equationOfState;
    const surface::ConstantPressureSurface condition{eos::PressureValue{0.25}};
    const auto constraint = surface::compilePressureSurfaceConstraint<Formulation>(condition, equationOfState);

    constexpr DensityCarrierState state{.density = 1.6, .entropy = 0.3};
    constexpr DensityCarrierState direction{.density = -0.2, .entropy = 0.4};

    CHECK(constraint.residual(state) == 0.0);
    CHECK(constraint.jacobianAction(state, direction) == -1.0);
}