// Ported from MeanField tests/quadrature/policy.cpp using the consolidated form/role API.

#include <array>
#include <cmath>
#include <format>
#include <stdexcept>
#include <string>

#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <mfem.hpp>

#include "serif/backend/mfem/discretization/quadrature/mfem_resolver.hpp"

#include "serif/tests/test_tags.hpp"

namespace quadrature = serif::discretization::quadrature;
namespace forms = serif::discretization::forms;
namespace blocks = serif::discretization::blocks;
using VectorMass = forms::mass<blocks::displacement::vector, blocks::displacement::vector>;
using DensityMass = forms::mass<blocks::density::mass, blocks::density::mass>;

TEST_CASE(
    "Quadrature Policy Computes Base Orders",
    tags::unit & tags::quadrature
) {
    const quadrature::Policy policy{quadrature::RuleSet{}};
    quadrature::Query query{.form = forms::form_descriptor<VectorMass>};
    query.trial_order = 3;
    query.test_order = 3;
    query.coefficient_order = 2;
    query.geometry_weight_order = 1;
    CHECK(policy.resolve(query).base_order == 9);
    query.trial_order_reduction = 1;
    query.test_order_reduction = 2;
    CHECK(policy.resolve(query).base_order == 6);
    query.trial_order_reduction = 5;
    CHECK(policy.resolve(query).base_order == 4);
    query.base_order = 7;
    CHECK(policy.resolve(query).base_order == 7);
}

TEST_CASE(
    "Quadrature Policy Composes Global Role And Form Boosts",
    tags::unit & tags::quadrature
) {
    quadrature::RuleSet rules;
    rules.set_fallback_control({.boost = 1});
    rules.set_role_control(quadrature::QuadratureRole::preconditioner, {.boost = 2});
    rules.set_form_control<VectorMass>({.boost = 3});
    const quadrature::Policy policy{rules};
    const quadrature::Query query{
        .form = forms::form_descriptor<VectorMass>, .role = quadrature::QuadratureRole::preconditioner, .base_order = 6
    };
    const auto result = policy.resolve(query);
    CHECK(result.base_order == 6);
    CHECK(result.boost == 6);
    CHECK(result.order == 12);
    CHECK_FALSE(result.used_fixed_order);
    CHECK(policy.resolve({.form = forms::form_descriptor<DensityMass>, .base_order = 4}).order == 5);
}

TEST_CASE(
    "Quadrature Fixed Orders Have Defined Precedence",
    tags::unit & tags::quadrature
) {
    quadrature::RuleSet rules;
    const quadrature::Query query{.form = forms::form_descriptor<VectorMass>, .base_order = 6};
    rules.set_fallback_control({.fixed_order = 9, .boost = 100});
    CHECK(quadrature::Policy{rules}.resolve(query).order == 9);
    rules.set_role_control(quadrature::QuadratureRole::discretization, {.fixed_order = 11, .boost = 200});
    CHECK(quadrature::Policy{rules}.resolve(query).order == 11);
    rules.set_form_control<VectorMass>({.fixed_order = 13, .boost = 300});
    const auto result = quadrature::Policy{rules}.resolve(query);
    CHECK(result.order == 13);
    CHECK(result.used_fixed_order);
    CHECK(result.boost == 0);
}

TEST_CASE(
    "Quadrature Policy Rejects Invalid Orders",
    tags::unit & tags::quadrature
) {
    const quadrature::Policy policy{quadrature::RuleSet{}};
    quadrature::Query query{.form = forms::form_descriptor<VectorMass>};
    query.base_order = -1;
    CHECK_THROWS_AS(policy.resolve(query), std::invalid_argument);
    query.base_order.reset();
    query.trial_order = -1;
    CHECK_THROWS_AS(policy.resolve(query), std::invalid_argument);
    query.trial_order = 1;
    query.test_order_reduction = -1;
    CHECK_THROWS_AS(policy.resolve(query), std::invalid_argument);
    quadrature::RuleSet rules;
    rules.set_form_control<VectorMass>({.fixed_order = -1});
    query.test_order_reduction = 0;
    CHECK_THROWS_AS(quadrature::Policy{rules}.resolve(query), std::invalid_argument);
    rules.set_form_control<VectorMass>({.boost = -2});
    CHECK_THROWS_AS(quadrature::Policy{rules}.resolve(query), std::invalid_argument);
}

TEST_CASE(
    "MFEM Rule Factory Includes Geometry And Dynamic Orders",
    tags::unit & tags::quadrature
) {
    using namespace quadrature::backend::mfem;
    const MFEMResolver resolver{quadrature::Policy{quadrature::RuleSet{}}};
    const std::array dynamic_orders{2, 1};
    MFEMQuadratureContext context{
        .geometry = mfem::Geometry::CUBE, .geometry_weight_order = 1, .dynamic_orders = dynamic_orders
    };
    const auto rule = resolver.resolve<VectorMass>(context);
    CHECK(rule.resolution.base_order == 10);
    CHECK(rule.integration_rule == &mfem::IntRules.Get(mfem::Geometry::CUBE, 10));
    context.geometry_weight_order = -1;
    CHECK_THROWS_AS(resolver.resolve<VectorMass>(context), std::invalid_argument);
    context.geometry_weight_order = 0;
    const std::array invalid_orders{-1};
    context.dynamic_orders = invalid_orders;
    CHECK_THROWS_AS(resolver.resolve<VectorMass>(context), std::invalid_argument);
}

TEST_CASE(
    "MFEM Quadrature Rule Integrates Tensor Polynomial Exactly",
    tags::unit & tags::quadrature
) {
    const quadrature::backend::mfem::MFEMResolver resolver{quadrature::Policy{quadrature::RuleSet{}}};
    const auto rule = resolver.resolve<VectorMass>({.geometry = mfem::Geometry::CUBE});
    double integral = 0.0;
    for (int i = 0; i < rule.integration_rule->GetNPoints(); ++i) {
        const auto& point = rule.integration_rule->IntPoint(i);
        integral += point.weight * std::pow(point.x, 4) * std::pow(point.y, 2) * point.z;
    }
    CHECK_THAT(integral, Catch::Matchers::WithinAbs(1.0 / 30.0, 1.0e-14));
}

TEST_CASE(
    "Policy Controlled Vector Mass Assembly Matches Overintegrated Reference",
    tags::unit & tags::quadrature
) {
    auto mesh = mfem::Mesh::MakeCartesian3D(1, 1, 1, mfem::Element::HEXAHEDRON, 1.2, 0.9, 1.4);
    mfem::H1_FECollection collection(3, 3);
    mfem::FiniteElementSpace space(&mesh, &collection, 3);
    auto& transformation = *mesh.GetElementTransformation(0);
    const auto& element = *space.GetFE(0);
    const quadrature::backend::mfem::MFEMResolver resolver{quadrature::Policy{quadrature::RuleSet{}}};
    mfem::VectorMassIntegrator controlled;
    const auto resolution = resolver.configure<VectorMass>(
        controlled, {.geometry = element.GetGeomType(), .geometry_weight_order = transformation.OrderW()}
    );
    mfem::VectorMassIntegrator reference;
    reference.SetIntRule(&mfem::IntRules.Get(element.GetGeomType(), resolution.order + 8));
    mfem::DenseMatrix actual, expected;
    controlled.AssembleElementMatrix(element, transformation, actual);
    reference.AssembleElementMatrix(element, transformation, expected);
    actual -= expected;
    CHECK(actual.FNorm() / expected.FNorm() < 1.0e-12);
}

TEST_CASE(
    "Integration Rule Adapter Supports Each Advertised Signature",
    tags::unit & tags::quadrature
) {
    struct ReferenceIntegrator {
        const mfem::IntegrationRule* selected = nullptr;
        void SetIntegrationRule(const mfem::IntegrationRule& rule) {
            selected = &rule;
        }
    } reference;
    struct PointerIntegrator {
        const mfem::IntegrationRule* selected = nullptr;
        void SetIntRule(const mfem::IntegrationRule* rule) {
            selected = rule;
        }
    } pointer;
    struct IntRuleReferenceIntegrator {
        const mfem::IntegrationRule* selected = nullptr;
        void SetIntRule(const mfem::IntegrationRule& rule) {
            selected = &rule;
        }
    } int_reference;

    const quadrature::backend::mfem::MFEMResolver resolver{quadrature::Policy{quadrature::RuleSet{}}};
    const quadrature::backend::mfem::MFEMQuadratureContext context{.geometry = mfem::Geometry::CUBE};
    resolver.configure<VectorMass>(reference, context);
    resolver.configure<VectorMass>(pointer, context);
    resolver.configure<VectorMass>(int_reference, context);
    CHECK(reference.selected == &mfem::IntRules.Get(mfem::Geometry::CUBE, 6));
    CHECK(pointer.selected == reference.selected);
    CHECK(int_reference.selected == reference.selected);
}

TEST_CASE(
    "Discrete Functional Space Can Be Formatted",
    tags::unit & tags::quadrature
) {
    const blocks::displacement::vector::FunctionalSpace space;
    CHECK(std::format("{}", space).find("rank=Vector") != std::string::npos);
}