From cdcdb32090e00b6b337fb08decd80d7a00411653 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Tue, 25 Aug 2026 12:19:30 -0700 Subject: [PATCH 01/27] Add quadrature generator test support --- .../Quadrature/test_QuadratureGenerators.cpp | 241 ++++++++++++++++++ 1 file changed, 241 insertions(+) create mode 100644 tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp new file mode 100644 index 000000000..c92bb93ef --- /dev/null +++ b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp @@ -0,0 +1,241 @@ +// 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 Shared test support for bounded one-dimensional quadrature generators. + */ + +#include + +#include "FE/Common/FEException.h" +#include "FE/Quadrature/QuadratureRule.h" + +#include +#include +#include +#include +#include +#include + +using namespace svmp::FE; +using namespace svmp::FE::quadrature; + +namespace { + +constexpr double kStructureTolerance = 1.0e-12; +constexpr double kMomentTolerance = 2.0e-12; + +enum class LineEndpointPolicy { + Excluded, + Included, +}; + +double analytic_line_monomial_integral(std::size_t power) +{ + if (power % 2u != 0u) { + return 0.0; + } + return 2.0 / (static_cast(power) + 1.0); +} + +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 static_cast(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_DOUBLE_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_DOUBLE_EQ(rule.point(point_index)[1], 0.0); + EXPECT_DOUBLE_EQ(rule.point(point_index)[2], 0.0); + } +} + +void expect_line_rule_invariants( + const QuadratureRule& rule, + LineEndpointPolicy endpoint_policy, + double tolerance = kStructureTolerance) +{ + 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], + tolerance); + EXPECT_NEAR(weight, rule.weight(mirror_index), tolerance); + } + + if (endpoint_policy == LineEndpointPolicy::Included) { + ASSERT_GE(rule.num_points(), 2u); + EXPECT_DOUBLE_EQ(rule.point(0)[0], -1.0); + EXPECT_DOUBLE_EQ(rule.point(rule.num_points() - 1u)[0], 1.0); + } + + if (rule.num_points() % 2u == 1u) { + EXPECT_DOUBLE_EQ(rule.point(rule.num_points() / 2u)[0], 0.0); + } + + EXPECT_NEAR( + accumulate_line_moment(rule, 0u), + rule.reference_cell_measure(), + tolerance); +} + +void expect_advertised_line_exactness( + const QuadratureRule& rule, + double tolerance = kMomentTolerance) +{ + 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); + EXPECT_NEAR( + accumulate_line_moment(rule, nonnegative_power), + analytic_line_monomial_integral(nonnegative_power), + tolerance); + } +} + +template +void expect_exception_with_message( + Function&& function, + std::string_view expected_substring) +{ + static_assert(std::is_base_of_v); + ASSERT_FALSE(expected_substring.empty()); + + try { + std::forward(function)(); + FAIL() << "Expected requested exception containing: " + << expected_substring; + } catch (const ExceptionType& 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) +{ + 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_DOUBLE_EQ(analytic_line_monomial_integral(0u), 2.0); + EXPECT_DOUBLE_EQ(analytic_line_monomial_integral(1u), 0.0); + EXPECT_DOUBLE_EQ(analytic_line_monomial_integral(2u), 2.0 / 3.0); + EXPECT_DOUBLE_EQ(analytic_line_monomial_integral(255u), 0.0); + EXPECT_NEAR( + accumulate_line_moment(interior_rule, 255u), + 0.0, + 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_exception_with_message( + [] { + (void)QuadratureRule( + svmp::CellFamily::Line, 1, {}, {}); + }, + "at least one point"); +} From 5e069792a7e3d1dc6920285fa148c129510eda48 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Tue, 25 Aug 2026 13:10:21 -0700 Subject: [PATCH 02/27] Add Gauss-Legendre public contract --- .../solver/FE/Quadrature/GaussQuadrature.h | 53 +++++++++++++++++++ .../Quadrature/test_QuadratureGenerators.cpp | 7 +++ 2 files changed, 60 insertions(+) create mode 100644 Code/Source/solver/FE/Quadrature/GaussQuadrature.h diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.h b/Code/Source/solver/FE/Quadrature/GaussQuadrature.h new file mode 100644 index 000000000..052875ec0 --- /dev/null +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.h @@ -0,0 +1,53 @@ +// 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 point count. + * @return The inclusive point-count limit, 128. + */ +[[nodiscard]] constexpr int max_gauss_legendre_points() noexcept +{ + return 128; +} + +/** + * @brief Generate an @p num_points Gauss-Legendre rule on @f$[-1,1]@f$. + * + * @details The returned line rule contains the roots of @f$P_n@f$ in strictly + * increasing order, where @f$n@f$ is @p num_points. Every point lies strictly + * inside @f$(-1,1)@f$, so neither endpoint is included. The rule has polynomial + * exactness @f$2n-1@f$ and positive weights aligned with its points. + * + * @param num_points Number of quadrature points; must be in + * @f$[1,\texttt{max\_gauss\_legendre\_points()}]@f$. + * @return A complete QuadratureRule value for CellFamily::Line. + * @throws InvalidArgumentException If @p num_points is outside the supported + * range. + * @throws ConvergenceException If root refinement or final numerical + * validation fails. + */ +[[nodiscard]] QuadratureRule +make_gauss_legendre_rule(int num_points); + +/** @} */ + +} // namespace svmp::FE::quadrature + +#endif // SVMP_FE_GAUSS_QUADRATURE_H diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp index c92bb93ef..562980f77 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp @@ -9,6 +9,7 @@ #include #include "FE/Common/FEException.h" +#include "FE/Quadrature/GaussQuadrature.h" #include "FE/Quadrature/QuadratureRule.h" #include @@ -26,6 +27,12 @@ namespace { constexpr double kStructureTolerance = 1.0e-12; constexpr double kMomentTolerance = 2.0e-12; +static_assert(max_gauss_legendre_points() == 128); +static_assert(noexcept(max_gauss_legendre_points())); +static_assert( + std::is_same_v); + enum class LineEndpointPolicy { Excluded, Included, From 24d3ffdcdd6369c31af1950ca48f9d098436c4f3 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Tue, 25 Aug 2026 14:30:40 -0700 Subject: [PATCH 03/27] Document Gauss-Legendre point limit --- Code/Source/solver/FE/Quadrature/GaussQuadrature.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.h b/Code/Source/solver/FE/Quadrature/GaussQuadrature.h index 052875ec0..14cdd376d 100644 --- a/Code/Source/solver/FE/Quadrature/GaussQuadrature.h +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.h @@ -20,6 +20,9 @@ namespace svmp::FE::quadrature { /** * @brief Return the largest supported Gauss-Legendre point count. + * @details The 128-point ceiling is a project support bound that limits + * generator work and downstream product-rule growth while providing line + * exactness through degree 255. * @return The inclusive point-count limit, 128. */ [[nodiscard]] constexpr int max_gauss_legendre_points() noexcept From 2f1006a418906a6e3add2eb689edf83b8d183e1e Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Tue, 25 Aug 2026 14:47:45 -0700 Subject: [PATCH 04/27] Implement bounded Gauss-Legendre generator --- .../solver/FE/Quadrature/GaussQuadrature.cpp | 462 ++++++++++++++++++ .../Quadrature/test_QuadratureGenerators.cpp | 39 ++ 2 files changed, 501 insertions(+) create mode 100644 Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp new file mode 100644 index 000000000..be2f1cb75 --- /dev/null +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp @@ -0,0 +1,462 @@ +// 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 Bounded generation and validation of Gauss-Legendre line rules. + * @ingroup FE_Quadrature + */ + +#include "FE/Quadrature/GaussQuadrature.h" + +#include "FE/Common/FEException.h" + +#include +#include +#include +#include +#include +#include + +namespace svmp::FE::quadrature { +namespace { + +constexpr int kMaximumNewtonIterations = 100; +constexpr double kNewtonCorrectionTolerance = 1.0e-14; +constexpr double kRuleValidationTolerance = 1.0e-12; + +struct LegendreEvaluation { + double value; + double derivative; +}; + +LegendreEvaluation evaluate_legendre_with_derivative( + int degree, + double coordinate) noexcept +{ + double previous_value = 1.0; + double previous_derivative = 0.0; + + if (degree == 0) { + return {previous_value, previous_derivative}; + } + + double value = coordinate; + double derivative = 1.0; + for (int recurrence_degree = 2; + recurrence_degree <= degree; + ++recurrence_degree) { + const double degree_value = + static_cast(recurrence_degree); + const double recurrence_factor = + static_cast(2 * recurrence_degree - 1); + const double next_value = + (recurrence_factor * coordinate * value - + static_cast(recurrence_degree - 1) * previous_value) / + degree_value; + const double next_derivative = + (recurrence_factor * (value + coordinate * derivative) - + static_cast(recurrence_degree - 1) * + previous_derivative) / + degree_value; + + previous_value = value; + previous_derivative = derivative; + value = next_value; + derivative = next_derivative; + } + + return {value, derivative}; +} + +[[noreturn]] void raise_generation_failure( + int num_points, + int root_index, + int iteration, + std::string_view diagnostic_name, + double diagnostic_value, + std::string_view detail) +{ + std::ostringstream message; + message << "Gauss-Legendre generator: " << detail + << ", num_points=" << num_points + << ", root_index=" << root_index + << ", " << diagnostic_name << '=' << diagnostic_value; + + const double residual = std::isfinite(diagnostic_value) + ? std::abs(diagnostic_value) + : 0.0; + svmp::raise( + message.str(), iteration, residual); +} + +void validate_num_points(int num_points) +{ + if (num_points < 1 || + num_points > max_gauss_legendre_points()) { + std::ostringstream message; + message << "Gauss-Legendre generator: num_points must be in [1, " + << max_gauss_legendre_points() << ']'; + svmp::raise(message.str()); + } +} + +void validate_generated_rule( + int num_points, + const std::vector& points, + const std::vector& weights) +{ + const std::size_t expected_size = + static_cast(num_points); + if (points.size() != expected_size) { + raise_generation_failure( + num_points, + -1, + -1, + "generated_point_count", + static_cast(points.size()), + "generated point storage has the wrong size"); + } + if (weights.size() != expected_size) { + raise_generation_failure( + num_points, + -1, + -1, + "generated_weight_count", + static_cast(weights.size()), + "generated weight storage has the wrong size"); + } + + long double weight_sum = 0.0L; + long double weight_sum_correction = 0.0L; + for (std::size_t point_index = 0; + point_index < expected_size; + ++point_index) { + const QuadPoint& point = points[point_index]; + const double coordinate = point[0]; + const double weight = weights[point_index]; + const int root_index = static_cast(point_index); + + if (!std::isfinite(coordinate)) { + raise_generation_failure( + num_points, + root_index, + -1, + "coordinate", + coordinate, + "generated a non-finite point"); + } + if (point[1] != 0.0 || point[2] != 0.0) { + const double inactive_coordinate = + point[1] != 0.0 ? point[1] : point[2]; + raise_generation_failure( + num_points, + root_index, + -1, + "inactive_coordinate", + inactive_coordinate, + "generated a nonzero inactive coordinate"); + } + if (!std::isfinite(weight) || weight <= 0.0) { + raise_generation_failure( + num_points, + root_index, + -1, + "weight", + weight, + "generated a non-finite or non-positive weight"); + } + if (coordinate <= -1.0 || coordinate >= 1.0) { + raise_generation_failure( + num_points, + root_index, + -1, + "coordinate", + coordinate, + "generated a point outside the open reference interval"); + } + if (point_index > 0u) { + const double spacing = + coordinate - points[point_index - 1u][0]; + if (!std::isfinite(spacing) || spacing <= 0.0) { + raise_generation_failure( + num_points, + root_index, + -1, + "point_spacing", + spacing, + "generated points are not strictly increasing"); + } + } + + const std::size_t mirror_index = + expected_size - 1u - point_index; + const double point_symmetry_error = + std::abs(coordinate + points[mirror_index][0]); + if (!std::isfinite(point_symmetry_error) || + point_symmetry_error > kRuleValidationTolerance) { + raise_generation_failure( + num_points, + root_index, + -1, + "point_symmetry_error", + point_symmetry_error, + "generated points are not symmetric"); + } + const double weight_symmetry_error = + std::abs(weight - weights[mirror_index]); + if (!std::isfinite(weight_symmetry_error) || + weight_symmetry_error > kRuleValidationTolerance) { + raise_generation_failure( + num_points, + root_index, + -1, + "weight_symmetry_error", + weight_symmetry_error, + "generated weights are not symmetric"); + } + + const long double weight_term = static_cast(weight); + const long double next_weight_sum = weight_sum + weight_term; + if (std::abs(weight_sum) >= std::abs(weight_term)) { + weight_sum_correction += + (weight_sum - next_weight_sum) + weight_term; + } else { + weight_sum_correction += + (weight_term - next_weight_sum) + weight_sum; + } + weight_sum = next_weight_sum; + } + + if (expected_size % 2u == 1u) { + const std::size_t center_index = expected_size / 2u; + if (points[center_index][0] != 0.0) { + raise_generation_failure( + num_points, + static_cast(center_index), + -1, + "center_coordinate", + points[center_index][0], + "generated an inexact center point"); + } + } + + const long double corrected_weight_sum = + weight_sum + weight_sum_correction; + const long double measure_error = + std::abs(corrected_weight_sum - 2.0L); + if (!std::isfinite(corrected_weight_sum) || + measure_error > + static_cast(kRuleValidationTolerance)) { + raise_generation_failure( + num_points, + -1, + -1, + "measure_error", + static_cast(measure_error), + "generated weights do not reproduce the reference measure"); + } +} + +} // namespace + +QuadratureRule make_gauss_legendre_rule(int num_points) +{ + validate_num_points(num_points); + + const std::size_t point_count = + static_cast(num_points); + std::vector points( + point_count, QuadPoint::Zero()); + std::vector weights(point_count); + + const double pi = std::acos(-1.0); + const int roots_to_refine = (num_points + 1) / 2; + for (int root_index = 0; + root_index < roots_to_refine; + ++root_index) { + double root = std::cos( + pi * (static_cast(root_index) + 0.75) / + (static_cast(num_points) + 0.5)); + if (!std::isfinite(root)) { + raise_generation_failure( + num_points, + root_index, + -1, + "initial_root", + root, + "computed a non-finite asymptotic root seed"); + } + + bool converged = false; + double correction = 0.0; + int iterations_used = 0; + for (int iteration = 1; + iteration <= kMaximumNewtonIterations; + ++iteration) { + iterations_used = iteration; + const LegendreEvaluation evaluation = + evaluate_legendre_with_derivative(num_points, root); + if (!std::isfinite(evaluation.value)) { + raise_generation_failure( + num_points, + root_index, + iteration, + "polynomial_value", + evaluation.value, + "encountered a non-finite Legendre value"); + } + if (!std::isfinite(evaluation.derivative) || + evaluation.derivative == 0.0) { + raise_generation_failure( + num_points, + root_index, + iteration, + "polynomial_derivative", + evaluation.derivative, + "encountered an invalid Legendre derivative"); + } + + correction = evaluation.value / evaluation.derivative; + if (!std::isfinite(correction)) { + raise_generation_failure( + num_points, + root_index, + iteration, + "newton_correction", + correction, + "computed a non-finite Newton correction"); + } + + const double updated_root = root - correction; + if (!std::isfinite(updated_root)) { + raise_generation_failure( + num_points, + root_index, + iteration, + "updated_root", + updated_root, + "computed a non-finite Newton update"); + } + root = updated_root; + + if (std::abs(correction) <= + kNewtonCorrectionTolerance) { + converged = true; + break; + } + } + + if (!converged) { + raise_generation_failure( + num_points, + root_index, + kMaximumNewtonIterations, + "newton_correction", + correction, + "Newton refinement did not converge"); + } + + const std::size_t left_index = + static_cast(root_index); + const std::size_t right_index = + point_count - 1u - left_index; + if (left_index == right_index) { + root = 0.0; + } + if (!std::isfinite(root) || root < 0.0 || root >= 1.0 || + (left_index != right_index && root == 0.0)) { + raise_generation_failure( + num_points, + root_index, + iterations_used, + "refined_root", + root, + "refined root is outside the expected half interval"); + } + + const LegendreEvaluation final_evaluation = + evaluate_legendre_with_derivative(num_points, root); + if (!std::isfinite(final_evaluation.value)) { + raise_generation_failure( + num_points, + root_index, + iterations_used, + "final_polynomial_value", + final_evaluation.value, + "refined root produced a non-finite Legendre value"); + } + if (!std::isfinite(final_evaluation.derivative) || + final_evaluation.derivative == 0.0) { + raise_generation_failure( + num_points, + root_index, + iterations_used, + "final_polynomial_derivative", + final_evaluation.derivative, + "refined root produced an invalid Legendre derivative"); + } + + const double final_correction = + final_evaluation.value / final_evaluation.derivative; + if (!std::isfinite(final_correction) || + std::abs(final_correction) > + kNewtonCorrectionTolerance) { + raise_generation_failure( + num_points, + root_index, + iterations_used, + "final_correction", + final_correction, + "refined root failed final correction validation"); + } + + const double interval_factor = + (1.0 - root) * (1.0 + root); + const double denominator = + interval_factor * final_evaluation.derivative * + final_evaluation.derivative; + if (!std::isfinite(interval_factor) || interval_factor <= 0.0 || + !std::isfinite(denominator) || denominator <= 0.0) { + raise_generation_failure( + num_points, + root_index, + iterations_used, + "weight_denominator", + 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, + iterations_used, + "weight", + weight, + "refined root produced a non-finite or non-positive weight"); + } + + if (left_index == right_index) { + points[left_index][0] = 0.0; + weights[left_index] = weight; + } else { + points[left_index][0] = -root; + points[right_index][0] = root; + weights[left_index] = weight; + weights[right_index] = weight; + } + } + + validate_generated_rule(num_points, points, weights); + + 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/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp index 562980f77..dcb842ad1 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp @@ -12,6 +12,7 @@ #include "FE/Quadrature/GaussQuadrature.h" #include "FE/Quadrature/QuadratureRule.h" +#include #include #include #include @@ -246,3 +247,41 @@ TEST(QuadratureGeneratorTestSupport, ExercisesSharedLineRuleChecks) }, "at least one point"); } + +TEST(GaussLegendreImplementation, GeneratesRepresentativeSupportedRules) +{ + const std::array point_counts{ + 1, 2, 17, max_gauss_legendre_points()}; + + for (const int num_points : point_counts) { + SCOPED_TRACE(::testing::Message() + << "num_points=" << num_points); + const QuadratureRule rule = + make_gauss_legendre_rule(num_points); + + expect_common_line_metadata( + rule, + static_cast(num_points), + 2 * num_points - 1); + expect_line_rule_invariants( + rule, + LineEndpointPolicy::Excluded); + expect_advertised_line_exactness(rule); + } +} + +TEST(GaussLegendreImplementation, RejectsRequestsOutsideSupportedRange) +{ + constexpr std::string_view expected_message = + "num_points must be in [1, 128]"; + + expect_exception_with_message( + [] { (void)make_gauss_legendre_rule(0); }, + expected_message); + expect_exception_with_message( + [] { + (void)make_gauss_legendre_rule( + max_gauss_legendre_points() + 1); + }, + expected_message); +} From b92b88eea2c9c662e209d92276d01fcea7d81da1 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Tue, 25 Aug 2026 23:32:51 -0700 Subject: [PATCH 05/27] Refine Gauss-Legendre generator internals --- .../solver/FE/Quadrature/GaussQuadrature.cpp | 463 +++++------------- 1 file changed, 135 insertions(+), 328 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp index be2f1cb75..766228570 100644 --- a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp @@ -13,6 +13,9 @@ #include #include +#include +#include +#include #include #include #include @@ -22,15 +25,13 @@ namespace svmp::FE::quadrature { namespace { constexpr int kMaximumNewtonIterations = 100; -constexpr double kNewtonCorrectionTolerance = 1.0e-14; -constexpr double kRuleValidationTolerance = 1.0e-12; +// Accommodate rounding in the Legendre recurrence and Newton update. +constexpr double kNewtonCorrectionTolerance = 64.0 * std::numeric_limits::epsilon(); +// Allow accumulated rounding across the largest supported rule. +constexpr double kRuleValidationTolerance = 32.0 * static_cast(max_gauss_legendre_points()) * + std::numeric_limits::epsilon(); -struct LegendreEvaluation { - double value; - double derivative; -}; - -LegendreEvaluation evaluate_legendre_with_derivative( +std::pair evaluate_legendre_with_derivative( int degree, double coordinate) noexcept { @@ -43,9 +44,7 @@ LegendreEvaluation evaluate_legendre_with_derivative( double value = coordinate; double derivative = 1.0; - for (int recurrence_degree = 2; - recurrence_degree <= degree; - ++recurrence_degree) { + for (int recurrence_degree = 2; recurrence_degree <= degree; ++recurrence_degree) { const double degree_value = static_cast(recurrence_degree); const double recurrence_factor = @@ -73,7 +72,6 @@ LegendreEvaluation evaluate_legendre_with_derivative( int num_points, int root_index, int iteration, - std::string_view diagnostic_name, double diagnostic_value, std::string_view detail) { @@ -81,13 +79,26 @@ LegendreEvaluation evaluate_legendre_with_derivative( message << "Gauss-Legendre generator: " << detail << ", num_points=" << num_points << ", root_index=" << root_index - << ", " << diagnostic_name << '=' << diagnostic_value; + << ", diagnostic_value=" << diagnostic_value; const double residual = std::isfinite(diagnostic_value) ? std::abs(diagnostic_value) : 0.0; - svmp::raise( - message.str(), iteration, residual); + svmp::raise(message.str(), iteration, residual); +} + +void require_generation( + bool condition, + int num_points, + int root_index, + int iteration, + double diagnostic_value, + std::string_view detail) +{ + if (!condition) { + raise_generation_failure( + num_points, root_index, iteration, diagnostic_value, detail); + } } void validate_num_points(int num_points) @@ -101,161 +112,117 @@ void validate_num_points(int num_points) } } -void validate_generated_rule( +std::pair generate_root_and_weight( int num_points, - const std::vector& points, - const std::vector& weights) + int root_index, + bool is_center) { - const std::size_t expected_size = - static_cast(num_points); - if (points.size() != expected_size) { - raise_generation_failure( - num_points, - -1, - -1, - "generated_point_count", - static_cast(points.size()), - "generated point storage has the wrong size"); - } - if (weights.size() != expected_size) { - raise_generation_failure( - num_points, - -1, - -1, - "generated_weight_count", - static_cast(weights.size()), - "generated weight storage has the wrong size"); - } - - long double weight_sum = 0.0L; - long double weight_sum_correction = 0.0L; - for (std::size_t point_index = 0; - point_index < expected_size; - ++point_index) { - const QuadPoint& point = points[point_index]; - const double coordinate = point[0]; - const double weight = weights[point_index]; - const int root_index = static_cast(point_index); - - if (!std::isfinite(coordinate)) { - raise_generation_failure( - num_points, - root_index, - -1, - "coordinate", - coordinate, - "generated a non-finite point"); - } - if (point[1] != 0.0 || point[2] != 0.0) { - const double inactive_coordinate = - point[1] != 0.0 ? point[1] : point[2]; - raise_generation_failure( - num_points, - root_index, - -1, - "inactive_coordinate", - inactive_coordinate, - "generated a nonzero inactive coordinate"); - } - if (!std::isfinite(weight) || weight <= 0.0) { - raise_generation_failure( - num_points, - root_index, - -1, - "weight", - weight, - "generated a non-finite or non-positive weight"); - } - if (coordinate <= -1.0 || coordinate >= 1.0) { - raise_generation_failure( - num_points, - root_index, - -1, - "coordinate", - coordinate, - "generated a point outside the open reference interval"); - } - if (point_index > 0u) { - const double spacing = - coordinate - points[point_index - 1u][0]; - if (!std::isfinite(spacing) || spacing <= 0.0) { - raise_generation_failure( - num_points, - root_index, - -1, - "point_spacing", - spacing, - "generated points are not strictly increasing"); - } + const double pi = std::numbers::pi_v; + double root = std::cos( + pi * (static_cast(root_index) + 0.75) / + (static_cast(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); + require_generation( + std::isfinite(polynomial_value) && + std::isfinite(polynomial_derivative) && + polynomial_derivative != 0.0, + 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; + require_generation( + std::isfinite(correction) && std::isfinite(updated_root), + num_points, root_index, iteration, correction, + "computed an invalid Newton update"); + root = updated_root; + + if (std::abs(correction) > kNewtonCorrectionTolerance) { + continue; } - const std::size_t mirror_index = - expected_size - 1u - point_index; - const double point_symmetry_error = - std::abs(coordinate + points[mirror_index][0]); - if (!std::isfinite(point_symmetry_error) || - point_symmetry_error > kRuleValidationTolerance) { - raise_generation_failure( - num_points, - root_index, - -1, - "point_symmetry_error", - point_symmetry_error, - "generated points are not symmetric"); - } - const double weight_symmetry_error = - std::abs(weight - weights[mirror_index]); - if (!std::isfinite(weight_symmetry_error) || - weight_symmetry_error > kRuleValidationTolerance) { - raise_generation_failure( - num_points, - root_index, - -1, - "weight_symmetry_error", - weight_symmetry_error, - "generated weights are not symmetric"); + if (is_center) { + root = 0.0; } + require_generation( + root >= 0.0 && root < 1.0 && + (is_center || root > 0.0), + 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); + require_generation( + std::isfinite(final_polynomial_value) && + std::isfinite(final_polynomial_derivative) && + final_polynomial_derivative != 0.0, + num_points, root_index, iteration, final_polynomial_derivative, + "refined root produced an invalid Legendre value or derivative"); - const long double weight_term = static_cast(weight); - const long double next_weight_sum = weight_sum + weight_term; - if (std::abs(weight_sum) >= std::abs(weight_term)) { - weight_sum_correction += - (weight_sum - next_weight_sum) + weight_term; - } else { - weight_sum_correction += - (weight_term - next_weight_sum) + weight_sum; - } - weight_sum = next_weight_sum; + const double final_correction = + final_polynomial_value / final_polynomial_derivative; + require_generation( + std::isfinite(final_correction) && + std::abs(final_correction) <= + kNewtonCorrectionTolerance, + 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; + require_generation( + std::isfinite(denominator) && denominator > 0.0, + num_points, root_index, iteration, denominator, + "refined root produced an invalid weight denominator"); + + const double weight = 2.0 / denominator; + require_generation( + std::isfinite(weight) && weight > 0.0, + num_points, root_index, iteration, weight, + "refined root produced an invalid quadrature weight"); + + return {root, weight}; } - if (expected_size % 2u == 1u) { - const std::size_t center_index = expected_size / 2u; - if (points[center_index][0] != 0.0) { - raise_generation_failure( - num_points, - static_cast(center_index), - -1, - "center_coordinate", - points[center_index][0], - "generated an inexact center point"); - } + raise_generation_failure( + num_points, root_index, kMaximumNewtonIterations, correction, + "Newton refinement did not converge"); +} + +void validate_ordering_and_measure( + int num_points, + const std::vector& points, + const std::vector& weights) +{ + 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]; + require_generation( + spacing > 0.0, + num_points, static_cast(point_index), -1, spacing, + "generated points are not strictly increasing"); } - const long double corrected_weight_sum = - weight_sum + weight_sum_correction; + const long double weight_sum = + std::accumulate(weights.begin(), weights.end(), 0.0L); const long double measure_error = - std::abs(corrected_weight_sum - 2.0L); - if (!std::isfinite(corrected_weight_sum) || - measure_error > - static_cast(kRuleValidationTolerance)) { - raise_generation_failure( - num_points, - -1, - -1, - "measure_error", - static_cast(measure_error), - "generated weights do not reproduce the reference measure"); - } + std::abs(weight_sum - 2.0L); + require_generation( + std::isfinite(weight_sum) && + measure_error <= + static_cast(kRuleValidationTolerance), + num_points, -1, -1, static_cast(measure_error), + "generated weights do not reproduce the reference measure"); } } // namespace @@ -270,186 +237,26 @@ QuadratureRule make_gauss_legendre_rule(int num_points) point_count, QuadPoint::Zero()); std::vector weights(point_count); - const double pi = std::acos(-1.0); const int roots_to_refine = (num_points + 1) / 2; for (int root_index = 0; root_index < roots_to_refine; ++root_index) { - double root = std::cos( - pi * (static_cast(root_index) + 0.75) / - (static_cast(num_points) + 0.5)); - if (!std::isfinite(root)) { - raise_generation_failure( - num_points, - root_index, - -1, - "initial_root", - root, - "computed a non-finite asymptotic root seed"); - } - - bool converged = false; - double correction = 0.0; - int iterations_used = 0; - for (int iteration = 1; - iteration <= kMaximumNewtonIterations; - ++iteration) { - iterations_used = iteration; - const LegendreEvaluation evaluation = - evaluate_legendre_with_derivative(num_points, root); - if (!std::isfinite(evaluation.value)) { - raise_generation_failure( - num_points, - root_index, - iteration, - "polynomial_value", - evaluation.value, - "encountered a non-finite Legendre value"); - } - if (!std::isfinite(evaluation.derivative) || - evaluation.derivative == 0.0) { - raise_generation_failure( - num_points, - root_index, - iteration, - "polynomial_derivative", - evaluation.derivative, - "encountered an invalid Legendre derivative"); - } - - correction = evaluation.value / evaluation.derivative; - if (!std::isfinite(correction)) { - raise_generation_failure( - num_points, - root_index, - iteration, - "newton_correction", - correction, - "computed a non-finite Newton correction"); - } - - const double updated_root = root - correction; - if (!std::isfinite(updated_root)) { - raise_generation_failure( - num_points, - root_index, - iteration, - "updated_root", - updated_root, - "computed a non-finite Newton update"); - } - root = updated_root; - - if (std::abs(correction) <= - kNewtonCorrectionTolerance) { - converged = true; - break; - } - } - - if (!converged) { - raise_generation_failure( - num_points, - root_index, - kMaximumNewtonIterations, - "newton_correction", - correction, - "Newton refinement did not converge"); - } - const std::size_t left_index = static_cast(root_index); const std::size_t right_index = point_count - 1u - left_index; - if (left_index == right_index) { - root = 0.0; - } - if (!std::isfinite(root) || root < 0.0 || root >= 1.0 || - (left_index != right_index && root == 0.0)) { - raise_generation_failure( - num_points, - root_index, - iterations_used, - "refined_root", - root, - "refined root is outside the expected half interval"); - } - - const LegendreEvaluation final_evaluation = - evaluate_legendre_with_derivative(num_points, root); - if (!std::isfinite(final_evaluation.value)) { - raise_generation_failure( - num_points, - root_index, - iterations_used, - "final_polynomial_value", - final_evaluation.value, - "refined root produced a non-finite Legendre value"); - } - if (!std::isfinite(final_evaluation.derivative) || - final_evaluation.derivative == 0.0) { - raise_generation_failure( - num_points, - root_index, - iterations_used, - "final_polynomial_derivative", - final_evaluation.derivative, - "refined root produced an invalid Legendre derivative"); - } - - const double final_correction = - final_evaluation.value / final_evaluation.derivative; - if (!std::isfinite(final_correction) || - std::abs(final_correction) > - kNewtonCorrectionTolerance) { - raise_generation_failure( - num_points, - root_index, - iterations_used, - "final_correction", - final_correction, - "refined root failed final correction validation"); - } - - const double interval_factor = - (1.0 - root) * (1.0 + root); - const double denominator = - interval_factor * final_evaluation.derivative * - final_evaluation.derivative; - if (!std::isfinite(interval_factor) || interval_factor <= 0.0 || - !std::isfinite(denominator) || denominator <= 0.0) { - raise_generation_failure( - num_points, - root_index, - iterations_used, - "weight_denominator", - 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, - iterations_used, - "weight", - weight, - "refined root produced a non-finite or non-positive weight"); - } + const auto [root, weight] = generate_root_and_weight( + num_points, + root_index, + left_index == right_index); - if (left_index == right_index) { - points[left_index][0] = 0.0; - weights[left_index] = weight; - } else { - points[left_index][0] = -root; - points[right_index][0] = root; - weights[left_index] = weight; - weights[right_index] = weight; - } + points[left_index][0] = -root; + points[right_index][0] = root; + weights[left_index] = weight; + weights[right_index] = weight; } - validate_generated_rule(num_points, points, weights); + validate_ordering_and_measure(num_points, points, weights); const int polynomial_exactness = 2 * num_points - 1; return QuadratureRule( From 63b84230e71c00013f2cc15e9dc53fbc97a0e5e4 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Tue, 25 Aug 2026 23:33:20 -0700 Subject: [PATCH 06/27] Complete Gauss-Legendre generator tests --- .../Quadrature/test_QuadratureGenerators.cpp | 279 +++++++++++++++--- 1 file changed, 237 insertions(+), 42 deletions(-) diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp index dcb842ad1..069d1aa2b 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include @@ -27,6 +28,8 @@ namespace { constexpr double kStructureTolerance = 1.0e-12; constexpr double kMomentTolerance = 2.0e-12; +constexpr double kFixtureTolerance = + 64.0 * std::numeric_limits::epsilon(); static_assert(max_gauss_legendre_points() == 128); static_assert(noexcept(max_gauss_legendre_points())); @@ -39,15 +42,17 @@ enum class LineEndpointPolicy { Included, }; -double analytic_line_monomial_integral(std::size_t power) +long double analytic_line_monomial_integral(std::size_t power) { if (power % 2u != 0u) { - return 0.0; + return 0.0L; } - return 2.0 / (static_cast(power) + 1.0); + return 2.0L / (static_cast(power) + 1.0L); } -double accumulate_line_moment(const QuadratureRule& rule, std::size_t power) +long double accumulate_line_moment( + const QuadratureRule& rule, + std::size_t power) { long double sum = 0.0L; long double correction = 0.0L; @@ -70,7 +75,7 @@ double accumulate_line_moment(const QuadratureRule& rule, std::size_t power) sum = next_sum; } - return static_cast(sum + correction); + return sum + correction; } void expect_common_line_metadata( @@ -152,17 +157,23 @@ void expect_line_rule_invariants( EXPECT_DOUBLE_EQ(rule.point(rule.num_points() / 2u)[0], 0.0); } - EXPECT_NEAR( - accumulate_line_moment(rule, 0u), - rule.reference_cell_measure(), - tolerance); + const long double measure_error = std::abs( + accumulate_line_moment(rule, 0u) - + static_cast(rule.reference_cell_measure())); + EXPECT_LE(measure_error, static_cast(tolerance)); } -void expect_advertised_line_exactness( +std::pair expect_advertised_line_exactness( const QuadratureRule& rule, double tolerance = kMomentTolerance) { - ASSERT_GE(rule.polynomial_exactness(), 0); + EXPECT_GE(rule.polynomial_exactness(), 0); + if (rule.polynomial_exactness() < 0) { + return {std::numeric_limits::infinity(), -1}; + } + + long double worst_error = 0.0L; + int worst_power = 0; for (int power = 0; power <= rule.polynomial_exactness(); @@ -170,11 +181,17 @@ void expect_advertised_line_exactness( SCOPED_TRACE(::testing::Message() << "monomial power=" << power); const std::size_t nonnegative_power = static_cast(power); - EXPECT_NEAR( - accumulate_line_moment(rule, nonnegative_power), - analytic_line_monomial_integral(nonnegative_power), - tolerance); + const long double error = std::abs( + accumulate_line_moment(rule, nonnegative_power) - + analytic_line_monomial_integral(nonnegative_power)); + if (error > worst_error) { + worst_error = error; + worst_power = power; + } + EXPECT_LE(error, static_cast(tolerance)); } + + return {worst_error, worst_power}; } template @@ -220,14 +237,15 @@ TEST(QuadratureGeneratorTestSupport, ExercisesSharedLineRuleChecks) {{-1.0, 0.0, 0.0}, {1.0, 0.0, 0.0}}, {1.0, 1.0}); - EXPECT_DOUBLE_EQ(analytic_line_monomial_integral(0u), 2.0); - EXPECT_DOUBLE_EQ(analytic_line_monomial_integral(1u), 0.0); - EXPECT_DOUBLE_EQ(analytic_line_monomial_integral(2u), 2.0 / 3.0); - EXPECT_DOUBLE_EQ(analytic_line_monomial_integral(255u), 0.0); - EXPECT_NEAR( - accumulate_line_moment(interior_rule, 255u), - 0.0, - kMomentTolerance); + 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( @@ -248,16 +266,15 @@ TEST(QuadratureGeneratorTestSupport, ExercisesSharedLineRuleChecks) "at least one point"); } -TEST(GaussLegendreImplementation, GeneratesRepresentativeSupportedRules) +TEST(GaussLegendreImplementation, GeneratesCanonicalLowOrderRules) { - const std::array point_counts{ - 1, 2, 17, max_gauss_legendre_points()}; - - for (const int num_points : point_counts) { - SCOPED_TRACE(::testing::Message() - << "num_points=" << num_points); - const QuadratureRule rule = - make_gauss_legendre_rule(num_points); + const auto expect_canonical_rule = []( + int num_points, + const auto& expected_points, + const auto& expected_weights) { + SCOPED_TRACE( + ::testing::Message() << "num_points=" << num_points); + const QuadratureRule rule = make_gauss_legendre_rule(num_points); expect_common_line_metadata( rule, @@ -267,21 +284,199 @@ TEST(GaussLegendreImplementation, GeneratesRepresentativeSupportedRules) rule, LineEndpointPolicy::Excluded); expect_advertised_line_exactness(rule); + + ASSERT_EQ(rule.num_points(), expected_points.size()); + ASSERT_EQ(rule.num_points(), expected_weights.size()); + for (std::size_t point_index = 0; + point_index < rule.num_points(); + ++point_index) { + SCOPED_TRACE( + ::testing::Message() << "point index=" << point_index); + if (expected_points[point_index] == 0.0) { + EXPECT_DOUBLE_EQ(rule.point(point_index)[0], 0.0); + } else { + EXPECT_NEAR( + rule.point(point_index)[0], + expected_points[point_index], + kFixtureTolerance); + } + EXPECT_NEAR( + rule.weight(point_index), + expected_weights[point_index], + kFixtureTolerance); + } + + const std::size_t first_unadvertised_even_power = + static_cast(2 * num_points); + 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)); + }; + + expect_canonical_rule( + 1, + std::array{0.0}, + std::array{2.0}); + + const double two_point_abscissa = 1.0 / std::sqrt(3.0); + expect_canonical_rule( + 2, + 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( + 3, + 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) +{ + long double worst_structure_error = 0.0L; + int worst_structure_num_points = 0; + std::size_t worst_structure_point_index = 0u; + std::string_view worst_structure_component = "none"; + + long double worst_measure_error = 0.0L; + int worst_measure_num_points = 0; + + long double worst_moment_error = 0.0L; + int worst_moment_num_points = 0; + int worst_moment_power = 0; + + for (int num_points = 1; + num_points <= max_gauss_legendre_points(); + ++num_points) { + SCOPED_TRACE( + ::testing::Message() << "num_points=" << num_points); + const QuadratureRule rule = make_gauss_legendre_rule(num_points); + + expect_common_line_metadata( + rule, + static_cast(num_points), + 2 * num_points - 1); + expect_line_rule_invariants( + rule, + LineEndpointPolicy::Excluded); + const auto [rule_moment_error, rule_moment_power] = + expect_advertised_line_exactness(rule); + if (worst_moment_num_points == 0 || + rule_moment_error > worst_moment_error) { + worst_moment_error = rule_moment_error; + worst_moment_num_points = num_points; + worst_moment_power = rule_moment_power; + } + + const long double measure_error = std::abs( + accumulate_line_moment(rule, 0u) - + static_cast(rule.reference_cell_measure())); + if (worst_measure_num_points == 0 || + measure_error > worst_measure_error) { + worst_measure_error = measure_error; + worst_measure_num_points = num_points; + } + + const auto update_worst_structure = [ + &worst_structure_error, + &worst_structure_num_points, + &worst_structure_point_index, + &worst_structure_component, + num_points]( + long double error, + std::size_t point_index, + std::string_view component) { + if (worst_structure_num_points == 0 || + error > worst_structure_error) { + worst_structure_error = error; + worst_structure_num_points = num_points; + worst_structure_point_index = point_index; + worst_structure_component = component; + } + }; + + for (std::size_t point_index = 0; + point_index < rule.num_points(); + ++point_index) { + const std::size_t mirror_index = + rule.num_points() - 1u - point_index; + update_worst_structure( + std::abs(static_cast( + rule.point(point_index)[1])), + point_index, + "inactive y coordinate"); + update_worst_structure( + std::abs(static_cast( + rule.point(point_index)[2])), + point_index, + "inactive z coordinate"); + update_worst_structure( + std::abs( + static_cast( + rule.point(point_index)[0]) + + static_cast( + rule.point(mirror_index)[0])), + point_index, + "mirrored point"); + update_worst_structure( + std::abs( + static_cast(rule.weight(point_index)) - + static_cast(rule.weight(mirror_index))), + point_index, + "mirrored weight"); + } + + if (rule.num_points() % 2u == 1u) { + const std::size_t center_index = rule.num_points() / 2u; + update_worst_structure( + std::abs(static_cast( + rule.point(center_index)[0])), + center_index, + "odd-rule center"); + } } + + EXPECT_LE( + worst_structure_error, + static_cast(kStructureTolerance)) + << "worst num_points=" << worst_structure_num_points + << ", point index=" << worst_structure_point_index + << ", component=" << worst_structure_component; + EXPECT_LE( + worst_measure_error, + static_cast(kStructureTolerance)) + << "worst num_points=" << worst_measure_num_points; + EXPECT_LE( + worst_moment_error, + static_cast(kMomentTolerance)) + << "worst num_points=" << worst_moment_num_points + << ", power=" << worst_moment_power; } TEST(GaussLegendreImplementation, RejectsRequestsOutsideSupportedRange) { constexpr std::string_view expected_message = "num_points must be in [1, 128]"; - - expect_exception_with_message( - [] { (void)make_gauss_legendre_rule(0); }, - expected_message); - expect_exception_with_message( - [] { - (void)make_gauss_legendre_rule( - max_gauss_legendre_points() + 1); - }, - expected_message); + constexpr std::array invalid_point_counts{ + std::numeric_limits::min(), + -1, + 0, + 129, + std::numeric_limits::max()}; + + for (const int num_points : invalid_point_counts) { + SCOPED_TRACE( + ::testing::Message() << "num_points=" << num_points); + expect_exception_with_message( + [num_points] { + (void)make_gauss_legendre_rule(num_points); + }, + expected_message); + } } From 772bb3159f778d8f0c53f05a5681e8295d188503 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Wed, 26 Aug 2026 10:10:56 -0700 Subject: [PATCH 07/27] Add Gauss-Lobatto public contract --- .../FE/Quadrature/GaussLobattoQuadrature.h | 60 +++++++++++++++++++ .../Quadrature/test_QuadratureGenerators.cpp | 6 ++ 2 files changed, 66 insertions(+) create mode 100644 Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.h diff --git a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.h b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.h new file mode 100644 index 000000000..a12efddbc --- /dev/null +++ b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.h @@ -0,0 +1,60 @@ +// 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 point count. + * @details The 128 total points include both endpoints. This project support + * bound limits generator work and downstream product-rule growth while + * providing endpoint-inclusive line exactness through degree 253; it is not a + * mathematical or convergence limit. + * @return The inclusive point-count limit, 128. + */ +[[nodiscard]] constexpr int max_gauss_lobatto_points() noexcept +{ + return 128; +} + +/** + * @brief Generate an @p num_points Gauss-Lobatto-Legendre rule on + * @f$[-1,1]@f$. + * + * @details The returned line rule has exactly @f$-1@f$ as its first point and + * exactly @f$+1@f$ as its last point. When present, its @f$n-2@f$ interior + * points are the roots of @f$P'_{n-1}@f$, where @f$n@f$ is @p num_points. + * Points are strictly increasing, weights are positive and aligned with their + * points, and the rule has polynomial exactness @f$2n-3@f$. + * + * @param num_points Signed number of quadrature points; must be in the + * inclusive range + * @f$[2,\texttt{max\_gauss\_lobatto\_points()}]@f$. + * @return A complete QuadratureRule value for CellFamily::Line. + * @throws InvalidArgumentException If @p num_points is outside the supported + * range. + * @throws ConvergenceException If root refinement or final numerical + * validation fails. + */ +[[nodiscard]] QuadratureRule +make_gauss_lobatto_rule(int num_points); + +/** @} */ + +} // namespace svmp::FE::quadrature + +#endif // SVMP_FE_GAUSS_LOBATTO_QUADRATURE_H diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp index 069d1aa2b..4708d9959 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp @@ -9,6 +9,7 @@ #include #include "FE/Common/FEException.h" +#include "FE/Quadrature/GaussLobattoQuadrature.h" #include "FE/Quadrature/GaussQuadrature.h" #include "FE/Quadrature/QuadratureRule.h" @@ -36,6 +37,11 @@ static_assert(noexcept(max_gauss_legendre_points())); static_assert( std::is_same_v); +static_assert(max_gauss_lobatto_points() == 128); +static_assert(noexcept(max_gauss_lobatto_points())); +static_assert( + std::is_same_v); enum class LineEndpointPolicy { Excluded, From 21699dcf69f7e5896e4d08be4c6ff57aa96b7577 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Wed, 26 Aug 2026 11:01:19 -0700 Subject: [PATCH 08/27] Implement bounded Gauss-Lobatto generator --- .../FE/Quadrature/GaussLobattoQuadrature.cpp | 373 ++++++++++++++++++ .../Quadrature/test_QuadratureGenerators.cpp | 38 ++ 2 files changed, 411 insertions(+) create mode 100644 Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp diff --git a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp new file mode 100644 index 000000000..e4899b301 --- /dev/null +++ b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp @@ -0,0 +1,373 @@ +// 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 Bounded generation and validation of 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 + +namespace svmp::FE::quadrature { +namespace { + +constexpr int kMaximumNewtonIterations = 100; +// Accommodate rounding in the Legendre recurrence and Newton update. +constexpr double kNewtonCorrectionTolerance = + 64.0 * std::numeric_limits::epsilon(); +// Allow accumulated rounding across the largest supported rule. +constexpr double kRuleValidationTolerance = + 32.0 * static_cast(max_gauss_lobatto_points()) * + std::numeric_limits::epsilon(); + +std::pair evaluate_adjacent_legendre_values( + int degree, + double coordinate) noexcept +{ + double previous_value = 1.0; + double value = coordinate; + + for (int recurrence_degree = 2; + recurrence_degree <= degree; + ++recurrence_degree) { + const double next_value = + (static_cast(2 * recurrence_degree - 1) * + coordinate * value - + static_cast(recurrence_degree - 1) * previous_value) / + static_cast(recurrence_degree); + previous_value = value; + value = next_value; + } + + return {value, previous_value}; +} + +[[noreturn]] void raise_generation_failure( + int 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; + + const double residual = std::isfinite(diagnostic_value) + ? std::abs(diagnostic_value) + : 0.0; + svmp::raise(message.str(), iteration, residual); +} + +void require_generation( + bool condition, + int num_points, + int half_root_index, + int iteration, + double diagnostic_value, + std::string_view detail) +{ + if (!condition) { + raise_generation_failure( + num_points, + half_root_index, + iteration, + diagnostic_value, + detail); + } +} + +void validate_num_points(int num_points) +{ + if (num_points < 2 || + num_points > max_gauss_lobatto_points()) { + std::ostringstream message; + message << "Gauss-Lobatto-Legendre generator: " + << "num_points must be in [2, " + << max_gauss_lobatto_points() << ']'; + svmp::raise(message.str()); + } +} + +std::pair generate_interior_root_and_weight( + int num_points, + int polynomial_degree, + int half_root_index, + bool is_center, + double normalization) +{ + const double num_points_value = static_cast(num_points); + const double degree_value = static_cast(polynomial_degree); + const double pi = std::numbers::pi_v; + double root = std::cos( + pi * static_cast(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); + require_generation( + std::isfinite(polynomial_value), + num_points, half_root_index, iteration, polynomial_value, + "encountered a non-finite P_m value"); + require_generation( + std::isfinite(previous_polynomial_value), + num_points, half_root_index, iteration, + previous_polynomial_value, + "encountered a non-finite P_(m-1) value"); + + const double residual = + root * polynomial_value - previous_polynomial_value; + require_generation( + std::isfinite(residual), + num_points, half_root_index, iteration, residual, + "computed a non-finite root-function residual"); + + const double derivative = num_points_value * polynomial_value; + require_generation( + std::isfinite(derivative), + num_points, half_root_index, iteration, derivative, + "computed a non-finite root-function derivative"); + require_generation( + derivative != 0.0, + num_points, half_root_index, iteration, derivative, + "encountered a zero root-function derivative"); + + correction = residual / derivative; + require_generation( + std::isfinite(correction), + num_points, half_root_index, iteration, correction, + "computed a non-finite Newton correction"); + + const double updated_root = root - correction; + require_generation( + std::isfinite(updated_root), + num_points, half_root_index, iteration, updated_root, + "computed a non-finite updated root"); + root = updated_root; + + if (std::abs(correction) > kNewtonCorrectionTolerance) { + continue; + } + + if (is_center) { + root = 0.0; + } + + const auto [final_polynomial_value, + final_previous_polynomial_value] = + evaluate_adjacent_legendre_values(polynomial_degree, root); + require_generation( + std::isfinite(final_polynomial_value), + num_points, half_root_index, iteration, + final_polynomial_value, + "refined root produced a non-finite P_m value"); + require_generation( + std::isfinite(final_previous_polynomial_value), + num_points, half_root_index, iteration, + final_previous_polynomial_value, + "refined root produced a non-finite P_(m-1) value"); + + const double final_residual = + root * final_polynomial_value - + final_previous_polynomial_value; + require_generation( + std::isfinite(final_residual), + num_points, half_root_index, iteration, final_residual, + "refined root produced a non-finite residual"); + + const double final_derivative = + num_points_value * final_polynomial_value; + require_generation( + std::isfinite(final_derivative), + num_points, half_root_index, iteration, final_derivative, + "refined root produced a non-finite derivative"); + require_generation( + final_derivative != 0.0, + num_points, half_root_index, iteration, final_derivative, + "refined root produced a zero derivative"); + + const double final_correction = + final_residual / final_derivative; + require_generation( + std::isfinite(final_correction), + num_points, half_root_index, iteration, final_correction, + "refined root produced a non-finite final correction"); + require_generation( + std::abs(final_correction) <= + kNewtonCorrectionTolerance, + num_points, half_root_index, iteration, final_correction, + "refined root failed final correction validation"); + require_generation( + root >= 0.0 && root < 1.0 && + (is_center || root > 0.0), + num_points, half_root_index, iteration, root, + "refined root is outside the expected half interval"); + + const double denominator = + normalization * final_polynomial_value * + final_polynomial_value; + require_generation( + std::isfinite(denominator), + num_points, half_root_index, iteration, denominator, + "refined root produced a non-finite weight denominator"); + require_generation( + denominator > 0.0, + num_points, half_root_index, iteration, denominator, + "refined root produced a non-positive weight denominator"); + + const double weight = 2.0 / denominator; + require_generation( + std::isfinite(weight), + num_points, half_root_index, iteration, weight, + "refined root produced a non-finite quadrature weight"); + require_generation( + weight > 0.0, + num_points, half_root_index, iteration, weight, + "refined root produced a non-positive quadrature weight"); + + return {root, weight}; + } + + raise_generation_failure( + num_points, + half_root_index, + kMaximumNewtonIterations, + correction, + "Newton refinement did not converge"); +} + +void validate_ordering_and_measure( + int num_points, + const std::vector& points, + const std::vector& weights) +{ + 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]; + require_generation( + spacing > 0.0, + num_points, + static_cast(point_index), + -1, + spacing, + "generated points are not strictly increasing"); + } + + const long double weight_sum = + std::accumulate(weights.begin(), weights.end(), 0.0L); + require_generation( + std::isfinite(weight_sum), + num_points, + -1, + -1, + static_cast(weight_sum), + "generated weights produced a non-finite measure"); + + const long double measure_error = std::abs(weight_sum - 2.0L); + require_generation( + measure_error <= + static_cast(kRuleValidationTolerance), + num_points, + -1, + -1, + static_cast(measure_error), + "generated weights do not reproduce the reference measure"); +} + +} // namespace + +QuadratureRule make_gauss_lobatto_rule(int num_points) +{ + validate_num_points(num_points); + + const std::size_t point_count = + static_cast(num_points); + std::vector points( + point_count, QuadPoint::Zero()); + std::vector weights(point_count); + + points.front()[0] = -1.0; + points.back()[0] = 1.0; + + const double num_points_value = static_cast(num_points); + const double normalization = + num_points_value * (num_points_value - 1.0); + require_generation( + std::isfinite(normalization), + num_points, -1, -1, normalization, + "computed a non-finite weight normalization"); + require_generation( + normalization > 0.0, + num_points, -1, -1, normalization, + "computed a non-positive weight normalization"); + + const double endpoint_weight = 2.0 / normalization; + require_generation( + std::isfinite(endpoint_weight), + num_points, -1, -1, endpoint_weight, + "computed a non-finite endpoint weight"); + require_generation( + endpoint_weight > 0.0, + num_points, -1, -1, endpoint_weight, + "computed a non-positive endpoint weight"); + weights.front() = endpoint_weight; + weights.back() = endpoint_weight; + + const int polynomial_degree = num_points - 1; + const int interior_roots_to_refine = (num_points - 1) / 2; + for (int half_root_index = 0; + half_root_index < interior_roots_to_refine; + ++half_root_index) { + const std::size_t left_index = + 1u + static_cast(half_root_index); + const std::size_t right_index = + point_count - 2u - + static_cast(half_root_index); + const auto [root, weight] = + generate_interior_root_and_weight( + num_points, + polynomial_degree, + half_root_index, + left_index == right_index, + normalization); + + points[left_index][0] = -root; + points[right_index][0] = root; + weights[left_index] = weight; + weights[right_index] = weight; + } + + validate_ordering_and_measure(num_points, points, weights); + + 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/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp index 4708d9959..4aa78816d 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp @@ -486,3 +486,41 @@ TEST(GaussLegendreImplementation, RejectsRequestsOutsideSupportedRange) expected_message); } } + +TEST(GaussLobattoImplementation, GeneratesRepresentativeSupportedRules) +{ + const std::array point_counts{ + 2, 3, 17, max_gauss_lobatto_points()}; + + for (const int num_points : point_counts) { + SCOPED_TRACE( + ::testing::Message() << "num_points=" << num_points); + const QuadratureRule rule = + make_gauss_lobatto_rule(num_points); + + expect_common_line_metadata( + rule, + static_cast(num_points), + 2 * num_points - 3); + expect_line_rule_invariants( + rule, + LineEndpointPolicy::Included); + expect_advertised_line_exactness(rule); + } +} + +TEST(GaussLobattoImplementation, RejectsRequestsOutsideSupportedRange) +{ + constexpr std::string_view expected_message = + "num_points must be in [2, 128]"; + + expect_exception_with_message( + [] { (void)make_gauss_lobatto_rule(1); }, + expected_message); + expect_exception_with_message( + [] { + (void)make_gauss_lobatto_rule( + max_gauss_lobatto_points() + 1); + }, + expected_message); +} From 4c7b4d96bc78591f2e12ef25f048d34a72e569df Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Wed, 26 Aug 2026 14:56:38 -0700 Subject: [PATCH 09/27] Complete Gauss-Lobatto generator tests --- .../Quadrature/test_QuadratureGenerators.cpp | 368 +++++++++++------- 1 file changed, 235 insertions(+), 133 deletions(-) diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp index 4aa78816d..aad25f8fc 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp @@ -200,6 +200,148 @@ std::pair expect_advertised_line_exactness( return {worst_error, worst_power}; } +void expect_every_supported_line_rule( + QuadratureRule (*generator)(int), + int first_num_points, + int last_num_points, + int exactness_subtrahend, + LineEndpointPolicy endpoint_policy) +{ + long double worst_structure_error = 0.0L; + int worst_structure_num_points = 0; + std::size_t worst_structure_point_index = 0u; + std::string_view worst_structure_component = "none"; + + long double worst_measure_error = 0.0L; + int worst_measure_num_points = 0; + + long double worst_moment_error = 0.0L; + int worst_moment_num_points = 0; + int worst_moment_power = 0; + + for (int num_points = first_num_points; + num_points <= last_num_points; + ++num_points) { + SCOPED_TRACE( + ::testing::Message() << "num_points=" << num_points); + const QuadratureRule rule = generator(num_points); + + expect_common_line_metadata( + rule, + static_cast(num_points), + 2 * num_points - exactness_subtrahend); + expect_line_rule_invariants(rule, endpoint_policy); + const auto [rule_moment_error, rule_moment_power] = + expect_advertised_line_exactness(rule); + if (worst_moment_num_points == 0 || + rule_moment_error > worst_moment_error) { + worst_moment_error = rule_moment_error; + worst_moment_num_points = num_points; + worst_moment_power = rule_moment_power; + } + + const long double measure_error = std::abs( + accumulate_line_moment(rule, 0u) - + static_cast(rule.reference_cell_measure())); + if (worst_measure_num_points == 0 || + measure_error > worst_measure_error) { + worst_measure_error = measure_error; + worst_measure_num_points = num_points; + } + + const auto update_worst_structure = [ + &worst_structure_error, + &worst_structure_num_points, + &worst_structure_point_index, + &worst_structure_component, + num_points]( + long double error, + std::size_t point_index, + std::string_view component) { + if (worst_structure_num_points == 0 || + error > worst_structure_error) { + worst_structure_error = error; + worst_structure_num_points = num_points; + worst_structure_point_index = point_index; + worst_structure_component = component; + } + }; + + for (std::size_t point_index = 0; + point_index < rule.num_points(); + ++point_index) { + const std::size_t mirror_index = + rule.num_points() - 1u - point_index; + update_worst_structure( + std::abs(static_cast( + rule.point(point_index)[1])), + point_index, + "inactive y coordinate"); + update_worst_structure( + std::abs(static_cast( + rule.point(point_index)[2])), + point_index, + "inactive z coordinate"); + update_worst_structure( + std::abs( + static_cast( + rule.point(point_index)[0]) + + static_cast( + rule.point(mirror_index)[0])), + point_index, + "mirrored point"); + update_worst_structure( + std::abs( + static_cast(rule.weight(point_index)) - + static_cast(rule.weight(mirror_index))), + point_index, + "mirrored weight"); + } + + if (rule.num_points() % 2u == 1u) { + const std::size_t center_index = rule.num_points() / 2u; + update_worst_structure( + std::abs(static_cast( + rule.point(center_index)[0])), + center_index, + "odd-rule center"); + } + + if (endpoint_policy == LineEndpointPolicy::Included) { + update_worst_structure( + std::abs( + static_cast(rule.point(0u)[0]) + 1.0L), + 0u, + "left endpoint"); + const std::size_t right_endpoint_index = + rule.num_points() - 1u; + update_worst_structure( + std::abs( + static_cast( + rule.point(right_endpoint_index)[0]) - + 1.0L), + right_endpoint_index, + "right endpoint"); + } + } + + EXPECT_LE( + worst_structure_error, + static_cast(kStructureTolerance)) + << "worst num_points=" << worst_structure_num_points + << ", point index=" << worst_structure_point_index + << ", component=" << worst_structure_component; + EXPECT_LE( + worst_measure_error, + static_cast(kStructureTolerance)) + << "worst num_points=" << worst_measure_num_points; + EXPECT_LE( + worst_moment_error, + static_cast(kMomentTolerance)) + << "worst num_points=" << worst_moment_num_points + << ", power=" << worst_moment_power; +} + template void expect_exception_with_message( Function&& function, @@ -345,124 +487,12 @@ TEST(GaussLegendreImplementation, GeneratesCanonicalLowOrderRules) TEST(GaussLegendreImplementation, GeneratesEverySupportedRule) { - long double worst_structure_error = 0.0L; - int worst_structure_num_points = 0; - std::size_t worst_structure_point_index = 0u; - std::string_view worst_structure_component = "none"; - - long double worst_measure_error = 0.0L; - int worst_measure_num_points = 0; - - long double worst_moment_error = 0.0L; - int worst_moment_num_points = 0; - int worst_moment_power = 0; - - for (int num_points = 1; - num_points <= max_gauss_legendre_points(); - ++num_points) { - SCOPED_TRACE( - ::testing::Message() << "num_points=" << num_points); - const QuadratureRule rule = make_gauss_legendre_rule(num_points); - - expect_common_line_metadata( - rule, - static_cast(num_points), - 2 * num_points - 1); - expect_line_rule_invariants( - rule, - LineEndpointPolicy::Excluded); - const auto [rule_moment_error, rule_moment_power] = - expect_advertised_line_exactness(rule); - if (worst_moment_num_points == 0 || - rule_moment_error > worst_moment_error) { - worst_moment_error = rule_moment_error; - worst_moment_num_points = num_points; - worst_moment_power = rule_moment_power; - } - - const long double measure_error = std::abs( - accumulate_line_moment(rule, 0u) - - static_cast(rule.reference_cell_measure())); - if (worst_measure_num_points == 0 || - measure_error > worst_measure_error) { - worst_measure_error = measure_error; - worst_measure_num_points = num_points; - } - - const auto update_worst_structure = [ - &worst_structure_error, - &worst_structure_num_points, - &worst_structure_point_index, - &worst_structure_component, - num_points]( - long double error, - std::size_t point_index, - std::string_view component) { - if (worst_structure_num_points == 0 || - error > worst_structure_error) { - worst_structure_error = error; - worst_structure_num_points = num_points; - worst_structure_point_index = point_index; - worst_structure_component = component; - } - }; - - for (std::size_t point_index = 0; - point_index < rule.num_points(); - ++point_index) { - const std::size_t mirror_index = - rule.num_points() - 1u - point_index; - update_worst_structure( - std::abs(static_cast( - rule.point(point_index)[1])), - point_index, - "inactive y coordinate"); - update_worst_structure( - std::abs(static_cast( - rule.point(point_index)[2])), - point_index, - "inactive z coordinate"); - update_worst_structure( - std::abs( - static_cast( - rule.point(point_index)[0]) + - static_cast( - rule.point(mirror_index)[0])), - point_index, - "mirrored point"); - update_worst_structure( - std::abs( - static_cast(rule.weight(point_index)) - - static_cast(rule.weight(mirror_index))), - point_index, - "mirrored weight"); - } - - if (rule.num_points() % 2u == 1u) { - const std::size_t center_index = rule.num_points() / 2u; - update_worst_structure( - std::abs(static_cast( - rule.point(center_index)[0])), - center_index, - "odd-rule center"); - } - } - - EXPECT_LE( - worst_structure_error, - static_cast(kStructureTolerance)) - << "worst num_points=" << worst_structure_num_points - << ", point index=" << worst_structure_point_index - << ", component=" << worst_structure_component; - EXPECT_LE( - worst_measure_error, - static_cast(kStructureTolerance)) - << "worst num_points=" << worst_measure_num_points; - EXPECT_LE( - worst_moment_error, - static_cast(kMomentTolerance)) - << "worst num_points=" << worst_moment_num_points - << ", power=" << worst_moment_power; + expect_every_supported_line_rule( + &make_gauss_legendre_rule, + 1, + max_gauss_legendre_points(), + 1, + LineEndpointPolicy::Excluded); } TEST(GaussLegendreImplementation, RejectsRequestsOutsideSupportedRange) @@ -487,12 +517,12 @@ TEST(GaussLegendreImplementation, RejectsRequestsOutsideSupportedRange) } } -TEST(GaussLobattoImplementation, GeneratesRepresentativeSupportedRules) +TEST(GaussLobattoImplementation, GeneratesCanonicalLowOrderRules) { - const std::array point_counts{ - 2, 3, 17, max_gauss_lobatto_points()}; - - for (const int num_points : point_counts) { + const auto expect_canonical_rule = []( + int num_points, + const auto& expected_points, + const auto& expected_weights) { SCOPED_TRACE( ::testing::Message() << "num_points=" << num_points); const QuadratureRule rule = @@ -506,21 +536,93 @@ TEST(GaussLobattoImplementation, GeneratesRepresentativeSupportedRules) rule, LineEndpointPolicy::Included); expect_advertised_line_exactness(rule); - } + + ASSERT_EQ(rule.num_points(), expected_points.size()); + ASSERT_EQ(rule.num_points(), expected_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 expected_coordinate = + expected_points[point_index]; + if (expected_coordinate == -1.0 || + expected_coordinate == 0.0 || + expected_coordinate == 1.0) { + EXPECT_DOUBLE_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(2 * num_points - 2); + 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)); + }; + + expect_canonical_rule( + 2, + std::array{-1.0, 1.0}, + std::array{1.0, 1.0}); + expect_canonical_rule( + 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( + 4, + 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) +{ + expect_every_supported_line_rule( + &make_gauss_lobatto_rule, + 2, + max_gauss_lobatto_points(), + 3, + LineEndpointPolicy::Included); } TEST(GaussLobattoImplementation, RejectsRequestsOutsideSupportedRange) { constexpr std::string_view expected_message = + "Gauss-Lobatto-Legendre generator: " "num_points must be in [2, 128]"; + constexpr std::array invalid_point_counts{ + std::numeric_limits::min(), + -1, + 0, + 1, + 129, + std::numeric_limits::max()}; - expect_exception_with_message( - [] { (void)make_gauss_lobatto_rule(1); }, - expected_message); - expect_exception_with_message( - [] { - (void)make_gauss_lobatto_rule( - max_gauss_lobatto_points() + 1); - }, - expected_message); + for (const int num_points : invalid_point_counts) { + SCOPED_TRACE( + ::testing::Message() << "num_points=" << num_points); + expect_exception_with_message( + [num_points] { + (void)make_gauss_lobatto_rule(num_points); + }, + expected_message); + } } From c20d1fe90d174af32a85e5cb0fe4b05d46224e62 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Wed, 26 Aug 2026 17:36:17 -0700 Subject: [PATCH 10/27] Add Lobatto Basis consistency coverage --- .../Quadrature/test_QuadratureGenerators.cpp | 54 +++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp index aad25f8fc..9285a26bc 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp @@ -8,6 +8,7 @@ #include +#include "FE/Basis/NodeOrderingConventions.h" #include "FE/Common/FEException.h" #include "FE/Quadrature/GaussLobattoQuadrature.h" #include "FE/Quadrature/GaussQuadrature.h" @@ -31,6 +32,10 @@ constexpr double kStructureTolerance = 1.0e-12; constexpr double kMomentTolerance = 2.0e-12; 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_points() == 128); static_assert(noexcept(max_gauss_legendre_points())); @@ -626,3 +631,52 @@ TEST(GaussLobattoImplementation, RejectsRequestsOutsideSupportedRange) expected_message); } } + +TEST(GaussLobattoBasisConsistency, MatchesRepresentativeNodeDistributions) +{ + constexpr std::array point_counts{ + 2, 4, 65, max_gauss_lobatto_points()}; + + for (const int num_points : point_counts) { + SCOPED_TRACE( + ::testing::Message() << "num_points=" << num_points); + const QuadratureRule rule = + make_gauss_lobatto_rule(num_points); + ASSERT_EQ( + rule.num_points(), + static_cast(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), + num_points - 1); + + if (point_index == 0u) { + EXPECT_EQ(quadrature_coordinate, -1.0); + EXPECT_EQ(basis_coordinate, -1.0); + EXPECT_EQ(quadrature_coordinate, basis_coordinate); + } else if (point_index + 1u == rule.num_points()) { + EXPECT_EQ(quadrature_coordinate, 1.0); + EXPECT_EQ(basis_coordinate, 1.0); + EXPECT_EQ(quadrature_coordinate, basis_coordinate); + } else if (num_points % 2 == 1 && + point_index == rule.num_points() / 2u) { + EXPECT_EQ(quadrature_coordinate, 0.0); + EXPECT_EQ(basis_coordinate, 0.0); + EXPECT_EQ(quadrature_coordinate, basis_coordinate); + } else { + EXPECT_NEAR( + quadrature_coordinate, + basis_coordinate, + kBasisConsistencyTolerance); + } + } + } +} From 20d28691c8e24275061584c8857d319744271325 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Wed, 26 Aug 2026 22:43:13 -0700 Subject: [PATCH 11/27] Harden quadrature generator validation --- .../FE/Quadrature/GaussLobattoQuadrature.cpp | 126 ++++++++---------- .../solver/FE/Quadrature/GaussQuadrature.cpp | 78 +++++------ .../Quadrature/test_QuadratureGenerators.cpp | 4 + 3 files changed, 98 insertions(+), 110 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp index e4899b301..32b6845c0 100644 --- a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp +++ b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp @@ -24,11 +24,15 @@ namespace svmp::FE::quadrature { namespace { +// Defensively bound supported cosine-seeded Newton refinements for +// deterministic termination. constexpr int kMaximumNewtonIterations = 100; -// Accommodate rounding in the Legendre recurrence and Newton update. +// Guard recurrence and Newton-update rounding; exhaustive supported-size +// sweeps qualify this scale. constexpr double kNewtonCorrectionTolerance = 64.0 * std::numeric_limits::epsilon(); -// Allow accumulated rounding across the largest supported rule. +// Provide conservative O(n epsilon) accumulation headroom, qualified by those +// sweeps. constexpr double kRuleValidationTolerance = 32.0 * static_cast(max_gauss_lobatto_points()) * std::numeric_limits::epsilon(); @@ -92,24 +96,12 @@ void require_generation( } } -void validate_num_points(int num_points) -{ - if (num_points < 2 || - num_points > max_gauss_lobatto_points()) { - std::ostringstream message; - message << "Gauss-Lobatto-Legendre generator: " - << "num_points must be in [2, " - << max_gauss_lobatto_points() << ']'; - svmp::raise(message.str()); - } -} - std::pair generate_interior_root_and_weight( int num_points, int polynomial_degree, int half_root_index, bool is_center, - double normalization) + double weight_denominator_scale) { const double num_points_value = static_cast(num_points); const double degree_value = static_cast(polynomial_degree); @@ -225,7 +217,7 @@ std::pair generate_interior_root_and_weight( "refined root is outside the expected half interval"); const double denominator = - normalization * final_polynomial_value * + weight_denominator_scale * final_polynomial_value * final_polynomial_value; require_generation( std::isfinite(denominator), @@ -257,51 +249,18 @@ std::pair generate_interior_root_and_weight( "Newton refinement did not converge"); } -void validate_ordering_and_measure( - int num_points, - const std::vector& points, - const std::vector& weights) -{ - 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]; - require_generation( - spacing > 0.0, - num_points, - static_cast(point_index), - -1, - spacing, - "generated points are not strictly increasing"); - } - - const long double weight_sum = - std::accumulate(weights.begin(), weights.end(), 0.0L); - require_generation( - std::isfinite(weight_sum), - num_points, - -1, - -1, - static_cast(weight_sum), - "generated weights produced a non-finite measure"); - - const long double measure_error = std::abs(weight_sum - 2.0L); - require_generation( - measure_error <= - static_cast(kRuleValidationTolerance), - num_points, - -1, - -1, - static_cast(measure_error), - "generated weights do not reproduce the reference measure"); -} - } // namespace QuadratureRule make_gauss_lobatto_rule(int num_points) { - validate_num_points(num_points); + if (num_points < 2 || + num_points > max_gauss_lobatto_points()) { + std::ostringstream message; + message << "Gauss-Lobatto-Legendre generator: " + << "num_points must be in [2, " + << max_gauss_lobatto_points() << ']'; + svmp::raise(message.str()); + } const std::size_t point_count = static_cast(num_points); @@ -313,18 +272,18 @@ QuadratureRule make_gauss_lobatto_rule(int num_points) points.back()[0] = 1.0; const double num_points_value = static_cast(num_points); - const double normalization = + const double weight_denominator_scale = num_points_value * (num_points_value - 1.0); require_generation( - std::isfinite(normalization), - num_points, -1, -1, normalization, - "computed a non-finite weight normalization"); + std::isfinite(weight_denominator_scale), + num_points, -1, -1, weight_denominator_scale, + "computed a non-finite weight denominator scale"); require_generation( - normalization > 0.0, - num_points, -1, -1, normalization, - "computed a non-positive weight normalization"); + weight_denominator_scale > 0.0, + num_points, -1, -1, weight_denominator_scale, + "computed a non-positive weight denominator scale"); - const double endpoint_weight = 2.0 / normalization; + const double endpoint_weight = 2.0 / weight_denominator_scale; require_generation( std::isfinite(endpoint_weight), num_points, -1, -1, endpoint_weight, @@ -352,7 +311,7 @@ QuadratureRule make_gauss_lobatto_rule(int num_points) polynomial_degree, half_root_index, left_index == right_index, - normalization); + weight_denominator_scale); points[left_index][0] = -root; points[right_index][0] = root; @@ -360,7 +319,40 @@ QuadratureRule make_gauss_lobatto_rule(int num_points) weights[right_index] = weight; } - validate_ordering_and_measure(num_points, points, weights); + 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]; + require_generation( + spacing > 0.0, + 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); + require_generation( + std::isfinite(weight_sum), + num_points, + -1, + -1, + static_cast(weight_sum), + "generated weights produced a non-finite measure"); + + const long double measure_error = std::abs(weight_sum - 2.0L); + require_generation( + measure_error <= + static_cast(kRuleValidationTolerance), + num_points, + -1, + -1, + static_cast(measure_error), + "generated weights do not reproduce the reference measure"); const int polynomial_exactness = 2 * num_points - 3; return QuadratureRule( diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp index 766228570..a5b52e37d 100644 --- a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp @@ -24,10 +24,14 @@ namespace svmp::FE::quadrature { namespace { +// Defensively bound supported cosine-seeded Newton refinements for +// deterministic termination. constexpr int kMaximumNewtonIterations = 100; -// Accommodate rounding in the Legendre recurrence and Newton update. +// Guard recurrence and Newton-update rounding; exhaustive supported-size +// sweeps qualify this scale. constexpr double kNewtonCorrectionTolerance = 64.0 * std::numeric_limits::epsilon(); -// Allow accumulated rounding across the largest supported rule. +// Provide conservative O(n epsilon) accumulation headroom, qualified by those +// sweeps. constexpr double kRuleValidationTolerance = 32.0 * static_cast(max_gauss_legendre_points()) * std::numeric_limits::epsilon(); @@ -101,17 +105,6 @@ void require_generation( } } -void validate_num_points(int num_points) -{ - if (num_points < 1 || - num_points > max_gauss_legendre_points()) { - std::ostringstream message; - message << "Gauss-Legendre generator: num_points must be in [1, " - << max_gauss_legendre_points() << ']'; - svmp::raise(message.str()); - } -} - std::pair generate_root_and_weight( int num_points, int root_index, @@ -197,39 +190,17 @@ std::pair generate_root_and_weight( "Newton refinement did not converge"); } -void validate_ordering_and_measure( - int num_points, - const std::vector& points, - const std::vector& weights) -{ - 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]; - require_generation( - spacing > 0.0, - num_points, static_cast(point_index), -1, spacing, - "generated points are not strictly increasing"); - } - - const long double weight_sum = - std::accumulate(weights.begin(), weights.end(), 0.0L); - const long double measure_error = - std::abs(weight_sum - 2.0L); - require_generation( - std::isfinite(weight_sum) && - measure_error <= - static_cast(kRuleValidationTolerance), - num_points, -1, -1, static_cast(measure_error), - "generated weights do not reproduce the reference measure"); -} - } // namespace QuadratureRule make_gauss_legendre_rule(int num_points) { - validate_num_points(num_points); + if (num_points < 1 || + num_points > max_gauss_legendre_points()) { + std::ostringstream message; + message << "Gauss-Legendre generator: num_points must be in [1, " + << max_gauss_legendre_points() << ']'; + svmp::raise(message.str()); + } const std::size_t point_count = static_cast(num_points); @@ -256,7 +227,28 @@ QuadratureRule make_gauss_legendre_rule(int num_points) weights[right_index] = weight; } - validate_ordering_and_measure(num_points, points, weights); + 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]; + require_generation( + spacing > 0.0, + 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); + require_generation( + std::isfinite(weight_sum) && + measure_error <= + static_cast(kRuleValidationTolerance), + num_points, -1, -1, static_cast(measure_error), + "generated weights do not reproduce the reference measure"); const int polynomial_exactness = 2 * num_points - 1; return QuadratureRule( diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp index 9285a26bc..c1b99300a 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp @@ -28,8 +28,12 @@ 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 From faa1db054dfb47ed79eac3f972ebca4c44f1f64d Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 27 Aug 2026 01:47:29 -0700 Subject: [PATCH 12/27] Verify Phase 02 quadrature generators From f8a0bd5482e7867e7948c1679df62d14675ace96 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Mon, 31 Aug 2026 12:52:47 -0700 Subject: [PATCH 13/27] Streamline Gauss quadrature generators --- .../FE/Quadrature/GaussLobattoQuadrature.cpp | 224 +++++------------- .../solver/FE/Quadrature/GaussQuadrature.cpp | 89 +++---- 2 files changed, 85 insertions(+), 228 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp index 32b6845c0..45138fc51 100644 --- a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp +++ b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp @@ -44,9 +44,7 @@ std::pair evaluate_adjacent_legendre_values( double previous_value = 1.0; double value = coordinate; - for (int recurrence_degree = 2; - recurrence_degree <= degree; - ++recurrence_degree) { + for (int recurrence_degree = 2; recurrence_degree <= degree; ++recurrence_degree) { const double next_value = (static_cast(2 * recurrence_degree - 1) * coordinate * value - @@ -59,12 +57,8 @@ std::pair evaluate_adjacent_legendre_values( return {value, previous_value}; } -[[noreturn]] void raise_generation_failure( - int num_points, - int half_root_index, - int iteration, - double diagnostic_value, - std::string_view detail) +[[noreturn]] void raise_generation_failure(int num_points, int half_root_index, + int iteration, double diagnostic_value, std::string_view detail) { std::ostringstream message; message << "Gauss-Lobatto-Legendre generator: " << detail @@ -78,84 +72,53 @@ std::pair evaluate_adjacent_legendre_values( svmp::raise(message.str(), iteration, residual); } -void require_generation( - bool condition, - int num_points, - int half_root_index, - int iteration, - double diagnostic_value, - std::string_view detail) +void require_generation(bool condition, int num_points, int half_root_index, + int iteration, double diagnostic_value, std::string_view detail) { if (!condition) { raise_generation_failure( - num_points, - half_root_index, - iteration, - diagnostic_value, - detail); + num_points, half_root_index, iteration, diagnostic_value, detail); } } -std::pair generate_interior_root_and_weight( - int num_points, - int polynomial_degree, - int half_root_index, - bool is_center, - double weight_denominator_scale) +std::pair generate_interior_root_and_weight(int num_points, int half_root_index, + bool is_center, double weight_denominator_scale) { + const int polynomial_degree = num_points - 1; const double num_points_value = static_cast(num_points); const double degree_value = static_cast(polynomial_degree); const double pi = std::numbers::pi_v; - double root = std::cos( - pi * static_cast(half_root_index + 1) / - degree_value); + double root = std::cos(pi * static_cast(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) { + for (int iteration = 1; iteration <= kMaximumNewtonIterations; ++iteration) { const auto [polynomial_value, previous_polynomial_value] = evaluate_adjacent_legendre_values(polynomial_degree, root); require_generation( - std::isfinite(polynomial_value), - num_points, half_root_index, iteration, polynomial_value, - "encountered a non-finite P_m value"); - require_generation( - std::isfinite(previous_polynomial_value), + std::isfinite(polynomial_value) && + std::isfinite(previous_polynomial_value), num_points, half_root_index, iteration, - previous_polynomial_value, - "encountered a non-finite P_(m-1) value"); + polynomial_value, + "encountered invalid adjacent Legendre values"); const double residual = root * polynomial_value - previous_polynomial_value; - require_generation( - std::isfinite(residual), - num_points, half_root_index, iteration, residual, - "computed a non-finite root-function residual"); - const double derivative = num_points_value * polynomial_value; require_generation( - std::isfinite(derivative), - num_points, half_root_index, iteration, derivative, - "computed a non-finite root-function derivative"); - require_generation( - derivative != 0.0, + std::isfinite(residual) && + std::isfinite(derivative) && + derivative != 0.0, num_points, half_root_index, iteration, derivative, - "encountered a zero root-function derivative"); + "computed an invalid root-function residual or derivative"); correction = residual / derivative; - require_generation( - std::isfinite(correction), - num_points, half_root_index, iteration, correction, - "computed a non-finite Newton correction"); - const double updated_root = root - correction; require_generation( - std::isfinite(updated_root), - num_points, half_root_index, iteration, updated_root, - "computed a non-finite updated root"); + std::isfinite(correction) && std::isfinite(updated_root), + num_points, half_root_index, iteration, correction, + "computed an invalid Newton update"); root = updated_root; if (std::abs(correction) > kNewtonCorrectionTolerance) { @@ -165,87 +128,61 @@ std::pair generate_interior_root_and_weight( if (is_center) { root = 0.0; } + require_generation( + root >= 0.0 && root < 1.0 && (is_center || root > 0.0), + 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); require_generation( - std::isfinite(final_polynomial_value), + std::isfinite(final_polynomial_value) && + std::isfinite(final_previous_polynomial_value), num_points, half_root_index, iteration, final_polynomial_value, - "refined root produced a non-finite P_m value"); - require_generation( - std::isfinite(final_previous_polynomial_value), - num_points, half_root_index, iteration, - final_previous_polynomial_value, - "refined root produced a non-finite P_(m-1) value"); + "refined root produced invalid adjacent Legendre values"); const double final_residual = root * final_polynomial_value - final_previous_polynomial_value; - require_generation( - std::isfinite(final_residual), - num_points, half_root_index, iteration, final_residual, - "refined root produced a non-finite residual"); - const double final_derivative = num_points_value * final_polynomial_value; require_generation( - std::isfinite(final_derivative), - num_points, half_root_index, iteration, final_derivative, - "refined root produced a non-finite derivative"); - require_generation( - final_derivative != 0.0, + std::isfinite(final_residual) && + std::isfinite(final_derivative) && + final_derivative != 0.0, num_points, half_root_index, iteration, final_derivative, - "refined root produced a zero derivative"); + "refined root produced an invalid residual or derivative"); const double final_correction = final_residual / final_derivative; require_generation( - std::isfinite(final_correction), - num_points, half_root_index, iteration, final_correction, - "refined root produced a non-finite final correction"); - require_generation( - std::abs(final_correction) <= - kNewtonCorrectionTolerance, + std::isfinite(final_correction) && + std::abs(final_correction) <= + kNewtonCorrectionTolerance, num_points, half_root_index, iteration, final_correction, "refined root failed final correction validation"); - require_generation( - root >= 0.0 && root < 1.0 && - (is_center || root > 0.0), - num_points, half_root_index, iteration, root, - "refined root is outside the expected half interval"); const double denominator = weight_denominator_scale * final_polynomial_value * final_polynomial_value; require_generation( - std::isfinite(denominator), + std::isfinite(denominator) && denominator > 0.0, num_points, half_root_index, iteration, denominator, - "refined root produced a non-finite weight denominator"); - require_generation( - denominator > 0.0, - num_points, half_root_index, iteration, denominator, - "refined root produced a non-positive weight denominator"); + "refined root produced an invalid weight denominator"); const double weight = 2.0 / denominator; require_generation( - std::isfinite(weight), - num_points, half_root_index, iteration, weight, - "refined root produced a non-finite quadrature weight"); - require_generation( - weight > 0.0, + std::isfinite(weight) && weight > 0.0, num_points, half_root_index, iteration, weight, - "refined root produced a non-positive quadrature weight"); + "refined root produced an invalid quadrature weight"); return {root, weight}; } raise_generation_failure( - num_points, - half_root_index, - kMaximumNewtonIterations, - correction, + num_points, half_root_index, kMaximumNewtonIterations, correction, "Newton refinement did not converge"); } @@ -253,8 +190,7 @@ std::pair generate_interior_root_and_weight( QuadratureRule make_gauss_lobatto_rule(int num_points) { - if (num_points < 2 || - num_points > max_gauss_lobatto_points()) { + if (num_points < 2 || num_points > max_gauss_lobatto_points()) { std::ostringstream message; message << "Gauss-Lobatto-Legendre generator: " << "num_points must be in [2, " @@ -262,55 +198,30 @@ QuadratureRule make_gauss_lobatto_rule(int num_points) svmp::raise(message.str()); } - const std::size_t point_count = - static_cast(num_points); std::vector points( - point_count, QuadPoint::Zero()); - std::vector weights(point_count); + static_cast(num_points), QuadPoint::Zero()); + std::vector weights(points.size()); points.front()[0] = -1.0; points.back()[0] = 1.0; - const double num_points_value = static_cast(num_points); const double weight_denominator_scale = - num_points_value * (num_points_value - 1.0); - require_generation( - std::isfinite(weight_denominator_scale), - num_points, -1, -1, weight_denominator_scale, - "computed a non-finite weight denominator scale"); - require_generation( - weight_denominator_scale > 0.0, - num_points, -1, -1, weight_denominator_scale, - "computed a non-positive weight denominator scale"); - + static_cast(num_points * (num_points - 1)); const double endpoint_weight = 2.0 / weight_denominator_scale; - require_generation( - std::isfinite(endpoint_weight), - num_points, -1, -1, endpoint_weight, - "computed a non-finite endpoint weight"); - require_generation( - endpoint_weight > 0.0, - num_points, -1, -1, endpoint_weight, - "computed a non-positive endpoint weight"); weights.front() = endpoint_weight; weights.back() = endpoint_weight; - const int polynomial_degree = num_points - 1; const int interior_roots_to_refine = (num_points - 1) / 2; for (int half_root_index = 0; - half_root_index < interior_roots_to_refine; - ++half_root_index) { + half_root_index < interior_roots_to_refine; ++half_root_index) { const std::size_t left_index = 1u + static_cast(half_root_index); const std::size_t right_index = - point_count - 2u - + points.size() - 2u - static_cast(half_root_index); const auto [root, weight] = generate_interior_root_and_weight( - num_points, - polynomial_degree, - half_root_index, - left_index == right_index, + num_points, half_root_index, left_index == right_index, weight_denominator_scale); points[left_index][0] = -root; @@ -319,47 +230,28 @@ QuadratureRule make_gauss_lobatto_rule(int num_points) 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]; + 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]; require_generation( spacing > 0.0, - num_points, - static_cast(point_index), - -1, - spacing, + 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); - require_generation( - std::isfinite(weight_sum), - num_points, - -1, - -1, - static_cast(weight_sum), - "generated weights produced a non-finite measure"); - + const long double weight_sum = std::accumulate(weights.begin(), weights.end(), 0.0L); const long double measure_error = std::abs(weight_sum - 2.0L); require_generation( - measure_error <= - static_cast(kRuleValidationTolerance), - num_points, - -1, - -1, - static_cast(measure_error), + std::isfinite(weight_sum) && + measure_error <= + static_cast(kRuleValidationTolerance), + num_points, -1, -1, static_cast(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)); + svmp::CellFamily::Line, polynomial_exactness, + std::move(points), std::move(weights)); } } // namespace svmp::FE::quadrature diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp index a5b52e37d..49ba8c748 100644 --- a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp @@ -29,11 +29,13 @@ namespace { 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(); +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 * static_cast(max_gauss_legendre_points()) * - std::numeric_limits::epsilon(); +constexpr double kRuleValidationTolerance = + 32.0 * static_cast(max_gauss_legendre_points()) * + std::numeric_limits::epsilon(); std::pair evaluate_legendre_with_derivative( int degree, @@ -41,27 +43,17 @@ std::pair evaluate_legendre_with_derivative( { double previous_value = 1.0; double previous_derivative = 0.0; - - if (degree == 0) { - return {previous_value, previous_derivative}; - } - double value = coordinate; double derivative = 1.0; for (int recurrence_degree = 2; recurrence_degree <= degree; ++recurrence_degree) { - const double degree_value = - static_cast(recurrence_degree); - const double recurrence_factor = - static_cast(2 * recurrence_degree - 1); + const double degree_value = static_cast(recurrence_degree); + const double recurrence_factor = static_cast(2 * recurrence_degree - 1); const double next_value = (recurrence_factor * coordinate * value - - static_cast(recurrence_degree - 1) * previous_value) / - degree_value; + static_cast(recurrence_degree - 1) * previous_value) / degree_value; const double next_derivative = (recurrence_factor * (value + coordinate * derivative) - - static_cast(recurrence_degree - 1) * - previous_derivative) / - degree_value; + static_cast(recurrence_degree - 1) * previous_derivative) / degree_value; previous_value = value; previous_derivative = derivative; @@ -91,24 +83,15 @@ std::pair evaluate_legendre_with_derivative( svmp::raise(message.str(), iteration, residual); } -void require_generation( - bool condition, - int num_points, - int root_index, - int iteration, - double diagnostic_value, - std::string_view detail) +void require_generation(bool condition, int num_points, int root_index, int iteration, + double diagnostic_value, std::string_view detail) { if (!condition) { - raise_generation_failure( - num_points, root_index, iteration, diagnostic_value, detail); + raise_generation_failure(num_points, root_index, iteration, diagnostic_value, detail); } } -std::pair generate_root_and_weight( - int num_points, - int root_index, - bool is_center) +std::pair generate_root_and_weight(int num_points, int root_index, bool is_center) { const double pi = std::numbers::pi_v; double root = std::cos( @@ -116,9 +99,7 @@ std::pair generate_root_and_weight( (static_cast(num_points) + 0.5)); double correction = 0.0; - for (int iteration = 1; - iteration <= kMaximumNewtonIterations; - ++iteration) { + for (int iteration = 1; iteration <= kMaximumNewtonIterations; ++iteration) { const auto [polynomial_value, polynomial_derivative] = evaluate_legendre_with_derivative(num_points, root); require_generation( @@ -194,32 +175,23 @@ std::pair generate_root_and_weight( QuadratureRule make_gauss_legendre_rule(int num_points) { - if (num_points < 1 || - num_points > max_gauss_legendre_points()) { + if (num_points < 1 || num_points > max_gauss_legendre_points()) { std::ostringstream message; message << "Gauss-Legendre generator: num_points must be in [1, " << max_gauss_legendre_points() << ']'; svmp::raise(message.str()); } - const std::size_t point_count = - static_cast(num_points); std::vector points( - point_count, QuadPoint::Zero()); - std::vector weights(point_count); + static_cast(num_points), QuadPoint::Zero()); + std::vector weights(points.size()); const int roots_to_refine = (num_points + 1) / 2; - for (int root_index = 0; - root_index < roots_to_refine; - ++root_index) { - const std::size_t left_index = - static_cast(root_index); - const std::size_t right_index = - point_count - 1u - left_index; + for (int root_index = 0; root_index < roots_to_refine; ++root_index) { + const std::size_t left_index = static_cast(root_index); + const std::size_t right_index = points.size() - 1u - left_index; const auto [root, weight] = generate_root_and_weight( - num_points, - root_index, - left_index == right_index); + num_points, root_index, left_index == right_index); points[left_index][0] = -root; points[right_index][0] = root; @@ -227,11 +199,8 @@ QuadratureRule make_gauss_legendre_rule(int num_points) 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]; + 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]; require_generation( spacing > 0.0, num_points, static_cast(point_index), -1, spacing, @@ -239,10 +208,8 @@ QuadratureRule make_gauss_legendre_rule(int num_points) } // 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); + const long double weight_sum = std::accumulate(weights.begin(), weights.end(), 0.0L); + const long double measure_error = std::abs(weight_sum - 2.0L); require_generation( std::isfinite(weight_sum) && measure_error <= @@ -252,10 +219,8 @@ QuadratureRule make_gauss_legendre_rule(int num_points) const int polynomial_exactness = 2 * num_points - 1; return QuadratureRule( - svmp::CellFamily::Line, - polynomial_exactness, - std::move(points), - std::move(weights)); + svmp::CellFamily::Line, polynomial_exactness, + std::move(points), std::move(weights)); } } // namespace svmp::FE::quadrature From 519f7749f86a1d3ef6cb9c22dbb326b90cc170b6 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Sat, 12 Sep 2026 01:12:37 -0700 Subject: [PATCH 14/27] Simplify quadrature generator tests --- .../Quadrature/test_QuadratureGenerators.cpp | 296 +++--------------- 1 file changed, 46 insertions(+), 250 deletions(-) diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp index c1b99300a..149407f11 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp @@ -3,7 +3,7 @@ /** * @file test_QuadratureGenerators.cpp - * @brief Shared test support for bounded one-dimensional quadrature generators. + * @brief Tests for bounded one-dimensional quadrature generators. */ #include @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -100,7 +101,7 @@ void expect_common_line_metadata( { EXPECT_EQ(rule.cell_family(), svmp::CellFamily::Line); EXPECT_EQ(rule.dimension(), 1u); - EXPECT_DOUBLE_EQ(rule.reference_cell_measure(), 2.0); + 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); @@ -110,8 +111,8 @@ void expect_common_line_metadata( point_index < rule.num_points(); ++point_index) { SCOPED_TRACE(::testing::Message() << "point index=" << point_index); - EXPECT_DOUBLE_EQ(rule.point(point_index)[1], 0.0); - EXPECT_DOUBLE_EQ(rule.point(point_index)[2], 0.0); + EXPECT_EQ(rule.point(point_index)[1], 0.0); + EXPECT_EQ(rule.point(point_index)[2], 0.0); } } @@ -164,12 +165,12 @@ void expect_line_rule_invariants( if (endpoint_policy == LineEndpointPolicy::Included) { ASSERT_GE(rule.num_points(), 2u); - EXPECT_DOUBLE_EQ(rule.point(0)[0], -1.0); - EXPECT_DOUBLE_EQ(rule.point(rule.num_points() - 1u)[0], 1.0); + 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_DOUBLE_EQ(rule.point(rule.num_points() / 2u)[0], 0.0); + EXPECT_EQ(rule.point(rule.num_points() / 2u)[0], 0.0); } const long double measure_error = std::abs( @@ -178,17 +179,11 @@ void expect_line_rule_invariants( EXPECT_LE(measure_error, static_cast(tolerance)); } -std::pair expect_advertised_line_exactness( +void expect_advertised_line_exactness( const QuadratureRule& rule, double tolerance = kMomentTolerance) { - EXPECT_GE(rule.polynomial_exactness(), 0); - if (rule.polynomial_exactness() < 0) { - return {std::numeric_limits::infinity(), -1}; - } - - long double worst_error = 0.0L; - int worst_power = 0; + ASSERT_GE(rule.polynomial_exactness(), 0); for (int power = 0; power <= rule.polynomial_exactness(); @@ -199,14 +194,8 @@ std::pair expect_advertised_line_exactness( const long double error = std::abs( accumulate_line_moment(rule, nonnegative_power) - analytic_line_monomial_integral(nonnegative_power)); - if (error > worst_error) { - worst_error = error; - worst_power = power; - } EXPECT_LE(error, static_cast(tolerance)); } - - return {worst_error, worst_power}; } void expect_every_supported_line_rule( @@ -216,18 +205,6 @@ void expect_every_supported_line_rule( int exactness_subtrahend, LineEndpointPolicy endpoint_policy) { - long double worst_structure_error = 0.0L; - int worst_structure_num_points = 0; - std::size_t worst_structure_point_index = 0u; - std::string_view worst_structure_component = "none"; - - long double worst_measure_error = 0.0L; - int worst_measure_num_points = 0; - - long double worst_moment_error = 0.0L; - int worst_moment_num_points = 0; - int worst_moment_power = 0; - for (int num_points = first_num_points; num_points <= last_num_points; ++num_points) { @@ -240,115 +217,41 @@ void expect_every_supported_line_rule( static_cast(num_points), 2 * num_points - exactness_subtrahend); expect_line_rule_invariants(rule, endpoint_policy); - const auto [rule_moment_error, rule_moment_power] = - expect_advertised_line_exactness(rule); - if (worst_moment_num_points == 0 || - rule_moment_error > worst_moment_error) { - worst_moment_error = rule_moment_error; - worst_moment_num_points = num_points; - worst_moment_power = rule_moment_power; - } - - const long double measure_error = std::abs( - accumulate_line_moment(rule, 0u) - - static_cast(rule.reference_cell_measure())); - if (worst_measure_num_points == 0 || - measure_error > worst_measure_error) { - worst_measure_error = measure_error; - worst_measure_num_points = num_points; - } - - const auto update_worst_structure = [ - &worst_structure_error, - &worst_structure_num_points, - &worst_structure_point_index, - &worst_structure_component, - num_points]( - long double error, - std::size_t point_index, - std::string_view component) { - if (worst_structure_num_points == 0 || - error > worst_structure_error) { - worst_structure_error = error; - worst_structure_num_points = num_points; - worst_structure_point_index = point_index; - worst_structure_component = component; - } - }; - - for (std::size_t point_index = 0; - point_index < rule.num_points(); - ++point_index) { - const std::size_t mirror_index = - rule.num_points() - 1u - point_index; - update_worst_structure( - std::abs(static_cast( - rule.point(point_index)[1])), - point_index, - "inactive y coordinate"); - update_worst_structure( - std::abs(static_cast( - rule.point(point_index)[2])), - point_index, - "inactive z coordinate"); - update_worst_structure( - std::abs( - static_cast( - rule.point(point_index)[0]) + - static_cast( - rule.point(mirror_index)[0])), - point_index, - "mirrored point"); - update_worst_structure( - std::abs( - static_cast(rule.weight(point_index)) - - static_cast(rule.weight(mirror_index))), - point_index, - "mirrored weight"); - } + expect_advertised_line_exactness(rule); + } +} - if (rule.num_points() % 2u == 1u) { - const std::size_t center_index = rule.num_points() / 2u; - update_worst_structure( - std::abs(static_cast( - rule.point(center_index)[0])), - center_index, - "odd-rule center"); - } +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); - if (endpoint_policy == LineEndpointPolicy::Included) { - update_worst_structure( - std::abs( - static_cast(rule.point(0u)[0]) + 1.0L), - 0u, - "left endpoint"); - const std::size_t right_endpoint_index = - rule.num_points() - 1u; - update_worst_structure( - std::abs( - static_cast( - rule.point(right_endpoint_index)[0]) - - 1.0L), - right_endpoint_index, - "right endpoint"); + // 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); } - EXPECT_LE( - worst_structure_error, - static_cast(kStructureTolerance)) - << "worst num_points=" << worst_structure_num_points - << ", point index=" << worst_structure_point_index - << ", component=" << worst_structure_component; - EXPECT_LE( - worst_measure_error, - static_cast(kStructureTolerance)) - << "worst num_points=" << worst_measure_num_points; - EXPECT_LE( - worst_moment_error, - static_cast(kMomentTolerance)) - << "worst num_points=" << worst_moment_num_points - << ", power=" << worst_moment_power; + 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 @@ -425,70 +328,20 @@ TEST(QuadratureGeneratorTestSupport, ExercisesSharedLineRuleChecks) TEST(GaussLegendreImplementation, GeneratesCanonicalLowOrderRules) { - const auto expect_canonical_rule = []( - int num_points, - const auto& expected_points, - const auto& expected_weights) { - SCOPED_TRACE( - ::testing::Message() << "num_points=" << num_points); - const QuadratureRule rule = make_gauss_legendre_rule(num_points); - - expect_common_line_metadata( - rule, - static_cast(num_points), - 2 * num_points - 1); - expect_line_rule_invariants( - rule, - LineEndpointPolicy::Excluded); - expect_advertised_line_exactness(rule); - - ASSERT_EQ(rule.num_points(), expected_points.size()); - ASSERT_EQ(rule.num_points(), expected_weights.size()); - for (std::size_t point_index = 0; - point_index < rule.num_points(); - ++point_index) { - SCOPED_TRACE( - ::testing::Message() << "point index=" << point_index); - if (expected_points[point_index] == 0.0) { - EXPECT_DOUBLE_EQ(rule.point(point_index)[0], 0.0); - } else { - EXPECT_NEAR( - rule.point(point_index)[0], - expected_points[point_index], - kFixtureTolerance); - } - EXPECT_NEAR( - rule.weight(point_index), - expected_weights[point_index], - kFixtureTolerance); - } - - const std::size_t first_unadvertised_even_power = - static_cast(2 * num_points); - 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)); - }; - expect_canonical_rule( - 1, + make_gauss_legendre_rule(1), 1, std::array{0.0}, std::array{2.0}); const double two_point_abscissa = 1.0 / std::sqrt(3.0); expect_canonical_rule( - 2, + 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( - 3, + make_gauss_legendre_rule(3), 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}); @@ -528,75 +381,18 @@ TEST(GaussLegendreImplementation, RejectsRequestsOutsideSupportedRange) TEST(GaussLobattoImplementation, GeneratesCanonicalLowOrderRules) { - const auto expect_canonical_rule = []( - int num_points, - const auto& expected_points, - const auto& expected_weights) { - SCOPED_TRACE( - ::testing::Message() << "num_points=" << num_points); - const QuadratureRule rule = - make_gauss_lobatto_rule(num_points); - - expect_common_line_metadata( - rule, - static_cast(num_points), - 2 * num_points - 3); - expect_line_rule_invariants( - rule, - LineEndpointPolicy::Included); - expect_advertised_line_exactness(rule); - - ASSERT_EQ(rule.num_points(), expected_points.size()); - ASSERT_EQ(rule.num_points(), expected_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 expected_coordinate = - expected_points[point_index]; - if (expected_coordinate == -1.0 || - expected_coordinate == 0.0 || - expected_coordinate == 1.0) { - EXPECT_DOUBLE_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(2 * num_points - 2); - 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)); - }; - expect_canonical_rule( - 2, + make_gauss_lobatto_rule(2), 1, std::array{-1.0, 1.0}, std::array{1.0, 1.0}); expect_canonical_rule( - 3, + make_gauss_lobatto_rule(3), 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( - 4, + 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}); From 94f9ca43f09b2f5e313f2c65598a4db0fac8b3d4 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Sun, 13 Sep 2026 19:33:15 -0700 Subject: [PATCH 15/27] Simplify quadrature test helpers and comparisons --- .../Quadrature/test_QuadratureGenerators.cpp | 33 ++++++++----------- 1 file changed, 13 insertions(+), 20 deletions(-) diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp index 149407f11..4e81dc007 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp @@ -118,8 +118,7 @@ void expect_common_line_metadata( void expect_line_rule_invariants( const QuadratureRule& rule, - LineEndpointPolicy endpoint_policy, - double tolerance = kStructureTolerance) + LineEndpointPolicy endpoint_policy) { ASSERT_GT(rule.num_points(), 0u); ASSERT_EQ(rule.points().size(), rule.weights().size()); @@ -159,8 +158,8 @@ void expect_line_rule_invariants( EXPECT_NEAR( coordinate, -rule.point(mirror_index)[0], - tolerance); - EXPECT_NEAR(weight, rule.weight(mirror_index), tolerance); + kStructureTolerance); + EXPECT_NEAR(weight, rule.weight(mirror_index), kStructureTolerance); } if (endpoint_policy == LineEndpointPolicy::Included) { @@ -176,12 +175,10 @@ void expect_line_rule_invariants( const long double measure_error = std::abs( accumulate_line_moment(rule, 0u) - static_cast(rule.reference_cell_measure())); - EXPECT_LE(measure_error, static_cast(tolerance)); + EXPECT_LE(measure_error, static_cast(kStructureTolerance)); } -void expect_advertised_line_exactness( - const QuadratureRule& rule, - double tolerance = kMomentTolerance) +void expect_advertised_line_exactness(const QuadratureRule& rule) { ASSERT_GE(rule.polynomial_exactness(), 0); @@ -194,7 +191,7 @@ void expect_advertised_line_exactness( const long double error = std::abs( accumulate_line_moment(rule, nonnegative_power) - analytic_line_monomial_integral(nonnegative_power)); - EXPECT_LE(error, static_cast(tolerance)); + EXPECT_LE(error, static_cast(kMomentTolerance)); } } @@ -254,19 +251,18 @@ void expect_canonical_rule( EXPECT_GT(first_unadvertised_error, static_cast(kMomentTolerance)); } -template -void expect_exception_with_message( +template +void expect_invalid_argument_with_message( Function&& function, std::string_view expected_substring) { - static_assert(std::is_base_of_v); ASSERT_FALSE(expected_substring.empty()); try { std::forward(function)(); - FAIL() << "Expected requested exception containing: " + FAIL() << "Expected InvalidArgumentException containing: " << expected_substring; - } catch (const ExceptionType& exception) { + } catch (const InvalidArgumentException& exception) { const std::string_view actual_message{exception.what()}; EXPECT_NE( actual_message.find(expected_substring), @@ -318,7 +314,7 @@ TEST(QuadratureGeneratorTestSupport, ExercisesSharedLineRuleChecks) endpoint_rule, LineEndpointPolicy::Included); expect_advertised_line_exactness(endpoint_rule); - expect_exception_with_message( + expect_invalid_argument_with_message( [] { (void)QuadratureRule( svmp::CellFamily::Line, 1, {}, {}); @@ -371,7 +367,7 @@ TEST(GaussLegendreImplementation, RejectsRequestsOutsideSupportedRange) for (const int num_points : invalid_point_counts) { SCOPED_TRACE( ::testing::Message() << "num_points=" << num_points); - expect_exception_with_message( + expect_invalid_argument_with_message( [num_points] { (void)make_gauss_legendre_rule(num_points); }, @@ -424,7 +420,7 @@ TEST(GaussLobattoImplementation, RejectsRequestsOutsideSupportedRange) for (const int num_points : invalid_point_counts) { SCOPED_TRACE( ::testing::Message() << "num_points=" << num_points); - expect_exception_with_message( + expect_invalid_argument_with_message( [num_points] { (void)make_gauss_lobatto_rule(num_points); }, @@ -461,16 +457,13 @@ TEST(GaussLobattoBasisConsistency, MatchesRepresentativeNodeDistributions) if (point_index == 0u) { EXPECT_EQ(quadrature_coordinate, -1.0); EXPECT_EQ(basis_coordinate, -1.0); - EXPECT_EQ(quadrature_coordinate, basis_coordinate); } else if (point_index + 1u == rule.num_points()) { EXPECT_EQ(quadrature_coordinate, 1.0); EXPECT_EQ(basis_coordinate, 1.0); - EXPECT_EQ(quadrature_coordinate, basis_coordinate); } else if (num_points % 2 == 1 && point_index == rule.num_points() / 2u) { EXPECT_EQ(quadrature_coordinate, 0.0); EXPECT_EQ(basis_coordinate, 0.0); - EXPECT_EQ(quadrature_coordinate, basis_coordinate); } else { EXPECT_NEAR( quadrature_coordinate, From 311301cd0d3d4ca1e469e1a23f3f465f8aff6028 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Mon, 14 Sep 2026 15:02:39 -0700 Subject: [PATCH 16/27] Use requested exactness for Gaussian quadrature generators --- .../FE/Quadrature/GaussLobattoQuadrature.cpp | 16 ++-- .../FE/Quadrature/GaussLobattoQuadrature.h | 42 +++++----- .../solver/FE/Quadrature/GaussQuadrature.cpp | 16 ++-- .../solver/FE/Quadrature/GaussQuadrature.h | 37 +++++---- .../Quadrature/test_QuadratureGenerators.cpp | 80 ++++++++++--------- 5 files changed, 102 insertions(+), 89 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp index 45138fc51..4c7f8e27a 100644 --- a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp +++ b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp @@ -3,7 +3,7 @@ /** * @file GaussLobattoQuadrature.cpp - * @brief Bounded generation and validation of Gauss-Lobatto-Legendre line rules. + * @brief Exactness-requested generation of bounded Gauss-Lobatto-Legendre line rules. * @ingroup FE_Quadrature */ @@ -24,6 +24,9 @@ namespace svmp::FE::quadrature { namespace { +constexpr int kMaximumPoints = 128; +static_assert(max_gauss_lobatto_exactness() == 2 * kMaximumPoints - 3); + // Defensively bound supported cosine-seeded Newton refinements for // deterministic termination. constexpr int kMaximumNewtonIterations = 100; @@ -34,7 +37,7 @@ constexpr double kNewtonCorrectionTolerance = // Provide conservative O(n epsilon) accumulation headroom, qualified by those // sweeps. constexpr double kRuleValidationTolerance = - 32.0 * static_cast(max_gauss_lobatto_points()) * + 32.0 * static_cast(kMaximumPoints) * std::numeric_limits::epsilon(); std::pair evaluate_adjacent_legendre_values( @@ -188,16 +191,17 @@ std::pair generate_interior_root_and_weight(int num_points, int } // namespace -QuadratureRule make_gauss_lobatto_rule(int num_points) +QuadratureRule make_gauss_lobatto_rule(int requested_exactness) { - if (num_points < 2 || num_points > max_gauss_lobatto_points()) { + if (requested_exactness < 0 || requested_exactness > max_gauss_lobatto_exactness()) { std::ostringstream message; message << "Gauss-Lobatto-Legendre generator: " - << "num_points must be in [2, " - << max_gauss_lobatto_points() << ']'; + << "requested_exactness must be in [0, " + << max_gauss_lobatto_exactness() << ']'; svmp::raise(message.str()); } + const int num_points = requested_exactness / 2 + 2; std::vector points( static_cast(num_points), QuadPoint::Zero()); std::vector weights(points.size()); diff --git a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.h b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.h index a12efddbc..7337d69d3 100644 --- a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.h +++ b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.h @@ -19,39 +19,39 @@ namespace svmp::FE::quadrature { */ /** - * @brief Return the largest supported Gauss-Lobatto-Legendre point count. - * @details The 128 total points include both endpoints. This project support - * bound limits generator work and downstream product-rule growth while - * providing endpoint-inclusive line exactness through degree 253; it is not a - * mathematical or convergence limit. - * @return The inclusive point-count limit, 128. + * @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_points() noexcept +[[nodiscard]] constexpr int max_gauss_lobatto_exactness() noexcept { - return 128; + return 253; } /** - * @brief Generate an @p num_points Gauss-Lobatto-Legendre rule on - * @f$[-1,1]@f$. + * @brief Generate a Gauss-Lobatto-Legendre rule on @f$[-1,1]@f$ with at least + * the requested exactness. * - * @details The returned line rule has exactly @f$-1@f$ as its first point and - * exactly @f$+1@f$ as its last point. When present, its @f$n-2@f$ interior - * points are the roots of @f$P'_{n-1}@f$, where @f$n@f$ is @p num_points. - * Points are strictly increasing, weights are positive and aligned with their - * points, and the rule has polynomial exactness @f$2n-3@f$. + * @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$. Points are strictly increasing, and weights + * are positive and aligned with their points. * - * @param num_points Signed number of quadrature points; must be in the - * inclusive range - * @f$[2,\texttt{max\_gauss\_lobatto\_points()}]@f$. + * @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 num_points is outside the supported - * range. + * @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 num_points); +make_gauss_lobatto_rule(int requested_exactness); /** @} */ diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp index 49ba8c748..ca7db18b0 100644 --- a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp @@ -3,7 +3,7 @@ /** * @file GaussQuadrature.cpp - * @brief Bounded generation and validation of Gauss-Legendre line rules. + * @brief Exactness-requested generation of bounded Gauss-Legendre line rules. * @ingroup FE_Quadrature */ @@ -24,6 +24,9 @@ namespace svmp::FE::quadrature { namespace { +constexpr int 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; @@ -34,7 +37,7 @@ constexpr double kNewtonCorrectionTolerance = // Provide conservative O(n epsilon) accumulation headroom, qualified by those // sweeps. constexpr double kRuleValidationTolerance = - 32.0 * static_cast(max_gauss_legendre_points()) * + 32.0 * static_cast(kMaximumPoints) * std::numeric_limits::epsilon(); std::pair evaluate_legendre_with_derivative( @@ -173,15 +176,16 @@ std::pair generate_root_and_weight(int num_points, int root_inde } // namespace -QuadratureRule make_gauss_legendre_rule(int num_points) +QuadratureRule make_gauss_legendre_rule(int requested_exactness) { - if (num_points < 1 || num_points > max_gauss_legendre_points()) { + if (requested_exactness < 0 || requested_exactness > max_gauss_legendre_exactness()) { std::ostringstream message; - message << "Gauss-Legendre generator: num_points must be in [1, " - << max_gauss_legendre_points() << ']'; + message << "Gauss-Legendre generator: requested_exactness must be in [0, " + << max_gauss_legendre_exactness() << ']'; svmp::raise(message.str()); } + const int num_points = requested_exactness / 2 + 1; std::vector points( static_cast(num_points), QuadPoint::Zero()); std::vector weights(points.size()); diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.h b/Code/Source/solver/FE/Quadrature/GaussQuadrature.h index 14cdd376d..66acdd8a5 100644 --- a/Code/Source/solver/FE/Quadrature/GaussQuadrature.h +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.h @@ -19,35 +19,38 @@ namespace svmp::FE::quadrature { */ /** - * @brief Return the largest supported Gauss-Legendre point count. - * @details The 128-point ceiling is a project support bound that limits - * generator work and downstream product-rule growth while providing line - * exactness through degree 255. - * @return The inclusive point-count limit, 128. + * @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_points() noexcept +[[nodiscard]] constexpr int max_gauss_legendre_exactness() noexcept { - return 128; + return 255; } /** - * @brief Generate an @p num_points Gauss-Legendre rule on @f$[-1,1]@f$. + * @brief Generate a Gauss-Legendre rule on @f$[-1,1]@f$ with at least the + * requested exactness. * - * @details The returned line rule contains the roots of @f$P_n@f$ in strictly - * increasing order, where @f$n@f$ is @p num_points. Every point lies strictly - * inside @f$(-1,1)@f$, so neither endpoint is included. The rule has polynomial - * exactness @f$2n-1@f$ and positive weights aligned with its points. + * @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 + * @f$P_n@f$ in strictly increasing order inside @f$(-1,1)@f$; neither endpoint + * is included. Weights are positive and aligned with their points. * - * @param num_points Number of quadrature points; must be in - * @f$[1,\texttt{max\_gauss\_legendre\_points()}]@f$. + * @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 num_points is outside the supported - * range. + * @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 num_points); +make_gauss_legendre_rule(int requested_exactness); /** @} */ diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp index 4e81dc007..2b37a1800 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp @@ -3,7 +3,7 @@ /** * @file test_QuadratureGenerators.cpp - * @brief Tests for bounded one-dimensional quadrature generators. + * @brief Tests for exactness-requested one-dimensional quadrature generators. */ #include @@ -42,13 +42,13 @@ constexpr double kFixtureTolerance = constexpr double kBasisConsistencyTolerance = 128.0 * std::numeric_limits::epsilon(); -static_assert(max_gauss_legendre_points() == 128); -static_assert(noexcept(max_gauss_legendre_points())); +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_points() == 128); -static_assert(noexcept(max_gauss_lobatto_points())); +static_assert(max_gauss_lobatto_exactness() == 253); +static_assert(noexcept(max_gauss_lobatto_exactness())); static_assert( std::is_same_v); @@ -207,14 +207,20 @@ void expect_every_supported_line_rule( ++num_points) { SCOPED_TRACE( ::testing::Message() << "num_points=" << num_points); - const QuadratureRule rule = generator(num_points); - - expect_common_line_metadata( - rule, - static_cast(num_points), - 2 * num_points - exactness_subtrahend); - expect_line_rule_invariants(rule, endpoint_policy); - expect_advertised_line_exactness(rule); + const int expected_exactness = 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, static_cast(num_points), expected_exactness); + expect_line_rule_invariants(rule, endpoint_policy); + expect_advertised_line_exactness(rule); + } } } @@ -325,7 +331,7 @@ TEST(QuadratureGeneratorTestSupport, ExercisesSharedLineRuleChecks) TEST(GaussLegendreImplementation, GeneratesCanonicalLowOrderRules) { expect_canonical_rule( - make_gauss_legendre_rule(1), 1, + make_gauss_legendre_rule(0), 1, std::array{0.0}, std::array{2.0}); @@ -337,7 +343,7 @@ TEST(GaussLegendreImplementation, GeneratesCanonicalLowOrderRules) const double three_point_abscissa = std::sqrt(3.0 / 5.0); expect_canonical_rule( - make_gauss_legendre_rule(3), 5, + 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}); @@ -348,7 +354,7 @@ TEST(GaussLegendreImplementation, GeneratesEverySupportedRule) expect_every_supported_line_rule( &make_gauss_legendre_rule, 1, - max_gauss_legendre_points(), + 128, 1, LineEndpointPolicy::Excluded); } @@ -356,20 +362,19 @@ TEST(GaussLegendreImplementation, GeneratesEverySupportedRule) TEST(GaussLegendreImplementation, RejectsRequestsOutsideSupportedRange) { constexpr std::string_view expected_message = - "num_points must be in [1, 128]"; - constexpr std::array invalid_point_counts{ + "requested_exactness must be in [0, 255]"; + constexpr std::array invalid_exactness{ std::numeric_limits::min(), -1, - 0, - 129, + 256, std::numeric_limits::max()}; - for (const int num_points : invalid_point_counts) { + for (const int requested_exactness : invalid_exactness) { SCOPED_TRACE( - ::testing::Message() << "num_points=" << num_points); + ::testing::Message() << "requested_exactness=" << requested_exactness); expect_invalid_argument_with_message( - [num_points] { - (void)make_gauss_legendre_rule(num_points); + [requested_exactness] { + (void)make_gauss_legendre_rule(requested_exactness); }, expected_message); } @@ -378,11 +383,11 @@ TEST(GaussLegendreImplementation, RejectsRequestsOutsideSupportedRange) TEST(GaussLobattoImplementation, GeneratesCanonicalLowOrderRules) { expect_canonical_rule( - make_gauss_lobatto_rule(2), 1, + 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(3), 3, + 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}); @@ -399,7 +404,7 @@ TEST(GaussLobattoImplementation, GeneratesEverySupportedRule) expect_every_supported_line_rule( &make_gauss_lobatto_rule, 2, - max_gauss_lobatto_points(), + 128, 3, LineEndpointPolicy::Included); } @@ -408,21 +413,19 @@ TEST(GaussLobattoImplementation, RejectsRequestsOutsideSupportedRange) { constexpr std::string_view expected_message = "Gauss-Lobatto-Legendre generator: " - "num_points must be in [2, 128]"; - constexpr std::array invalid_point_counts{ + "requested_exactness must be in [0, 253]"; + constexpr std::array invalid_exactness{ std::numeric_limits::min(), -1, - 0, - 1, - 129, + 254, std::numeric_limits::max()}; - for (const int num_points : invalid_point_counts) { + for (const int requested_exactness : invalid_exactness) { SCOPED_TRACE( - ::testing::Message() << "num_points=" << num_points); + ::testing::Message() << "requested_exactness=" << requested_exactness); expect_invalid_argument_with_message( - [num_points] { - (void)make_gauss_lobatto_rule(num_points); + [requested_exactness] { + (void)make_gauss_lobatto_rule(requested_exactness); }, expected_message); } @@ -430,14 +433,13 @@ TEST(GaussLobattoImplementation, RejectsRequestsOutsideSupportedRange) TEST(GaussLobattoBasisConsistency, MatchesRepresentativeNodeDistributions) { - constexpr std::array point_counts{ - 2, 4, 65, max_gauss_lobatto_points()}; + constexpr std::array point_counts{2, 4, 65, 128}; for (const int num_points : point_counts) { SCOPED_TRACE( ::testing::Message() << "num_points=" << num_points); const QuadratureRule rule = - make_gauss_lobatto_rule(num_points); + make_gauss_lobatto_rule(2 * num_points - 3); ASSERT_EQ( rule.num_points(), static_cast(num_points)); From 7b2fab2ea918b9741e7c01a244f0de856508142f Mon Sep 17 00:00:00 2001 From: Zachary Sexton <47196674+zasexton@users.noreply.github.com> Date: Tue, 22 Sep 2026 14:30:14 -0700 Subject: [PATCH 17/27] Update Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp Co-authored-by: Michele Bucelli --- .../solver/FE/Quadrature/GaussQuadrature.cpp | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp index ca7db18b0..79c739693 100644 --- a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp @@ -178,12 +178,14 @@ std::pair generate_root_and_weight(int num_points, int root_inde QuadratureRule make_gauss_legendre_rule(int requested_exactness) { - if (requested_exactness < 0 || requested_exactness > max_gauss_legendre_exactness()) { - std::ostringstream message; - message << "Gauss-Legendre generator: requested_exactness must be in [0, " - << max_gauss_legendre_exactness() << ']'; - svmp::raise(message.str()); - } + svmp::check( + requested_exactness >= 0, + "requested_exactness cannot be negative."); + + svmp::check( + requested_exactness <= max_gauss_legendre_exactness(), + "requested exactness cannot be greater than " + + std::to_string(max_gauss_legendre_exactness() + "."); const int num_points = requested_exactness / 2 + 1; std::vector points( From 4bfc923755b56f30689948a254b5134485692c8a Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Wed, 23 Sep 2026 16:40:51 -0700 Subject: [PATCH 18/27] Fix Gauss-Legendre diagnostic string concatenation --- Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp index 79c739693..d42458a2d 100644 --- a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp @@ -185,7 +185,7 @@ QuadratureRule make_gauss_legendre_rule(int requested_exactness) svmp::check( requested_exactness <= max_gauss_legendre_exactness(), "requested exactness cannot be greater than " + - std::to_string(max_gauss_legendre_exactness() + "."); + std::to_string(max_gauss_legendre_exactness()) + "."); const int num_points = requested_exactness / 2 + 1; std::vector points( From 963ba69cccac4ed4e8c43dda7303175064aafcc5 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Wed, 23 Sep 2026 16:46:54 -0700 Subject: [PATCH 19/27] Restore Gauss-Legendre supported-range diagnostics --- Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp index d42458a2d..053669290 100644 --- a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp @@ -179,13 +179,10 @@ std::pair generate_root_and_weight(int num_points, int root_inde QuadratureRule make_gauss_legendre_rule(int requested_exactness) { svmp::check( - requested_exactness >= 0, - "requested_exactness cannot be negative."); - - svmp::check( - requested_exactness <= max_gauss_legendre_exactness(), - "requested exactness cannot be greater than " + - std::to_string(max_gauss_legendre_exactness()) + "."); + 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 int num_points = requested_exactness / 2 + 1; std::vector points( From ce871ec318dfa014767f83837dde0c858195953a Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 24 Sep 2026 13:02:53 -0700 Subject: [PATCH 20/27] Document Gaussian quadrature test coverage Addresses https://github.com/SimVascular/svMultiPhysics/pull/645#discussion_r4039200488 --- Documentation/Doxyfile | 1 + .../Quadrature/test_QuadratureGenerators.cpp | 40 +++++++++++++++++++ 2 files changed, 41 insertions(+) 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 index 2b37a1800..01bd3b23f 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp @@ -283,6 +283,11 @@ void expect_invalid_argument_with_message( } // 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); @@ -328,6 +333,11 @@ TEST(QuadratureGeneratorTestSupport, ExercisesSharedLineRuleChecks) "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( @@ -349,6 +359,11 @@ TEST(GaussLegendreImplementation, GeneratesCanonicalLowOrderRules) 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( @@ -359,6 +374,11 @@ TEST(GaussLegendreImplementation, GeneratesEverySupportedRule) 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 = @@ -380,6 +400,11 @@ TEST(GaussLegendreImplementation, RejectsRequestsOutsideSupportedRange) } } +/** + * @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( @@ -399,6 +424,11 @@ TEST(GaussLobattoImplementation, GeneratesCanonicalLowOrderRules) 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( @@ -409,6 +439,11 @@ TEST(GaussLobattoImplementation, GeneratesEverySupportedRule) 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 = @@ -431,6 +466,11 @@ TEST(GaussLobattoImplementation, RejectsRequestsOutsideSupportedRange) } } +/** + * @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}; From a6ca9bf56587cf37cea6f36c8efd5757a26c2056 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 24 Sep 2026 13:03:40 -0700 Subject: [PATCH 21/27] Define the Legendre polynomial in Gauss documentation Addresses https://github.com/SimVascular/svMultiPhysics/pull/645#discussion_r4039226389 --- Code/Source/solver/FE/Quadrature/GaussQuadrature.h | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.h b/Code/Source/solver/FE/Quadrature/GaussQuadrature.h index 66acdd8a5..b9d594a0e 100644 --- a/Code/Source/solver/FE/Quadrature/GaussQuadrature.h +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.h @@ -38,8 +38,11 @@ namespace svmp::FE::quadrature { * 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 - * @f$P_n@f$ in strictly increasing order inside @f$(-1,1)@f$; neither endpoint - * is included. Weights are positive and aligned with their points. + * 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()). From cfd0658ac52895f09ef9dd396944cd18fb2e559d Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 24 Sep 2026 13:03:40 -0700 Subject: [PATCH 22/27] Document Gauss-Legendre construction helpers Addresses https://github.com/SimVascular/svMultiPhysics/pull/645#discussion_r4039343915 --- Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp index 053669290..3c92853c6 100644 --- a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp @@ -40,6 +40,9 @@ constexpr double kRuleValidationTolerance = 32.0 * static_cast(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( int degree, double coordinate) noexcept @@ -67,6 +70,8 @@ std::pair evaluate_legendre_with_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( int num_points, int root_index, @@ -94,6 +99,10 @@ void require_generation(bool condition, int num_points, int root_index, int iter } } +// 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(int num_points, int root_index, bool is_center) { const double pi = std::numbers::pi_v; From b751f2d3420ca0aa252232dbcb67b2f7f8e569b0 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 24 Sep 2026 13:03:40 -0700 Subject: [PATCH 23/27] Simplify Gauss-Legendre floating-point conversions Addresses https://github.com/SimVascular/svMultiPhysics/pull/645#discussion_r4039357866 --- .../Source/solver/FE/Quadrature/GaussQuadrature.cpp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp index 3c92853c6..babac458e 100644 --- a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp @@ -37,7 +37,7 @@ constexpr double kNewtonCorrectionTolerance = // Provide conservative O(n epsilon) accumulation headroom, qualified by those // sweeps. constexpr double kRuleValidationTolerance = - 32.0 * static_cast(kMaximumPoints) * + 32.0 * kMaximumPoints * std::numeric_limits::epsilon(); // Return (P_degree(x), P'_degree(x)) for degree >= 1. Advance the Legendre @@ -52,14 +52,14 @@ std::pair evaluate_legendre_with_derivative( double value = coordinate; double derivative = 1.0; for (int recurrence_degree = 2; recurrence_degree <= degree; ++recurrence_degree) { - const double degree_value = static_cast(recurrence_degree); - const double recurrence_factor = static_cast(2 * recurrence_degree - 1); + const double degree_value = recurrence_degree; + const double recurrence_factor = 2 * recurrence_degree - 1; const double next_value = (recurrence_factor * coordinate * value - - static_cast(recurrence_degree - 1) * previous_value) / degree_value; + (recurrence_degree - 1) * previous_value) / degree_value; const double next_derivative = (recurrence_factor * (value + coordinate * derivative) - - static_cast(recurrence_degree - 1) * previous_derivative) / degree_value; + (recurrence_degree - 1) * previous_derivative) / degree_value; previous_value = value; previous_derivative = derivative; @@ -107,8 +107,7 @@ std::pair generate_root_and_weight(int num_points, int root_inde { const double pi = std::numbers::pi_v; double root = std::cos( - pi * (static_cast(root_index) + 0.75) / - (static_cast(num_points) + 0.5)); + pi * (root_index + 0.75) / (num_points + 0.5)); double correction = 0.0; for (int iteration = 1; iteration <= kMaximumNewtonIterations; ++iteration) { From a177daf3b072e611746b341edf39197757225ad8 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 24 Sep 2026 13:03:41 -0700 Subject: [PATCH 24/27] Preserve nonfinite Gauss-Legendre failure residuals Addresses https://github.com/SimVascular/svMultiPhysics/pull/645#discussion_r4039381701 --- Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp index babac458e..de3b0e1f1 100644 --- a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp @@ -85,10 +85,8 @@ std::pair evaluate_legendre_with_derivative( << ", root_index=" << root_index << ", diagnostic_value=" << diagnostic_value; - const double residual = std::isfinite(diagnostic_value) - ? std::abs(diagnostic_value) - : 0.0; - svmp::raise(message.str(), iteration, residual); + svmp::raise( + message.str(), iteration, std::abs(diagnostic_value)); } void require_generation(bool condition, int num_points, int root_index, int iteration, From eef3a4a500a2b1864327ea524b434c6c19d4a2f1 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 24 Sep 2026 13:03:41 -0700 Subject: [PATCH 25/27] Make Gauss-Legendre failure checks explicit Addresses https://github.com/SimVascular/svMultiPhysics/pull/645#discussion_r4039390307 --- .../solver/FE/Quadrature/GaussQuadrature.cpp | 107 +++++++++--------- 1 file changed, 54 insertions(+), 53 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp index de3b0e1f1..fa8527505 100644 --- a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp @@ -89,14 +89,6 @@ std::pair evaluate_legendre_with_derivative( message.str(), iteration, std::abs(diagnostic_value)); } -void require_generation(bool condition, int num_points, int root_index, int iteration, - double diagnostic_value, std::string_view detail) -{ - if (!condition) { - raise_generation_failure(num_points, root_index, iteration, diagnostic_value, detail); - } -} - // 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 @@ -111,19 +103,21 @@ std::pair generate_root_and_weight(int num_points, int root_inde for (int iteration = 1; iteration <= kMaximumNewtonIterations; ++iteration) { const auto [polynomial_value, polynomial_derivative] = evaluate_legendre_with_derivative(num_points, root); - require_generation( - std::isfinite(polynomial_value) && - std::isfinite(polynomial_derivative) && - polynomial_derivative != 0.0, - num_points, root_index, iteration, polynomial_derivative, - "encountered an invalid Legendre value or derivative"); + 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; - require_generation( - std::isfinite(correction) && std::isfinite(updated_root), - num_points, root_index, iteration, correction, - "computed an invalid Newton update"); + 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) { @@ -133,44 +127,49 @@ std::pair generate_root_and_weight(int num_points, int root_inde if (is_center) { root = 0.0; } - require_generation( - root >= 0.0 && root < 1.0 && - (is_center || root > 0.0), - num_points, root_index, iteration, root, - "refined root is outside the expected half interval"); + 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); - require_generation( - std::isfinite(final_polynomial_value) && - std::isfinite(final_polynomial_derivative) && - final_polynomial_derivative != 0.0, - num_points, root_index, iteration, final_polynomial_derivative, - "refined root produced an invalid Legendre value or derivative"); + 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; - require_generation( - std::isfinite(final_correction) && - std::abs(final_correction) <= - kNewtonCorrectionTolerance, - num_points, root_index, iteration, final_correction, - "refined root failed final correction validation"); + 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; - require_generation( - std::isfinite(denominator) && denominator > 0.0, - num_points, root_index, iteration, denominator, - "refined root produced an invalid weight denominator"); + 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; - require_generation( - std::isfinite(weight) && weight > 0.0, - num_points, root_index, iteration, weight, - "refined root produced an invalid quadrature weight"); + 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}; } @@ -210,21 +209,23 @@ QuadratureRule make_gauss_legendre_rule(int requested_exactness) 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]; - require_generation( - spacing > 0.0, - num_points, static_cast(point_index), -1, spacing, - "generated points are not strictly increasing"); + 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); - require_generation( - std::isfinite(weight_sum) && - measure_error <= - static_cast(kRuleValidationTolerance), - num_points, -1, -1, static_cast(measure_error), - "generated weights do not reproduce the reference measure"); + 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 = 2 * num_points - 1; return QuadratureRule( From eaa558559f2828e0d214915692c2ce4c1f4230f6 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 24 Sep 2026 13:03:42 -0700 Subject: [PATCH 26/27] Apply quadrature review updates to Lobatto rules Addresses https://github.com/SimVascular/svMultiPhysics/pull/645#discussion_r4040077818 --- .../FE/Quadrature/GaussLobattoQuadrature.cpp | 173 +++++++++--------- .../FE/Quadrature/GaussLobattoQuadrature.h | 7 +- 2 files changed, 93 insertions(+), 87 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp index 4c7f8e27a..ca1f5f1f4 100644 --- a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp +++ b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include @@ -37,9 +38,11 @@ constexpr double kNewtonCorrectionTolerance = // Provide conservative O(n epsilon) accumulation headroom, qualified by those // sweeps. constexpr double kRuleValidationTolerance = - 32.0 * static_cast(kMaximumPoints) * + 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( int degree, double coordinate) noexcept @@ -49,10 +52,8 @@ std::pair evaluate_adjacent_legendre_values( for (int recurrence_degree = 2; recurrence_degree <= degree; ++recurrence_degree) { const double next_value = - (static_cast(2 * recurrence_degree - 1) * - coordinate * value - - static_cast(recurrence_degree - 1) * previous_value) / - static_cast(recurrence_degree); + ((2 * recurrence_degree - 1) * coordinate * value - + (recurrence_degree - 1) * previous_value) / recurrence_degree; previous_value = value; value = next_value; } @@ -60,6 +61,8 @@ std::pair evaluate_adjacent_legendre_values( 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(int num_points, int half_root_index, int iteration, double diagnostic_value, std::string_view detail) { @@ -69,29 +72,22 @@ std::pair evaluate_adjacent_legendre_values( << ", half_root_index=" << half_root_index << ", diagnostic_value=" << diagnostic_value; - const double residual = std::isfinite(diagnostic_value) - ? std::abs(diagnostic_value) - : 0.0; - svmp::raise(message.str(), iteration, residual); -} - -void require_generation(bool condition, int num_points, int half_root_index, - int iteration, double diagnostic_value, std::string_view detail) -{ - if (!condition) { - raise_generation_failure( - num_points, half_root_index, iteration, diagnostic_value, detail); - } + 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(int num_points, int half_root_index, bool is_center, double weight_denominator_scale) { const int polynomial_degree = num_points - 1; - const double num_points_value = static_cast(num_points); - const double degree_value = static_cast(polynomial_degree); + 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 * static_cast(half_root_index + 1) / degree_value); + 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 @@ -99,29 +95,31 @@ std::pair generate_interior_root_and_weight(int num_points, int for (int iteration = 1; iteration <= kMaximumNewtonIterations; ++iteration) { const auto [polynomial_value, previous_polynomial_value] = evaluate_adjacent_legendre_values(polynomial_degree, root); - require_generation( - std::isfinite(polynomial_value) && - std::isfinite(previous_polynomial_value), - num_points, half_root_index, iteration, - polynomial_value, - "encountered invalid adjacent Legendre values"); + 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; - require_generation( - std::isfinite(residual) && - std::isfinite(derivative) && - derivative != 0.0, - num_points, half_root_index, iteration, derivative, - "computed an invalid root-function residual or derivative"); + 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; - require_generation( - std::isfinite(correction) && std::isfinite(updated_root), - num_points, half_root_index, iteration, correction, - "computed an invalid Newton update"); + 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) { @@ -131,55 +129,60 @@ std::pair generate_interior_root_and_weight(int num_points, int if (is_center) { root = 0.0; } - require_generation( - root >= 0.0 && root < 1.0 && (is_center || root > 0.0), - num_points, half_root_index, iteration, root, - "refined root is outside the expected half interval"); + 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); - require_generation( - std::isfinite(final_polynomial_value) && - std::isfinite(final_previous_polynomial_value), - num_points, half_root_index, iteration, - final_polynomial_value, - "refined root produced invalid adjacent Legendre values"); + 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; - require_generation( - std::isfinite(final_residual) && - std::isfinite(final_derivative) && - final_derivative != 0.0, - num_points, half_root_index, iteration, final_derivative, - "refined root produced an invalid residual or derivative"); + 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; - require_generation( - std::isfinite(final_correction) && - std::abs(final_correction) <= - kNewtonCorrectionTolerance, - num_points, half_root_index, iteration, final_correction, - "refined root failed final correction validation"); + 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; - require_generation( - std::isfinite(denominator) && denominator > 0.0, - num_points, half_root_index, iteration, denominator, - "refined root produced an invalid weight denominator"); + 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; - require_generation( - std::isfinite(weight) && weight > 0.0, - num_points, half_root_index, iteration, weight, - "refined root produced an invalid quadrature weight"); + 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}; } @@ -193,13 +196,11 @@ std::pair generate_interior_root_and_weight(int num_points, int QuadratureRule make_gauss_lobatto_rule(int requested_exactness) { - if (requested_exactness < 0 || requested_exactness > max_gauss_lobatto_exactness()) { - std::ostringstream message; - message << "Gauss-Lobatto-Legendre generator: " - << "requested_exactness must be in [0, " - << max_gauss_lobatto_exactness() << ']'; - svmp::raise(message.str()); - } + 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 int num_points = requested_exactness / 2 + 2; std::vector points( @@ -210,7 +211,7 @@ QuadratureRule make_gauss_lobatto_rule(int requested_exactness) points.back()[0] = 1.0; const double weight_denominator_scale = - static_cast(num_points * (num_points - 1)); + num_points * (num_points - 1); const double endpoint_weight = 2.0 / weight_denominator_scale; weights.front() = endpoint_weight; weights.back() = endpoint_weight; @@ -236,21 +237,23 @@ QuadratureRule make_gauss_lobatto_rule(int requested_exactness) 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]; - require_generation( - spacing > 0.0, - num_points, static_cast(point_index), -1, spacing, - "generated points are not strictly increasing"); + 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); - require_generation( - std::isfinite(weight_sum) && - measure_error <= - static_cast(kRuleValidationTolerance), - num_points, -1, -1, static_cast(measure_error), - "generated weights do not reproduce the reference measure"); + 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 = 2 * num_points - 3; return QuadratureRule( diff --git a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.h b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.h index 7337d69d3..d00b453be 100644 --- a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.h +++ b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.h @@ -39,8 +39,11 @@ namespace svmp::FE::quadrature { * 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$. Points are strictly increasing, and weights - * are positive and aligned with their 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()). From 0dcf5f9359b28c5884ede6295d1af8e2b6ec91be Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 24 Sep 2026 16:44:05 -0700 Subject: [PATCH 27/27] Use size_t for Gaussian quadrature point counts Propagate nonnegative counts through both generators and test helpers, keeping signed exactness and diagnostic conversions explicit. Addresses https://github.com/SimVascular/svMultiPhysics/pull/645#discussion_r4039325122 --- .../FE/Quadrature/GaussLobattoQuadrature.cpp | 35 +++++++++---------- .../solver/FE/Quadrature/GaussQuadrature.cpp | 27 +++++++------- .../Quadrature/test_QuadratureGenerators.cpp | 23 ++++++------ 3 files changed, 40 insertions(+), 45 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp index ca1f5f1f4..8065f7529 100644 --- a/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp +++ b/Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp @@ -25,8 +25,8 @@ namespace svmp::FE::quadrature { namespace { -constexpr int kMaximumPoints = 128; -static_assert(max_gauss_lobatto_exactness() == 2 * kMaximumPoints - 3); +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. @@ -44,13 +44,13 @@ constexpr double kRuleValidationTolerance = // 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( - int degree, + std::size_t degree, double coordinate) noexcept { double previous_value = 1.0; double value = coordinate; - for (int recurrence_degree = 2; recurrence_degree <= degree; ++recurrence_degree) { + 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; @@ -63,7 +63,7 @@ std::pair evaluate_adjacent_legendre_values( // 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(int num_points, int half_root_index, +[[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; @@ -80,10 +80,10 @@ std::pair evaluate_adjacent_legendre_values( // 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(int num_points, int half_root_index, - bool is_center, double weight_denominator_scale) +std::pair generate_interior_root_and_weight(std::size_t num_points, + int half_root_index, bool is_center, double weight_denominator_scale) { - const int polynomial_degree = num_points - 1; + 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; @@ -202,9 +202,8 @@ QuadratureRule make_gauss_lobatto_rule(int requested_exactness) "Gauss-Lobatto-Legendre generator: requested_exactness must be in [0, " + std::to_string(max_gauss_lobatto_exactness()) + ']'); - const int num_points = requested_exactness / 2 + 2; - std::vector points( - static_cast(num_points), QuadPoint::Zero()); + 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; @@ -216,17 +215,15 @@ QuadratureRule make_gauss_lobatto_rule(int requested_exactness) weights.front() = endpoint_weight; weights.back() = endpoint_weight; - const int interior_roots_to_refine = (num_points - 1) / 2; - for (int half_root_index = 0; + 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 + static_cast(half_root_index); + const std::size_t left_index = 1u + half_root_index; const std::size_t right_index = - points.size() - 2u - - static_cast(half_root_index); + points.size() - 2u - half_root_index; const auto [root, weight] = generate_interior_root_and_weight( - num_points, half_root_index, left_index == right_index, + num_points, static_cast(half_root_index), left_index == right_index, weight_denominator_scale); points[left_index][0] = -root; @@ -255,7 +252,7 @@ QuadratureRule make_gauss_lobatto_rule(int requested_exactness) "generated weights do not reproduce the reference measure"); } - const int polynomial_exactness = 2 * num_points - 3; + const int polynomial_exactness = static_cast(2 * num_points - 3); return QuadratureRule( svmp::CellFamily::Line, polynomial_exactness, std::move(points), std::move(weights)); diff --git a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp index fa8527505..0750365cf 100644 --- a/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp +++ b/Code/Source/solver/FE/Quadrature/GaussQuadrature.cpp @@ -24,8 +24,8 @@ namespace svmp::FE::quadrature { namespace { -constexpr int kMaximumPoints = 128; -static_assert(max_gauss_legendre_exactness() == 2 * kMaximumPoints - 1); +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. @@ -44,14 +44,14 @@ constexpr double kRuleValidationTolerance = // 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( - int degree, + 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 (int recurrence_degree = 2; recurrence_degree <= degree; ++recurrence_degree) { + 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 = @@ -73,7 +73,7 @@ std::pair evaluate_legendre_with_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( - int num_points, + std::size_t num_points, int root_index, int iteration, double diagnostic_value, @@ -93,7 +93,8 @@ std::pair evaluate_legendre_with_derivative( // 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(int num_points, int root_index, bool is_center) +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( @@ -189,17 +190,15 @@ QuadratureRule make_gauss_legendre_rule(int requested_exactness) "Gauss-Legendre generator: requested_exactness must be in [0, " + std::to_string(max_gauss_legendre_exactness()) + ']'); - const int num_points = requested_exactness / 2 + 1; - std::vector points( - static_cast(num_points), QuadPoint::Zero()); + 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 int roots_to_refine = (num_points + 1) / 2; - for (int root_index = 0; root_index < roots_to_refine; ++root_index) { - const std::size_t left_index = static_cast(root_index); + 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, root_index, left_index == right_index); + num_points, static_cast(left_index), left_index == right_index); points[left_index][0] = -root; points[right_index][0] = root; @@ -227,7 +226,7 @@ QuadratureRule make_gauss_legendre_rule(int requested_exactness) "generated weights do not reproduce the reference measure"); } - const int polynomial_exactness = 2 * num_points - 1; + const int polynomial_exactness = static_cast(2 * num_points - 1); return QuadratureRule( svmp::CellFamily::Line, polynomial_exactness, std::move(points), std::move(weights)); diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp index 01bd3b23f..53d100f62 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureGenerators.cpp @@ -197,17 +197,18 @@ void expect_advertised_line_exactness(const QuadratureRule& rule) void expect_every_supported_line_rule( QuadratureRule (*generator)(int), - int first_num_points, - int last_num_points, + std::size_t first_num_points, + std::size_t last_num_points, int exactness_subtrahend, LineEndpointPolicy endpoint_policy) { - for (int num_points = first_num_points; + 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 = 2 * num_points - exactness_subtrahend; + 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 : @@ -217,7 +218,7 @@ void expect_every_supported_line_rule( const QuadratureRule rule = generator(requested_exactness); expect_common_line_metadata( - rule, static_cast(num_points), expected_exactness); + rule, num_points, expected_exactness); expect_line_rule_invariants(rule, endpoint_policy); expect_advertised_line_exactness(rule); } @@ -473,16 +474,14 @@ TEST(GaussLobattoImplementation, RejectsRequestsOutsideSupportedRange) */ TEST(GaussLobattoBasisConsistency, MatchesRepresentativeNodeDistributions) { - constexpr std::array point_counts{2, 4, 65, 128}; + constexpr std::array point_counts{2, 4, 65, 128}; - for (const int num_points : point_counts) { + for (const std::size_t num_points : point_counts) { SCOPED_TRACE( ::testing::Message() << "num_points=" << num_points); const QuadratureRule rule = - make_gauss_lobatto_rule(2 * num_points - 3); - ASSERT_EQ( - rule.num_points(), - static_cast(num_points)); + 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(); @@ -494,7 +493,7 @@ TEST(GaussLobattoBasisConsistency, MatchesRepresentativeNodeDistributions) const double basis_coordinate = svmp::FE::basis::line_coord_pm_one( static_cast(point_index), - num_points - 1); + static_cast(num_points - 1)); if (point_index == 0u) { EXPECT_EQ(quadrature_coordinate, -1.0);