]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add wrappers for ArborX 2.0
authorBruno Turcksin <bruno.turcksin@gmail.com>
Tue, 17 Jun 2025 20:47:57 +0000 (16:47 -0400)
committerBruno Turcksin <bruno.turcksin@gmail.com>
Wed, 18 Jun 2025 18:12:14 +0000 (14:12 -0400)
include/deal.II/arborx/access_traits.h
include/deal.II/arborx/bvh.h
include/deal.II/arborx/distributed_tree.h
source/arborx/access_traits.cc
source/arborx/access_traits.inst.in
tests/arborx/bounding_box_intersect.cc
tests/arborx/bounding_box_nearest.cc
tests/arborx/point_intersect.cc
tests/arborx/point_nearest.cc
tests/arborx/sphere_intersect.cc
tests/arborx/sphere_nearest.cc

index 88cb92f09f6cc9b3fbe92db3d25e1bd2de6922e2..9ba7212d3aa3a8ce19a55e1297a9ff1c1e313aa5 100644 (file)
@@ -29,6 +29,7 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace ArborXWrappers
 {
+#  if ARBORX_VERSION_MAJOR < 2
   /**
    * Base class for Point-based predicates providing basic functionality for
    * derived classes, not supposed to be used on its own.
@@ -300,6 +301,384 @@ namespace ArborXWrappers
   private:
     unsigned int n_nearest_neighbors;
   };
+#  else
+  namespace internal
+  {
+    /**
+     * This structure returns the ArborX object associated with the deal.II
+     * object stored in a PairValueIndex.
+     */
+    struct IndexableGetter
+    {
+      template <int dim, typename Number>
+      ArborX::Point<dim, Number>
+      operator()(const ArborX::PairValueIndex<dealii::Point<dim, Number>,
+                                              unsigned int> &pair) const
+      {
+        if constexpr (dim == 1)
+          {
+            return {pair.value[0]};
+          }
+        if constexpr (dim == 2)
+          {
+            return {pair.value[0], pair.value[1]};
+          }
+        if constexpr (dim == 3)
+          {
+            return {pair.value[0], pair.value[1], pair.value[2]};
+          }
+      }
+
+
+
+      template <int dim, typename Number>
+      ArborX::Box<dim, Number>
+      operator()(const ArborX::PairValueIndex<dealii::BoundingBox<dim, Number>,
+                                              unsigned int> &pair) const
+      {
+        const auto boundary_points = pair.value.get_boundary_points();
+        const dealii::Point<dim, Number> min_corner = boundary_points.first;
+        const dealii::Point<dim, Number> max_corner = boundary_points.second;
+        if constexpr (dim == 1)
+          {
+            return {{min_corner[0]}, {max_corner[0]}};
+          }
+        if constexpr (dim == 2)
+          {
+            return {{min_corner[0], min_corner[1]},
+                    {max_corner[0], max_corner[1]}};
+          }
+        if constexpr (dim == 3)
+          {
+            return {{min_corner[0], min_corner[1], min_corner[2]},
+                    {max_corner[0], max_corner[1], max_corner[2]}};
+          }
+      }
+
+
+
+      template <int dim, typename Number>
+      ArborX::Sphere<dim, Number>
+      operator()(const ArborX::PairValueIndex<
+                 std::pair<dealii::Point<dim, Number>, Number>,
+                 unsigned int> &pair) const
+      {
+        if constexpr (dim == 1)
+          {
+            return {{pair.value.first[0]}, pair.value.second};
+          }
+        if constexpr (dim == 2)
+          {
+            return {{pair.value.first[0], pair.value.first[1]},
+                    pair.value.second};
+          }
+        if constexpr (dim == 3)
+          {
+            return {{pair.value.first[0],
+                     pair.value.first[1],
+                     pair.value.first[2]},
+                    pair.value.second};
+          }
+      }
+    };
+
+
+
+    /**
+     * Callback to extract the index of each primitive that satisfies a query.
+     */
+    struct ExtractIndex
+    {
+      template <typename Query, typename Value, typename Output>
+      KOKKOS_FUNCTION void
+      operator()(const Query &, const Value &value, const Output &out) const
+      {
+        out(value.index);
+      }
+    };
+
+
+
+    /**
+     * Callback to extract the index and the rank of each primitive that
+     * satisfies a query.
+     */
+    struct ExtractIndexRank
+    {
+      unsigned int rank;
+
+      template <typename Predicate, typename Value, typename Output>
+      KOKKOS_FUNCTION void
+      operator()(const Predicate &,
+                 const ArborX::PairValueIndex<Value> &value,
+                 const Output                        &out) const
+      {
+        out({value.index, rank});
+      }
+    };
+  } // namespace internal
+
+
+
+  /**
+   * This class defines a predicate used by ArborXWrappers::BVH to determine
+   * for given points which of the primitives used to build the
+   * ArborXWrappers::BVH intersect with them.
+   */
+  template <int dim, typename Number>
+  class PointIntersectPredicate
+  {
+  public:
+    /**
+     * Constructor. @p points is a list of points which we are interested in
+     * knowing if they intersect ArborXWrappers::BVH primitives.
+     */
+    PointIntersectPredicate(
+      const std::vector<dealii::Point<dim, Number>> &points);
+
+    /**
+     * The number of points stored in the structure.
+     */
+    std::size_t
+    size() const;
+
+    /**
+     * Return the `i`th Point stored in the object.
+     */
+    const dealii::Point<dim, Number> &
+    get(unsigned int i) const;
+
+    /**
+     * A flag that specifies if the predicate is nearest neighbors search.
+     */
+    static constexpr bool is_nearest = false;
+
+  private:
+    std::vector<dealii::Point<dim, Number>> points;
+  };
+
+
+
+  /**
+   * This class defines a predicate used by ArborXWrappers::BVH to determine
+   * for given points which are the nearest primitives among the ones
+   * used to build the ArborXWrappers::BVH.
+   */
+  template <int dim, typename Number>
+  class PointNearestPredicate
+  {
+  public:
+    /**
+     * Constructor. @p points is a list of points for which we are interested in
+     * the @p n_nearest_neighbors in the ArborXWrappers::BVH primitives.
+     */
+    PointNearestPredicate(const std::vector<dealii::Point<dim, Number>> &points,
+                          const unsigned int n_nearest_neighbors);
+
+    /**
+     * Return the number of nearest neighbors we are looking for.
+     */
+    unsigned int
+    get_n_nearest_neighbors() const;
+
+    /**
+     * The number of points stored in the structure.
+     */
+    std::size_t
+    size() const;
+
+    /**
+     * Return the `i`th Point stored in the object.
+     */
+    const dealii::Point<dim, Number> &
+    get(unsigned int i) const;
+
+
+    /**
+     * A flag that specifies if the predicate is nearest neighbors search.
+     */
+    static constexpr bool is_nearest = true;
+
+  private:
+    std::vector<dealii::Point<dim, Number>> points;
+    unsigned int                            n_nearest_neighbors;
+  };
+
+
+
+  /**
+   * This class is used by ArborXWrappers::BVH to determine for given bounding
+   * boxes which of the primitives used to build the ArborXWrappers::BVH
+   * intersect with them.
+   */
+  template <int dim, typename Number>
+  class BoundingBoxIntersectPredicate
+  {
+  public:
+    /**
+     * Constructor. @p bounding_boxes is a list of bounding boxes which we are interested in
+     * knowing if they intersect ArborXWrappers::BVH primitives.
+     */
+    BoundingBoxIntersectPredicate(
+      const std::vector<dealii::BoundingBox<dim, Number>> &bounding_boxes);
+
+    /**
+     * The number of points stored in the structure.
+     */
+    std::size_t
+    size() const;
+
+    /**
+     * Return the `i`th Point stored in the object.
+     */
+    const dealii::BoundingBox<dim, Number> &
+    get(unsigned int i) const;
+
+    /**
+     * A flag that specifies if the predicate is nearest neighbors search.
+     */
+    static constexpr bool is_nearest = false;
+
+  private:
+    std::vector<dealii::BoundingBox<dim, Number>> bounding_boxes;
+  };
+
+
+
+  /**
+   * This class is used by ArborXWrappers::BVH to determine for given bounding
+   * boxes which are the nearest primitives among the ones used to
+   * build the ArborXWrappers::BVH.
+   */
+  template <int dim, typename Number>
+  class BoundingBoxNearestPredicate
+  {
+  public:
+    /**
+     * Constructor. @p bounding_boxes is a list of bounding boxes for which are interested in
+     * knowing the @p n_nearest_neighbors nearest primitives used to build the
+     * ArborXWrappers::BVH.
+     */
+    BoundingBoxNearestPredicate(
+      const std::vector<dealii::BoundingBox<dim, Number>> &bounding_boxes,
+      const unsigned int                                   n_nearest_neighbors);
+
+    /**
+     * Return the number of nearest neighbors we are looking for.
+     */
+    unsigned int
+    get_n_nearest_neighbors() const;
+
+    /**
+     * The number of points stored in the structure.
+     */
+    std::size_t
+    size() const;
+
+    /**
+     * Return the `i`th Point stored in the object.
+     */
+    const dealii::BoundingBox<dim, Number> &
+    get(unsigned int i) const;
+
+    /**
+     * A flag that specifies if the predicate is nearest neighbors search.
+     */
+    static constexpr bool is_nearest = true;
+
+  private:
+    std::vector<dealii::BoundingBox<dim, Number>> bounding_boxes;
+    unsigned int                                  n_nearest_neighbors;
+  };
+
+
+
+  /**
+   * This class defines a predicate used by ArborXWrappers::BVH to determine
+   * for given spheres which of the primitives used to build the
+   * ArborXWrappers::BVH intersect with them.
+   */
+  template <int dim, typename Number>
+  class SphereIntersectPredicate
+  {
+  public:
+    /**
+     * Constructor. @p spheres is a list of spheres which we are interested in
+     * knowing if they intersect ArborXWrappers::BVH primitives.
+     */
+    SphereIntersectPredicate(
+      const std::vector<std::pair<dealii::Point<dim, Number>, Number>>
+        &spheres);
+
+    /**
+     * The number of points stored in the structure.
+     */
+    std::size_t
+    size() const;
+
+    /**
+     * Return the `i`th Point stored in the object.
+     */
+    const std::pair<dealii::Point<dim, Number>, Number> &
+    get(unsigned int) const;
+
+    /**
+     * A flag that specifies if the predicate is nearest neighbors search.
+     */
+    static constexpr bool is_nearest = false;
+
+  private:
+    std::vector<std::pair<dealii::Point<dim, Number>, Number>> spheres;
+  };
+
+
+
+  /**
+   * This class defines a predicate used by ArborXWrappers::BVH to determine,
+   * for the given spheres, which are the nearest bounding primitives among
+   * the ones used to build the ArborXWrappers::BVH.
+   */
+  template <int dim, typename Number>
+  class SphereNearestPredicate
+  {
+  public:
+    /**
+     * Constructor. @p spheres is a list of spheres for which we are
+     * interested in the @p n_nearest_neighbors in the ArborXWrappers::BVH
+     * primitives.
+     */
+    SphereNearestPredicate(
+      const std::vector<std::pair<dealii::Point<dim, Number>, Number>> &spheres,
+      const unsigned int n_nearest_neighbors);
+
+    /**
+     * Return the number of nearest neighbors we are looking for.
+     */
+    unsigned int
+    get_n_nearest_neighbors() const;
+
+    /**
+     * The number of points stored in the structure.
+     */
+    std::size_t
+    size() const;
+
+    /**
+     * Return the `i`th Point stored in the object.
+     */
+    const std::pair<dealii::Point<dim, Number>, Number> &
+    get(unsigned int) const;
+
+    /**
+     * A flag that specifies if the predicate is nearest neighbors search.
+     */
+    static constexpr bool is_nearest = true;
+
+  private:
+    std::vector<std::pair<dealii::Point<dim, Number>, Number>> spheres;
+    unsigned int n_nearest_neighbors;
+  };
+#  endif
 } // namespace ArborXWrappers
 
 DEAL_II_NAMESPACE_CLOSE
@@ -309,6 +688,7 @@ DEAL_II_NAMESPACE_CLOSE
  */
 namespace ArborX
 {
+#  if ARBORX_VERSION_MAJOR < 2
   /**
    * This struct allows ArborX to use std::vector<dealii::Point> as
    * primitive.
@@ -383,9 +763,32 @@ namespace ArborX
     get(const std::vector<std::pair<dealii::Point<dim, Number>, Number>> &v,
         std::size_t                                                       i);
   };
+#  else
+  /**
+   * This struct allows ArborX to use std::vector<T> as primitive.
+   */
+  template <typename T>
+  struct AccessTraits<std::vector<T>>
+  {
+    using memory_space = Kokkos::HostSpace;
+
+    /**
+     * Return the size of the vector @p v.
+     */
+    static std::size_t
+    size(const std::vector<T> &v);
+
+    /**
+     * Return the ith element of the vector @p v.
+     */
+    static T
+    get(const std::vector<T> &v, std::size_t i);
+  };
+#  endif
 
 
 
+#  if ARBORX_VERSION_MAJOR < 2
   /**
    * This struct allows ArborX to use PointIntersectPredicate as a predicate.
    */
@@ -485,11 +888,9 @@ namespace ArborX
     size(const dealii::ArborXWrappers::BoundingBoxNearestPredicate &bb_nearest);
 
     /**
-     * Return an
-     * Arbox::nearest(ArborX::Box,
+     * Return an Arbox::nearest(ArborX::Box,
      * BoundingBoxtNearestPredicate::get_n_nearest_neighbors) object constructed
-     * from the
-     * `i`th dealii::BoundingBox stored in @p bb_nearest.
+     * from the `i`th dealii::BoundingBox stored in @p bb_nearest.
      */
     static auto
     get(const dealii::ArborXWrappers::BoundingBoxNearestPredicate &bb_nearest,
@@ -499,7 +900,8 @@ namespace ArborX
 
 
   /**
-   * This struct allows ArborX to use SphereIntersectPredicate as a predicate.
+   * This struct allows ArborX to use SphereIntersectPredicate as a
+   predicate.
    */
   template <>
   struct AccessTraits<dealii::ArborXWrappers::SphereIntersectPredicate,
@@ -548,6 +950,180 @@ namespace ArborX
     get(const dealii::ArborXWrappers::SphereNearestPredicate &sph_nearest,
         std::size_t                                           i);
   };
+#  else
+  /**
+   * This struct allows ArborX to use PointIntersectPredicate as a predicate.
+   */
+  template <int dim, typename Number>
+  struct AccessTraits<
+    dealii::ArborXWrappers::PointIntersectPredicate<dim, Number>>
+  {
+    using memory_space = Kokkos::HostSpace;
+
+    /**
+     * The number of points stored in @p pt_intersect.
+     */
+    static std::size_t
+    size(const dealii::ArborXWrappers::PointIntersectPredicate<dim, Number>
+           &pt_intersect);
+
+    /**
+     * Return an ArborX::intersects(ArborX::Point) object constructed from the
+     * `i`th dealii::Point stored in @p pt_intersect.
+     */
+    static auto
+    get(const dealii::ArborXWrappers::PointIntersectPredicate<dim, Number>
+                   &pt_intersect,
+        std::size_t i);
+  };
+
+
+
+  /**
+   * This struct allows ArborX to use PointNearestPredicate as a predicate.
+   */
+  template <int dim, typename Number>
+  struct AccessTraits<
+    dealii::ArborXWrappers::PointNearestPredicate<dim, Number>>
+  {
+    using memory_space = Kokkos::HostSpace;
+
+    /**
+     * The number of points stored in @p pt_nearest.
+     */
+    static std::size_t
+    size(const dealii::ArborXWrappers::PointNearestPredicate<dim, Number>
+           &pt_nearest);
+
+    /**
+     * Return an ArborX::nearest(ArborX::Point,
+     * PointNearestPredicate::get_n_nearest_neighbors) object constructed from
+     * the `i`th dealii::Point stored in @p pt_nearest.
+     */
+    static auto
+    get(const dealii::ArborXWrappers::PointNearestPredicate<dim, Number>
+                   &pt_nearest,
+        std::size_t i);
+  };
+
+
+
+  /**
+   * This struct allows ArborX to use BoundingBoxIntersectPredicate as a
+   * predicate.
+   */
+  template <int dim, typename Number>
+  struct AccessTraits<
+    dealii::ArborXWrappers::BoundingBoxIntersectPredicate<dim, Number>>
+  {
+    using memory_space = Kokkos::HostSpace;
+
+    /**
+     * The number of bounding boxes stored in @p bb_intersect.
+     */
+    static std::size_t
+    size(
+      const dealii::ArborXWrappers::BoundingBoxIntersectPredicate<dim, Number>
+        &bb_intersect);
+
+    /**
+     * Return an Arbox::intersects(ArborX::Box) object constructed from the
+     * `i`th dealii::BoundingBox stored in @p bb_intersect.
+     */
+    static auto
+    get(const dealii::ArborXWrappers::BoundingBoxIntersectPredicate<dim, Number>
+                   &bb_intersect,
+        std::size_t i);
+  };
+
+
+
+  /**
+   * This struct allows ArborX to use BoundingBoxNearstPredicate as a
+   * predicate.
+   */
+  template <int dim, typename Number>
+  struct AccessTraits<
+    dealii::ArborXWrappers::BoundingBoxNearestPredicate<dim, Number>>
+  {
+    using memory_space = Kokkos::HostSpace;
+
+    /**
+     * The number of bounding boxes stored in @p bb_nearest.
+     */
+    static std::size_t
+    size(const dealii::ArborXWrappers::BoundingBoxNearestPredicate<dim, Number>
+           &bb_nearest);
+
+    /**
+     * Return an Arbox::nearest(ArborX::Box,
+     * BoundingBoxtNearestPredicate::get_n_nearest_neighbors) object constructed
+     * from the `i`th dealii::BoundingBox stored in @p bb_nearest.
+     */
+    static auto
+    get(const dealii::ArborXWrappers::BoundingBoxNearestPredicate<dim, Number>
+                   &bb_nearest,
+        std::size_t i);
+  };
+
+
+
+  /**
+   * This struct allows ArborX to use SphereIntersectPredicate as a
+   predicate.
+   */
+  template <int dim, typename Number>
+  struct AccessTraits<
+    dealii::ArborXWrappers::SphereIntersectPredicate<dim, Number>>
+  {
+    using memory_space = Kokkos::HostSpace;
+
+    /**
+     * The number of points stored in @p sph_intersect.
+     */
+    static std::size_t
+    size(const dealii::ArborXWrappers::SphereIntersectPredicate<dim, Number>
+           &sph_intersect);
+
+    /**
+     * Return an ArborX::intersects(ArborX::Sphere) object constructed from the
+     * `i`th sphere stored in @p sph_intersect.
+     */
+    static auto
+    get(const dealii::ArborXWrappers::SphereIntersectPredicate<dim, Number>
+                   &sph_intersect,
+        std::size_t i);
+  };
+
+
+
+  /**
+   * This struct allows ArborX to use SphereNearestPredicate as a predicate.
+   */
+  template <int dim, typename Number>
+  struct AccessTraits<
+    dealii::ArborXWrappers::SphereNearestPredicate<dim, Number>>
+  {
+    using memory_space = Kokkos::HostSpace;
+
+    /**
+     * The number of spheres stored in @p sph_nearest.
+     */
+    static std::size_t
+    size(const dealii::ArborXWrappers::SphereNearestPredicate<dim, Number>
+           &sph_nearest);
+
+    /**
+     * Return an ArborX::nearest(ArborX::Sphere,
+     * SphereNearestPredicate::get_n_nearest_neighbors) object constructed from
+     * the `i`th sphere stored in @p sph_nearest.
+     */
+    static auto
+    get(const dealii::ArborXWrappers::SphereNearestPredicate<dim, Number>
+                   &sph_nearest,
+        std::size_t i);
+  };
+#  endif
 
   // ------------------------------- Inline ----------------------------------//
 
@@ -555,6 +1131,7 @@ namespace ArborX
   // header file otherwise the return type of auto get() cannot be determined.
   // We use auto because ArborX does not expose the type of intersects
 
+#  if ARBORX_VERSION_MAJOR < 2
   inline std::size_t
   AccessTraits<dealii::ArborXWrappers::PointIntersectPredicate, PredicatesTag>::
     size(const dealii::ArborXWrappers::PointIntersectPredicate &pt_intersect)
@@ -695,6 +1272,247 @@ namespace ArborX
                    sph_nearest.get_n_nearest_neighbors());
   }
 } // namespace ArborX
+#  else
+  template <int dim, typename Number>
+  inline std::size_t
+  AccessTraits<dealii::ArborXWrappers::PointIntersectPredicate<dim, Number>>::
+    size(const dealii::ArborXWrappers::PointIntersectPredicate<dim, Number>
+           &pt_intersect)
+  {
+    return pt_intersect.size();
+  }
+
+
+
+  template <int dim, typename Number>
+  inline auto
+  AccessTraits<dealii::ArborXWrappers::PointIntersectPredicate<dim, Number>>::
+    get(const dealii::ArborXWrappers::PointIntersectPredicate<dim, Number>
+                   &pt_intersect,
+        std::size_t i)
+  {
+    const auto dealii_point = pt_intersect.get(i);
+    if constexpr (dim == 1)
+      {
+        return intersects(Point{dealii_point[0]});
+      }
+    if constexpr (dim == 2)
+      {
+        return intersects(Point{dealii_point[0], dealii_point[1]});
+      }
+    if constexpr (dim == 3)
+      {
+        return intersects(
+          Point{dealii_point[0], dealii_point[1], dealii_point[2]});
+      }
+  }
+
+
+
+  template <int dim, typename Number>
+  inline std::size_t
+  AccessTraits<dealii::ArborXWrappers::PointNearestPredicate<dim, Number>>::
+    size(const dealii::ArborXWrappers::PointNearestPredicate<dim, Number>
+           &pt_nearest)
+  {
+    return pt_nearest.size();
+  }
+
+
+
+  template <int dim, typename Number>
+  inline auto
+  AccessTraits<dealii::ArborXWrappers::PointNearestPredicate<dim, Number>>::get(
+    const dealii::ArborXWrappers::PointNearestPredicate<dim, Number>
+               &pt_nearest,
+    std::size_t i)
+  {
+    const auto dealii_point = pt_nearest.get(i);
+    if constexpr (dim == 1)
+      {
+        return nearest(Point{dealii_point[0]},
+                       pt_nearest.get_n_nearest_neighbors());
+      }
+    if constexpr (dim == 2)
+      {
+        return nearest(Point{dealii_point[0], dealii_point[1]},
+                       pt_nearest.get_n_nearest_neighbors());
+      }
+    if constexpr (dim == 3)
+      {
+        return nearest(Point{dealii_point[0], dealii_point[1], dealii_point[2]},
+                       pt_nearest.get_n_nearest_neighbors());
+      }
+  }
+
+
+
+  template <int dim, typename Number>
+  inline std::size_t
+  AccessTraits<
+    dealii::ArborXWrappers::BoundingBoxIntersectPredicate<dim, Number>>::
+    size(
+      const dealii::ArborXWrappers::BoundingBoxIntersectPredicate<dim, Number>
+        &bb_intersect)
+  {
+    return bb_intersect.size();
+  }
+
+
+
+  template <int dim, typename Number>
+  inline auto
+  AccessTraits<
+    dealii::ArborXWrappers::BoundingBoxIntersectPredicate<dim, Number>>::
+    get(const dealii::ArborXWrappers::BoundingBoxIntersectPredicate<dim, Number>
+                   &bb_intersect,
+        std::size_t i)
+  {
+    const auto boundary_points = bb_intersect.get(i).get_boundary_points();
+    const dealii::Point<dim, Number> min_corner = boundary_points.first;
+    const dealii::Point<dim, Number> max_corner = boundary_points.second;
+
+    if constexpr (dim == 1)
+      {
+        return intersects(Box{{min_corner[0]}, {max_corner[0]}});
+      }
+    if constexpr (dim == 2)
+      {
+        return intersects(
+          Box{{min_corner[0], min_corner[1]}, {max_corner[0], max_corner[1]}});
+      }
+    if constexpr (dim == 3)
+      {
+        return intersects(Box{{min_corner[0], min_corner[1], min_corner[2]},
+                              {max_corner[0], max_corner[1], max_corner[2]}});
+      }
+  }
+
+
+
+  template <int dim, typename Number>
+  inline std::size_t
+  AccessTraits<
+    dealii::ArborXWrappers::BoundingBoxNearestPredicate<dim, Number>>::
+    size(const dealii::ArborXWrappers::BoundingBoxNearestPredicate<dim, Number>
+           &bb_nearest)
+  {
+    return bb_nearest.size();
+  }
+
+
+
+  template <int dim, typename Number>
+  inline auto
+  AccessTraits<
+    dealii::ArborXWrappers::BoundingBoxNearestPredicate<dim, Number>>::
+    get(const dealii::ArborXWrappers::BoundingBoxNearestPredicate<dim, Number>
+                   &bb_nearest,
+        std::size_t i)
+  {
+    const auto boundary_points = bb_nearest.get(i).get_boundary_points();
+    const dealii::Point<dim, Number> min_corner = boundary_points.first;
+    const dealii::Point<dim, Number> max_corner = boundary_points.second;
+
+    if constexpr (dim == 1)
+      {
+        return nearest(Box{{min_corner[0]}, {max_corner[0]}},
+                       bb_nearest.get_n_nearest_neighbors());
+      }
+    if constexpr (dim == 2)
+      {
+        return nearest(Box{{min_corner[0], min_corner[1]},
+                           {max_corner[0], max_corner[1]}},
+                       bb_nearest.get_n_nearest_neighbors());
+      }
+    if constexpr (dim == 3)
+      {
+        return nearest(Box{{min_corner[0], min_corner[1], min_corner[2]},
+                           {max_corner[0], max_corner[1], max_corner[2]}},
+                       bb_nearest.get_n_nearest_neighbors());
+      }
+  }
+
+
+
+  template <int dim, typename Number>
+  inline std::size_t
+  AccessTraits<dealii::ArborXWrappers::SphereIntersectPredicate<dim, Number>>::
+    size(const dealii::ArborXWrappers::SphereIntersectPredicate<dim, Number>
+           &sph_intersect)
+  {
+    return sph_intersect.size();
+  }
+
+
+
+  template <int dim, typename Number>
+  inline auto
+  AccessTraits<dealii::ArborXWrappers::SphereIntersectPredicate<dim, Number>>::
+    get(const dealii::ArborXWrappers::SphereIntersectPredicate<dim, Number>
+                   &sph_intersect,
+        std::size_t i)
+  {
+    const auto sphere = sph_intersect.get(i);
+    if constexpr (dim == 1)
+      {
+        return intersects(Sphere{{sphere.first[0]}, sphere.second});
+      }
+    if constexpr (dim == 2)
+      {
+        return intersects(
+          Sphere{{sphere.first[0], sphere.first[1]}, sphere.second});
+      }
+    if constexpr (dim == 3)
+      {
+        return intersects(
+          Sphere{{sphere.first[0], sphere.first[1], sphere.first[2]},
+                 sphere.second});
+      }
+  }
+
+
+
+  template <int dim, typename Number>
+  inline std::size_t
+  AccessTraits<dealii::ArborXWrappers::SphereNearestPredicate<dim, Number>>::
+    size(const dealii::ArborXWrappers::SphereNearestPredicate<dim, Number>
+           &sph_nearest)
+  {
+    return sph_nearest.size();
+  }
+
+
+
+  template <int dim, typename Number>
+  inline auto
+  AccessTraits<dealii::ArborXWrappers::SphereNearestPredicate<dim, Number>>::
+    get(const dealii::ArborXWrappers::SphereNearestPredicate<dim, Number>
+                   &sph_nearest,
+        std::size_t i)
+  {
+    const auto sphere = sph_nearest.get(i);
+    if constexpr (dim == 1)
+      {
+        return nearest(Sphere{{sphere.first[0]}, sphere.second},
+                       sph_nearest.get_n_nearest_neighbors());
+      }
+    if constexpr (dim == 2)
+      {
+        return nearest(Sphere{{sphere.first[0], sphere.first[1]},
+                              sphere.second},
+                       sph_nearest.get_n_nearest_neighbors());
+      }
+    if constexpr (dim == 3)
+      {
+        return nearest(
+          Sphere{{sphere.first[0], sphere.first[1], sphere.first[2]},
+                 sphere.second},
+          sph_nearest.get_n_nearest_neighbors());
+      }
+  }
+#  endif
+} // namespace ArborX
 
 #else
 
index 6819cdd6ecad76a0e1e7868a70b57fa03dc099da..ff639ee7e766725ccbf2f1cf405b888b57baceb6 100644 (file)
@@ -20,6 +20,7 @@
 #ifdef DEAL_II_WITH_ARBORX
 #  include <deal.II/arborx/access_traits.h>
 
+#  include <ArborX_Config.hpp>
 #  include <ArborX_LinearBVH.hpp>
 #  include <Kokkos_Core.hpp>
 
@@ -46,10 +47,8 @@ namespace ArborXWrappers
    * several operations on sets of geometric objects efficiently, such as in
    * collision detection and ray tracing.
    * </blockquote>
-   *
-   * Because ArborX uses Kokkos, Kokkos needs to be initialized and finalized
-   * before using this class.
    */
+#  if ARBORX_VERSION_MAJOR < 2
   class BVH
   {
   public:
@@ -181,6 +180,135 @@ namespace ArborXWrappers
 
     return {indices_vector, offset_vector};
   }
+#  else
+  template <typename Value>
+  class BVH
+  {
+  public:
+    /**
+     * Constructor. Use a vector of Point, BoundingBox, or sphere (std::pair<Point,Number>) @p values as primitives.
+     */
+    BVH(const std::vector<Value> &values);
+
+    /**
+     * Return the indices of those primitive objects that satisfy the @p queries.
+     * Because @p queries can contain multiple queries, the function returns a pair
+     * of indices and offsets.
+     *
+     * Below is an example piece of code that does the following: Let us assume
+     * that we have a set of bounding boxes for objects -- say, the bounding
+     * boxes of each of the cells in a triangulation, or the bounding boxes for
+     * each of the parts of a triangulation in a parallel triangulation. We will
+     * store those in the `bvh_bounding_boxes` array below.
+     *
+     * Let us then also assume that we have a set of other bounding boxes, let's
+     * say for small objects that are moving around in our domain. We will put
+     * these bounding boxes into the `bb_intersect` array. The question we would
+     * then like to answer is the following: with which of the BVH bounding
+     * box(es) do each of the bb_intersect bounding boxes intersect? In other
+     * words, in which cell(s) or partition(s) are the particles?
+     *
+     * This query is answered by the following piece of code:
+     *
+     * @code
+     * const std::vector<BoundingBox<dim>> query_bounding_boxes = ...
+     * ArborXWrappers::BoundingBoxIntersectPredicate
+     * bb_intersect(query_bounding_boxes);
+     *
+     * const std::vector<BoundingBox<dim>> bvh_bounding_boxes = ...
+     * ArborxWrappers::BVH<BoundingBox<dim>> bvh(bvh_bounding_boxes);
+     *
+     * auto [indices, offset] = bvh.query(bb_intersect);
+     * @endcode
+     *
+     * The elements of `bvh_bounding_boxes` that intersect the `j`th BoundingBox
+     * of `query_bounding_boxes` are given by:
+     *
+     * @code
+     * std::vector<int> bvh_bounding_box_indices;
+     * for (int i = offset[j]; i < offset[j+1]; ++i)
+     *   bvh_bounding_box_indices.push_back(indices[i]);
+     * @endcode
+     *
+     * In many other applications, we are interested not only in finding which
+     * bounding boxes another bounding box lies in, but in fact which bounding
+     * boxes individual points lie in -- say, if instead of objects we have
+     * point-like particles moving around. In that case, we would need to answer
+     * a query for points, and this can be done as follows:
+     *
+     * @code
+     * const std::vector<Point<dim>> query_points = ...
+     * ArborXWrappers::PointIntersectPredicate pt_intersect(query_points);
+     *
+     * const std::vector<BoundingBox<dim>> bvh_bounding_boxes = ...
+     * ArborxWrappers::BVH<BoundingBox<dim>> bvh(bvh_bounding_boxes);
+     *
+     * auto [indices, offset] = bvh.query(pt_intersect);
+     * @endcode
+     *
+     * As a final example, we want to show how to find the five nearest points
+     * of a given set of points. This can done as follows:
+     *
+     * @code
+     * const std::vector<Point<dim>> query_points = ...
+     * ArborXWrappers::PointNearestPredicate pt_nearest(query_points, 5);
+     *
+     * const std::vector<Point<dim>> bvh_points = ...
+     * ArborxWrappers::BVH<Point<dim>> bvh(bvh_points);
+     *
+     * auto [indices, offset] = bvh.query(pt_nearest);
+     * @endcode
+     */
+    template <typename QueryType>
+    std::pair<std::vector<int>, std::vector<int>>
+    query(const QueryType &queries);
+
+  private:
+    /**
+     * Underlying ArborX object.
+     */
+    ArborX::BVH<Kokkos::HostSpace,
+                ArborX::PairValueIndex<Value, unsigned int>,
+                internal::IndexableGetter>
+      bvh;
+  };
+
+
+
+  template <typename Value>
+  BVH<Value>::BVH(const std::vector<Value> &values)
+    : bvh(Kokkos::DefaultHostExecutionSpace{},
+          ArborX::Experimental::attach_indices(values),
+          internal::IndexableGetter{})
+  {}
+
+
+
+  template <typename Value>
+  template <typename QueryType>
+  std::pair<std::vector<int>, std::vector<int>>
+  BVH<Value>::query(const QueryType &queries)
+  {
+    Kokkos::View<int *, Kokkos::HostSpace> indices("indices", 0);
+
+    Kokkos::View<int *, Kokkos::HostSpace> offset("offset", 0);
+    bvh.query(Kokkos::DefaultHostExecutionSpace{},
+              queries,
+              internal::ExtractIndex{},
+              indices,
+              offset);
+    std::vector<int> indices_vector;
+    indices_vector.insert(indices_vector.begin(),
+                          indices.data(),
+                          indices.data() + indices.extent(0));
+    std::vector<int> offset_vector;
+    offset_vector.insert(offset_vector.begin(),
+                         offset.data(),
+                         offset.data() + offset.extent(0));
+
+    return {indices_vector, offset_vector};
+  }
+#  endif
 } // namespace ArborXWrappers
 
 DEAL_II_NAMESPACE_CLOSE
index c02bd9d2feff821c33164e98080d60c169e5f83c..728a0bf6b8e2ed593393dd8aaac19a35c4d5446d 100644 (file)
@@ -20,6 +20,8 @@
 #if defined(DEAL_II_ARBORX_WITH_MPI) && defined(DEAL_II_WITH_MPI)
 #  include <deal.II/arborx/access_traits.h>
 
+#  include <deal.II/base/mpi.h>
+
 #  include <ArborX_DistributedTree.hpp>
 #  include <Kokkos_Core.hpp>
 
@@ -30,10 +32,8 @@ namespace ArborXWrappers
   /**
    * This class implements a wrapper around ArborX::DistributedTree, the
    * distributed version of ArborX::BVH.
-   *
-   * Because ArborX uses Kokkos, Kokkos needs to be initialized and finalized
-   * before using this class.
    */
+#  if ARBORX_VERSION_MAJOR < 2
   class DistributedTree
   {
   public:
@@ -63,8 +63,10 @@ namespace ArborXWrappers
      * Valid `QueryType` classes are
      * ArborXWrappers::BoundingBoxIntersectPredicate,
      * ArborXWrappers::BoundingBoxNearestPredicate,
-     * ArborXWrappers::PointIntersectPredicate, and
+     * ArborXWrappers::PointIntersectPredicate,
      * ArborXWrappers::PointNearestPredicate,
+     * ArborXWrappers::SphereIntersectPredicate, and
+     * ArborXWrappers::SphereNearestPredicate.
      */
     template <typename QueryType>
     std::pair<std::vector<std::pair<int, int>>, std::vector<int>>
@@ -125,6 +127,105 @@ namespace ArborXWrappers
 
     return {indices_ranks_vector, offsets_vector};
   }
+#  else
+  template <typename Value>
+  class DistributedTree
+  {
+  public:
+    /**
+     * Constructor. Use a vector of Point, BoundingBox, or sphere
+     * (std::pair<Point,Number>) @p values as primitives. The objects in @p values
+     * are local to the MPI process.
+     */
+    DistributedTree(const MPI_Comm comm, const std::vector<Value> &values);
+
+    /**
+     * Return the indices and the MPI ranks of those BoundingBox objects that
+     * satisfy the @p queries.
+     * Because @p queries can contain multiple queries, the function returns a pair
+     * of indices and ranks and the associated offsets.
+     *
+     * Valid `QueryType` classes are
+     * ArborXWrappers::BoundingBoxIntersectPredicate,
+     * ArborXWrappers::BoundingBoxNearestPredicate,
+     * ArborXWrappers::PointIntersectPredicate,
+     * ArborXWrappers::PointNearestPredicate,
+     * ArborXWrappers::SphereIntersectPredicate, and
+     * ArborXWrappers::SphereNearestPredicate.
+     */
+    template <typename QueryType>
+    std::pair<std::vector<std::pair<int, int>>, std::vector<int>>
+    query(const QueryType &queries);
+
+  private:
+    /**
+     * Underlying ArborX object.
+     */
+    ArborX::DistributedTree<Kokkos::HostSpace,
+                            ArborX::PairValueIndex<Value, unsigned int>,
+                            internal::IndexableGetter>
+      distributed_tree;
+    /**
+     * Callback to extract the index and the rank of each primitive that
+     * satisfies a query.
+     */
+    internal::ExtractIndexRank callback;
+  };
+
+
+
+  template <typename Value>
+  DistributedTree<Value>::DistributedTree(const MPI_Comm            comm,
+                                          const std::vector<Value> &values)
+    : distributed_tree(comm,
+                       Kokkos::DefaultHostExecutionSpace{},
+                       ArborX::Experimental::attach_indices(values),
+                       internal::IndexableGetter{})
+    , callback{Utilities::MPI::this_mpi_process(comm)}
+  {}
+
+  template <typename Value>
+  template <typename QueryType>
+  std::pair<std::vector<std::pair<int, int>>, std::vector<int>>
+  DistributedTree<Value>::query(const QueryType &queries)
+  {
+    Kokkos::View<int *, Kokkos::HostSpace> offsets("offsets", 0);
+    Kokkos::View<Kokkos::pair<int, int> *, Kokkos::HostSpace> indices_ranks(
+      "indices_ranks", 0);
+    if constexpr (QueryType::is_nearest)
+      {
+        distributed_tree.query(
+          Kokkos::DefaultHostExecutionSpace{},
+          queries,
+          ArborX::Experimental::declare_callback_constrained(callback),
+          indices_ranks,
+          offsets);
+      }
+    else
+      {
+        distributed_tree.query(Kokkos::DefaultHostExecutionSpace{},
+                               queries,
+                               callback,
+                               indices_ranks,
+                               offsets);
+      }
+
+    std::vector<std::pair<int, int>> indices_ranks_vector;
+    for (unsigned int i = 0; i < indices_ranks.extent(0); ++i)
+      {
+        indices_ranks_vector.emplace_back(indices_ranks(i).first,
+                                          indices_ranks(i).second);
+      }
+
+    std::vector<int> offsets_vector;
+    offsets_vector.insert(offsets_vector.begin(),
+                          offsets.data(),
+                          offsets.data() + offsets.extent(0));
+
+    return {indices_ranks_vector, offsets_vector};
+  }
+
+#  endif
 } // namespace ArborXWrappers
 
 DEAL_II_NAMESPACE_CLOSE
index b191ca27540b787187803b8bab8ce346e4ba33ba..bead71f741b2178192d09c5e08a90806b5ee6567 100644 (file)
@@ -20,6 +20,8 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace ArborXWrappers
 {
+#  if ARBORX_VERSION_MAJOR < 2
+
   namespace
   {
     template <int dim, typename Number>
@@ -213,12 +215,200 @@ namespace ArborXWrappers
   {
     return n_nearest_neighbors;
   }
+#  else
+  template <int dim, typename Number>
+  PointIntersectPredicate<dim, Number>::PointIntersectPredicate(
+    const std::vector<dealii::Point<dim, Number>> &points)
+    : points(points)
+  {}
+
+
+
+  template <int dim, typename Number>
+  std::size_t
+  PointIntersectPredicate<dim, Number>::size() const
+  {
+    return points.size();
+  }
+
+
+
+  template <int dim, typename Number>
+  const dealii::Point<dim, Number> &
+  PointIntersectPredicate<dim, Number>::get(unsigned int i) const
+  {
+    return points[i];
+  }
+
+
+  template <int dim, typename Number>
+  PointNearestPredicate<dim, Number>::PointNearestPredicate(
+    const std::vector<dealii::Point<dim, Number>> &points,
+    const unsigned int                             n_nearest_neighbors)
+    : points(points)
+    , n_nearest_neighbors(n_nearest_neighbors)
+  {}
+
+
+
+  template <int dim, typename Number>
+  unsigned int
+  PointNearestPredicate<dim, Number>::get_n_nearest_neighbors() const
+  {
+    return n_nearest_neighbors;
+  }
+
+
+
+  template <int dim, typename Number>
+  std::size_t
+  PointNearestPredicate<dim, Number>::size() const
+  {
+    return points.size();
+  }
+
+
+
+  template <int dim, typename Number>
+  const dealii::Point<dim, Number> &
+  PointNearestPredicate<dim, Number>::get(unsigned int i) const
+  {
+    return points[i];
+  }
+
+
+  // ------------------- BoundingBoxPredicate ------------------- //
+  template <int dim, typename Number>
+  BoundingBoxIntersectPredicate<dim, Number>::BoundingBoxIntersectPredicate(
+    const std::vector<dealii::BoundingBox<dim, Number>> &bb)
+    : bounding_boxes(bb)
+  {}
+
+
+
+  template <int dim, typename Number>
+  std::size_t
+  BoundingBoxIntersectPredicate<dim, Number>::size() const
+  {
+    return bounding_boxes.size();
+  }
+
+
+
+  template <int dim, typename Number>
+  const dealii::BoundingBox<dim, Number> &
+  BoundingBoxIntersectPredicate<dim, Number>::get(unsigned int i) const
+  {
+    return bounding_boxes[i];
+  }
+
+
+
+  template <int dim, typename Number>
+  BoundingBoxNearestPredicate<dim, Number>::BoundingBoxNearestPredicate(
+    const std::vector<dealii::BoundingBox<dim, Number>> &bounding_boxes,
+    const unsigned int                                   n_nearest_neighbors)
+    : bounding_boxes(bounding_boxes)
+    , n_nearest_neighbors(n_nearest_neighbors)
+  {}
+
+
+
+  template <int dim, typename Number>
+  unsigned int
+  BoundingBoxNearestPredicate<dim, Number>::get_n_nearest_neighbors() const
+  {
+    return n_nearest_neighbors;
+  }
+
+
+
+  template <int dim, typename Number>
+  std::size_t
+  BoundingBoxNearestPredicate<dim, Number>::size() const
+  {
+    return bounding_boxes.size();
+  }
+
+
+
+  template <int dim, typename Number>
+  const dealii::BoundingBox<dim, Number> &
+  BoundingBoxNearestPredicate<dim, Number>::get(unsigned int i) const
+  {
+    return bounding_boxes[i];
+  }
+
+  // ------------------- SpherePredicate ------------------- //
+  template <int dim, typename Number>
+  SphereIntersectPredicate<dim, Number>::SphereIntersectPredicate(
+    const std::vector<std::pair<dealii::Point<dim, Number>, Number>> &spheres)
+    : spheres(spheres)
+  {}
+
+
+
+  template <int dim, typename Number>
+  std::size_t
+  SphereIntersectPredicate<dim, Number>::size() const
+  {
+    return spheres.size();
+  }
+
+
+
+  template <int dim, typename Number>
+  const std::pair<dealii::Point<dim, Number>, Number> &
+  SphereIntersectPredicate<dim, Number>::get(unsigned int i) const
+  {
+    return spheres[i];
+  }
+
+
+
+  template <int dim, typename Number>
+  SphereNearestPredicate<dim, Number>::SphereNearestPredicate(
+    const std::vector<std::pair<dealii::Point<dim, Number>, Number>> &spheres,
+    const unsigned int n_nearest_neighbors)
+    : spheres(spheres)
+    , n_nearest_neighbors(n_nearest_neighbors)
+  {}
+
+
+
+  template <int dim, typename Number>
+  unsigned int
+  SphereNearestPredicate<dim, Number>::get_n_nearest_neighbors() const
+  {
+    return n_nearest_neighbors;
+  }
+
+
+
+  template <int dim, typename Number>
+  std::size_t
+  SphereNearestPredicate<dim, Number>::size() const
+  {
+    return spheres.size();
+  }
+
+
+
+  template <int dim, typename Number>
+  const std::pair<dealii::Point<dim, Number>, Number> &
+  SphereNearestPredicate<dim, Number>::get(unsigned int i) const
+  {
+    return spheres[i];
+  }
+#  endif
 } // namespace ArborXWrappers
 
+
 DEAL_II_NAMESPACE_CLOSE
 
 namespace ArborX
 {
+#  if ARBORX_VERSION_MAJOR < 2
   namespace
   {
     template <int dim, typename Number>
@@ -310,6 +500,24 @@ namespace ArborX
     // the sphere is 3d
     return {to_arborx_point(v[i].first), static_cast<float>(v[i].second)};
   }
+
+#  else
+  template <typename T>
+  std::size_t
+  AccessTraits<std::vector<T>>::size(const std::vector<T> &v)
+  {
+    return v.size();
+  }
+
+
+
+  template <typename T>
+  T
+  AccessTraits<std::vector<T>>::get(const std::vector<T> &v, std::size_t i)
+  {
+    return v[i];
+  }
+#  endif
 } // namespace ArborX
 
 // ----------------------- Instantiations --------------------//
index 902325d33a15f9980cec967faf62e6c10187eb11..e64a50d80b55e5be4498f5fdd4fa61c71b920d55 100644 (file)
@@ -18,6 +18,7 @@ for (DIM : DIMENSIONS; SCALAR : REAL_SCALARS)
     DEAL_II_NAMESPACE_OPEN
     namespace ArborXWrappers
     {
+#if ARBORX_VERSION_MAJOR < 2
       template PointIntersectPredicate::PointIntersectPredicate(
         const std::vector<dealii::Point<DIM, SCALAR>> &points);
       template PointNearestPredicate::PointNearestPredicate(
@@ -37,11 +38,23 @@ for (DIM : DIMENSIONS; SCALAR : REAL_SCALARS)
         const std::vector<std::pair<dealii::Point<DIM, SCALAR>, SCALAR>>
                           &spheres,
         const unsigned int n_nearest_neighbors);
+#else
+      template class PointIntersectPredicate<DIM, SCALAR>;
+      template class PointNearestPredicate<DIM, SCALAR>;
+
+      template class BoundingBoxIntersectPredicate<DIM, SCALAR>;
+      template class BoundingBoxNearestPredicate<DIM, SCALAR>;
+
+      template class SphereIntersectPredicate<DIM, SCALAR>;
+      template class SphereNearestPredicate<DIM, SCALAR>;
+#endif
     \}
     DEAL_II_NAMESPACE_CLOSE
 
     namespace ArborX
     {
+
+#if ARBORX_VERSION_MAJOR < 2
       template struct AccessTraits<std::vector<dealii::Point<DIM, SCALAR>>,
                                    PrimitivesTag>;
 
@@ -52,5 +65,14 @@ for (DIM : DIMENSIONS; SCALAR : REAL_SCALARS)
       template struct AccessTraits<
         std::vector<std::pair<dealii::Point<DIM, SCALAR>, SCALAR>>,
         PrimitivesTag>;
+#else
+      template struct AccessTraits<std::vector<dealii::Point<DIM, SCALAR>>>;
+
+      template struct AccessTraits<
+        std::vector<dealii::BoundingBox<DIM, SCALAR>>>;
+
+      template struct AccessTraits<
+        std::vector<std::pair<dealii::Point<DIM, SCALAR>, SCALAR>>>;
+#endif
     \}
   }
index 644eb0a754ceae87864d2c3d1b1b9ef78c03bca5..b49a7bc199df14b0d04ac34d52615364da8d3d99 100644 (file)
@@ -45,7 +45,11 @@ test_1d()
   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);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(bounding_boxes);
+#else
+  ArborXWrappers::BVH<BoundingBox<1>> bvh(bounding_boxes);
+#endif
   ArborXWrappers::BoundingBoxIntersectPredicate bb_intersect(
     query_bounding_boxes);
   auto             indices_offset = bvh.query(bb_intersect);
@@ -115,7 +119,11 @@ test_2d()
   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);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(bounding_boxes);
+#else
+  ArborXWrappers::BVH<BoundingBox<2>> bvh(bounding_boxes);
+#endif
   ArborXWrappers::BoundingBoxIntersectPredicate bb_intersect(
     query_bounding_boxes);
   auto             indices_offset = bvh.query(bb_intersect);
@@ -193,7 +201,11 @@ test_3d()
   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);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(bounding_boxes);
+#else
+  ArborXWrappers::BVH<BoundingBox<3>> bvh(bounding_boxes);
+#endif
   ArborXWrappers::BoundingBoxIntersectPredicate bb_intersect(
     query_bounding_boxes);
   auto             indices_offset = bvh.query(bb_intersect);
index fc220c04ba70aceda7fdc712ec97bf1be33e8c42..6d3dc6c1d4e9155185cef9870994c110b207d9cb 100644 (file)
@@ -59,7 +59,11 @@ test_bounding_box_2d()
   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);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(bounding_boxes);
+#else
+  ArborXWrappers::BVH<BoundingBox<2>> bvh(bounding_boxes);
+#endif
   ArborXWrappers::BoundingBoxNearestPredicate bb_nearest(query_bounding_boxes,
                                                          1);
   auto             indices_offset = bvh.query(bb_nearest);
@@ -108,16 +112,30 @@ test_point_2d()
         }
     }
 
+#if ARBORX_VERSION_MAJOR < 2
   Point<3> point_a(0.5, 0.5, 0.5);
   Point<3> point_b(1.5, 1.5, 1.5);
   Point<3> point_c(2.2, 2.2, 2.2);
   Point<3> point_d(2.6, 2.6, 2.6);
 
   std::vector<BoundingBox<3>> query_bounding_boxes;
+#else
+  Point<2>                            point_a(0.5, 0.5);
+  Point<2>                            point_b(1.5, 1.5);
+  Point<2>                            point_c(2.2, 2.2);
+  Point<2>                            point_d(2.6, 2.6);
+
+  std::vector<BoundingBox<2>>         query_bounding_boxes;
+#endif
+
   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(points);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(points);
+#else
+  ArborXWrappers::BVH<Point<2>>       bvh(points);
+#endif
   ArborXWrappers::BoundingBoxNearestPredicate bb_nearest(query_bounding_boxes,
                                                          1);
   auto             indices_offset = bvh.query(bb_nearest);
@@ -195,7 +213,11 @@ test_bounding_box_3d()
   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);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(bounding_boxes);
+#else
+  ArborXWrappers::BVH<BoundingBox<3>> bvh(bounding_boxes);
+#endif
   ArborXWrappers::BoundingBoxNearestPredicate bb_nearest(query_bounding_boxes,
                                                          1);
   auto             indices_offset = bvh.query(bb_nearest);
@@ -255,7 +277,11 @@ test_point_3d()
   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(points);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(points);
+#else
+  ArborXWrappers::BVH<Point<3>>       bvh(points);
+#endif
   ArborXWrappers::BoundingBoxNearestPredicate bb_nearest(query_bounding_boxes,
                                                          1);
   auto             indices_offset = bvh.query(bb_nearest);
index 2735972b248ac627fc7b5f1bc13508c11f8a4631..a2c30a9a40db87d761312ce6d7aa813c010c443b 100644 (file)
@@ -44,7 +44,11 @@ test_1d()
   query_points.emplace_back(2.6);
 
 
-  ArborXWrappers::BVH                     bvh(bounding_boxes);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(bounding_boxes);
+#else
+  ArborXWrappers::BVH<BoundingBox<1>> bvh(bounding_boxes);
+#endif
   ArborXWrappers::PointIntersectPredicate pt_intersect(query_points);
   auto             indices_offset = bvh.query(pt_intersect);
   std::vector<int> indices        = indices_offset.first;
@@ -112,7 +116,11 @@ test_2d()
   query_points.emplace_back(2.6, 2.6);
 
 
-  ArborXWrappers::BVH                     bvh(bounding_boxes);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(bounding_boxes);
+#else
+  ArborXWrappers::BVH<BoundingBox<2>> bvh(bounding_boxes);
+#endif
   ArborXWrappers::PointIntersectPredicate pt_intersect(query_points);
   auto             indices_offset = bvh.query(pt_intersect);
   std::vector<int> indices        = indices_offset.first;
@@ -187,7 +195,11 @@ test_3d()
   query_points.emplace_back(2.6, 2.6, 2.6);
 
 
-  ArborXWrappers::BVH                     bvh(bounding_boxes);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(bounding_boxes);
+#else
+  ArborXWrappers::BVH<BoundingBox<3>> bvh(bounding_boxes);
+#endif
   ArborXWrappers::PointIntersectPredicate pt_intersect(query_points);
   auto             indices_offset = bvh.query(pt_intersect);
   std::vector<int> indices        = indices_offset.first;
index c1c232f010e851a97fe64aca96226bfc3572e46c..7360ad76b5a93ead4969e030faa8fb7dffba4659 100644 (file)
@@ -57,7 +57,11 @@ test_bounding_box_2d()
   query_points.emplace_back(2.6, 2.6);
 
 
-  ArborXWrappers::BVH                   bvh(bounding_boxes);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(bounding_boxes);
+#else
+  ArborXWrappers::BVH<BoundingBox<2>> bvh(bounding_boxes);
+#endif
   ArborXWrappers::PointNearestPredicate pt_nearest(query_points, 1);
   auto                                  indices_offset = bvh.query(pt_nearest);
   std::vector<int>                      indices        = indices_offset.first;
@@ -113,7 +117,11 @@ test_points_2d()
   query_points.emplace_back(2.6, 2.6);
 
 
-  ArborXWrappers::BVH                   bvh(points);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(points);
+#else
+  ArborXWrappers::BVH<Point<2>>       bvh(points);
+#endif
   ArborXWrappers::PointNearestPredicate pt_nearest(query_points, 1);
   auto                                  indices_offset = bvh.query(pt_nearest);
   std::vector<int>                      indices        = indices_offset.first;
@@ -188,7 +196,11 @@ test_bounding_box_3d()
   query_points.emplace_back(2.6, 2.6, 2.6);
 
 
-  ArborXWrappers::BVH                   bvh(bounding_boxes);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(bounding_boxes);
+#else
+  ArborXWrappers::BVH<BoundingBox<3>> bvh(bounding_boxes);
+#endif
   ArborXWrappers::PointNearestPredicate pt_nearest(query_points, 1);
   auto                                  indices_offset = bvh.query(pt_nearest);
   std::vector<int>                      indices        = indices_offset.first;
index ff056a4eb562b417ea04e90ad97adff089d5cffd..6f88c1d061ebfa898e24202b13c70bf6b4534f51 100644 (file)
@@ -36,11 +36,15 @@ test_1d()
   std::vector<std::pair<Point<1>, double>> query_spheres;
   query_spheres.emplace_back(Point<1>(0.5), 0.8);
   query_spheres.emplace_back(Point<1>(1.5), 0.8);
-  query_spheres.emplace_back(Point<1>(2.2), 1.2);
+  query_spheres.emplace_back(Point<1>(2.2), 1.2001);
   query_spheres.emplace_back(Point<1>(2.6), 0.9);
 
 
-  ArborXWrappers::BVH                      bvh(points);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(points);
+#else
+  ArborXWrappers::BVH<Point<1>> bvh(points);
+#endif
   ArborXWrappers::SphereIntersectPredicate sph_intersect(query_spheres);
   auto             indices_offsets = bvh.query(sph_intersect);
   std::vector<int> indices         = indices_offsets.first;
@@ -49,6 +53,9 @@ test_1d()
   std::vector<int> indices_ref = {0, 1, 1, 2, 1, 2, 3, 2, 3};
   std::vector<int> offsets_ref = {0, 2, 4, 7, 9};
 
+  std::cout << indices.size() << " " << indices_ref.size() << std::endl;
+  for (auto &val : indices)
+    std::cout << val << std::endl;
   AssertThrow(indices.size() == indices_ref.size(), ExcInternalError());
   AssertThrow(offsets.size() == offsets_ref.size(), ExcInternalError());
   for (unsigned int i = 0; i < offsets.size() - 1; ++i)
@@ -95,7 +102,11 @@ test_2d()
   query_spheres.push_back(std::make_pair(Point<2>(2.6, 2.6), 0.9));
 
 
-  ArborXWrappers::BVH                      bvh(points);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(points);
+#else
+  ArborXWrappers::BVH<Point<2>> bvh(points);
+#endif
   ArborXWrappers::SphereIntersectPredicate sph_intersect(query_spheres);
   auto             indices_offsets = bvh.query(sph_intersect);
   std::vector<int> indices         = indices_offsets.first;
@@ -154,7 +165,11 @@ test_3d()
   query_spheres.push_back(std::make_pair(Point<3>(2.6, 2.6, 2.6), 1.1));
 
 
-  ArborXWrappers::BVH                      bvh(points);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(points);
+#else
+  ArborXWrappers::BVH<Point<3>> bvh(points);
+#endif
   ArborXWrappers::SphereIntersectPredicate sph_intersect(query_spheres);
   auto             indices_offsets = bvh.query(sph_intersect);
   std::vector<int> indices         = indices_offsets.first;
@@ -198,7 +213,14 @@ main(int argc, char **argv)
   // Initialize ArborX
   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
 
-  // tests
+  // The 1D test hits a bug in clang:
+  // https://github.com/llvm/llvm-project/issues/18060
+#if defined(DEAL_II_HAVE_FP_EXCEPTIONS)
+  {
+    const int current_fe_except = fegetexcept();
+    fedisableexcept(current_fe_except);
+  }
+#endif
   test_1d();
   test_2d();
   test_3d();
index d6ae73d052f514c6384b646cb8fc1e9967023578..028023a87ea3737b3ed8bb01069fc14c395ff9b9 100644 (file)
 void
 test_1d()
 {
-  std::vector<Point<1>> points;
+  std::vector<Point<1, float>> points;
 
   unsigned int n_points_1d = 5;
   for (unsigned int i = 0; i < n_points_1d; ++i)
     points.emplace_back(i);
 
-  std::vector<std::pair<Point<1>, double>> query_spheres;
-  query_spheres.emplace_back(Point<1>(0.5), 0.1);
-  query_spheres.emplace_back(Point<1>(1.5), 0.1);
-  query_spheres.emplace_back(Point<1>(2.2), 0.1);
-  query_spheres.emplace_back(Point<1>(2.6), 0.1);
+  std::vector<std::pair<Point<1, float>, float>> query_spheres;
+  query_spheres.emplace_back(Point<1, float>(0.5), 0.1);
+  query_spheres.emplace_back(Point<1, float>(1.5), 0.1);
+  query_spheres.emplace_back(Point<1, float>(2.2), 0.1);
+  query_spheres.emplace_back(Point<1, float>(2.6), 0.1);
 
 
-  ArborXWrappers::BVH                    bvh(points);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(points);
+#else
+  std::vector<BoundingBox<1, float>> bounding_boxes;
+  for (const auto &p : points)
+    {
+      bounding_boxes.emplace_back(std::pair(p, p));
+    }
+  ArborXWrappers::BVH<BoundingBox<1, float>> bvh(bounding_boxes);
+#endif
   ArborXWrappers::SphereNearestPredicate sph_nearest(query_spheres, 1);
   auto             indices_offsets = bvh.query(sph_nearest);
   std::vector<int> indices         = indices_offsets.first;
@@ -79,7 +88,7 @@ test_1d()
 void
 test_2d()
 {
-  std::vector<Point<2>> points;
+  std::vector<Point<2, float>> points;
 
   unsigned int n_points_1d = 5;
   for (unsigned int i = 0; i < n_points_1d; ++i)
@@ -91,14 +100,23 @@ test_2d()
     }
 
 
-  std::vector<std::pair<Point<2>, double>> query_spheres;
-  query_spheres.push_back(std::make_pair(Point<2>(0.5, 0.5), 0.1));
-  query_spheres.push_back(std::make_pair(Point<2>(1.5, 1.5), 0.1));
-  query_spheres.push_back(std::make_pair(Point<2>(2.2, 2.2), 0.1));
-  query_spheres.push_back(std::make_pair(Point<2>(2.6, 2.6), 0.1));
+  std::vector<std::pair<Point<2, float>, float>> query_spheres;
+  query_spheres.push_back(std::make_pair(Point<2, float>(0.5, 0.5), 0.1));
+  query_spheres.push_back(std::make_pair(Point<2, float>(1.5, 1.5), 0.1));
+  query_spheres.push_back(std::make_pair(Point<2, float>(2.2, 2.2), 0.1));
+  query_spheres.push_back(std::make_pair(Point<2, float>(2.6, 2.6), 0.1));
 
 
-  ArborXWrappers::BVH                    bvh(points);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(points);
+#else
+  std::vector<BoundingBox<2, float>>         bounding_boxes;
+  for (const auto &p : points)
+    {
+      bounding_boxes.emplace_back(std::pair(p, p));
+    }
+  ArborXWrappers::BVH<BoundingBox<2, float>> bvh(bounding_boxes);
+#endif
   ArborXWrappers::SphereNearestPredicate sph_nearest(query_spheres, 1);
   auto             indices_offsets = bvh.query(sph_nearest);
   std::vector<int> indices         = indices_offsets.first;
@@ -136,7 +154,7 @@ test_2d()
 void
 test_3d()
 {
-  std::vector<Point<3>> points;
+  std::vector<Point<3, float>> points;
 
   unsigned int n_points_1d = 5;
   for (unsigned int i = 0; i < n_points_1d; ++i)
@@ -150,14 +168,23 @@ test_3d()
         }
     }
 
-  std::vector<std::pair<Point<3>, double>> query_spheres;
-  query_spheres.push_back(std::make_pair(Point<3>(0.5, 0.5, 0.5), 0.1));
-  query_spheres.push_back(std::make_pair(Point<3>(1.5, 1.5, 1.5), 0.1));
-  query_spheres.push_back(std::make_pair(Point<3>(2.2, 2.2, 2.2), 0.1));
-  query_spheres.push_back(std::make_pair(Point<3>(2.6, 2.6, 2.6), 0.1));
+  std::vector<std::pair<Point<3, float>, float>> query_spheres;
+  query_spheres.push_back(std::make_pair(Point<3, float>(0.5, 0.5, 0.5), 0.1));
+  query_spheres.push_back(std::make_pair(Point<3, float>(1.5, 1.5, 1.5), 0.1));
+  query_spheres.push_back(std::make_pair(Point<3, float>(2.2, 2.2, 2.2), 0.1));
+  query_spheres.push_back(std::make_pair(Point<3, float>(2.6, 2.6, 2.6), 0.1));
 
 
-  ArborXWrappers::BVH                    bvh(points);
+#if ARBORX_VERSION_MAJOR < 2
+  ArborXWrappers::BVH bvh(points);
+#else
+  std::vector<BoundingBox<3, float>>         bounding_boxes;
+  for (const auto &p : points)
+    {
+      bounding_boxes.emplace_back(std::pair(p, p));
+    }
+  ArborXWrappers::BVH<BoundingBox<3, float>> bvh(bounding_boxes);
+#endif
   ArborXWrappers::SphereNearestPredicate sph_nearest(query_spheres, 1);
   auto             indices_offsets = bvh.query(sph_nearest);
   std::vector<int> indices         = indices_offsets.first;
@@ -198,7 +225,14 @@ main(int argc, char **argv)
   // Initialize ArborX
   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
 
-  // tests
+  // The 1D test hits a bug in clang:
+  // https://github.com/llvm/llvm-project/issues/18060
+#if defined(DEAL_II_HAVE_FP_EXCEPTIONS)
+  {
+    const int current_fe_except = fegetexcept();
+    fedisableexcept(current_fe_except);
+  }
+#endif
   test_1d();
   test_2d();
   test_3d();

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.