Conservative transfer 3D

CMakeLists.txt

@@ -49,6 +49,15 @@ option(PCMS_ENABLE_Python "Enable pcms Python api" OFF)
 
 option(PCMS_ENABLE_PRINT "PCMS print statements enabled" ON)
 
+# Per-thread on-stack ring-buffer capacities for the adjacency-based
+# intersection/adjacency BFS (queue_visited.hpp). Defaults are sized for 3D
+# tetrahedral meshes on the serial/OpenMP backends. Lower them for device
+# (GPU) builds, where each search thread's stack must hold these arrays.
+set(PCMS_INTERSECTION_QUEUE_SIZE 1024 CACHE STRING
+    "Max BFS queue size per target element in the intersection search")
+set(PCMS_INTERSECTION_TRACK_SIZE 2048 CACHE STRING
+    "Max visited-set size per target element in the intersection search")
+
 option(PETSC_LINK_STATIC "Use pkg-config --static results for PETSc"
        ${_pcms_link_petsc_static_default})
 

src/pcms/configuration.h.in

@@ -7,3 +7,6 @@
 #cmakedefine PCMS_ENABLE_Fortran
 #cmakedefine PCMS_ENABLE_MESHFIELDS
 #cmakedefine PCMS_ENABLE_PETSC
+
+#cmakedefine PCMS_INTERSECTION_QUEUE_SIZE @PCMS_INTERSECTION_QUEUE_SIZE@
+#cmakedefine PCMS_INTERSECTION_TRACK_SIZE @PCMS_INTERSECTION_TRACK_SIZE@

src/pcms/field/layout/omega_h_lagrange.cpp

@@ -171,6 +171,7 @@ OmegaHLagrangeLayout::OmegaHLagrangeLayout(
   owned_ = BuildOwned(mesh_, entity_dim, owned_mask);
   owned_host_ =
     Kokkos::View<bool*, HostMemorySpace>("owned_host", owned_.size());
+  Kokkos::deep_copy(owned_host_, owned_);
 
   class_ids_ = Omega_h::Read<Omega_h::ClassId>(
     mesh_.get_array<Omega_h::ClassId>(entity_dim, "class_id"));

src/pcms/localization/queue_visited.hpp

@@ -7,8 +7,23 @@
 #include <Omega_h_library.hpp>
 #include <Omega_h_mesh.hpp>
 #include <Omega_h_reduce.hpp>
-#define MAX_SIZE_QUEUE 500
-#define MAX_SIZE_TRACK 800
+#include "pcms/configuration.h"
+// Per-thread BFS ring-buffer capacities for the adjacency-based intersection
+// search. Sized for 3D: a target tetrahedron can overlap many more source
+// elements than a 2D triangle, and the visited set additionally holds the
+// non-intersecting neighbors probed along the way. These are on-stack arrays,
+// so raising them raises per-thread stack usage; lower them for device (GPU)
+// builds via the PCMS_INTERSECTION_{QUEUE,TRACK}_SIZE CMake cache variables,
+// which flow in through pcms/configuration.h. The fallbacks below apply only if
+// the configured header is unavailable (e.g. a header-only consumer).
+#ifndef PCMS_INTERSECTION_QUEUE_SIZE
+#define PCMS_INTERSECTION_QUEUE_SIZE 1024
+#endif
+#ifndef PCMS_INTERSECTION_TRACK_SIZE
+#define PCMS_INTERSECTION_TRACK_SIZE 2048
+#endif
+#define MAX_SIZE_QUEUE PCMS_INTERSECTION_QUEUE_SIZE
+#define MAX_SIZE_TRACK PCMS_INTERSECTION_TRACK_SIZE
 
 namespace pcms
 {

src/pcms/transfer/mass_matrix_integrator.hpp

@@ -17,16 +17,19 @@ template <typename FieldElement>
 class MassMatrixIntegrator : public MeshField::Integrator
 {
 public:
+  // Linear simplex: numNodes = spatial dim + 1 (3 for triangles, 4 for tets).
+  static constexpr int numNodes = FieldElement::MeshEntDim + 1;
+
   MassMatrixIntegrator(Omega_h::Mesh& mesh_in, FieldElement& fe_in,
                        int order = 2)
     : mesh(mesh_in),
       fe(fe_in),
-      subMatrixSize(3 * 3), // FIXME remove hard coded size
-      elmMassMatrix("elmMassMatrix", mesh_in.nelems() * 3 * 3),
+      subMatrixSize(numNodes * numNodes),
+      elmMassMatrix("elmMassMatrix", mesh_in.nelems() * numNodes * numNodes),
       Integrator(order)
   {
     Kokkos::deep_copy(elmMassMatrix, 0);
-    assert(mesh.dim() == 2); // TODO support 1d,2d,3d
+    assert(mesh.dim() == 2 || mesh.dim() == 3);
     assert(mesh.family() == OMEGA_H_SIMPLEX);
   }
   void atPoints(Kokkos::View<MeshField::Real**> p,
@@ -58,12 +61,18 @@ class MassMatrixIntegrator : public MeshField::Integrator
           }
           const auto N = shapeFn.getValues(localCoord);
           const auto wPt = w(pt);
-          const auto dVPt = dV(pt);
+          // Use the unsigned volume element: MeshField returns a signed
+          // Jacobian determinant, which is negative for tetrahedra whose vertex
+          // ordering has negative orientation. A mass matrix integrates against
+          // the positive volume measure, so take the magnitude (a no-op in 2D
+          // where the differential area is already positive).
+          const auto dVPt = Kokkos::fabs(dV(pt));
           //       printf("Shape Functions: %f, %f, %f \n", N[0], N[1], N[2]);
           //       printf("wPt, dVPt: %f, %f \n", wPt, dVPt);
           for (auto i = 0; i < N.size(); i++) {
             for (auto j = 0; j < N.size(); j++) {
-              massMatrix(elm * subMat + i * 3 + j) += N[i] * N[j] * wPt * dVPt;
+              massMatrix(elm * subMat + i * numNodes + j) +=
+                N[i] * N[j] * wPt * dVPt;
             }
           }
         }

src/pcms/transfer/mesh_intersection.cpp

@@ -4,51 +4,66 @@
 
 namespace pcms
 {
+namespace
+{
+// Construct the source-mesh containing-element search appropriate to the
+// spatial dimension: a uniform 20^Dim background grid over the source mesh.
+template <int Dim>
+auto MakeGridPointSearch(Omega_h::Mesh& source_mesh)
+{
+  if constexpr (Dim == 3) {
+    return pcms::GridPointSearch3D(source_mesh, 20, 20, 20);
+  } else {
+    return pcms::GridPointSearch2D(source_mesh, 20, 20);
+  }
+}
+} // namespace
+
+template <int Dim>
 void FindIntersections::adjBasedIntersectSearch(
   const Omega_h::LOs& tgt2src_offsets,
   Omega_h::Write<Omega_h::LO>& nIntersections,
   Omega_h::Write<Omega_h::LO>& tgt2src_indices, bool is_count_only)
 {
-
+  // Element entity dimension equals the spatial dimension (FACE for 2D, REGION
+  // for 3D); measures are triangle areas (2D) or tet volumes (3D).
   const auto& tgt_coords = target_mesh_.coords();
   const auto& src_coords = source_mesh_.coords();
-  const auto& tgt_faces2nodes =
-    target_mesh_.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b;
-  const auto& src_faces2nodes =
-    source_mesh_.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b;
-  const auto& src_elem_areas = measure_elements_real(&source_mesh_);
-  const auto& tgt_elem_areas = measure_elements_real(&target_mesh_);
+  const auto& tgt_elems2nodes = target_mesh_.ask_down(Dim, Omega_h::VERT).ab2b;
+  const auto& src_elems2nodes = source_mesh_.ask_down(Dim, Omega_h::VERT).ab2b;
+  const auto& src_elem_measures = measure_elements_real(&source_mesh_);
+  const auto& tgt_elem_measures = measure_elements_real(&target_mesh_);
   const auto& t2t =
     source_mesh_.ask_dual(); // gives connected element neighbors
   const auto& t2tt = t2t.a2ab;
   const auto& tt2t = t2t.ab2b;
 
-  const auto flat_centroids =
-    pcms::get_entity_centroids(target_mesh_, Omega_h::FACE);
+  const auto flat_centroids = pcms::get_entity_centroids(target_mesh_, Dim);
   // Convert layout_right 1D Omega_h array to 2D Kokkos view with correct layout
-  auto centroids = ConvertCoordsTo2D(flat_centroids, target_mesh_.nfaces(), 2);
+  auto centroids =
+    ConvertCoordsTo2D(flat_centroids, target_mesh_.nelems(), Dim);
 
-  pcms::GridPointSearch2D search_cell(source_mesh_, 20, 20);
+  auto search_cell = MakeGridPointSearch<Dim>(source_mesh_);
   auto results = search_cell(centroids);
   auto owning_cell_ids = search_cell.GetOwningElementIds(results);
 
-  auto nfaces_target = target_mesh_.nfaces();
+  auto nelems_target = target_mesh_.nelems();
   Omega_h::parallel_for(
-    nfaces_target,
+    nelems_target,
     OMEGA_H_LAMBDA(const Omega_h::LO id) {
       Queue queue;
       Track visited;
 
       auto current_cell_id = owning_cell_ids(id);
-      auto current_tgt_elm_area = tgt_elem_areas[id];
+      auto current_tgt_elm_measure = tgt_elem_measures[id];
 
       OMEGA_H_CHECK_PRINTF(current_cell_id >= 0,
                            "ERROR: source cell id not found for given target "
-                           "centroid %d (%f, %f)\n",
-                           id, centroids(id, 0), centroids(id, 1));
+                           "centroid %d\n",
+                           id);
 
       auto tgt_elm_vert_coords =
-        get_vert_coords_of_elem(tgt_coords, tgt_faces2nodes, id);
+        get_vert_coords_of_elem<Dim>(tgt_coords, tgt_elems2nodes, id);
 
       Omega_h::LO start_counter;
       if (!is_count_only) {
@@ -76,22 +91,42 @@ void FindIntersections::adjBasedIntersectSearch(
           auto neighborElmId = tt2t[i];
 
           if (visited.notVisited(neighborElmId)) {
-            visited.push_back(neighborElmId);
-            auto elm_vert_coords = get_vert_coords_of_elem(
-              src_coords, src_faces2nodes, neighborElmId);
-            r3d::Polytope<2> intersection;
+            // If the visited buffer is full, skip this neighbor so the BFS
+            // terminates. Without this, an unrecorded neighbor stays "not
+            // visited" and is re-queued forever (infinite loop). Mirrors the
+            // guard in adj_search.cpp. Raise MAX_SIZE_TRACK if this fires.
+            if (!visited.push_back(neighborElmId)) {
+              printf("ERROR: visited buffer full (MAX_SIZE_TRACK=%d) for "
+                     "target %d; some intersections may be missed\n",
+                     MAX_SIZE_TRACK, id);
+              continue;
+            }
+            auto elm_vert_coords = get_vert_coords_of_elem<Dim>(
+              src_coords, src_elems2nodes, neighborElmId);
+            r3d::Polytope<Dim> intersection;
             r3d::intersect_simplices(intersection, tgt_elm_vert_coords,
                                      elm_vert_coords);
-            auto intersected_area = r3d::measure(intersection);
-            auto current_src_elm_area = src_elem_areas[neighborElmId];
+            // Take the magnitude: r3d::measure is signed by the orientation of
+            // the target simplex used to initialize the polytope, which can be
+            // negative for tetrahedra. This mirrors the fabs applied in the
+            // sub-simplex decomposition and mass assembly; without it a
+            // negatively-oriented target element would reject all of its real
+            // overlaps and break conservation.
+            auto intersected_measure = Kokkos::fabs(r3d::measure(intersection));
+            auto current_src_elm_measure = src_elem_measures[neighborElmId];
             auto scale =
-              Kokkos::fmax(current_tgt_elm_area, current_src_elm_area);
+              Kokkos::fmax(current_tgt_elm_measure, current_src_elm_measure);
             auto eps = Kokkos::fmax(abs_tol, rel_tol * scale);
-            if (intersection.nverts >= 3 && intersected_area >= eps) {
+            // A valid intersection is a non-degenerate simplex-simplex overlap:
+            // at least Dim+1 vertices (a polygon in 2D, a polyhedron in 3D).
+            if (intersection.nverts >= Dim + 1 && intersected_measure >= eps) {
               count++;
 
               OMEGA_H_CHECK_PRINTF(
-                count < 500, "WARNING: count exceeds 500 for target %d", id);
+                count < kMaxIntersectionsPerTarget,
+                "intersection count for target %d reached the cap %d; raise "
+                "kMaxIntersectionsPerTarget/MAX_SIZE_QUEUE",
+                id, kMaxIntersectionsPerTarget);
 
               queue.push_back(neighborElmId);
 
@@ -113,20 +148,32 @@ void FindIntersections::adjBasedIntersectSearch(
     }, // end of lambda
     "count the number of intersections for each target element");
 }
-IntersectionResults intersectTargets(Omega_h::Mesh& source_mesh,
-                                     Omega_h::Mesh& target_mesh)
+
+// Explicit instantiations for the supported spatial dimensions.
+template void FindIntersections::adjBasedIntersectSearch<2>(
+  const Omega_h::LOs&, Omega_h::Write<Omega_h::LO>&,
+  Omega_h::Write<Omega_h::LO>&, bool);
+template void FindIntersections::adjBasedIntersectSearch<3>(
+  const Omega_h::LOs&, Omega_h::Write<Omega_h::LO>&,
+  Omega_h::Write<Omega_h::LO>&, bool);
+
+namespace
+{
+template <int Dim>
+IntersectionResults intersectTargetsImpl(Omega_h::Mesh& source_mesh,
+                                         Omega_h::Mesh& target_mesh)
 {
   FindIntersections intersect(source_mesh, target_mesh);
 
-  auto nfaces_target = target_mesh.nfaces();
+  auto nelems_target = target_mesh.nelems();
 
   Omega_h::Write<Omega_h::LO> nIntersections(
-    nfaces_target, 0, "number of intersections in each target vertex");
+    nelems_target, 0, "number of intersections in each target element");
 
   Omega_h::Write<Omega_h::LO> tgt2src_indices;
 
-  intersect.adjBasedIntersectSearch(Omega_h::LOs(), nIntersections,
-                                    tgt2src_indices, true);
+  intersect.adjBasedIntersectSearch<Dim>(Omega_h::LOs(), nIntersections,
+                                         tgt2src_indices, true);
 
   Kokkos::fence();
   auto tgt2src_offsets = Omega_h::offset_scan(Omega_h::Read(nIntersections),
@@ -139,9 +186,20 @@ IntersectionResults intersectTargets(Omega_h::Mesh& source_mesh,
     ntotal_intersections, 0,
     "indices of the source elements that intersect the given target element");
 
-  intersect.adjBasedIntersectSearch(tgt2src_offsets, nIntersections,
-                                    tgt2src_indices, false);
+  intersect.adjBasedIntersectSearch<Dim>(tgt2src_offsets, nIntersections,
+                                         tgt2src_indices, false);
   return {.tgt2src_offsets = tgt2src_offsets,
           .tgt2src_indices = Omega_h::read(tgt2src_indices)};
 }
+} // namespace
+
+IntersectionResults intersectTargets(Omega_h::Mesh& source_mesh,
+                                     Omega_h::Mesh& target_mesh)
+{
+  OMEGA_H_CHECK(source_mesh.dim() == target_mesh.dim());
+  if (source_mesh.dim() == 3) {
+    return intersectTargetsImpl<3>(source_mesh, target_mesh);
+  }
+  return intersectTargetsImpl<2>(source_mesh, target_mesh);
+}
 } // namespace pcms

src/pcms/transfer/mesh_intersection.hpp

@@ -15,19 +15,29 @@ namespace pcms
 constexpr static double abs_tol = 1e-18; /// abs tolerance
 constexpr static double rel_tol = 1e-12; /// rel tolerance
 
-[[nodiscard]] OMEGA_H_INLINE r3d::Few<r3d::Vector<2>, 3>
+// Upper bound on intersecting source elements recorded per target element.
+// Must not exceed the BFS queue capacity (MAX_SIZE_QUEUE); raised for 3D where
+// a single target tet can overlap many source tets.
+constexpr static int kMaxIntersectionsPerTarget = MAX_SIZE_QUEUE;
+
+// Gather the vertex coordinates of a simplex element (triangle for Dim==2,
+// tetrahedron for Dim==3) into an r3d simplex, ready for
+// r3d::intersect_simplices.
+template <int Dim>
+[[nodiscard]] OMEGA_H_INLINE r3d::Few<r3d::Vector<Dim>, Dim + 1>
 get_vert_coords_of_elem(const Omega_h::Reals& coords,
-                        const Omega_h::LOs& faces2nodes, const int id)
+                        const Omega_h::LOs& elems2nodes, const int id)
 {
-  const auto elm_verts = Omega_h::gather_verts<3>(faces2nodes, id);
+  const auto elm_verts = Omega_h::gather_verts<Dim + 1>(elems2nodes, id);
 
-  const Omega_h::Matrix<2, 3> elm_vert_coords =
-    Omega_h::gather_vectors<3, 2>(coords, elm_verts);
+  const Omega_h::Matrix<Dim, Dim + 1> elm_vert_coords =
+    Omega_h::gather_vectors<Dim + 1, Dim>(coords, elm_verts);
 
-  r3d::Few<r3d::Vector<2>, 3> r3d_vector;
-  for (int i = 0; i < 3; ++i) {
-    r3d_vector[i][0] = elm_vert_coords[i][0];
-    r3d_vector[i][1] = elm_vert_coords[i][1];
+  r3d::Few<r3d::Vector<Dim>, Dim + 1> r3d_vector;
+  for (int i = 0; i < Dim + 1; ++i) {
+    for (int d = 0; d < Dim; ++d) {
+      r3d_vector[i][d] = elm_vert_coords[i][d];
+    }
   }
 
   return r3d_vector;
@@ -78,11 +88,13 @@ class FindIntersections
    * @param is_count_only If true, only counts intersections; if false, also
    * fills tgt2src_indices.
    *
-   * @note This method assumes 2D linear triangles and uses
-   * `r3d::intersect_simplices` for geometric intersection.
+   * @note Templated on spatial dimension `Dim`: linear triangles (Dim==2) or
+   * linear tetrahedra (Dim==3), using `r3d::intersect_simplices` for geometric
+   * intersection.
    *
    * @see r3d::intersect_simplices, intersectTargets
    */
+  template <int Dim>
   void adjBasedIntersectSearch(const Omega_h::LOs& tgt2src_offsets,
                                Omega_h::Write<Omega_h::LO>& nIntersections,
                                Omega_h::Write<Omega_h::LO>& tgt2src_indices,