FE Quadrature Phase 01 abstract contract

Code/CMake/SimVascularExternals.cmake

@@ -6,11 +6,13 @@ find_package(Doxygen)
 if(DOXYGEN_FOUND)
   file(TO_NATIVE_PATH "${SV_BINARY_DIR}/Doxygen/" SV_DOCS_DIR_WORK)
   set(SV_DOCS_DIR ${SV_DOCS_DIR_WORK} CACHE PATH "Location to place docs")
-  configure_file(${SV_SOURCE_DIR}/../Documentation/Doxyfile
+  configure_file(${SV_SOURCE_DIR}/../Documentation/Doxyfile.cmake.in
     ${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_DOCS_DIR}"
+    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

@@ -267,11 +267,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,131 @@
+// 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 <cmath>
+#include <string>
+#include <utility>
+
+namespace svmp::FE::quadrature {
+namespace {
+
+std::size_t reference_dimension(svmp::CellFamily family)
+{
+    switch (family) {
+        case svmp::CellFamily::Point:
+            return 0u;
+        case svmp::CellFamily::Line:
+            return 1u;
+        case svmp::CellFamily::Triangle:
+        case svmp::CellFamily::Quad:
+            return 2u;
+        case svmp::CellFamily::Tetra:
+        case svmp::CellFamily::Hex:
+        case svmp::CellFamily::Wedge:
+            return 3u;
+        default:
+            svmp::raise<InvalidArgumentException>(
+                "QuadratureRule: unsupported reference-cell family");
+    }
+}
+
+double reference_measure(svmp::CellFamily family)
+{
+    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:
+            svmp::raise<InvalidArgumentException>(
+                "QuadratureRule: unsupported reference-cell family");
+    }
+}
+
+} // namespace
+
+QuadratureRule::~QuadratureRule() = default;
+
+std::size_t QuadratureRule::dimension() const
+{
+    return reference_dimension(cell_family_);
+}
+
+double QuadratureRule::reference_cell_measure() const
+{
+    return reference_measure(cell_family_);
+}
+
+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 std::size_t dimension = reference_dimension(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) {
+        const QuadPoint& point = points_[point_index];
+        for (std::size_t component = 0;
+             component < point.size();
+             ++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 >= 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));
+            }
+        }
+        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 svmp::FE::quadrature

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

@@ -0,0 +1,152 @@
+// 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 Abstract base for reference-space quadrature rules.
+ * @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 use these canonical reference cells:
+ *
+ * | Family   | Reference cell                            | Measure   |
+ * |----------|-------------------------------------------|-----------|
+ * | Point    | @f$(0,0,0)@f$                             | @f$1@f$   |
+ * | Line     | @f$[-1,1]@f$                              | @f$2@f$   |
+ * | Triangle | @f$\{(x,y):x,y\geq0,\ x+y\leq1\}@f$       | @f$1/2@f$ |
+ * | Quad     | @f$[-1,1]^2@f$                            | @f$4@f$   |
+ * | Tetra    | @f$\{(x,y,z):x,y,z\geq0,\ x+y+z\leq1\}@f$ | @f$1/6@f$ |
+ * | Hex      | @f$[-1,1]^3@f$                            | @f$8@f$   |
+ * | Wedge    | unit triangle @f$\times[-1,1]@f$          | @f$1@f$   |
+ *
+ * The family therefore determines both reference dimension and cell measure.
+ * Quadrature points are not required to lie inside the reference cell.
+ * 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 abstract base for quadrature rules on canonical reference cells.
+ *
+ * Construction requires:
+ *
+ * - a supported cell family and non-negative polynomial exactness;
+ * - a vector of points, with at least one element;
+ * - for all points, coordinates must be finite;
+ * - for all points, coordinates beyond the reference dimension of the cell
+ *   family must be equal to zero; and
+ * - a vector of finite weights, with as many elements as the points.
+ *
+ * The constructor checks each requirement and throws InvalidArgumentException
+ * when one is violated.
+ *
+ * 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 {
+public:
+    /** @brief Enable polymorphic destruction of concrete quadrature rules. */
+    virtual ~QuadratureRule() = 0;
+
+    /** @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.
+     */
+    std::size_t dimension() const;
+
+    /** @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;
+
+protected:
+    /**
+     * @brief Initialize a concrete 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);
+
+    QuadratureRule(const QuadratureRule&) = default;
+    QuadratureRule(QuadratureRule&&) noexcept = default;
+    QuadratureRule& operator=(const QuadratureRule&) = default;
+    QuadratureRule& operator=(QuadratureRule&&) noexcept = default;
+
+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

Documentation/Doxyfile.cmake.in

@@ -0,0 +1,5 @@
+# SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others.
+# SPDX-License-Identifier: BSD-3-Clause
+
+@INCLUDE = "@SV_SOURCE_DIR@/../Documentation/Doxyfile"
+OUTPUT_DIRECTORY = "@SV_DOCS_DIR@"