#include <array>
#include <memory>
#include <span>

#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>

#include "mfem/linalg/densemat.hpp"
#include "serif/backend/mfem/solver/preconditioner/coupling.hpp"
#include "serif/backend/mfem/solver/preconditioner/inverse/diagonal.hpp"

namespace
{
    namespace backend = serif::solver::preconditioner::backend::mfem;
    using CouplingOrder = serif::solver::preconditioner::strategies::coupling::CouplingOrder;

    /** @brief Build a small independent physical block from row-major entries. */
    [[nodiscard]] backend::OperatorHandle material_block(int rows, int columns, std::span<const double> entries)
    {
        auto matrix = std::make_shared<::mfem::DenseMatrix>(rows, columns);
        for (int row = 0; row < rows; ++row)
        {
            for (int column = 0; column < columns; ++column)
            {
                (*matrix)(row, column) = entries[row * columns + column];
            }
        }
        return matrix;
    }

    /** @brief Use the production diagonal inverse for each independently scaled material component. */
    [[nodiscard]] backend::Inverse scalar_inverse(double diagonal)
    {
        const auto inverse = backend::try_diagonal_inverse(std::span<const double>{&diagonal, 1});
        REQUIRE(inverse);
        return *inverse;
    }
} // namespace

TEST_CASE("Material surface policies compose the same two-block inverse with signed triangular algebra", "[unit][preconditioner][material]")
{
    const backend::TwoBlockCouplings material_couplings{.first_from_second = material_block(1, 1, std::array{4.0})};
    const backend::TwoBlockCouplings surface_couplings{.first_from_second = material_block(2, 1, std::array{3.0, 8.0}), .second_from_first = material_block(1, 2, std::array{5.0, 7.0})};
    const auto material = backend::try_coupled_inverse(scalar_inverse(2.0), scalar_inverse(9.0), material_couplings, CouplingOrder::second_then_first);
    REQUIRE(material);
    const auto diagonal_material = backend::try_coupled_inverse(scalar_inverse(2.0), scalar_inverse(9.0), {}, CouplingOrder::independent);
    REQUIRE(diagonal_material);

    // Each subsystem keeps its own native order: material=(density, enthalpy), followed by surface.
    ::mfem::Vector right_hand_side(3), actual(3);
    right_hand_side[0] = 29.0;
    right_hand_side[1] = 43.0;
    right_hand_side[2] = 44.0;
    const double enthalpy = 43.0 / 9.0;
    const double density = (29.0 - 4.0 * enthalpy) / 2.0;
    const double surface = 44.0 / 6.0;
    const double corrected_enthalpy = (43.0 - 8.0 * surface) / 9.0;
    const double ldu_surface = (44.0 - 5.0 * density - 7.0 * enthalpy) / 6.0;
    const double ldu_enthalpy = (43.0 - 8.0 * ldu_surface) / 9.0;
    struct PolicyCase
    {
        backend::Inverse material;
        CouplingOrder order;
        std::array<double, 3> expected;
    };
    const std::array cases{
        PolicyCase{*diagonal_material, CouplingOrder::independent, {14.5, enthalpy, surface}},
        PolicyCase{*material, CouplingOrder::independent, {density, enthalpy, surface}},
        PolicyCase{*material, CouplingOrder::second_then_first, {(29.0 - 3.0 * surface - 4.0 * corrected_enthalpy) / 2.0, corrected_enthalpy, surface}},
        PolicyCase{*material, CouplingOrder::lower_diagonal_upper, {(29.0 - 3.0 * ldu_surface - 4.0 * ldu_enthalpy) / 2.0, ldu_enthalpy, ldu_surface}}};
    for (const PolicyCase& policy : cases)
    {
        const auto inverse = backend::try_coupled_inverse(policy.material, scalar_inverse(6.0), surface_couplings, policy.order);
        REQUIRE(inverse);
        actual = 7.0;
        inverse->action->Mult(right_hand_side, actual);
        for (int index = 0; index < 3; ++index)
        {
            CHECK(actual[index] == Catch::Approx(policy.expected[index]).margin(2e-13));
        }
        ::mfem::Vector zero(3);
        zero = 0.0;
        inverse->action->Mult(zero, actual);
        CHECK(actual.Norml2() == 0.0);
    }

    const auto schur = backend::try_schur_operator(material_block(1, 1, std::array{6.0}), *material, surface_couplings);
    REQUIRE(schur);
}