diff --git a/Code/CMake/SimVascularExternals.cmake b/Code/CMake/SimVascularExternals.cmake index 50772c5b9..880882e3e 100644 --- a/Code/CMake/SimVascularExternals.cmake +++ b/Code/CMake/SimVascularExternals.cmake @@ -9,8 +9,10 @@ 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 "${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 ) endif(DOXYGEN_FOUND) diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index c5ab81146..30baa33b3 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -266,11 +266,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). */ /** diff --git a/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp new file mode 100644 index 000000000..9a07c41fb --- /dev/null +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.cpp @@ -0,0 +1,444 @@ +// 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 +#include +#include +#include +#include +#include +#include + +namespace svmp::FE::quadrature { +namespace { + +struct ReferenceCellTraits { + int dimension; + double reference_cell_measure; +}; + +constexpr double coordinate_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. +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(); + + std::string_view reason{}; + std::size_t sample{no_sample}; + + constexpr bool valid() const noexcept + { + return reason.empty(); + } +}; + +constexpr std::optional reference_cell_traits( + svmp::CellFamily family) noexcept +{ + switch (family) { + case svmp::CellFamily::Point: + return ReferenceCellTraits{0, 1.0}; + case svmp::CellFamily::Line: + return ReferenceCellTraits{1, 2.0}; + case svmp::CellFamily::Triangle: + return ReferenceCellTraits{2, 0.5}; + case svmp::CellFamily::Quad: + return ReferenceCellTraits{2, 4.0}; + case svmp::CellFamily::Tetra: + return ReferenceCellTraits{3, 1.0 / 6.0}; + case svmp::CellFamily::Hex: + return ReferenceCellTraits{3, 8.0}; + case svmp::CellFamily::Wedge: + return ReferenceCellTraits{3, 1.0}; + default: + return std::nullopt; + } +} + +ValidationResult validate_point( + const QuadPoint& point, + svmp::CellFamily family, + const ReferenceCellTraits& traits, + std::size_t sample) noexcept +{ + for (std::size_t component = 0; component < 3u; ++component) { + if (!std::isfinite(point[component])) { + return {"quadrature point contains a non-finite coordinate", sample}; + } + if (component >= + static_cast(traits.dimension) && + std::abs(point[component]) > coordinate_validation_tolerance) { + return {"quadrature point has a nonzero inactive coordinate", sample}; + } + } + + 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]; + const double y = point[1]; + const double z = point[2]; + + bool is_contained = false; + switch (family) { + case svmp::CellFamily::Point: + is_contained = true; + break; + case svmp::CellFamily::Line: + is_contained = in_interval(x, -1.0, 1.0); + break; + case svmp::CellFamily::Triangle: + 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 >= -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) && + in_interval(y, -1.0, 1.0) && + in_interval(z, -1.0, 1.0); + break; + case svmp::CellFamily::Wedge: + 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"}; + } + + if (!is_contained) { + return { + "quadrature point lies outside the canonical reference cell", + sample}; + } + return {}; +} + +ValidationResult validate_weights( + const std::vector& weights, + double reference_cell_measure) noexcept +{ + ExactBinary64Sum exact_sum; + + for (std::size_t sample = 0; sample < weights.size(); ++sample) { + if (!exact_sum.add(weights[sample])) { + return {"quadrature weight must be finite", sample}; + } + } + + 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 {}; +} + +ValidationResult validate_rule_data( + svmp::CellFamily family, + const ReferenceCellTraits& traits, + int polynomial_exactness, + const std::vector& points, + const std::vector& weights) noexcept +{ + if (polynomial_exactness < 0) { + return {"polynomial exactness must be non-negative"}; + } + if (points.empty()) { + return {"a rule must contain at least one sample"}; + } + if (points.size() != weights.size()) { + return {"points/weights size mismatch"}; + } + + for (std::size_t sample = 0; sample < points.size(); ++sample) { + const auto result = + validate_point(points[sample], family, traits, sample); + if (!result.valid()) { + return result; + } + } + + return validate_weights(weights, traits.reference_cell_measure); +} + +std::string validation_failure_message(const ValidationResult& result) +{ + if (result.valid()) { + return {}; + } + + 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 message; +} + +} // 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_cell_measure_(state.reference_cell_measure), + points_(std::move(state.points)), + weights_(std::move(state.weights)) +{ +} + +QuadratureRule::ValidatedState QuadratureRule::validate( + svmp::CellFamily family, + RuleData data) +{ + const auto traits = reference_cell_traits(family); + if (!traits) { + svmp::raise( + validation_failure_message( + {"unsupported reference-cell family"})); + } + + const auto validation = validate_rule_data( + family, + *traits, + data.polynomial_exactness, + data.points, + data.weights); + if (!validation.valid()) { + svmp::raise( + validation_failure_message(validation)); + } + + return { + family, + traits->dimension, + data.polynomial_exactness, + traits->reference_cell_measure, + std::move(data.points), + std::move(data.weights), + }; +} + +} // 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..3cdb42c2a --- /dev/null +++ b/Code/Source/solver/FE/Quadrature/QuadratureRule.h @@ -0,0 +1,304 @@ +// 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 unweighted 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] + * The reference-cell measure satisfies + * @f[ + * |\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 + * 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, + * exact 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. + * 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 + * + * 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 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 + * 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 + * 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 + * 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 | Reference-cell 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. Concurrent const access + * is safe while an owning handle keeps the rule alive. + */ + +#include "FE/Common/Types.h" +#include "FE/Math/Vector.h" + +#include +#include + +namespace svmp::FE::quadrature { + +/** @addtogroup 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. + */ +using QuadPoint = math::Vector; + +/** + * @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. 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: + /** @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. + * + * 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. + */ + int polynomial_exactness() const noexcept { return polynomial_exactness_; } + + /** + * @brief Return the dimension of the canonical reference cell. + * + * This is also the number of active components in each QuadPoint. + * + * @return Reference dimension, 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 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[ + * |\hat K| = \int_{\hat K} 1\,d\hat x = \sum_q w_q. + * @f] + * + * @return Geometric measure of the canonical reference cell. + */ + double reference_cell_measure() const noexcept { return reference_cell_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-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. + * @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, or the weights do not reproduce the + * reference-cell measure within the scaled measure tolerance. + */ + 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_cell_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_cell_measure_; ///< Canonical reference-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 diff --git a/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp new file mode 100644 index 000000000..f0f3465c6 --- /dev/null +++ b/tests/unitTests/FE/Quadrature/test_QuadratureRules.cpp @@ -0,0 +1,540 @@ +/** + * @file test_QuadratureRules.cpp + * @brief Unit tests for the core quadrature rule infrastructure. + */ + +#include + +#include "FE/Common/FEException.h" +#include "FE/Quadrature/QuadratureRule.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace svmp::FE; +using namespace svmp::FE::quadrature; + +namespace { + +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; + +using ExpectedPoint = std::array; + +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"; + } +} + +double weight_sum(const QuadratureRule& rule) +{ + double sum = 0.0; + for (const double weight : rule.weights()) { + sum += weight; + } + return sum; +} + +} // namespace + +TEST(QuadPointContract, UsesFixedSizeFEVectorRepresentation) +{ + static_assert(std::is_same_v>); + static_assert(QuadPoint::RowsAtCompileTime == 3); + static_assert(QuadPoint::ColsAtCompileTime == 1); + + const QuadPoint origin = QuadPoint::Zero(); + for (std::size_t component = 0; component < 3u; ++component) { + EXPECT_DOUBLE_EQ(origin[component], 0.0); + } + + 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, 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 mutable_point = QuadPoint::Zero(); + mutable_point[2] = 0.75; + EXPECT_DOUBLE_EQ(mutable_point[2], 0.75); + + 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); + } + } +} + +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_cell_measure(), + c.expected_measure); + } +} + +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"); +} + +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, EnforcesReferenceCellMeasureButAllowsNegativeWeights) +{ + 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_LT(rule.weight(0), 0.0); + EXPECT_DOUBLE_EQ( + weight_sum(rule), + rule.reference_cell_measure()); +} + +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, RejectsIncorrectMeasureHiddenByCancellation) +{ + 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_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 reference-cell measure"); +} + +TEST(QuadratureRuleValidation, ExactMeasureValidationIsOrderIndependent) +{ + 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 reference-cell measure"); + ++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()), + { + 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 reference-cell measure"); +} + +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 reference-cell measure"); +} + +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, sample_weight); + + const RuleProbe rule( + svmp::CellFamily::Line, + 0, + std::move(points), + std::move(weights)); + EXPECT_EQ(rule.num_points(), sample_count); + EXPECT_DOUBLE_EQ(rule.weight(0), sample_weight); +} + +TEST(QuadratureRuleValidation, AcceptsLargeSignedRule) +{ + constexpr std::size_t sample_count = 8192u; + constexpr double cancelling_weight = 1.0 / 1048576.0; + std::vector points(sample_count, QuadPoint::Zero()); + 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_LT(rule.weight(0), 0.0); + EXPECT_NEAR( + weight_sum(rule), + rule.reference_cell_measure(), + kTol); +} + +TEST(QuadratureRuleValidation, AppliesConstructionCoordinateTolerance) +{ + constexpr double accepted_offset = 1.0e-14; + const RuleProbe rule( + 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 = 1.0e-10; + 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, AppliesConstructionWeightTolerance) +{ + constexpr double accepted_offset = 1.0e-14; + const RuleProbe rule( + svmp::CellFamily::Triangle, + 0, + {{0.25, 0.25, 0.0}}, + {0.5 + accepted_offset}); + EXPECT_DOUBLE_EQ(rule.weight(0), 0.5 + accepted_offset); + + constexpr double rejected_offset = 1.0e-10; + 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 reference-cell measure"); +} + +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_cell_measure(), 2.0); + 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); +}