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 <int dim, typename Number>
- PointIntersect(const std::vector<dealii::Point<dim, Number>> &dim_points);
+ PointPredicate(const std::vector<dealii::Point<dim, Number>> &points);
/**
* Number of points stored in the structure.
/**
- * 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 <int dim, typename Number>
- BoundingBoxIntersect(
- const std::vector<dealii::BoundingBox<dim, Number>> &bb);
+ PointIntersectPredicate(
+ const 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 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 <int dim, typename Number>
+ 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;
+
+ 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 <int dim, typename Number>
+ BoundingBoxPredicate(
+ const std::vector<dealii::BoundingBox<dim, Number>> &bounding_boxes);
/**
* Number of bounding boxes stored in the structure.
private:
std::vector<dealii::BoundingBox<3, float>> 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 <int dim, typename Number>
+ BoundingBoxIntersectPredicate(
+ const 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 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 <int dim, typename Number>
+ 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;
+
+ private:
+ unsigned int n_nearest_neighbors;
+ };
} // namespace ArborXWrappers
DEAL_II_NAMESPACE_CLOSE
/**
- * This struct allows ArborX to use PointIntersect in a predicate.
+ * This struct allows ArborX to use PointIntersectPredicate as a predicate.
*/
template <>
- struct AccessTraits<dealii::ArborXWrappers::PointIntersect, PredicatesTag>
+ struct AccessTraits<dealii::ArborXWrappers::PointIntersectPredicate,
+ PredicatesTag>
{
using memory_space = Kokkos::HostSpace;
* 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<dealii::ArborXWrappers::PointNearestPredicate,
+ PredicatesTag>
+ {
+ 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<dealii::ArborXWrappers::BoundingBoxIntersect,
+ struct AccessTraits<dealii::ArborXWrappers::BoundingBoxIntersectPredicate,
PredicatesTag>
{
using memory_space = Kokkos::HostSpace;
* 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<dealii::ArborXWrappers::BoundingBoxNearestPredicate,
+ PredicatesTag>
+ {
+ 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 ----------------------------------//
// We use auto because ArborX does not expose the type of intersects
inline std::size_t
- AccessTraits<dealii::ArborXWrappers::PointIntersect, PredicatesTag>::size(
- const dealii::ArborXWrappers::PointIntersect &pt_intersect)
+ AccessTraits<dealii::ArborXWrappers::PointIntersectPredicate, PredicatesTag>::
+ size(const dealii::ArborXWrappers::PointIntersectPredicate &pt_intersect)
{
return pt_intersect.size();
}
inline auto
- AccessTraits<dealii::ArborXWrappers::PointIntersect, PredicatesTag>::get(
- const dealii::ArborXWrappers::PointIntersect &pt_intersect,
- std::size_t i)
+ AccessTraits<dealii::ArborXWrappers::PointIntersectPredicate, PredicatesTag>::
+ 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<dealii::ArborXWrappers::BoundingBoxIntersect, PredicatesTag>::
- size(const dealii::ArborXWrappers::BoundingBoxIntersect &bb_intersect)
+ AccessTraits<dealii::ArborXWrappers::PointNearestPredicate, PredicatesTag>::
+ size(const dealii::ArborXWrappers::PointNearestPredicate &pt_nearest)
+ {
+ return pt_nearest.size();
+ }
+
+
+
+ inline auto
+ AccessTraits<dealii::ArborXWrappers::PointNearestPredicate, PredicatesTag>::
+ 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<dealii::ArborXWrappers::BoundingBoxIntersectPredicate,
+ PredicatesTag>::
+ size(
+ const dealii::ArborXWrappers::BoundingBoxIntersectPredicate &bb_intersect)
{
return bb_intersect.size();
}
inline auto
- AccessTraits<dealii::ArborXWrappers::BoundingBoxIntersect, PredicatesTag>::
- get(const dealii::ArborXWrappers::BoundingBoxIntersect &bb_intersect,
- std::size_t i)
+ AccessTraits<dealii::ArborXWrappers::BoundingBoxIntersectPredicate,
+ PredicatesTag>::
+ 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;
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<dealii::ArborXWrappers::BoundingBoxNearestPredicate,
+ PredicatesTag>::
+ size(const dealii::ArborXWrappers::BoundingBoxNearestPredicate &bb_nearest)
+ {
+ return bb_nearest.size();
+ }
+
+
+
+ inline auto
+ AccessTraits<dealii::ArborXWrappers::BoundingBoxNearestPredicate,
+ PredicatesTag>::
+ 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
template <int dim, typename Number>
BVH(const std::vector<BoundingBox<dim, Number>> &bounding_boxes);
+ /**
+ * Constructor. Use a vector of @p points as primitives.
+ */
+ template <int dim, typename Number>
+ BVH(const std::vector<Point<dim, Number>> &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
*
* @code
* const std::vector<BoundingBox<dim>> query_bounding_boxes = ...
- * ArborXWrappers::BoundingBoxIntersect bb_intersect(query_bounding_boxes);
+ * ArborXWrappers::BoundingBoxIntersectPredicate
+ * bb_intersect(query_bounding_boxes);
*
* const std::vector<BoundingBox<dim>> bvh_bounding_boxes = ...
* ArborxWrappers::BVH bvh(bvh_bounding_boxes);
*
* @code
* const std::vector<Point<dim>> query_points = ...
- * ArborXWrappers::PointIntersect pt_intersect(query_points);
+ * ArborXWrappers::PointIntersectPredicate pt_intersect(query_points);
*
* const std::vector<BoundingBox<dim>> 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<Point<dim>> query_points = ...
+ * ArborXWrappers::PointNearestPredicate pt_nearest(query_points, 5);
+ *
+ * const std::vector<Point<dim>> bvh_points = ...
+ * ArborxWrappers::BVH bvh(bvh_points);
+ *
+ * auto [indices, offset] = bvh.query(pt_nearest);
+ * @endcode
*/
template <typename QueryType>
std::pair<std::vector<int>, std::vector<int>>
+ template <int dim, typename Number>
+ BVH::BVH(const std::vector<Point<dim, Number>> &points)
+ : bvh(Kokkos::DefaultHostExecutionSpace{}, points)
+ {}
+
+
+
template <typename QueryType>
std::pair<std::vector<int>, std::vector<int>>
BVH::query(const QueryType &queries)
namespace ArborXWrappers
{
- // ------------------- PointIntersect ------------------- //
+ // ------------------- PointPredicate ------------------- //
template <int dim, typename Number>
- PointIntersect::PointIntersect(
+ PointPredicate::PointPredicate(
const std::vector<dealii::Point<dim, Number>> &dim_points)
{
static_assert(dim != 1, "dim equal to one is not supported.");
std::size_t
- PointIntersect::size() const
+ PointPredicate::size() const
{
return points.size();
}
const dealii::Point<3, float> &
- PointIntersect::get(unsigned int i) const
+ PointPredicate::get(unsigned int i) const
{
return points[i];
}
- // ------------------- BoundingBoxIntersect ------------------- //
template <int dim, typename Number>
- BoundingBoxIntersect::BoundingBoxIntersect(
+ PointIntersectPredicate::PointIntersectPredicate(
+ const std::vector<dealii::Point<dim, Number>> &points)
+ : PointPredicate(points)
+ {}
+
+
+
+ template <int dim, typename Number>
+ PointNearestPredicate::PointNearestPredicate(
+ const std::vector<dealii::Point<dim, Number>> &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 <int dim, typename Number>
+ BoundingBoxPredicate::BoundingBoxPredicate(
const std::vector<dealii::BoundingBox<dim, Number>> &bb)
{
const unsigned int size = bb.size();
std::size_t
- BoundingBoxIntersect::size() const
+ BoundingBoxPredicate::size() const
{
return bounding_boxes.size();
}
const dealii::BoundingBox<3, float> &
- BoundingBoxIntersect::get(unsigned int i) const
+ BoundingBoxPredicate::get(unsigned int i) const
{
return bounding_boxes[i];
}
+
+
+
+ template <int dim, typename Number>
+ BoundingBoxIntersectPredicate::BoundingBoxIntersectPredicate(
+ const std::vector<dealii::BoundingBox<dim, Number>> &bounding_boxes)
+ : BoundingBoxPredicate(bounding_boxes)
+ {}
+
+
+
+ template <int dim, typename Number>
+ BoundingBoxNearestPredicate::BoundingBoxNearestPredicate(
+ const std::vector<dealii::BoundingBox<dim, Number>> &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
DEAL_II_NAMESPACE_OPEN
namespace ArborXWrappers
{
- template PointIntersect::PointIntersect(
- const std::vector<dealii::Point<DIM, SCALAR>> &dim_points);
- template BoundingBoxIntersect::BoundingBoxIntersect(
+ template PointIntersectPredicate::PointIntersectPredicate(
+ const std::vector<dealii::Point<DIM, SCALAR>> &points);
+ template PointNearestPredicate::PointNearestPredicate(
+ const std::vector<dealii::Point<DIM, SCALAR>> &points,
+ const unsigned int n_nearest_neighbors);
+
+ template BoundingBoxIntersectPredicate::BoundingBoxIntersectPredicate(
const std::vector<dealii::BoundingBox<DIM, SCALAR>> &bb);
+ template BoundingBoxNearestPredicate::BoundingBoxNearestPredicate(
+ const std::vector<dealii::BoundingBox<DIM, SCALAR>> &bb,
+ const unsigned int n_nearest_neighbors);
\}
DEAL_II_NAMESPACE_CLOSE
-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)
// ---------------------------------------------------------------------
-// Check ArborX wrapper: intersection of boundary boxes
+// Check ArborX wrapper: intersection of bounding boxes
#include <deal.II/arborx/access_traits.h>
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<int> indices = indices_offset.first;
- std::vector<int> 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<int> indices = indices_offset.first;
+ std::vector<int> offset = indices_offset.second;
std::vector<int> indices_ref = {0, 1, 4, 5, 10};
std::vector<int> offset_ref = {0, 4, 5};
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<int> indices = indices_offset.first;
- std::vector<int> 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<int> indices = indices_offset.first;
+ std::vector<int> offset = indices_offset.second;
std::vector<int> indices_ref = {0, 1, 4, 5, 16, 17, 20, 21, 42};
std::vector<int> offset_ref = {0, 8, 9};
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/arborx/access_traits.h>
+#include <deal.II/arborx/bvh.h>
+
+#include <deal.II/base/bounding_box.h>
+
+#include "../tests.h"
+
+
+void
+test_bounding_box_2d()
+{
+ std::vector<BoundingBox<2>> bounding_boxes;
+ std::vector<Point<2>> 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<BoundingBox<2>> 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<int> indices = indices_offset.first;
+ std::vector<int> offset = indices_offset.second;
+
+ std::vector<int> indices_ref = {0, 15};
+ std::vector<int> 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<Point<2>> 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<BoundingBox<3>> 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<int> indices = indices_offset.first;
+ std::vector<int> offset = indices_offset.second;
+
+ std::vector<int> indices_ref = {6, 12};
+ std::vector<int> 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<BoundingBox<3>> bounding_boxes;
+ std::vector<Point<3>> 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<BoundingBox<3>> 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<int> indices = indices_offset.first;
+ std::vector<int> offset = indices_offset.second;
+
+ std::vector<int> indices_ref = {0, 63};
+ std::vector<int> 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<Point<3>> 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<BoundingBox<3>> 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<int> indices = indices_offset.first;
+ std::vector<int> offset = indices_offset.second;
+
+ std::vector<int> indices_ref = {31, 62};
+ std::vector<int> 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();
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
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<int> indices = indices_offset.first;
- std::vector<int> offset = indices_offset.second;
+ ArborXWrappers::BVH bvh(bounding_boxes);
+ ArborXWrappers::PointIntersectPredicate pt_intersect(query_points);
+ auto indices_offset = bvh.query(pt_intersect);
+ std::vector<int> indices = indices_offset.first;
+ std::vector<int> offset = indices_offset.second;
std::vector<int> indices_ref = {0, 5, 10, 10};
std::vector<int> offset_ref = {0, 1, 2, 3, 4};
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<int> indices = indices_offset.first;
- std::vector<int> offset = indices_offset.second;
+ ArborXWrappers::BVH bvh(bounding_boxes);
+ ArborXWrappers::PointIntersectPredicate pt_intersect(query_points);
+ auto indices_offset = bvh.query(pt_intersect);
+ std::vector<int> indices = indices_offset.first;
+ std::vector<int> offset = indices_offset.second;
std::vector<int> indices_ref = {0, 21, 42, 42};
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/arborx/access_traits.h>
+#include <deal.II/arborx/bvh.h>
+
+#include <deal.II/base/point.h>
+
+#include "../tests.h"
+
+
+void
+test_bounding_box_2d()
+{
+ std::vector<BoundingBox<2>> bounding_boxes;
+ std::vector<Point<2>> 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<Point<2>> 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<int> indices = indices_offset.first;
+ std::vector<int> offset = indices_offset.second;
+
+ std::vector<int> indices_ref = {0, 5, 10, 10};
+ std::vector<int> offset_ref = {0, 1, 2, 3, 4};
+
+ AssertThrow(indices.size() == indices_ref.size(), ExcInternalError());
+ AssertThrow(offset.size() == offset_ref.size(), ExcInternalError());
+ for (unsigned int i = 0; i < offset.size() - 1; ++i)
+ {
+ for (unsigned int j = offset[i]; j < offset[i + 1]; ++j)
+ {
+ // The indices associated to each query are not ordered.
+ bool found = false;
+ for (unsigned int k = offset[i]; k < offset[i + 1]; ++k)
+ {
+ if (indices[j] == indices_ref[k])
+ {
+ found = true;
+ break;
+ }
+ }
+ AssertThrow(found, ExcInternalError());
+ }
+ }
+ for (unsigned int i = 0; i < offset.size(); ++i)
+ AssertThrow(offset[i] == offset_ref[i], ExcInternalError());
+
+ deallog << "OK" << std::endl;
+}
+
+void
+test_points_2d()
+{
+ std::vector<Point<2>> 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<Point<2>> 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<int> indices = indices_offset.first;
+ std::vector<int> offset = indices_offset.second;
+
+ std::vector<int> indices_ref = {0, 6, 12, 18};
+ std::vector<int> offset_ref = {0, 1, 2, 3, 4};
+
+ AssertThrow(indices.size() == indices_ref.size(), ExcInternalError());
+ AssertThrow(offset.size() == offset_ref.size(), ExcInternalError());
+ for (unsigned int i = 0; i < offset.size() - 1; ++i)
+ {
+ for (unsigned int j = offset[i]; j < offset[i + 1]; ++j)
+ {
+ // The indices associated to each query are not ordered.
+ bool found = false;
+ for (unsigned int k = offset[i]; k < offset[i + 1]; ++k)
+ {
+ if (indices[j] == indices_ref[k])
+ {
+ found = true;
+ break;
+ }
+ }
+ AssertThrow(found, ExcInternalError());
+ }
+ }
+ for (unsigned int i = 0; i < offset.size(); ++i)
+ AssertThrow(offset[i] == offset_ref[i], ExcInternalError());
+
+ deallog << "OK" << std::endl;
+}
+
+void
+test_bounding_box_3d()
+{
+ std::vector<BoundingBox<3>> bounding_boxes;
+ std::vector<Point<3>> 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<Point<3>> 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<int> indices = indices_offset.first;
+ std::vector<int> offset = indices_offset.second;
+
+
+ std::vector<int> indices_ref = {0, 21, 42, 42};
+ std::vector<int> offset_ref = {0, 1, 2, 3, 4};
+
+ AssertThrow(indices.size() == indices_ref.size(), ExcInternalError());
+ AssertThrow(offset.size() == offset_ref.size(), ExcInternalError());
+ for (unsigned int i = 0; i < offset.size() - 1; ++i)
+ {
+ for (unsigned int j = offset[i]; j < offset[i + 1]; ++j)
+ {
+ // The indices associated to each query are not ordered.
+ bool found = false;
+ for (unsigned int k = offset[i]; k < offset[i + 1]; ++k)
+ {
+ if (indices[j] == indices_ref[k])
+ {
+ found = true;
+ break;
+ }
+ }
+ AssertThrow(found, ExcInternalError());
+ }
+ }
+ for (unsigned int i = 0; i < offset.size(); ++i)
+ AssertThrow(offset[i] == offset_ref[i], ExcInternalError());
+ deallog << "OK" << std::endl;
+}
+
+void
+test_points_3d()
+{
+ std::vector<Point<3>> 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<Point<3>> 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<int> indices = indices_offset.first;
+ std::vector<int> offset = indices_offset.second;
+
+ std::vector<int> indices_ref = {0, 31, 62, 93};
+ std::vector<int> offset_ref = {0, 1, 2, 3, 4};
+
+ AssertThrow(indices.size() == indices_ref.size(), ExcInternalError());
+ AssertThrow(offset.size() == offset_ref.size(), ExcInternalError());
+ for (unsigned int i = 0; i < offset.size() - 1; ++i)
+ {
+ for (unsigned int j = offset[i]; j < offset[i + 1]; ++j)
+ {
+ // The indices associated to each query are not ordered.
+ bool found = false;
+ for (unsigned int k = offset[i]; k < offset[i + 1]; ++k)
+ {
+ if (indices[j] == indices_ref[k])
+ {
+ found = true;
+ break;
+ }
+ }
+ AssertThrow(found, ExcInternalError());
+ }
+ }
+ for (unsigned int i = 0; i < offset.size(); ++i)
+ AssertThrow(offset[i] == offset_ref[i], ExcInternalError());
+ deallog << "OK" << std::endl;
+}
+
+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();
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
projects / dealii.git / commitdiff