FE Quadrature Phase 01 abstract contract

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)

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.

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
  */

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).
  */
 
 /**

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

@@ -0,0 +1,156 @@
+// 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 <cassert>
+#include <cmath>
+#include <string>
+#include <utility>
+
+namespace svmp::FE::quadrature {
+namespace {
+
+constexpr int reference_dimension(svmp::CellFamily family) noexcept
+{
+    switch (family) {
+        case svmp::CellFamily::Point:
+            return 0;
+        case svmp::CellFamily::Line:
+            return 1;
+        case svmp::CellFamily::Triangle:
+        case svmp::CellFamily::Quad:
+            return 2;
+        case svmp::CellFamily::Tetra:
+        case svmp::CellFamily::Hex:
+        case svmp::CellFamily::Wedge:
+            return 3;
+        default:
+            return -1;
+    }
+}
+
+constexpr double reference_measure(svmp::CellFamily family) noexcept
+{
+    switch (family) {
+        case svmp::CellFamily::Point:
+            return 1.0;
+        case svmp::CellFamily::Line:
+            return 2.0;
+        case svmp::CellFamily::Triangle:
+            return 0.5;
+        case svmp::CellFamily::Quad:
+            return 4.0;
+        case svmp::CellFamily::Tetra:
+            return 1.0 / 6.0;
+        case svmp::CellFamily::Hex:
+            return 8.0;
+        case svmp::CellFamily::Wedge:
+            return 1.0;
+        default:
+            return -1.0;
+    }
+}
+
+void validate_point(
+    const QuadPoint& point,
+    int dimension,
+    std::size_t point_index)
+{
+    for (std::size_t component = 0; component < 3u; ++component) {
+        if (!std::isfinite(point[component])) {
+            svmp::raise<InvalidArgumentException>(
+                std::string{
+                    "QuadratureRule: quadrature point contains a non-finite "
+                    "coordinate at point index "} +
+                std::to_string(point_index));
+        }
+        if (component >=
+            static_cast<std::size_t>(dimension) &&
+            point[component] != 0.0) {
+            svmp::raise<InvalidArgumentException>(
+                std::string{
+                    "QuadratureRule: quadrature point has a nonzero inactive "
+                    "coordinate at point index "} +
+                std::to_string(point_index));
+        }
+    }
+}
+
+void validate_weights(const std::vector<double>& weights)
+{
+    for (std::size_t point_index = 0;
+         point_index < weights.size();
+         ++point_index) {
+        if (!std::isfinite(weights[point_index])) {
+            svmp::raise<InvalidArgumentException>(
+                std::string{
+                    "QuadratureRule: quadrature weight must be finite at point "
+                    "index "} +
+                std::to_string(point_index));
+        }
+    }
+}
+
+} // namespace
+
+int QuadratureRule::dimension() const noexcept
+{
+    const int dimension = reference_dimension(cell_family_);
+    assert(dimension >= 0);
+    return dimension;
+}
+
+double QuadratureRule::reference_cell_measure() const noexcept
+{
+    const double measure = reference_measure(cell_family_);
+    assert(measure > 0.0);
+    return measure;
+}
+
+QuadratureRule::QuadratureRule(
+    svmp::CellFamily family,
+    int polynomial_exactness,
+    std::vector<QuadPoint> points,
+    std::vector<double> weights)
+    : cell_family_(family),
+      polynomial_exactness_(polynomial_exactness),
+      points_(std::move(points)),
+      weights_(std::move(weights))
+{
+    const int dimension = reference_dimension(cell_family_);
+    svmp::check<InvalidArgumentException>(
+        dimension >= 0,
+        "QuadratureRule: unsupported reference-cell family");
+
+    svmp::check<InvalidArgumentException>(
+        polynomial_exactness_ >= 0,
+        "QuadratureRule: polynomial exactness must be non-negative");
+    svmp::check<InvalidArgumentException>(
+        !points_.empty(),
+        "QuadratureRule: a rule must contain at least one point");
+    svmp::check<InvalidArgumentException>(
+        points_.size() == weights_.size(),
+        "QuadratureRule: points/weights size mismatch");
+
+    for (std::size_t point_index = 0;
+         point_index < points_.size();
+         ++point_index) {
+        validate_point(
+            points_[point_index],
+            dimension,
+            point_index);
+    }
+
+    validate_weights(weights_);
+}
+
+} // namespace svmp::FE::quadrature

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

@@ -0,0 +1,126 @@
+// 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 Reference-space quadrature rule value type.
+ * @ingroup FE_Quadrature
+ */
+
+/**
+ * @defgroup FE_Quadrature Quadrature
+ * @ingroup FE
+ * @brief Integration rules on canonical reference cells.
+ *
+ * @details
+ * A QuadratureRule owns ordered reference coordinates and weights for
+ * @f[
+ *   \int_{\hat K} f(\hat x)\,d\hat x
+ *   \approx \sum_q w_q f(\hat x_q).
+ * @f]
+ * Supported families are Point, Line, Triangle, Quad, Tetra, Hex, and Wedge.
+ * The family determines the reference dimension and cell measure. Generating
+ * code is responsible for verifying weight normalization and declared
+ * polynomial exactness through analytic moment tests.
+ */
+
+#include "FE/Common/Types.h"
+#include "FE/Math/Vector.h"
+
+#include <cstddef>
+#include <vector>
+
+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; remaining
+ * components must be zero.
+ */
+using QuadPoint = math::Vector<double, 3>;
+
+/**
+ * @brief Owning quadrature rule on a canonical reference cell.
+ *
+ * Construction requires:
+ *
+ * - a supported cell family and non-negative polynomial exactness;
+ * - at least one point and the same number of points and weights;
+ * - finite coordinates and weights; and
+ * - inactive coordinates equal to zero.
+ *
+ * Points may be duplicate or outside the reference cell, and weights may be
+ * zero or negative. Construction does not verify weight normalization or
+ * polynomial exactness.
+ */
+class QuadratureRule final {
+public:
+    /**
+     * @brief Construct a rule from complete point and weight data.
+     * @param family Reference-cell family; also determines dimension and measure.
+     * @param polynomial_exactness Declared total-degree polynomial exactness.
+     * @param points Ordered reference coordinates.
+     * @param weights Weights paired with @p points.
+     * @throws InvalidArgumentException If a construction requirement is violated.
+     */
+    explicit QuadratureRule(
+        svmp::CellFamily family,
+        int polynomial_exactness,
+        std::vector<QuadPoint> points,
+        std::vector<double> weights);
+
+    /** @brief Return the number of point/weight pairs. */
+    std::size_t num_points() const noexcept { return points_.size(); }
+
+    /** @brief Return the declared total-degree polynomial exactness. */
+    int polynomial_exactness() const noexcept { return polynomial_exactness_; }
+
+    /**
+     * @brief Return the reference dimension and active QuadPoint component count.
+     */
+    int dimension() const noexcept;
+
+    /** @brief Return the canonical reference-cell family. */
+    svmp::CellFamily cell_family() const noexcept { return cell_family_; }
+
+    /**
+     * @brief Return point @p i without bounds checking.
+     * @pre @p i is less than num_points().
+     */
+    const QuadPoint& point(std::size_t i) const noexcept { return points_[i]; }
+
+    /**
+     * @brief Return the weight paired with point @p i without bounds checking.
+     * @pre @p i is less than num_points().
+     */
+    double weight(std::size_t i) const noexcept { return weights_[i]; }
+
+    /** @brief Return all points in integration order. */
+    const std::vector<QuadPoint>& points() const noexcept { return points_; }
+
+    /** @brief Return all weights in point order. */
+    const std::vector<double>& weights() const noexcept { return weights_; }
+
+    /** @brief Return the reference-cell measure derived from cell_family(). */
+    double reference_cell_measure() const noexcept;
+
+private:
+    svmp::CellFamily cell_family_;          ///< Canonical reference topology.
+    int polynomial_exactness_;              ///< Exactness declared by the generator.
+    std::vector<QuadPoint> points_;          ///< Ordered reference coordinates.
+    std::vector<double> weights_;            ///< Weights paired with points_.
+};
+
+/** @} */
+
+} // namespace svmp::FE::quadrature
+
+#endif // SVMP_FE_QUADRATURE_RULE_H