Added BVH accelerated point localization and updated localization API

README.md

@@ -91,6 +91,14 @@ cmake -S redev -B build-redev \
   -DADIOS2_DIR=$PWD/build-ADIOS2/install/lib64/cmake/adios2 \
   -Dperfstubs_DIR=$PWD/build-perfstubs/install/lib/cmake
 cmake --build build-redev --target install
+
+git clone git@github.com:arborx/ArborX.git
+cmake -S ArborX -B build-ArborX -D CMAKE_INSTALL_PREFIX="$PWD/build-ArborX/install" \
+-D ARBORX_ENABLE_MPI=ON \
+-D Kokkos_ROOT="$PWD/build-kokkos/install" \
+-D CMAKE_CXX_COMPILER=mpicxx \
+-D CMAKE_CXX_EXTENSIONS=OFF 
+cmake --build build-ArborX --target install
 ```
 ### Build dependencies with CUDA
 ```
@@ -164,6 +172,14 @@ git clone git@github.com:catchorg/Catch2
 cmake -S Catch2 -B build-Catch2 \
   -DCMAKE_INSTALL_PREFIX=$PWD/build-Catch2/install
 cmake --build build-Catch2 --target install
+
+git clone git@github.com:arborx/ArborX.git
+cmake -S ArborX -B build-ArborX -D CMAKE_INSTALL_PREFIX="$PWD/build-ArborX/install" \
+-D ARBORX_ENABLE_MPI=ON \
+-D Kokkos_ROOT="$PWD/build-kokkos/install" \
+-D CMAKE_CXX_COMPILER="$PWD/build-kokkos/install/bin/nvcc_wrapper" \
+-D CMAKE_CXX_EXTENSIONS=OFF
+cmake --build build-ArborX --target install
 ```
 ### Build, install, and test pcms
 ```
@@ -181,6 +197,7 @@ cmake -S pcms -B build-pcms \
   -Dperfstubs_DIR=$PWD/build-perfstubs/install/lib/cmake \
   -DCatch2_DIR=$PWD/build-Catch2/install/lib64/cmake/Catch2/ \
   -DKokkosKernels_DIR=$PWD/build-kokkos-kernels/install/lib64/cmake/KokkosKernels/ \
+  -DArborX_DIR=$PWD/build-ArborX/install/share/cmake/ArborX \
   -DPCMS_TEST_DATA_DIR=$PWD/pcms_testcases
 cmake --build build-pcms -j 8
 

src/pcms/field/evaluator/omega_h_lagrange.h

@@ -65,21 +65,19 @@ struct CopyCoordsFunctor
 template <int Dim>
 OmegaHLagrangeLocHint BuildLagrangeLocHint(
   Omega_h::Mesh& mesh, int mesh_dim,
-  Kokkos::View<typename PointLocalizationSearch<Dim>::Result*,
-               DeviceMemorySpace>
-    results,
+  PointSearch::Results& results,
   Kokkos::View<Real**, DeviceMemorySpace> coords_d,
   Kokkos::View<LO*, DeviceMemorySpace> owning_ids, OutOfBoundsMode mode)
 {
-  LO n = static_cast<LO>(results.size());
+  LO n = static_cast<LO>(results.dimensionalities.size());
 
   // First pass: count valid and missing
   Kokkos::View<int*, DeviceMemorySpace> is_valid("is_valid", n);
   Kokkos::parallel_for(
     "CheckValidity",
     Kokkos::RangePolicy<typename DeviceMemorySpace::execution_space>(0, n),
     KOKKOS_LAMBDA(LO i) {
-      bool out = (owning_ids(i) < 0) || (results(i).element_id < 0);
+      bool out = (owning_ids(i) < 0) || (results.element_ids(i) < 0);
       is_valid(i) = out ? 0 : 1;
     });
 
@@ -344,7 +342,7 @@ class OmegaHLagrangeEvaluatorFactory : public FieldEvaluatorFactory<T>
         detail::CopyCoordsFunctor<Dim> copy_functor(coords_d, raw_coords);
         Kokkos::parallel_for("copy_coords", n_pts, copy_functor);
 
-        auto results_d = search(coords_d);
+        auto results_d = search.apply(coords);
         auto owning_ids = search.GetOwningElementIds(results_d);
 
         return detail::BuildLagrangeLocHint<Dim>(

src/pcms/field/uniform_grid_binary_field.h

@@ -47,7 +47,7 @@ CreateUniformGridBinaryField(Omega_h::Mesh& mesh, const UniformGrid<Dim>& grid)
   auto coords = coord_view.GetValues();
   LO n = layout->GetNumOwnedDofHolder();
 
-  Kokkos::View<Real* [Dim]> coords_d("coords_d", n);
+  Kokkos::View<Real**> coords_d("coords_d", n, Dim);
   Kokkos::parallel_for(
     "CopyCoords", Kokkos::RangePolicy<>(0, n), KOKKOS_LAMBDA(int i) {
       for (int d = 0; d < Dim; ++d) {
@@ -56,21 +56,21 @@ CreateUniformGridBinaryField(Omega_h::Mesh& mesh, const UniformGrid<Dim>& grid)
     });
 
   // Run point-in-mesh search
-  Kokkos::View<typename PointLocalizationSearch<Dim>::Result*> results_d;
+  PointSearch::Results results_d;
   if constexpr (Dim == 2) {
     GridPointSearch2D search(mesh, grid.divisions[0], grid.divisions[1]);
-    results_d = search(coords_d);
+    results_d = search.apply(pcms::CoordinateView(pcms::CoordinateSystem::Cartesian, pcms::MakeConstRank2View(coords_d)));
   } else {
     GridPointSearch3D search(mesh, grid.divisions[0], grid.divisions[1],
                              grid.divisions[2]);
-    results_d = search(coords_d);
+    results_d = search.apply(pcms::CoordinateView(pcms::CoordinateSystem::Cartesian, pcms::MakeConstRank2View(coords_d)));
   }
-  auto results_h =
-    Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, results_d);
+  auto result_ids_h =
+    Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, results_d.element_ids);
 
   Kokkos::View<Real*, HostMemorySpace> data("binary_mask", n);
   for (LO i = 0; i < n; ++i)
-    data(i) = (results_h(i).element_id >= 0) ? 1.0 : 0.0;
+    data(i) = (result_ids_h(i) >= 0) ? 1.0 : 0.0;
 
   field.SetDOFHolderDataHost(
     Rank2View<const Real, HostMemorySpace>(data.data(), n, 1));

src/pcms/localization/CMakeLists.txt

@@ -7,6 +7,7 @@ set(
         localization_factory.h
         point_cloud_localization.h
         mesh_localization.h
+        point_localization.h
 )
 
 set(
@@ -15,7 +16,9 @@ set(
         adj_search.cpp
         mls_support_helpers.cpp
         mesh_localization.cpp
+        point_localization.cpp
 )
+find_package(ArborX REQUIRED)
 add_library(pcms_localization ${PCMS_LOCALIZATION_SOURCES})
 add_library(pcms::localization ALIAS pcms_localization)
 target_include_directories(
@@ -32,6 +35,7 @@ target_link_libraries(
         pcms::utility
         pcms::discretization
         Omega_h::omega_h
+        ArborX::ArborX
 )
 target_compile_features(pcms_localization PUBLIC cxx_std_20)
 

src/pcms/localization/adj_search.cpp

@@ -7,7 +7,7 @@ namespace pcms
 // Debug helper retained intentionally: useful for diagnosing point-localization
 // failures while developing support-search logic.
 [[maybe_unused]] static void checkTargetPoints(
-  const Kokkos::View<pcms::GridPointSearch2D::Result*>& results,
+  const GridPointSearch2D::Results& results,
   const Kokkos::View<pcms::LO*>& owning_cell_ids)
 {
   Kokkos::fence();
@@ -22,7 +22,7 @@ namespace pcms
       printf("%d, ", i);
     }
   };
-  Omega_h::parallel_for(results.size(), check_target_points,
+  Omega_h::parallel_for(results.dimensionalities.size(), check_target_points,
                         "check_target_points");
   Kokkos::fence();
   pcms::printInfo("\n");
@@ -81,17 +81,19 @@ static Omega_h::Write<Omega_h::LO> locate_target_cells(
     Omega_h::Write<Omega_h::LO>(nvertices_target, -1, "source cell ids");
 
   if (dim == 2) {
-    Kokkos::View<Omega_h::Real* [2]> target_points("target_points",
-                                                   nvertices_target);
+    Kokkos::View<pcms::Real**> target_points("target_points",
+                                                   nvertices_target, 2);
     Omega_h::parallel_for(
       nvertices_target, OMEGA_H_LAMBDA(const Omega_h::LO i) {
         target_points(i, 0) = target_coords[i * dim];
         target_points(i, 1) = target_coords[i * dim + 1];
       });
     Kokkos::fence();
 
-    pcms::GridPointSearch2D search_cell(source_mesh, 10, 10);
-    auto results = search_cell(target_points);
+    pcms::GridPointSearch2D search_cell(source_mesh, 10, 10); //
+    auto results = search_cell.apply(pcms::CoordinateView(
+		pcms::CoordinateSystem::Cartesian,
+		pcms::MakeConstRank2View(target_points)));
     auto owning_cell_ids = search_cell.GetOwningElementIds(results);
     Omega_h::parallel_for(
       nvertices_target, OMEGA_H_LAMBDA(const Omega_h::LO i) {
@@ -103,8 +105,8 @@ static Omega_h::Write<Omega_h::LO> locate_target_cells(
         source_cell_ids[i] = source_cell_id;
       });
   } else if (dim == 3) {
-    Kokkos::View<Omega_h::Real* [3]> target_points("target_points",
-                                                   nvertices_target);
+    Kokkos::View<pcms::Real**> target_points("target_points",
+                                                   nvertices_target, 3);
     Omega_h::parallel_for(
       nvertices_target, OMEGA_H_LAMBDA(const Omega_h::LO i) {
         target_points(i, 0) = target_coords[i * dim];
@@ -114,7 +116,9 @@ static Omega_h::Write<Omega_h::LO> locate_target_cells(
     Kokkos::fence();
 
     pcms::GridPointSearch3D search_cell(source_mesh, 10, 10, 10);
-    auto results = search_cell(target_points);
+    auto results = search_cell.apply(pcms::CoordinateView(
+		pcms::CoordinateSystem::Cartesian,
+		pcms::MakeConstRank2View(target_points)));
     auto owning_cell_ids = search_cell.GetOwningElementIds(results);
     Omega_h::parallel_for(
       nvertices_target, OMEGA_H_LAMBDA(const Omega_h::LO i) {