Merge pull request #15720 from masterleinad/arborx_1d

Add 1d support for ArborXWrappers

source/arborx/access_traits.cc

@@ -21,22 +21,31 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace ArborXWrappers
 {
+  namespace
+  {
+    template <int dim, typename Number>
+    dealii::Point<3, float>
+    to_3d_float_point(dealii::Point<dim, Number> p)
+    {
+      static_assert(dim >= 1 && dim <= 3);
+      if (dim == 1)
+        return Point<3, float>(p[0], 0.f, 0.f);
+      else if (dim == 2)
+        return Point<3, float>(p[0], p[1], 0.f);
+      else
+        return Point<3, float>(p[0], p[1], p[2]);
+    }
+  } // namespace
+
   // ------------------- PointPredicate ------------------- //
   template <int dim, typename Number>
   PointPredicate::PointPredicate(
     const std::vector<dealii::Point<dim, Number>> &dim_points)
   {
-    static_assert(dim != 1, "dim equal to one is not supported.");
-
     const unsigned int size = dim_points.size();
     points.reserve(size);
     for (unsigned int i = 0; i < size; ++i)
-      {
-        points.emplace_back(static_cast<float>(dim_points[i][0]),
-                            static_cast<float>(dim_points[i][1]),
-                            dim == 2 ? 0.f :
-                                       static_cast<float>(dim_points[i][2]));
-      }
+      points.emplace_back(to_3d_float_point(dim_points[i]));
   }
 
 
@@ -153,20 +162,14 @@ namespace ArborXWrappers
     const std::vector<std::pair<dealii::Point<dim, Number>, Number>>
       &dim_spheres)
   {
-    static_assert(dim != 1, "dim equal to one is not supported.");
-
     const unsigned int size = dim_spheres.size();
     spheres.reserve(size);
     for (unsigned int i = 0; i < size; ++i)
       {
         // ArborX assumes that the center coordinates and the radius use float
         // and the sphere is 3d
-        spheres.emplace_back(std::make_pair(
-          dealii::Point<3, float>(
-            static_cast<float>(dim_spheres[i].first[0]),
-            static_cast<float>(dim_spheres[i].first[1]),
-            dim == 2 ? 0.f : static_cast<float>(dim_spheres[i].first[2])),
-          static_cast<float>(dim_spheres[i].second)));
+        spheres.emplace_back(to_3d_float_point(dim_spheres[i].first),
+                             static_cast<float>(dim_spheres[i].second));
       }
   }
 
@@ -217,6 +220,24 @@ DEAL_II_NAMESPACE_CLOSE
 
 namespace ArborX
 {
+  namespace
+  {
+    template <int dim, typename Number>
+    Point
+    to_arborx_point(dealii::Point<dim, Number> p)
+    {
+      static_assert(dim >= 1 && dim <= 3);
+      if (dim == 1)
+        return {float(p[0]), 0.f, 0.f};
+      else if (dim == 2)
+        return {float(p[0]), float(p[1]), 0.f};
+      else
+        return {float(p[0]), float(p[1]), float(p[2])};
+    }
+  } // namespace
+
+
+
   // ------------------- Point Primitives AccessTraits ------------------- //
   template <int dim, typename Number>
   std::size_t
@@ -236,9 +257,7 @@ namespace ArborX
   {
     // ArborX assumes that the point coordinates use float and that the point
     // is 3d
-    return {static_cast<float>(v[i][0]),
-            static_cast<float>(v[i][1]),
-            dim == 2 ? 0 : static_cast<float>(v[i][2])};
+    return to_arborx_point(v[i]);
   }
 
 
@@ -264,12 +283,7 @@ namespace ArborX
     const dealii::Point<dim, Number> max_corner = boundary_points.second;
     // ArborX assumes that the bounding box coordinates use float and that the
     // bounding box is 3d
-    return {{static_cast<float>(min_corner[0]),
-             static_cast<float>(min_corner[1]),
-             dim == 2 ? 0.f : static_cast<float>(min_corner[2])},
-            {static_cast<float>(max_corner[0]),
-             static_cast<float>(max_corner[1]),
-             dim == 2 ? 0.f : static_cast<float>(max_corner[2])}};
+    return {to_arborx_point(min_corner), to_arborx_point(max_corner)};
   }
 
 
@@ -295,10 +309,7 @@ namespace ArborX
   {
     // ArborX assumes that the center coordinates and the radius use float and
     // the sphere is 3d
-    return {{static_cast<float>(v[i].first[0]),
-             static_cast<float>(v[i].first[1]),
-             dim == 2 ? 0 : static_cast<float>(v[i].first[2])},
-            static_cast<float>(v[i].second)};
+    return {to_arborx_point(v[i].first), static_cast<float>(v[i].second)};
   }
 } // namespace ArborX
 

source/arborx/access_traits.inst.in

@@ -16,7 +16,6 @@
 
 for (DIM : DIMENSIONS; SCALAR : REAL_SCALARS)
   {
-#if DIM > 1
     DEAL_II_NAMESPACE_OPEN
     namespace ArborXWrappers
     {
@@ -55,5 +54,4 @@ for (DIM : DIMENSIONS; SCALAR : REAL_SCALARS)
         std::vector<std::pair<dealii::Point<DIM, SCALAR>, SCALAR>>,
         PrimitivesTag>;
     \}
-#endif
   }

tests/arborx/bounding_box_intersect.cc

@@ -24,6 +24,62 @@
 
 #include "../tests.h"
 
+void
+test_1d()
+{
+  std::vector<BoundingBox<1>> bounding_boxes;
+  std::vector<Point<1>>       points;
+
+  unsigned int n_points_1d = 5;
+  for (unsigned int i = 0; i < n_points_1d; ++i)
+    points.emplace_back(i);
+
+  for (unsigned int i = 0; i < n_points_1d - 1; ++i)
+    bounding_boxes.emplace_back(std::make_pair(points[i], points[i + 1]));
+
+  Point<1> point_a(0.5);
+  Point<1> point_b(1.5);
+  Point<1> point_c(2.2);
+  Point<1> point_d(2.6);
+
+  std::vector<BoundingBox<1>> query_bounding_boxes;
+  query_bounding_boxes.emplace_back(std::make_pair(point_a, point_b));
+  query_bounding_boxes.emplace_back(std::make_pair(point_c, point_d));
+
+  ArborXWrappers::BVH                           bvh(bounding_boxes);
+  ArborXWrappers::BoundingBoxIntersectPredicate bb_intersect(
+    query_bounding_boxes);
+  auto             indices_offset = bvh.query(bb_intersect);
+  std::vector<int> indices        = indices_offset.first;
+  std::vector<int> offset         = indices_offset.second;
+
+  std::vector<int> indices_ref = {0, 1, 2};
+  std::vector<int> offset_ref  = {0, 2, 3};
+
+  AssertDimension(indices.size(), indices_ref.size());
+  AssertDimension(offset.size(), offset_ref.size());
+  for (unsigned int i = 0; i < offset.size() - 1; ++i)
+    {
+      for (unsigned int j = offset[i]; j < offset[i + 1]; ++j)
+        {
+          // The indices associated to each query are not ordered.
+          bool found = false;
+          for (unsigned int k = offset[i]; k < offset[i + 1]; ++k)
+            {
+              if (indices[j] == indices_ref[k])
+                {
+                  found = true;
+                  break;
+                }
+            }
+          AssertThrow(found, ExcInternalError());
+        }
+    }
+  for (unsigned int i = 0; i < offset.size(); ++i)
+    AssertThrow(offset[i] == offset_ref[i], ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
 
 void
 test_2d()
@@ -181,6 +237,7 @@ main(int argc, char **argv)
   initlog();
 
   // tests
+  test_1d();
   test_2d();
   test_3d();
 

tests/arborx/bounding_box_intersect.output

@@ -1,3 +1,4 @@
 
 DEAL::OK
 DEAL::OK
+DEAL::OK

tests/arborx/distributed_tree_intersect.cc

@@ -25,6 +25,70 @@
 #include "../tests.h"
 
 
+void
+test_1d()
+{
+  std::vector<BoundingBox<1>> bounding_boxes;
+  std::vector<Point<1>>       points;
+
+  const unsigned int n_points_1d = 5;
+  const int          rank    = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  const int          n_procs = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+  const int          offset_bb    = 2 * rank * n_points_1d;
+  const int          offset_query = 2 * ((rank + 1) % n_procs) * n_points_1d;
+  for (unsigned int i = 0; i < n_points_1d; ++i)
+    points.emplace_back(i + offset_bb);
+
+  for (unsigned int i = 0; i < n_points_1d - 1; ++i)
+    bounding_boxes.emplace_back(std::make_pair(points[i], points[i + 1]));
+
+  Point<1> point_a(0.5 + offset_query);
+  Point<1> point_b(1.5 + offset_query);
+  Point<1> point_c(2.2 + offset_query);
+  Point<1> point_d(2.6 + offset_query);
+
+  std::vector<BoundingBox<1>> query_bounding_boxes;
+  query_bounding_boxes.emplace_back(std::make_pair(point_a, point_b));
+  query_bounding_boxes.emplace_back(std::make_pair(point_c, point_d));
+
+  ArborXWrappers::DistributedTree               distributed_tree(MPI_COMM_WORLD,
+                                                   bounding_boxes);
+  ArborXWrappers::BoundingBoxIntersectPredicate bb_intersect(
+    query_bounding_boxes);
+  auto indices_ranks_offset = distributed_tree.query(bb_intersect);
+  auto indices_ranks        = indices_ranks_offset.first;
+  auto offsets              = indices_ranks_offset.second;
+
+  std::vector<int> indices_ref = {1, 0, 2};
+  std::vector<int> offsets_ref = {0, 2, 3};
+
+  AssertThrow(indices_ranks.size() == indices_ref.size(), ExcInternalError());
+  AssertThrow(offsets.size() == offsets_ref.size(), ExcInternalError());
+  for (unsigned int i = 0; i < offsets.size() - 1; ++i)
+    {
+      for (int j = offsets[i]; j < offsets[i + 1]; ++j)
+        {
+          // The indices associated to each query are not ordered.
+          bool found = false;
+          for (int k = offsets[i]; k < offsets[i + 1]; ++k)
+            {
+              if ((indices_ranks[j].first == indices_ref[k]) &&
+                  (indices_ranks[j].second == (rank + 1) % n_procs))
+                {
+                  found = true;
+                  break;
+                }
+            }
+          AssertThrow(found, ExcInternalError());
+        }
+    }
+  for (unsigned int i = 0; i < offsets.size(); ++i)
+    AssertThrow(offsets[i] == offsets_ref[i], ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
+
+
 void
 test_2d()
 {
@@ -192,6 +256,7 @@ main(int argc, char **argv)
   initlog();
 
   // tests
+  test_1d();
   test_2d();
   test_3d();
 }

tests/arborx/distributed_tree_intersect.with_arborx_with_mpi=on.mpirun=2.output

@@ -1,3 +1,4 @@
 
-DEAL::OK                                                                                                                                                                                                                                     
+DEAL::OK
+DEAL::OK
 DEAL::OK

tests/arborx/point_intersect.cc

@@ -25,6 +25,61 @@
 #include "../tests.h"
 
 
+void
+test_1d()
+{
+  std::vector<BoundingBox<1>> bounding_boxes;
+  std::vector<Point<1>>       points;
+
+  unsigned int n_points_1d = 5;
+  for (unsigned int i = 0; i < n_points_1d; ++i)
+    points.emplace_back(i);
+
+  for (unsigned int i = 0; i < n_points_1d - 1; ++i)
+    bounding_boxes.emplace_back(std::make_pair(points[i], points[i + 1]));
+
+  std::vector<Point<1>> query_points;
+  query_points.emplace_back(0.5);
+  query_points.emplace_back(1.5);
+  query_points.emplace_back(2.2);
+  query_points.emplace_back(2.6);
+
+
+  ArborXWrappers::BVH                     bvh(bounding_boxes);
+  ArborXWrappers::PointIntersectPredicate pt_intersect(query_points);
+  auto             indices_offset = bvh.query(pt_intersect);
+  std::vector<int> indices        = indices_offset.first;
+  std::vector<int> offset         = indices_offset.second;
+
+  std::vector<int> indices_ref = {0, 1, 2, 2};
+  std::vector<int> offset_ref  = {0, 1, 2, 3, 4};
+
+  AssertThrow(indices.size() == indices_ref.size(), ExcInternalError());
+  AssertThrow(offset.size() == offset_ref.size(), ExcInternalError());
+  for (unsigned int i = 0; i < offset.size() - 1; ++i)
+    {
+      for (unsigned int j = offset[i]; j < offset[i + 1]; ++j)
+        {
+          // The indices associated to each query are not ordered.
+          bool found = false;
+          for (unsigned int k = offset[i]; k < offset[i + 1]; ++k)
+            {
+              if (indices[j] == indices_ref[k])
+                {
+                  found = true;
+                  break;
+                }
+            }
+          AssertThrow(found, ExcInternalError());
+        }
+    }
+  for (unsigned int i = 0; i < offset.size(); ++i)
+    AssertThrow(offset[i] == offset_ref[i], ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
+
+
 void
 test_2d()
 {
@@ -176,6 +231,7 @@ main(int argc, char **argv)
   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
 
   // tests
+  test_1d();
   test_2d();
   test_3d();
 }

tests/arborx/point_intersect.output

@@ -1,3 +1,4 @@
 
 DEAL::OK
 DEAL::OK
+DEAL::OK