Issue 525 FE Basis refactor implementation

.github/workflows/documentation.yml

@@ -10,7 +10,7 @@ jobs:
       - uses: actions/checkout@v4
       - name: Build doxygen documentation
         run: |
-          sudo apt install -y doxygen
+          sudo apt install -y doxygen graphviz
           doxygen Documentation/Doxyfile
       - name: Save documentation
         uses: actions/upload-artifact@v4

Code/Source/solver/CMakeLists.txt

@@ -23,15 +23,18 @@
 # LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 # NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 
-set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED True)
 
 include_directories(${SV_SOURCE_DIR}/ThirdParty/eigen/include)
+include_directories(${SV_SOURCE_DIR}/ThirdParty/eigen/include/eigen3)
 include_directories(${SV_SOURCE_DIR}/ThirdParty/parmetis_internal/simvascular_parmetis_internal/ParMETISLib)
 include_directories(${SV_SOURCE_DIR}/ThirdParty/tetgen/simvascular_tetgen)
 include_directories(${SV_SOURCE_DIR}/ThirdParty/tinyxml/simvascular_tinyxml)
 include_directories(${MPI_C_INCLUDE_PATH})
+include_directories(${CMAKE_CURRENT_SOURCE_DIR})
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/Core)
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/FE)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/FE/Common)
 
 # Find Trilinos package if requested
@@ -86,7 +89,7 @@ endif()
 # add trilinos flags and defines
 if(USE_TRILINOS)
   ADD_DEFINITIONS(-DWITH_TRILINOS)
-  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++20")
 endif()
 
 # Build with the PETSc linear algebra package.
@@ -245,9 +248,21 @@ file(GLOB SOLVER_FE_COMMON_SRCS CONFIGURE_DEPENDS
   FE/Common/*.h
 )
 
+file(GLOB SOLVER_FE_BASIS_SRCS CONFIGURE_DEPENDS
+  FE/Basis/*.cpp
+  FE/Basis/*.h
+)
+
+file(GLOB SOLVER_FE_MATH_SRCS CONFIGURE_DEPENDS
+  FE/Math/*.cpp
+  FE/Math/*.h
+)
+
 list(APPEND CSRCS
   ${SOLVER_CORE_SRCS}
   ${SOLVER_FE_COMMON_SRCS}
+  ${SOLVER_FE_BASIS_SRCS}
+  ${SOLVER_FE_MATH_SRCS}
 )
 
   # Set PETSc interace code.
@@ -324,16 +339,24 @@ if(ENABLE_UNIT_TEST)
 
   # link pthread on ubuntu20
   find_package(Threads REQUIRED)
-
+  include(FetchContent)
   # install Google Test
   #if(NOT TARGET gtest_main AND NOT TARGET gtest)
-  include(FetchContent)
   FetchContent_Declare(
-    googletest
-    URL https://github.com/google/googletest/archive/refs/heads/main.zip
-    DOWNLOAD_EXTRACT_TIMESTAMP TRUE
+          googletest
+          GIT_REPOSITORY https://github.com/google/googletest.git
+          GIT_TAG v1.17.0
+          DOWNLOAD_EXTRACT_TIMESTAMP TRUE
   )
   FetchContent_MakeAvailable(googletest)
+
+  if(CMAKE_CXX_COMPILER_ID MATCHES "Clang" AND CMAKE_CXX_STANDARD GREATER_EQUAL 20)
+    foreach(GTEST_TARGET gtest gtest_main gmock gmock_main)
+      if(TARGET ${GTEST_TARGET})
+        target_compile_options(${GTEST_TARGET} PRIVATE -std=gnu++17)
+      endif()
+    endforeach()
+  endif()
   #endif()
 
   enable_testing()

Code/Source/solver/FE/Basis/BasisExceptions.h

@@ -0,0 +1,90 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others.
+// SPDX-License-Identifier: BSD-3-Clause
+
+#ifndef SVMP_FE_BASIS_BASISEXCEPTIONS_H
+#define SVMP_FE_BASIS_BASISEXCEPTIONS_H
+
+#include "FEException.h"
+
+namespace svmp {
+namespace FE {
+namespace basis {
+
+/**
+ * @brief Base exception type for errors originating in the Basis module
+ */
+class BasisException : public FEException {
+public:
+    BasisException(const std::string& message,
+                   const char* file = "",
+                   int line = 0,
+                   const char* function = "",
+                   StatusCode status = StatusCode::Unknown)
+        : FEException(message, status, file, line, function) {}
+};
+
+/**
+ * @brief Invalid Basis request or configuration
+ */
+class BasisConfigurationException : public BasisException {
+public:
+    BasisConfigurationException(const std::string& message,
+                                const char* file = "",
+                                int line = 0,
+                                const char* function = "")
+        : BasisException(message, file, line, function, StatusCode::InvalidArgument) {}
+};
+
+/**
+ * @brief Requested element topology is incompatible with the basis family
+ */
+class BasisElementCompatibilityException : public BasisException {
+public:
+    BasisElementCompatibilityException(const std::string& message,
+                                       const char* file = "",
+                                       int line = 0,
+                                       const char* function = "")
+        : BasisException(message, file, line, function, StatusCode::InvalidArgument) {}
+};
+
+/**
+ * @brief Basis evaluation request cannot be satisfied
+ */
+class BasisEvaluationException : public BasisException {
+public:
+    BasisEvaluationException(const std::string& message,
+                             const char* file = "",
+                             int line = 0,
+                             const char* function = "")
+        : BasisException(message, file, line, function, StatusCode::InvalidArgument) {}
+};
+
+/**
+ * @brief Public-to-canonical node ordering or coordinate lookup failure
+ */
+class BasisNodeOrderingException : public BasisException {
+public:
+    BasisNodeOrderingException(const std::string& message,
+                               const char* file = "",
+                               int line = 0,
+                               const char* function = "")
+        : BasisException(message, file, line, function, StatusCode::InvalidArgument) {}
+};
+
+/**
+ * @brief Internal basis construction or transform setup failure
+ */
+class BasisConstructionException : public BasisException {
+public:
+    BasisConstructionException(const std::string& message,
+                               const char* file = "",
+                               int line = 0,
+                               const char* function = "")
+        : BasisException(message, file, line, function, StatusCode::InternalError) {}
+};
+
+} // namespace basis
+} // namespace FE
+} // namespace svmp
+
+#endif // SVMP_FE_BASIS_BASISEXCEPTIONS_H

Code/Source/solver/FE/Basis/BasisFactory.cpp

@@ -0,0 +1,98 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others.
+// SPDX-License-Identifier: BSD-3-Clause
+
+#include "BasisFactory.h"
+
+#include "BasisTraits.h"
+#include "LagrangeBasis.h"
+#include "SerendipityBasis.h"
+
+namespace svmp {
+namespace FE {
+namespace basis {
+
+namespace {
+
+int require_basis_order(const BasisRequest& req,
+                        const char* missing_message,
+                        const char* negative_message) {
+    FE::throw_if<BasisConfigurationException>(!req.order.has_value(), SVMP_HERE,
+                                              missing_message);
+    FE::throw_if<BasisConfigurationException>(*req.order < 0, SVMP_HERE,
+                                              negative_message);
+    return *req.order;
+}
+
+void require_scalar_c0_request(const BasisRequest& req) {
+    FE::throw_if<BasisConfigurationException>(
+        req.field_type != FieldType::Scalar, SVMP_HERE,
+        "BasisFactory: Lagrange/Serendipity bases support scalar fields only");
+    FE::throw_if<BasisConfigurationException>(
+        req.continuity != Continuity::C0, SVMP_HERE,
+        "BasisFactory: Lagrange/Serendipity bases support C0 continuity only");
+}
+
+std::shared_ptr<BasisFunction> create_lagrange(const BasisRequest& req) {
+    require_scalar_c0_request(req);
+    const int order = require_basis_order(
+        req,
+        "BasisFactory: Lagrange creation requires an explicit order",
+        "BasisFactory: Lagrange requires non-negative order");
+    return std::make_shared<LagrangeBasis>(req.element_type, order);
+}
+
+std::shared_ptr<BasisFunction> create_serendipity(const BasisRequest& req) {
+    require_scalar_c0_request(req);
+    const int order = require_basis_order(
+        req,
+        "BasisFactory: Serendipity creation requires an explicit order",
+        "BasisFactory: Serendipity requires non-negative order");
+    return std::make_shared<SerendipityBasis>(req.element_type, order);
+}
+
+} // namespace
+
+namespace basis_factory {
+
+std::shared_ptr<BasisFunction> create(const BasisRequest& req) {
+    switch (req.basis_type) {
+        case BasisType::Lagrange:
+            return create_lagrange(req);
+        case BasisType::Serendipity:
+            return create_serendipity(req);
+        default:
+            FE::raise<BasisConfigurationException>(SVMP_HERE,
+                "BasisFactory: requested basis family is outside the scalar Lagrange/Serendipity scope");
+    }
+}
+
+BasisRequest default_basis_request(ElementType element_type) {
+    switch (element_type) {
+        // Reduced serendipity node layouts have no complete Lagrange basis at
+        // their node count; they always use the quadratic serendipity space.
+        case ElementType::Quad8:
+        case ElementType::Hex20:
+        case ElementType::Wedge15:
+            return BasisRequest{element_type, BasisType::Serendipity, 2};
+        case ElementType::Point1:
+            return BasisRequest{element_type, BasisType::Lagrange, 0};
+        default: {
+            const int order = complete_lagrange_alias_order(element_type);
+            if (order >= 0) {
+                return BasisRequest{element_type, BasisType::Lagrange, order};
+            }
+            FE::raise<BasisElementCompatibilityException>(SVMP_HERE,
+                "BasisFactory: no default basis is defined for the requested element type");
+        }
+    }
+}
+
+std::shared_ptr<BasisFunction> create_default_for(ElementType element_type) {
+    return create(default_basis_request(element_type));
+}
+
+} // namespace basis_factory
+
+} // namespace basis
+} // namespace FE
+} // namespace svmp

Code/Source/solver/FE/Basis/BasisFactory.h

@@ -0,0 +1,71 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others.
+// SPDX-License-Identifier: BSD-3-Clause
+
+#ifndef SVMP_FE_BASIS_BASISFACTORY_H
+#define SVMP_FE_BASIS_BASISFACTORY_H
+
+/**
+ * @file BasisFactory.h
+ * @brief Runtime creation of basis families
+ */
+
+#include "BasisFunction.h"
+#include <memory>
+#include <optional>
+#include <string>
+#include <vector>
+
+namespace svmp {
+namespace FE {
+namespace basis {
+
+struct BasisRequest {
+    ElementType element_type;
+    BasisType basis_type;
+    std::optional<int> order{};
+    Continuity continuity{Continuity::C0};
+    FieldType field_type{FieldType::Scalar};
+    std::vector<Real> knot_vector{};
+    std::vector<Real> weights{};
+    std::vector<int> axis_orders{};
+    std::vector<std::vector<Real>> axis_knot_vectors{};
+    std::vector<std::vector<Real>> axis_weights{};
+    std::vector<int> tensor_extents{};
+    std::string custom_id{};
+};
+
+namespace basis_factory {
+
+[[nodiscard]] std::shared_ptr<BasisFunction> create(const BasisRequest& req);
+
+/// \brief Return the default basis request (family and order) for an element type.
+///
+/// \details This is the single source of truth for which basis family and
+/// polynomial order a given element type uses by default: serendipity node
+/// layouts (Quad8, Hex20, Wedge15) select the quadratic serendipity family,
+/// and every complete Lagrange element selects the Lagrange family at the
+/// order given by its node layout. Solver-facing adapters should translate
+/// their element names to ElementType and delegate the basis choice here
+/// rather than tabulating family/order themselves.
+///
+/// \param element_type Element type to select a default basis for.
+/// \return Basis request suitable for create().
+/// \throws BasisElementCompatibilityException If no default basis is defined
+///         for the element type.
+[[nodiscard]] BasisRequest default_basis_request(ElementType element_type);
+
+/// \brief Create the default basis for an element type.
+///
+/// \details Equivalent to create(default_basis_request(element_type)).
+///
+/// \param element_type Element type to create a default basis for.
+/// \return Shared basis instance.
+[[nodiscard]] std::shared_ptr<BasisFunction> create_default_for(ElementType element_type);
+
+} // namespace basis_factory
+
+} // namespace basis
+} // namespace FE
+} // namespace svmp
+
+#endif // SVMP_FE_BASIS_BASISFACTORY_H