#include <catch2/catch_test_macros.hpp>

#include "serif/tests/equations/fixture.hpp"
#include "serif/tests/test_tags.hpp"

namespace serif::tests::equations {
    TEST_CASE("Equation Closure Geometry Has The Analytic Dilation Factor", tags::unit & tags::equations) {
        Fixture f;
        const auto rho = constant(f.density_binding.size(), 0.55), h = constant(f.enthalpy_binding.size(), 0.9);
        const auto zero_rho = constant(rho.Size(), 0.0), zero_h = constant(h.Size(), 0.0), zero_u = f.displacement(0.0);
        constexpr double dilation = 0.2, direction_scale = 0.07;
        ::mfem::VectorFunctionCoefficient coefficient(3, [](const ::mfem::Vector& x, ::mfem::Vector& u) {
            u = x;
            u *= dilation;
        });
        ::mfem::GridFunction displacement(&f.displacement_space);
        displacement.ProjectCoefficient(coefficient);
        const auto u = f.displacement_binding.gather(displacement);
        auto du = u;
        du *= direction_scale / dilation;
        PreparedClosure prepared(
            Closure{eos::models::Polytrope{3.0, 1.5}}, f.density_binding, f.displacement_binding, f.compactification,
            f.mapper
        );
        ::mfem::Vector identity, deformed, action;
        prepared.prepare(f.state(rho, h), zero_u, f.workspace);
        prepared.residual(f.workspace, identity);
        prepared.prepare(f.state(rho, h), u, f.workspace);
        prepared.residual(f.workspace, deformed);
        prepared.apply_jacobian(f.state(zero_rho, zero_h), du, f.workspace, action);
        identity *= std::pow(1.0 + dilation, 3.0);
        CHECK(relative_error(deformed, identity) < 5.0e-12);
        deformed *= 3.0 * direction_scale / (1.0 + dilation);
        CHECK(relative_error(action, deformed) < 5.0e-12);
    }
} // namespace serif::tests::equations