diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index 3f57b4daf..060fbbe09 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -174,8 +174,8 @@ endif() add_subdirectory(Core) add_subdirectory(FE/Common) add_subdirectory(FE/Math) -add_subdirectory(FE/Basis) add_subdirectory(FE/Quadrature) +add_subdirectory(FE/Basis) set(SVMP_ARCHITECTURE_TARGETS svmp_core diff --git a/Code/Source/solver/FE/Basis/CMakeLists.txt b/Code/Source/solver/FE/Basis/CMakeLists.txt index 166985ab8..3d0ee337f 100644 --- a/Code/Source/solver/FE/Basis/CMakeLists.txt +++ b/Code/Source/solver/FE/Basis/CMakeLists.txt @@ -20,4 +20,5 @@ target_include_directories(svmp_fe_basis PUBLIC target_link_libraries(svmp_fe_basis PUBLIC svmp_fe_common svmp_fe_math + svmp_fe_quadrature ) diff --git a/Code/Source/solver/FE/Basis/NodeOrderingConventions.cpp b/Code/Source/solver/FE/Basis/NodeOrderingConventions.cpp index b0262551a..09cf5c47c 100644 --- a/Code/Source/solver/FE/Basis/NodeOrderingConventions.cpp +++ b/Code/Source/solver/FE/Basis/NodeOrderingConventions.cpp @@ -4,9 +4,10 @@ #include "NodeOrderingConventions.h" #include "BasisExceptions.h" #include "BasisTraits.h" +#include "FE/Quadrature/GaussLobattoNodes.h" #include -#include +#include #include #include #include @@ -22,11 +23,8 @@ namespace { using Point = math::Vector; using Lattice = std::array; -// Gauss-Lobatto-Legendre nodes on [-1, 1] for a degree-`order` distribution -// (order + 1 nodes). The endpoints are -1 and +1; the interior nodes are the -// roots of P'_order, found by Newton iteration on f(x) = x P_order(x) - -// P_{order-1}(x) -- whose roots are exactly the GLL nodes -- from the -// Chebyshev-Gauss-Lobatto seed. +// Cache GLL tensor-axis nodes by interpolation order; Quadrature owns the +// shared refinement, while Basis retains its correction tolerance and errors. const std::vector& gll_points(int order) { thread_local std::map> cache; const auto found = cache.find(order); @@ -34,62 +32,22 @@ const std::vector& gll_points(int order) { return found->second; } - // Newton converges quadratically from the Chebyshev-Gauss-Lobatto seed, so a - // few iterations suffice for any practical order; the cap is only a safety - // bound. Reaching it without meeting the tolerance signals a real failure - constexpr int kMaxNewtonIterations = 100; - constexpr double kNewtonTolerance = double(1e-15); - - std::vector pts(static_cast(order + 1), double(0)); - if (order >= 1) { - pts.front() = double(-1); - pts.back() = double(1); - } - const double pi = std::acos(double(-1)); - const int half = order / 2; - for (int j = 1; j <= half; ++j) { - if (2 * j == order) { - pts[static_cast(j)] = double(0); // exact center, even order - continue; + constexpr double kNewtonTolerance = 1.0e-15; + const std::size_t num_points = static_cast(order) + 1; + std::vector pts(num_points, 0.0); + pts.front() = -1.0; + pts.back() = 1.0; + try { + for (std::size_t half_root_index = 0; + half_root_index < (num_points - 1) / 2; ++half_root_index) { + const auto node = quadrature::detail::refine_gauss_lobatto_node( + num_points, half_root_index, kNewtonTolerance); + pts[1u + half_root_index] = -node.coordinate; + pts[num_points - 2u - half_root_index] = node.coordinate; } - double x = -std::cos(pi * static_cast(j) / static_cast(order)); - bool converged = false; - for (int iter = 0; iter < kMaxNewtonIterations; ++iter) { - // Legendre P_k and P'_k up to k = order at x, by the three-term - // recurrences (regular at x = +/-1). - double p_km1 = double(1); // P_0 - double p_k = x; // P_1 - double d_km1 = double(0); // P'_0 - double d_k = double(1); // P'_1 - for (int k = 1; k < order; ++k) { - const double kk = static_cast(k); - const double inv = double(1) / (kk + double(1)); - const double p_kp1 = - ((double(2) * kk + double(1)) * x * p_k - kk * p_km1) * inv; - const double d_kp1 = - ((double(2) * kk + double(1)) * (p_k + x * d_k) - kk * d_km1) * inv; - p_km1 = p_k; - p_k = p_kp1; - d_km1 = d_k; - d_k = d_kp1; - } - // p_k = P_order, p_km1 = P_{order-1}, d_k = P'_order, d_km1 = P'_{order-1}. - const double f = x * p_k - p_km1; - const double f_prime = p_k + x * d_k - d_km1; - const double dx = f / f_prime; - x -= dx; - if (std::abs(dx) <= kNewtonTolerance) { - converged = true; - break; - } - } - svmp::throw_if( - !converged, "ReferenceNodeLayout: Gauss-Lobatto-Legendre Newton iteration did not converge " - "(order outside the trustworthy range)"); - pts[static_cast(j)] = x; - } - for (int j = half + 1; j < order; ++j) { - pts[static_cast(j)] = -pts[static_cast(order - j)]; + } catch (const ConvergenceException& error) { + svmp::raise( + "ReferenceNodeLayout: " + error.message()); } auto inserted = cache.emplace(order, std::move(pts)); diff --git a/Code/Source/solver/FE/Quadrature/GaussLobattoNodes.cpp b/Code/Source/solver/FE/Quadrature/GaussLobattoNodes.cpp new file mode 100644 index 000000000..ef3b9b115 --- /dev/null +++ b/Code/Source/solver/FE/Quadrature/GaussLobattoNodes.cpp @@ -0,0 +1,165 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +#include "FE/Quadrature/GaussLobattoNodes.h" + +#include "FE/Common/FEException.h" + +#include +#include +#include +#include +#include + +namespace svmp::FE::quadrature::detail { +namespace { + +// Bound cosine-seeded Newton refinement for deterministic termination. +constexpr int kMaximumNewtonIterations = 100; + +// Return (P_degree(x), P_(degree-1)(x)) for degree >= 1 using the Legendre +// three-term recurrence, starting from P_0 = 1 and P_1 = x. +std::pair evaluate_adjacent_legendre_values( + std::size_t degree, + double coordinate) noexcept +{ + double previous_value = 1.0; + double value = coordinate; + + for (std::size_t recurrence_degree = 2; recurrence_degree <= degree; ++recurrence_degree) { + const double next_value = + ((2 * recurrence_degree - 1) * coordinate * value - + (recurrence_degree - 1) * previous_value) / recurrence_degree; + previous_value = value; + value = next_value; + } + + return {value, previous_value}; +} + +// Preserve the failed quantity and generator context in a convergence error. +[[noreturn]] void raise_generation_failure(std::size_t num_points, std::size_t half_root_index, + int iteration, double diagnostic_value, std::string_view detail) +{ + std::ostringstream message; + message << "Gauss-Lobatto-Legendre generator: " << detail + << ", num_points=" << num_points + << ", half_root_index=" << half_root_index + << ", diagnostic_value=" << diagnostic_value; + + svmp::raise( + message.str(), iteration, std::abs(diagnostic_value)); +} + +} // namespace + +// Refine f(x) = x*P_(n-1)(x) - P_(n-2)(x), then recheck its correction at +// the returned node. Callers retain their existing correction tolerances. +GaussLobattoNode refine_gauss_lobatto_node(std::size_t num_points, + std::size_t half_root_index, double correction_tolerance) +{ + if (num_points < 3 || half_root_index >= (num_points - 1) / 2) { + svmp::raise( + "Gauss-Lobatto-Legendre nodes: invalid interior node request"); + } + if (!(std::isfinite(correction_tolerance) && correction_tolerance > 0.0)) { + svmp::raise( + "Gauss-Lobatto-Legendre nodes: correction tolerance must be finite and positive"); + } + + const bool is_center = + 1u + half_root_index == num_points - 2u - half_root_index; + const std::size_t polynomial_degree = num_points - 1; + const double num_points_value = num_points; + const double degree_value = polynomial_degree; + const double pi = std::numbers::pi_v; + double root = std::cos(pi * (half_root_index + 1) / degree_value); + double correction = 0.0; + + // For f = x*P_m - P_(m-1), Legendre identities give + // f' = (m+1)*P_m = n*P_m exactly. + for (int iteration = 1; iteration <= kMaximumNewtonIterations; ++iteration) { + const auto [polynomial_value, previous_polynomial_value] = + evaluate_adjacent_legendre_values(polynomial_degree, root); + if (!(std::isfinite(polynomial_value) && + std::isfinite(previous_polynomial_value))) { + raise_generation_failure( + num_points, half_root_index, iteration, polynomial_value, + "encountered invalid adjacent Legendre values"); + } + + const double residual = + root * polynomial_value - previous_polynomial_value; + const double derivative = num_points_value * polynomial_value; + if (!(std::isfinite(residual) && + std::isfinite(derivative) && + derivative != 0.0)) { + raise_generation_failure( + num_points, half_root_index, iteration, derivative, + "computed an invalid root-function residual or derivative"); + } + + correction = residual / derivative; + const double updated_root = root - correction; + if (!(std::isfinite(correction) && std::isfinite(updated_root))) { + raise_generation_failure( + num_points, half_root_index, iteration, correction, + "computed an invalid Newton update"); + } + root = updated_root; + + if (std::abs(correction) > correction_tolerance) { + continue; + } + + if (is_center) { + root = 0.0; + } + if (!(root >= 0.0 && root < 1.0 && (is_center || root > 0.0))) { + raise_generation_failure( + num_points, half_root_index, iteration, root, + "refined root is outside the expected half interval"); + } + + const auto [final_polynomial_value, + final_previous_polynomial_value] = + evaluate_adjacent_legendre_values(polynomial_degree, root); + if (!(std::isfinite(final_polynomial_value) && + std::isfinite(final_previous_polynomial_value))) { + raise_generation_failure( + num_points, half_root_index, iteration, final_polynomial_value, + "refined root produced invalid adjacent Legendre values"); + } + + const double final_residual = + root * final_polynomial_value - + final_previous_polynomial_value; + const double final_derivative = + num_points_value * final_polynomial_value; + if (!(std::isfinite(final_residual) && + std::isfinite(final_derivative) && + final_derivative != 0.0)) { + raise_generation_failure( + num_points, half_root_index, iteration, final_derivative, + "refined root produced an invalid residual or derivative"); + } + + const double final_correction = + final_residual / final_derivative; + if (!(std::isfinite(final_correction) && + std::abs(final_correction) <= + correction_tolerance)) { + raise_generation_failure( + num_points, half_root_index, iteration, final_correction, + "refined root failed final correction validation"); + } + + return {root, final_polynomial_value, iteration}; + } + + raise_generation_failure( + num_points, half_root_index, kMaximumNewtonIterations, correction, + "Newton refinement did not converge"); +} + +} // namespace svmp::FE::quadrature::detail diff --git a/Code/Source/solver/FE/Quadrature/GaussLobattoNodes.h b/Code/Source/solver/FE/Quadrature/GaussLobattoNodes.h new file mode 100644 index 000000000..b117d0f73 --- /dev/null +++ b/Code/Source/solver/FE/Quadrature/GaussLobattoNodes.h @@ -0,0 +1,35 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +#ifndef SVMP_FE_GAUSS_LOBATTO_NODES_H +#define SVMP_FE_GAUSS_LOBATTO_NODES_H + +#include + +namespace svmp::FE::quadrature::detail { + +/** @cond INTERNAL */ + +// Shared node refinement for interpolation and integration. The Legendre +// value is retained so Quadrature can form a weight without another recurrence. +struct GaussLobattoNode { + double coordinate; + double legendre_value; // P_(num_points-1) at coordinate. + int iterations; +}; + +// Refine a nonnegative interior node of an n-point GLL distribution on [-1, 1]. +// half_root_index runs from the positive endpoint toward the center, starting +// at zero; an odd distribution's center is returned as exactly zero. +// Requires n >= 3, index < (n-1)/2, and a finite positive correction tolerance. +// Invalid requests throw InvalidArgumentException; unsuccessful refinement +// throws ConvergenceException. Callers own endpoints, mirroring, and size caps. +[[nodiscard]] GaussLobattoNode refine_gauss_lobatto_node( + std::size_t num_points, std::size_t half_root_index, + double correction_tolerance); + +/** @endcond */ + +} // namespace svmp::FE::quadrature::detail + +#endif // SVMP_FE_GAUSS_LOBATTO_NODES_H diff --git a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp new file mode 100644 index 000000000..2046e9f3c --- /dev/null +++ b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp @@ -0,0 +1,143 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +/** + * @file GaussLobattoQuadrature.cpp + * @brief Exactness-requested generation of bounded Gauss-Lobatto-Legendre line rules. + * @ingroup FE_Quadrature + */ + +#include "FE/Quadrature/GaussLobattoQuadrature.h" + +#include "FE/Common/FEException.h" +#include "FE/Quadrature/GaussLobattoNodes.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace svmp::FE::quadrature { +namespace { + +constexpr std::size_t kMaximumPoints = 128; +static_assert(max_gauss_lobatto_exactness() == 2 * kMaximumPoints - 3); + +// Guard recurrence and Newton-update rounding; exhaustive supported-size +// sweeps qualify this scale. +constexpr double kNewtonCorrectionTolerance = + 64.0 * std::numeric_limits::epsilon(); +// Provide conservative O(n epsilon) accumulation headroom, qualified by those +// sweeps. +constexpr double kRuleValidationTolerance = + 32.0 * kMaximumPoints * + std::numeric_limits::epsilon(); + +// Preserve the failed quantity and generator context in a convergence error. +// A half-root index or iteration of -1 identifies a rule-wide validation check. +[[noreturn]] void raise_generation_failure(std::size_t num_points, int half_root_index, + int iteration, double diagnostic_value, std::string_view detail) +{ + std::ostringstream message; + message << "Gauss-Lobatto-Legendre generator: " << detail + << ", num_points=" << num_points + << ", half_root_index=" << half_root_index + << ", diagnostic_value=" << diagnostic_value; + + svmp::raise( + message.str(), iteration, std::abs(diagnostic_value)); +} + +// Form a weight from the polynomial value already computed during refinement. +std::pair generate_interior_root_and_weight(std::size_t num_points, + std::size_t half_root_index, double weight_denominator_scale) +{ + const auto node = detail::refine_gauss_lobatto_node( + num_points, half_root_index, kNewtonCorrectionTolerance); + const double denominator = + weight_denominator_scale * node.legendre_value * node.legendre_value; + if (!(std::isfinite(denominator) && denominator > 0.0)) { + raise_generation_failure( + num_points, half_root_index, node.iterations, denominator, + "refined root produced an invalid weight denominator"); + } + + const double weight = 2.0 / denominator; + if (!(std::isfinite(weight) && weight > 0.0)) { + raise_generation_failure( + num_points, half_root_index, node.iterations, weight, + "refined root produced an invalid quadrature weight"); + } + return {node.coordinate, weight}; +} + +} // namespace + +QuadratureRule make_gauss_lobatto_rule(int requested_exactness) +{ + svmp::check( + requested_exactness >= 0 && + requested_exactness <= max_gauss_lobatto_exactness(), + "Gauss-Lobatto-Legendre generator: requested_exactness must be in [0, " + + std::to_string(max_gauss_lobatto_exactness()) + ']'); + + const std::size_t num_points = requested_exactness / 2 + 2; + std::vector points(num_points, QuadPoint::Zero()); + std::vector weights(points.size()); + + points.front()[0] = -1.0; + points.back()[0] = 1.0; + + const double weight_denominator_scale = + num_points * (num_points - 1); + const double endpoint_weight = 2.0 / weight_denominator_scale; + weights.front() = endpoint_weight; + weights.back() = endpoint_weight; + + const std::size_t interior_roots_to_refine = (num_points - 1) / 2; + for (std::size_t half_root_index = 0; + half_root_index < interior_roots_to_refine; ++half_root_index) { + const std::size_t left_index = 1u + half_root_index; + const std::size_t right_index = + points.size() - 2u - half_root_index; + const auto [root, weight] = + generate_interior_root_and_weight( + num_points, half_root_index, weight_denominator_scale); + + points[left_index][0] = -root; + points[right_index][0] = root; + weights[left_index] = weight; + weights[right_index] = weight; + } + + for (std::size_t point_index = 1; point_index < points.size(); ++point_index) { + const double spacing = points[point_index][0] - points[point_index - 1u][0]; + if (!(spacing > 0.0)) { + raise_generation_failure( + num_points, point_index, -1, spacing, + "generated points are not strictly increasing"); + } + } + + // Report a failed measure instead of repairing or rescaling the weights. + const long double weight_sum = std::accumulate(weights.begin(), weights.end(), 0.0L); + const long double measure_error = std::abs(weight_sum - 2.0L); + if (!(std::isfinite(weight_sum) && + measure_error <= kRuleValidationTolerance)) { + raise_generation_failure( + num_points, -1, -1, measure_error, + "generated weights do not reproduce the reference measure"); + } + + const int polynomial_exactness = 2 * num_points - 3; + return QuadratureRule( + svmp::CellFamily::Line, polynomial_exactness, + std::move(points), std::move(weights)); +} + +} // namespace svmp::FE::quadrature diff --git a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.h b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.h new file mode 100644 index 000000000..d00b453be --- /dev/null +++ b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.h @@ -0,0 +1,63 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +#ifndef SVMP_FE_GAUSS_LOBATTO_QUADRATURE_H +#define SVMP_FE_GAUSS_LOBATTO_QUADRATURE_H + +/** + * @file GaussLobattoQuadrature.h + * @brief Bounded Gauss-Lobatto-Legendre quadrature generation on the reference line. + * @ingroup FE_Quadrature + */ + +#include "FE/Quadrature/QuadratureRule.h" + +namespace svmp::FE::quadrature { + +/** @addtogroup FE_Quadrature + * @{ + */ + +/** + * @brief Return the largest supported Gauss-Lobatto-Legendre exactness request. + * @details Degree 253 corresponds to the internal 128-point project support + * bound, including both endpoints, which limits generator work and downstream + * product-rule growth; it is not a mathematical or convergence limit. + * @return The inclusive requested-exactness limit, 253. + */ +[[nodiscard]] constexpr int max_gauss_lobatto_exactness() noexcept +{ + return 253; +} + +/** + * @brief Generate a Gauss-Lobatto-Legendre rule on @f$[-1,1]@f$ with at least + * the requested exactness. + * + * @details For @f$d@f$ = @p requested_exactness, integer division gives the + * minimum point count @f$n=\lfloor d/2\rfloor+2@f$. The returned metadata reports + * the actual polynomial exactness @f$2n-3@f$, which exceeds even requests by one. + * Degree zero produces two points with exactness one. The first and last points + * are exactly @f$-1@f$ and @f$+1@f$. When present, the @f$n-2@f$ interior points + * are the roots of @f$P'_{n-1}@f$, the derivative of the Legendre polynomial + * of degree @f$n-1@f$. Points are strictly increasing, and weights are positive + * and aligned with their points. + * + * @see [NIST DLMF: Legendre polynomials](https://dlmf.nist.gov/18.3) + * + * @param requested_exactness Minimum polynomial degree to integrate exactly; + * must be in [0, 253], inclusive (see max_gauss_lobatto_exactness()). + * @return A complete QuadratureRule value for CellFamily::Line. + * @throws InvalidArgumentException If @p requested_exactness is outside the + * supported range; checked before point-count conversion or allocation. + * @throws ConvergenceException If root refinement or final numerical + * validation fails. + */ +[[nodiscard]] QuadratureRule +make_gauss_lobatto_rule(int requested_exactness); + +/** @} */ + +} // namespace svmp::FE::quadrature + +#endif // SVMP_FE_GAUSS_LOBATTO_QUADRATURE_H diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp new file mode 100644 index 000000000..b57ffe8dc --- /dev/null +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp @@ -0,0 +1,234 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +/** + * @file GaussQuadrature.cpp + * @brief Exactness-requested generation of bounded Gauss-Legendre line rules. + * @ingroup FE_Quadrature + */ + +#include "FE/Quadrature/GaussQuadrature.h" + +#include "FE/Common/FEException.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace svmp::FE::quadrature { +namespace { + +constexpr std::size_t kMaximumPoints = 128; +static_assert(max_gauss_legendre_exactness() == 2 * kMaximumPoints - 1); + +// Defensively bound supported cosine-seeded Newton refinements for +// deterministic termination. +constexpr int kMaximumNewtonIterations = 100; +// Guard recurrence and Newton-update rounding; exhaustive supported-size +// sweeps qualify this scale. +constexpr double kNewtonCorrectionTolerance = + 64.0 * std::numeric_limits::epsilon(); +// Provide conservative O(n epsilon) accumulation headroom, qualified by those +// sweeps. +constexpr double kRuleValidationTolerance = + 32.0 * kMaximumPoints * + std::numeric_limits::epsilon(); + +// Return (P_degree(x), P'_degree(x)) for degree >= 1. Advance the Legendre +// three-term recurrence and its derivative together, starting from P_0 = 1 +// and P_1 = x, without numerical differentiation. +std::pair evaluate_legendre_with_derivative( + std::size_t degree, + double coordinate) noexcept +{ + double previous_value = 1.0; + double previous_derivative = 0.0; + double value = coordinate; + double derivative = 1.0; + for (std::size_t recurrence_degree = 2; recurrence_degree <= degree; ++recurrence_degree) { + const double degree_value = recurrence_degree; + const double recurrence_factor = 2 * recurrence_degree - 1; + const double next_value = + (recurrence_factor * coordinate * value - + (recurrence_degree - 1) * previous_value) / degree_value; + const double next_derivative = + (recurrence_factor * (value + coordinate * derivative) - + (recurrence_degree - 1) * previous_derivative) / degree_value; + + previous_value = value; + previous_derivative = derivative; + value = next_value; + derivative = next_derivative; + } + + return {value, derivative}; +} + +// Preserve the failed quantity and generator context in a convergence error. +// A root index or iteration of -1 identifies a rule-wide validation check. +[[noreturn]] void raise_generation_failure( + std::size_t num_points, + int root_index, + int iteration, + double diagnostic_value, + std::string_view detail) +{ + std::ostringstream message; + message << "Gauss-Legendre generator: " << detail + << ", num_points=" << num_points + << ", root_index=" << root_index + << ", diagnostic_value=" << diagnostic_value; + + svmp::raise( + message.str(), iteration, std::abs(diagnostic_value)); +} + +// Refine a nonnegative root of P_n with cosine-seeded Newton iteration, bounded +// by the iteration and correction limits above. Recheck P_n/P'_n before forming +// w = 2 / ((1 - x*x) * P'_n(x)^2). The caller mirrors (x, w); is_center assigns +// the odd rule's center exactly to zero before the final validation. +std::pair generate_root_and_weight( + std::size_t num_points, std::size_t root_index, bool is_center) +{ + const double pi = std::numbers::pi_v; + double root = std::cos( + pi * (root_index + 0.75) / (num_points + 0.5)); + double correction = 0.0; + + for (int iteration = 1; iteration <= kMaximumNewtonIterations; ++iteration) { + const auto [polynomial_value, polynomial_derivative] = + evaluate_legendre_with_derivative(num_points, root); + if (!(std::isfinite(polynomial_value) && + std::isfinite(polynomial_derivative) && + polynomial_derivative != 0.0)) { + raise_generation_failure( + num_points, root_index, iteration, polynomial_derivative, + "encountered an invalid Legendre value or derivative"); + } + + correction = polynomial_value / polynomial_derivative; + const double updated_root = root - correction; + if (!(std::isfinite(correction) && std::isfinite(updated_root))) { + raise_generation_failure( + num_points, root_index, iteration, correction, + "computed an invalid Newton update"); + } + root = updated_root; + + if (std::abs(correction) > kNewtonCorrectionTolerance) { + continue; + } + + if (is_center) { + root = 0.0; + } + if (!(root >= 0.0 && root < 1.0 && + (is_center || root > 0.0))) { + raise_generation_failure( + num_points, root_index, iteration, root, + "refined root is outside the expected half interval"); + } + + const auto [final_polynomial_value, + final_polynomial_derivative] = + evaluate_legendre_with_derivative(num_points, root); + if (!(std::isfinite(final_polynomial_value) && + std::isfinite(final_polynomial_derivative) && + final_polynomial_derivative != 0.0)) { + raise_generation_failure( + num_points, root_index, iteration, final_polynomial_derivative, + "refined root produced an invalid Legendre value or derivative"); + } + + const double final_correction = + final_polynomial_value / final_polynomial_derivative; + if (!(std::isfinite(final_correction) && + std::abs(final_correction) <= + kNewtonCorrectionTolerance)) { + raise_generation_failure( + num_points, root_index, iteration, final_correction, + "refined root failed final correction validation"); + } + + const double denominator = + (1.0 - root) * (1.0 + root) * + final_polynomial_derivative * final_polynomial_derivative; + if (!(std::isfinite(denominator) && denominator > 0.0)) { + raise_generation_failure( + num_points, root_index, iteration, denominator, + "refined root produced an invalid weight denominator"); + } + + const double weight = 2.0 / denominator; + if (!(std::isfinite(weight) && weight > 0.0)) { + raise_generation_failure( + num_points, root_index, iteration, weight, + "refined root produced an invalid quadrature weight"); + } + + return {root, weight}; + } + + raise_generation_failure( + num_points, root_index, kMaximumNewtonIterations, correction, + "Newton refinement did not converge"); +} + +} // namespace + +QuadratureRule make_gauss_legendre_rule(int requested_exactness) +{ + svmp::check( + requested_exactness >= 0 && + requested_exactness <= max_gauss_legendre_exactness(), + "Gauss-Legendre generator: requested_exactness must be in [0, " + + std::to_string(max_gauss_legendre_exactness()) + ']'); + + const std::size_t num_points = requested_exactness / 2 + 1; + std::vector points(num_points, QuadPoint::Zero()); + std::vector weights(points.size()); + + const std::size_t roots_to_refine = (num_points + 1) / 2; + for (std::size_t left_index = 0; left_index < roots_to_refine; ++left_index) { + const std::size_t right_index = points.size() - 1u - left_index; + const auto [root, weight] = generate_root_and_weight( + num_points, left_index, left_index == right_index); + + points[left_index][0] = -root; + points[right_index][0] = root; + weights[left_index] = weight; + weights[right_index] = weight; + } + + for (std::size_t point_index = 1; point_index < points.size(); ++point_index) { + const double spacing = points[point_index][0] - points[point_index - 1u][0]; + if (!(spacing > 0.0)) { + raise_generation_failure( + num_points, point_index, -1, spacing, + "generated points are not strictly increasing"); + } + } + + // Report a failed measure instead of repairing or rescaling the weights. + const long double weight_sum = std::accumulate(weights.begin(), weights.end(), 0.0L); + const long double measure_error = std::abs(weight_sum - 2.0L); + if (!(std::isfinite(weight_sum) && + measure_error <= kRuleValidationTolerance)) { + raise_generation_failure( + num_points, -1, -1, measure_error, + "generated weights do not reproduce the reference measure"); + } + + const int polynomial_exactness = 2 * num_points - 1; + return QuadratureRule( + svmp::CellFamily::Line, polynomial_exactness, + std::move(points), std::move(weights)); +} + +} // namespace svmp::FE::quadrature diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.h b/Code/Source/solver/FE/Quadrature/GaussQuadrature.h new file mode 100644 index 000000000..b9d594a0e --- /dev/null +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.h @@ -0,0 +1,62 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +#ifndef SVMP_FE_GAUSS_QUADRATURE_H +#define SVMP_FE_GAUSS_QUADRATURE_H + +/** + * @file GaussQuadrature.h + * @brief Bounded Gauss-Legendre quadrature generation on the reference line. + * @ingroup FE_Quadrature + */ + +#include "FE/Quadrature/QuadratureRule.h" + +namespace svmp::FE::quadrature { + +/** @addtogroup FE_Quadrature + * @{ + */ + +/** + * @brief Return the largest supported Gauss-Legendre exactness request. + * @details Degree 255 corresponds to the internal 128-point project support + * bound, which limits generator work and downstream product-rule growth; it + * is not a mathematical or convergence limit. + * @return The inclusive requested-exactness limit, 255. + */ +[[nodiscard]] constexpr int max_gauss_legendre_exactness() noexcept +{ + return 255; +} + +/** + * @brief Generate a Gauss-Legendre rule on @f$[-1,1]@f$ with at least the + * requested exactness. + * + * @details For @f$d@f$ = @p requested_exactness, integer division gives the + * minimum point count @f$n=\lfloor d/2\rfloor+1@f$. The returned metadata reports + * the actual polynomial exactness @f$2n-1@f$, which exceeds even requests by one. + * Degree zero produces one point with exactness one. Points are the roots of + * the Legendre polynomial @f$P_n@f$ of degree @f$n@f$, in strictly increasing + * order inside @f$(-1,1)@f$; neither endpoint is included. Weights are positive + * and aligned with their points. + * + * @see [NIST DLMF: Legendre polynomials](https://dlmf.nist.gov/18.3) + * + * @param requested_exactness Minimum polynomial degree to integrate exactly; + * must be in [0, 255], inclusive (see max_gauss_legendre_exactness()). + * @return A complete QuadratureRule value for CellFamily::Line. + * @throws InvalidArgumentException If @p requested_exactness is outside the + * supported range; checked before point-count conversion or allocation. + * @throws ConvergenceException If root refinement or final numerical + * validation fails. + */ +[[nodiscard]] QuadratureRule +make_gauss_legendre_rule(int requested_exactness); + +/** @} */ + +} // namespace svmp::FE::quadrature + +#endif // SVMP_FE_GAUSS_QUADRATURE_H diff --git a/Documentation/Doxyfile b/Documentation/Doxyfile index a4502c58c..48803a491 100644 --- a/Documentation/Doxyfile +++ b/Documentation/Doxyfile @@ -108,6 +108,7 @@ WARN_LOGFILE = INPUT = \ Documentation/internal.md \ Code/Source \ + tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp \ tests/unitTests/linear_solver_tests \ tests/unitTests/ionic_model_tests \ tests/unitTests/active_stress_tests \ diff --git a/tests/unitTests/FE/Basis/test_LagrangeBasis.cpp b/tests/unitTests/FE/Basis/test_LagrangeBasis.cpp index 48c5897a1..a24bcb96a 100644 --- a/tests/unitTests/FE/Basis/test_LagrangeBasis.cpp +++ b/tests/unitTests/FE/Basis/test_LagrangeBasis.cpp @@ -869,7 +869,8 @@ TEST(LagrangeBasis, TensorAxisNodesAreGaussLobattoLegendre) { EXPECT_EQ(line_coord_pm_one(1, 2), double(0)); EXPECT_EQ(line_coord_pm_one(2, 2), double(1)); - for (int order = 1; order <= 12; ++order) { + // Higher orders also guard against inheriting Quadrature's 128-point cap. + for (const int order : {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 128, 129, 256}) { EXPECT_EQ(line_coord_pm_one(0, order), double(-1)) << "order=" << order; EXPECT_EQ(line_coord_pm_one(order, order), double(1)) << "order=" << order; @@ -889,4 +890,5 @@ TEST(LagrangeBasis, TensorAxisNodesAreGaussLobattoLegendre) { // the endpoints. At order 4 the first interior node is -sqrt(3/7) ~ -0.6547, // past the equispaced position -0.5. EXPECT_LT(line_coord_pm_one(1, 4), double(-0.5) - double(1e-6)); + EXPECT_EQ(line_coord_pm_one(0, 0), double(0)); } diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp new file mode 100644 index 000000000..28bbd3859 --- /dev/null +++ b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp @@ -0,0 +1,549 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +/** + * @file test_QuadratureGenerators.cpp + * @brief Tests for exactness-requested one-dimensional quadrature generators. + */ + +#include + +#include "FE/Basis/NodeOrderingConventions.h" +#include "FE/Common/FEException.h" +#include "FE/Quadrature/GaussLobattoNodes.h" +#include "FE/Quadrature/GaussLobattoQuadrature.h" +#include "FE/Quadrature/GaussQuadrature.h" +#include "FE/Quadrature/QuadratureRule.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace svmp::FE; +using namespace svmp::FE::quadrature; + +namespace { + +// Exhaustive supported-domain sweeps observed structure, measure, and moment +// errors below 1.2e-15; these envelopes retain cross-toolchain headroom +// through degree 255. +constexpr double kStructureTolerance = 1.0e-12; +constexpr double kMomentTolerance = 2.0e-12; +// Computed fixtures follow the production 64*epsilon refinement scale. +constexpr double kFixtureTolerance = + 64.0 * std::numeric_limits::epsilon(); +// Basis and Quadrature use different correction tolerances, so terminal +// binary64 rounding can differ; this remains far below node spacing. +constexpr double kBasisConsistencyTolerance = + 128.0 * std::numeric_limits::epsilon(); + +static_assert(max_gauss_legendre_exactness() == 255); +static_assert(noexcept(max_gauss_legendre_exactness())); +static_assert( + std::is_same_v); +static_assert(max_gauss_lobatto_exactness() == 253); +static_assert(noexcept(max_gauss_lobatto_exactness())); +static_assert( + std::is_same_v); + +enum class LineEndpointPolicy { + Excluded, + Included, +}; + +long double analytic_line_monomial_integral(std::size_t power) +{ + if (power % 2u != 0u) { + return 0.0L; + } + return 2.0L / (static_cast(power) + 1.0L); +} + +long double accumulate_line_moment( + const QuadratureRule& rule, + std::size_t power) +{ + long double sum = 0.0L; + long double correction = 0.0L; + + for (std::size_t point_index = 0; + point_index < rule.num_points(); + ++point_index) { + const long double coordinate = + static_cast(rule.point(point_index)[0]); + const long double term = + static_cast(rule.weight(point_index)) * + std::pow(coordinate, static_cast(power)); + const long double next_sum = sum + term; + + if (std::abs(sum) >= std::abs(term)) { + correction += (sum - next_sum) + term; + } else { + correction += (term - next_sum) + sum; + } + sum = next_sum; + } + + return sum + correction; +} + +void expect_common_line_metadata( + const QuadratureRule& rule, + std::size_t expected_num_points, + int expected_exactness) +{ + EXPECT_EQ(rule.cell_family(), svmp::CellFamily::Line); + EXPECT_EQ(rule.dimension(), 1u); + EXPECT_EQ(rule.reference_cell_measure(), 2.0); + EXPECT_EQ(rule.polynomial_exactness(), expected_exactness); + ASSERT_EQ(rule.num_points(), expected_num_points); + ASSERT_EQ(rule.points().size(), expected_num_points); + ASSERT_EQ(rule.weights().size(), expected_num_points); + + for (std::size_t point_index = 0; + point_index < rule.num_points(); + ++point_index) { + SCOPED_TRACE(::testing::Message() << "point index=" << point_index); + EXPECT_EQ(rule.point(point_index)[1], 0.0); + EXPECT_EQ(rule.point(point_index)[2], 0.0); + } +} + +void expect_line_rule_invariants( + const QuadratureRule& rule, + LineEndpointPolicy endpoint_policy) +{ + ASSERT_GT(rule.num_points(), 0u); + ASSERT_EQ(rule.points().size(), rule.weights().size()); + + for (std::size_t point_index = 0; + point_index < rule.num_points(); + ++point_index) { + SCOPED_TRACE(::testing::Message() << "point index=" << point_index); + + const double coordinate = rule.point(point_index)[0]; + const double weight = rule.weight(point_index); + EXPECT_TRUE(std::isfinite(coordinate)); + EXPECT_TRUE(std::isfinite(weight)); + EXPECT_GT(weight, 0.0); + + if (endpoint_policy == LineEndpointPolicy::Included) { + EXPECT_GE(coordinate, -1.0); + EXPECT_LE(coordinate, 1.0); + if (point_index > 0u && + point_index + 1u < rule.num_points()) { + EXPECT_GT(coordinate, -1.0); + EXPECT_LT(coordinate, 1.0); + } + } else { + EXPECT_GT(coordinate, -1.0); + EXPECT_LT(coordinate, 1.0); + } + + if (point_index > 0u) { + EXPECT_LT( + rule.point(point_index - 1u)[0], + coordinate); + } + + const std::size_t mirror_index = + rule.num_points() - 1u - point_index; + EXPECT_NEAR( + coordinate, + -rule.point(mirror_index)[0], + kStructureTolerance); + EXPECT_NEAR(weight, rule.weight(mirror_index), kStructureTolerance); + } + + if (endpoint_policy == LineEndpointPolicy::Included) { + ASSERT_GE(rule.num_points(), 2u); + EXPECT_EQ(rule.point(0)[0], -1.0); + EXPECT_EQ(rule.point(rule.num_points() - 1u)[0], 1.0); + } + + if (rule.num_points() % 2u == 1u) { + EXPECT_EQ(rule.point(rule.num_points() / 2u)[0], 0.0); + } + + const long double measure_error = std::abs( + accumulate_line_moment(rule, 0u) - + static_cast(rule.reference_cell_measure())); + EXPECT_LE(measure_error, static_cast(kStructureTolerance)); +} + +void expect_advertised_line_exactness(const QuadratureRule& rule) +{ + ASSERT_GE(rule.polynomial_exactness(), 0); + + for (int power = 0; + power <= rule.polynomial_exactness(); + ++power) { + SCOPED_TRACE(::testing::Message() << "monomial power=" << power); + const std::size_t nonnegative_power = + static_cast(power); + const long double error = std::abs( + accumulate_line_moment(rule, nonnegative_power) - + analytic_line_monomial_integral(nonnegative_power)); + EXPECT_LE(error, static_cast(kMomentTolerance)); + } +} + +void expect_every_supported_line_rule( + QuadratureRule (*generator)(int), + std::size_t first_num_points, + std::size_t last_num_points, + int exactness_subtrahend, + LineEndpointPolicy endpoint_policy) +{ + for (std::size_t num_points = first_num_points; + num_points <= last_num_points; + ++num_points) { + SCOPED_TRACE( + ::testing::Message() << "num_points=" << num_points); + const int expected_exactness = + static_cast(2 * num_points) - exactness_subtrahend; + // Both requests must select this minimum point count and report its + // actual exactness, not just echo the requested degree. + for (const int requested_exactness : + {expected_exactness - 1, expected_exactness}) { + SCOPED_TRACE( + ::testing::Message() << "requested_exactness=" << requested_exactness); + const QuadratureRule rule = generator(requested_exactness); + + expect_common_line_metadata( + rule, num_points, expected_exactness); + expect_line_rule_invariants(rule, endpoint_policy); + expect_advertised_line_exactness(rule); + } + } +} + +void expect_canonical_rule( + const QuadratureRule& rule, + int expected_exactness, + std::span expected_points, + std::span expected_weights) +{ + SCOPED_TRACE(::testing::Message() << "num_points=" << expected_points.size()); + ASSERT_EQ(rule.num_points(), expected_points.size()); + ASSERT_EQ(rule.num_points(), expected_weights.size()); + EXPECT_EQ(rule.polynomial_exactness(), expected_exactness); + + // The exhaustive sweeps cover metadata, invariants, and advertised moments. + // These fixtures independently anchor the samples and the exactness limit. + for (std::size_t point_index = 0; point_index < rule.num_points(); ++point_index) { + SCOPED_TRACE(::testing::Message() << "point index=" << point_index); + const double expected_coordinate = expected_points[point_index]; + if (expected_coordinate == -1.0 || expected_coordinate == 0.0 || + expected_coordinate == 1.0) { + EXPECT_EQ(rule.point(point_index)[0], expected_coordinate); + } else { + EXPECT_NEAR(rule.point(point_index)[0], expected_coordinate, kFixtureTolerance); + } + EXPECT_NEAR(rule.weight(point_index), expected_weights[point_index], kFixtureTolerance); + } + + const std::size_t first_unadvertised_even_power = + static_cast(expected_exactness + 1); + const long double first_unadvertised_error = std::abs( + accumulate_line_moment(rule, first_unadvertised_even_power) - + analytic_line_monomial_integral(first_unadvertised_even_power)); + EXPECT_GT(first_unadvertised_error, static_cast(kMomentTolerance)); +} + +template +void expect_invalid_argument_with_message( + Function&& function, + std::string_view expected_substring) +{ + ASSERT_FALSE(expected_substring.empty()); + + try { + std::forward(function)(); + FAIL() << "Expected InvalidArgumentException containing: " + << expected_substring; + } catch (const InvalidArgumentException& exception) { + const std::string_view actual_message{exception.what()}; + EXPECT_NE( + actual_message.find(expected_substring), + std::string_view::npos) + << "actual message: " << actual_message; + } catch (const std::exception& exception) { + FAIL() << "Received a different exception type: " << exception.what(); + } catch (...) { + FAIL() << "Received an unknown exception type"; + } +} + +} // namespace + +/** + * @test QuadratureGeneratorTestSupport.ExercisesSharedLineRuleChecks exercises + * the analytic moments, metadata, invariants, and exception-message helpers + * against explicit interior and endpoint-inclusive line rules. + */ +TEST(QuadratureGeneratorTestSupport, ExercisesSharedLineRuleChecks) +{ + const double abscissa = std::sqrt(3.0 / 5.0); + const QuadratureRule interior_rule( + svmp::CellFamily::Line, + 5, + {{-abscissa, 0.0, 0.0}, + {0.0, 0.0, 0.0}, + {abscissa, 0.0, 0.0}}, + {5.0 / 9.0, 8.0 / 9.0, 5.0 / 9.0}); + const QuadratureRule endpoint_rule( + svmp::CellFamily::Line, + 1, + {{-1.0, 0.0, 0.0}, {1.0, 0.0, 0.0}}, + {1.0, 1.0}); + + EXPECT_EQ(analytic_line_monomial_integral(0u), 2.0L); + EXPECT_EQ(analytic_line_monomial_integral(1u), 0.0L); + EXPECT_EQ( + analytic_line_monomial_integral(2u), + 2.0L / 3.0L); + EXPECT_EQ(analytic_line_monomial_integral(255u), 0.0L); + EXPECT_LE( + std::abs(accumulate_line_moment(interior_rule, 255u)), + static_cast(kMomentTolerance)); + + expect_common_line_metadata(interior_rule, 3u, 5); + expect_line_rule_invariants( + interior_rule, + LineEndpointPolicy::Excluded); + expect_advertised_line_exactness(interior_rule); + + expect_common_line_metadata(endpoint_rule, 2u, 1); + expect_line_rule_invariants( + endpoint_rule, + LineEndpointPolicy::Included); + expect_advertised_line_exactness(endpoint_rule); + expect_invalid_argument_with_message( + [] { + (void)QuadratureRule( + svmp::CellFamily::Line, 1, {}, {}); + }, + "at least one point"); +} + +/** + * @test GaussLegendreImplementation.GeneratesCanonicalLowOrderRules compares + * requests 0, 2, and 4 with analytic one-, two-, and three-point rules, including + * their weights and the actual exactness exceeding each request by one. + */ +TEST(GaussLegendreImplementation, GeneratesCanonicalLowOrderRules) +{ + expect_canonical_rule( + make_gauss_legendre_rule(0), 1, + std::array{0.0}, + std::array{2.0}); + + const double two_point_abscissa = 1.0 / std::sqrt(3.0); + expect_canonical_rule( + make_gauss_legendre_rule(2), 3, + std::array{-two_point_abscissa, two_point_abscissa}, + std::array{1.0, 1.0}); + + const double three_point_abscissa = std::sqrt(3.0 / 5.0); + expect_canonical_rule( + make_gauss_legendre_rule(4), 5, + std::array{ + -three_point_abscissa, 0.0, three_point_abscissa}, + std::array{5.0 / 9.0, 8.0 / 9.0, 5.0 / 9.0}); +} + +/** + * @test GaussLegendreImplementation.GeneratesEverySupportedRule checks every + * request from 0 through 255 for minimum point count, actual exactness, ordered + * symmetric interior samples, positive weights, and all advertised moments. + */ +TEST(GaussLegendreImplementation, GeneratesEverySupportedRule) +{ + expect_every_supported_line_rule( + &make_gauss_legendre_rule, + 1, + 128, + 1, + LineEndpointPolicy::Excluded); +} + +/** + * @test GaussLegendreImplementation.RejectsRequestsOutsideSupportedRange checks + * that negative requests, the first unsupported degree, and integer extremes + * raise InvalidArgumentException with the supported-range diagnostic. + */ +TEST(GaussLegendreImplementation, RejectsRequestsOutsideSupportedRange) +{ + constexpr std::string_view expected_message = + "requested_exactness must be in [0, 255]"; + constexpr std::array invalid_exactness{ + std::numeric_limits::min(), + -1, + 256, + std::numeric_limits::max()}; + + for (const int requested_exactness : invalid_exactness) { + SCOPED_TRACE( + ::testing::Message() << "requested_exactness=" << requested_exactness); + expect_invalid_argument_with_message( + [requested_exactness] { + (void)make_gauss_legendre_rule(requested_exactness); + }, + expected_message); + } +} + +/** + * @test GaussLobattoImplementation.GeneratesCanonicalLowOrderRules compares + * requests 0, 2, and 4 with analytic two-, three-, and four-point rules, including + * exact endpoints, weights, and the actual exactness exceeding each request. + */ +TEST(GaussLobattoImplementation, GeneratesCanonicalLowOrderRules) +{ + expect_canonical_rule( + make_gauss_lobatto_rule(0), 1, + std::array{-1.0, 1.0}, + std::array{1.0, 1.0}); + expect_canonical_rule( + make_gauss_lobatto_rule(2), 3, + std::array{-1.0, 0.0, 1.0}, + std::array{1.0 / 3.0, 4.0 / 3.0, 1.0 / 3.0}); + + const double four_point_abscissa = 1.0 / std::sqrt(5.0); + expect_canonical_rule( + make_gauss_lobatto_rule(4), 5, + std::array{ + -1.0, -four_point_abscissa, four_point_abscissa, 1.0}, + std::array{1.0 / 6.0, 5.0 / 6.0, 5.0 / 6.0, 1.0 / 6.0}); +} + +/** + * @test GaussLobattoImplementation.GeneratesEverySupportedRule checks every + * request from 0 through 253 for minimum point count, actual exactness, ordered + * symmetric samples, exact endpoints, positive weights, and advertised moments. + */ +TEST(GaussLobattoImplementation, GeneratesEverySupportedRule) +{ + expect_every_supported_line_rule( + &make_gauss_lobatto_rule, + 2, + 128, + 3, + LineEndpointPolicy::Included); +} + +/** + * @test GaussLobattoImplementation.RejectsRequestsOutsideSupportedRange checks + * that negative requests, the first unsupported degree, and integer extremes + * raise InvalidArgumentException with the supported-range diagnostic. + */ +TEST(GaussLobattoImplementation, RejectsRequestsOutsideSupportedRange) +{ + constexpr std::string_view expected_message = + "Gauss-Lobatto-Legendre generator: " + "requested_exactness must be in [0, 253]"; + constexpr std::array invalid_exactness{ + std::numeric_limits::min(), + -1, + 254, + std::numeric_limits::max()}; + + for (const int requested_exactness : invalid_exactness) { + SCOPED_TRACE( + ::testing::Message() << "requested_exactness=" << requested_exactness); + expect_invalid_argument_with_message( + [requested_exactness] { + (void)make_gauss_lobatto_rule(requested_exactness); + }, + expected_message); + } +} + +/** + * @test GaussLobattoBasisConsistency.MatchesRepresentativeNodeDistributions + * compares every supported Lobatto point count with Basis nodes, checking the + * two callers' indexing, symmetry, endpoints, and exact centers. + */ +TEST(GaussLobattoBasisConsistency, MatchesRepresentativeNodeDistributions) +{ + for (std::size_t num_points = 2; num_points <= 128; ++num_points) { + SCOPED_TRACE( + ::testing::Message() << "num_points=" << num_points); + const QuadratureRule rule = + make_gauss_lobatto_rule(static_cast(2 * num_points - 3)); + ASSERT_EQ(rule.num_points(), num_points); + + for (std::size_t point_index = 0; + point_index < rule.num_points(); + ++point_index) { + SCOPED_TRACE( + ::testing::Message() << "point index=" << point_index); + const double quadrature_coordinate = + rule.point(point_index)[0]; + const double basis_coordinate = + svmp::FE::basis::line_coord_pm_one( + static_cast(point_index), + static_cast(num_points - 1)); + + if (point_index == 0u) { + EXPECT_EQ(quadrature_coordinate, -1.0); + EXPECT_EQ(basis_coordinate, -1.0); + } else if (point_index + 1u == rule.num_points()) { + EXPECT_EQ(quadrature_coordinate, 1.0); + EXPECT_EQ(basis_coordinate, 1.0); + } else if (num_points % 2 == 1 && + point_index == rule.num_points() / 2u) { + EXPECT_EQ(quadrature_coordinate, 0.0); + EXPECT_EQ(basis_coordinate, 0.0); + } else { + EXPECT_NEAR( + quadrature_coordinate, + basis_coordinate, + kBasisConsistencyTolerance); + } + } + } +} + +TEST(GaussLobattoNodes, RefinesAnalyticNodesAndLegendreValues) +{ + const auto center = detail::refine_gauss_lobatto_node(3, 0, 1.0e-15); + EXPECT_EQ(center.coordinate, 0.0); + EXPECT_EQ(center.legendre_value, -0.5); + + const auto cubic = detail::refine_gauss_lobatto_node(4, 0, 1.0e-15); + EXPECT_NEAR(cubic.coordinate, 1.0 / std::sqrt(5.0), kFixtureTolerance); + EXPECT_NEAR(cubic.legendre_value, -1.0 / std::sqrt(5.0), kFixtureTolerance); + + const auto quartic = detail::refine_gauss_lobatto_node(5, 0, 1.0e-15); + EXPECT_NEAR(quartic.coordinate, std::sqrt(3.0 / 7.0), kFixtureTolerance); + EXPECT_NEAR(quartic.legendre_value, -3.0 / 7.0, kFixtureTolerance); + const auto quartic_center = detail::refine_gauss_lobatto_node(5, 1, 1.0e-15); + EXPECT_EQ(quartic_center.coordinate, 0.0); + EXPECT_EQ(quartic_center.legendre_value, 3.0 / 8.0); +} + +TEST(GaussLobattoNodes, RejectsInvalidRefinementRequests) +{ + for (const std::size_t num_points : {0u, 1u, 2u}) { + EXPECT_THROW((void)detail::refine_gauss_lobatto_node(num_points, 0, 1.0e-15), + InvalidArgumentException); + } + EXPECT_THROW((void)detail::refine_gauss_lobatto_node(4, 1, 1.0e-15), + InvalidArgumentException); + for (const double tolerance : {0.0, -1.0, + std::numeric_limits::infinity(), + std::numeric_limits::quiet_NaN()}) { + EXPECT_THROW((void)detail::refine_gauss_lobatto_node(4, 0, tolerance), + InvalidArgumentException); + } +}