From a2f1677bc504758bbf9bd95f2a1a8b30a32a87b7 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Wed, 22 Jul 2026 10:55:04 -0700 Subject: [PATCH 01/21] Define the FE Quadrature module boundary --- Code/Source/solver/CMakeLists.txt | 6 ++++++ Code/Source/solver/FE/Common/Types.h | 28 +++++++--------------------- Code/Source/solver/FE/FE.h | 3 ++- 3 files changed, 15 insertions(+), 22 deletions(-) diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index bac65c976..691c030ba 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -258,11 +258,17 @@ file(GLOB SOLVER_FE_MATH_SRCS CONFIGURE_DEPENDS FE/Math/*.h ) +file(GLOB SOLVER_FE_QUADRATURE_SRCS CONFIGURE_DEPENDS + FE/Quadrature/*.cpp + FE/Quadrature/*.h +) + list(APPEND CSRCS ${SOLVER_CORE_SRCS} ${SOLVER_FE_COMMON_SRCS} ${SOLVER_FE_BASIS_SRCS} ${SOLVER_FE_MATH_SRCS} + ${SOLVER_FE_QUADRATURE_SRCS} ) # Set PETSc interace code. diff --git a/Code/Source/solver/FE/Common/Types.h b/Code/Source/solver/FE/Common/Types.h index f443d0225..74669f010 100644 --- a/Code/Source/solver/FE/Common/Types.h +++ b/Code/Source/solver/FE/Common/Types.h @@ -69,10 +69,9 @@ enum class CellFamily { * @brief Shared vocabulary types, constants, and exception infrastructure used by every FE module. * * @details The Common module collects the foundational definitions that the - * rest of the FE library builds on: index and scalar type aliases; element, - * basis, quadrature, and field enumerations; sentinel constants and strong - * type wrappers; and the FE exception hierarchy together with its - * argument-checking helpers. + * rest of the FE library builds on: index and scalar type aliases; shared + * enumerations and strong types; sentinel constants; and the FE exception + * hierarchy together with its argument-checking helpers. */ namespace svmp::FE { @@ -83,10 +82,10 @@ namespace svmp::FE { * @brief Core type aliases, enumerations, constants, geometric types, and compile-time traits. * * @details This group documents the index and identifier types used for - * element-local and global numbering, the element/basis/quadrature/field - * enumerations shared across modules, sentinel constants, reference- and - * physical-space geometric aliases, and the strong-type utilities that - * prevent accidental mixing of conceptually distinct values. + * element-local and global numbering, the enumerations shared across modules, + * sentinel constants, reference- and physical-space geometric aliases, and + * the strong-type utilities that prevent accidental mixing of conceptually + * distinct values. * @{ */ @@ -242,19 +241,6 @@ enum class ElementType : std::uint8_t { Unknown ///< Unrecognized or uninitialized element type }; -/** - * @brief Quadrature rule types - */ -enum class QuadratureType : std::uint8_t { - GaussLegendre, ///< Standard Gaussian quadrature - GaussLobatto, ///< Includes endpoints (for spectral elements) - Newton, ///< Newton-Cotes rules - Reduced, ///< Order-based reduced integration for locking - PositionBased, ///< Position-based reduced integration (legacy compatible) - Composite, ///< Composite rules for adaptivity - Custom ///< User-defined quadrature points -}; - /** * @brief Basis function families */ diff --git a/Code/Source/solver/FE/FE.h b/Code/Source/solver/FE/FE.h index 125660942..c51c4a787 100644 --- a/Code/Source/solver/FE/FE.h +++ b/Code/Source/solver/FE/FE.h @@ -11,7 +11,8 @@ * This header intentionally contains no declarations. It gives Doxygen a * header-based home for the top-level FE group; submodule groups attach to it * from their own headers, including FE_Basis (Basis/BasisFunction.h), - * FE_Common (Common/Types.h), and FE_Math (Math/Vector.h). + * FE_Common (Common/Types.h), FE_Math (Math/Vector.h), and FE_Quadrature + * (Quadrature/QuadratureRule.h). */ /** From a07b87b522f772e6b589360ed8ce657f167428c1 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Wed, 22 Jul 2026 10:55:10 -0700 Subject: [PATCH 02/21] Add immutable reference-space quadrature rules --- .../solver/FE/Quadrature/QuadratureRule.cpp | 386 +++++++++++++++++ .../solver/FE/Quadrature/QuadratureRule.h | 389 ++++++++++++++++++ 2 files changed, 775 insertions(+) create mode 100644 Code/Source/solver/FE/Quadrature/QuadratureRule.cpp create mode 100644 Code/Source/solver/FE/Quadrature/QuadratureRule.h diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp new file mode 100644 index 000000000..c99d72f2d --- /dev/null +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp @@ -0,0 +1,386 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +/** + * @file QuadratureRule.cpp + * @brief Internal construction and structural validation for quadrature rules. + * @ingroup FE_Quadrature + */ + +#include "FE/Quadrature/QuadratureRule.h" + +#include "FE/Common/FEException.h" + +#include +#include +#include +#include +#include + +namespace svmp::FE::quadrature { +namespace { + +enum class ReferenceDomain { + Unsupported, + Point, + Line, + Triangle, + Quad, + Tetra, + Hex, + Wedge, +}; + +struct ReferenceCellTraits { + ReferenceDomain domain; + int dimension; + double measure; +}; + +enum class ValidationFailure { + None, + UnsupportedCellFamily, + NegativePolynomialExactness, + InvalidTolerance, + EmptyRule, + PointWeightSizeMismatch, + NonFiniteCoordinate, + NonzeroInactiveCoordinate, + PointOutsideReferenceCell, + NonFiniteWeight, + NonFiniteWeightAccumulation, + IncorrectZerothMoment, + UnstableStoredOrderAccumulation, + IllConditionedWeightCancellation, +}; + +struct ValidationResult { + ValidationFailure failure{ValidationFailure::None}; + std::size_t sample{std::numeric_limits::max()}; +}; + +constexpr ReferenceCellTraits unsupported_traits() noexcept +{ + return { + ReferenceDomain::Unsupported, + -1, + std::numeric_limits::quiet_NaN(), + }; +} + +constexpr ReferenceCellTraits reference_cell_traits(svmp::CellFamily family) noexcept +{ + switch (family) { + case svmp::CellFamily::Point: + return {ReferenceDomain::Point, 0, 1.0}; + case svmp::CellFamily::Line: + return {ReferenceDomain::Line, 1, 2.0}; + case svmp::CellFamily::Triangle: + return {ReferenceDomain::Triangle, 2, 0.5}; + case svmp::CellFamily::Quad: + return {ReferenceDomain::Quad, 2, 4.0}; + case svmp::CellFamily::Tetra: + return {ReferenceDomain::Tetra, 3, 1.0 / 6.0}; + case svmp::CellFamily::Hex: + return {ReferenceDomain::Hex, 3, 8.0}; + case svmp::CellFamily::Wedge: + return {ReferenceDomain::Wedge, 3, 1.0}; + default: + return unsupported_traits(); + } +} + +ValidationResult validate_point( + const QuadPoint& point, + const ReferenceCellTraits& traits, + double tolerance, + std::size_t sample) noexcept +{ + for (std::size_t component = 0; component < 3u; ++component) { + if (!std::isfinite(point[component])) { + return {ValidationFailure::NonFiniteCoordinate, sample}; + } + if (component >= static_cast(traits.dimension) && + std::abs(point[component]) > tolerance) { + return {ValidationFailure::NonzeroInactiveCoordinate, sample}; + } + } + + const auto in_interval = [tolerance](double value, double lower, double upper) { + return value >= lower - tolerance && value <= upper + tolerance; + }; + + const double x = point[0]; + const double y = point[1]; + const double z = point[2]; + + bool is_contained = false; + switch (traits.domain) { + case ReferenceDomain::Point: + is_contained = true; + break; + case ReferenceDomain::Line: + is_contained = in_interval(x, -1.0, 1.0); + break; + case ReferenceDomain::Triangle: + is_contained = x >= -tolerance && y >= -tolerance && + x + y <= 1.0 + tolerance; + break; + case ReferenceDomain::Quad: + is_contained = in_interval(x, -1.0, 1.0) && + in_interval(y, -1.0, 1.0); + break; + case ReferenceDomain::Tetra: + is_contained = x >= -tolerance && y >= -tolerance && + z >= -tolerance && + x + y + z <= 1.0 + tolerance; + break; + case ReferenceDomain::Hex: + is_contained = in_interval(x, -1.0, 1.0) && + in_interval(y, -1.0, 1.0) && + in_interval(z, -1.0, 1.0); + break; + case ReferenceDomain::Wedge: + is_contained = x >= -tolerance && y >= -tolerance && + x + y <= 1.0 + tolerance && + in_interval(z, -1.0, 1.0); + break; + case ReferenceDomain::Unsupported: + break; + } + + if (!is_contained) { + return {ValidationFailure::PointOutsideReferenceCell, sample}; + } + return {}; +} + +double add_as_binary64(double left, double right) noexcept +{ + volatile double rounded_sum = left + right; + return rounded_sum; +} + +ValidationResult validate_weights( + const std::vector& weights, + double reference_measure, + double tolerance) noexcept +{ + double ordered_sum = 0.0; + long double compensated_sum = 0.0L; + long double correction = 0.0L; + long double absolute_sum = 0.0L; + bool has_negative_weight = false; + + for (std::size_t sample = 0; sample < weights.size(); ++sample) { + const double weight = weights[sample]; + if (!std::isfinite(weight)) { + return {ValidationFailure::NonFiniteWeight, sample}; + } + + ordered_sum = add_as_binary64(ordered_sum, weight); + if (!std::isfinite(ordered_sum)) { + return {ValidationFailure::NonFiniteWeightAccumulation, sample}; + } + has_negative_weight = has_negative_weight || weight < 0.0; + + const long double value = static_cast(weight); + const long double next = compensated_sum + value; + if (std::abs(compensated_sum) >= std::abs(value)) { + correction += (compensated_sum - next) + value; + } else { + correction += (value - next) + compensated_sum; + } + compensated_sum = next; + absolute_sum += std::abs(value); + if (!std::isfinite(compensated_sum) || + !std::isfinite(correction) || + !std::isfinite(absolute_sum)) { + return {ValidationFailure::NonFiniteWeightAccumulation, sample}; + } + } + + const long double total = compensated_sum + correction; + if (!std::isfinite(total)) { + return {ValidationFailure::NonFiniteWeightAccumulation, weights.size() - 1u}; + } + + const long double expected = static_cast(reference_measure); + const long double scale = std::max(1.0L, std::abs(expected)); + const long double requested_error_budget = + static_cast(tolerance) * scale; + const long double compensated_error = std::abs(total - expected); + + if (compensated_error > requested_error_budget) { + return {ValidationFailure::IncorrectZerothMoment}; + } + + const long double ordered_error = + std::abs(static_cast(ordered_sum) - expected); + if (ordered_error > requested_error_budget) { + return {ValidationFailure::UnstableStoredOrderAccumulation}; + } + + if (has_negative_weight) { + // The actual stored-order sum is checked above. This independent guard + // limits sensitivity to signed-weight cancellation without imposing a + // sample-count ceiling on otherwise well-conditioned rules. + constexpr long double cancellation_safety_factor = 8.0L; + const long double epsilon = + static_cast(std::numeric_limits::epsilon()); + const long double stability_error_bound = + compensated_error + + cancellation_safety_factor * epsilon * absolute_sum; + const long double stability_budget = + static_cast( + QuadratureRule::default_validation_tolerance()) * scale; + if (!std::isfinite(stability_error_bound) || + stability_error_bound > stability_budget) { + return {ValidationFailure::IllConditionedWeightCancellation}; + } + } + + return {}; +} + +ValidationResult validate_rule_data( + svmp::CellFamily family, + int polynomial_exactness, + const std::vector& points, + const std::vector& weights, + double tolerance) noexcept +{ + const auto traits = reference_cell_traits(family); + if (traits.domain == ReferenceDomain::Unsupported) { + return {ValidationFailure::UnsupportedCellFamily}; + } + if (polynomial_exactness < 0) { + return {ValidationFailure::NegativePolynomialExactness}; + } + if (!std::isfinite(tolerance) || tolerance < 0.0) { + return {ValidationFailure::InvalidTolerance}; + } + if (points.empty()) { + return {ValidationFailure::EmptyRule}; + } + if (points.size() != weights.size()) { + return {ValidationFailure::PointWeightSizeMismatch}; + } + + for (std::size_t sample = 0; sample < points.size(); ++sample) { + const auto result = + validate_point(points[sample], traits, tolerance, sample); + if (result.failure != ValidationFailure::None) { + return result; + } + } + + return validate_weights(weights, traits.measure, tolerance); +} + +std::string validation_failure_message(const ValidationResult& result) +{ + const auto sample_suffix = [&result]() { + if (result.sample == std::numeric_limits::max()) { + return std::string{}; + } + return std::string{" at sample "} + std::to_string(result.sample); + }; + + switch (result.failure) { + case ValidationFailure::None: + return {}; + case ValidationFailure::UnsupportedCellFamily: + return "QuadratureRule: unsupported reference-cell family"; + case ValidationFailure::NegativePolynomialExactness: + return "QuadratureRule: polynomial exactness must be non-negative"; + case ValidationFailure::InvalidTolerance: + return "QuadratureRule: validation tolerance must be finite and non-negative"; + case ValidationFailure::EmptyRule: + return "QuadratureRule: a rule must contain at least one sample"; + case ValidationFailure::PointWeightSizeMismatch: + return "QuadratureRule: points/weights size mismatch"; + case ValidationFailure::NonFiniteCoordinate: + return "QuadratureRule: quadrature point contains a non-finite coordinate" + + sample_suffix(); + case ValidationFailure::NonzeroInactiveCoordinate: + return "QuadratureRule: quadrature point has a nonzero inactive coordinate" + + sample_suffix(); + case ValidationFailure::PointOutsideReferenceCell: + return "QuadratureRule: quadrature point lies outside the canonical reference cell" + + sample_suffix(); + case ValidationFailure::NonFiniteWeight: + return "QuadratureRule: quadrature weight must be finite" + + sample_suffix(); + case ValidationFailure::NonFiniteWeightAccumulation: + return "QuadratureRule: weight accumulation is not finite" + + sample_suffix(); + case ValidationFailure::IncorrectZerothMoment: + return "QuadratureRule: weights do not reproduce the reference-cell measure"; + case ValidationFailure::UnstableStoredOrderAccumulation: + return "QuadratureRule: stored-order double-precision weight accumulation " + "does not reproduce the reference-cell measure"; + case ValidationFailure::IllConditionedWeightCancellation: + return "QuadratureRule: weight cancellation is too ill-conditioned for " + "stable double-precision integration"; + } + return "QuadratureRule: unknown validation failure"; +} + +} // namespace + +QuadratureRule::~QuadratureRule() = default; + +QuadratureRule::QuadratureRule(svmp::CellFamily family, RuleData data) + : QuadratureRule(validate(family, std::move(data))) +{ +} + +QuadratureRule::QuadratureRule(ValidatedState state) + : cell_family_(state.cell_family), + dimension_(state.dimension), + polynomial_exactness_(state.polynomial_exactness), + reference_measure_(state.reference_measure), + points_(std::move(state.points)), + weights_(std::move(state.weights)) +{ +} + +QuadratureRule::ValidatedState QuadratureRule::validate( + svmp::CellFamily family, + RuleData data) +{ + const auto validation = validate_rule_data( + family, + data.polynomial_exactness, + data.points, + data.weights, + default_validation_tolerance()); + if (validation.failure != ValidationFailure::None) { + svmp::raise( + validation_failure_message(validation)); + } + + const auto traits = reference_cell_traits(family); + return { + family, + traits.dimension, + data.polynomial_exactness, + traits.measure, + std::move(data.points), + std::move(data.weights), + }; +} + +bool QuadratureRule::is_structurally_valid(double tolerance) const noexcept +{ + return validate_rule_data( + cell_family_, + polynomial_exactness_, + points_, + weights_, + tolerance) + .failure == ValidationFailure::None; +} + +} // namespace svmp::FE::quadrature diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.h b/Code/Source/solver/FE/Quadrature/QuadratureRule.h new file mode 100644 index 000000000..7c2c15441 --- /dev/null +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.h @@ -0,0 +1,389 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others. +// SPDX-License-Identifier: BSD-3-Clause + +#ifndef SVMP_FE_QUADRATURE_RULE_H +#define SVMP_FE_QUADRATURE_RULE_H + +/** + * @file QuadratureRule.h + * @brief Immutable reference-space quadrature rule contract. + * @ingroup FE_Quadrature + * + * This header defines the consumer-facing representation of a finite-element + * quadrature rule. Rule construction and validation are implemented separately + * so consumers depend only on the stable query interface. + */ + +/** + * @defgroup FE_Quadrature Quadrature + * @ingroup FE + * @brief Immutable integration rules on canonical finite-element reference cells. + * + * @details + * ## Scope + * + * The Quadrature module owns ordered reference coordinates and weights used to + * approximate an integral on a canonical reference cell: + * @f[ + * \int_{\hat K} f(\hat x)\,d\hat x + * \approx \sum_q w_q f(\hat x_q). + * @f] + * A rule identifies its reference-cell family, reports its intrinsic dimension + * and declared polynomial exactness, and keeps every point paired with its + * corresponding weight. + * + * The module does not choose the exactness required by an equation term, apply + * reduced-integration policy, select a basis, own mesh storage, embed or orient + * a face in its parent element, or map a reference integral into physical + * space. Those operations require solver, basis, mesh, and geometry context and + * remain at the caller's integration boundary. + * + * ## Public API + * + * @ref svmp::FE::quadrature::QuadPoint "QuadPoint" and the const query surface + * of @ref svmp::FE::quadrature::QuadratureRule "QuadratureRule" form the public + * API. Integration consumers read rule metadata and ordered samples; they do + * not derive new rules or modify rule storage. + * + * Complete-data construction, reference-cell traits, point-containment checks, + * compensated weight summation, concrete generators, caches, and rule-selection + * facilities are module implementation details. + * + * ## Rule-provider contract + * + * Concrete rule providers are the only supported subclasses. A provider builds + * one complete RuleData payload before invoking the protected constructor and + * exposes no mutation afterward. It must advertise only exactness established + * for every rule it supplies through analytic moment tests. Derivation is a + * provider extension seam, not an integration-consumer customization point. + * + * ## Rule contract + * + * A rule is complete and structurally valid when construction returns. The + * constructor rejects unsupported cells, negative exactness, empty or mismatched + * storage, non-finite coordinates or weights, points outside the declared + * reference cell, and weights whose zeroth moment does not equal the reference + * measure in both compensated arithmetic and ordinary stored-order double + * accumulation. A condition estimate based on the absolute weight sum rejects + * signed rules whose cancellation is too sensitive for double precision. + * Negative individual weights remain valid because some quadrature families + * require them. + * + * polynomial_exactness() is declared by the concrete construction algorithm. + * Structural validation verifies metadata, containment, finiteness, and the + * zeroth moment; it does not prove higher-order polynomial moments. Every + * concrete rule provider is responsible for making its declaration a guarantee + * by establishing the advertised exactness with analytic moment tests. + * + * Points use one zero-initialized three-component representation. Constructors + * that receive fewer than three coordinates zero-fill the omitted components. + * Only the first dimension() components are active, and every inactive + * component is zero within the construction tolerance. The supported canonical + * domains are: + * + * | Cell family | Canonical reference domain | Measure | + * | ----------- | -------------------------- | ------- | + * | Point | @f$(0,0,0)@f$ | @f$1@f$ | + * | Line | @f$[-1,1]@f$ | @f$2@f$ | + * | Triangle | @f$\xi,\eta\geq0;\ \xi+\eta\leq1@f$ | @f$1/2@f$ | + * | Quad | @f$[-1,1]^2@f$ | @f$4@f$ | + * | Tetra | @f$\xi,\eta,\zeta\geq0;\ \xi+\eta+\zeta\leq1@f$ | @f$1/6@f$ | + * | Hex | @f$[-1,1]^3@f$ | @f$8@f$ | + * | Wedge | reference triangle @f$\times[-1,1]@f$ | @f$1@f$ | + * + * Pyramid, Polygon, Polyhedron, and unknown cell families are intentionally + * unsupported. + * + * ## Ownership and lifetime + * + * Rule objects are non-copyable and non-movable. They are intended to be built + * once, retained through a const owning handle, and shared across integrations. + * References returned by points() and weights() remain valid for the lifetime + * of the rule. Consumers should retain that rule rather than copying its point + * count or samples into parallel authoritative state. + */ + +#include "FE/Common/Types.h" +#include "Math/Vector.h" + +#include +#include + +namespace svmp::FE::quadrature { + +/** @addtogroup FE_Quadrature + * @{ + */ + +/** + * @brief Zero-initialized three-component reference coordinate. + * + * Only the first QuadratureRule::dimension() components are active. Remaining + * components are zero, giving point, line, surface, and volume rules a uniform + * representation. + */ +class QuadPoint { +public: + /** @brief Fixed-size FE vector used as the coordinate storage type. */ + using CoordinateVector = math::Vector; + + /** @brief Construct the origin, with all three components initialized to zero. */ + QuadPoint() : coordinates_(CoordinateVector::Zero()) {} + + /** + * @brief Construct a line coordinate and zero-fill the remaining components. + * @param x First reference coordinate. + */ + explicit QuadPoint(double x) + : coordinates_(x, 0.0, 0.0) + { + } + + /** + * @brief Construct a surface coordinate and zero-fill the third component. + * @param x First reference coordinate. + * @param y Second reference coordinate. + */ + QuadPoint(double x, double y) + : coordinates_(x, y, 0.0) + { + } + + /** + * @brief Construct a three-component coordinate. + * @param x First reference coordinate. + * @param y Second reference coordinate. + * @param z Third reference coordinate. + */ + QuadPoint(double x, double y, double z) + : coordinates_(x, y, z) + { + } + + /** + * @brief Construct from an initialized fixed-size FE vector. + * @param coordinates Three reference-coordinate components. + */ + explicit QuadPoint(const CoordinateVector& coordinates) + : coordinates_(coordinates) + { + } + + /** + * @brief Return the zero coordinate. + * @return Three-component coordinate initialized to the origin. + */ + static QuadPoint Zero() { return {}; } + + /** + * @brief Access one coordinate component without bounds checking. + * @param component Component index in the half-open range `[0, 3)`. + * @return Mutable component reference for rule construction. + */ + double& operator[](std::size_t component) noexcept + { + return coordinates_[static_cast(component)]; + } + + /** + * @brief Access one coordinate component without bounds checking. + * @param component Component index in the half-open range `[0, 3)`. + * @return Immutable component reference. + */ + const double& operator[](std::size_t component) const noexcept + { + return coordinates_[static_cast(component)]; + } + + /** + * @brief Return mutable fixed-size FE coordinate storage for rule generation. + * @return Mutable three-component coordinate storage. + */ + CoordinateVector& coordinates() noexcept { return coordinates_; } + + /** + * @brief Return immutable fixed-size FE coordinate storage. + * @return Immutable three-component coordinate storage. + */ + const CoordinateVector& coordinates() const noexcept { return coordinates_; } + +private: + CoordinateVector coordinates_; +}; + +/** + * @brief Immutable consumer interface for a quadrature rule on a reference cell. + * + * Concrete rule providers initialize the base with one complete RuleData + * payload. The payload is validated before construction returns, and the object + * exposes no mutation or assignment path afterward. General solver consumers + * use only the public const query interface. + */ +class QuadratureRule { +public: + /** @brief Destroy a quadrature rule through the abstract interface. */ + virtual ~QuadratureRule() = 0; + + /** @brief Rule objects cannot be copied through the abstract interface. */ + QuadratureRule(const QuadratureRule&) = delete; + + /** @brief Rule objects cannot be moved through the abstract interface. */ + QuadratureRule(QuadratureRule&&) = delete; + + /** @brief Rule objects cannot be replaced through base assignment. */ + QuadratureRule& operator=(const QuadratureRule&) = delete; + + /** @brief Rule objects cannot be replaced through base move assignment. */ + QuadratureRule& operator=(QuadratureRule&&) = delete; + + /** + * @brief Return the number of ordered point/weight pairs. + * @return Quadrature sample count. + */ + std::size_t num_points() const noexcept { return points_.size(); } + + /** + * @brief Return the total-degree polynomial exactness declared by the rule. + * + * Structural validation does not independently prove this degree. + * Concrete rule providers establish it through analytic moment tests. + * + * @return Declared total degree that a conforming rule integrates exactly. + */ + int polynomial_exactness() const noexcept { return polynomial_exactness_; } + + /** + * @brief Return the intrinsic dimension of the reference integration domain. + * @return Active coordinate count, from zero for Point through three for volume cells. + */ + int dimension() const noexcept { return dimension_; } + + /** + * @brief Return the canonical reference-cell family. + * @return Reference topology integrated by this rule. + */ + svmp::CellFamily cell_family() const noexcept { return cell_family_; } + + /** + * @brief Return one reference coordinate without bounds checking. + * @param i Point index in the half-open range `[0, num_points())`. + * @return Immutable reference to the indexed quadrature point, valid for + * the lifetime of this rule. + * @pre @p i is less than num_points(). + */ + const QuadPoint& point(std::size_t i) const noexcept { return points_[i]; } + + /** + * @brief Return one reference weight without bounds checking. + * @param i Weight index in the half-open range `[0, num_points())`. + * @return Weight paired with point(@p i). + * @pre @p i is less than num_points(). + */ + double weight(std::size_t i) const noexcept { return weights_[i]; } + + /** + * @brief Return all reference coordinates in integration order. + * @return Immutable point storage, valid for the lifetime of this rule. + */ + const std::vector& points() const noexcept { return points_; } + + /** + * @brief Return all reference weights in point order. + * @return Immutable weight storage, valid for the lifetime of this rule. + */ + const std::vector& weights() const noexcept { return weights_; } + + /** + * @brief Return the default tolerance used during structural validation. + * + * Coordinate bounds and inactive components use this as an absolute + * tolerance. Compensated and ordinary stored-order zeroth-moment checks + * scale it by the larger of one and the reference-cell measure. For signed + * rules, a separate cancellation-sensitivity bound based on double epsilon + * and the absolute weight sum always uses this default tolerance, so callers + * cannot waive the minimum double-precision stability requirement. + * + * @return Default absolute validation tolerance. + */ + static constexpr double default_validation_tolerance() noexcept { return 1.0e-12; } + + /** + * @brief Recheck the rule's structural invariants using a caller tolerance. + * + * This diagnostic verifies metadata, storage, finite values, point + * containment, and the zeroth moment. It does not prove + * polynomial_exactness(). Construction already performs this check with + * default_validation_tolerance(). + * + * @param tol Finite, non-negative absolute coordinate and scaled measure tolerance. + * @return True when the stored rule satisfies every structural invariant. + */ + bool is_structurally_valid( + double tol = default_validation_tolerance()) const noexcept; + + /** + * @brief Return the measure of the canonical reference cell. + * @return Reference length, area, volume, or unit point measure. + */ + double reference_measure() const noexcept { return reference_measure_; } + +protected: + /** + * @brief Complete construction payload used by concrete rule providers. + * + * Integration consumers do not construct this payload. A provider computes + * all three fields before invoking the protected base constructor and must + * substantiate its declared exactness with analytic moment tests. + */ + struct RuleData { + int polynomial_exactness{-1}; ///< Declared total-degree polynomial exactness. + std::vector points; ///< Ordered canonical reference coordinates. + std::vector weights; ///< Weights paired with points in the same order. + }; + + /** + * @brief Construct and validate one complete immutable rule. + * + * Dimension and reference measure are derived from @p family; callers cannot + * supply redundant topology metadata. + * + * @param family Supported canonical reference-cell family. + * @param data Complete exactness, point, and weight payload. + * @throws InvalidArgumentException If the family is unsupported, exactness + * is negative, storage is empty or mismatched, a value is non-finite, a point + * is outside the reference cell, the weights do not reproduce its measure in + * compensated and stored-order arithmetic, or their cancellation is too + * ill-conditioned for stable double-precision integration. + */ + explicit QuadratureRule(svmp::CellFamily family, RuleData data); + +private: + /** @brief Fully checked state used by the delegating constructor. */ + struct ValidatedState { + svmp::CellFamily cell_family; + int dimension; + int polynomial_exactness; + double reference_measure; + std::vector points; + std::vector weights; + }; + + /** @brief Validate a local payload before initializing immutable members. */ + static ValidatedState validate(svmp::CellFamily family, RuleData data); + + /** @brief Initialize members from state already checked by validate(). */ + explicit QuadratureRule(ValidatedState state); + + const svmp::CellFamily cell_family_; ///< Canonical reference topology. + const int dimension_; ///< Number of active coordinate components. + const int polynomial_exactness_; ///< Exactness declared by the concrete generator. + const double reference_measure_; ///< Canonical cell measure. + const std::vector points_; ///< Ordered immutable reference coordinates. + const std::vector weights_; ///< Immutable weights paired with points_. +}; + +/** @} */ + +} // namespace svmp::FE::quadrature + +#endif // SVMP_FE_QUADRATURE_RULE_H From fbf363094c84faa403fc5d905cd9bb7317f1ae1f Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Wed, 22 Jul 2026 10:55:17 -0700 Subject: [PATCH 03/21] Establish Quadrature Phase 01 baselines --- .../FE/Quadrature/test_QuadratureRules.cpp | 1141 +++++++++++++++++ 1 file changed, 1141 insertions(+) create mode 100644 tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp new file mode 100644 index 000000000..7e8341ec9 --- /dev/null +++ b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp @@ -0,0 +1,1141 @@ +/** + * @file test_QuadratureRules.cpp + * @brief Baseline solver data and core quadrature rule contract tests. + */ + +#include + +#include "FE/Common/FEException.h" +#include "FE/Quadrature/QuadratureRule.h" +#include "nn.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace svmp::FE; +using namespace svmp::FE::quadrature; + +namespace { + +class RuleProbe final : public QuadratureRule { +public: + RuleProbe(svmp::CellFamily family, + int polynomial_exactness, + std::vector points, + std::vector weights) + : QuadratureRule( + family, + RuleData{ + polynomial_exactness, + std::move(points), + std::move(weights)}) + { + } +}; + +constexpr double kTol = 1.0e-12; +constexpr double kExactnessTol = 3.0e-12; + +using ExpectedPoint = std::array; + +struct ExpectedSample { + ExpectedPoint point; + double weight; +}; + +using ExpectedSamples = std::vector; + +struct SolverRuleCase { + const char* name; + consts::ElementType solver_type; + ElementType fe_type; + int dimension; + int requested_exactness; + int advertised_exactness; + double reference_measure; + ExpectedSamples samples; +}; + +template +void expect_invalid_argument_with_message( + Function&& function, + const std::string& expected_message) +{ + try { + std::forward(function)(); + FAIL() << "Expected InvalidArgumentException containing: " + << expected_message; + } catch (const InvalidArgumentException& exception) { + const std::string actual_message = exception.what(); + EXPECT_NE(actual_message.find(expected_message), std::string::npos) + << "actual message: " << actual_message; + } catch (const std::exception& exception) { + FAIL() << "Expected InvalidArgumentException, received: " + << exception.what(); + } catch (...) { + FAIL() << "Expected InvalidArgumentException, received an unknown exception"; + } +} + +ExpectedSamples point_samples() +{ + return {{{0.0, 0.0, 0.0}, 1.0}}; +} + +ExpectedSamples line2_samples() +{ + const double a = 1.0 / std::sqrt(3.0); + return { + {{-a, 0.0, 0.0}, 1.0}, + {{ a, 0.0, 0.0}, 1.0}, + }; +} + +ExpectedSamples line3_samples() +{ + const double a = std::sqrt(0.6); + return { + {{-a, 0.0, 0.0}, 5.0 / 9.0}, + {{ a, 0.0, 0.0}, 5.0 / 9.0}, + {{0.0, 0.0, 0.0}, 8.0 / 9.0}, + }; +} + +ExpectedSamples triangle3_samples() +{ + const double s = 2.0 / 3.0; + const double t = 1.0 / 6.0; + return { + {{t, t, 0.0}, 1.0 / 6.0}, + {{s, t, 0.0}, 1.0 / 6.0}, + {{t, s, 0.0}, 1.0 / 6.0}, + }; +} + +ExpectedSamples triangle6_samples() +{ + const double centroid = 0.333333333333333; + const double a1 = 0.797426985353087; + const double b1 = 0.101286507323456; + const double a2 = 0.059715871789770; + const double b2 = 0.470142064105115; + const double w0 = 0.225000000000000 * 0.5; + const double w1 = 0.125939180544827 * 0.5; + const double w2 = 0.132394152788506 * 0.5; + return { + {{centroid, centroid, 0.0}, w0}, + {{a1, b1, 0.0}, w1}, + {{b1, a1, 0.0}, w1}, + {{b1, b1, 0.0}, w1}, + {{a2, b2, 0.0}, w2}, + {{b2, a2, 0.0}, w2}, + {{b2, b2, 0.0}, w2}, + }; +} + +ExpectedSamples quad4_samples() +{ + const double a = 1.0 / std::sqrt(3.0); + return { + {{-a, -a, 0.0}, 1.0}, + {{ a, -a, 0.0}, 1.0}, + {{ a, a, 0.0}, 1.0}, + {{-a, a, 0.0}, 1.0}, + }; +} + +ExpectedSamples quad9_samples(double a = std::sqrt(0.6)) +{ + const double corner_weight = 25.0 / 81.0; + const double edge_weight = 40.0 / 81.0; + return { + {{-a, -a, 0.0}, corner_weight}, + {{ a, -a, 0.0}, corner_weight}, + {{ a, a, 0.0}, corner_weight}, + {{-a, a, 0.0}, corner_weight}, + {{0.0, -a, 0.0}, edge_weight}, + {{ a, 0.0, 0.0}, edge_weight}, + {{0.0, a, 0.0}, edge_weight}, + {{-a, 0.0, 0.0}, edge_weight}, + {{0.0, 0.0, 0.0}, 64.0 / 81.0}, + }; +} + +ExpectedSamples tetra4_samples() +{ + const double a = (5.0 + 3.0 * std::sqrt(5.0)) / 20.0; + const double b = (1.0 - a) / 3.0; + return { + {{a, b, b}, 1.0 / 24.0}, + {{b, a, b}, 1.0 / 24.0}, + {{b, b, a}, 1.0 / 24.0}, + {{b, b, b}, 1.0 / 24.0}, + }; +} + +ExpectedSamples tetra10_samples() +{ + const double one_third = 0.3333333333333330; + const double a1 = 0.0909090909090910; + const double b1 = 0.7272727272727270; + const double a2 = 0.0665501535736640; + const double b2 = 0.4334498464263360; + const double w0 = 0.0302836780970890; + const double w1 = 0.0060267857142860; + const double w2 = 0.0116452490860290; + const double w3 = 0.0109491415613860; + return { + {{0.25, 0.25, 0.25}, w0}, + {{0.0, one_third, one_third}, w1}, + {{one_third, 0.0, one_third}, w1}, + {{one_third, one_third, 0.0}, w1}, + {{one_third, one_third, one_third}, w1}, + {{b1, a1, a1}, w2}, + {{a1, b1, a1}, w2}, + {{a1, a1, b1}, w2}, + {{a1, a1, a1}, w2}, + {{a2, a2, b2}, w3}, + {{a2, b2, a2}, w3}, + {{a2, b2, b2}, w3}, + {{b2, b2, a2}, w3}, + {{b2, a2, b2}, w3}, + {{b2, a2, a2}, w3}, + }; +} + +ExpectedSamples hex8_samples() +{ + const double a = 1.0 / std::sqrt(3.0); + return { + {{-a, -a, -a}, 1.0}, + {{ a, -a, -a}, 1.0}, + {{ a, a, -a}, 1.0}, + {{-a, a, -a}, 1.0}, + {{-a, -a, a}, 1.0}, + {{ a, -a, a}, 1.0}, + {{ a, a, a}, 1.0}, + {{-a, a, a}, 1.0}, + }; +} + +ExpectedSamples hex27_samples() +{ + const double a = std::sqrt(0.6); + const double w0 = 125.0 / 729.0; + const double w1 = 200.0 / 729.0; + const double w2 = 320.0 / 729.0; + return { + {{-a, -a, -a}, w0}, {{ a, -a, -a}, w0}, + {{ a, a, -a}, w0}, {{-a, a, -a}, w0}, + {{-a, -a, a}, w0}, {{ a, -a, a}, w0}, + {{ a, a, a}, w0}, {{-a, a, a}, w0}, + + {{0.0, -a, -a}, w1}, {{ a, 0.0, -a}, w1}, + {{0.0, a, -a}, w1}, {{-a, 0.0, -a}, w1}, + {{0.0, -a, a}, w1}, {{ a, 0.0, a}, w1}, + {{0.0, a, a}, w1}, {{-a, 0.0, a}, w1}, + {{-a, -a, 0.0}, w1}, {{ a, -a, 0.0}, w1}, + {{ a, a, 0.0}, w1}, {{-a, a, 0.0}, w1}, + + {{-a, 0.0, 0.0}, w2}, {{ a, 0.0, 0.0}, w2}, + {{0.0, -a, 0.0}, w2}, {{0.0, a, 0.0}, w2}, + {{0.0, 0.0, -a}, w2}, {{0.0, 0.0, a}, w2}, + {{0.0, 0.0, 0.0}, 512.0 / 729.0}, + }; +} + +ExpectedSamples wedge6_samples() +{ + const double s = 2.0 / 3.0; + const double t = 1.0 / 6.0; + const double z = 1.0 / std::sqrt(3.0); + return { + {{s, t, -z}, 1.0 / 6.0}, + {{t, s, -z}, 1.0 / 6.0}, + {{t, t, -z}, 1.0 / 6.0}, + {{s, t, z}, 1.0 / 6.0}, + {{t, s, z}, 1.0 / 6.0}, + {{t, t, z}, 1.0 / 6.0}, + }; +} + +const std::vector& standard_solver_cases() +{ + static const std::vector cases = { + {"PNT", consts::ElementType::PNT, ElementType::Point1, 0, 1, 1, 1.0, point_samples()}, + {"LIN1", consts::ElementType::LIN1, ElementType::Line2, 1, 2, 3, 2.0, line2_samples()}, + {"LIN2", consts::ElementType::LIN2, ElementType::Line3, 1, 4, 5, 2.0, line3_samples()}, + {"TRI3", consts::ElementType::TRI3, ElementType::Triangle3, 2, 2, 2, 0.5, triangle3_samples()}, + {"TRI6", consts::ElementType::TRI6, ElementType::Triangle6, 2, 5, 5, 0.5, triangle6_samples()}, + {"QUD4", consts::ElementType::QUD4, ElementType::Quad4, 2, 2, 3, 4.0, quad4_samples()}, + {"QUD8", consts::ElementType::QUD8, ElementType::Quad8, 2, 4, 5, 4.0, quad9_samples()}, + {"QUD9", consts::ElementType::QUD9, ElementType::Quad9, 2, 4, 5, 4.0, quad9_samples()}, + {"TET4", consts::ElementType::TET4, ElementType::Tetra4, 3, 2, 2, 1.0 / 6.0, tetra4_samples()}, + {"TET10", consts::ElementType::TET10, ElementType::Tetra10, 3, 5, 5, 1.0 / 6.0, tetra10_samples()}, + {"HEX8", consts::ElementType::HEX8, ElementType::Hex8, 3, 2, 3, 8.0, hex8_samples()}, + {"HEX20", consts::ElementType::HEX20, ElementType::Hex20, 3, 4, 5, 8.0, hex27_samples()}, + {"HEX27", consts::ElementType::HEX27, ElementType::Hex27, 3, 4, 5, 8.0, hex27_samples()}, + {"WDG", consts::ElementType::WDG, ElementType::Wedge6, 3, 2, 2, 1.0, wedge6_samples()}, + }; + return cases; +} + +std::vector rule_points(const ExpectedSamples& samples) +{ + std::vector points; + points.reserve(samples.size()); + for (const auto& sample : samples) { + points.push_back({ + sample.point[0], + sample.point[1], + sample.point[2], + }); + } + return points; +} + +std::vector rule_weights(const ExpectedSamples& samples) +{ + std::vector weights; + weights.reserve(samples.size()); + for (const auto& sample : samples) { + weights.push_back(sample.weight); + } + return weights; +} + +std::vector rule_points(const Array& points, int dimension) +{ + std::vector result; + result.reserve(static_cast(points.ncols())); + for (int q = 0; q < points.ncols(); ++q) { + QuadPoint point = QuadPoint::Zero(); + for (int component = 0; component < dimension; ++component) { + point[static_cast(component)] = points(component, q); + } + result.push_back(point); + } + return result; +} + +std::vector rule_weights(const Vector& weights) +{ + std::vector result; + result.reserve(static_cast(weights.size())); + for (int q = 0; q < weights.size(); ++q) { + result.push_back(weights[q]); + } + return result; +} + +double factorial(int n) +{ + double value = 1.0; + for (int factor = 2; factor <= n; ++factor) { + value *= static_cast(factor); + } + return value; +} + +double interval_moment(int power) +{ + return power % 2 == 0 + ? 2.0 / static_cast(power + 1) + : 0.0; +} + +double simplex_moment(int px, int py, int pz, int dimension) +{ + const int total = px + + (dimension >= 2 ? py : 0) + + (dimension >= 3 ? pz : 0); + double numerator = factorial(px); + if (dimension >= 2) { + numerator *= factorial(py); + } + if (dimension >= 3) { + numerator *= factorial(pz); + } + return numerator / factorial(total + dimension); +} + +double exact_reference_moment( + svmp::CellFamily family, + int px, + int py, + int pz) +{ + switch (family) { + case svmp::CellFamily::Point: + return px == 0 && py == 0 && pz == 0 ? 1.0 : 0.0; + case svmp::CellFamily::Line: + return interval_moment(px); + case svmp::CellFamily::Triangle: + return simplex_moment(px, py, 0, 2); + case svmp::CellFamily::Quad: + return interval_moment(px) * interval_moment(py); + case svmp::CellFamily::Tetra: + return simplex_moment(px, py, pz, 3); + case svmp::CellFamily::Hex: + return interval_moment(px) * + interval_moment(py) * + interval_moment(pz); + case svmp::CellFamily::Wedge: + return simplex_moment(px, py, 0, 2) * interval_moment(pz); + default: + svmp::raise( + "Unsupported reference cell in exactness test"); + } +} + +double integer_power(double base, int exponent) +{ + double value = 1.0; + for (int factor = 0; factor < exponent; ++factor) { + value *= base; + } + return value; +} + +double integrate_monomial( + const QuadratureRule& rule, + int px, + int py, + int pz) +{ + double result = 0.0; + for (std::size_t q = 0; q < rule.num_points(); ++q) { + const auto point = rule.point(q); + result += rule.weight(q) * + integer_power(point[0], px) * + integer_power(point[1], py) * + integer_power(point[2], pz); + } + return result; +} + +void expect_total_degree_exact( + const QuadratureRule& rule, + int degree, + double tolerance = kExactnessTol) +{ + const int max_px = rule.dimension() >= 1 ? degree : 0; + for (int px = 0; px <= max_px; ++px) { + const int max_py = rule.dimension() >= 2 ? degree - px : 0; + for (int py = 0; py <= max_py; ++py) { + const int max_pz = + rule.dimension() >= 3 ? degree - px - py : 0; + for (int pz = 0; pz <= max_pz; ++pz) { + EXPECT_NEAR( + integrate_monomial(rule, px, py, pz), + exact_reference_moment( + rule.cell_family(), px, py, pz), + tolerance) + << "powers=(" << px << ',' << py << ',' << pz << ')'; + } + } + } +} + +void expect_samples_in_order( + const QuadratureRule& rule, + const ExpectedSamples& expected, + double tolerance = kTol) +{ + ASSERT_EQ(rule.num_points(), expected.size()); + for (std::size_t q = 0; q < expected.size(); ++q) { + EXPECT_NEAR(rule.weight(q), expected[q].weight, tolerance) + << "sample=" << q; + const auto point = rule.point(q); + for (std::size_t component = 0; component < 3u; ++component) { + EXPECT_NEAR( + point[component], + expected[q].point[component], + tolerance) + << "sample=" << q << " component=" << component; + } + } +} + +RuleProbe make_rule(const SolverRuleCase& c) +{ + return RuleProbe( + to_mesh_family(c.fe_type), + c.advertised_exactness, + rule_points(c.samples), + rule_weights(c.samples)); +} + +void expect_legacy_arrays_in_order(const Vector& weights, + const Array& points, + int dimension, + const ExpectedSamples& expected, + double tol = kTol) +{ + ASSERT_EQ(weights.size(), static_cast(expected.size())); + ASSERT_EQ(points.nrows(), dimension); + if (dimension > 0) { + ASSERT_EQ(points.ncols(), static_cast(expected.size())); + } + + for (std::size_t q = 0; q < expected.size(); ++q) { + EXPECT_NEAR(weights[static_cast(q)], expected[q].weight, tol) + << "sample=" << q; + for (int d = 0; d < dimension; ++d) { + EXPECT_NEAR(points(d, static_cast(q)), + expected[q].point[static_cast(d)], tol) + << "sample=" << q << " component=" << d; + } + } +} + +bool generic_context_supports(consts::ElementType type) +{ + switch (type) { + case consts::ElementType::LIN1: + case consts::ElementType::TRI3: + case consts::ElementType::QUD4: + case consts::ElementType::QUD9: + case consts::ElementType::TET4: + case consts::ElementType::HEX8: + case consts::ElementType::HEX20: + case consts::ElementType::HEX27: + case consts::ElementType::WDG: + return true; + default: + return false; + } +} + +bool volume_context_supports(consts::ElementType type) +{ + switch (type) { + case consts::ElementType::LIN1: + case consts::ElementType::LIN2: + case consts::ElementType::TRI3: + case consts::ElementType::TRI6: + case consts::ElementType::QUD4: + case consts::ElementType::QUD9: + case consts::ElementType::TET4: + case consts::ElementType::TET10: + case consts::ElementType::HEX8: + case consts::ElementType::HEX20: + case consts::ElementType::HEX27: + case consts::ElementType::WDG: + return true; + default: + return false; + } +} + +bool face_context_supports(consts::ElementType type) +{ + switch (type) { + case consts::ElementType::PNT: + case consts::ElementType::LIN1: + case consts::ElementType::LIN2: + case consts::ElementType::TRI3: + case consts::ElementType::TRI6: + case consts::ElementType::QUD4: + case consts::ElementType::QUD8: + case consts::ElementType::QUD9: + return true; + default: + return false; + } +} + +} // namespace + +TEST(QuadraturePhase01Baseline, StandardSelectionTableHasOrderedPointWeightData) +{ + const auto& cases = standard_solver_cases(); + ASSERT_EQ(cases.size(), 14u); + + for (const auto& c : cases) { + SCOPED_TRACE(c.name); + EXPECT_FALSE(c.samples.empty()); + EXPECT_GT(c.requested_exactness, 0); + EXPECT_GE(c.advertised_exactness, c.requested_exactness); + EXPECT_GT(c.reference_measure, 0.0); + + double weight_sum = 0.0; + for (const auto& sample : c.samples) { + weight_sum += sample.weight; + } + EXPECT_NEAR(weight_sum, c.reference_measure, kTol); + } +} + +TEST(QuadraturePhase01Baseline, CanonicalFixturesSatisfyTheRuleAndExactnessContracts) +{ + for (const auto& c : standard_solver_cases()) { + SCOPED_TRACE(c.name); + const auto rule = make_rule(c); + + EXPECT_EQ(rule.cell_family(), to_mesh_family(c.fe_type)); + EXPECT_EQ(rule.dimension(), c.dimension); + EXPECT_EQ(rule.polynomial_exactness(), c.advertised_exactness); + EXPECT_EQ(rule.num_points(), c.samples.size()); + EXPECT_DOUBLE_EQ(rule.reference_measure(), c.reference_measure); + EXPECT_TRUE(rule.is_structurally_valid()); + expect_samples_in_order(rule, c.samples); + expect_total_degree_exact(rule, c.advertised_exactness); + } +} + +TEST(QuadraturePhase01Baseline, GenericFunctionSpaceRulesMatchOrderedLegacyData) +{ + for (const auto& c : standard_solver_cases()) { + if (!generic_context_supports(c.solver_type) || + c.solver_type == consts::ElementType::QUD9) { + // QUD9 is characterized separately because the legacy generic + // table contains an out-of-domain sqrt(6) coordinate defect. + continue; + } + + SCOPED_TRACE(c.name); + Vector weights(static_cast(c.samples.size())); + Array points(c.dimension, static_cast(c.samples.size())); + nn::get_gip(c.dimension, c.solver_type, + static_cast(c.samples.size()), weights, points); + expect_legacy_arrays_in_order(weights, points, c.dimension, c.samples); + } +} + +TEST(QuadraturePhase01Baseline, GenericQud9SqrtSixCoordinatesAreAnExplicitLegacyDefect) +{ + const auto legacy_defect = quad9_samples(std::sqrt(6.0)); + Vector weights(static_cast(legacy_defect.size())); + Array points(2, static_cast(legacy_defect.size())); + + nn::get_gip( + 2, + consts::ElementType::QUD9, + static_cast(legacy_defect.size()), + weights, + points); + + expect_legacy_arrays_in_order(weights, points, 2, legacy_defect); + EXPECT_GT(std::abs(points(0, 0)), 1.0); + + double x_squared_moment = 0.0; + for (int q = 0; q < weights.size(); ++q) { + x_squared_moment += weights[q] * points(0, q) * points(0, q); + } + EXPECT_NEAR(x_squared_moment, 40.0 / 3.0, kExactnessTol); + EXPECT_NEAR( + exact_reference_moment(svmp::CellFamily::Quad, 2, 0, 0), + 4.0 / 3.0, + kTol); + + EXPECT_THROW( + (void)RuleProbe( + svmp::CellFamily::Quad, + 5, + rule_points(points, 2), + rule_weights(weights)), + InvalidArgumentException); +} + +TEST(QuadraturePhase01Baseline, VolumeMeshRulesMatchOrderedLegacyData) +{ + for (const auto& c : standard_solver_cases()) { + if (!volume_context_supports(c.solver_type)) { + continue; + } + + SCOPED_TRACE(c.name); + mshType mesh; + mesh.eType = c.solver_type; + mesh.nG = static_cast(c.samples.size()); + mesh.w = Vector(mesh.nG); + mesh.xi = Array(c.dimension, mesh.nG); + nn::get_gip(mesh); + expect_legacy_arrays_in_order(mesh.w, mesh.xi, c.dimension, c.samples); + } +} + +TEST(QuadraturePhase01Baseline, Qud8VolumeSelectionHasNoLegacyPopulationEntry) +{ + const auto& qud8_case = standard_solver_cases()[6]; + ASSERT_EQ(qud8_case.solver_type, consts::ElementType::QUD8); + + mshType mesh; + mesh.eType = qud8_case.solver_type; + mesh.nG = static_cast(qud8_case.samples.size()); + mesh.w = Vector(mesh.nG); + mesh.xi = Array(qud8_case.dimension, mesh.nG); + + EXPECT_THROW(nn::get_gip(mesh), InvalidElementException); + + const auto canonical_rule = make_rule(qud8_case); + expect_samples_in_order(canonical_rule, qud8_case.samples); + expect_total_degree_exact( + canonical_rule, + qud8_case.advertised_exactness); +} + +TEST(QuadraturePhase01Baseline, GenericFunctionSpaceMissingEntriesAreExplicit) +{ + for (const auto& c : standard_solver_cases()) { + if (generic_context_supports(c.solver_type)) { + continue; + } + + SCOPED_TRACE(c.name); + Vector weights(static_cast(c.samples.size())); + Array points(c.dimension, static_cast(c.samples.size())); + EXPECT_THROW( + nn::get_gip(c.dimension, c.solver_type, + static_cast(c.samples.size()), weights, points), + InvalidElementException); + } +} + +TEST(QuadraturePhase01Baseline, FaceRulesMatchOrderedLegacyData) +{ + for (const auto& c : standard_solver_cases()) { + if (!face_context_supports(c.solver_type)) { + continue; + } + + SCOPED_TRACE(c.name); + faceType face; + face.eType = c.solver_type; + face.nG = static_cast(c.samples.size()); + face.w = Vector(face.nG); + face.xi = Array(c.dimension, face.nG); + nn::get_gip(nullptr, face); + expect_legacy_arrays_in_order(face.w, face.xi, c.dimension, c.samples); + } +} + +TEST(QuadraturePhase01Baseline, LegacyPositionModifiersAreOnlyDegreeOneExact) +{ + mshType tetra_mesh; + tetra_mesh.eType = consts::ElementType::TET4; + tetra_mesh.nG = 4; + tetra_mesh.qmTET4 = 0.25; + tetra_mesh.w = Vector(tetra_mesh.nG); + tetra_mesh.xi = Array(3, tetra_mesh.nG); + nn::get_gip(tetra_mesh); + + const RuleProbe tetra_rule( + svmp::CellFamily::Tetra, + 1, + rule_points(tetra_mesh.xi, 3), + rule_weights(tetra_mesh.w)); + expect_total_degree_exact(tetra_rule, 1); + EXPECT_NEAR( + integrate_monomial(tetra_rule, 2, 0, 0), + 1.0 / 96.0, + kExactnessTol); + EXPECT_GT( + std::abs( + integrate_monomial(tetra_rule, 2, 0, 0) - + exact_reference_moment(svmp::CellFamily::Tetra, 2, 0, 0)), + 1.0e-3); + + faceType triangle_face; + triangle_face.eType = consts::ElementType::TRI3; + triangle_face.nG = 3; + triangle_face.qmTRI3 = 1.0 / 3.0; + triangle_face.w = Vector(triangle_face.nG); + triangle_face.xi = Array(2, triangle_face.nG); + nn::get_gip(nullptr, triangle_face); + + const RuleProbe triangle_rule( + svmp::CellFamily::Triangle, + 1, + rule_points(triangle_face.xi, 2), + rule_weights(triangle_face.w)); + expect_total_degree_exact(triangle_rule, 1); + EXPECT_NEAR( + integrate_monomial(triangle_rule, 2, 0, 0), + 1.0 / 18.0, + kExactnessTol); + EXPECT_GT( + std::abs( + integrate_monomial(triangle_rule, 2, 0, 0) - + exact_reference_moment(svmp::CellFamily::Triangle, 2, 0, 0)), + 1.0e-3); +} + +TEST(QuadPointContract, DefaultAndPartialConstructionZeroFillCoordinates) +{ + static_assert(!std::is_convertible_v); + + const QuadPoint origin; + for (std::size_t component = 0; component < 3u; ++component) { + EXPECT_DOUBLE_EQ(origin[component], 0.0); + EXPECT_DOUBLE_EQ( + origin.coordinates()[static_cast(component)], 0.0); + } + + const QuadPoint line_point{0.25}; + EXPECT_DOUBLE_EQ(line_point[0], 0.25); + EXPECT_DOUBLE_EQ(line_point[1], 0.0); + EXPECT_DOUBLE_EQ(line_point[2], 0.0); + + const QuadPoint surface_point{0.25, 0.5}; + EXPECT_DOUBLE_EQ(surface_point[0], 0.25); + EXPECT_DOUBLE_EQ(surface_point[1], 0.5); + EXPECT_DOUBLE_EQ(surface_point[2], 0.0); + + QuadPoint::CoordinateVector coordinates = + QuadPoint::CoordinateVector::Zero(); + coordinates[0] = 0.2; + coordinates[1] = 0.3; + const QuadPoint vector_point{coordinates}; + EXPECT_DOUBLE_EQ(vector_point[0], 0.2); + EXPECT_DOUBLE_EQ(vector_point[1], 0.3); + EXPECT_DOUBLE_EQ(vector_point[2], 0.0); + + QuadPoint mutable_point; + mutable_point.coordinates()[2] = 0.75; + EXPECT_DOUBLE_EQ(mutable_point[2], 0.75); + + const std::vector points(2); + for (const auto& point : points) { + for (std::size_t component = 0; component < 3u; ++component) { + EXPECT_DOUBLE_EQ(point[component], 0.0); + } + } +} + +TEST(QuadratureRuleValidation, AcceptsEverySupportedReferenceCell) +{ + struct Case { + svmp::CellFamily family; + int expected_dimension; + double expected_measure; + ExpectedPoint point; + }; + + const std::vector cases = { + {svmp::CellFamily::Point, 0, 1.0, {0.0, 0.0, 0.0}}, + {svmp::CellFamily::Line, 1, 2.0, {0.0, 0.0, 0.0}}, + {svmp::CellFamily::Triangle, 2, 0.5, {0.25, 0.25, 0.0}}, + {svmp::CellFamily::Quad, 2, 4.0, {0.0, 0.0, 0.0}}, + {svmp::CellFamily::Tetra, 3, 1.0 / 6.0, {0.25, 0.25, 0.25}}, + {svmp::CellFamily::Hex, 3, 8.0, {0.0, 0.0, 0.0}}, + {svmp::CellFamily::Wedge, 3, 1.0, {0.25, 0.25, 0.0}}, + }; + + for (const auto& c : cases) { + const RuleProbe rule( + c.family, + 0, + {{c.point[0], c.point[1], c.point[2]}}, + {c.expected_measure}); + EXPECT_EQ(rule.dimension(), c.expected_dimension); + EXPECT_DOUBLE_EQ(rule.reference_measure(), c.expected_measure); + EXPECT_TRUE(rule.is_structurally_valid()); + } +} + +TEST(QuadratureRuleValidation, RejectsInvalidMetadata) +{ + expect_invalid_argument_with_message( + [] { + (void)RuleProbe( + svmp::CellFamily::Triangle, -1, {{0.0, 0.0, 0.0}}, {0.5}); + }, + "polynomial exactness must be non-negative"); + + const std::array unsupported_families = { + svmp::CellFamily::Pyramid, + svmp::CellFamily::Polygon, + svmp::CellFamily::Polyhedron, + }; + for (const auto family : unsupported_families) { + SCOPED_TRACE(static_cast(family)); + expect_invalid_argument_with_message( + [family] { + (void)RuleProbe(family, 1, {{0.0, 0.0, 0.0}}, {1.0}); + }, + "unsupported reference-cell family"); + } + + expect_invalid_argument_with_message( + [] { + (void)RuleProbe( + static_cast(255), + 1, + {{0.0, 0.0, 0.0}}, + {1.0}); + }, + "unsupported reference-cell family"); +} + +TEST(QuadratureRuleValidation, RejectsMalformedStorageAndNonfiniteValues) +{ + const double nan = std::numeric_limits::quiet_NaN(); + const double inf = std::numeric_limits::infinity(); + + expect_invalid_argument_with_message( + [] { (void)RuleProbe(svmp::CellFamily::Line, 1, {}, {}); }, + "at least one sample"); + expect_invalid_argument_with_message( + [] { + (void)RuleProbe( + svmp::CellFamily::Line, 1, {{0.0, 0.0, 0.0}}, {}); + }, + "points/weights size mismatch"); + expect_invalid_argument_with_message( + [nan] { + (void)RuleProbe( + svmp::CellFamily::Line, 1, {{nan, 0.0, 0.0}}, {2.0}); + }, + "non-finite coordinate at sample 0"); + expect_invalid_argument_with_message( + [inf] { + (void)RuleProbe( + svmp::CellFamily::Line, 1, {{0.0, 0.0, 0.0}}, {inf}); + }, + "quadrature weight must be finite at sample 0"); + expect_invalid_argument_with_message( + [nan] { + (void)RuleProbe( + svmp::CellFamily::Line, 1, {{0.0, 0.0, 0.0}}, {nan}); + }, + "quadrature weight must be finite at sample 0"); + expect_invalid_argument_with_message( + [] { + (void)RuleProbe( + svmp::CellFamily::Line, + 1, + {{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}}, + {std::numeric_limits::max(), + std::numeric_limits::max()}); + }, + "weight accumulation is not finite at sample 1"); +} + +TEST(QuadratureRuleValidation, RejectsInactiveCoordinatesAndOutOfCellPoints) +{ + expect_invalid_argument_with_message( + [] { + (void)RuleProbe( + svmp::CellFamily::Point, 0, {{1.0e-4, 0.0, 0.0}}, {1.0}); + }, + "nonzero inactive coordinate at sample 0"); + expect_invalid_argument_with_message( + [] { + (void)RuleProbe( + svmp::CellFamily::Line, 1, {{0.0, 1.0e-4, 0.0}}, {2.0}); + }, + "nonzero inactive coordinate at sample 0"); + expect_invalid_argument_with_message( + [] { + (void)RuleProbe( + svmp::CellFamily::Triangle, 1, {{0.8, 0.3, 0.0}}, {0.5}); + }, + "outside the canonical reference cell at sample 0"); + EXPECT_THROW( + (void)RuleProbe(svmp::CellFamily::Quad, 1, {{0.0, 1.1, 0.0}}, {4.0}), + InvalidArgumentException); + EXPECT_THROW( + (void)RuleProbe(svmp::CellFamily::Tetra, 1, {{0.4, 0.4, 0.3}}, {1.0 / 6.0}), + InvalidArgumentException); + EXPECT_THROW( + (void)RuleProbe(svmp::CellFamily::Hex, 1, {{0.0, 0.0, -1.1}}, {8.0}), + InvalidArgumentException); + EXPECT_THROW( + (void)RuleProbe(svmp::CellFamily::Wedge, 1, {{0.6, 0.5, 0.0}}, {1.0}), + InvalidArgumentException); +} + +TEST(QuadratureRuleValidation, EnforcesMeasureButAllowsNegativeWeights) +{ + expect_invalid_argument_with_message( + [] { + (void)RuleProbe( + svmp::CellFamily::Triangle, + 1, + {{0.25, 0.25, 0.0}}, + {1.0}); + }, + "weights do not reproduce the reference-cell measure"); + + const RuleProbe rule( + svmp::CellFamily::Triangle, + 0, + {{1.0 / 3.0, 1.0 / 3.0, 0.0}, {0.2, 0.2, 0.0}}, + {-0.25, 0.75}); + EXPECT_TRUE(rule.is_structurally_valid()); + EXPECT_TRUE(rule.is_structurally_valid(0.0)); + EXPECT_LT(rule.weight(0), 0.0); + EXPECT_DOUBLE_EQ( + integrate_monomial(rule, 0, 0, 0), + rule.reference_measure()); + EXPECT_FALSE(rule.is_structurally_valid(-1.0)); + EXPECT_FALSE(rule.is_structurally_valid( + std::numeric_limits::quiet_NaN())); + EXPECT_FALSE(rule.is_structurally_valid( + std::numeric_limits::infinity())); +} + +TEST(QuadratureRuleValidation, RejectsStoredOrderCancellation) +{ + expect_invalid_argument_with_message( + [] { + (void)RuleProbe( + svmp::CellFamily::Line, + 0, + {{-0.75, 0.0, 0.0}, + {-0.25, 0.0, 0.0}, + {0.25, 0.0, 0.0}, + {0.75, 0.0, 0.0}}, + {1.0e16, 1.0, -1.0e16, 1.0}); + }, + "stored-order double-precision weight accumulation"); +} + +TEST(QuadratureRuleValidation, RejectsIllConditionedCancellationDespiteExactOrderedSum) +{ + expect_invalid_argument_with_message( + [] { + (void)RuleProbe( + svmp::CellFamily::Line, + 0, + {{-0.75, 0.0, 0.0}, + {-0.25, 0.0, 0.0}, + {0.25, 0.0, 0.0}, + {0.75, 0.0, 0.0}}, + {1.0e16, -1.0e16, 1.0, 1.0}); + }, + "too ill-conditioned for stable double-precision integration"); +} + +TEST(QuadratureRuleValidation, RejectsIncorrectMeasureDespiteLargeCancellation) +{ + expect_invalid_argument_with_message( + [] { + (void)RuleProbe( + svmp::CellFamily::Line, + 0, + {{-0.5, 0.0, 0.0}, + {0.0, 0.0, 0.0}, + {0.5, 0.0, 0.0}}, + {1.0e20, 0.0, -1.0e20}); + }, + "weights do not reproduce the reference-cell measure"); +} + +TEST(QuadratureRuleValidation, AcceptsWellConditionedLargePositiveRule) +{ + constexpr std::size_t sample_count = 8192u; + std::vector points(sample_count); + std::vector weights( + sample_count, + 2.0 / static_cast(sample_count)); + + const RuleProbe rule( + svmp::CellFamily::Line, + 0, + std::move(points), + std::move(weights)); + EXPECT_TRUE(rule.is_structurally_valid()); + EXPECT_DOUBLE_EQ( + integrate_monomial(rule, 0, 0, 0), + rule.reference_measure()); +} + +TEST(QuadratureRuleValidation, AcceptsWellConditionedLargeSignedRule) +{ + constexpr std::size_t sample_count = 8192u; + constexpr double cancelling_weight = 1.0 / 1048576.0; + std::vector points(sample_count); + std::vector weights( + sample_count, + 2.0 / static_cast(sample_count - 2u)); + weights[0] = -cancelling_weight; + weights[1] = cancelling_weight; + + const RuleProbe rule( + svmp::CellFamily::Line, + 0, + std::move(points), + std::move(weights)); + EXPECT_TRUE(rule.is_structurally_valid()); + EXPECT_LT(rule.weight(0), 0.0); + EXPECT_NEAR( + integrate_monomial(rule, 0, 0, 0), + rule.reference_measure(), + QuadratureRule::default_validation_tolerance()); +} + +TEST(QuadratureRuleValidation, UsesDocumentedCoordinateTolerance) +{ + constexpr double offset = 0.5 * QuadratureRule::default_validation_tolerance(); + const RuleProbe rule( + svmp::CellFamily::Line, 0, {{1.0 + offset, offset, 0.0}}, {2.0}); + EXPECT_TRUE(rule.is_structurally_valid()); + EXPECT_FALSE(rule.is_structurally_valid(offset / 2.0)); +} + +TEST(QuadratureRuleValidation, UsesDocumentedWeightTolerance) +{ + constexpr double offset = + 0.5 * QuadratureRule::default_validation_tolerance(); + const RuleProbe rule( + svmp::CellFamily::Triangle, + 0, + {{0.25, 0.25, 0.0}}, + {0.5 + offset}); + EXPECT_TRUE(rule.is_structurally_valid()); + EXPECT_FALSE(rule.is_structurally_valid(offset / 2.0)); +} + +TEST(QuadratureRuleContract, PublishesOnlyACompleteImmutableQueryInterface) +{ + static_assert(std::is_abstract_v); + static_assert(!std::is_copy_constructible_v); + static_assert(!std::is_move_constructible_v); + static_assert(!std::is_copy_assignable_v); + static_assert(!std::is_move_assignable_v); + static_assert(!std::is_copy_constructible_v); + static_assert(!std::is_move_constructible_v); + static_assert(!std::is_copy_assignable_v); + static_assert(!std::is_move_assignable_v); + static_assert( + std::is_same().point(0)), + const QuadPoint&>::value, + "A quadrature point must be exposed through an immutable reference"); + static_assert( + std::is_same().points()), + const std::vector&>::value, + "Quadrature points must be exposed through an immutable view"); + static_assert( + std::is_same().weights()), + const std::vector&>::value, + "Quadrature weights must be exposed through an immutable view"); + + const double a = 1.0 / std::sqrt(3.0); + const RuleProbe rule( + svmp::CellFamily::Line, + 3, + {{-a, 0.0, 0.0}, {a, 0.0, 0.0}}, + {1.0, 1.0}); + + EXPECT_EQ(rule.cell_family(), svmp::CellFamily::Line); + EXPECT_EQ(rule.dimension(), 1); + EXPECT_EQ(rule.polynomial_exactness(), 3); + EXPECT_DOUBLE_EQ(rule.reference_measure(), 2.0); + EXPECT_TRUE(rule.is_structurally_valid()); + ASSERT_EQ(rule.num_points(), 2u); + ASSERT_EQ(rule.points().size(), 2u); + ASSERT_EQ(rule.weights().size(), 2u); + EXPECT_DOUBLE_EQ(rule.point(0)[0], -a); + EXPECT_DOUBLE_EQ(rule.point(1)[0], a); + EXPECT_DOUBLE_EQ(rule.weight(0), 1.0); + EXPECT_DOUBLE_EQ(rule.weight(1), 1.0); +} From 6110849ff07d625bdf00eca180e449181fe0ec93 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Wed, 22 Jul 2026 10:55:22 -0700 Subject: [PATCH 04/21] Run API documentation from the repository root --- Code/CMake/SimVascularExternals.cmake | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Code/CMake/SimVascularExternals.cmake b/Code/CMake/SimVascularExternals.cmake index 50772c5b9..ed9e078ec 100644 --- a/Code/CMake/SimVascularExternals.cmake +++ b/Code/CMake/SimVascularExternals.cmake @@ -9,8 +9,8 @@ if(DOXYGEN_FOUND) configure_file(${SV_SOURCE_DIR}/../Documentation/Doxyfile ${SV_BINARY_DIR}/Doxyfile @ONLY) add_custom_target(doc - ${DOXYGEN_EXECUTABLE} ${SV_BINARY_DIR}/Doxyfile - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + COMMAND "${DOXYGEN_EXECUTABLE}" "${SV_BINARY_DIR}/Doxyfile" + WORKING_DIRECTORY "${SV_SOURCE_DIR}/.." COMMENT "Generating API documentation with Doxygen" VERBATIM ) endif(DOXYGEN_FOUND) From 593f765baf517cb82ce6afa86dc243f9dfd0a834 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Wed, 22 Jul 2026 13:46:03 -0700 Subject: [PATCH 05/21] Polish quadrature validation result handling --- .../solver/FE/Quadrature/QuadratureRule.cpp | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp index c99d72f2d..cac6882c7 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp @@ -55,8 +55,16 @@ enum class ValidationFailure { }; struct ValidationResult { + static constexpr std::size_t no_sample = + std::numeric_limits::max(); + ValidationFailure failure{ValidationFailure::None}; - std::size_t sample{std::numeric_limits::max()}; + std::size_t sample{no_sample}; + + constexpr bool valid() const noexcept + { + return failure == ValidationFailure::None; + } }; constexpr ReferenceCellTraits unsupported_traits() noexcept @@ -270,7 +278,7 @@ ValidationResult validate_rule_data( for (std::size_t sample = 0; sample < points.size(); ++sample) { const auto result = validate_point(points[sample], traits, tolerance, sample); - if (result.failure != ValidationFailure::None) { + if (!result.valid()) { return result; } } @@ -281,7 +289,7 @@ ValidationResult validate_rule_data( std::string validation_failure_message(const ValidationResult& result) { const auto sample_suffix = [&result]() { - if (result.sample == std::numeric_limits::max()) { + if (result.sample == ValidationResult::no_sample) { return std::string{}; } return std::string{" at sample "} + std::to_string(result.sample); @@ -356,7 +364,7 @@ QuadratureRule::ValidatedState QuadratureRule::validate( data.points, data.weights, default_validation_tolerance()); - if (validation.failure != ValidationFailure::None) { + if (!validation.valid()) { svmp::raise( validation_failure_message(validation)); } @@ -380,7 +388,7 @@ bool QuadratureRule::is_structurally_valid(double tolerance) const noexcept points_, weights_, tolerance) - .failure == ValidationFailure::None; + .valid(); } } // namespace svmp::FE::quadrature From 444988542de6dd9aa3cacfd9ac9f3718cfa7b188 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 23 Jul 2026 09:36:02 -0700 Subject: [PATCH 06/21] Simplify reference-cell validation --- .../solver/FE/Quadrature/QuadratureRule.cpp | 87 +++++++++---------- 1 file changed, 39 insertions(+), 48 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp index cac6882c7..dca6822c7 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp @@ -14,25 +14,14 @@ #include #include #include +#include #include #include namespace svmp::FE::quadrature { namespace { -enum class ReferenceDomain { - Unsupported, - Point, - Line, - Triangle, - Quad, - Tetra, - Hex, - Wedge, -}; - struct ReferenceCellTraits { - ReferenceDomain domain; int dimension; double measure; }; @@ -67,39 +56,32 @@ struct ValidationResult { } }; -constexpr ReferenceCellTraits unsupported_traits() noexcept -{ - return { - ReferenceDomain::Unsupported, - -1, - std::numeric_limits::quiet_NaN(), - }; -} - -constexpr ReferenceCellTraits reference_cell_traits(svmp::CellFamily family) noexcept +constexpr std::optional reference_cell_traits( + svmp::CellFamily family) noexcept { switch (family) { case svmp::CellFamily::Point: - return {ReferenceDomain::Point, 0, 1.0}; + return ReferenceCellTraits{0, 1.0}; case svmp::CellFamily::Line: - return {ReferenceDomain::Line, 1, 2.0}; + return ReferenceCellTraits{1, 2.0}; case svmp::CellFamily::Triangle: - return {ReferenceDomain::Triangle, 2, 0.5}; + return ReferenceCellTraits{2, 0.5}; case svmp::CellFamily::Quad: - return {ReferenceDomain::Quad, 2, 4.0}; + return ReferenceCellTraits{2, 4.0}; case svmp::CellFamily::Tetra: - return {ReferenceDomain::Tetra, 3, 1.0 / 6.0}; + return ReferenceCellTraits{3, 1.0 / 6.0}; case svmp::CellFamily::Hex: - return {ReferenceDomain::Hex, 3, 8.0}; + return ReferenceCellTraits{3, 8.0}; case svmp::CellFamily::Wedge: - return {ReferenceDomain::Wedge, 3, 1.0}; + return ReferenceCellTraits{3, 1.0}; default: - return unsupported_traits(); + return std::nullopt; } } ValidationResult validate_point( const QuadPoint& point, + svmp::CellFamily family, const ReferenceCellTraits& traits, double tolerance, std::size_t sample) noexcept @@ -123,38 +105,38 @@ ValidationResult validate_point( const double z = point[2]; bool is_contained = false; - switch (traits.domain) { - case ReferenceDomain::Point: + switch (family) { + case svmp::CellFamily::Point: is_contained = true; break; - case ReferenceDomain::Line: + case svmp::CellFamily::Line: is_contained = in_interval(x, -1.0, 1.0); break; - case ReferenceDomain::Triangle: + case svmp::CellFamily::Triangle: is_contained = x >= -tolerance && y >= -tolerance && x + y <= 1.0 + tolerance; break; - case ReferenceDomain::Quad: + case svmp::CellFamily::Quad: is_contained = in_interval(x, -1.0, 1.0) && in_interval(y, -1.0, 1.0); break; - case ReferenceDomain::Tetra: + case svmp::CellFamily::Tetra: is_contained = x >= -tolerance && y >= -tolerance && z >= -tolerance && x + y + z <= 1.0 + tolerance; break; - case ReferenceDomain::Hex: + case svmp::CellFamily::Hex: is_contained = in_interval(x, -1.0, 1.0) && in_interval(y, -1.0, 1.0) && in_interval(z, -1.0, 1.0); break; - case ReferenceDomain::Wedge: + case svmp::CellFamily::Wedge: is_contained = x >= -tolerance && y >= -tolerance && x + y <= 1.0 + tolerance && in_interval(z, -1.0, 1.0); break; - case ReferenceDomain::Unsupported: - break; + default: + return {ValidationFailure::UnsupportedCellFamily}; } if (!is_contained) { @@ -253,15 +235,12 @@ ValidationResult validate_weights( ValidationResult validate_rule_data( svmp::CellFamily family, + const ReferenceCellTraits& traits, int polynomial_exactness, const std::vector& points, const std::vector& weights, double tolerance) noexcept { - const auto traits = reference_cell_traits(family); - if (traits.domain == ReferenceDomain::Unsupported) { - return {ValidationFailure::UnsupportedCellFamily}; - } if (polynomial_exactness < 0) { return {ValidationFailure::NegativePolynomialExactness}; } @@ -277,7 +256,7 @@ ValidationResult validate_rule_data( for (std::size_t sample = 0; sample < points.size(); ++sample) { const auto result = - validate_point(points[sample], traits, tolerance, sample); + validate_point(points[sample], family, traits, tolerance, sample); if (!result.valid()) { return result; } @@ -358,8 +337,16 @@ QuadratureRule::ValidatedState QuadratureRule::validate( svmp::CellFamily family, RuleData data) { + const auto traits = reference_cell_traits(family); + if (!traits) { + svmp::raise( + validation_failure_message( + {ValidationFailure::UnsupportedCellFamily})); + } + const auto validation = validate_rule_data( family, + *traits, data.polynomial_exactness, data.points, data.weights, @@ -369,12 +356,11 @@ QuadratureRule::ValidatedState QuadratureRule::validate( validation_failure_message(validation)); } - const auto traits = reference_cell_traits(family); return { family, - traits.dimension, + traits->dimension, data.polynomial_exactness, - traits.measure, + traits->measure, std::move(data.points), std::move(data.weights), }; @@ -382,8 +368,13 @@ QuadratureRule::ValidatedState QuadratureRule::validate( bool QuadratureRule::is_structurally_valid(double tolerance) const noexcept { + const auto traits = reference_cell_traits(cell_family_); + if (!traits) { + return false; + } return validate_rule_data( cell_family_, + *traits, polynomial_exactness_, points_, weights_, From 454c7dd260e303a3f1a1bad7d09662b1d6d83118 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 23 Jul 2026 09:55:27 -0700 Subject: [PATCH 07/21] Use FE vectors for quadrature points --- .../solver/FE/Quadrature/QuadratureRule.h | 107 ++---------------- .../FE/Quadrature/test_QuadratureRules.cpp | 33 ++---- 2 files changed, 22 insertions(+), 118 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.h b/Code/Source/solver/FE/Quadrature/QuadratureRule.h index 7c2c15441..8bbb12dda 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.h +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.h @@ -75,11 +75,11 @@ * concrete rule provider is responsible for making its declaration a guarantee * by establishing the advertised exactness with analytic moment tests. * - * Points use one zero-initialized three-component representation. Constructors - * that receive fewer than three coordinates zero-fill the omitted components. - * Only the first dimension() components are active, and every inactive - * component is zero within the construction tolerance. The supported canonical - * domains are: + * Points use one fixed-size three-component representation. Providers initialize + * all three coordinates explicitly because the Eigen-backed vector is not + * zero-initialized by default. Only the first dimension() components are active, + * and every inactive component is zero within the construction tolerance. The + * supported canonical domains are: * * | Cell family | Canonical reference domain | Measure | * | ----------- | -------------------------- | ------- | @@ -116,100 +116,13 @@ namespace svmp::FE::quadrature { */ /** - * @brief Zero-initialized three-component reference coordinate. + * @brief Three-component coordinate used for every reference quadrature point. * - * Only the first QuadratureRule::dimension() components are active. Remaining - * components are zero, giving point, line, surface, and volume rules a uniform - * representation. + * Only the first QuadratureRule::dimension() components are active. Providers + * explicitly zero remaining components, giving point, line, surface, and volume + * rules a uniform representation directly compatible with FE math consumers. */ -class QuadPoint { -public: - /** @brief Fixed-size FE vector used as the coordinate storage type. */ - using CoordinateVector = math::Vector; - - /** @brief Construct the origin, with all three components initialized to zero. */ - QuadPoint() : coordinates_(CoordinateVector::Zero()) {} - - /** - * @brief Construct a line coordinate and zero-fill the remaining components. - * @param x First reference coordinate. - */ - explicit QuadPoint(double x) - : coordinates_(x, 0.0, 0.0) - { - } - - /** - * @brief Construct a surface coordinate and zero-fill the third component. - * @param x First reference coordinate. - * @param y Second reference coordinate. - */ - QuadPoint(double x, double y) - : coordinates_(x, y, 0.0) - { - } - - /** - * @brief Construct a three-component coordinate. - * @param x First reference coordinate. - * @param y Second reference coordinate. - * @param z Third reference coordinate. - */ - QuadPoint(double x, double y, double z) - : coordinates_(x, y, z) - { - } - - /** - * @brief Construct from an initialized fixed-size FE vector. - * @param coordinates Three reference-coordinate components. - */ - explicit QuadPoint(const CoordinateVector& coordinates) - : coordinates_(coordinates) - { - } - - /** - * @brief Return the zero coordinate. - * @return Three-component coordinate initialized to the origin. - */ - static QuadPoint Zero() { return {}; } - - /** - * @brief Access one coordinate component without bounds checking. - * @param component Component index in the half-open range `[0, 3)`. - * @return Mutable component reference for rule construction. - */ - double& operator[](std::size_t component) noexcept - { - return coordinates_[static_cast(component)]; - } - - /** - * @brief Access one coordinate component without bounds checking. - * @param component Component index in the half-open range `[0, 3)`. - * @return Immutable component reference. - */ - const double& operator[](std::size_t component) const noexcept - { - return coordinates_[static_cast(component)]; - } - - /** - * @brief Return mutable fixed-size FE coordinate storage for rule generation. - * @return Mutable three-component coordinate storage. - */ - CoordinateVector& coordinates() noexcept { return coordinates_; } - - /** - * @brief Return immutable fixed-size FE coordinate storage. - * @return Immutable three-component coordinate storage. - */ - const CoordinateVector& coordinates() const noexcept { return coordinates_; } - -private: - CoordinateVector coordinates_; -}; +using QuadPoint = math::Vector; /** * @brief Immutable consumer interface for a quadrature rule on a reference cell. diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp index 7e8341ec9..5d12bb347 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp @@ -768,41 +768,32 @@ TEST(QuadraturePhase01Baseline, LegacyPositionModifiersAreOnlyDegreeOneExact) 1.0e-3); } -TEST(QuadPointContract, DefaultAndPartialConstructionZeroFillCoordinates) +TEST(QuadPointContract, UsesFixedSizeFEVectorRepresentation) { - static_assert(!std::is_convertible_v); + static_assert(std::is_same_v>); + static_assert(QuadPoint::RowsAtCompileTime == 3); + static_assert(QuadPoint::ColsAtCompileTime == 1); - const QuadPoint origin; + const QuadPoint origin = QuadPoint::Zero(); for (std::size_t component = 0; component < 3u; ++component) { EXPECT_DOUBLE_EQ(origin[component], 0.0); - EXPECT_DOUBLE_EQ( - origin.coordinates()[static_cast(component)], 0.0); } - const QuadPoint line_point{0.25}; + const QuadPoint line_point{0.25, 0.0, 0.0}; EXPECT_DOUBLE_EQ(line_point[0], 0.25); EXPECT_DOUBLE_EQ(line_point[1], 0.0); EXPECT_DOUBLE_EQ(line_point[2], 0.0); - const QuadPoint surface_point{0.25, 0.5}; + const QuadPoint surface_point{0.25, 0.5, 0.0}; EXPECT_DOUBLE_EQ(surface_point[0], 0.25); EXPECT_DOUBLE_EQ(surface_point[1], 0.5); EXPECT_DOUBLE_EQ(surface_point[2], 0.0); - QuadPoint::CoordinateVector coordinates = - QuadPoint::CoordinateVector::Zero(); - coordinates[0] = 0.2; - coordinates[1] = 0.3; - const QuadPoint vector_point{coordinates}; - EXPECT_DOUBLE_EQ(vector_point[0], 0.2); - EXPECT_DOUBLE_EQ(vector_point[1], 0.3); - EXPECT_DOUBLE_EQ(vector_point[2], 0.0); - - QuadPoint mutable_point; - mutable_point.coordinates()[2] = 0.75; + QuadPoint mutable_point = QuadPoint::Zero(); + mutable_point[2] = 0.75; EXPECT_DOUBLE_EQ(mutable_point[2], 0.75); - const std::vector points(2); + const std::vector points(2, QuadPoint::Zero()); for (const auto& point : points) { for (std::size_t component = 0; component < 3u; ++component) { EXPECT_DOUBLE_EQ(point[component], 0.0); @@ -1033,7 +1024,7 @@ TEST(QuadratureRuleValidation, RejectsIncorrectMeasureDespiteLargeCancellation) TEST(QuadratureRuleValidation, AcceptsWellConditionedLargePositiveRule) { constexpr std::size_t sample_count = 8192u; - std::vector points(sample_count); + std::vector points(sample_count, QuadPoint::Zero()); std::vector weights( sample_count, 2.0 / static_cast(sample_count)); @@ -1053,7 +1044,7 @@ TEST(QuadratureRuleValidation, AcceptsWellConditionedLargeSignedRule) { constexpr std::size_t sample_count = 8192u; constexpr double cancelling_weight = 1.0 / 1048576.0; - std::vector points(sample_count); + std::vector points(sample_count, QuadPoint::Zero()); std::vector weights( sample_count, 2.0 / static_cast(sample_count - 2u)); From f9ec988af27244e279e78bd43205f753e7c82495 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 23 Jul 2026 11:20:39 -0700 Subject: [PATCH 08/21] Separate quadrature dimensions and zeroth moment --- .../solver/FE/Quadrature/QuadratureRule.cpp | 42 ++++--- .../solver/FE/Quadrature/QuadratureRule.h | 108 ++++++++++++------ .../FE/Quadrature/test_QuadratureRules.cpp | 42 ++++--- 3 files changed, 125 insertions(+), 67 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp index dca6822c7..9939d4194 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp @@ -22,8 +22,9 @@ namespace svmp::FE::quadrature { namespace { struct ReferenceCellTraits { - int dimension; - double measure; + int integration_dimension; + int coordinate_dimension; + double zeroth_moment; }; enum class ValidationFailure { @@ -61,19 +62,19 @@ constexpr std::optional reference_cell_traits( { switch (family) { case svmp::CellFamily::Point: - return ReferenceCellTraits{0, 1.0}; + return ReferenceCellTraits{0, 0, 1.0}; case svmp::CellFamily::Line: - return ReferenceCellTraits{1, 2.0}; + return ReferenceCellTraits{1, 1, 2.0}; case svmp::CellFamily::Triangle: - return ReferenceCellTraits{2, 0.5}; + return ReferenceCellTraits{2, 2, 0.5}; case svmp::CellFamily::Quad: - return ReferenceCellTraits{2, 4.0}; + return ReferenceCellTraits{2, 2, 4.0}; case svmp::CellFamily::Tetra: - return ReferenceCellTraits{3, 1.0 / 6.0}; + return ReferenceCellTraits{3, 3, 1.0 / 6.0}; case svmp::CellFamily::Hex: - return ReferenceCellTraits{3, 8.0}; + return ReferenceCellTraits{3, 3, 8.0}; case svmp::CellFamily::Wedge: - return ReferenceCellTraits{3, 1.0}; + return ReferenceCellTraits{3, 3, 1.0}; default: return std::nullopt; } @@ -90,7 +91,8 @@ ValidationResult validate_point( if (!std::isfinite(point[component])) { return {ValidationFailure::NonFiniteCoordinate, sample}; } - if (component >= static_cast(traits.dimension) && + if (component >= + static_cast(traits.coordinate_dimension) && std::abs(point[component]) > tolerance) { return {ValidationFailure::NonzeroInactiveCoordinate, sample}; } @@ -153,7 +155,7 @@ double add_as_binary64(double left, double right) noexcept ValidationResult validate_weights( const std::vector& weights, - double reference_measure, + double zeroth_moment, double tolerance) noexcept { double ordered_sum = 0.0; @@ -195,7 +197,7 @@ ValidationResult validate_weights( return {ValidationFailure::NonFiniteWeightAccumulation, weights.size() - 1u}; } - const long double expected = static_cast(reference_measure); + const long double expected = static_cast(zeroth_moment); const long double scale = std::max(1.0L, std::abs(expected)); const long double requested_error_budget = static_cast(tolerance) * scale; @@ -262,7 +264,7 @@ ValidationResult validate_rule_data( } } - return validate_weights(weights, traits.measure, tolerance); + return validate_weights(weights, traits.zeroth_moment, tolerance); } std::string validation_failure_message(const ValidationResult& result) @@ -303,10 +305,10 @@ std::string validation_failure_message(const ValidationResult& result) return "QuadratureRule: weight accumulation is not finite" + sample_suffix(); case ValidationFailure::IncorrectZerothMoment: - return "QuadratureRule: weights do not reproduce the reference-cell measure"; + return "QuadratureRule: weights do not reproduce the zeroth moment"; case ValidationFailure::UnstableStoredOrderAccumulation: return "QuadratureRule: stored-order double-precision weight accumulation " - "does not reproduce the reference-cell measure"; + "does not reproduce the zeroth moment"; case ValidationFailure::IllConditionedWeightCancellation: return "QuadratureRule: weight cancellation is too ill-conditioned for " "stable double-precision integration"; @@ -325,9 +327,10 @@ QuadratureRule::QuadratureRule(svmp::CellFamily family, RuleData data) QuadratureRule::QuadratureRule(ValidatedState state) : cell_family_(state.cell_family), - dimension_(state.dimension), + integration_dimension_(state.integration_dimension), + coordinate_dimension_(state.coordinate_dimension), polynomial_exactness_(state.polynomial_exactness), - reference_measure_(state.reference_measure), + zeroth_moment_(state.zeroth_moment), points_(std::move(state.points)), weights_(std::move(state.weights)) { @@ -358,9 +361,10 @@ QuadratureRule::ValidatedState QuadratureRule::validate( return { family, - traits->dimension, + traits->integration_dimension, + traits->coordinate_dimension, data.polynomial_exactness, - traits->measure, + traits->zeroth_moment, std::move(data.points), std::move(data.weights), }; diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.h b/Code/Source/solver/FE/Quadrature/QuadratureRule.h index 8bbb12dda..426ab91d4 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.h +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.h @@ -28,9 +28,9 @@ * \int_{\hat K} f(\hat x)\,d\hat x * \approx \sum_q w_q f(\hat x_q). * @f] - * A rule identifies its reference-cell family, reports its intrinsic dimension - * and declared polynomial exactness, and keeps every point paired with its - * corresponding weight. + * A rule identifies its reference-cell family, reports its integration and + * coordinate dimensions and declared polynomial exactness, and keeps every + * point paired with its corresponding weight. * * The module does not choose the exactness required by an equation term, apply * reduced-integration policy, select a basis, own mesh storage, embed or orient @@ -62,10 +62,10 @@ * A rule is complete and structurally valid when construction returns. The * constructor rejects unsupported cells, negative exactness, empty or mismatched * storage, non-finite coordinates or weights, points outside the declared - * reference cell, and weights whose zeroth moment does not equal the reference - * measure in both compensated arithmetic and ordinary stored-order double - * accumulation. A condition estimate based on the absolute weight sum rejects - * signed rules whose cancellation is too sensitive for double precision. + * reference cell, and weights whose sum does not equal the canonical rule's + * zeroth moment in both compensated arithmetic and ordinary stored-order + * double accumulation. A condition estimate based on the absolute weight sum + * rejects signed rules whose cancellation is too sensitive for double precision. * Negative individual weights remain valid because some quadrature families * require them. * @@ -77,9 +77,9 @@ * * Points use one fixed-size three-component representation. Providers initialize * all three coordinates explicitly because the Eigen-backed vector is not - * zero-initialized by default. Only the first dimension() components are active, - * and every inactive component is zero within the construction tolerance. The - * supported canonical domains are: + * zero-initialized by default. Only the first coordinate_dimension() components + * are active, and every inactive component is zero within the construction + * tolerance. The supported canonical domains are: * * | Cell family | Canonical reference domain | Measure | * | ----------- | -------------------------- | ------- | @@ -118,9 +118,10 @@ namespace svmp::FE::quadrature { /** * @brief Three-component coordinate used for every reference quadrature point. * - * Only the first QuadratureRule::dimension() components are active. Providers - * explicitly zero remaining components, giving point, line, surface, and volume - * rules a uniform representation directly compatible with FE math consumers. + * Only the first QuadratureRule::coordinate_dimension() components are active. + * Providers explicitly zero remaining components, giving point, line, surface, + * and volume rules a uniform representation directly compatible with FE math + * consumers. */ using QuadPoint = math::Vector; @@ -166,10 +167,36 @@ class QuadratureRule { int polynomial_exactness() const noexcept { return polynomial_exactness_; } /** - * @brief Return the intrinsic dimension of the reference integration domain. - * @return Active coordinate count, from zero for Point through three for volume cells. + * @brief Return the intrinsic dimension of the integration domain. + * + * This is the dimension of the measure being integrated. It is distinct + * from coordinate_dimension() so future embedded rules can integrate a + * lower-dimensional entity represented in parent-reference coordinates. + * + * @return Integration dimension, from zero for Point through three for volume cells. + */ + int integration_dimension() const noexcept { return integration_dimension_; } + + /** + * @brief Return the number of active reference-coordinate components. + * + * Canonical reference-cell rules currently have the same integration and + * coordinate dimensions. Keeping the concepts distinct permits future + * embedded rules to use parent-reference coordinates without changing the + * meaning of integration_dimension(). + * + * @return Active coordinate count in each QuadPoint. */ - int dimension() const noexcept { return dimension_; } + int coordinate_dimension() const noexcept { return coordinate_dimension_; } + + /** + * @brief Return the intrinsic integration dimension. + * + * Compatibility alias for integration_dimension(). + * + * @return Integration dimension. + */ + int dimension() const noexcept { return integration_dimension(); } /** * @brief Return the canonical reference-cell family. @@ -211,7 +238,7 @@ class QuadratureRule { * * Coordinate bounds and inactive components use this as an absolute * tolerance. Compensated and ordinary stored-order zeroth-moment checks - * scale it by the larger of one and the reference-cell measure. For signed + * scale it by the larger of one and the zeroth moment. For signed * rules, a separate cancellation-sensitivity bound based on double epsilon * and the absolute weight sum always uses this default tolerance, so callers * cannot waive the minimum double-precision stability requirement. @@ -228,17 +255,32 @@ class QuadratureRule { * polynomial_exactness(). Construction already performs this check with * default_validation_tolerance(). * - * @param tol Finite, non-negative absolute coordinate and scaled measure tolerance. + * @param tol Finite, non-negative absolute coordinate and scaled zeroth-moment tolerance. * @return True when the stored rule satisfies every structural invariant. */ bool is_structurally_valid( double tol = default_validation_tolerance()) const noexcept; + /** + * @brief Return the rule's zeroth moment. + * + * This is the integral of the constant function one under the rule's + * integration measure. For current canonical unweighted rules it equals + * the geometric measure of the reference cell. + * + * @return Expected sum of the quadrature weights. + */ + double zeroth_moment() const noexcept { return zeroth_moment_; } + /** * @brief Return the measure of the canonical reference cell. + * + * Compatibility alias for zeroth_moment(). All currently supported rules + * are unweighted rules on complete canonical reference cells. + * * @return Reference length, area, volume, or unit point measure. */ - double reference_measure() const noexcept { return reference_measure_; } + double reference_measure() const noexcept { return zeroth_moment(); } protected: /** @@ -257,16 +299,16 @@ class QuadratureRule { /** * @brief Construct and validate one complete immutable rule. * - * Dimension and reference measure are derived from @p family; callers cannot - * supply redundant topology metadata. + * Integration dimension, coordinate dimension, and zeroth moment are + * derived from @p family; callers cannot supply redundant topology metadata. * * @param family Supported canonical reference-cell family. * @param data Complete exactness, point, and weight payload. * @throws InvalidArgumentException If the family is unsupported, exactness * is negative, storage is empty or mismatched, a value is non-finite, a point - * is outside the reference cell, the weights do not reproduce its measure in - * compensated and stored-order arithmetic, or their cancellation is too - * ill-conditioned for stable double-precision integration. + * is outside the reference cell, the weights do not reproduce its zeroth + * moment in compensated and stored-order arithmetic, or their cancellation + * is too ill-conditioned for stable double-precision integration. */ explicit QuadratureRule(svmp::CellFamily family, RuleData data); @@ -274,9 +316,10 @@ class QuadratureRule { /** @brief Fully checked state used by the delegating constructor. */ struct ValidatedState { svmp::CellFamily cell_family; - int dimension; + int integration_dimension; + int coordinate_dimension; int polynomial_exactness; - double reference_measure; + double zeroth_moment; std::vector points; std::vector weights; }; @@ -287,12 +330,13 @@ class QuadratureRule { /** @brief Initialize members from state already checked by validate(). */ explicit QuadratureRule(ValidatedState state); - const svmp::CellFamily cell_family_; ///< Canonical reference topology. - const int dimension_; ///< Number of active coordinate components. - const int polynomial_exactness_; ///< Exactness declared by the concrete generator. - const double reference_measure_; ///< Canonical cell measure. - const std::vector points_; ///< Ordered immutable reference coordinates. - const std::vector weights_; ///< Immutable weights paired with points_. + const svmp::CellFamily cell_family_; ///< Canonical reference topology. + const int integration_dimension_; ///< Dimension of the integration measure. + const int coordinate_dimension_; ///< Number of active coordinate components. + const int polynomial_exactness_; ///< Exactness declared by the concrete generator. + const double zeroth_moment_; ///< Integral of one under the rule's measure. + const std::vector points_; ///< Ordered immutable reference coordinates. + const std::vector weights_; ///< Immutable weights paired with points_. }; /** @} */ diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp index 5d12bb347..959238dcc 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp @@ -425,12 +425,13 @@ void expect_total_degree_exact( int degree, double tolerance = kExactnessTol) { - const int max_px = rule.dimension() >= 1 ? degree : 0; + const int max_px = rule.integration_dimension() >= 1 ? degree : 0; for (int px = 0; px <= max_px; ++px) { - const int max_py = rule.dimension() >= 2 ? degree - px : 0; + const int max_py = + rule.integration_dimension() >= 2 ? degree - px : 0; for (int py = 0; py <= max_py; ++py) { const int max_pz = - rule.dimension() >= 3 ? degree - px - py : 0; + rule.integration_dimension() >= 3 ? degree - px - py : 0; for (int pz = 0; pz <= max_pz; ++pz) { EXPECT_NEAR( integrate_monomial(rule, px, py, pz), @@ -580,10 +581,13 @@ TEST(QuadraturePhase01Baseline, CanonicalFixturesSatisfyTheRuleAndExactnessContr const auto rule = make_rule(c); EXPECT_EQ(rule.cell_family(), to_mesh_family(c.fe_type)); - EXPECT_EQ(rule.dimension(), c.dimension); + EXPECT_EQ(rule.integration_dimension(), c.dimension); + EXPECT_EQ(rule.coordinate_dimension(), c.dimension); + EXPECT_EQ(rule.dimension(), rule.integration_dimension()); EXPECT_EQ(rule.polynomial_exactness(), c.advertised_exactness); EXPECT_EQ(rule.num_points(), c.samples.size()); - EXPECT_DOUBLE_EQ(rule.reference_measure(), c.reference_measure); + EXPECT_DOUBLE_EQ(rule.zeroth_moment(), c.reference_measure); + EXPECT_DOUBLE_EQ(rule.reference_measure(), rule.zeroth_moment()); EXPECT_TRUE(rule.is_structurally_valid()); expect_samples_in_order(rule, c.samples); expect_total_degree_exact(rule, c.advertised_exactness); @@ -826,8 +830,11 @@ TEST(QuadratureRuleValidation, AcceptsEverySupportedReferenceCell) 0, {{c.point[0], c.point[1], c.point[2]}}, {c.expected_measure}); - EXPECT_EQ(rule.dimension(), c.expected_dimension); - EXPECT_DOUBLE_EQ(rule.reference_measure(), c.expected_measure); + EXPECT_EQ(rule.integration_dimension(), c.expected_dimension); + EXPECT_EQ(rule.coordinate_dimension(), c.expected_dimension); + EXPECT_EQ(rule.dimension(), rule.integration_dimension()); + EXPECT_DOUBLE_EQ(rule.zeroth_moment(), c.expected_measure); + EXPECT_DOUBLE_EQ(rule.reference_measure(), rule.zeroth_moment()); EXPECT_TRUE(rule.is_structurally_valid()); } } @@ -944,7 +951,7 @@ TEST(QuadratureRuleValidation, RejectsInactiveCoordinatesAndOutOfCellPoints) InvalidArgumentException); } -TEST(QuadratureRuleValidation, EnforcesMeasureButAllowsNegativeWeights) +TEST(QuadratureRuleValidation, EnforcesZerothMomentButAllowsNegativeWeights) { expect_invalid_argument_with_message( [] { @@ -954,7 +961,7 @@ TEST(QuadratureRuleValidation, EnforcesMeasureButAllowsNegativeWeights) {{0.25, 0.25, 0.0}}, {1.0}); }, - "weights do not reproduce the reference-cell measure"); + "weights do not reproduce the zeroth moment"); const RuleProbe rule( svmp::CellFamily::Triangle, @@ -966,7 +973,7 @@ TEST(QuadratureRuleValidation, EnforcesMeasureButAllowsNegativeWeights) EXPECT_LT(rule.weight(0), 0.0); EXPECT_DOUBLE_EQ( integrate_monomial(rule, 0, 0, 0), - rule.reference_measure()); + rule.zeroth_moment()); EXPECT_FALSE(rule.is_structurally_valid(-1.0)); EXPECT_FALSE(rule.is_structurally_valid( std::numeric_limits::quiet_NaN())); @@ -1006,7 +1013,7 @@ TEST(QuadratureRuleValidation, RejectsIllConditionedCancellationDespiteExactOrde "too ill-conditioned for stable double-precision integration"); } -TEST(QuadratureRuleValidation, RejectsIncorrectMeasureDespiteLargeCancellation) +TEST(QuadratureRuleValidation, RejectsIncorrectZerothMomentDespiteLargeCancellation) { expect_invalid_argument_with_message( [] { @@ -1018,7 +1025,7 @@ TEST(QuadratureRuleValidation, RejectsIncorrectMeasureDespiteLargeCancellation) {0.5, 0.0, 0.0}}, {1.0e20, 0.0, -1.0e20}); }, - "weights do not reproduce the reference-cell measure"); + "weights do not reproduce the zeroth moment"); } TEST(QuadratureRuleValidation, AcceptsWellConditionedLargePositiveRule) @@ -1037,7 +1044,7 @@ TEST(QuadratureRuleValidation, AcceptsWellConditionedLargePositiveRule) EXPECT_TRUE(rule.is_structurally_valid()); EXPECT_DOUBLE_EQ( integrate_monomial(rule, 0, 0, 0), - rule.reference_measure()); + rule.zeroth_moment()); } TEST(QuadratureRuleValidation, AcceptsWellConditionedLargeSignedRule) @@ -1060,7 +1067,7 @@ TEST(QuadratureRuleValidation, AcceptsWellConditionedLargeSignedRule) EXPECT_LT(rule.weight(0), 0.0); EXPECT_NEAR( integrate_monomial(rule, 0, 0, 0), - rule.reference_measure(), + rule.zeroth_moment(), QuadratureRule::default_validation_tolerance()); } @@ -1118,9 +1125,12 @@ TEST(QuadratureRuleContract, PublishesOnlyACompleteImmutableQueryInterface) {1.0, 1.0}); EXPECT_EQ(rule.cell_family(), svmp::CellFamily::Line); - EXPECT_EQ(rule.dimension(), 1); + EXPECT_EQ(rule.integration_dimension(), 1); + EXPECT_EQ(rule.coordinate_dimension(), 1); + EXPECT_EQ(rule.dimension(), rule.integration_dimension()); EXPECT_EQ(rule.polynomial_exactness(), 3); - EXPECT_DOUBLE_EQ(rule.reference_measure(), 2.0); + EXPECT_DOUBLE_EQ(rule.zeroth_moment(), 2.0); + EXPECT_DOUBLE_EQ(rule.reference_measure(), rule.zeroth_moment()); EXPECT_TRUE(rule.is_structurally_valid()); ASSERT_EQ(rule.num_points(), 2u); ASSERT_EQ(rule.points().size(), 2u); From 2e48f4c5e772ef0ddd0fd72cdc1ee2eb2256616c Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 23 Jul 2026 11:39:58 -0700 Subject: [PATCH 09/21] Inline binary64 weight accumulation --- Code/Source/solver/FE/Quadrature/QuadratureRule.cpp | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp index 9939d4194..dc7dc8484 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp @@ -147,12 +147,6 @@ ValidationResult validate_point( return {}; } -double add_as_binary64(double left, double right) noexcept -{ - volatile double rounded_sum = left + right; - return rounded_sum; -} - ValidationResult validate_weights( const std::vector& weights, double zeroth_moment, @@ -170,7 +164,10 @@ ValidationResult validate_weights( return {ValidationFailure::NonFiniteWeight, sample}; } - ordered_sum = add_as_binary64(ordered_sum, weight); + // The volatile store makes each stored-order addition round to + // binary64 even if intermediate expressions retain excess precision. + volatile double rounded_sum = ordered_sum + weight; + ordered_sum = rounded_sum; if (!std::isfinite(ordered_sum)) { return {ValidationFailure::NonFiniteWeightAccumulation, sample}; } From 0b4d7bf8e9677764d2289efde18e63deb4e8ac73 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 23 Jul 2026 12:01:13 -0700 Subject: [PATCH 10/21] Simplify quadrature validation diagnostics --- .../solver/FE/Quadrature/QuadratureRule.cpp | 112 ++++++------------ 1 file changed, 34 insertions(+), 78 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp index dc7dc8484..be37cac3c 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp @@ -16,6 +16,7 @@ #include #include #include +#include #include namespace svmp::FE::quadrature { @@ -27,33 +28,16 @@ struct ReferenceCellTraits { double zeroth_moment; }; -enum class ValidationFailure { - None, - UnsupportedCellFamily, - NegativePolynomialExactness, - InvalidTolerance, - EmptyRule, - PointWeightSizeMismatch, - NonFiniteCoordinate, - NonzeroInactiveCoordinate, - PointOutsideReferenceCell, - NonFiniteWeight, - NonFiniteWeightAccumulation, - IncorrectZerothMoment, - UnstableStoredOrderAccumulation, - IllConditionedWeightCancellation, -}; - struct ValidationResult { static constexpr std::size_t no_sample = std::numeric_limits::max(); - ValidationFailure failure{ValidationFailure::None}; + std::string_view reason{}; std::size_t sample{no_sample}; constexpr bool valid() const noexcept { - return failure == ValidationFailure::None; + return reason.empty(); } }; @@ -89,12 +73,12 @@ ValidationResult validate_point( { for (std::size_t component = 0; component < 3u; ++component) { if (!std::isfinite(point[component])) { - return {ValidationFailure::NonFiniteCoordinate, sample}; + return {"quadrature point contains a non-finite coordinate", sample}; } if (component >= static_cast(traits.coordinate_dimension) && std::abs(point[component]) > tolerance) { - return {ValidationFailure::NonzeroInactiveCoordinate, sample}; + return {"quadrature point has a nonzero inactive coordinate", sample}; } } @@ -138,11 +122,13 @@ ValidationResult validate_point( in_interval(z, -1.0, 1.0); break; default: - return {ValidationFailure::UnsupportedCellFamily}; + return {"unsupported reference-cell family"}; } if (!is_contained) { - return {ValidationFailure::PointOutsideReferenceCell, sample}; + return { + "quadrature point lies outside the canonical reference cell", + sample}; } return {}; } @@ -161,7 +147,7 @@ ValidationResult validate_weights( for (std::size_t sample = 0; sample < weights.size(); ++sample) { const double weight = weights[sample]; if (!std::isfinite(weight)) { - return {ValidationFailure::NonFiniteWeight, sample}; + return {"quadrature weight must be finite", sample}; } // The volatile store makes each stored-order addition round to @@ -169,7 +155,7 @@ ValidationResult validate_weights( volatile double rounded_sum = ordered_sum + weight; ordered_sum = rounded_sum; if (!std::isfinite(ordered_sum)) { - return {ValidationFailure::NonFiniteWeightAccumulation, sample}; + return {"weight accumulation is not finite", sample}; } has_negative_weight = has_negative_weight || weight < 0.0; @@ -185,13 +171,13 @@ ValidationResult validate_weights( if (!std::isfinite(compensated_sum) || !std::isfinite(correction) || !std::isfinite(absolute_sum)) { - return {ValidationFailure::NonFiniteWeightAccumulation, sample}; + return {"weight accumulation is not finite", sample}; } } const long double total = compensated_sum + correction; if (!std::isfinite(total)) { - return {ValidationFailure::NonFiniteWeightAccumulation, weights.size() - 1u}; + return {"weight accumulation is not finite", weights.size() - 1u}; } const long double expected = static_cast(zeroth_moment); @@ -201,13 +187,15 @@ ValidationResult validate_weights( const long double compensated_error = std::abs(total - expected); if (compensated_error > requested_error_budget) { - return {ValidationFailure::IncorrectZerothMoment}; + return {"weights do not reproduce the zeroth moment"}; } const long double ordered_error = std::abs(static_cast(ordered_sum) - expected); if (ordered_error > requested_error_budget) { - return {ValidationFailure::UnstableStoredOrderAccumulation}; + return { + "stored-order double-precision weight accumulation does not " + "reproduce the zeroth moment"}; } if (has_negative_weight) { @@ -225,7 +213,9 @@ ValidationResult validate_weights( QuadratureRule::default_validation_tolerance()) * scale; if (!std::isfinite(stability_error_bound) || stability_error_bound > stability_budget) { - return {ValidationFailure::IllConditionedWeightCancellation}; + return { + "weight cancellation is too ill-conditioned for stable " + "double-precision integration"}; } } @@ -241,16 +231,16 @@ ValidationResult validate_rule_data( double tolerance) noexcept { if (polynomial_exactness < 0) { - return {ValidationFailure::NegativePolynomialExactness}; + return {"polynomial exactness must be non-negative"}; } if (!std::isfinite(tolerance) || tolerance < 0.0) { - return {ValidationFailure::InvalidTolerance}; + return {"validation tolerance must be finite and non-negative"}; } if (points.empty()) { - return {ValidationFailure::EmptyRule}; + return {"a rule must contain at least one sample"}; } if (points.size() != weights.size()) { - return {ValidationFailure::PointWeightSizeMismatch}; + return {"points/weights size mismatch"}; } for (std::size_t sample = 0; sample < points.size(); ++sample) { @@ -266,51 +256,17 @@ ValidationResult validate_rule_data( std::string validation_failure_message(const ValidationResult& result) { - const auto sample_suffix = [&result]() { - if (result.sample == ValidationResult::no_sample) { - return std::string{}; - } - return std::string{" at sample "} + std::to_string(result.sample); - }; + if (result.valid()) { + return {}; + } - switch (result.failure) { - case ValidationFailure::None: - return {}; - case ValidationFailure::UnsupportedCellFamily: - return "QuadratureRule: unsupported reference-cell family"; - case ValidationFailure::NegativePolynomialExactness: - return "QuadratureRule: polynomial exactness must be non-negative"; - case ValidationFailure::InvalidTolerance: - return "QuadratureRule: validation tolerance must be finite and non-negative"; - case ValidationFailure::EmptyRule: - return "QuadratureRule: a rule must contain at least one sample"; - case ValidationFailure::PointWeightSizeMismatch: - return "QuadratureRule: points/weights size mismatch"; - case ValidationFailure::NonFiniteCoordinate: - return "QuadratureRule: quadrature point contains a non-finite coordinate" + - sample_suffix(); - case ValidationFailure::NonzeroInactiveCoordinate: - return "QuadratureRule: quadrature point has a nonzero inactive coordinate" + - sample_suffix(); - case ValidationFailure::PointOutsideReferenceCell: - return "QuadratureRule: quadrature point lies outside the canonical reference cell" + - sample_suffix(); - case ValidationFailure::NonFiniteWeight: - return "QuadratureRule: quadrature weight must be finite" + - sample_suffix(); - case ValidationFailure::NonFiniteWeightAccumulation: - return "QuadratureRule: weight accumulation is not finite" + - sample_suffix(); - case ValidationFailure::IncorrectZerothMoment: - return "QuadratureRule: weights do not reproduce the zeroth moment"; - case ValidationFailure::UnstableStoredOrderAccumulation: - return "QuadratureRule: stored-order double-precision weight accumulation " - "does not reproduce the zeroth moment"; - case ValidationFailure::IllConditionedWeightCancellation: - return "QuadratureRule: weight cancellation is too ill-conditioned for " - "stable double-precision integration"; + std::string message{"QuadratureRule: "}; + message.append(result.reason); + if (result.sample != ValidationResult::no_sample) { + message.append(" at sample "); + message.append(std::to_string(result.sample)); } - return "QuadratureRule: unknown validation failure"; + return message; } } // namespace @@ -341,7 +297,7 @@ QuadratureRule::ValidatedState QuadratureRule::validate( if (!traits) { svmp::raise( validation_failure_message( - {ValidationFailure::UnsupportedCellFamily})); + {"unsupported reference-cell family"})); } const auto validation = validate_rule_data( From 02aff3077ccee46bdc556ae5695a3e9422513c99 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 23 Jul 2026 12:13:26 -0700 Subject: [PATCH 11/21] Keep quadrature validation policy internal --- .../solver/FE/Quadrature/QuadratureRule.cpp | 26 ++------- .../solver/FE/Quadrature/QuadratureRule.h | 28 --------- .../FE/Quadrature/test_QuadratureRules.cpp | 58 +++++++++++-------- 3 files changed, 38 insertions(+), 74 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp index be37cac3c..6930640d9 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp @@ -28,6 +28,8 @@ struct ReferenceCellTraits { double zeroth_moment; }; +constexpr double construction_validation_tolerance = 1.0e-12; + struct ValidationResult { static constexpr std::size_t no_sample = std::numeric_limits::max(); @@ -209,8 +211,7 @@ ValidationResult validate_weights( compensated_error + cancellation_safety_factor * epsilon * absolute_sum; const long double stability_budget = - static_cast( - QuadratureRule::default_validation_tolerance()) * scale; + static_cast(construction_validation_tolerance) * scale; if (!std::isfinite(stability_error_bound) || stability_error_bound > stability_budget) { return { @@ -233,9 +234,6 @@ ValidationResult validate_rule_data( if (polynomial_exactness < 0) { return {"polynomial exactness must be non-negative"}; } - if (!std::isfinite(tolerance) || tolerance < 0.0) { - return {"validation tolerance must be finite and non-negative"}; - } if (points.empty()) { return {"a rule must contain at least one sample"}; } @@ -306,7 +304,7 @@ QuadratureRule::ValidatedState QuadratureRule::validate( data.polynomial_exactness, data.points, data.weights, - default_validation_tolerance()); + construction_validation_tolerance); if (!validation.valid()) { svmp::raise( validation_failure_message(validation)); @@ -323,20 +321,4 @@ QuadratureRule::ValidatedState QuadratureRule::validate( }; } -bool QuadratureRule::is_structurally_valid(double tolerance) const noexcept -{ - const auto traits = reference_cell_traits(cell_family_); - if (!traits) { - return false; - } - return validate_rule_data( - cell_family_, - *traits, - polynomial_exactness_, - points_, - weights_, - tolerance) - .valid(); -} - } // namespace svmp::FE::quadrature diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.h b/Code/Source/solver/FE/Quadrature/QuadratureRule.h index 426ab91d4..949c6d004 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.h +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.h @@ -233,34 +233,6 @@ class QuadratureRule { */ const std::vector& weights() const noexcept { return weights_; } - /** - * @brief Return the default tolerance used during structural validation. - * - * Coordinate bounds and inactive components use this as an absolute - * tolerance. Compensated and ordinary stored-order zeroth-moment checks - * scale it by the larger of one and the zeroth moment. For signed - * rules, a separate cancellation-sensitivity bound based on double epsilon - * and the absolute weight sum always uses this default tolerance, so callers - * cannot waive the minimum double-precision stability requirement. - * - * @return Default absolute validation tolerance. - */ - static constexpr double default_validation_tolerance() noexcept { return 1.0e-12; } - - /** - * @brief Recheck the rule's structural invariants using a caller tolerance. - * - * This diagnostic verifies metadata, storage, finite values, point - * containment, and the zeroth moment. It does not prove - * polynomial_exactness(). Construction already performs this check with - * default_validation_tolerance(). - * - * @param tol Finite, non-negative absolute coordinate and scaled zeroth-moment tolerance. - * @return True when the stored rule satisfies every structural invariant. - */ - bool is_structurally_valid( - double tol = default_validation_tolerance()) const noexcept; - /** * @brief Return the rule's zeroth moment. * diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp index 959238dcc..ebca61c2c 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp @@ -588,7 +588,6 @@ TEST(QuadraturePhase01Baseline, CanonicalFixturesSatisfyTheRuleAndExactnessContr EXPECT_EQ(rule.num_points(), c.samples.size()); EXPECT_DOUBLE_EQ(rule.zeroth_moment(), c.reference_measure); EXPECT_DOUBLE_EQ(rule.reference_measure(), rule.zeroth_moment()); - EXPECT_TRUE(rule.is_structurally_valid()); expect_samples_in_order(rule, c.samples); expect_total_degree_exact(rule, c.advertised_exactness); } @@ -835,7 +834,6 @@ TEST(QuadratureRuleValidation, AcceptsEverySupportedReferenceCell) EXPECT_EQ(rule.dimension(), rule.integration_dimension()); EXPECT_DOUBLE_EQ(rule.zeroth_moment(), c.expected_measure); EXPECT_DOUBLE_EQ(rule.reference_measure(), rule.zeroth_moment()); - EXPECT_TRUE(rule.is_structurally_valid()); } } @@ -968,17 +966,10 @@ TEST(QuadratureRuleValidation, EnforcesZerothMomentButAllowsNegativeWeights) 0, {{1.0 / 3.0, 1.0 / 3.0, 0.0}, {0.2, 0.2, 0.0}}, {-0.25, 0.75}); - EXPECT_TRUE(rule.is_structurally_valid()); - EXPECT_TRUE(rule.is_structurally_valid(0.0)); EXPECT_LT(rule.weight(0), 0.0); EXPECT_DOUBLE_EQ( integrate_monomial(rule, 0, 0, 0), rule.zeroth_moment()); - EXPECT_FALSE(rule.is_structurally_valid(-1.0)); - EXPECT_FALSE(rule.is_structurally_valid( - std::numeric_limits::quiet_NaN())); - EXPECT_FALSE(rule.is_structurally_valid( - std::numeric_limits::infinity())); } TEST(QuadratureRuleValidation, RejectsStoredOrderCancellation) @@ -1041,7 +1032,6 @@ TEST(QuadratureRuleValidation, AcceptsWellConditionedLargePositiveRule) 0, std::move(points), std::move(weights)); - EXPECT_TRUE(rule.is_structurally_valid()); EXPECT_DOUBLE_EQ( integrate_monomial(rule, 0, 0, 0), rule.zeroth_moment()); @@ -1063,34 +1053,55 @@ TEST(QuadratureRuleValidation, AcceptsWellConditionedLargeSignedRule) 0, std::move(points), std::move(weights)); - EXPECT_TRUE(rule.is_structurally_valid()); EXPECT_LT(rule.weight(0), 0.0); EXPECT_NEAR( integrate_monomial(rule, 0, 0, 0), rule.zeroth_moment(), - QuadratureRule::default_validation_tolerance()); + kTol); } -TEST(QuadratureRuleValidation, UsesDocumentedCoordinateTolerance) +TEST(QuadratureRuleValidation, AppliesConstructionCoordinateTolerance) { - constexpr double offset = 0.5 * QuadratureRule::default_validation_tolerance(); + constexpr double accepted_offset = 0.5 * kTol; const RuleProbe rule( - svmp::CellFamily::Line, 0, {{1.0 + offset, offset, 0.0}}, {2.0}); - EXPECT_TRUE(rule.is_structurally_valid()); - EXPECT_FALSE(rule.is_structurally_valid(offset / 2.0)); + svmp::CellFamily::Line, + 0, + {{1.0 + accepted_offset, accepted_offset, 0.0}}, + {2.0}); + EXPECT_DOUBLE_EQ(rule.point(0)[0], 1.0 + accepted_offset); + + constexpr double rejected_offset = 2.0 * kTol; + expect_invalid_argument_with_message( + [rejected_offset] { + (void)RuleProbe( + svmp::CellFamily::Line, + 0, + {{1.0 + rejected_offset, 0.0, 0.0}}, + {2.0}); + }, + "outside the canonical reference cell"); } -TEST(QuadratureRuleValidation, UsesDocumentedWeightTolerance) +TEST(QuadratureRuleValidation, AppliesConstructionWeightTolerance) { - constexpr double offset = - 0.5 * QuadratureRule::default_validation_tolerance(); + constexpr double accepted_offset = 0.5 * kTol; const RuleProbe rule( svmp::CellFamily::Triangle, 0, {{0.25, 0.25, 0.0}}, - {0.5 + offset}); - EXPECT_TRUE(rule.is_structurally_valid()); - EXPECT_FALSE(rule.is_structurally_valid(offset / 2.0)); + {0.5 + accepted_offset}); + EXPECT_DOUBLE_EQ(rule.weight(0), 0.5 + accepted_offset); + + constexpr double rejected_offset = 2.0 * kTol; + expect_invalid_argument_with_message( + [rejected_offset] { + (void)RuleProbe( + svmp::CellFamily::Triangle, + 0, + {{0.25, 0.25, 0.0}}, + {0.5 + rejected_offset}); + }, + "weights do not reproduce the zeroth moment"); } TEST(QuadratureRuleContract, PublishesOnlyACompleteImmutableQueryInterface) @@ -1131,7 +1142,6 @@ TEST(QuadratureRuleContract, PublishesOnlyACompleteImmutableQueryInterface) EXPECT_EQ(rule.polynomial_exactness(), 3); EXPECT_DOUBLE_EQ(rule.zeroth_moment(), 2.0); EXPECT_DOUBLE_EQ(rule.reference_measure(), rule.zeroth_moment()); - EXPECT_TRUE(rule.is_structurally_valid()); ASSERT_EQ(rule.num_points(), 2u); ASSERT_EQ(rule.points().size(), 2u); ASSERT_EQ(rule.weights().size(), 2u); From bad918840b89f335e824b85483be59a529919c7a Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 23 Jul 2026 12:28:09 -0700 Subject: [PATCH 12/21] Remove redundant reference measure alias --- Code/Source/solver/FE/Quadrature/QuadratureRule.h | 10 ---------- .../unitTests/FE/Quadrature/test_QuadratureRules.cpp | 11 ++++------- 2 files changed, 4 insertions(+), 17 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.h b/Code/Source/solver/FE/Quadrature/QuadratureRule.h index 949c6d004..e46d0fd64 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.h +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.h @@ -244,16 +244,6 @@ class QuadratureRule { */ double zeroth_moment() const noexcept { return zeroth_moment_; } - /** - * @brief Return the measure of the canonical reference cell. - * - * Compatibility alias for zeroth_moment(). All currently supported rules - * are unweighted rules on complete canonical reference cells. - * - * @return Reference length, area, volume, or unit point measure. - */ - double reference_measure() const noexcept { return zeroth_moment(); } - protected: /** * @brief Complete construction payload used by concrete rule providers. diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp index ebca61c2c..2e948bfe6 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp @@ -58,7 +58,7 @@ struct SolverRuleCase { int dimension; int requested_exactness; int advertised_exactness; - double reference_measure; + double zeroth_moment; ExpectedSamples samples; }; @@ -564,13 +564,13 @@ TEST(QuadraturePhase01Baseline, StandardSelectionTableHasOrderedPointWeightData) EXPECT_FALSE(c.samples.empty()); EXPECT_GT(c.requested_exactness, 0); EXPECT_GE(c.advertised_exactness, c.requested_exactness); - EXPECT_GT(c.reference_measure, 0.0); + EXPECT_GT(c.zeroth_moment, 0.0); double weight_sum = 0.0; for (const auto& sample : c.samples) { weight_sum += sample.weight; } - EXPECT_NEAR(weight_sum, c.reference_measure, kTol); + EXPECT_NEAR(weight_sum, c.zeroth_moment, kTol); } } @@ -586,8 +586,7 @@ TEST(QuadraturePhase01Baseline, CanonicalFixturesSatisfyTheRuleAndExactnessContr EXPECT_EQ(rule.dimension(), rule.integration_dimension()); EXPECT_EQ(rule.polynomial_exactness(), c.advertised_exactness); EXPECT_EQ(rule.num_points(), c.samples.size()); - EXPECT_DOUBLE_EQ(rule.zeroth_moment(), c.reference_measure); - EXPECT_DOUBLE_EQ(rule.reference_measure(), rule.zeroth_moment()); + EXPECT_DOUBLE_EQ(rule.zeroth_moment(), c.zeroth_moment); expect_samples_in_order(rule, c.samples); expect_total_degree_exact(rule, c.advertised_exactness); } @@ -833,7 +832,6 @@ TEST(QuadratureRuleValidation, AcceptsEverySupportedReferenceCell) EXPECT_EQ(rule.coordinate_dimension(), c.expected_dimension); EXPECT_EQ(rule.dimension(), rule.integration_dimension()); EXPECT_DOUBLE_EQ(rule.zeroth_moment(), c.expected_measure); - EXPECT_DOUBLE_EQ(rule.reference_measure(), rule.zeroth_moment()); } } @@ -1141,7 +1139,6 @@ TEST(QuadratureRuleContract, PublishesOnlyACompleteImmutableQueryInterface) EXPECT_EQ(rule.dimension(), rule.integration_dimension()); EXPECT_EQ(rule.polynomial_exactness(), 3); EXPECT_DOUBLE_EQ(rule.zeroth_moment(), 2.0); - EXPECT_DOUBLE_EQ(rule.reference_measure(), rule.zeroth_moment()); ASSERT_EQ(rule.num_points(), 2u); ASSERT_EQ(rule.points().size(), 2u); ASSERT_EQ(rule.weights().size(), 2u); From 752bbee171fefffa0e7ee576768d82a4f18e9fce Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 23 Jul 2026 13:00:56 -0700 Subject: [PATCH 13/21] Clarify the quadrature rule contract --- Code/CMake/SimVascularExternals.cmake | 2 + .../solver/FE/Quadrature/QuadratureRule.h | 77 +++++++++++-------- .../FE/Quadrature/test_QuadratureRules.cpp | 3 - 3 files changed, 49 insertions(+), 33 deletions(-) diff --git a/Code/CMake/SimVascularExternals.cmake b/Code/CMake/SimVascularExternals.cmake index ed9e078ec..880882e3e 100644 --- a/Code/CMake/SimVascularExternals.cmake +++ b/Code/CMake/SimVascularExternals.cmake @@ -9,6 +9,8 @@ if(DOXYGEN_FOUND) configure_file(${SV_SOURCE_DIR}/../Documentation/Doxyfile ${SV_BINARY_DIR}/Doxyfile @ONLY) add_custom_target(doc + COMMAND "${CMAKE_COMMAND}" -E remove_directory + "${SV_SOURCE_DIR}/../Documentation/build" COMMAND "${DOXYGEN_EXECUTABLE}" "${SV_BINARY_DIR}/Doxyfile" WORKING_DIRECTORY "${SV_SOURCE_DIR}/.." COMMENT "Generating API documentation with Doxygen" VERBATIM diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.h b/Code/Source/solver/FE/Quadrature/QuadratureRule.h index e46d0fd64..f0e86b0d1 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.h +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.h @@ -25,9 +25,14 @@ * The Quadrature module owns ordered reference coordinates and weights used to * approximate an integral on a canonical reference cell: * @f[ - * \int_{\hat K} f(\hat x)\,d\hat x + * \int_{\hat K} f(\hat x)\,d\mu(\hat x) * \approx \sum_q w_q f(\hat x_q). * @f] + * For the current unweighted rules, @f$d\mu=d\hat x@f$. The zeroth + * moment is the analytic normalization + * @f[ + * M_0 = \int_{\hat K} 1\,d\mu = \sum_q w_q. + * @f] * A rule identifies its reference-cell family, reports its integration and * coordinate dimensions and declared polynomial exactness, and keeps every * point paired with its corresponding weight. @@ -56,24 +61,30 @@ * exposes no mutation afterward. It must advertise only exactness established * for every rule it supplies through analytic moment tests. Derivation is a * provider extension seam, not an integration-consumer customization point. + * The protected constructor is the only provider extension point; public + * metadata and sample queries are fixed, nonvirtual operations on validated + * base-class state. * * ## Rule contract * - * A rule is complete and structurally valid when construction returns. The - * constructor rejects unsupported cells, negative exactness, empty or mismatched - * storage, non-finite coordinates or weights, points outside the declared - * reference cell, and weights whose sum does not equal the canonical rule's - * zeroth moment in both compensated arithmetic and ordinary stored-order - * double accumulation. A condition estimate based on the absolute weight sum - * rejects signed rules whose cancellation is too sensitive for double precision. - * Negative individual weights remain valid because some quadrature families - * require them. + * Construction is the sole validity boundary: a rule is complete and + * structurally valid when its constructor returns, and consumers do not perform + * a separate revalidation step. The constructor rejects unsupported cells, + * negative exactness, empty or mismatched storage, non-finite coordinates or + * weights, points outside the declared reference cell, and weights whose sum + * does not equal the canonical rule's zeroth moment in both compensated + * arithmetic and ordinary stored-order double accumulation. A condition + * estimate based on the absolute weight sum rejects signed rules whose + * cancellation is too sensitive for double precision. * - * polynomial_exactness() is declared by the concrete construction algorithm. - * Structural validation verifies metadata, containment, finiteness, and the - * zeroth moment; it does not prove higher-order polynomial moments. Every - * concrete rule provider is responsible for making its declaration a guarantee - * by establishing the advertised exactness with analytic moment tests. + * Structural validation does not require unique points or nonzero, positive + * individual weights. It verifies metadata, containment, finiteness, and the + * zeroth moment; it does not prove higher-order polynomial moments. A + * polynomial exactness of @f$p@f$ guarantees every polynomial of total degree + * at most @f$p@f$. A rule can integrate selected higher-degree polynomials + * without increasing that common guarantee. Every concrete provider is + * responsible for establishing its advertised exactness with analytic moment + * tests. * * Points use one fixed-size three-component representation. Providers initialize * all three coordinates explicitly because the Eigen-backed vector is not @@ -81,8 +92,8 @@ * are active, and every inactive component is zero within the construction * tolerance. The supported canonical domains are: * - * | Cell family | Canonical reference domain | Measure | - * | ----------- | -------------------------- | ------- | + * | Cell family | Canonical reference domain | Zeroth moment | + * | ----------- | -------------------------- | ------------- | * | Point | @f$(0,0,0)@f$ | @f$1@f$ | * | Line | @f$[-1,1]@f$ | @f$2@f$ | * | Triangle | @f$\xi,\eta\geq0;\ \xi+\eta\leq1@f$ | @f$1/2@f$ | @@ -100,7 +111,8 @@ * once, retained through a const owning handle, and shared across integrations. * References returned by points() and weights() remain valid for the lifetime * of the rule. Consumers should retain that rule rather than copying its point - * count or samples into parallel authoritative state. + * count or samples into parallel authoritative state. Concurrent const access + * is safe while an owning handle keeps the rule alive. */ #include "FE/Common/Types.h" @@ -131,7 +143,11 @@ using QuadPoint = math::Vector; * Concrete rule providers initialize the base with one complete RuleData * payload. The payload is validated before construction returns, and the object * exposes no mutation or assignment path afterward. General solver consumers - * use only the public const query interface. + * use only the public const query interface. Providers supply data through the + * protected constructor; they do not override the fixed metadata or sample + * queries. Successful construction establishes the rule invariant; no public + * revalidation operation is required or exposed. Concurrent const access is + * safe while an owning handle keeps the rule alive. */ class QuadratureRule { public: @@ -159,8 +175,11 @@ class QuadratureRule { /** * @brief Return the total-degree polynomial exactness declared by the rule. * - * Structural validation does not independently prove this degree. - * Concrete rule providers establish it through analytic moment tests. + * A value @f$p@f$ guarantees exact integration of every polynomial with + * total degree at most @f$p@f$. The rule can also integrate selected + * higher-degree polynomials. Structural validation does not independently + * prove this guarantee; concrete providers establish it through analytic + * moment tests. * * @return Declared total degree that a conforming rule integrates exactly. */ @@ -172,6 +191,8 @@ class QuadratureRule { * This is the dimension of the measure being integrated. It is distinct * from coordinate_dimension() so future embedded rules can integrate a * lower-dimensional entity represented in parent-reference coordinates. + * For example, a face rule can have integration dimension two while using + * three active parent-reference coordinates. * * @return Integration dimension, from zero for Point through three for volume cells. */ @@ -189,15 +210,6 @@ class QuadratureRule { */ int coordinate_dimension() const noexcept { return coordinate_dimension_; } - /** - * @brief Return the intrinsic integration dimension. - * - * Compatibility alias for integration_dimension(). - * - * @return Integration dimension. - */ - int dimension() const noexcept { return integration_dimension(); } - /** * @brief Return the canonical reference-cell family. * @return Reference topology integrated by this rule. @@ -239,6 +251,9 @@ class QuadratureRule { * This is the integral of the constant function one under the rule's * integration measure. For current canonical unweighted rules it equals * the geometric measure of the reference cell. + * @f[ + * M_0 = \int_{\hat K} 1\,d\mu = \sum_q w_q. + * @f] * * @return Expected sum of the quadrature weights. */ @@ -266,6 +281,8 @@ class QuadratureRule { * * @param family Supported canonical reference-cell family. * @param data Complete exactness, point, and weight payload. + * @note Duplicate points and zero or negative weights remain admissible + * when every other rule invariant is satisfied. * @throws InvalidArgumentException If the family is unsupported, exactness * is negative, storage is empty or mismatched, a value is non-finite, a point * is outside the reference cell, the weights do not reproduce its zeroth diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp index 2e948bfe6..287ef5828 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp @@ -583,7 +583,6 @@ TEST(QuadraturePhase01Baseline, CanonicalFixturesSatisfyTheRuleAndExactnessContr EXPECT_EQ(rule.cell_family(), to_mesh_family(c.fe_type)); EXPECT_EQ(rule.integration_dimension(), c.dimension); EXPECT_EQ(rule.coordinate_dimension(), c.dimension); - EXPECT_EQ(rule.dimension(), rule.integration_dimension()); EXPECT_EQ(rule.polynomial_exactness(), c.advertised_exactness); EXPECT_EQ(rule.num_points(), c.samples.size()); EXPECT_DOUBLE_EQ(rule.zeroth_moment(), c.zeroth_moment); @@ -830,7 +829,6 @@ TEST(QuadratureRuleValidation, AcceptsEverySupportedReferenceCell) {c.expected_measure}); EXPECT_EQ(rule.integration_dimension(), c.expected_dimension); EXPECT_EQ(rule.coordinate_dimension(), c.expected_dimension); - EXPECT_EQ(rule.dimension(), rule.integration_dimension()); EXPECT_DOUBLE_EQ(rule.zeroth_moment(), c.expected_measure); } } @@ -1136,7 +1134,6 @@ TEST(QuadratureRuleContract, PublishesOnlyACompleteImmutableQueryInterface) EXPECT_EQ(rule.cell_family(), svmp::CellFamily::Line); EXPECT_EQ(rule.integration_dimension(), 1); EXPECT_EQ(rule.coordinate_dimension(), 1); - EXPECT_EQ(rule.dimension(), rule.integration_dimension()); EXPECT_EQ(rule.polynomial_exactness(), 3); EXPECT_DOUBLE_EQ(rule.zeroth_moment(), 2.0); ASSERT_EQ(rule.num_points(), 2u); From fba6733ad4c42347c2184384540360bad8810f0e Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 23 Jul 2026 15:46:48 -0700 Subject: [PATCH 14/21] Simplify quadrature rule validation --- .../solver/FE/Quadrature/QuadratureRule.cpp | 120 +++++++----------- .../solver/FE/Quadrature/QuadratureRule.h | 89 +++++-------- .../FE/Quadrature/test_QuadratureRules.cpp | 100 +++++---------- 3 files changed, 116 insertions(+), 193 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp index 6930640d9..bf26044dd 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp @@ -23,12 +23,12 @@ namespace svmp::FE::quadrature { namespace { struct ReferenceCellTraits { - int integration_dimension; - int coordinate_dimension; + int dimension; double zeroth_moment; }; -constexpr double construction_validation_tolerance = 1.0e-12; +constexpr double coordinate_validation_tolerance = 1.0e-12; +constexpr double moment_validation_tolerance = 1.0e-12; struct ValidationResult { static constexpr std::size_t no_sample = @@ -48,19 +48,19 @@ constexpr std::optional reference_cell_traits( { switch (family) { case svmp::CellFamily::Point: - return ReferenceCellTraits{0, 0, 1.0}; + return ReferenceCellTraits{0, 1.0}; case svmp::CellFamily::Line: - return ReferenceCellTraits{1, 1, 2.0}; + return ReferenceCellTraits{1, 2.0}; case svmp::CellFamily::Triangle: - return ReferenceCellTraits{2, 2, 0.5}; + return ReferenceCellTraits{2, 0.5}; case svmp::CellFamily::Quad: - return ReferenceCellTraits{2, 2, 4.0}; + return ReferenceCellTraits{2, 4.0}; case svmp::CellFamily::Tetra: - return ReferenceCellTraits{3, 3, 1.0 / 6.0}; + return ReferenceCellTraits{3, 1.0 / 6.0}; case svmp::CellFamily::Hex: - return ReferenceCellTraits{3, 3, 8.0}; + return ReferenceCellTraits{3, 8.0}; case svmp::CellFamily::Wedge: - return ReferenceCellTraits{3, 3, 1.0}; + return ReferenceCellTraits{3, 1.0}; default: return std::nullopt; } @@ -70,7 +70,6 @@ ValidationResult validate_point( const QuadPoint& point, svmp::CellFamily family, const ReferenceCellTraits& traits, - double tolerance, std::size_t sample) noexcept { for (std::size_t component = 0; component < 3u; ++component) { @@ -78,14 +77,15 @@ ValidationResult validate_point( return {"quadrature point contains a non-finite coordinate", sample}; } if (component >= - static_cast(traits.coordinate_dimension) && - std::abs(point[component]) > tolerance) { + static_cast(traits.dimension) && + std::abs(point[component]) > coordinate_validation_tolerance) { return {"quadrature point has a nonzero inactive coordinate", sample}; } } - const auto in_interval = [tolerance](double value, double lower, double upper) { - return value >= lower - tolerance && value <= upper + tolerance; + const auto in_interval = [](double value, double lower, double upper) { + return value >= lower - coordinate_validation_tolerance && + value <= upper + coordinate_validation_tolerance; }; const double x = point[0]; @@ -101,17 +101,21 @@ ValidationResult validate_point( is_contained = in_interval(x, -1.0, 1.0); break; case svmp::CellFamily::Triangle: - is_contained = x >= -tolerance && y >= -tolerance && - x + y <= 1.0 + tolerance; + is_contained = + x >= -coordinate_validation_tolerance && + y >= -coordinate_validation_tolerance && + x + y <= 1.0 + coordinate_validation_tolerance; break; case svmp::CellFamily::Quad: is_contained = in_interval(x, -1.0, 1.0) && in_interval(y, -1.0, 1.0); break; case svmp::CellFamily::Tetra: - is_contained = x >= -tolerance && y >= -tolerance && - z >= -tolerance && - x + y + z <= 1.0 + tolerance; + is_contained = + x >= -coordinate_validation_tolerance && + y >= -coordinate_validation_tolerance && + z >= -coordinate_validation_tolerance && + x + y + z <= 1.0 + coordinate_validation_tolerance; break; case svmp::CellFamily::Hex: is_contained = in_interval(x, -1.0, 1.0) && @@ -119,9 +123,11 @@ ValidationResult validate_point( in_interval(z, -1.0, 1.0); break; case svmp::CellFamily::Wedge: - is_contained = x >= -tolerance && y >= -tolerance && - x + y <= 1.0 + tolerance && - in_interval(z, -1.0, 1.0); + is_contained = + x >= -coordinate_validation_tolerance && + y >= -coordinate_validation_tolerance && + x + y <= 1.0 + coordinate_validation_tolerance && + in_interval(z, -1.0, 1.0); break; default: return {"unsupported reference-cell family"}; @@ -137,14 +143,11 @@ ValidationResult validate_point( ValidationResult validate_weights( const std::vector& weights, - double zeroth_moment, - double tolerance) noexcept + double zeroth_moment) noexcept { - double ordered_sum = 0.0; long double compensated_sum = 0.0L; long double correction = 0.0L; long double absolute_sum = 0.0L; - bool has_negative_weight = false; for (std::size_t sample = 0; sample < weights.size(); ++sample) { const double weight = weights[sample]; @@ -152,15 +155,6 @@ ValidationResult validate_weights( return {"quadrature weight must be finite", sample}; } - // The volatile store makes each stored-order addition round to - // binary64 even if intermediate expressions retain excess precision. - volatile double rounded_sum = ordered_sum + weight; - ordered_sum = rounded_sum; - if (!std::isfinite(ordered_sum)) { - return {"weight accumulation is not finite", sample}; - } - has_negative_weight = has_negative_weight || weight < 0.0; - const long double value = static_cast(weight); const long double next = compensated_sum + value; if (std::abs(compensated_sum) >= std::abs(value)) { @@ -185,39 +179,25 @@ ValidationResult validate_weights( const long double expected = static_cast(zeroth_moment); const long double scale = std::max(1.0L, std::abs(expected)); const long double requested_error_budget = - static_cast(tolerance) * scale; + static_cast(moment_validation_tolerance) * scale; const long double compensated_error = std::abs(total - expected); if (compensated_error > requested_error_budget) { return {"weights do not reproduce the zeroth moment"}; } - const long double ordered_error = - std::abs(static_cast(ordered_sum) - expected); - if (ordered_error > requested_error_budget) { + // A single compensated accumulator is not exact for unrestricted signed + // data. Use a conservative uncertainty estimate in the precision used by + // validation before accepting the computed moment. + constexpr long double validation_safety_factor = 8.0L; + const long double validation_error_bound = + compensated_error + + validation_safety_factor * + std::numeric_limits::epsilon() * absolute_sum; + if (!std::isfinite(validation_error_bound) || + validation_error_bound > requested_error_budget) { return { - "stored-order double-precision weight accumulation does not " - "reproduce the zeroth moment"}; - } - - if (has_negative_weight) { - // The actual stored-order sum is checked above. This independent guard - // limits sensitivity to signed-weight cancellation without imposing a - // sample-count ceiling on otherwise well-conditioned rules. - constexpr long double cancellation_safety_factor = 8.0L; - const long double epsilon = - static_cast(std::numeric_limits::epsilon()); - const long double stability_error_bound = - compensated_error + - cancellation_safety_factor * epsilon * absolute_sum; - const long double stability_budget = - static_cast(construction_validation_tolerance) * scale; - if (!std::isfinite(stability_error_bound) || - stability_error_bound > stability_budget) { - return { - "weight cancellation is too ill-conditioned for stable " - "double-precision integration"}; - } + "weight sum is too ill-conditioned to validate the zeroth moment"}; } return {}; @@ -228,8 +208,7 @@ ValidationResult validate_rule_data( const ReferenceCellTraits& traits, int polynomial_exactness, const std::vector& points, - const std::vector& weights, - double tolerance) noexcept + const std::vector& weights) noexcept { if (polynomial_exactness < 0) { return {"polynomial exactness must be non-negative"}; @@ -243,13 +222,13 @@ ValidationResult validate_rule_data( for (std::size_t sample = 0; sample < points.size(); ++sample) { const auto result = - validate_point(points[sample], family, traits, tolerance, sample); + validate_point(points[sample], family, traits, sample); if (!result.valid()) { return result; } } - return validate_weights(weights, traits.zeroth_moment, tolerance); + return validate_weights(weights, traits.zeroth_moment); } std::string validation_failure_message(const ValidationResult& result) @@ -278,8 +257,7 @@ QuadratureRule::QuadratureRule(svmp::CellFamily family, RuleData data) QuadratureRule::QuadratureRule(ValidatedState state) : cell_family_(state.cell_family), - integration_dimension_(state.integration_dimension), - coordinate_dimension_(state.coordinate_dimension), + dimension_(state.dimension), polynomial_exactness_(state.polynomial_exactness), zeroth_moment_(state.zeroth_moment), points_(std::move(state.points)), @@ -303,8 +281,7 @@ QuadratureRule::ValidatedState QuadratureRule::validate( *traits, data.polynomial_exactness, data.points, - data.weights, - construction_validation_tolerance); + data.weights); if (!validation.valid()) { svmp::raise( validation_failure_message(validation)); @@ -312,8 +289,7 @@ QuadratureRule::ValidatedState QuadratureRule::validate( return { family, - traits->integration_dimension, - traits->coordinate_dimension, + traits->dimension, data.polynomial_exactness, traits->zeroth_moment, std::move(data.points), diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.h b/Code/Source/solver/FE/Quadrature/QuadratureRule.h index f0e86b0d1..517fd74f6 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.h +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.h @@ -23,19 +23,18 @@ * ## Scope * * The Quadrature module owns ordered reference coordinates and weights used to - * approximate an integral on a canonical reference cell: + * approximate an unweighted integral on a canonical reference cell: * @f[ - * \int_{\hat K} f(\hat x)\,d\mu(\hat x) + * \int_{\hat K} f(\hat x)\,d\hat x * \approx \sum_q w_q f(\hat x_q). * @f] - * For the current unweighted rules, @f$d\mu=d\hat x@f$. The zeroth - * moment is the analytic normalization + * The zeroth moment is the analytic normalization * @f[ - * M_0 = \int_{\hat K} 1\,d\mu = \sum_q w_q. + * M_0 = \int_{\hat K} 1\,d\hat x = \sum_q w_q. * @f] - * A rule identifies its reference-cell family, reports its integration and - * coordinate dimensions and declared polynomial exactness, and keeps every - * point paired with its corresponding weight. + * A rule identifies its reference-cell family, reports its dimension and + * declared polynomial exactness, and keeps every point paired with its + * corresponding weight. * * The module does not choose the exactness required by an equation term, apply * reduced-integration policy, select a basis, own mesh storage, embed or orient @@ -72,10 +71,8 @@ * a separate revalidation step. The constructor rejects unsupported cells, * negative exactness, empty or mismatched storage, non-finite coordinates or * weights, points outside the declared reference cell, and weights whose sum - * does not equal the canonical rule's zeroth moment in both compensated - * arithmetic and ordinary stored-order double accumulation. A condition - * estimate based on the absolute weight sum rejects signed rules whose - * cancellation is too sensitive for double precision. + * does not equal the canonical rule's zeroth moment in compensated arithmetic + * or cannot be established within the validator's numerical precision. * * Structural validation does not require unique points or nonzero, positive * individual weights. It verifies metadata, containment, finiteness, and the @@ -88,8 +85,8 @@ * * Points use one fixed-size three-component representation. Providers initialize * all three coordinates explicitly because the Eigen-backed vector is not - * zero-initialized by default. Only the first coordinate_dimension() components - * are active, and every inactive component is zero within the construction + * zero-initialized by default. Only the first dimension() components + * are active, and every inactive component is zero within the coordinate * tolerance. The supported canonical domains are: * * | Cell family | Canonical reference domain | Zeroth moment | @@ -130,10 +127,9 @@ namespace svmp::FE::quadrature { /** * @brief Three-component coordinate used for every reference quadrature point. * - * Only the first QuadratureRule::coordinate_dimension() components are active. - * Providers explicitly zero remaining components, giving point, line, surface, - * and volume rules a uniform representation directly compatible with FE math - * consumers. + * Only the first QuadratureRule::dimension() components are active. Providers + * explicitly zero remaining components, giving point, line, surface, and volume + * rules a uniform representation directly compatible with FE math consumers. */ using QuadPoint = math::Vector; @@ -186,29 +182,13 @@ class QuadratureRule { int polynomial_exactness() const noexcept { return polynomial_exactness_; } /** - * @brief Return the intrinsic dimension of the integration domain. + * @brief Return the dimension of the canonical reference cell. * - * This is the dimension of the measure being integrated. It is distinct - * from coordinate_dimension() so future embedded rules can integrate a - * lower-dimensional entity represented in parent-reference coordinates. - * For example, a face rule can have integration dimension two while using - * three active parent-reference coordinates. + * This is also the number of active components in each QuadPoint. * - * @return Integration dimension, from zero for Point through three for volume cells. + * @return Reference dimension, from zero for Point through three for volume cells. */ - int integration_dimension() const noexcept { return integration_dimension_; } - - /** - * @brief Return the number of active reference-coordinate components. - * - * Canonical reference-cell rules currently have the same integration and - * coordinate dimensions. Keeping the concepts distinct permits future - * embedded rules to use parent-reference coordinates without changing the - * meaning of integration_dimension(). - * - * @return Active coordinate count in each QuadPoint. - */ - int coordinate_dimension() const noexcept { return coordinate_dimension_; } + int dimension() const noexcept { return dimension_; } /** * @brief Return the canonical reference-cell family. @@ -248,11 +228,12 @@ class QuadratureRule { /** * @brief Return the rule's zeroth moment. * - * This is the integral of the constant function one under the rule's - * integration measure. For current canonical unweighted rules it equals - * the geometric measure of the reference cell. + * This is the integral of the constant function one. All supported rules + * are unweighted rules on complete canonical reference cells, so the + * constructor derives this value from cell_family() and it equals the + * geometric measure of that cell. * @f[ - * M_0 = \int_{\hat K} 1\,d\mu = \sum_q w_q. + * M_0 = \int_{\hat K} 1\,d\hat x = \sum_q w_q. * @f] * * @return Expected sum of the quadrature weights. @@ -276,8 +257,8 @@ class QuadratureRule { /** * @brief Construct and validate one complete immutable rule. * - * Integration dimension, coordinate dimension, and zeroth moment are - * derived from @p family; callers cannot supply redundant topology metadata. + * Dimension and zeroth moment are derived from @p family; callers cannot + * supply redundant topology metadata. * * @param family Supported canonical reference-cell family. * @param data Complete exactness, point, and weight payload. @@ -286,8 +267,8 @@ class QuadratureRule { * @throws InvalidArgumentException If the family is unsupported, exactness * is negative, storage is empty or mismatched, a value is non-finite, a point * is outside the reference cell, the weights do not reproduce its zeroth - * moment in compensated and stored-order arithmetic, or their cancellation - * is too ill-conditioned for stable double-precision integration. + * moment in compensated arithmetic, or cancellation prevents validating + * that moment reliably. */ explicit QuadratureRule(svmp::CellFamily family, RuleData data); @@ -295,8 +276,7 @@ class QuadratureRule { /** @brief Fully checked state used by the delegating constructor. */ struct ValidatedState { svmp::CellFamily cell_family; - int integration_dimension; - int coordinate_dimension; + int dimension; int polynomial_exactness; double zeroth_moment; std::vector points; @@ -309,13 +289,12 @@ class QuadratureRule { /** @brief Initialize members from state already checked by validate(). */ explicit QuadratureRule(ValidatedState state); - const svmp::CellFamily cell_family_; ///< Canonical reference topology. - const int integration_dimension_; ///< Dimension of the integration measure. - const int coordinate_dimension_; ///< Number of active coordinate components. - const int polynomial_exactness_; ///< Exactness declared by the concrete generator. - const double zeroth_moment_; ///< Integral of one under the rule's measure. - const std::vector points_; ///< Ordered immutable reference coordinates. - const std::vector weights_; ///< Immutable weights paired with points_. + const svmp::CellFamily cell_family_; ///< Canonical reference topology. + const int dimension_; ///< Number of active coordinate components. + const int polynomial_exactness_; ///< Exactness declared by the concrete generator. + const double zeroth_moment_; ///< Canonical reference-cell measure. + const std::vector points_; ///< Ordered immutable reference coordinates. + const std::vector weights_; ///< Immutable weights paired with points_. }; /** @} */ diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp index 287ef5828..5586b55e3 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp @@ -425,13 +425,11 @@ void expect_total_degree_exact( int degree, double tolerance = kExactnessTol) { - const int max_px = rule.integration_dimension() >= 1 ? degree : 0; + const int max_px = rule.dimension() >= 1 ? degree : 0; for (int px = 0; px <= max_px; ++px) { - const int max_py = - rule.integration_dimension() >= 2 ? degree - px : 0; + const int max_py = rule.dimension() >= 2 ? degree - px : 0; for (int py = 0; py <= max_py; ++py) { - const int max_pz = - rule.integration_dimension() >= 3 ? degree - px - py : 0; + const int max_pz = rule.dimension() >= 3 ? degree - px - py : 0; for (int pz = 0; pz <= max_pz; ++pz) { EXPECT_NEAR( integrate_monomial(rule, px, py, pz), @@ -581,8 +579,7 @@ TEST(QuadraturePhase01Baseline, CanonicalFixturesSatisfyTheRuleAndExactnessContr const auto rule = make_rule(c); EXPECT_EQ(rule.cell_family(), to_mesh_family(c.fe_type)); - EXPECT_EQ(rule.integration_dimension(), c.dimension); - EXPECT_EQ(rule.coordinate_dimension(), c.dimension); + EXPECT_EQ(rule.dimension(), c.dimension); EXPECT_EQ(rule.polynomial_exactness(), c.advertised_exactness); EXPECT_EQ(rule.num_points(), c.samples.size()); EXPECT_DOUBLE_EQ(rule.zeroth_moment(), c.zeroth_moment); @@ -827,8 +824,7 @@ TEST(QuadratureRuleValidation, AcceptsEverySupportedReferenceCell) 0, {{c.point[0], c.point[1], c.point[2]}}, {c.expected_measure}); - EXPECT_EQ(rule.integration_dimension(), c.expected_dimension); - EXPECT_EQ(rule.coordinate_dimension(), c.expected_dimension); + EXPECT_EQ(rule.dimension(), c.expected_dimension); EXPECT_DOUBLE_EQ(rule.zeroth_moment(), c.expected_measure); } } @@ -899,16 +895,6 @@ TEST(QuadratureRuleValidation, RejectsMalformedStorageAndNonfiniteValues) svmp::CellFamily::Line, 1, {{0.0, 0.0, 0.0}}, {nan}); }, "quadrature weight must be finite at sample 0"); - expect_invalid_argument_with_message( - [] { - (void)RuleProbe( - svmp::CellFamily::Line, - 1, - {{0.0, 0.0, 0.0}, {0.5, 0.0, 0.0}}, - {std::numeric_limits::max(), - std::numeric_limits::max()}); - }, - "weight accumulation is not finite at sample 1"); } TEST(QuadratureRuleValidation, RejectsInactiveCoordinatesAndOutOfCellPoints) @@ -968,38 +954,6 @@ TEST(QuadratureRuleValidation, EnforcesZerothMomentButAllowsNegativeWeights) rule.zeroth_moment()); } -TEST(QuadratureRuleValidation, RejectsStoredOrderCancellation) -{ - expect_invalid_argument_with_message( - [] { - (void)RuleProbe( - svmp::CellFamily::Line, - 0, - {{-0.75, 0.0, 0.0}, - {-0.25, 0.0, 0.0}, - {0.25, 0.0, 0.0}, - {0.75, 0.0, 0.0}}, - {1.0e16, 1.0, -1.0e16, 1.0}); - }, - "stored-order double-precision weight accumulation"); -} - -TEST(QuadratureRuleValidation, RejectsIllConditionedCancellationDespiteExactOrderedSum) -{ - expect_invalid_argument_with_message( - [] { - (void)RuleProbe( - svmp::CellFamily::Line, - 0, - {{-0.75, 0.0, 0.0}, - {-0.25, 0.0, 0.0}, - {0.25, 0.0, 0.0}, - {0.75, 0.0, 0.0}}, - {1.0e16, -1.0e16, 1.0, 1.0}); - }, - "too ill-conditioned for stable double-precision integration"); -} - TEST(QuadratureRuleValidation, RejectsIncorrectZerothMomentDespiteLargeCancellation) { expect_invalid_argument_with_message( @@ -1015,25 +969,40 @@ TEST(QuadratureRuleValidation, RejectsIncorrectZerothMomentDespiteLargeCancellat "weights do not reproduce the zeroth moment"); } -TEST(QuadratureRuleValidation, AcceptsWellConditionedLargePositiveRule) +TEST(QuadratureRuleValidation, RejectsSignedDataWithUnreliableMoment) { - constexpr std::size_t sample_count = 8192u; + const double medium = std::ldexp(1.0, 93); + const double large = std::ldexp(1.0, 233); + const double residual = std::ldexp(1.0, 14); + + EXPECT_THROW( + (void)RuleProbe( + svmp::CellFamily::Line, + 0, + std::vector(6u, QuadPoint::Zero()), + {-medium, -large, residual, medium, large, 2.0}), + InvalidArgumentException); +} + +TEST(QuadratureRuleValidation, AcceptsLargePositiveRule) +{ + // Repeating a non-dyadic weight exposes ordinary left-fold drift without + // changing the mathematical normalization of the rule. + constexpr std::size_t sample_count = 100000u; + const double sample_weight = 2.0 / static_cast(sample_count); std::vector points(sample_count, QuadPoint::Zero()); - std::vector weights( - sample_count, - 2.0 / static_cast(sample_count)); + std::vector weights(sample_count, sample_weight); const RuleProbe rule( svmp::CellFamily::Line, 0, std::move(points), std::move(weights)); - EXPECT_DOUBLE_EQ( - integrate_monomial(rule, 0, 0, 0), - rule.zeroth_moment()); + EXPECT_EQ(rule.num_points(), sample_count); + EXPECT_DOUBLE_EQ(rule.weight(0), sample_weight); } -TEST(QuadratureRuleValidation, AcceptsWellConditionedLargeSignedRule) +TEST(QuadratureRuleValidation, AcceptsLargeSignedRule) { constexpr std::size_t sample_count = 8192u; constexpr double cancelling_weight = 1.0 / 1048576.0; @@ -1058,7 +1027,7 @@ TEST(QuadratureRuleValidation, AcceptsWellConditionedLargeSignedRule) TEST(QuadratureRuleValidation, AppliesConstructionCoordinateTolerance) { - constexpr double accepted_offset = 0.5 * kTol; + constexpr double accepted_offset = 1.0e-14; const RuleProbe rule( svmp::CellFamily::Line, 0, @@ -1066,7 +1035,7 @@ TEST(QuadratureRuleValidation, AppliesConstructionCoordinateTolerance) {2.0}); EXPECT_DOUBLE_EQ(rule.point(0)[0], 1.0 + accepted_offset); - constexpr double rejected_offset = 2.0 * kTol; + constexpr double rejected_offset = 1.0e-10; expect_invalid_argument_with_message( [rejected_offset] { (void)RuleProbe( @@ -1080,7 +1049,7 @@ TEST(QuadratureRuleValidation, AppliesConstructionCoordinateTolerance) TEST(QuadratureRuleValidation, AppliesConstructionWeightTolerance) { - constexpr double accepted_offset = 0.5 * kTol; + constexpr double accepted_offset = 1.0e-14; const RuleProbe rule( svmp::CellFamily::Triangle, 0, @@ -1088,7 +1057,7 @@ TEST(QuadratureRuleValidation, AppliesConstructionWeightTolerance) {0.5 + accepted_offset}); EXPECT_DOUBLE_EQ(rule.weight(0), 0.5 + accepted_offset); - constexpr double rejected_offset = 2.0 * kTol; + constexpr double rejected_offset = 1.0e-10; expect_invalid_argument_with_message( [rejected_offset] { (void)RuleProbe( @@ -1132,8 +1101,7 @@ TEST(QuadratureRuleContract, PublishesOnlyACompleteImmutableQueryInterface) {1.0, 1.0}); EXPECT_EQ(rule.cell_family(), svmp::CellFamily::Line); - EXPECT_EQ(rule.integration_dimension(), 1); - EXPECT_EQ(rule.coordinate_dimension(), 1); + EXPECT_EQ(rule.dimension(), 1); EXPECT_EQ(rule.polynomial_exactness(), 3); EXPECT_DOUBLE_EQ(rule.zeroth_moment(), 2.0); ASSERT_EQ(rule.num_points(), 2u); From b1433cb586fe3282028d67c426584071e0311b20 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 23 Jul 2026 16:24:16 -0700 Subject: [PATCH 15/21] Use exact summation for quadrature validation --- .../solver/FE/Quadrature/QuadratureRule.cpp | 237 ++++++++++++++---- .../solver/FE/Quadrature/QuadratureRule.h | 10 +- .../FE/Quadrature/test_QuadratureRules.cpp | 124 ++++++++- 3 files changed, 316 insertions(+), 55 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp index bf26044dd..d4db4ec64 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp @@ -12,7 +12,11 @@ #include "FE/Common/FEException.h" #include +#include +#include +#include #include +#include #include #include #include @@ -30,6 +34,188 @@ struct ReferenceCellTraits { constexpr double coordinate_validation_tolerance = 1.0e-12; constexpr double moment_validation_tolerance = 1.0e-12; +static_assert( + std::numeric_limits::is_iec559 && + std::numeric_limits::radix == 2 && + std::numeric_limits::digits == 53 && + std::numeric_limits::min_exponent == -1021 && + std::numeric_limits::max_exponent == 1024 && + sizeof(double) == sizeof(std::uint64_t), + "Quadrature validation requires IEEE 754 binary64 doubles"); +static_assert( + std::bit_cast(1.0) == 0x3ff0000000000000ULL && + std::bit_cast(-0.0) == 0x8000000000000000ULL && + std::bit_cast( + std::numeric_limits::denorm_min()) == 1u, + "Quadrature validation requires the standard binary64 bit layout"); + +/** + * @brief Exact, order-independent sum of finite binary64 values. + * + * Every finite binary64 value is an integer multiple of 2^-1074. Positive and + * negative coefficients are accumulated separately in fixed-width unsigned + * integers, so neither cancellation nor an intermediate floating-point + * overflow can affect the result. + */ +class ExactBinary64Sum { +public: + /** + * @brief Add one value exactly. + * @return False when @p value is infinite or NaN. + */ + bool add(double value) noexcept + { + const std::uint64_t bits = std::bit_cast(value); + const std::uint64_t exponent = + (bits >> fraction_bit_count) & exponent_mask; + if (exponent == exponent_mask) { + return false; + } + + const std::uint64_t fraction = bits & fraction_mask; + const std::uint64_t coefficient = + exponent == 0u ? fraction : hidden_bit | fraction; + const std::size_t shift = + exponent == 0u ? 0u : static_cast(exponent - 1u); + + auto& magnitude = + (bits & sign_mask) == 0u ? positive_ : negative_; + add_shifted(magnitude, coefficient, shift); + return true; + } + + /** + * @brief Test the exact sum against a binary64 target and error budget. + */ + bool is_within(double expected, double error_budget) const noexcept + { + if (!(error_budget >= 0.0)) { + return false; + } + + ExactBinary64Sum residual = *this; + if (!residual.add(-expected)) { + return false; + } + + ExactBinary64Sum budget; + if (!budget.add(error_budget)) { + return false; + } + + return compare( + residual.absolute_magnitude(), + budget.absolute_magnitude()) <= 0; + } + +private: + static constexpr std::size_t limb_bit_count = + std::numeric_limits::digits; + static constexpr std::size_t fraction_bit_count = 52u; + static constexpr std::uint64_t hidden_bit = + std::uint64_t{1} << fraction_bit_count; + static constexpr std::uint64_t fraction_mask = hidden_bit - 1u; + static constexpr std::uint64_t exponent_mask = 0x7ffu; + static constexpr std::uint64_t sign_mask = + std::uint64_t{1} << (limb_bit_count - 1u); + + // A finite binary64 coefficient occupies at most this many bits when + // scaled by 2^1074. The size_t term covers every possible vector length + // plus the expected moment without dynamic allocation. + static constexpr std::size_t maximum_scaled_value_bits = + static_cast( + std::numeric_limits::max_exponent - + std::numeric_limits::min_exponent + + std::numeric_limits::digits); + static constexpr std::size_t accumulator_bit_count = + maximum_scaled_value_bits + + std::numeric_limits::digits; + static constexpr std::size_t accumulator_limb_count = + (accumulator_bit_count + limb_bit_count - 1u) / limb_bit_count; + + using Magnitude = + std::array; + + static void add_word( + Magnitude& magnitude, + std::size_t index, + std::uint64_t word) noexcept + { + while (word != 0u && index < magnitude.size()) { + const std::uint64_t previous = magnitude[index]; + magnitude[index] += word; + word = magnitude[index] < previous ? 1u : 0u; + ++index; + } + assert(word == 0u); + } + + static void add_shifted( + Magnitude& magnitude, + std::uint64_t coefficient, + std::size_t shift) noexcept + { + if (coefficient == 0u) { + return; + } + + const std::size_t index = shift / limb_bit_count; + const std::size_t offset = shift % limb_bit_count; + add_word(magnitude, index, coefficient << offset); + if (offset != 0u) { + add_word( + magnitude, + index + 1u, + coefficient >> (limb_bit_count - offset)); + } + } + + static int compare( + const Magnitude& left, + const Magnitude& right) noexcept + { + for (std::size_t index = left.size(); index-- > 0u;) { + if (left[index] < right[index]) { + return -1; + } + if (left[index] > right[index]) { + return 1; + } + } + return 0; + } + + static Magnitude subtract( + const Magnitude& larger, + const Magnitude& smaller) noexcept + { + Magnitude difference{}; + std::uint64_t borrow = 0u; + for (std::size_t index = 0u; index < difference.size(); ++index) { + const std::uint64_t after_subtraction = + larger[index] - smaller[index]; + const bool subtraction_borrow = + larger[index] < smaller[index]; + difference[index] = after_subtraction - borrow; + const bool incoming_borrow = after_subtraction < borrow; + borrow = + subtraction_borrow || incoming_borrow ? 1u : 0u; + } + assert(borrow == 0u); + return difference; + } + + Magnitude absolute_magnitude() const noexcept + { + return compare(positive_, negative_) >= 0 + ? subtract(positive_, negative_) + : subtract(negative_, positive_); + } + + Magnitude positive_{}; + Magnitude negative_{}; +}; + struct ValidationResult { static constexpr std::size_t no_sample = std::numeric_limits::max(); @@ -145,61 +331,20 @@ ValidationResult validate_weights( const std::vector& weights, double zeroth_moment) noexcept { - long double compensated_sum = 0.0L; - long double correction = 0.0L; - long double absolute_sum = 0.0L; + ExactBinary64Sum exact_sum; for (std::size_t sample = 0; sample < weights.size(); ++sample) { - const double weight = weights[sample]; - if (!std::isfinite(weight)) { + if (!exact_sum.add(weights[sample])) { return {"quadrature weight must be finite", sample}; } - - const long double value = static_cast(weight); - const long double next = compensated_sum + value; - if (std::abs(compensated_sum) >= std::abs(value)) { - correction += (compensated_sum - next) + value; - } else { - correction += (value - next) + compensated_sum; - } - compensated_sum = next; - absolute_sum += std::abs(value); - if (!std::isfinite(compensated_sum) || - !std::isfinite(correction) || - !std::isfinite(absolute_sum)) { - return {"weight accumulation is not finite", sample}; - } } - const long double total = compensated_sum + correction; - if (!std::isfinite(total)) { - return {"weight accumulation is not finite", weights.size() - 1u}; - } - - const long double expected = static_cast(zeroth_moment); - const long double scale = std::max(1.0L, std::abs(expected)); - const long double requested_error_budget = - static_cast(moment_validation_tolerance) * scale; - const long double compensated_error = std::abs(total - expected); - - if (compensated_error > requested_error_budget) { + const double scale = std::max(1.0, std::abs(zeroth_moment)); + const double error_budget = moment_validation_tolerance * scale; + if (!exact_sum.is_within(zeroth_moment, error_budget)) { return {"weights do not reproduce the zeroth moment"}; } - // A single compensated accumulator is not exact for unrestricted signed - // data. Use a conservative uncertainty estimate in the precision used by - // validation before accepting the computed moment. - constexpr long double validation_safety_factor = 8.0L; - const long double validation_error_bound = - compensated_error + - validation_safety_factor * - std::numeric_limits::epsilon() * absolute_sum; - if (!std::isfinite(validation_error_bound) || - validation_error_bound > requested_error_budget) { - return { - "weight sum is too ill-conditioned to validate the zeroth moment"}; - } - return {}; } diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.h b/Code/Source/solver/FE/Quadrature/QuadratureRule.h index 517fd74f6..3f3c4277f 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.h +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.h @@ -50,7 +50,7 @@ * not derive new rules or modify rule storage. * * Complete-data construction, reference-cell traits, point-containment checks, - * compensated weight summation, concrete generators, caches, and rule-selection + * exact weight summation, concrete generators, caches, and rule-selection * facilities are module implementation details. * * ## Rule-provider contract @@ -71,8 +71,9 @@ * a separate revalidation step. The constructor rejects unsupported cells, * negative exactness, empty or mismatched storage, non-finite coordinates or * weights, points outside the declared reference cell, and weights whose sum - * does not equal the canonical rule's zeroth moment in compensated arithmetic - * or cannot be established within the validator's numerical precision. + * does not equal the canonical rule's zeroth moment within the scaled moment + * tolerance. The sum of the stored binary64 weights is evaluated exactly and + * independently of their order. * * Structural validation does not require unique points or nonzero, positive * individual weights. It verifies metadata, containment, finiteness, and the @@ -267,8 +268,7 @@ class QuadratureRule { * @throws InvalidArgumentException If the family is unsupported, exactness * is negative, storage is empty or mismatched, a value is non-finite, a point * is outside the reference cell, the weights do not reproduce its zeroth - * moment in compensated arithmetic, or cancellation prevents validating - * that moment reliably. + * moment within the scaled moment tolerance. */ explicit QuadratureRule(svmp::CellFamily family, RuleData data); diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp index 5586b55e3..2d95e3296 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp @@ -9,6 +9,7 @@ #include "FE/Quadrature/QuadratureRule.h" #include "nn.h" +#include #include #include #include @@ -969,19 +970,134 @@ TEST(QuadratureRuleValidation, RejectsIncorrectZerothMomentDespiteLargeCancellat "weights do not reproduce the zeroth moment"); } -TEST(QuadratureRuleValidation, RejectsSignedDataWithUnreliableMoment) +TEST(QuadratureRuleValidation, RejectsIncorrectMomentHiddenByCancellation) { const double medium = std::ldexp(1.0, 93); const double large = std::ldexp(1.0, 233); const double residual = std::ldexp(1.0, 14); - EXPECT_THROW( + expect_invalid_argument_with_message( + [medium, large, residual] { + (void)RuleProbe( + svmp::CellFamily::Line, + 0, + std::vector(6u, QuadPoint::Zero()), + {-medium, -large, residual, medium, large, 2.0}); + }, + "weights do not reproduce the zeroth moment"); +} + +TEST(QuadratureRuleValidation, ExactMomentValidationIsOrderIndependent) +{ + const double maximum = std::numeric_limits::max(); + std::vector valid_weights{ + -maximum, -maximum, 2.0, maximum, maximum}; + std::sort(valid_weights.begin(), valid_weights.end()); + + std::size_t permutation_count = 0u; + do { + SCOPED_TRACE(permutation_count); + EXPECT_NO_THROW( + (void)RuleProbe( + svmp::CellFamily::Line, + 0, + std::vector( + valid_weights.size(), + QuadPoint::Zero()), + valid_weights)); + ++permutation_count; + } while (std::next_permutation( + valid_weights.begin(), + valid_weights.end())); + EXPECT_EQ(permutation_count, 30u); + + const double excessive_residual = std::ldexp(1.0, -38); + std::vector invalid_weights{ + -maximum, + -maximum, + 2.0, + excessive_residual, + maximum, + maximum}; + std::sort(invalid_weights.begin(), invalid_weights.end()); + + permutation_count = 0u; + do { + SCOPED_TRACE(permutation_count); + expect_invalid_argument_with_message( + [&invalid_weights] { + (void)RuleProbe( + svmp::CellFamily::Line, + 0, + std::vector( + invalid_weights.size(), + QuadPoint::Zero()), + invalid_weights); + }, + "weights do not reproduce the zeroth moment"); + ++permutation_count; + } while (std::next_permutation( + invalid_weights.begin(), + invalid_weights.end())); + EXPECT_EQ(permutation_count, 180u); +} + +TEST(QuadratureRuleValidation, HandlesExtremeAndSubnormalCancellation) +{ + const double maximum = std::numeric_limits::max(); + const double minimum_normal = std::numeric_limits::min(); + const double minimum_subnormal = + std::numeric_limits::denorm_min(); + const double next_normal = + std::nextafter(minimum_normal, std::numeric_limits::infinity()); + + EXPECT_NO_THROW( (void)RuleProbe( svmp::CellFamily::Line, 0, std::vector(6u, QuadPoint::Zero()), - {-medium, -large, residual, medium, large, 2.0}), - InvalidArgumentException); + { + maximum, + minimum_normal, + minimum_subnormal, + 2.0, + -next_normal, + -maximum, + })); + + expect_invalid_argument_with_message( + [maximum] { + (void)RuleProbe( + svmp::CellFamily::Line, + 0, + std::vector(2u, QuadPoint::Zero()), + {maximum, maximum}); + }, + "weights do not reproduce the zeroth moment"); +} + +TEST(QuadratureRuleValidation, AppliesToleranceToTheExactWeightSum) +{ + const double maximum = std::numeric_limits::max(); + const double accepted_residual = std::ldexp(1.0, -42); + const double rejected_residual = std::ldexp(1.0, -38); + + EXPECT_NO_THROW( + (void)RuleProbe( + svmp::CellFamily::Line, + 0, + std::vector(4u, QuadPoint::Zero()), + {maximum, accepted_residual, -maximum, 2.0})); + + expect_invalid_argument_with_message( + [maximum, rejected_residual] { + (void)RuleProbe( + svmp::CellFamily::Line, + 0, + std::vector(4u, QuadPoint::Zero()), + {maximum, rejected_residual, -maximum, 2.0}); + }, + "weights do not reproduce the zeroth moment"); } TEST(QuadratureRuleValidation, AcceptsLargePositiveRule) From ce8ca0781a0cb99c09e7de09b78d609e6da23e43 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 23 Jul 2026 16:36:59 -0700 Subject: [PATCH 16/21] fixing spacing for exact binary summation --- Code/Source/solver/FE/Quadrature/QuadratureRule.cpp | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp index d4db4ec64..5f345b6f9 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp @@ -109,15 +109,12 @@ class ExactBinary64Sum { } private: - static constexpr std::size_t limb_bit_count = - std::numeric_limits::digits; + static constexpr std::size_t limb_bit_count = std::numeric_limits::digits; static constexpr std::size_t fraction_bit_count = 52u; - static constexpr std::uint64_t hidden_bit = - std::uint64_t{1} << fraction_bit_count; + static constexpr std::uint64_t hidden_bit = std::uint64_t{1} << fraction_bit_count; static constexpr std::uint64_t fraction_mask = hidden_bit - 1u; static constexpr std::uint64_t exponent_mask = 0x7ffu; - static constexpr std::uint64_t sign_mask = - std::uint64_t{1} << (limb_bit_count - 1u); + static constexpr std::uint64_t sign_mask = std::uint64_t{1} << (limb_bit_count - 1u); // A finite binary64 coefficient occupies at most this many bits when // scaled by 2^1074. The size_t term covers every possible vector length From d66e30ccdba30676dca042624ef086496c012472 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Thu, 23 Jul 2026 17:20:04 -0700 Subject: [PATCH 17/21] Document binary64 validation assumptions --- Code/Source/solver/FE/Quadrature/QuadratureRule.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp index 5f345b6f9..e506dc8fa 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp @@ -34,6 +34,8 @@ struct ReferenceCellTraits { constexpr double coordinate_validation_tolerance = 1.0e-12; constexpr double moment_validation_tolerance = 1.0e-12; +// Exact summation decodes double bit patterns directly, so fail at compile +// time unless both the binary64 value model and object layout are supported. static_assert( std::numeric_limits::is_iec559 && std::numeric_limits::radix == 2 && From e8d368046bd6943a9e28d2816634669dcd4bc914 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Fri, 24 Jul 2026 11:37:04 -0700 Subject: [PATCH 18/21] Qualify quadrature vector include --- Code/Source/solver/FE/Quadrature/QuadratureRule.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.h b/Code/Source/solver/FE/Quadrature/QuadratureRule.h index 3f3c4277f..965b19914 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.h +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.h @@ -114,7 +114,7 @@ */ #include "FE/Common/Types.h" -#include "Math/Vector.h" +#include "FE/Math/Vector.h" #include #include From 42c9e9ec966ab21b1a56b91f879d4c1d69570a94 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Fri, 24 Jul 2026 12:00:37 -0700 Subject: [PATCH 19/21] Use reference cell measure terminology --- .../solver/FE/Quadrature/QuadratureRule.cpp | 20 ++++----- .../solver/FE/Quadrature/QuadratureRule.h | 30 ++++++------- .../FE/Quadrature/test_QuadratureRules.cpp | 42 ++++++++++--------- 3 files changed, 48 insertions(+), 44 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp index e506dc8fa..9a07c41fb 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp @@ -28,11 +28,11 @@ namespace { struct ReferenceCellTraits { int dimension; - double zeroth_moment; + double reference_cell_measure; }; constexpr double coordinate_validation_tolerance = 1.0e-12; -constexpr double moment_validation_tolerance = 1.0e-12; +constexpr double measure_validation_tolerance = 1.0e-12; // Exact summation decodes double bit patterns directly, so fail at compile // time unless both the binary64 value model and object layout are supported. @@ -328,7 +328,7 @@ ValidationResult validate_point( ValidationResult validate_weights( const std::vector& weights, - double zeroth_moment) noexcept + double reference_cell_measure) noexcept { ExactBinary64Sum exact_sum; @@ -338,10 +338,10 @@ ValidationResult validate_weights( } } - const double scale = std::max(1.0, std::abs(zeroth_moment)); - const double error_budget = moment_validation_tolerance * scale; - if (!exact_sum.is_within(zeroth_moment, error_budget)) { - return {"weights do not reproduce the zeroth moment"}; + const double scale = std::max(1.0, std::abs(reference_cell_measure)); + const double error_budget = measure_validation_tolerance * scale; + if (!exact_sum.is_within(reference_cell_measure, error_budget)) { + return {"weights do not reproduce the reference-cell measure"}; } return {}; @@ -372,7 +372,7 @@ ValidationResult validate_rule_data( } } - return validate_weights(weights, traits.zeroth_moment); + return validate_weights(weights, traits.reference_cell_measure); } std::string validation_failure_message(const ValidationResult& result) @@ -403,7 +403,7 @@ QuadratureRule::QuadratureRule(ValidatedState state) : cell_family_(state.cell_family), dimension_(state.dimension), polynomial_exactness_(state.polynomial_exactness), - zeroth_moment_(state.zeroth_moment), + reference_cell_measure_(state.reference_cell_measure), points_(std::move(state.points)), weights_(std::move(state.weights)) { @@ -435,7 +435,7 @@ QuadratureRule::ValidatedState QuadratureRule::validate( family, traits->dimension, data.polynomial_exactness, - traits->zeroth_moment, + traits->reference_cell_measure, std::move(data.points), std::move(data.weights), }; diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.h b/Code/Source/solver/FE/Quadrature/QuadratureRule.h index 965b19914..72e28c4f0 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.h +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.h @@ -28,9 +28,9 @@ * \int_{\hat K} f(\hat x)\,d\hat x * \approx \sum_q w_q f(\hat x_q). * @f] - * The zeroth moment is the analytic normalization + * The reference-cell measure satisfies * @f[ - * M_0 = \int_{\hat K} 1\,d\hat x = \sum_q w_q. + * |\hat K| = \int_{\hat K} 1\,d\hat x = \sum_q w_q. * @f] * A rule identifies its reference-cell family, reports its dimension and * declared polynomial exactness, and keeps every point paired with its @@ -71,13 +71,13 @@ * a separate revalidation step. The constructor rejects unsupported cells, * negative exactness, empty or mismatched storage, non-finite coordinates or * weights, points outside the declared reference cell, and weights whose sum - * does not equal the canonical rule's zeroth moment within the scaled moment - * tolerance. The sum of the stored binary64 weights is evaluated exactly and - * independently of their order. + * does not equal the canonical reference-cell measure within the scaled + * measure tolerance. The sum of the stored binary64 weights is evaluated + * exactly and independently of their order. * * Structural validation does not require unique points or nonzero, positive * individual weights. It verifies metadata, containment, finiteness, and the - * zeroth moment; it does not prove higher-order polynomial moments. A + * reference-cell measure; it does not prove higher-order polynomial moments. A * polynomial exactness of @f$p@f$ guarantees every polynomial of total degree * at most @f$p@f$. A rule can integrate selected higher-degree polynomials * without increasing that common guarantee. Every concrete provider is @@ -90,7 +90,7 @@ * are active, and every inactive component is zero within the coordinate * tolerance. The supported canonical domains are: * - * | Cell family | Canonical reference domain | Zeroth moment | + * | Cell family | Canonical reference domain | Reference-cell measure | * | ----------- | -------------------------- | ------------- | * | Point | @f$(0,0,0)@f$ | @f$1@f$ | * | Line | @f$[-1,1]@f$ | @f$2@f$ | @@ -227,19 +227,19 @@ class QuadratureRule { const std::vector& weights() const noexcept { return weights_; } /** - * @brief Return the rule's zeroth moment. + * @brief Return the measure of the canonical reference cell. * * This is the integral of the constant function one. All supported rules * are unweighted rules on complete canonical reference cells, so the * constructor derives this value from cell_family() and it equals the * geometric measure of that cell. * @f[ - * M_0 = \int_{\hat K} 1\,d\hat x = \sum_q w_q. + * |\hat K| = \int_{\hat K} 1\,d\hat x = \sum_q w_q. * @f] * - * @return Expected sum of the quadrature weights. + * @return Geometric measure of the canonical reference cell. */ - double zeroth_moment() const noexcept { return zeroth_moment_; } + double reference_cell_measure() const noexcept { return reference_cell_measure_; } protected: /** @@ -258,8 +258,8 @@ class QuadratureRule { /** * @brief Construct and validate one complete immutable rule. * - * Dimension and zeroth moment are derived from @p family; callers cannot - * supply redundant topology metadata. + * Dimension and reference-cell measure are derived from @p family; callers + * cannot supply redundant topology metadata. * * @param family Supported canonical reference-cell family. * @param data Complete exactness, point, and weight payload. @@ -278,7 +278,7 @@ class QuadratureRule { svmp::CellFamily cell_family; int dimension; int polynomial_exactness; - double zeroth_moment; + double reference_cell_measure; std::vector points; std::vector weights; }; @@ -292,7 +292,7 @@ class QuadratureRule { const svmp::CellFamily cell_family_; ///< Canonical reference topology. const int dimension_; ///< Number of active coordinate components. const int polynomial_exactness_; ///< Exactness declared by the concrete generator. - const double zeroth_moment_; ///< Canonical reference-cell measure. + const double reference_cell_measure_; ///< Canonical reference-cell measure. const std::vector points_; ///< Ordered immutable reference coordinates. const std::vector weights_; ///< Immutable weights paired with points_. }; diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp index 2d95e3296..7208cd25a 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp @@ -59,7 +59,7 @@ struct SolverRuleCase { int dimension; int requested_exactness; int advertised_exactness; - double zeroth_moment; + double reference_cell_measure; ExpectedSamples samples; }; @@ -563,13 +563,13 @@ TEST(QuadraturePhase01Baseline, StandardSelectionTableHasOrderedPointWeightData) EXPECT_FALSE(c.samples.empty()); EXPECT_GT(c.requested_exactness, 0); EXPECT_GE(c.advertised_exactness, c.requested_exactness); - EXPECT_GT(c.zeroth_moment, 0.0); + EXPECT_GT(c.reference_cell_measure, 0.0); double weight_sum = 0.0; for (const auto& sample : c.samples) { weight_sum += sample.weight; } - EXPECT_NEAR(weight_sum, c.zeroth_moment, kTol); + EXPECT_NEAR(weight_sum, c.reference_cell_measure, kTol); } } @@ -583,7 +583,9 @@ TEST(QuadraturePhase01Baseline, CanonicalFixturesSatisfyTheRuleAndExactnessContr EXPECT_EQ(rule.dimension(), c.dimension); EXPECT_EQ(rule.polynomial_exactness(), c.advertised_exactness); EXPECT_EQ(rule.num_points(), c.samples.size()); - EXPECT_DOUBLE_EQ(rule.zeroth_moment(), c.zeroth_moment); + EXPECT_DOUBLE_EQ( + rule.reference_cell_measure(), + c.reference_cell_measure); expect_samples_in_order(rule, c.samples); expect_total_degree_exact(rule, c.advertised_exactness); } @@ -826,7 +828,9 @@ TEST(QuadratureRuleValidation, AcceptsEverySupportedReferenceCell) {{c.point[0], c.point[1], c.point[2]}}, {c.expected_measure}); EXPECT_EQ(rule.dimension(), c.expected_dimension); - EXPECT_DOUBLE_EQ(rule.zeroth_moment(), c.expected_measure); + EXPECT_DOUBLE_EQ( + rule.reference_cell_measure(), + c.expected_measure); } } @@ -932,7 +936,7 @@ TEST(QuadratureRuleValidation, RejectsInactiveCoordinatesAndOutOfCellPoints) InvalidArgumentException); } -TEST(QuadratureRuleValidation, EnforcesZerothMomentButAllowsNegativeWeights) +TEST(QuadratureRuleValidation, EnforcesReferenceCellMeasureButAllowsNegativeWeights) { expect_invalid_argument_with_message( [] { @@ -942,7 +946,7 @@ TEST(QuadratureRuleValidation, EnforcesZerothMomentButAllowsNegativeWeights) {{0.25, 0.25, 0.0}}, {1.0}); }, - "weights do not reproduce the zeroth moment"); + "weights do not reproduce the reference-cell measure"); const RuleProbe rule( svmp::CellFamily::Triangle, @@ -952,10 +956,10 @@ TEST(QuadratureRuleValidation, EnforcesZerothMomentButAllowsNegativeWeights) EXPECT_LT(rule.weight(0), 0.0); EXPECT_DOUBLE_EQ( integrate_monomial(rule, 0, 0, 0), - rule.zeroth_moment()); + rule.reference_cell_measure()); } -TEST(QuadratureRuleValidation, RejectsIncorrectZerothMomentDespiteLargeCancellation) +TEST(QuadratureRuleValidation, RejectsIncorrectMeasureDespiteLargeCancellation) { expect_invalid_argument_with_message( [] { @@ -967,10 +971,10 @@ TEST(QuadratureRuleValidation, RejectsIncorrectZerothMomentDespiteLargeCancellat {0.5, 0.0, 0.0}}, {1.0e20, 0.0, -1.0e20}); }, - "weights do not reproduce the zeroth moment"); + "weights do not reproduce the reference-cell measure"); } -TEST(QuadratureRuleValidation, RejectsIncorrectMomentHiddenByCancellation) +TEST(QuadratureRuleValidation, RejectsIncorrectMeasureHiddenByCancellation) { const double medium = std::ldexp(1.0, 93); const double large = std::ldexp(1.0, 233); @@ -984,10 +988,10 @@ TEST(QuadratureRuleValidation, RejectsIncorrectMomentHiddenByCancellation) std::vector(6u, QuadPoint::Zero()), {-medium, -large, residual, medium, large, 2.0}); }, - "weights do not reproduce the zeroth moment"); + "weights do not reproduce the reference-cell measure"); } -TEST(QuadratureRuleValidation, ExactMomentValidationIsOrderIndependent) +TEST(QuadratureRuleValidation, ExactMeasureValidationIsOrderIndependent) { const double maximum = std::numeric_limits::max(); std::vector valid_weights{ @@ -1034,7 +1038,7 @@ TEST(QuadratureRuleValidation, ExactMomentValidationIsOrderIndependent) QuadPoint::Zero()), invalid_weights); }, - "weights do not reproduce the zeroth moment"); + "weights do not reproduce the reference-cell measure"); ++permutation_count; } while (std::next_permutation( invalid_weights.begin(), @@ -1073,7 +1077,7 @@ TEST(QuadratureRuleValidation, HandlesExtremeAndSubnormalCancellation) std::vector(2u, QuadPoint::Zero()), {maximum, maximum}); }, - "weights do not reproduce the zeroth moment"); + "weights do not reproduce the reference-cell measure"); } TEST(QuadratureRuleValidation, AppliesToleranceToTheExactWeightSum) @@ -1097,7 +1101,7 @@ TEST(QuadratureRuleValidation, AppliesToleranceToTheExactWeightSum) std::vector(4u, QuadPoint::Zero()), {maximum, rejected_residual, -maximum, 2.0}); }, - "weights do not reproduce the zeroth moment"); + "weights do not reproduce the reference-cell measure"); } TEST(QuadratureRuleValidation, AcceptsLargePositiveRule) @@ -1137,7 +1141,7 @@ TEST(QuadratureRuleValidation, AcceptsLargeSignedRule) EXPECT_LT(rule.weight(0), 0.0); EXPECT_NEAR( integrate_monomial(rule, 0, 0, 0), - rule.zeroth_moment(), + rule.reference_cell_measure(), kTol); } @@ -1182,7 +1186,7 @@ TEST(QuadratureRuleValidation, AppliesConstructionWeightTolerance) {{0.25, 0.25, 0.0}}, {0.5 + rejected_offset}); }, - "weights do not reproduce the zeroth moment"); + "weights do not reproduce the reference-cell measure"); } TEST(QuadratureRuleContract, PublishesOnlyACompleteImmutableQueryInterface) @@ -1219,7 +1223,7 @@ TEST(QuadratureRuleContract, PublishesOnlyACompleteImmutableQueryInterface) EXPECT_EQ(rule.cell_family(), svmp::CellFamily::Line); EXPECT_EQ(rule.dimension(), 1); EXPECT_EQ(rule.polynomial_exactness(), 3); - EXPECT_DOUBLE_EQ(rule.zeroth_moment(), 2.0); + EXPECT_DOUBLE_EQ(rule.reference_cell_measure(), 2.0); ASSERT_EQ(rule.num_points(), 2u); ASSERT_EQ(rule.points().size(), 2u); ASSERT_EQ(rule.weights().size(), 2u); From 4feee7874ceaefc4ffebddeded803cbf11864f17 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Fri, 24 Jul 2026 12:17:08 -0700 Subject: [PATCH 20/21] Complete reference cell measure documentation --- Code/Source/solver/FE/Quadrature/QuadratureRule.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.h b/Code/Source/solver/FE/Quadrature/QuadratureRule.h index 72e28c4f0..3cdb42c2a 100644 --- a/Code/Source/solver/FE/Quadrature/QuadratureRule.h +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.h @@ -267,8 +267,8 @@ class QuadratureRule { * when every other rule invariant is satisfied. * @throws InvalidArgumentException If the family is unsupported, exactness * is negative, storage is empty or mismatched, a value is non-finite, a point - * is outside the reference cell, the weights do not reproduce its zeroth - * moment within the scaled moment tolerance. + * is outside the reference cell, or the weights do not reproduce the + * reference-cell measure within the scaled measure tolerance. */ explicit QuadratureRule(svmp::CellFamily family, RuleData data); From 6fb1dad649706d6fdac67383312095b80aa4ded6 Mon Sep 17 00:00:00 2001 From: Zachary Sexton Date: Fri, 24 Jul 2026 12:48:04 -0700 Subject: [PATCH 21/21] Focus quadrature tests on new infrastructure --- .../FE/Quadrature/test_QuadratureRules.cpp | 710 +----------------- 1 file changed, 8 insertions(+), 702 deletions(-) diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp index 7208cd25a..f0f3465c6 100644 --- a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp +++ b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp @@ -1,13 +1,12 @@ /** * @file test_QuadratureRules.cpp - * @brief Baseline solver data and core quadrature rule contract tests. + * @brief Unit tests for the core quadrature rule infrastructure. */ #include #include "FE/Common/FEException.h" #include "FE/Quadrature/QuadratureRule.h" -#include "nn.h" #include #include @@ -41,28 +40,9 @@ class RuleProbe final : public QuadratureRule { }; constexpr double kTol = 1.0e-12; -constexpr double kExactnessTol = 3.0e-12; using ExpectedPoint = std::array; -struct ExpectedSample { - ExpectedPoint point; - double weight; -}; - -using ExpectedSamples = std::vector; - -struct SolverRuleCase { - const char* name; - consts::ElementType solver_type; - ElementType fe_type; - int dimension; - int requested_exactness; - int advertised_exactness; - double reference_cell_measure; - ExpectedSamples samples; -}; - template void expect_invalid_argument_with_message( Function&& function, @@ -84,691 +64,17 @@ void expect_invalid_argument_with_message( } } -ExpectedSamples point_samples() -{ - return {{{0.0, 0.0, 0.0}, 1.0}}; -} - -ExpectedSamples line2_samples() -{ - const double a = 1.0 / std::sqrt(3.0); - return { - {{-a, 0.0, 0.0}, 1.0}, - {{ a, 0.0, 0.0}, 1.0}, - }; -} - -ExpectedSamples line3_samples() -{ - const double a = std::sqrt(0.6); - return { - {{-a, 0.0, 0.0}, 5.0 / 9.0}, - {{ a, 0.0, 0.0}, 5.0 / 9.0}, - {{0.0, 0.0, 0.0}, 8.0 / 9.0}, - }; -} - -ExpectedSamples triangle3_samples() -{ - const double s = 2.0 / 3.0; - const double t = 1.0 / 6.0; - return { - {{t, t, 0.0}, 1.0 / 6.0}, - {{s, t, 0.0}, 1.0 / 6.0}, - {{t, s, 0.0}, 1.0 / 6.0}, - }; -} - -ExpectedSamples triangle6_samples() -{ - const double centroid = 0.333333333333333; - const double a1 = 0.797426985353087; - const double b1 = 0.101286507323456; - const double a2 = 0.059715871789770; - const double b2 = 0.470142064105115; - const double w0 = 0.225000000000000 * 0.5; - const double w1 = 0.125939180544827 * 0.5; - const double w2 = 0.132394152788506 * 0.5; - return { - {{centroid, centroid, 0.0}, w0}, - {{a1, b1, 0.0}, w1}, - {{b1, a1, 0.0}, w1}, - {{b1, b1, 0.0}, w1}, - {{a2, b2, 0.0}, w2}, - {{b2, a2, 0.0}, w2}, - {{b2, b2, 0.0}, w2}, - }; -} - -ExpectedSamples quad4_samples() -{ - const double a = 1.0 / std::sqrt(3.0); - return { - {{-a, -a, 0.0}, 1.0}, - {{ a, -a, 0.0}, 1.0}, - {{ a, a, 0.0}, 1.0}, - {{-a, a, 0.0}, 1.0}, - }; -} - -ExpectedSamples quad9_samples(double a = std::sqrt(0.6)) -{ - const double corner_weight = 25.0 / 81.0; - const double edge_weight = 40.0 / 81.0; - return { - {{-a, -a, 0.0}, corner_weight}, - {{ a, -a, 0.0}, corner_weight}, - {{ a, a, 0.0}, corner_weight}, - {{-a, a, 0.0}, corner_weight}, - {{0.0, -a, 0.0}, edge_weight}, - {{ a, 0.0, 0.0}, edge_weight}, - {{0.0, a, 0.0}, edge_weight}, - {{-a, 0.0, 0.0}, edge_weight}, - {{0.0, 0.0, 0.0}, 64.0 / 81.0}, - }; -} - -ExpectedSamples tetra4_samples() -{ - const double a = (5.0 + 3.0 * std::sqrt(5.0)) / 20.0; - const double b = (1.0 - a) / 3.0; - return { - {{a, b, b}, 1.0 / 24.0}, - {{b, a, b}, 1.0 / 24.0}, - {{b, b, a}, 1.0 / 24.0}, - {{b, b, b}, 1.0 / 24.0}, - }; -} - -ExpectedSamples tetra10_samples() -{ - const double one_third = 0.3333333333333330; - const double a1 = 0.0909090909090910; - const double b1 = 0.7272727272727270; - const double a2 = 0.0665501535736640; - const double b2 = 0.4334498464263360; - const double w0 = 0.0302836780970890; - const double w1 = 0.0060267857142860; - const double w2 = 0.0116452490860290; - const double w3 = 0.0109491415613860; - return { - {{0.25, 0.25, 0.25}, w0}, - {{0.0, one_third, one_third}, w1}, - {{one_third, 0.0, one_third}, w1}, - {{one_third, one_third, 0.0}, w1}, - {{one_third, one_third, one_third}, w1}, - {{b1, a1, a1}, w2}, - {{a1, b1, a1}, w2}, - {{a1, a1, b1}, w2}, - {{a1, a1, a1}, w2}, - {{a2, a2, b2}, w3}, - {{a2, b2, a2}, w3}, - {{a2, b2, b2}, w3}, - {{b2, b2, a2}, w3}, - {{b2, a2, b2}, w3}, - {{b2, a2, a2}, w3}, - }; -} - -ExpectedSamples hex8_samples() -{ - const double a = 1.0 / std::sqrt(3.0); - return { - {{-a, -a, -a}, 1.0}, - {{ a, -a, -a}, 1.0}, - {{ a, a, -a}, 1.0}, - {{-a, a, -a}, 1.0}, - {{-a, -a, a}, 1.0}, - {{ a, -a, a}, 1.0}, - {{ a, a, a}, 1.0}, - {{-a, a, a}, 1.0}, - }; -} - -ExpectedSamples hex27_samples() -{ - const double a = std::sqrt(0.6); - const double w0 = 125.0 / 729.0; - const double w1 = 200.0 / 729.0; - const double w2 = 320.0 / 729.0; - return { - {{-a, -a, -a}, w0}, {{ a, -a, -a}, w0}, - {{ a, a, -a}, w0}, {{-a, a, -a}, w0}, - {{-a, -a, a}, w0}, {{ a, -a, a}, w0}, - {{ a, a, a}, w0}, {{-a, a, a}, w0}, - - {{0.0, -a, -a}, w1}, {{ a, 0.0, -a}, w1}, - {{0.0, a, -a}, w1}, {{-a, 0.0, -a}, w1}, - {{0.0, -a, a}, w1}, {{ a, 0.0, a}, w1}, - {{0.0, a, a}, w1}, {{-a, 0.0, a}, w1}, - {{-a, -a, 0.0}, w1}, {{ a, -a, 0.0}, w1}, - {{ a, a, 0.0}, w1}, {{-a, a, 0.0}, w1}, - - {{-a, 0.0, 0.0}, w2}, {{ a, 0.0, 0.0}, w2}, - {{0.0, -a, 0.0}, w2}, {{0.0, a, 0.0}, w2}, - {{0.0, 0.0, -a}, w2}, {{0.0, 0.0, a}, w2}, - {{0.0, 0.0, 0.0}, 512.0 / 729.0}, - }; -} - -ExpectedSamples wedge6_samples() -{ - const double s = 2.0 / 3.0; - const double t = 1.0 / 6.0; - const double z = 1.0 / std::sqrt(3.0); - return { - {{s, t, -z}, 1.0 / 6.0}, - {{t, s, -z}, 1.0 / 6.0}, - {{t, t, -z}, 1.0 / 6.0}, - {{s, t, z}, 1.0 / 6.0}, - {{t, s, z}, 1.0 / 6.0}, - {{t, t, z}, 1.0 / 6.0}, - }; -} - -const std::vector& standard_solver_cases() -{ - static const std::vector cases = { - {"PNT", consts::ElementType::PNT, ElementType::Point1, 0, 1, 1, 1.0, point_samples()}, - {"LIN1", consts::ElementType::LIN1, ElementType::Line2, 1, 2, 3, 2.0, line2_samples()}, - {"LIN2", consts::ElementType::LIN2, ElementType::Line3, 1, 4, 5, 2.0, line3_samples()}, - {"TRI3", consts::ElementType::TRI3, ElementType::Triangle3, 2, 2, 2, 0.5, triangle3_samples()}, - {"TRI6", consts::ElementType::TRI6, ElementType::Triangle6, 2, 5, 5, 0.5, triangle6_samples()}, - {"QUD4", consts::ElementType::QUD4, ElementType::Quad4, 2, 2, 3, 4.0, quad4_samples()}, - {"QUD8", consts::ElementType::QUD8, ElementType::Quad8, 2, 4, 5, 4.0, quad9_samples()}, - {"QUD9", consts::ElementType::QUD9, ElementType::Quad9, 2, 4, 5, 4.0, quad9_samples()}, - {"TET4", consts::ElementType::TET4, ElementType::Tetra4, 3, 2, 2, 1.0 / 6.0, tetra4_samples()}, - {"TET10", consts::ElementType::TET10, ElementType::Tetra10, 3, 5, 5, 1.0 / 6.0, tetra10_samples()}, - {"HEX8", consts::ElementType::HEX8, ElementType::Hex8, 3, 2, 3, 8.0, hex8_samples()}, - {"HEX20", consts::ElementType::HEX20, ElementType::Hex20, 3, 4, 5, 8.0, hex27_samples()}, - {"HEX27", consts::ElementType::HEX27, ElementType::Hex27, 3, 4, 5, 8.0, hex27_samples()}, - {"WDG", consts::ElementType::WDG, ElementType::Wedge6, 3, 2, 2, 1.0, wedge6_samples()}, - }; - return cases; -} - -std::vector rule_points(const ExpectedSamples& samples) -{ - std::vector points; - points.reserve(samples.size()); - for (const auto& sample : samples) { - points.push_back({ - sample.point[0], - sample.point[1], - sample.point[2], - }); - } - return points; -} - -std::vector rule_weights(const ExpectedSamples& samples) -{ - std::vector weights; - weights.reserve(samples.size()); - for (const auto& sample : samples) { - weights.push_back(sample.weight); - } - return weights; -} - -std::vector rule_points(const Array& points, int dimension) -{ - std::vector result; - result.reserve(static_cast(points.ncols())); - for (int q = 0; q < points.ncols(); ++q) { - QuadPoint point = QuadPoint::Zero(); - for (int component = 0; component < dimension; ++component) { - point[static_cast(component)] = points(component, q); - } - result.push_back(point); - } - return result; -} - -std::vector rule_weights(const Vector& weights) -{ - std::vector result; - result.reserve(static_cast(weights.size())); - for (int q = 0; q < weights.size(); ++q) { - result.push_back(weights[q]); - } - return result; -} - -double factorial(int n) -{ - double value = 1.0; - for (int factor = 2; factor <= n; ++factor) { - value *= static_cast(factor); - } - return value; -} - -double interval_moment(int power) -{ - return power % 2 == 0 - ? 2.0 / static_cast(power + 1) - : 0.0; -} - -double simplex_moment(int px, int py, int pz, int dimension) -{ - const int total = px + - (dimension >= 2 ? py : 0) + - (dimension >= 3 ? pz : 0); - double numerator = factorial(px); - if (dimension >= 2) { - numerator *= factorial(py); - } - if (dimension >= 3) { - numerator *= factorial(pz); - } - return numerator / factorial(total + dimension); -} - -double exact_reference_moment( - svmp::CellFamily family, - int px, - int py, - int pz) -{ - switch (family) { - case svmp::CellFamily::Point: - return px == 0 && py == 0 && pz == 0 ? 1.0 : 0.0; - case svmp::CellFamily::Line: - return interval_moment(px); - case svmp::CellFamily::Triangle: - return simplex_moment(px, py, 0, 2); - case svmp::CellFamily::Quad: - return interval_moment(px) * interval_moment(py); - case svmp::CellFamily::Tetra: - return simplex_moment(px, py, pz, 3); - case svmp::CellFamily::Hex: - return interval_moment(px) * - interval_moment(py) * - interval_moment(pz); - case svmp::CellFamily::Wedge: - return simplex_moment(px, py, 0, 2) * interval_moment(pz); - default: - svmp::raise( - "Unsupported reference cell in exactness test"); - } -} - -double integer_power(double base, int exponent) -{ - double value = 1.0; - for (int factor = 0; factor < exponent; ++factor) { - value *= base; - } - return value; -} - -double integrate_monomial( - const QuadratureRule& rule, - int px, - int py, - int pz) -{ - double result = 0.0; - for (std::size_t q = 0; q < rule.num_points(); ++q) { - const auto point = rule.point(q); - result += rule.weight(q) * - integer_power(point[0], px) * - integer_power(point[1], py) * - integer_power(point[2], pz); - } - return result; -} - -void expect_total_degree_exact( - const QuadratureRule& rule, - int degree, - double tolerance = kExactnessTol) -{ - const int max_px = rule.dimension() >= 1 ? degree : 0; - for (int px = 0; px <= max_px; ++px) { - const int max_py = rule.dimension() >= 2 ? degree - px : 0; - for (int py = 0; py <= max_py; ++py) { - const int max_pz = rule.dimension() >= 3 ? degree - px - py : 0; - for (int pz = 0; pz <= max_pz; ++pz) { - EXPECT_NEAR( - integrate_monomial(rule, px, py, pz), - exact_reference_moment( - rule.cell_family(), px, py, pz), - tolerance) - << "powers=(" << px << ',' << py << ',' << pz << ')'; - } - } - } -} - -void expect_samples_in_order( - const QuadratureRule& rule, - const ExpectedSamples& expected, - double tolerance = kTol) -{ - ASSERT_EQ(rule.num_points(), expected.size()); - for (std::size_t q = 0; q < expected.size(); ++q) { - EXPECT_NEAR(rule.weight(q), expected[q].weight, tolerance) - << "sample=" << q; - const auto point = rule.point(q); - for (std::size_t component = 0; component < 3u; ++component) { - EXPECT_NEAR( - point[component], - expected[q].point[component], - tolerance) - << "sample=" << q << " component=" << component; - } - } -} - -RuleProbe make_rule(const SolverRuleCase& c) -{ - return RuleProbe( - to_mesh_family(c.fe_type), - c.advertised_exactness, - rule_points(c.samples), - rule_weights(c.samples)); -} - -void expect_legacy_arrays_in_order(const Vector& weights, - const Array& points, - int dimension, - const ExpectedSamples& expected, - double tol = kTol) -{ - ASSERT_EQ(weights.size(), static_cast(expected.size())); - ASSERT_EQ(points.nrows(), dimension); - if (dimension > 0) { - ASSERT_EQ(points.ncols(), static_cast(expected.size())); - } - - for (std::size_t q = 0; q < expected.size(); ++q) { - EXPECT_NEAR(weights[static_cast(q)], expected[q].weight, tol) - << "sample=" << q; - for (int d = 0; d < dimension; ++d) { - EXPECT_NEAR(points(d, static_cast(q)), - expected[q].point[static_cast(d)], tol) - << "sample=" << q << " component=" << d; - } - } -} - -bool generic_context_supports(consts::ElementType type) -{ - switch (type) { - case consts::ElementType::LIN1: - case consts::ElementType::TRI3: - case consts::ElementType::QUD4: - case consts::ElementType::QUD9: - case consts::ElementType::TET4: - case consts::ElementType::HEX8: - case consts::ElementType::HEX20: - case consts::ElementType::HEX27: - case consts::ElementType::WDG: - return true; - default: - return false; - } -} - -bool volume_context_supports(consts::ElementType type) -{ - switch (type) { - case consts::ElementType::LIN1: - case consts::ElementType::LIN2: - case consts::ElementType::TRI3: - case consts::ElementType::TRI6: - case consts::ElementType::QUD4: - case consts::ElementType::QUD9: - case consts::ElementType::TET4: - case consts::ElementType::TET10: - case consts::ElementType::HEX8: - case consts::ElementType::HEX20: - case consts::ElementType::HEX27: - case consts::ElementType::WDG: - return true; - default: - return false; - } -} - -bool face_context_supports(consts::ElementType type) +double weight_sum(const QuadratureRule& rule) { - switch (type) { - case consts::ElementType::PNT: - case consts::ElementType::LIN1: - case consts::ElementType::LIN2: - case consts::ElementType::TRI3: - case consts::ElementType::TRI6: - case consts::ElementType::QUD4: - case consts::ElementType::QUD8: - case consts::ElementType::QUD9: - return true; - default: - return false; + double sum = 0.0; + for (const double weight : rule.weights()) { + sum += weight; } + return sum; } } // namespace -TEST(QuadraturePhase01Baseline, StandardSelectionTableHasOrderedPointWeightData) -{ - const auto& cases = standard_solver_cases(); - ASSERT_EQ(cases.size(), 14u); - - for (const auto& c : cases) { - SCOPED_TRACE(c.name); - EXPECT_FALSE(c.samples.empty()); - EXPECT_GT(c.requested_exactness, 0); - EXPECT_GE(c.advertised_exactness, c.requested_exactness); - EXPECT_GT(c.reference_cell_measure, 0.0); - - double weight_sum = 0.0; - for (const auto& sample : c.samples) { - weight_sum += sample.weight; - } - EXPECT_NEAR(weight_sum, c.reference_cell_measure, kTol); - } -} - -TEST(QuadraturePhase01Baseline, CanonicalFixturesSatisfyTheRuleAndExactnessContracts) -{ - for (const auto& c : standard_solver_cases()) { - SCOPED_TRACE(c.name); - const auto rule = make_rule(c); - - EXPECT_EQ(rule.cell_family(), to_mesh_family(c.fe_type)); - EXPECT_EQ(rule.dimension(), c.dimension); - EXPECT_EQ(rule.polynomial_exactness(), c.advertised_exactness); - EXPECT_EQ(rule.num_points(), c.samples.size()); - EXPECT_DOUBLE_EQ( - rule.reference_cell_measure(), - c.reference_cell_measure); - expect_samples_in_order(rule, c.samples); - expect_total_degree_exact(rule, c.advertised_exactness); - } -} - -TEST(QuadraturePhase01Baseline, GenericFunctionSpaceRulesMatchOrderedLegacyData) -{ - for (const auto& c : standard_solver_cases()) { - if (!generic_context_supports(c.solver_type) || - c.solver_type == consts::ElementType::QUD9) { - // QUD9 is characterized separately because the legacy generic - // table contains an out-of-domain sqrt(6) coordinate defect. - continue; - } - - SCOPED_TRACE(c.name); - Vector weights(static_cast(c.samples.size())); - Array points(c.dimension, static_cast(c.samples.size())); - nn::get_gip(c.dimension, c.solver_type, - static_cast(c.samples.size()), weights, points); - expect_legacy_arrays_in_order(weights, points, c.dimension, c.samples); - } -} - -TEST(QuadraturePhase01Baseline, GenericQud9SqrtSixCoordinatesAreAnExplicitLegacyDefect) -{ - const auto legacy_defect = quad9_samples(std::sqrt(6.0)); - Vector weights(static_cast(legacy_defect.size())); - Array points(2, static_cast(legacy_defect.size())); - - nn::get_gip( - 2, - consts::ElementType::QUD9, - static_cast(legacy_defect.size()), - weights, - points); - - expect_legacy_arrays_in_order(weights, points, 2, legacy_defect); - EXPECT_GT(std::abs(points(0, 0)), 1.0); - - double x_squared_moment = 0.0; - for (int q = 0; q < weights.size(); ++q) { - x_squared_moment += weights[q] * points(0, q) * points(0, q); - } - EXPECT_NEAR(x_squared_moment, 40.0 / 3.0, kExactnessTol); - EXPECT_NEAR( - exact_reference_moment(svmp::CellFamily::Quad, 2, 0, 0), - 4.0 / 3.0, - kTol); - - EXPECT_THROW( - (void)RuleProbe( - svmp::CellFamily::Quad, - 5, - rule_points(points, 2), - rule_weights(weights)), - InvalidArgumentException); -} - -TEST(QuadraturePhase01Baseline, VolumeMeshRulesMatchOrderedLegacyData) -{ - for (const auto& c : standard_solver_cases()) { - if (!volume_context_supports(c.solver_type)) { - continue; - } - - SCOPED_TRACE(c.name); - mshType mesh; - mesh.eType = c.solver_type; - mesh.nG = static_cast(c.samples.size()); - mesh.w = Vector(mesh.nG); - mesh.xi = Array(c.dimension, mesh.nG); - nn::get_gip(mesh); - expect_legacy_arrays_in_order(mesh.w, mesh.xi, c.dimension, c.samples); - } -} - -TEST(QuadraturePhase01Baseline, Qud8VolumeSelectionHasNoLegacyPopulationEntry) -{ - const auto& qud8_case = standard_solver_cases()[6]; - ASSERT_EQ(qud8_case.solver_type, consts::ElementType::QUD8); - - mshType mesh; - mesh.eType = qud8_case.solver_type; - mesh.nG = static_cast(qud8_case.samples.size()); - mesh.w = Vector(mesh.nG); - mesh.xi = Array(qud8_case.dimension, mesh.nG); - - EXPECT_THROW(nn::get_gip(mesh), InvalidElementException); - - const auto canonical_rule = make_rule(qud8_case); - expect_samples_in_order(canonical_rule, qud8_case.samples); - expect_total_degree_exact( - canonical_rule, - qud8_case.advertised_exactness); -} - -TEST(QuadraturePhase01Baseline, GenericFunctionSpaceMissingEntriesAreExplicit) -{ - for (const auto& c : standard_solver_cases()) { - if (generic_context_supports(c.solver_type)) { - continue; - } - - SCOPED_TRACE(c.name); - Vector weights(static_cast(c.samples.size())); - Array points(c.dimension, static_cast(c.samples.size())); - EXPECT_THROW( - nn::get_gip(c.dimension, c.solver_type, - static_cast(c.samples.size()), weights, points), - InvalidElementException); - } -} - -TEST(QuadraturePhase01Baseline, FaceRulesMatchOrderedLegacyData) -{ - for (const auto& c : standard_solver_cases()) { - if (!face_context_supports(c.solver_type)) { - continue; - } - - SCOPED_TRACE(c.name); - faceType face; - face.eType = c.solver_type; - face.nG = static_cast(c.samples.size()); - face.w = Vector(face.nG); - face.xi = Array(c.dimension, face.nG); - nn::get_gip(nullptr, face); - expect_legacy_arrays_in_order(face.w, face.xi, c.dimension, c.samples); - } -} - -TEST(QuadraturePhase01Baseline, LegacyPositionModifiersAreOnlyDegreeOneExact) -{ - mshType tetra_mesh; - tetra_mesh.eType = consts::ElementType::TET4; - tetra_mesh.nG = 4; - tetra_mesh.qmTET4 = 0.25; - tetra_mesh.w = Vector(tetra_mesh.nG); - tetra_mesh.xi = Array(3, tetra_mesh.nG); - nn::get_gip(tetra_mesh); - - const RuleProbe tetra_rule( - svmp::CellFamily::Tetra, - 1, - rule_points(tetra_mesh.xi, 3), - rule_weights(tetra_mesh.w)); - expect_total_degree_exact(tetra_rule, 1); - EXPECT_NEAR( - integrate_monomial(tetra_rule, 2, 0, 0), - 1.0 / 96.0, - kExactnessTol); - EXPECT_GT( - std::abs( - integrate_monomial(tetra_rule, 2, 0, 0) - - exact_reference_moment(svmp::CellFamily::Tetra, 2, 0, 0)), - 1.0e-3); - - faceType triangle_face; - triangle_face.eType = consts::ElementType::TRI3; - triangle_face.nG = 3; - triangle_face.qmTRI3 = 1.0 / 3.0; - triangle_face.w = Vector(triangle_face.nG); - triangle_face.xi = Array(2, triangle_face.nG); - nn::get_gip(nullptr, triangle_face); - - const RuleProbe triangle_rule( - svmp::CellFamily::Triangle, - 1, - rule_points(triangle_face.xi, 2), - rule_weights(triangle_face.w)); - expect_total_degree_exact(triangle_rule, 1); - EXPECT_NEAR( - integrate_monomial(triangle_rule, 2, 0, 0), - 1.0 / 18.0, - kExactnessTol); - EXPECT_GT( - std::abs( - integrate_monomial(triangle_rule, 2, 0, 0) - - exact_reference_moment(svmp::CellFamily::Triangle, 2, 0, 0)), - 1.0e-3); -} - TEST(QuadPointContract, UsesFixedSizeFEVectorRepresentation) { static_assert(std::is_same_v>); @@ -955,7 +261,7 @@ TEST(QuadratureRuleValidation, EnforcesReferenceCellMeasureButAllowsNegativeWeig {-0.25, 0.75}); EXPECT_LT(rule.weight(0), 0.0); EXPECT_DOUBLE_EQ( - integrate_monomial(rule, 0, 0, 0), + weight_sum(rule), rule.reference_cell_measure()); } @@ -1140,7 +446,7 @@ TEST(QuadratureRuleValidation, AcceptsLargeSignedRule) std::move(weights)); EXPECT_LT(rule.weight(0), 0.0); EXPECT_NEAR( - integrate_monomial(rule, 0, 0, 0), + weight_sum(rule), rule.reference_cell_measure(), kTol); }