From 70ec02b82837746ce87b83ffbc488ce1de49be03 Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Fri, 19 Mar 2021 21:51:34 +0000 Subject: [PATCH] Refactor ArborX access traits and add nearest neighbor predicate --- include/deal.II/arborx/access_traits.h | 268 +++++++++++++++++--- include/deal.II/arborx/bvh.h | 31 ++- source/arborx/access_traits.cc | 66 ++++- source/arborx/access_traits.inst.in | 13 +- tests/arborx/CMakeLists.txt | 2 +- tests/arborx/bounding_box_intersect.cc | 24 +- tests/arborx/bounding_box_nearest.cc | 308 +++++++++++++++++++++++ tests/arborx/bounding_box_nearest.output | 5 + tests/arborx/point_intersect.cc | 20 +- tests/arborx/point_nearest.cc | 301 ++++++++++++++++++++++ tests/arborx/point_nearest.output | 5 + 11 files changed, 976 insertions(+), 67 deletions(-) create mode 100644 tests/arborx/bounding_box_nearest.cc create mode 100644 tests/arborx/bounding_box_nearest.output create mode 100644 tests/arborx/point_nearest.cc create mode 100644 tests/arborx/point_nearest.output diff --git a/include/deal.II/arborx/access_traits.h b/include/deal.II/arborx/access_traits.h index ce9204c1f8..7c02fd85c2 100644 --- a/include/deal.II/arborx/access_traits.h +++ b/include/deal.II/arborx/access_traits.h @@ -29,18 +29,16 @@ DEAL_II_NAMESPACE_OPEN namespace ArborXWrappers { /** - * This class is used by ArborXWrappers::BVH to determine which points - * intersect the bounding boxes used to build the ArborXWrappers::BVH. + * Base class for Point-based predicates. */ - class PointIntersect + class PointPredicate { public: /** - * Constructor. @p dim_points is a list of points which we are interested in - * knowing if they intersect ArborXWrappers::BVH bounding boxes. + * Constructor. @p points is a list of points used by the predicate. */ template - PointIntersect(const std::vector> &dim_points); + PointPredicate(const std::vector> &points); /** * Number of points stored in the structure. @@ -61,19 +59,66 @@ namespace ArborXWrappers /** - * This class is used by ArborXWrappers::BVH to determine which bounding - * boxes intersect the bounding boxes used to build the ArborXWrappers::BVH. + * This class defines a predicate used by ArborXWrappers::BVH to determine + * for given points which of the bounding boxes used to build the + * ArborXWrappers::BVH intersect with them. */ - class BoundingBoxIntersect + class PointIntersectPredicate : public PointPredicate { public: /** - * Constructor. @p bb is a list of bounding boxes which we are interested in + * Constructor. @p points is a list of points which we are interested in * knowing if they intersect ArborXWrappers::BVH bounding boxes. */ template - BoundingBoxIntersect( - const std::vector> &bb); + PointIntersectPredicate( + const std::vector> &points); + }; + + + + /** + * This class defines a predicate used by ArborXWrappers::BVH to determine + * for given points which are the nearest bounding boxes/points among the ones + * used to build the ArborXWrappers::BVH. + */ + class PointNearestPredicate : public PointPredicate + { + public: + /** + * Constructor. @p points is a list of points for which we are interested in + * the @p n_nearest_neighbors in the ArborXWrappers::BVH bounding + * boxes/points. + */ + template + 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; + + private: + unsigned int n_nearest_neighbors; + }; + + + + /** + * Base class for BoundingBox predicates. + */ + class BoundingBoxPredicate + { + public: + /** + * Constructor. @p bounding_boxes is a list of bounding boxes used by the + * predicate. + */ + template + BoundingBoxPredicate( + const std::vector> &bounding_boxes); /** * Number of bounding boxes stored in the structure. @@ -90,6 +135,54 @@ namespace ArborXWrappers private: std::vector> bounding_boxes; }; + + + + /** + * This class is used by ArborXWrappers::BVH to determine for given bounding + * boxes which of the bounding boxes used to build the ArborXWrappers::BVH + * intersect with them. + */ + class BoundingBoxIntersectPredicate : public BoundingBoxPredicate + { + public: + /** + * Constructor. @p bounding_boxes is a list of bounding boxes which we are interested in + * knowing if they intersect ArborXWrappers::BVH bounding boxes. + */ + template + BoundingBoxIntersectPredicate( + const std::vector> &bounding_boxes); + }; + + + /** + * This class is used by ArborXWrappers::BVH to determine for given bounding + * boxes which are the nearest bounding boxes/points among the ones used to + * build the ArborXWrappers::BVH. + */ + class BoundingBoxNearestPredicate : public BoundingBoxPredicate + { + public: + /** + * Constructor. @p bounding_boxes is a list of bounding boxes for which are interested in + * knowing the @p n_nearest_neighbors nearest bounding boxes used to build the + * ArborXWrappers::BVH. + */ + template + 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; + + private: + unsigned int n_nearest_neighbors; + }; } // namespace ArborXWrappers DEAL_II_NAMESPACE_CLOSE @@ -149,10 +242,11 @@ namespace ArborX /** - * This struct allows ArborX to use PointIntersect in a predicate. + * This struct allows ArborX to use PointIntersectPredicate as a predicate. */ template <> - struct AccessTraits + struct AccessTraits { using memory_space = Kokkos::HostSpace; @@ -160,24 +254,50 @@ namespace ArborX * Number of Point stored in @p pt_intersect. */ static std::size_t - size(const dealii::ArborXWrappers::PointIntersect &pt_intersect); + size(const dealii::ArborXWrappers::PointIntersectPredicate &pt_intersect); /** * Return an Arbox::intersects(ArborX::Point) object constructed from the * `i`th dealii::Point stored in @p pt_intersect. */ static auto - get(const dealii::ArborXWrappers::PointIntersect &pt_intersect, - std::size_t i); + get(const dealii::ArborXWrappers::PointIntersectPredicate &pt_intersect, + std::size_t i); }; + /** + * This struct allows ArborX to use PointNearestPredicate as a predicate. + */ + template <> + struct AccessTraits + { + using memory_space = Kokkos::HostSpace; + + /** + * Number of Point stored in @p pt_nearest. + */ + static std::size_t + size(const dealii::ArborXWrappers::PointNearestPredicate &pt_nearest); + + /** + * Return an Arbox::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 BoundingBoxIntersect in a predicate. + * This struct allows ArborX to use BoundingBoxIntersectPredicate as a + * predicate. */ template <> - struct AccessTraits { using memory_space = Kokkos::HostSpace; @@ -186,15 +306,46 @@ namespace ArborX * Number of BoundingBox stored in @p bb_intersect. */ static std::size_t - size(const dealii::ArborXWrappers::BoundingBoxIntersect &bb_intersect); + 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::BoundingBoxIntersect &bb_intersect, - std::size_t i); + get( + const dealii::ArborXWrappers::BoundingBoxIntersectPredicate &bb_intersect, + std::size_t i); + }; + + + /** + * This struct allows ArborX to use BoundingBoxNearstPredicate as a + * predicate. + */ + template <> + struct AccessTraits + { + using memory_space = Kokkos::HostSpace; + + /** + * Number of BoundingBox 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); }; // ------------------------------- Inline ----------------------------------// @@ -204,8 +355,8 @@ namespace ArborX // We use auto because ArborX does not expose the type of intersects inline std::size_t - AccessTraits::size( - const dealii::ArborXWrappers::PointIntersect &pt_intersect) + AccessTraits:: + size(const dealii::ArborXWrappers::PointIntersectPredicate &pt_intersect) { return pt_intersect.size(); } @@ -213,18 +364,42 @@ namespace ArborX inline auto - AccessTraits::get( - const dealii::ArborXWrappers::PointIntersect &pt_intersect, - std::size_t i) + AccessTraits:: + get(const dealii::ArborXWrappers::PointIntersectPredicate &pt_intersect, + std::size_t i) { const auto dealii_point = pt_intersect.get(i); return intersects(Point{dealii_point[0], dealii_point[1], dealii_point[2]}); } + inline std::size_t - AccessTraits:: - size(const dealii::ArborXWrappers::BoundingBoxIntersect &bb_intersect) + AccessTraits:: + size(const dealii::ArborXWrappers::PointNearestPredicate &pt_nearest) + { + return pt_nearest.size(); + } + + + + inline auto + AccessTraits:: + get(const dealii::ArborXWrappers::PointNearestPredicate &pt_nearest, + std::size_t i) + { + const auto dealii_point = pt_nearest.get(i); + return nearest(Point{dealii_point[0], dealii_point[1], dealii_point[2]}, + pt_nearest.get_n_nearest_neighbors()); + } + + + + inline std::size_t + AccessTraits:: + size( + const dealii::ArborXWrappers::BoundingBoxIntersectPredicate &bb_intersect) { return bb_intersect.size(); } @@ -232,9 +407,11 @@ namespace ArborX inline auto - AccessTraits:: - get(const dealii::ArborXWrappers::BoundingBoxIntersect &bb_intersect, - std::size_t i) + AccessTraits:: + 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<3, float> min_corner = boundary_points.first; @@ -243,6 +420,33 @@ namespace ArborX return intersects(Box{{min_corner[0], min_corner[1], min_corner[2]}, {max_corner[0], max_corner[1], max_corner[2]}}); } + + + + inline std::size_t + AccessTraits:: + size(const dealii::ArborXWrappers::BoundingBoxNearestPredicate &bb_nearest) + { + return bb_nearest.size(); + } + + + + inline auto + AccessTraits:: + 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<3, float> min_corner = boundary_points.first; + const dealii::Point<3, float> max_corner = boundary_points.second; + + 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()); + } } // namespace ArborX #endif diff --git a/include/deal.II/arborx/bvh.h b/include/deal.II/arborx/bvh.h index 98fec0d674..dbd895bb4f 100644 --- a/include/deal.II/arborx/bvh.h +++ b/include/deal.II/arborx/bvh.h @@ -60,6 +60,12 @@ namespace ArborXWrappers template BVH(const std::vector> &bounding_boxes); + /** + * Constructor. Use a vector of @p points as primitives. + */ + template + BVH(const std::vector> &points); + /** * Return the indices of those BoundingBox objects that satisfy the @p queries. * Because @p queries can contain multiple queries, the function returns a pair @@ -82,7 +88,8 @@ namespace ArborXWrappers * * @code * const std::vector> query_bounding_boxes = ... - * ArborXWrappers::BoundingBoxIntersect bb_intersect(query_bounding_boxes); + * ArborXWrappers::BoundingBoxIntersectPredicate + * bb_intersect(query_bounding_boxes); * * const std::vector> bvh_bounding_boxes = ... * ArborxWrappers::BVH bvh(bvh_bounding_boxes); @@ -107,13 +114,26 @@ namespace ArborXWrappers * * @code * const std::vector> query_points = ... - * ArborXWrappers::PointIntersect pt_intersect(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> @@ -135,6 +155,13 @@ namespace ArborXWrappers + template + BVH::BVH(const std::vector> &points) + : bvh(Kokkos::DefaultHostExecutionSpace{}, points) + {} + + + template std::pair, std::vector> BVH::query(const QueryType &queries) diff --git a/source/arborx/access_traits.cc b/source/arborx/access_traits.cc index b6612bd9a8..08fbbc0a60 100644 --- a/source/arborx/access_traits.cc +++ b/source/arborx/access_traits.cc @@ -21,9 +21,9 @@ DEAL_II_NAMESPACE_OPEN namespace ArborXWrappers { - // ------------------- PointIntersect ------------------- // + // ------------------- PointPredicate ------------------- // template - PointIntersect::PointIntersect( + PointPredicate::PointPredicate( const std::vector> &dim_points) { static_assert(dim != 1, "dim equal to one is not supported."); @@ -42,7 +42,7 @@ namespace ArborXWrappers std::size_t - PointIntersect::size() const + PointPredicate::size() const { return points.size(); } @@ -50,16 +50,40 @@ namespace ArborXWrappers const dealii::Point<3, float> & - PointIntersect::get(unsigned int i) const + PointPredicate::get(unsigned int i) const { return points[i]; } - // ------------------- BoundingBoxIntersect ------------------- // template - BoundingBoxIntersect::BoundingBoxIntersect( + PointIntersectPredicate::PointIntersectPredicate( + const std::vector> &points) + : PointPredicate(points) + {} + + + + template + PointNearestPredicate::PointNearestPredicate( + const std::vector> &points, + const unsigned int n_nearest_neighbors) + : PointPredicate(points) + , n_nearest_neighbors(n_nearest_neighbors) + {} + + + + unsigned int + PointNearestPredicate::get_n_nearest_neighbors() const + { + return n_nearest_neighbors; + } + + // ------------------- BoundingBoxPredicate ------------------- // + template + BoundingBoxPredicate::BoundingBoxPredicate( const std::vector> &bb) { const unsigned int size = bb.size(); @@ -84,7 +108,7 @@ namespace ArborXWrappers std::size_t - BoundingBoxIntersect::size() const + BoundingBoxPredicate::size() const { return bounding_boxes.size(); } @@ -92,10 +116,36 @@ namespace ArborXWrappers const dealii::BoundingBox<3, float> & - BoundingBoxIntersect::get(unsigned int i) const + BoundingBoxPredicate::get(unsigned int i) const { return bounding_boxes[i]; } + + + + template + BoundingBoxIntersectPredicate::BoundingBoxIntersectPredicate( + const std::vector> &bounding_boxes) + : BoundingBoxPredicate(bounding_boxes) + {} + + + + template + BoundingBoxNearestPredicate::BoundingBoxNearestPredicate( + const std::vector> &bounding_boxes, + const unsigned int n_nearest_neighbors) + : BoundingBoxPredicate(bounding_boxes) + , n_nearest_neighbors(n_nearest_neighbors) + {} + + + + unsigned int + BoundingBoxNearestPredicate::get_n_nearest_neighbors() const + { + return n_nearest_neighbors; + } } // namespace ArborXWrappers DEAL_II_NAMESPACE_CLOSE diff --git a/source/arborx/access_traits.inst.in b/source/arborx/access_traits.inst.in index c1d094aed3..cc949b0f12 100644 --- a/source/arborx/access_traits.inst.in +++ b/source/arborx/access_traits.inst.in @@ -20,10 +20,17 @@ for (DIM : DIMENSIONS; SCALAR : REAL_SCALARS) DEAL_II_NAMESPACE_OPEN namespace ArborXWrappers { - template PointIntersect::PointIntersect( - const std::vector> &dim_points); - template BoundingBoxIntersect::BoundingBoxIntersect( + template PointIntersectPredicate::PointIntersectPredicate( + const std::vector> &points); + template PointNearestPredicate::PointNearestPredicate( + const std::vector> &points, + const unsigned int n_nearest_neighbors); + + template BoundingBoxIntersectPredicate::BoundingBoxIntersectPredicate( const std::vector> &bb); + template BoundingBoxNearestPredicate::BoundingBoxNearestPredicate( + const std::vector> &bb, + const unsigned int n_nearest_neighbors); \} DEAL_II_NAMESPACE_CLOSE diff --git a/tests/arborx/CMakeLists.txt b/tests/arborx/CMakeLists.txt index 288e920d4c..50b43e346b 100644 --- a/tests/arborx/CMakeLists.txt +++ b/tests/arborx/CMakeLists.txt @@ -1,4 +1,4 @@ -CMAKE_MINIMUM_REQUIRED(VERSION 2.8.9) +CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0) INCLUDE(../setup_testsubproject.cmake) PROJECT(testsuite CXX) IF(DEAL_II_WITH_ARBORX) diff --git a/tests/arborx/bounding_box_intersect.cc b/tests/arborx/bounding_box_intersect.cc index 741d727dba..30aec63c20 100644 --- a/tests/arborx/bounding_box_intersect.cc +++ b/tests/arborx/bounding_box_intersect.cc @@ -14,7 +14,7 @@ // --------------------------------------------------------------------- -// Check ArborX wrapper: intersection of boundary boxes +// Check ArborX wrapper: intersection of bounding boxes #include @@ -60,11 +60,12 @@ 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); - ArborXWrappers::BoundingBoxIntersect bb_intersect(query_bounding_boxes); - auto indices_offset = bvh.query(bb_intersect); - std::vector indices = indices_offset.first; - std::vector offset = indices_offset.second; + ArborXWrappers::BVH bvh(bounding_boxes); + ArborXWrappers::BoundingBoxIntersectPredicate bb_intersect( + query_bounding_boxes); + auto indices_offset = bvh.query(bb_intersect); + std::vector indices = indices_offset.first; + std::vector offset = indices_offset.second; std::vector indices_ref = {0, 1, 4, 5, 10}; std::vector offset_ref = {0, 4, 5}; @@ -137,11 +138,12 @@ 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); - ArborXWrappers::BoundingBoxIntersect bb_intersect(query_bounding_boxes); - auto indices_offset = bvh.query(bb_intersect); - std::vector indices = indices_offset.first; - std::vector offset = indices_offset.second; + ArborXWrappers::BVH bvh(bounding_boxes); + ArborXWrappers::BoundingBoxIntersectPredicate bb_intersect( + query_bounding_boxes); + auto indices_offset = bvh.query(bb_intersect); + std::vector indices = indices_offset.first; + std::vector offset = indices_offset.second; std::vector indices_ref = {0, 1, 4, 5, 16, 17, 20, 21, 42}; std::vector offset_ref = {0, 8, 9}; diff --git a/tests/arborx/bounding_box_nearest.cc b/tests/arborx/bounding_box_nearest.cc new file mode 100644 index 0000000000..e1b01b1730 --- /dev/null +++ b/tests/arborx/bounding_box_nearest.cc @@ -0,0 +1,308 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Check ArborX wrapper: nearest bounding boxes + + +#include +#include + +#include + +#include "../tests.h" + + +void +test_bounding_box_2d() +{ + std::vector> bounding_boxes; + std::vector> points; + + unsigned int n_points_1d = 5; + for (unsigned int i = 0; i < n_points_1d; ++i) + { + for (unsigned int j = 0; j < n_points_1d; ++j) + { + points.emplace_back(j, i); + } + } + + for (unsigned int i = 0; i < n_points_1d - 1; ++i) + { + for (unsigned int j = 0; j < n_points_1d - 1; ++j) + { + unsigned int point_index = j + i * n_points_1d; + bounding_boxes.push_back( + std::make_pair(points[point_index], + points[point_index + n_points_1d + 1])); + } + } + + Point<2> point_a(-1.0, -1.0); + Point<2> point_b(-0.5, -0.5); + Point<2> point_c(5.2, 5.2); + Point<2> point_d(5.6, 5.6); + + std::vector> query_bounding_boxes; + query_bounding_boxes.emplace_back(std::make_pair(point_a, point_b)); + query_bounding_boxes.emplace_back(std::make_pair(point_c, point_d)); + + ArborXWrappers::BVH bvh(bounding_boxes); + ArborXWrappers::BoundingBoxNearestPredicate bb_nearest(query_bounding_boxes, + 1); + auto indices_offset = bvh.query(bb_nearest); + std::vector indices = indices_offset.first; + std::vector offset = indices_offset.second; + + std::vector indices_ref = {0, 15}; + std::vector offset_ref = {0, 1, 2}; + + AssertThrow(indices.size() == indices_ref.size(), ExcInternalError()); + AssertThrow(offset.size() == offset_ref.size(), ExcInternalError()); + for (unsigned int i = 0; i < offset.size() - 1; ++i) + { + for (unsigned int j = offset[i]; j < offset[i + 1]; ++j) + { + // The indices associated to each query are not ordered. + bool found = false; + for (unsigned int k = offset[i]; k < offset[i + 1]; ++k) + { + if (indices[j] == indices_ref[k]) + { + found = true; + break; + } + } + AssertThrow(found, ExcInternalError()); + } + } + for (unsigned int i = 0; i < offset.size(); ++i) + AssertThrow(offset[i] == offset_ref[i], ExcInternalError()); + + deallog << "OK" << std::endl; +} + +void +test_point_2d() +{ + std::vector> points; + + unsigned int n_points_1d = 5; + for (unsigned int i = 0; i < n_points_1d; ++i) + { + for (unsigned int j = 0; j < n_points_1d; ++j) + { + points.emplace_back(j, i); + } + } + + 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; + 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); + ArborXWrappers::BoundingBoxNearestPredicate bb_nearest(query_bounding_boxes, + 1); + auto indices_offset = bvh.query(bb_nearest); + std::vector indices = indices_offset.first; + std::vector offset = indices_offset.second; + + std::vector indices_ref = {6, 12}; + std::vector offset_ref = {0, 1, 2}; + + AssertThrow(indices.size() == indices_ref.size(), ExcInternalError()); + AssertThrow(offset.size() == offset_ref.size(), ExcInternalError()); + for (unsigned int i = 0; i < offset.size() - 1; ++i) + { + for (unsigned int j = offset[i]; j < offset[i + 1]; ++j) + { + // The indices associated to each query are not ordered. + bool found = false; + for (unsigned int k = offset[i]; k < offset[i + 1]; ++k) + { + if (indices[j] == indices_ref[k]) + { + found = true; + break; + } + } + AssertThrow(found, ExcInternalError()); + } + } + for (unsigned int i = 0; i < offset.size(); ++i) + AssertThrow(offset[i] == offset_ref[i], ExcInternalError()); + + deallog << "OK" << std::endl; +} + +void +test_bounding_box_3d() +{ + std::vector> bounding_boxes; + std::vector> points; + + unsigned int n_points_1d = 5; + for (unsigned int i = 0; i < n_points_1d; ++i) + { + for (unsigned int j = 0; j < n_points_1d; ++j) + { + for (unsigned int k = 0; k < n_points_1d; ++k) + { + points.emplace_back(k, j, i); + } + } + } + + for (unsigned int i = 0; i < n_points_1d - 1; ++i) + { + for (unsigned int j = 0; j < n_points_1d - 1; ++j) + { + for (unsigned int k = 0; k < n_points_1d - 1; ++k) + { + unsigned int point_index = + k + j * n_points_1d + i * n_points_1d * n_points_1d; + bounding_boxes.push_back( + std::make_pair(points[point_index], + points[point_index + n_points_1d * n_points_1d + + n_points_1d + 1])); + } + } + } + + Point<3> point_a(-1.0, -1.0, -1.0); + Point<3> point_b(-0.5, -0.5, -0.5); + Point<3> point_c(10.2, 10.2, 10.2); + Point<3> point_d(10.6, 10.6, 10.6); + + std::vector> query_bounding_boxes; + query_bounding_boxes.emplace_back(std::make_pair(point_a, point_b)); + query_bounding_boxes.emplace_back(std::make_pair(point_c, point_d)); + + ArborXWrappers::BVH bvh(bounding_boxes); + ArborXWrappers::BoundingBoxNearestPredicate bb_nearest(query_bounding_boxes, + 1); + auto indices_offset = bvh.query(bb_nearest); + std::vector indices = indices_offset.first; + std::vector offset = indices_offset.second; + + std::vector indices_ref = {0, 63}; + std::vector offset_ref = {0, 1, 2}; + + AssertThrow(indices.size() == indices_ref.size(), ExcInternalError()); + AssertThrow(offset.size() == offset_ref.size(), ExcInternalError()); + for (unsigned int i = 0; i < offset.size() - 1; ++i) + { + for (unsigned int j = offset[i]; j < offset[i + 1]; ++j) + { + // The indices associated to each query are not ordered. + bool found = false; + for (unsigned int k = offset[i]; k < offset[i + 1]; ++k) + { + if (indices[j] == indices_ref[k]) + { + found = true; + break; + } + } + AssertThrow(found, ExcInternalError()); + } + } + for (unsigned int i = 0; i < offset.size(); ++i) + AssertThrow(offset[i] == offset_ref[i], ExcInternalError()); + deallog << "OK" << std::endl; +} + +void +test_point_3d() +{ + std::vector> points; + + unsigned int n_points_1d = 5; + for (unsigned int i = 0; i < n_points_1d; ++i) + { + for (unsigned int j = 0; j < n_points_1d; ++j) + { + for (unsigned int k = 0; k < n_points_1d; ++k) + { + points.emplace_back(k, j, i); + } + } + } + + 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; + 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); + ArborXWrappers::BoundingBoxNearestPredicate bb_nearest(query_bounding_boxes, + 1); + auto indices_offset = bvh.query(bb_nearest); + std::vector indices = indices_offset.first; + std::vector offset = indices_offset.second; + + std::vector indices_ref = {31, 62}; + std::vector offset_ref = {0, 1, 2}; + + AssertThrow(indices.size() == indices_ref.size(), ExcInternalError()); + AssertThrow(offset.size() == offset_ref.size(), ExcInternalError()); + for (unsigned int i = 0; i < offset.size() - 1; ++i) + { + for (unsigned int j = offset[i]; j < offset[i + 1]; ++j) + { + // The indices associated to each query are not ordered. + bool found = false; + for (unsigned int k = offset[i]; k < offset[i + 1]; ++k) + { + if (indices[j] == indices_ref[k]) + { + found = true; + break; + } + } + AssertThrow(found, ExcInternalError()); + } + } + for (unsigned int i = 0; i < offset.size(); ++i) + AssertThrow(offset[i] == offset_ref[i], ExcInternalError()); + deallog << "OK" << std::endl; +} + +int +main(int argc, char **argv) +{ + // Initialize Kokkos + Kokkos::initialize(argc, argv); + + initlog(); + + // tests + test_bounding_box_2d(); + test_bounding_box_3d(); + test_point_2d(); + test_point_3d(); + + Kokkos::finalize(); +} diff --git a/tests/arborx/bounding_box_nearest.output b/tests/arborx/bounding_box_nearest.output new file mode 100644 index 0000000000..5f42bb230f --- /dev/null +++ b/tests/arborx/bounding_box_nearest.output @@ -0,0 +1,5 @@ + +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK diff --git a/tests/arborx/point_intersect.cc b/tests/arborx/point_intersect.cc index 6093454953..ca972f7bb6 100644 --- a/tests/arborx/point_intersect.cc +++ b/tests/arborx/point_intersect.cc @@ -58,11 +58,11 @@ test_2d() query_points.emplace_back(2.6, 2.6); - ArborXWrappers::BVH bvh(bounding_boxes); - ArborXWrappers::PointIntersect pt_intersect(query_points); - auto indices_offset = bvh.query(pt_intersect); - std::vector indices = indices_offset.first; - std::vector offset = indices_offset.second; + ArborXWrappers::BVH bvh(bounding_boxes); + ArborXWrappers::PointIntersectPredicate pt_intersect(query_points); + auto indices_offset = bvh.query(pt_intersect); + std::vector indices = indices_offset.first; + std::vector offset = indices_offset.second; std::vector indices_ref = {0, 5, 10, 10}; std::vector offset_ref = {0, 1, 2, 3, 4}; @@ -133,11 +133,11 @@ test_3d() query_points.emplace_back(2.6, 2.6, 2.6); - ArborXWrappers::BVH bvh(bounding_boxes); - ArborXWrappers::PointIntersect pt_intersect(query_points); - auto indices_offset = bvh.query(pt_intersect); - std::vector indices = indices_offset.first; - std::vector offset = indices_offset.second; + ArborXWrappers::BVH bvh(bounding_boxes); + ArborXWrappers::PointIntersectPredicate pt_intersect(query_points); + auto indices_offset = bvh.query(pt_intersect); + std::vector indices = indices_offset.first; + std::vector offset = indices_offset.second; std::vector indices_ref = {0, 21, 42, 42}; diff --git a/tests/arborx/point_nearest.cc b/tests/arborx/point_nearest.cc new file mode 100644 index 0000000000..1e9af5ee6e --- /dev/null +++ b/tests/arborx/point_nearest.cc @@ -0,0 +1,301 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Check ArborX wrapper: nearest points from bounding boxes and others points + + +#include +#include + +#include + +#include "../tests.h" + + +void +test_bounding_box_2d() +{ + std::vector> bounding_boxes; + std::vector> points; + + unsigned int n_points_1d = 5; + for (unsigned int i = 0; i < n_points_1d; ++i) + { + for (unsigned int j = 0; j < n_points_1d; ++j) + { + points.emplace_back(j, i); + } + } + + for (unsigned int i = 0; i < n_points_1d - 1; ++i) + { + for (unsigned int j = 0; j < n_points_1d - 1; ++j) + { + unsigned int point_index = j + i * n_points_1d; + bounding_boxes.push_back( + std::make_pair(points[point_index], + points[point_index + n_points_1d + 1])); + } + } + + std::vector> query_points; + query_points.emplace_back(0.5, 0.5); + query_points.emplace_back(1.5, 1.5); + query_points.emplace_back(2.2, 2.2); + query_points.emplace_back(2.6, 2.6); + + + ArborXWrappers::BVH bvh(bounding_boxes); + ArborXWrappers::PointNearestPredicate pt_nearest(query_points, 1); + auto indices_offset = bvh.query(pt_nearest); + std::vector indices = indices_offset.first; + std::vector offset = indices_offset.second; + + std::vector indices_ref = {0, 5, 10, 10}; + std::vector offset_ref = {0, 1, 2, 3, 4}; + + AssertThrow(indices.size() == indices_ref.size(), ExcInternalError()); + AssertThrow(offset.size() == offset_ref.size(), ExcInternalError()); + for (unsigned int i = 0; i < offset.size() - 1; ++i) + { + for (unsigned int j = offset[i]; j < offset[i + 1]; ++j) + { + // The indices associated to each query are not ordered. + bool found = false; + for (unsigned int k = offset[i]; k < offset[i + 1]; ++k) + { + if (indices[j] == indices_ref[k]) + { + found = true; + break; + } + } + AssertThrow(found, ExcInternalError()); + } + } + for (unsigned int i = 0; i < offset.size(); ++i) + AssertThrow(offset[i] == offset_ref[i], ExcInternalError()); + + deallog << "OK" << std::endl; +} + +void +test_points_2d() +{ + std::vector> points; + + unsigned int n_points_1d = 5; + for (unsigned int i = 0; i < n_points_1d; ++i) + { + for (unsigned int j = 0; j < n_points_1d; ++j) + { + points.emplace_back(j, i); + } + } + + + std::vector> query_points; + query_points.emplace_back(0.5, 0.5); + query_points.emplace_back(1.5, 1.5); + query_points.emplace_back(2.2, 2.2); + query_points.emplace_back(2.6, 2.6); + + + ArborXWrappers::BVH bvh(points); + ArborXWrappers::PointNearestPredicate pt_nearest(query_points, 1); + auto indices_offset = bvh.query(pt_nearest); + std::vector indices = indices_offset.first; + std::vector offset = indices_offset.second; + + std::vector indices_ref = {0, 6, 12, 18}; + std::vector offset_ref = {0, 1, 2, 3, 4}; + + AssertThrow(indices.size() == indices_ref.size(), ExcInternalError()); + AssertThrow(offset.size() == offset_ref.size(), ExcInternalError()); + for (unsigned int i = 0; i < offset.size() - 1; ++i) + { + for (unsigned int j = offset[i]; j < offset[i + 1]; ++j) + { + // The indices associated to each query are not ordered. + bool found = false; + for (unsigned int k = offset[i]; k < offset[i + 1]; ++k) + { + if (indices[j] == indices_ref[k]) + { + found = true; + break; + } + } + AssertThrow(found, ExcInternalError()); + } + } + for (unsigned int i = 0; i < offset.size(); ++i) + AssertThrow(offset[i] == offset_ref[i], ExcInternalError()); + + deallog << "OK" << std::endl; +} + +void +test_bounding_box_3d() +{ + std::vector> bounding_boxes; + std::vector> points; + + unsigned int n_points_1d = 5; + for (unsigned int i = 0; i < n_points_1d; ++i) + { + for (unsigned int j = 0; j < n_points_1d; ++j) + { + for (unsigned int k = 0; k < n_points_1d; ++k) + { + points.emplace_back(k, j, i); + } + } + } + + for (unsigned int i = 0; i < n_points_1d - 1; ++i) + { + for (unsigned int j = 0; j < n_points_1d - 1; ++j) + { + for (unsigned int k = 0; k < n_points_1d - 1; ++k) + { + unsigned int point_index = + k + j * n_points_1d + i * n_points_1d * n_points_1d; + bounding_boxes.push_back( + std::make_pair(points[point_index], + points[point_index + n_points_1d * n_points_1d + + n_points_1d + 1])); + } + } + } + + std::vector> query_points; + query_points.emplace_back(0.5, 0.5, 0.5); + query_points.emplace_back(1.5, 1.5, 1.5); + query_points.emplace_back(2.2, 2.2, 2.2); + query_points.emplace_back(2.6, 2.6, 2.6); + + + ArborXWrappers::BVH bvh(bounding_boxes); + ArborXWrappers::PointNearestPredicate pt_nearest(query_points, 1); + auto indices_offset = bvh.query(pt_nearest); + std::vector indices = indices_offset.first; + std::vector offset = indices_offset.second; + + + std::vector indices_ref = {0, 21, 42, 42}; + std::vector offset_ref = {0, 1, 2, 3, 4}; + + AssertThrow(indices.size() == indices_ref.size(), ExcInternalError()); + AssertThrow(offset.size() == offset_ref.size(), ExcInternalError()); + for (unsigned int i = 0; i < offset.size() - 1; ++i) + { + for (unsigned int j = offset[i]; j < offset[i + 1]; ++j) + { + // The indices associated to each query are not ordered. + bool found = false; + for (unsigned int k = offset[i]; k < offset[i + 1]; ++k) + { + if (indices[j] == indices_ref[k]) + { + found = true; + break; + } + } + AssertThrow(found, ExcInternalError()); + } + } + for (unsigned int i = 0; i < offset.size(); ++i) + AssertThrow(offset[i] == offset_ref[i], ExcInternalError()); + deallog << "OK" << std::endl; +} + +void +test_points_3d() +{ + std::vector> points; + + unsigned int n_points_1d = 5; + for (unsigned int i = 0; i < n_points_1d; ++i) + { + for (unsigned int j = 0; j < n_points_1d; ++j) + { + for (unsigned int k = 0; k < n_points_1d; ++k) + { + points.emplace_back(k, j, i); + } + } + } + + std::vector> query_points; + query_points.emplace_back(0.5, 0.5, 0.5); + query_points.emplace_back(1.5, 1.5, 1.5); + query_points.emplace_back(2.2, 2.2, 2.2); + query_points.emplace_back(2.6, 2.6, 2.6); + + + ArborXWrappers::BVH bvh(points); + ArborXWrappers::PointNearestPredicate pt_nearest(query_points, 1); + auto indices_offset = bvh.query(pt_nearest); + std::vector indices = indices_offset.first; + std::vector offset = indices_offset.second; + + std::vector indices_ref = {0, 31, 62, 93}; + std::vector offset_ref = {0, 1, 2, 3, 4}; + + AssertThrow(indices.size() == indices_ref.size(), ExcInternalError()); + AssertThrow(offset.size() == offset_ref.size(), ExcInternalError()); + for (unsigned int i = 0; i < offset.size() - 1; ++i) + { + for (unsigned int j = offset[i]; j < offset[i + 1]; ++j) + { + // The indices associated to each query are not ordered. + bool found = false; + for (unsigned int k = offset[i]; k < offset[i + 1]; ++k) + { + if (indices[j] == indices_ref[k]) + { + found = true; + break; + } + } + AssertThrow(found, ExcInternalError()); + } + } + for (unsigned int i = 0; i < offset.size(); ++i) + AssertThrow(offset[i] == offset_ref[i], ExcInternalError()); + deallog << "OK" << std::endl; +} + +int +main(int argc, char **argv) +{ + // Initialize Kokkos + Kokkos::initialize(argc, argv); + + initlog(); + + // Initialize ArborX + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); + + // tests + test_bounding_box_2d(); + test_bounding_box_3d(); + test_points_2d(); + test_points_3d(); + + Kokkos::finalize(); +} diff --git a/tests/arborx/point_nearest.output b/tests/arborx/point_nearest.output new file mode 100644 index 0000000000..5f42bb230f --- /dev/null +++ b/tests/arborx/point_nearest.output @@ -0,0 +1,5 @@ + +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK -- 2.39.5