FE Quadrature Phase 02: bounded Gaussian line generators

Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.cpp

@@ -0,0 +1,261 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others.
+// SPDX-License-Identifier: BSD-3-Clause
+
+/**
+ * @file GaussLobattoQuadrature.cpp
+ * @brief Exactness-requested generation of bounded Gauss-Lobatto-Legendre line rules.
+ * @ingroup FE_Quadrature
+ */
+
+#include "FE/Quadrature/GaussLobattoQuadrature.h"
+
+#include "FE/Common/FEException.h"
+
+#include <cmath>
+#include <cstddef>
+#include <limits>
+#include <numbers>
+#include <numeric>
+#include <sstream>
+#include <string>
+#include <string_view>
+#include <utility>
+#include <vector>
+
+namespace svmp::FE::quadrature {
+namespace {
+
+constexpr std::size_t kMaximumPoints = 128;
+static_assert(max_gauss_lobatto_exactness() == static_cast<int>(2 * kMaximumPoints - 3));
+
+// Defensively bound supported cosine-seeded Newton refinements for
+// deterministic termination.
+constexpr int kMaximumNewtonIterations = 100;
+// Guard recurrence and Newton-update rounding; exhaustive supported-size
+// sweeps qualify this scale.
+constexpr double kNewtonCorrectionTolerance =
+    64.0 * std::numeric_limits<double>::epsilon();
+// Provide conservative O(n epsilon) accumulation headroom, qualified by those
+// sweeps.
+constexpr double kRuleValidationTolerance =
+    32.0 * kMaximumPoints *
+    std::numeric_limits<double>::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<double, double> evaluate_adjacent_legendre_values(
+    std::size_t degree,
+    double coordinate) noexcept
+{
+    double previous_value = 1.0;
+    double value = coordinate;
+
+    for (std::size_t recurrence_degree = 2; recurrence_degree <= degree; ++recurrence_degree) {
+        const double next_value =
+            ((2 * recurrence_degree - 1) * coordinate * value -
+             (recurrence_degree - 1) * previous_value) / recurrence_degree;
+        previous_value = value;
+        value = next_value;
+    }
+
+    return {value, previous_value};
+}
+
+// Preserve the failed quantity and generator context in a convergence error.
+// A half-root index or iteration of -1 identifies a rule-wide validation check.
+[[noreturn]] void raise_generation_failure(std::size_t num_points, int half_root_index,
+    int iteration, double diagnostic_value, std::string_view detail)
+{
+    std::ostringstream message;
+    message << "Gauss-Lobatto-Legendre generator: " << detail
+            << ", num_points=" << num_points
+            << ", half_root_index=" << half_root_index
+            << ", diagnostic_value=" << diagnostic_value;
+
+    svmp::raise<ConvergenceException>(
+        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<double, double> generate_interior_root_and_weight(std::size_t num_points,
+    int half_root_index, bool is_center, double weight_denominator_scale)
+{
+    const std::size_t polynomial_degree = num_points - 1;
+    const double num_points_value = num_points;
+    const double degree_value = polynomial_degree;
+    const double pi = std::numbers::pi_v<double>;
+    double root = std::cos(pi * (half_root_index + 1) / degree_value);
+    double correction = 0.0;
+
+    // For f = x*P_m - P_(m-1), Legendre identities give
+    // f' = (m+1)*P_m = n*P_m exactly.
+    for (int iteration = 1; iteration <= kMaximumNewtonIterations; ++iteration) {
+        const auto [polynomial_value, previous_polynomial_value] =
+            evaluate_adjacent_legendre_values(polynomial_degree, root);
+        if (!(std::isfinite(polynomial_value) &&
+              std::isfinite(previous_polynomial_value))) {
+            raise_generation_failure(
+                num_points, half_root_index, iteration, polynomial_value,
+                "encountered invalid adjacent Legendre values");
+        }
+
+        const double residual =
+            root * polynomial_value - previous_polynomial_value;
+        const double derivative = num_points_value * polynomial_value;
+        if (!(std::isfinite(residual) &&
+              std::isfinite(derivative) &&
+              derivative != 0.0)) {
+            raise_generation_failure(
+                num_points, half_root_index, iteration, derivative,
+                "computed an invalid root-function residual or derivative");
+        }
+
+        correction = residual / derivative;
+        const double updated_root = root - correction;
+        if (!(std::isfinite(correction) && std::isfinite(updated_root))) {
+            raise_generation_failure(
+                num_points, half_root_index, iteration, correction,
+                "computed an invalid Newton update");
+        }
+        root = updated_root;
+
+        if (std::abs(correction) > kNewtonCorrectionTolerance) {
+            continue;
+        }
+
+        if (is_center) {
+            root = 0.0;
+        }
+        if (!(root >= 0.0 && root < 1.0 && (is_center || root > 0.0))) {
+            raise_generation_failure(
+                num_points, half_root_index, iteration, root,
+                "refined root is outside the expected half interval");
+        }
+
+        const auto [final_polynomial_value,
+                    final_previous_polynomial_value] =
+            evaluate_adjacent_legendre_values(polynomial_degree, root);
+        if (!(std::isfinite(final_polynomial_value) &&
+              std::isfinite(final_previous_polynomial_value))) {
+            raise_generation_failure(
+                num_points, half_root_index, iteration, final_polynomial_value,
+                "refined root produced invalid adjacent Legendre values");
+        }
+
+        const double final_residual =
+            root * final_polynomial_value -
+            final_previous_polynomial_value;
+        const double final_derivative =
+            num_points_value * final_polynomial_value;
+        if (!(std::isfinite(final_residual) &&
+              std::isfinite(final_derivative) &&
+              final_derivative != 0.0)) {
+            raise_generation_failure(
+                num_points, half_root_index, iteration, final_derivative,
+                "refined root produced an invalid residual or derivative");
+        }
+
+        const double final_correction =
+            final_residual / final_derivative;
+        if (!(std::isfinite(final_correction) &&
+              std::abs(final_correction) <=
+                  kNewtonCorrectionTolerance)) {
+            raise_generation_failure(
+                num_points, half_root_index, iteration, final_correction,
+                "refined root failed final correction validation");
+        }
+
+        const double denominator =
+            weight_denominator_scale * final_polynomial_value *
+            final_polynomial_value;
+        if (!(std::isfinite(denominator) && denominator > 0.0)) {
+            raise_generation_failure(
+                num_points, half_root_index, iteration, denominator,
+                "refined root produced an invalid weight denominator");
+        }
+
+        const double weight = 2.0 / denominator;
+        if (!(std::isfinite(weight) && weight > 0.0)) {
+            raise_generation_failure(
+                num_points, half_root_index, iteration, weight,
+                "refined root produced an invalid quadrature weight");
+        }
+
+        return {root, weight};
+    }
+
+    raise_generation_failure(
+        num_points, half_root_index, kMaximumNewtonIterations, correction,
+        "Newton refinement did not converge");
+}
+
+} // namespace
+
+QuadratureRule make_gauss_lobatto_rule(int requested_exactness)
+{
+    svmp::check<InvalidArgumentException>(
+        requested_exactness >= 0 &&
+            requested_exactness <= max_gauss_lobatto_exactness(),
+        "Gauss-Lobatto-Legendre generator: requested_exactness must be in [0, " +
+            std::to_string(max_gauss_lobatto_exactness()) + ']');
+
+    const std::size_t num_points = static_cast<std::size_t>(requested_exactness) / 2 + 2;
+    std::vector<QuadPoint> points(num_points, QuadPoint::Zero());
+    std::vector<double> weights(points.size());
+
+    points.front()[0] = -1.0;
+    points.back()[0] = 1.0;
+
+    const double weight_denominator_scale =
+        num_points * (num_points - 1);
+    const double endpoint_weight = 2.0 / weight_denominator_scale;
+    weights.front() = endpoint_weight;
+    weights.back() = endpoint_weight;
+
+    const std::size_t interior_roots_to_refine = (num_points - 1) / 2;
+    for (std::size_t half_root_index = 0;
+         half_root_index < interior_roots_to_refine; ++half_root_index) {
+        const std::size_t left_index = 1u + half_root_index;
+        const std::size_t right_index =
+            points.size() - 2u - half_root_index;
+        const auto [root, weight] =
+            generate_interior_root_and_weight(
+                num_points, static_cast<int>(half_root_index), left_index == right_index,
+                weight_denominator_scale);
+
+        points[left_index][0] = -root;
+        points[right_index][0] = root;
+        weights[left_index] = weight;
+        weights[right_index] = weight;
+    }
+
+    for (std::size_t point_index = 1; point_index < points.size(); ++point_index) {
+        const double spacing = points[point_index][0] - points[point_index - 1u][0];
+        if (!(spacing > 0.0)) {
+            raise_generation_failure(
+                num_points, static_cast<int>(point_index), -1, spacing,
+                "generated points are not strictly increasing");
+        }
+    }
+
+    // Report a failed measure instead of repairing or rescaling the weights.
+    const long double weight_sum = std::accumulate(weights.begin(), weights.end(), 0.0L);
+    const long double measure_error = std::abs(weight_sum - 2.0L);
+    if (!(std::isfinite(weight_sum) &&
+          measure_error <=
+              static_cast<long double>(kRuleValidationTolerance))) {
+        raise_generation_failure(
+            num_points, -1, -1, static_cast<double>(measure_error),
+            "generated weights do not reproduce the reference measure");
+    }
+
+    const int polynomial_exactness = static_cast<int>(2 * num_points - 3);
+    return QuadratureRule(
+        svmp::CellFamily::Line, polynomial_exactness,
+        std::move(points), std::move(weights));
+}
+
+} // namespace svmp::FE::quadrature

Code/Source/solver/FE/Quadrature/GaussLobattoQuadrature.h

@@ -0,0 +1,63 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others.
+// SPDX-License-Identifier: BSD-3-Clause
+
+#ifndef SVMP_FE_GAUSS_LOBATTO_QUADRATURE_H
+#define SVMP_FE_GAUSS_LOBATTO_QUADRATURE_H
+
+/**
+ * @file GaussLobattoQuadrature.h
+ * @brief Bounded Gauss-Lobatto-Legendre quadrature generation on the reference line.
+ * @ingroup FE_Quadrature
+ */
+
+#include "FE/Quadrature/QuadratureRule.h"
+
+namespace svmp::FE::quadrature {
+
+/** @addtogroup FE_Quadrature
+ * @{
+ */
+
+/**
+ * @brief Return the largest supported Gauss-Lobatto-Legendre exactness request.
+ * @details Degree 253 corresponds to the internal 128-point project support
+ * bound, including both endpoints, which limits generator work and downstream
+ * product-rule growth; it is not a mathematical or convergence limit.
+ * @return The inclusive requested-exactness limit, 253.
+ */
+[[nodiscard]] constexpr int max_gauss_lobatto_exactness() noexcept
+{
+    return 253;
+}
+
+/**
+ * @brief Generate a Gauss-Lobatto-Legendre rule on @f$[-1,1]@f$ with at least
+ *        the requested exactness.
+ *
+ * @details For @f$d@f$ = @p requested_exactness, integer division gives the
+ * minimum point count @f$n=\lfloor d/2\rfloor+2@f$. The returned metadata reports
+ * the actual polynomial exactness @f$2n-3@f$, which exceeds even requests by one.
+ * Degree zero produces two points with exactness one. The first and last points
+ * are exactly @f$-1@f$ and @f$+1@f$. When present, the @f$n-2@f$ interior points
+ * are the roots of @f$P'_{n-1}@f$, the derivative of the Legendre polynomial
+ * of degree @f$n-1@f$. Points are strictly increasing, and weights are positive
+ * and aligned with their points.
+ *
+ * @see [NIST DLMF: Legendre polynomials](https://dlmf.nist.gov/18.3)
+ *
+ * @param requested_exactness Minimum polynomial degree to integrate exactly;
+ *        must be in [0, 253], inclusive (see max_gauss_lobatto_exactness()).
+ * @return A complete QuadratureRule value for CellFamily::Line.
+ * @throws InvalidArgumentException If @p requested_exactness is outside the
+ *         supported range; checked before point-count conversion or allocation.
+ * @throws ConvergenceException If root refinement or final numerical
+ *         validation fails.
+ */
+[[nodiscard]] QuadratureRule
+make_gauss_lobatto_rule(int requested_exactness);
+
+/** @} */
+
+} // namespace svmp::FE::quadrature
+
+#endif // SVMP_FE_GAUSS_LOBATTO_QUADRATURE_H