From e1fad5eff43d7db37fbbccf0390c8f1423f9fa7a Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Tue, 17 Jun 2025 16:47:57 -0400 Subject: [PATCH] Add wrappers for ArborX 2.0 --- include/deal.II/arborx/access_traits.h | 828 +++++++++++++++++++++- include/deal.II/arborx/bvh.h | 134 +++- include/deal.II/arborx/distributed_tree.h | 109 ++- source/arborx/access_traits.cc | 208 ++++++ source/arborx/access_traits.inst.in | 22 + tests/arborx/bounding_box_intersect.cc | 18 +- tests/arborx/bounding_box_nearest.cc | 34 +- tests/arborx/point_intersect.cc | 18 +- tests/arborx/point_nearest.cc | 18 +- tests/arborx/sphere_intersect.cc | 32 +- tests/arborx/sphere_nearest.cc | 78 +- 11 files changed, 1447 insertions(+), 52 deletions(-) diff --git a/include/deal.II/arborx/access_traits.h b/include/deal.II/arborx/access_traits.h index 88cb92f09f..9ba7212d3a 100644 --- a/include/deal.II/arborx/access_traits.h +++ b/include/deal.II/arborx/access_traits.h @@ -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 + ArborX::Point + operator()(const ArborX::PairValueIndex, + 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 + ArborX::Box + operator()(const ArborX::PairValueIndex, + unsigned int> &pair) const + { + const auto boundary_points = pair.value.get_boundary_points(); + const dealii::Point min_corner = boundary_points.first; + const dealii::Point 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 + ArborX::Sphere + operator()(const ArborX::PairValueIndex< + std::pair, 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 + 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 + KOKKOS_FUNCTION void + operator()(const Predicate &, + const ArborX::PairValueIndex &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 + 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> &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 & + 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> 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 + 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> &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 & + 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> 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 + 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> &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 & + 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> 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 + 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> &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 & + 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> 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 + 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, 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, 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, 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 + 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, 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, 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, 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 as * primitive. @@ -383,9 +763,32 @@ namespace ArborX get(const std::vector, Number>> &v, std::size_t i); }; +# else + /** + * This struct allows ArborX to use std::vector as primitive. + */ + template + struct AccessTraits> + { + using memory_space = Kokkos::HostSpace; + + /** + * Return the size of the vector @p v. + */ + static std::size_t + size(const std::vector &v); + + /** + * Return the ith element of the vector @p v. + */ + static T + get(const std::vector &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 + struct AccessTraits< + dealii::ArborXWrappers::PointIntersectPredicate> + { + using memory_space = Kokkos::HostSpace; + + /** + * The number of points stored in @p pt_intersect. + */ + static std::size_t + size(const dealii::ArborXWrappers::PointIntersectPredicate + &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 + &pt_intersect, + std::size_t i); + }; + + + + /** + * This struct allows ArborX to use PointNearestPredicate as a predicate. + */ + template + struct AccessTraits< + dealii::ArborXWrappers::PointNearestPredicate> + { + using memory_space = Kokkos::HostSpace; + + /** + * The number of points stored in @p pt_nearest. + */ + static std::size_t + size(const dealii::ArborXWrappers::PointNearestPredicate + &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 + &pt_nearest, + std::size_t i); + }; + + + + /** + * This struct allows ArborX to use BoundingBoxIntersectPredicate as a + * predicate. + */ + template + struct AccessTraits< + dealii::ArborXWrappers::BoundingBoxIntersectPredicate> + { + using memory_space = Kokkos::HostSpace; + + /** + * The number of bounding boxes stored in @p bb_intersect. + */ + static std::size_t + size( + const dealii::ArborXWrappers::BoundingBoxIntersectPredicate + &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 + &bb_intersect, + std::size_t i); + }; + + + + /** + * This struct allows ArborX to use BoundingBoxNearstPredicate as a + * predicate. + */ + template + struct AccessTraits< + dealii::ArborXWrappers::BoundingBoxNearestPredicate> + { + using memory_space = Kokkos::HostSpace; + + /** + * The number of bounding boxes stored in @p bb_nearest. + */ + static std::size_t + size(const dealii::ArborXWrappers::BoundingBoxNearestPredicate + &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 + &bb_nearest, + std::size_t i); + }; + + + + /** + * This struct allows ArborX to use SphereIntersectPredicate as a + predicate. + */ + template + struct AccessTraits< + dealii::ArborXWrappers::SphereIntersectPredicate> + { + using memory_space = Kokkos::HostSpace; + + /** + * The number of points stored in @p sph_intersect. + */ + static std::size_t + size(const dealii::ArborXWrappers::SphereIntersectPredicate + &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 + &sph_intersect, + std::size_t i); + }; + + + + /** + * This struct allows ArborX to use SphereNearestPredicate as a predicate. + */ + template + struct AccessTraits< + dealii::ArborXWrappers::SphereNearestPredicate> + { + using memory_space = Kokkos::HostSpace; + + /** + * The number of spheres stored in @p sph_nearest. + */ + static std::size_t + size(const dealii::ArborXWrappers::SphereNearestPredicate + &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 + &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:: size(const dealii::ArborXWrappers::PointIntersectPredicate &pt_intersect) @@ -695,6 +1272,247 @@ namespace ArborX sph_nearest.get_n_nearest_neighbors()); } } // namespace ArborX +# else + template + inline std::size_t + AccessTraits>:: + size(const dealii::ArborXWrappers::PointIntersectPredicate + &pt_intersect) + { + return pt_intersect.size(); + } + + + + template + inline auto + AccessTraits>:: + get(const dealii::ArborXWrappers::PointIntersectPredicate + &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 + inline std::size_t + AccessTraits>:: + size(const dealii::ArborXWrappers::PointNearestPredicate + &pt_nearest) + { + return pt_nearest.size(); + } + + + + template + inline auto + AccessTraits>::get( + const dealii::ArborXWrappers::PointNearestPredicate + &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 + inline std::size_t + AccessTraits< + dealii::ArborXWrappers::BoundingBoxIntersectPredicate>:: + size( + const dealii::ArborXWrappers::BoundingBoxIntersectPredicate + &bb_intersect) + { + return bb_intersect.size(); + } + + + + template + inline auto + AccessTraits< + dealii::ArborXWrappers::BoundingBoxIntersectPredicate>:: + get(const dealii::ArborXWrappers::BoundingBoxIntersectPredicate + &bb_intersect, + std::size_t i) + { + const auto boundary_points = bb_intersect.get(i).get_boundary_points(); + const dealii::Point min_corner = boundary_points.first; + const dealii::Point 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 + inline std::size_t + AccessTraits< + dealii::ArborXWrappers::BoundingBoxNearestPredicate>:: + size(const dealii::ArborXWrappers::BoundingBoxNearestPredicate + &bb_nearest) + { + return bb_nearest.size(); + } + + + + template + inline auto + AccessTraits< + dealii::ArborXWrappers::BoundingBoxNearestPredicate>:: + get(const dealii::ArborXWrappers::BoundingBoxNearestPredicate + &bb_nearest, + std::size_t i) + { + const auto boundary_points = bb_nearest.get(i).get_boundary_points(); + const dealii::Point min_corner = boundary_points.first; + const dealii::Point 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 + inline std::size_t + AccessTraits>:: + size(const dealii::ArborXWrappers::SphereIntersectPredicate + &sph_intersect) + { + return sph_intersect.size(); + } + + + + template + inline auto + AccessTraits>:: + get(const dealii::ArborXWrappers::SphereIntersectPredicate + &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 + inline std::size_t + AccessTraits>:: + size(const dealii::ArborXWrappers::SphereNearestPredicate + &sph_nearest) + { + return sph_nearest.size(); + } + + + + template + inline auto + AccessTraits>:: + get(const dealii::ArborXWrappers::SphereNearestPredicate + &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 diff --git a/include/deal.II/arborx/bvh.h b/include/deal.II/arborx/bvh.h index 6819cdd6ec..ff639ee7e7 100644 --- a/include/deal.II/arborx/bvh.h +++ b/include/deal.II/arborx/bvh.h @@ -20,6 +20,7 @@ #ifdef DEAL_II_WITH_ARBORX # include +# include # include # include @@ -46,10 +47,8 @@ namespace ArborXWrappers * several operations on sets of geometric objects efficiently, such as in * collision detection and ray tracing. * - * - * 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 + class BVH + { + public: + /** + * Constructor. Use a vector of Point, BoundingBox, or sphere (std::pair) @p values as primitives. + */ + BVH(const std::vector &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> query_bounding_boxes = ... + * ArborXWrappers::BoundingBoxIntersectPredicate + * bb_intersect(query_bounding_boxes); + * + * const std::vector> bvh_bounding_boxes = ... + * ArborxWrappers::BVH> 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 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> query_points = ... + * ArborXWrappers::PointIntersectPredicate pt_intersect(query_points); + * + * const std::vector> bvh_bounding_boxes = ... + * ArborxWrappers::BVH> 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> query_points = ... + * ArborXWrappers::PointNearestPredicate pt_nearest(query_points, 5); + * + * const std::vector> bvh_points = ... + * ArborxWrappers::BVH> bvh(bvh_points); + * + * auto [indices, offset] = bvh.query(pt_nearest); + * @endcode + */ + template + std::pair, std::vector> + query(const QueryType &queries); + + private: + /** + * Underlying ArborX object. + */ + ArborX::BVH, + internal::IndexableGetter> + bvh; + }; + + + + template + BVH::BVH(const std::vector &values) + : bvh(Kokkos::DefaultHostExecutionSpace{}, + ArborX::Experimental::attach_indices(values), + internal::IndexableGetter{}) + {} + + + + template + template + std::pair, std::vector> + BVH::query(const QueryType &queries) + { + Kokkos::View indices("indices", 0); + + Kokkos::View offset("offset", 0); + bvh.query(Kokkos::DefaultHostExecutionSpace{}, + queries, + internal::ExtractIndex{}, + indices, + offset); + std::vector indices_vector; + indices_vector.insert(indices_vector.begin(), + indices.data(), + indices.data() + indices.extent(0)); + std::vector 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 diff --git a/include/deal.II/arborx/distributed_tree.h b/include/deal.II/arborx/distributed_tree.h index c02bd9d2fe..728a0bf6b8 100644 --- a/include/deal.II/arborx/distributed_tree.h +++ b/include/deal.II/arborx/distributed_tree.h @@ -20,6 +20,8 @@ #if defined(DEAL_II_ARBORX_WITH_MPI) && defined(DEAL_II_WITH_MPI) # include +# include + # include # include @@ -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 std::pair>, std::vector> @@ -125,6 +127,105 @@ namespace ArborXWrappers return {indices_ranks_vector, offsets_vector}; } +# else + template + class DistributedTree + { + public: + /** + * Constructor. Use a vector of Point, BoundingBox, or sphere + * (std::pair) @p values as primitives. The objects in @p values + * are local to the MPI process. + */ + DistributedTree(const MPI_Comm comm, const std::vector &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 + std::pair>, std::vector> + query(const QueryType &queries); + + private: + /** + * Underlying ArborX object. + */ + ArborX::DistributedTree, + internal::IndexableGetter> + distributed_tree; + /** + * Callback to extract the index and the rank of each primitive that + * satisfies a query. + */ + internal::ExtractIndexRank callback; + }; + + + + template + DistributedTree::DistributedTree(const MPI_Comm comm, + const std::vector &values) + : distributed_tree(comm, + Kokkos::DefaultHostExecutionSpace{}, + ArborX::Experimental::attach_indices(values), + internal::IndexableGetter{}) + , callback{Utilities::MPI::this_mpi_process(comm)} + {} + + template + template + std::pair>, std::vector> + DistributedTree::query(const QueryType &queries) + { + Kokkos::View offsets("offsets", 0); + Kokkos::View *, 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> 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 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 diff --git a/source/arborx/access_traits.cc b/source/arborx/access_traits.cc index b191ca2754..bead71f741 100644 --- a/source/arborx/access_traits.cc +++ b/source/arborx/access_traits.cc @@ -20,6 +20,8 @@ DEAL_II_NAMESPACE_OPEN namespace ArborXWrappers { +# if ARBORX_VERSION_MAJOR < 2 + namespace { template @@ -213,12 +215,200 @@ namespace ArborXWrappers { return n_nearest_neighbors; } +# else + template + PointIntersectPredicate::PointIntersectPredicate( + const std::vector> &points) + : points(points) + {} + + + + template + std::size_t + PointIntersectPredicate::size() const + { + return points.size(); + } + + + + template + const dealii::Point & + PointIntersectPredicate::get(unsigned int i) const + { + return points[i]; + } + + + template + PointNearestPredicate::PointNearestPredicate( + const std::vector> &points, + const unsigned int n_nearest_neighbors) + : points(points) + , n_nearest_neighbors(n_nearest_neighbors) + {} + + + + template + unsigned int + PointNearestPredicate::get_n_nearest_neighbors() const + { + return n_nearest_neighbors; + } + + + + template + std::size_t + PointNearestPredicate::size() const + { + return points.size(); + } + + + + template + const dealii::Point & + PointNearestPredicate::get(unsigned int i) const + { + return points[i]; + } + + + // ------------------- BoundingBoxPredicate ------------------- // + template + BoundingBoxIntersectPredicate::BoundingBoxIntersectPredicate( + const std::vector> &bb) + : bounding_boxes(bb) + {} + + + + template + std::size_t + BoundingBoxIntersectPredicate::size() const + { + return bounding_boxes.size(); + } + + + + template + const dealii::BoundingBox & + BoundingBoxIntersectPredicate::get(unsigned int i) const + { + return bounding_boxes[i]; + } + + + + template + BoundingBoxNearestPredicate::BoundingBoxNearestPredicate( + const std::vector> &bounding_boxes, + const unsigned int n_nearest_neighbors) + : bounding_boxes(bounding_boxes) + , n_nearest_neighbors(n_nearest_neighbors) + {} + + + + template + unsigned int + BoundingBoxNearestPredicate::get_n_nearest_neighbors() const + { + return n_nearest_neighbors; + } + + + + template + std::size_t + BoundingBoxNearestPredicate::size() const + { + return bounding_boxes.size(); + } + + + + template + const dealii::BoundingBox & + BoundingBoxNearestPredicate::get(unsigned int i) const + { + return bounding_boxes[i]; + } + + // ------------------- SpherePredicate ------------------- // + template + SphereIntersectPredicate::SphereIntersectPredicate( + const std::vector, Number>> &spheres) + : spheres(spheres) + {} + + + + template + std::size_t + SphereIntersectPredicate::size() const + { + return spheres.size(); + } + + + + template + const std::pair, Number> & + SphereIntersectPredicate::get(unsigned int i) const + { + return spheres[i]; + } + + + + template + SphereNearestPredicate::SphereNearestPredicate( + const std::vector, Number>> &spheres, + const unsigned int n_nearest_neighbors) + : spheres(spheres) + , n_nearest_neighbors(n_nearest_neighbors) + {} + + + + template + unsigned int + SphereNearestPredicate::get_n_nearest_neighbors() const + { + return n_nearest_neighbors; + } + + + + template + std::size_t + SphereNearestPredicate::size() const + { + return spheres.size(); + } + + + + template + const std::pair, Number> & + SphereNearestPredicate::get(unsigned int i) const + { + return spheres[i]; + } +# endif } // namespace ArborXWrappers + DEAL_II_NAMESPACE_CLOSE namespace ArborX { +# if ARBORX_VERSION_MAJOR < 2 namespace { template @@ -310,6 +500,24 @@ namespace ArborX // the sphere is 3d return {to_arborx_point(v[i].first), static_cast(v[i].second)}; } + +# else + template + std::size_t + AccessTraits>::size(const std::vector &v) + { + return v.size(); + } + + + + template + T + AccessTraits>::get(const std::vector &v, std::size_t i) + { + return v[i]; + } +# endif } // namespace ArborX // ----------------------- Instantiations --------------------// diff --git a/source/arborx/access_traits.inst.in b/source/arborx/access_traits.inst.in index 902325d33a..e64a50d80b 100644 --- a/source/arborx/access_traits.inst.in +++ b/source/arborx/access_traits.inst.in @@ -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> &points); template PointNearestPredicate::PointNearestPredicate( @@ -37,11 +38,23 @@ for (DIM : DIMENSIONS; SCALAR : REAL_SCALARS) const std::vector, SCALAR>> &spheres, const unsigned int n_nearest_neighbors); +#else + template class PointIntersectPredicate; + template class PointNearestPredicate; + + template class BoundingBoxIntersectPredicate; + template class BoundingBoxNearestPredicate; + + template class SphereIntersectPredicate; + template class SphereNearestPredicate; +#endif \} DEAL_II_NAMESPACE_CLOSE namespace ArborX { + +#if ARBORX_VERSION_MAJOR < 2 template struct AccessTraits>, PrimitivesTag>; @@ -52,5 +65,14 @@ for (DIM : DIMENSIONS; SCALAR : REAL_SCALARS) template struct AccessTraits< std::vector, SCALAR>>, PrimitivesTag>; +#else + template struct AccessTraits>>; + + template struct AccessTraits< + std::vector>>; + + template struct AccessTraits< + std::vector, SCALAR>>>; +#endif \} } diff --git a/tests/arborx/bounding_box_intersect.cc b/tests/arborx/bounding_box_intersect.cc index 644eb0a754..b49a7bc199 100644 --- a/tests/arborx/bounding_box_intersect.cc +++ b/tests/arborx/bounding_box_intersect.cc @@ -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> 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> 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> bvh(bounding_boxes); +#endif ArborXWrappers::BoundingBoxIntersectPredicate bb_intersect( query_bounding_boxes); auto indices_offset = bvh.query(bb_intersect); diff --git a/tests/arborx/bounding_box_nearest.cc b/tests/arborx/bounding_box_nearest.cc index fc220c04ba..6d3dc6c1d4 100644 --- a/tests/arborx/bounding_box_nearest.cc +++ b/tests/arborx/bounding_box_nearest.cc @@ -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> 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> 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> 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> 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> 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> bvh(points); +#endif ArborXWrappers::BoundingBoxNearestPredicate bb_nearest(query_bounding_boxes, 1); auto indices_offset = bvh.query(bb_nearest); diff --git a/tests/arborx/point_intersect.cc b/tests/arborx/point_intersect.cc index 2735972b24..a2c30a9a40 100644 --- a/tests/arborx/point_intersect.cc +++ b/tests/arborx/point_intersect.cc @@ -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> bvh(bounding_boxes); +#endif ArborXWrappers::PointIntersectPredicate pt_intersect(query_points); auto indices_offset = bvh.query(pt_intersect); std::vector 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> bvh(bounding_boxes); +#endif ArborXWrappers::PointIntersectPredicate pt_intersect(query_points); auto indices_offset = bvh.query(pt_intersect); std::vector 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> bvh(bounding_boxes); +#endif ArborXWrappers::PointIntersectPredicate pt_intersect(query_points); auto indices_offset = bvh.query(pt_intersect); std::vector indices = indices_offset.first; diff --git a/tests/arborx/point_nearest.cc b/tests/arborx/point_nearest.cc index c1c232f010..7360ad76b5 100644 --- a/tests/arborx/point_nearest.cc +++ b/tests/arborx/point_nearest.cc @@ -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> bvh(bounding_boxes); +#endif ArborXWrappers::PointNearestPredicate pt_nearest(query_points, 1); auto indices_offset = bvh.query(pt_nearest); std::vector 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> bvh(points); +#endif ArborXWrappers::PointNearestPredicate pt_nearest(query_points, 1); auto indices_offset = bvh.query(pt_nearest); std::vector 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> bvh(bounding_boxes); +#endif ArborXWrappers::PointNearestPredicate pt_nearest(query_points, 1); auto indices_offset = bvh.query(pt_nearest); std::vector indices = indices_offset.first; diff --git a/tests/arborx/sphere_intersect.cc b/tests/arborx/sphere_intersect.cc index ff056a4eb5..6f88c1d061 100644 --- a/tests/arborx/sphere_intersect.cc +++ b/tests/arborx/sphere_intersect.cc @@ -36,11 +36,15 @@ test_1d() std::vector, 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> bvh(points); +#endif ArborXWrappers::SphereIntersectPredicate sph_intersect(query_spheres); auto indices_offsets = bvh.query(sph_intersect); std::vector indices = indices_offsets.first; @@ -49,6 +53,9 @@ test_1d() std::vector indices_ref = {0, 1, 1, 2, 1, 2, 3, 2, 3}; std::vector 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> bvh(points); +#endif ArborXWrappers::SphereIntersectPredicate sph_intersect(query_spheres); auto indices_offsets = bvh.query(sph_intersect); std::vector 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> bvh(points); +#endif ArborXWrappers::SphereIntersectPredicate sph_intersect(query_spheres); auto indices_offsets = bvh.query(sph_intersect); std::vector 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(); diff --git a/tests/arborx/sphere_nearest.cc b/tests/arborx/sphere_nearest.cc index d6ae73d052..028023a87e 100644 --- a/tests/arborx/sphere_nearest.cc +++ b/tests/arborx/sphere_nearest.cc @@ -28,20 +28,29 @@ void test_1d() { - std::vector> points; + std::vector> points; unsigned int n_points_1d = 5; for (unsigned int i = 0; i < n_points_1d; ++i) points.emplace_back(i); - std::vector, 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, 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> bounding_boxes; + for (const auto &p : points) + { + bounding_boxes.emplace_back(std::pair(p, p)); + } + ArborXWrappers::BVH> bvh(bounding_boxes); +#endif ArborXWrappers::SphereNearestPredicate sph_nearest(query_spheres, 1); auto indices_offsets = bvh.query(sph_nearest); std::vector indices = indices_offsets.first; @@ -79,7 +88,7 @@ test_1d() void test_2d() { - std::vector> points; + std::vector> 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, 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, 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> bounding_boxes; + for (const auto &p : points) + { + bounding_boxes.emplace_back(std::pair(p, p)); + } + ArborXWrappers::BVH> bvh(bounding_boxes); +#endif ArborXWrappers::SphereNearestPredicate sph_nearest(query_spheres, 1); auto indices_offsets = bvh.query(sph_nearest); std::vector indices = indices_offsets.first; @@ -136,7 +154,7 @@ test_2d() void test_3d() { - std::vector> points; + std::vector> 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, 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, 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> bounding_boxes; + for (const auto &p : points) + { + bounding_boxes.emplace_back(std::pair(p, p)); + } + ArborXWrappers::BVH> bvh(bounding_boxes); +#endif ArborXWrappers::SphereNearestPredicate sph_nearest(query_spheres, 1); auto indices_offsets = bvh.query(sph_nearest); std::vector 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(); -- 2.39.5