diff --git a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp new file mode 100644 index 000000000..8065f7529 --- /dev/null +++ b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp @@ -0,0 +1,261 @@ +// 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 +#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() == static_cast(2 * kMaximumPoints - 3)); + +// 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-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. +// 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)); +} + +// Refine a nonnegative interior node with bounded, cosine-seeded Newton +// iteration on f(x) = x*P_(n-1)(x) - P_(n-2)(x). Recheck the final correction +// before forming w = 2 / (n*(n-1)*P_(n-1)(x)^2). The caller mirrors the node and +// handles endpoints; is_center assigns an odd rule's center exactly to zero. +std::pair generate_interior_root_and_weight(std::size_t num_points, + int half_root_index, bool is_center, double weight_denominator_scale) +{ + 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) > 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, 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) <= + kNewtonCorrectionTolerance)) { + raise_generation_failure( + num_points, half_root_index, iteration, final_correction, + "refined root failed final correction validation"); + } + + const double denominator = + weight_denominator_scale * final_polynomial_value * + final_polynomial_value; + if (!(std::isfinite(denominator) && denominator > 0.0)) { + raise_generation_failure( + num_points, half_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, half_root_index, iteration, weight, + "refined root produced an invalid quadrature weight"); + } + + return {root, weight}; + } + + raise_generation_failure( + num_points, half_root_index, kMaximumNewtonIterations, correction, + "Newton refinement did not converge"); +} + +} // 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 = static_cast(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, static_cast(half_root_index), left_index == right_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, static_cast(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 <= + static_cast(kRuleValidationTolerance))) { + raise_generation_failure( + num_points, -1, -1, static_cast(measure_error), + "generated weights do not reproduce the reference measure"); + } + + const int polynomial_exactness = static_cast(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..0750365cf --- /dev/null +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp @@ -0,0 +1,235 @@ +// 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() == static_cast(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, int 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 = static_cast(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, static_cast(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, static_cast(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 <= + static_cast(kRuleValidationTolerance))) { + raise_generation_failure( + num_points, -1, -1, static_cast(measure_error), + "generated weights do not reproduce the reference measure"); + } + + const int polynomial_exactness = static_cast(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/Quadrature/test_QuadratureGenerators.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp new file mode 100644 index 000000000..53d100f62 --- /dev/null +++ b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp @@ -0,0 +1,516 @@ +// 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/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 independently refine computed interior GLL roots, 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 2-, 4-, 65-, and 128-point Lobatto rules with the independently + * computed Basis nodes, requiring exact endpoints and close interior values. + */ +TEST(GaussLobattoBasisConsistency, MatchesRepresentativeNodeDistributions) +{ + constexpr std::array point_counts{2, 4, 65, 128}; + + for (const std::size_t num_points : point_counts) { + 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); + } + } + } +}