#include <fesa/fem/gauss_rule.hpp>

#include <array>
#include <limits>
#include <stdexcept>

#include <gtest/gtest.h>

namespace {

TEST(Quadrature, OnePointRuleIntegratesDegreeZeroAndOneExactly) {
    const auto rule = fesa::gauss_rule_1d(1);

    ASSERT_EQ(rule.size(), 1);
    EXPECT_DOUBLE_EQ(rule[0].weight, 2.0);
    EXPECT_DOUBLE_EQ(rule[0].weight * rule[0].xi, 0.0);
}

TEST(Quadrature, TwoPointRuleIntegratesThroughDegreeThreeWithinRoundoff) {
    const auto rule = fesa::gauss_rule_1d(2);

    ASSERT_EQ(rule.size(), 2);

    double degree_zero = 0.0;
    double degree_one = 0.0;
    double degree_two = 0.0;
    double degree_three = 0.0;
    for (const auto& point : rule) {
        const double xi_squared = point.xi * point.xi;
        degree_zero += point.weight;
        degree_one += point.weight * point.xi;
        degree_two += point.weight * xi_squared;
        degree_three += point.weight * xi_squared * point.xi;
    }

    // Sixteen ulps at unit scale cover only the rounding in these short sums.
    constexpr double tolerance =
        16.0 * std::numeric_limits<double>::epsilon();
    EXPECT_NEAR(degree_zero, 2.0, tolerance);
    EXPECT_NEAR(degree_one, 0.0, tolerance);
    EXPECT_NEAR(degree_two, 2.0 / 3.0, tolerance);
    EXPECT_NEAR(degree_three, 0.0, tolerance);
}

TEST(Quadrature, RejectsUnsupportedOrders) {
    for (const int order : std::array{0, 3, -1}) {
        EXPECT_THROW(
            static_cast<void>(fesa::gauss_rule_1d(order)),
            std::invalid_argument);
    }
}

}  // namespace