vertex based patch creation

src/CMakeLists.txt

@@ -71,6 +71,7 @@ set(Omega_h_SOURCES
   Omega_h_owners.cpp
   Omega_h_parser.cpp
   Omega_h_parser_graph.cpp
+  Omega_h_patches.cpp
   Omega_h_pool.cpp
   Omega_h_print.cpp
   Omega_h_profile.cpp
@@ -394,6 +395,12 @@ if(BUILD_TESTING)
   if (Omega_h_USE_Kokkos)
     osh_add_exe(bbox_reduce_test)
     osh_add_exe(initKokkosAndLib)
+    osh_add_exe(test_patches)
+    if(Omega_h_USE_MPI)
+      test_func(run_test_patches_par 4 ./test_patches)
+    else()
+      test_func(run_test_patches 1 ./test_patches)
+    endif()
   endif()
   set(TEST_EXES ${TEST_EXES} unit_math)
   test_basefunc(run_unit_math 1 ./unit_math)

src/Omega_h_mesh.hpp

@@ -167,6 +167,23 @@ class Mesh {
       std::string const& name, Read<T> array);
   friend class ScopedChangeRCFieldsToMesh;
 
+  #if defined(OMEGA_H_USE_KOKKOS)
+  /**
+   * \brief form a patch of at least minPatchSize elements surrounding each mesh
+   *        vertex
+   * \remark the patch is expanded via 2nd order adjacencies using meshDim-1 as
+   *         the bridge entity (e.g., faces for 3d, edges for 2d)
+   * \param m (in) mesh of simplices
+   * \param minPatchSize (in) the minimum number of elements in each patch
+   * \return a graph whose source nodes are mesh vertices, and
+   *         edges are connecting to elements in the patch 
+   *         OR
+   *         an empty graph upon failure
+   */
+  [[nodiscard]] Graph get_vtx_patches(Int minPatchSize);
+  #endif
+
+
  private:
   bool change_all_rcFieldsToMesh();
   bool change_all_rcFieldsTorc();

src/Omega_h_patches.cpp

@@ -0,0 +1,119 @@
+#include <Omega_h_mesh.hpp>
+#include <Omega_h_for.hpp> //parallel_for
+#include <Omega_h_array_ops.hpp> //get_min
+#include <Omega_h_int_scan.hpp> //offset_scan
+
+using namespace Omega_h;
+
+#if defined(OMEGA_H_USE_KOKKOS)
+#include <Kokkos_NestedSort.hpp> //sort_team
+
+namespace {
+[[nodiscard]] Graph adj_segment_sort(Graph& g) {
+  using ExecSpace = Kokkos::DefaultExecutionSpace;
+  using TeamPol = Kokkos::TeamPolicy<ExecSpace>;
+  using TeamMem = typename TeamPol::member_type;
+  auto offsets = g.a2ab;
+  auto elms_r = g.ab2b.view(); //read only
+  Kokkos::View<LO*, ExecSpace> elms("elms", elms_r.size());
+  Kokkos::deep_copy(elms, elms_r);
+  auto segment_sort = KOKKOS_LAMBDA(const TeamMem& t) {
+    auto i = t.league_rank();
+    auto patch = Kokkos::subview(elms, Kokkos::make_pair(offsets[i], offsets[i+1]));
+    Kokkos::Experimental::sort_team(t, patch);
+  };
+  Kokkos::parallel_for(TeamPol(g.nnodes(), Kokkos::AUTO()), segment_sort);
+  return Graph(offsets,Write<LO>(elms));
+}
+
+[[nodiscard]] Graph remove_duplicate_edges(Graph g) {
+  auto offsets = g.a2ab;
+  auto values = g.ab2b;
+  Write<I8> keep(g.nedges());
+  auto markDups = OMEGA_H_LAMBDA(LO i) {
+    keep[offsets[i]] = 1; //always keep the first edge in the segment
+    for(int j=offsets[i]+1; j<offsets[i+1]; j++) {
+      keep[j] = (values[j-1] != values[j]);
+    }
+  };
+  parallel_for(g.nnodes(), markDups);
+  auto filtered = filter_graph_edges(g,keep);
+  return filtered;
+}
+
+[[nodiscard]] Read<I8> patchSufficient(Graph patches, Int minPatchSize) {
+  const auto num_patches = patches.nnodes();
+  auto offsets = patches.a2ab;
+  Write<I8> done(num_patches);
+  parallel_for(num_patches, OMEGA_H_LAMBDA(LO i) {
+    done[i] = ((offsets[i+1]-offsets[i]) >= minPatchSize);
+  });
+  return read(done);
+}
+
+/**
+ * \brief expand the patches
+ * \param m (in) mesh of simplices
+ * \param patches (in) graph of key entities to elements
+ * \param adjElms (in) second order element-to-element adjacencies
+ *                     used for expansion
+ * \return an expanded graph from key entities to elements
+*/
+//TODO use Omega_h_map and Omega_h_graph functions
+[[nodiscard]] Graph expandPatches(Mesh& m, Graph patches, Graph adjElms, Read<I8> patchDone) {
+  auto adjElms_offsets = adjElms.a2ab;
+  auto adjElms_elms = adjElms.ab2b;
+  const auto num_patches = patches.nnodes();
+  auto patch_offsets = patches.a2ab;
+  auto patch_elms = patches.ab2b;
+  Write<LO> degree(num_patches);
+  parallel_for(num_patches, OMEGA_H_LAMBDA(LO patch) {
+    degree[patch] = patch_offsets[patch+1] - patch_offsets[patch];
+    if(!patchDone[patch]) {
+      for(int j=patch_offsets[patch]; j<patch_offsets[patch+1]; j++) {
+        auto elm = patch_elms[j];
+        degree[patch] += adjElms_offsets[elm+1]-adjElms_offsets[elm]; //counts duplicates
+      }
+    }
+  });
+  auto patchExpDup_offsets = offset_scan(read(degree));
+  Write<LO> patchExpDup_elms(patchExpDup_offsets.last());
+  parallel_for(num_patches, OMEGA_H_LAMBDA(LO patch) {
+    auto idx = patchExpDup_offsets[patch];
+    for(int j=patch_offsets[patch]; j<patch_offsets[patch+1]; j++) {
+      patchExpDup_elms[idx++] = patch_elms[j];
+      if(!patchDone[patch]) {
+        auto elm = patch_elms[j];
+        for(int k=adjElms_offsets[elm]; k<adjElms_offsets[elm+1]; k++) {
+          patchExpDup_elms[idx++] = adjElms_elms[k];
+        }
+      }
+    }
+  });
+  Graph patchExpDup(patchExpDup_offsets,patchExpDup_elms);
+  auto sorted = adj_segment_sort(patchExpDup);
+  auto dedup = remove_duplicate_edges(sorted);
+  return dedup;
+}
+}//end anonymous namespace
+
+[[nodiscard]] Graph Mesh::get_vtx_patches(Int minPatchSize) {
+  OMEGA_H_CHECK(minPatchSize > 0);
+  auto patches = ask_up(VERT,dim());
+  auto patchDone = patchSufficient(patches, minPatchSize);
+  if( get_min(patchDone) == 1 )
+    return patches;
+  auto adjElms = ask_dual();
+  //assuming each iteration adds at least one element then
+  //the minPatchSize is a conservative upper bound on the
+  //iteration count
+  for(Int iter = 0; iter < minPatchSize; iter++) {
+    patches = expandPatches(*this, patches, adjElms, patchDone);
+    patchDone = patchSufficient(patches, minPatchSize);
+    if( get_min(patchDone) == 1 ) {
+      return patches;
+    }
+  }
+  return Graph();
+}
+#endif

src/test_patches.cpp

@@ -0,0 +1,111 @@
+#include <Omega_h_file.hpp>
+#include <Omega_h_library.hpp>
+#include <Omega_h_mesh.hpp>
+#include <Omega_h_build.hpp> //build_box
+
+using namespace Omega_h;
+
+void test2x2(Omega_h::CommPtr comm) {
+  OMEGA_H_CHECK(comm->size() == 1);
+  const auto x = 2.0;
+  const auto y = 2.0;
+  const auto z = 0.0;
+  const auto nx = 2;
+  const auto ny = 2;
+  const auto nz = 0;
+  const auto symmetric = false;
+  auto mesh = Omega_h::build_box(comm, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric);
+  const auto minPatchSize = 3;
+  auto patches = mesh.get_vtx_patches(minPatchSize);
+  Graph expected(
+   {0,4,7,11,17,20,24,27,30,34},
+   {0,1,2,6,1,2,4,1,2,4,5,0,1,2,3,5,6,1,4,5,1,3,5,6,3,6,7,0,6,7,0,3,6,7});
+  OMEGA_H_CHECK(patches == expected);
+}
+
+void test1x5(Omega_h::CommPtr comm) {
+  OMEGA_H_CHECK(comm->size() == 1);
+  const auto x = 1.0;
+  const auto y = 5.0;
+  const auto z = 0.0;
+  const auto nx = 1;
+  const auto ny = 5;
+  const auto nz = 0;
+  const auto symmetric = false;
+  auto mesh = Omega_h::build_box(comm, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric);
+  {
+    const auto minPatchSize = 3;
+    auto patches = mesh.get_vtx_patches(minPatchSize);
+    Graph expected(
+      {0,3,6,9,12,15,18,21,24,27,30,33,36},
+      {0,1,2,1,2,3,0,1,2,0,1,2,1,3,4,3,4,6,5,6,9,4,5,6,7,8,9,7,8,9,7,8,9,5,7,9});
+    OMEGA_H_CHECK(patches == expected);
+  }
+  {
+    const auto minPatchSize = 4;
+    auto patches = mesh.get_vtx_patches(minPatchSize);
+    Graph expected(
+       {0,4,9,13,17,22,27,32,37,41,45,49,54},
+       {0,1,2,3,0,1,2,3,4,0,1,2,3,0,1,2,3,1,2,3,4,6,1,3,
+        4,5,6,4,5,6,7,9,3,4,5,6,9,5,7,8,9,5,7,8,9,5,7,8,9,5,6,7,8,9});
+    OMEGA_H_CHECK(patches == expected);
+  }
+}
+
+void test3D(Omega_h::CommPtr comm) {
+  OMEGA_H_CHECK(comm->size() == 1);
+  const auto x = 1.0;
+  const auto y = 1.0;
+  const auto z = 1.0;
+  const auto nx = 2;
+  const auto ny = 2;
+  const auto nz = 2;
+  const auto symmetric = false;
+  auto mesh = Omega_h::build_box(comm, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric);
+  const auto minPatchSize = 3;
+  auto patches = mesh.get_vtx_patches(minPatchSize);
+  OMEGA_H_CHECK(patches.nnodes() != -1);
+}
+
+void testPar(Omega_h::CommPtr comm) {
+  OMEGA_H_CHECK(comm->size() == 4);
+  const auto x = 1.0;
+  const auto y = 1.0;
+  const auto z = 0.0;
+  const auto nx = 4;
+  const auto ny = 4;
+  const auto nz = 0;
+  const auto symmetric = false;
+  auto mesh = Omega_h::build_box(comm, OMEGA_H_SIMPLEX, x, y, z, nx, ny, nz, symmetric);
+  mesh.set_parting(OMEGA_H_GHOSTED);
+  {
+    const auto min_elems = comm->allreduce(mesh.nelems(), OMEGA_H_MIN);
+    const auto minPatchSize = min_elems+1;
+    auto patches = mesh.get_vtx_patches(minPatchSize);
+    if( mesh.nelems() < minPatchSize ) {
+      OMEGA_H_CHECK(patches.nnodes() == -1); //expected to fail
+    } else {
+      OMEGA_H_CHECK(patches.nnodes() != -1); //expected to pass
+    }
+  }
+  {
+    const auto minPatchSize = 6;
+    auto patches = mesh.get_vtx_patches(minPatchSize);
+    OMEGA_H_CHECK(patches.nnodes() != -1);
+  }
+}
+
+int main(int argc, char** argv) {
+  auto lib = Omega_h::Library(&argc, &argv);
+  auto world = lib.world();
+  OMEGA_H_CHECK(argc == 1);
+  if (world->rank() == 0) {
+    test2x2(lib.self());
+    test1x5(lib.self());
+    test3D(lib.self());
+  }
+  if (world->size() == 4) {
+    testPar(world);
+  }
+  return 0;
+}