Merge pull request #1614 from PrincetonUniversity/issue-1580

Issue 1580 - Implementing ChunkedFaceIterator & assembly::face_types

Author: Rohit-Kakodkar

core/specfem/assembly.hpp

@@ -5,6 +5,8 @@
 #include "assembly/boundary_values.hpp"
 #include "assembly/conforming_interfaces.hpp"
 #include "assembly/element_types.hpp"
+#include "assembly/edge_types.hpp"
+#include "assembly/face_types.hpp"
 #include "assembly/fields.hpp"
 #include "assembly/jacobian_matrix.hpp"
 #include "assembly/kernels.hpp"

core/specfem/assembly/edge_types/dim2/edge_types.hpp

@@ -9,40 +9,14 @@
 namespace specfem::assembly {
 
 /**
- * @brief 2D spectral element edge classification and coupling management
- *
- * This template specialization provides storage and management for edge
- * information in 2D spectral element meshes. It handles edge connectivity,
- * interface types, and boundary conditions for coupling between different
- * media types in 2D wave propagation problems.
+ * @brief Individual 2D edge representation with quadrature point access
  *
- * @code
- * // Construct 2D edge types from mesh data
- * specfem::assembly::edge_types<specfem::element::dimension_tag::dim2> edges(
- *     ngllx, ngllz, mesh, element_types);
+ * This structure represents a single edge in the 2D spectral element mesh,
+ * providing access to quadrature points along the edge for coupling
+ * computations, boundary condition enforcement, and flux calculations.
  *
- * // Get elastic-acoustic coupling edges on device
- * auto [self_edges, coupled_edges] = edges.get_edges_on_device(
- *     specfem::connections::type::weakly_conforming,
- *     specfem::interface::interface_tag::elastic_acoustic,
- *     specfem::element::boundary_tag::none);
- * @endcode
+ * @tparam ExecutionSpace Kokkos execution space (host or device)
  */
-template <> struct edge_types<specfem::element::dimension_tag::dim2> {
-
-public:
-  constexpr static auto dimension_tag =
-      specfem::element::dimension_tag::dim2; ///< Dimension tag
-
-  /**
-   * @brief Individual 2D edge representation with quadrature point access
-   *
-   * This structure represents a single edge in the 2D spectral element mesh,
-   * providing access to quadrature points along the edge for coupling
-   * computations, boundary condition enforcement, and flux calculations.
-   *
-   * @tparam ExecutionSpace Kokkos execution space (host or device)
-   */
   template <typename ExecutionSpace> struct Edge {
     int n_points; ///< Number of quadrature points on this edge
     using IndexView = Kokkos::View<int *, Kokkos::LayoutStride,
@@ -88,114 +62,142 @@ template <> struct edge_types<specfem::element::dimension_tag::dim2> {
     }
   };
 
+/**
+ * @brief Collection of 2D edges with parallel access capabilities
+ *
+ * This structure manages collections of edges for efficient parallel
+ * processing of edge-based operations such as coupling computations,
+ * boundary condition enforcement, and flux calculations in 2D.
+ *
+ * @tparam ExecutionSpace Kokkos execution space (host or device)
+ * @tparam Layout Memory layout for Kokkos views
+ */
+template <typename ExecutionSpace,
+          typename Layout = typename ExecutionSpace::array_layout>
+struct EdgeView {
+  int n_edges;  ///< Number of edges in this view
+  int n_points; ///< Number of quadrature points per edge
+  using IndexView =
+      Kokkos::View<int *, Layout, ExecutionSpace>; ///< View type for integer
+                                                    ///< indices
+  using QPView =
+      Kokkos::View<int **, Layout, ExecutionSpace>; ///< View type for
+                                                    ///< quadrature point
+                                                    ///< arrays
+  using EdgeTypeView = ///< View type for 2D edge classifications
+      Kokkos::View<specfem::mesh_entity::dim2::type *, ExecutionSpace>;
+
+  using HostMirror = std::conditional_t<
+      std::is_same<typename ExecutionSpace::memory_space,
+                    Kokkos::HostSpace>::value,
+      EdgeView, EdgeView<Kokkos::DefaultHostExecutionSpace, Layout> >;
+
   /**
-   * @brief Collection of 2D edges with parallel access capabilities
-   *
-   * This structure manages collections of edges for efficient parallel
-   * processing of edge-based operations such as coupling computations,
-   * boundary condition enforcement, and flux calculations in 2D.
+   * @brief Default constructor creating empty edge view.
+   */
+  EdgeView() : n_edges(0), n_points(0) {}
+
+  /**
+   * @brief Construct edge view with allocated storage.
    *
-   * @tparam ExecutionSpace Kokkos execution space (host or device)
-   * @tparam Layout Memory layout for Kokkos views
+   * @param label Base label for Kokkos view names
+   * @param n_edges Number of edges to allocate
+   * @param n_points Number of quadrature points per edge
    */
-  template <typename ExecutionSpace,
-            typename Layout = typename ExecutionSpace::array_layout>
-  struct EdgeView {
-    int n_edges;  ///< Number of edges in this view
-    int n_points; ///< Number of quadrature points per edge
-    using IndexView =
-        Kokkos::View<int *, Layout, ExecutionSpace>; ///< View type for integer
-                                                     ///< indices
-    using QPView =
-        Kokkos::View<int **, Layout, ExecutionSpace>; ///< View type for
-                                                      ///< quadrature point
-                                                      ///< arrays
-    using EdgeTypeView = ///< View type for 2D edge classifications
-        Kokkos::View<specfem::mesh_entity::dim2::type *, ExecutionSpace>;
+  EdgeView(const std::string &label, const int n_edges, const int n_points)
+      : n_edges(n_edges), n_points(n_points),
+        element_index(label + "_element_index", n_edges),
+        edge_index(label + "_edge_index", n_edges),
+        edge_types(label + "_edge_types", n_edges),
+        iz(label + "_iz", n_edges, n_points),
+        ix(label + "_ix", n_edges, n_points) {}
 
-    using HostMirror = std::conditional_t<
-        std::is_same<typename ExecutionSpace::memory_space,
-                     Kokkos::HostSpace>::value,
-        EdgeView, EdgeView<Kokkos::DefaultHostExecutionSpace, Layout> >;
+  IndexView element_index; ///< Element indices for each edge
+  IndexView edge_index;    ///< Local edge indices within elements
+  EdgeTypeView edge_types; ///< 2D edge type classifications
+  QPView iz;               ///< Z-direction quadrature indices for all edges
+  QPView ix;               ///< X-direction quadrature indices for all edges
 
-    /**
-     * @brief Default constructor creating empty edge view.
-     */
-    EdgeView() : n_edges(0), n_points(0) {}
+  /**
+   * @brief Device-side constructor from existing views.
+   *
+   * @param n_edges Number of edges
+   * @param n_points Number of quadrature points per edge
+   * @param element_index Element indices view
+   * @param edge_index Edge indices view
+   * @param edge_types Edge types view
+   * @param iz Z-direction quadrature indices
+   * @param ix X-direction quadrature indices
+   */
+  KOKKOS_INLINE_FUNCTION
+  EdgeView(const int n_edges, const int n_points,
+            const IndexView &element_index, const IndexView &edge_index,
+            const EdgeTypeView &edge_types, const QPView &iz, const QPView &ix)
+      : n_edges(n_edges), n_points(n_points), element_index(element_index),
+        edge_index(edge_index), edge_types(edge_types), iz(iz), ix(ix) {}
 
-    /**
-     * @brief Construct edge view with allocated storage.
-     *
-     * @param label Base label for Kokkos view names
-     * @param n_edges Number of edges to allocate
-     * @param n_points Number of quadrature points per edge
-     */
-    EdgeView(const std::string &label, const int n_edges, const int n_points)
-        : n_edges(n_edges), n_points(n_points),
-          element_index(label + "_element_index", n_edges),
-          edge_index(label + "_edge_index", n_edges),
-          edge_types(label + "_edge_types", n_edges),
-          iz(label + "_iz", n_edges, n_points),
-          ix(label + "_ix", n_edges, n_points) {}
+  /**
+   * @brief Access individual edge by index.
+   *
+   * @param edge_id Index of the edge to access
+   * @return Edge structure for the specified edge
+   */
+  KOKKOS_INLINE_FUNCTION
+  Edge<ExecutionSpace> operator()(const int edge_id) const {
+    return { n_points,
+              element_index(edge_id),
+              edge_index(edge_id),
+              edge_types(edge_id),
+              Kokkos::subview(iz, edge_id, Kokkos::ALL()),
+              Kokkos::subview(ix, edge_id, Kokkos::ALL()) };
+  }
 
-    IndexView element_index; ///< Element indices for each edge
-    IndexView edge_index;    ///< Local edge indices within elements
-    EdgeTypeView edge_types; ///< 2D edge type classifications
-    QPView iz;               ///< Z-direction quadrature indices for all edges
-    QPView ix;               ///< X-direction quadrature indices for all edges
+  /**
+   * @brief Access subrange of edges.
+   *
+   * @param edge_range Pair specifying start and end indices
+   * @return EdgeView containing the specified range of edges
+   */
+  KOKKOS_INLINE_FUNCTION
+  EdgeView<ExecutionSpace>
+  operator()(const Kokkos::pair<int, int> &edge_range) const {
+    return { edge_range.second - edge_range.first,
+              n_points,
+              Kokkos::subview(element_index, edge_range),
+              Kokkos::subview(edge_index, edge_range),
+              Kokkos::subview(edge_types, edge_range),
+              Kokkos::subview(iz, edge_range, Kokkos::ALL()),
+              Kokkos::subview(ix, edge_range, Kokkos::ALL()) };
+  }
+};
 
-    /**
-     * @brief Device-side constructor from existing views.
-     *
-     * @param n_edges Number of edges
-     * @param n_points Number of quadrature points per edge
-     * @param element_index Element indices view
-     * @param edge_index Edge indices view
-     * @param edge_types Edge types view
-     * @param iz Z-direction quadrature indices
-     * @param ix X-direction quadrature indices
-     */
-    KOKKOS_INLINE_FUNCTION
-    EdgeView(const int n_edges, const int n_points,
-             const IndexView &element_index, const IndexView &edge_index,
-             const EdgeTypeView &edge_types, const QPView &iz, const QPView &ix)
-        : n_edges(n_edges), n_points(n_points), element_index(element_index),
-          edge_index(edge_index), edge_types(edge_types), iz(iz), ix(ix) {}
 
-    /**
-     * @brief Access individual edge by index.
-     *
-     * @param edge_id Index of the edge to access
-     * @return Edge structure for the specified edge
-     */
-    KOKKOS_INLINE_FUNCTION
-    Edge<ExecutionSpace> operator()(const int edge_id) const {
-      return { n_points,
-               element_index(edge_id),
-               edge_index(edge_id),
-               edge_types(edge_id),
-               Kokkos::subview(iz, edge_id, Kokkos::ALL()),
-               Kokkos::subview(ix, edge_id, Kokkos::ALL()) };
-    }
+/**
+ * @brief 2D spectral element edge classification and coupling management
+ *
+ * This template specialization provides storage and management for edge
+ * information in 2D spectral element meshes. It handles edge connectivity,
+ * interface types, and boundary conditions for coupling between different
+ * media types in 2D wave propagation problems.
+ *
+ * @code
+ * // Construct 2D edge types from mesh data
+ * specfem::assembly::edge_types<specfem::element::dimension_tag::dim2> edges(
+ *     ngllx, ngllz, mesh, element_types);
+ *
+ * // Get elastic-acoustic coupling edges on device
+ * auto [self_edges, coupled_edges] = edges.get_edges_on_device(
+ *     specfem::connections::type::weakly_conforming,
+ *     specfem::interface::interface_tag::elastic_acoustic,
+ *     specfem::element::boundary_tag::none);
+ * @endcode
+ */
+template <> struct edge_types<specfem::element::dimension_tag::dim2> {
+
+public:
+  constexpr static auto dimension_tag =
+      specfem::element::dimension_tag::dim2; ///< Dimension tag
 
-    /**
-     * @brief Access subrange of edges.
-     *
-     * @param edge_range Pair specifying start and end indices
-     * @return EdgeView containing the specified range of edges
-     */
-    KOKKOS_INLINE_FUNCTION
-    EdgeView<ExecutionSpace>
-    operator()(const Kokkos::pair<int, int> &edge_range) const {
-      return { edge_range.second - edge_range.first,
-               n_points,
-               Kokkos::subview(element_index, edge_range),
-               Kokkos::subview(edge_index, edge_range),
-               Kokkos::subview(edge_types, edge_range),
-               Kokkos::subview(iz, edge_range, Kokkos::ALL()),
-               Kokkos::subview(ix, edge_range, Kokkos::ALL()) };
-    }
-  };
 
 public:
   /**

core/specfem/assembly/face_types.hpp

@@ -0,0 +1,28 @@
+#pragma once
+
+#include "enumerations/interface.hpp"
+
+namespace specfem::assembly {
+
+/**
+ * @brief Spectral element face classification and coupling management
+ *
+ * This class provides storage and management for face information in
+ * spectral element meshes, including face connectivity, interface types,
+ * and boundary conditions. It handles coupling between different media
+ * types and manages face-based operations essential for discontinuous
+ * Galerkin (DG) and coupled field formulations in 3D.
+ *
+ * The class manages faces that connect spectral elements, including:
+ * - Conforming and non-conforming interfaces
+ * - Elastic-acoustic coupling interfaces
+ * - Free surface and absorbing boundary conditions
+ * - Mortar element connections for non-matching meshes
+ *
+ * @tparam DimensionTag Spatial dimension (3D)
+ */
+template <specfem::element::dimension_tag DimensionTag> struct face_types;
+
+} // namespace specfem::assembly
+
+#include "face_types/dim3/face_types.hpp"