#include <array>
#include <cmath>

#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
#include "mfem/fem/fe/fe_h1.hpp"
#include "serif/backend/mfem/discretization/basis_evaluation.hpp"

TEST_CASE("Surface power retains nodal density and its analytic surface limit", "[unit][reconstruction]") {
    namespace reconstruction = serif::discretization::reconstruction;
    const std::array nodes{0.0, 0.2, 0.7, 1.0};
    std::array<double, 4> values, slopes;
    for (const double exponent : {1.0, 4.0 / 3.0, 5.0 / 3.0, 2.0, 2.8}) {
        const reconstruction::SurfacePowerProfile profile(exponent);
        for (std::size_t node = 0; node < nodes.size(); ++node) {
            profile.evaluate(nodes, nodes[node], values, slopes);
            for (std::size_t basis = 0; basis < nodes.size(); ++basis) {
                CHECK(values[basis] == Catch::Approx(node == basis ? 1.0 : 0.0).margin(2e-14));
                CHECK(std::isfinite(slopes[basis]));
            }
        }
        for (const double point : {0.0, 1e-12, 0.01, 0.35, 0.9, 1.0}) {
            profile.evaluate(nodes, point, values, slopes);
            double density = 0.0, derivative = 0.0, constant = 0.0, constant_derivative = 0.0;
            for (std::size_t basis = 0; basis < nodes.size(); ++basis) {
                density += values[basis] * std::pow(nodes[basis], exponent);
                derivative += slopes[basis] * std::pow(nodes[basis], exponent);
                constant += values[basis];
                constant_derivative += slopes[basis];
            }
            CHECK(density == Catch::Approx(std::pow(point, exponent)).margin(3e-14));
            CHECK(derivative == Catch::Approx(exponent * std::pow(point, exponent - 1.0)).margin(3e-13));
            CHECK(constant == Catch::Approx(1.0).margin(3e-14));
            CHECK(constant_derivative == Catch::Approx(0.0).margin(3e-13));
        }
    }
    CHECK_THROWS(reconstruction::SurfacePowerProfile(0.8));
    CHECK_THROWS(reconstruction::SurfacePowerProfile(NAN));
}

TEST_CASE("Reconstructed tensor values gradients and transposes agree in every surface orientation", "[unit][reconstruction]") {
    namespace backend = serif::discretization::backend::mfem;
    namespace reconstruction = serif::discretization::reconstruction;
    ::mfem::H1_HexahedronElement element(3);
    const auto& rule = ::mfem::IntRules.Get(::mfem::Geometry::CUBE, 12);
    backend::TensorBasisCache cache;
    backend::BasisWorkspace workspace;
    for (const double exponent : {1.0, 4.0 / 3.0, 2.8}) {
        const reconstruction::SurfacePowerProfile profile(exponent);
        for (int axis = 0; axis < 3; ++axis) {
            for (const bool high : {false, true}) {
                const backend::ElementReconstruction reconstruction{&profile, axis, high};
                for (const auto operation : {backend::ReferenceOperation::value, backend::ReferenceOperation::gradient, backend::ReferenceOperation::value_and_gradient}) {
                    backend::BasisEvaluation dense, tensor;
                    dense.sample(element, rule, operation, 1, cache, backend::ContractionMode::dense, &reconstruction);
                    tensor.sample(element, rule, operation, 1, cache, backend::ContractionMode::automatic, &reconstruction);
                    REQUIRE(tensor.uses_tensor());
                    ::mfem::Vector coefficients(element.GetDof()), dense_values, tensor_values;
                    for (int i = 0; i < coefficients.Size(); ++i) { coefficients[i] = std::sin(0.3 * i); }
                    dense.evaluate(coefficients, dense_values, workspace);
                    tensor.evaluate(coefficients, tensor_values, workspace);
                    tensor_values -= dense_values;
                    CHECK(tensor_values.Normlinf() < 2e-12);
                    ::mfem::Vector dual(dense.height()), dense_transpose, tensor_transpose;
                    for (int i = 0; i < dual.Size(); ++i) { dual[i] = std::cos(0.13 * i); }
                    dense.integrate(dual, dense_transpose, workspace);
                    tensor.integrate(dual, tensor_transpose, workspace);
                    CHECK((coefficients * tensor_transpose) == Catch::Approx(dense_values * dual).margin(1e-9));
                    tensor_transpose -= dense_transpose;
                    CHECK(tensor_transpose.Normlinf() < 2e-10);
                }
            }
        }
        // The cache borrows the immutable profile, just as it borrows finite elements.
        cache = backend::TensorBasisCache{};
    }
}