--- /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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_arborx_access_traits_h
+#define dealii_arborx_access_traits_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_ARBORX
+# include <deal.II/base/bounding_box.h>
+
+# include <ArborX.hpp>
+
+
+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.
+ */
+ class PointIntersect
+ {
+ public:
+ /**
+ * Constructor. @p dim_points is a list of points which we are interested in
+ * knowing if they intersect ArborXWrappers::BVH bounding boxes.
+ */
+ template <int dim, typename Number>
+ PointIntersect(const std::vector<dealii::Point<dim, Number>> &dim_points);
+
+ /**
+ * Number of points stored in the structure.
+ */
+ std::size_t
+ size() const;
+
+ /**
+ * Return the `i`th Point stored in the object.
+ */
+ const dealii::Point<3, float> &
+ get(unsigned int i) const;
+
+ private:
+ std::vector<dealii::Point<3, float>> points;
+ };
+
+
+
+ /**
+ * This class is used by ArborXWrappers::BVH to determine which bounding
+ * boxes intersect the bounding boxes used to build the ArborXWrappers::BVH.
+ */
+ class BoundingBoxIntersect
+ {
+ public:
+ /**
+ * Constructor. @p bb is a list of bounding boxes 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);
+
+ /**
+ * Number of bounding boxes stored in the structure.
+ */
+ std::size_t
+ size() const;
+
+ /**
+ * Return the `i`th BoundingBox stored in the object.
+ */
+ const dealii::BoundingBox<3, float> &
+ get(unsigned int i) const;
+
+ private:
+ std::vector<dealii::BoundingBox<3, float>> bounding_boxes;
+ };
+} // namespace ArborXWrappers
+
+DEAL_II_NAMESPACE_CLOSE
+
+/**
+ * This namespace contains the implementation of AccessTraits used by ArborX.
+ */
+namespace ArborX
+{
+ /**
+ * This struct allows ArborX to use std::vector<dealii::Point> as
+ * primitive.
+ */
+ template <int dim, typename Number>
+ struct AccessTraits<std::vector<dealii::Point<dim, Number>>, PrimitivesTag>
+ {
+ using memory_space = Kokkos::HostSpace;
+
+ /**
+ * Return the size of the vector @p v.
+ */
+ static std::size_t
+ size(const std::vector<dealii::Point<dim, Number>> &v);
+
+ /**
+ * Return an ArborX::Point from the dealii::Point `v[i]`.
+ */
+ static Point
+ get(const std::vector<dealii::Point<dim, Number>> &v, std::size_t i);
+ };
+
+
+
+ /**
+ * This struct allows ArborX to use std::vector<dealii::BoundingBox> as
+ * primitive.
+ */
+ template <int dim, typename Number>
+ struct AccessTraits<std::vector<dealii::BoundingBox<dim, Number>>,
+ PrimitivesTag>
+ {
+ using memory_space = Kokkos::HostSpace;
+
+ /**
+ * Return the size of the vector @p v.
+ */
+ static std::size_t
+ size(const std::vector<dealii::BoundingBox<dim, Number>> &v);
+
+ /**
+ * Return an ArborX::Box from the dealii::BoundingBox `v[i]`.
+ */
+ static Box
+ get(const std::vector<dealii::BoundingBox<dim, Number>> &v, std::size_t i);
+ };
+
+
+
+ /**
+ * This struct allows ArborX to use PointIntersect in a predicate.
+ */
+ template <>
+ struct AccessTraits<dealii::ArborXWrappers::PointIntersect, 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);
+
+ /**
+ * 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);
+ };
+
+
+
+ /**
+ * This struct allows ArborX to use BoundingBoxIntersect in a predicate.
+ */
+ template <>
+ struct AccessTraits<dealii::ArborXWrappers::BoundingBoxIntersect,
+ 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);
+
+ /**
+ * 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);
+ };
+
+ // ------------------------------- Inline ----------------------------------//
+
+ // The implementation of AccessTraits<..., PredicatesTag> needs to be in the
+ // header file otherwise the return type of auto get() cannot be determined.
+ // 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)
+ {
+ return pt_intersect.size();
+ }
+
+
+
+ inline auto
+ AccessTraits<dealii::ArborXWrappers::PointIntersect, PredicatesTag>::get(
+ const dealii::ArborXWrappers::PointIntersect &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)
+ {
+ return bb_intersect.size();
+ }
+
+
+
+ inline auto
+ AccessTraits<dealii::ArborXWrappers::BoundingBoxIntersect, PredicatesTag>::
+ get(const dealii::ArborXWrappers::BoundingBoxIntersect &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;
+ const dealii::Point<3, float> max_corner = boundary_points.second;
+
+ return intersects(Box{{min_corner[0], min_corner[1], min_corner[2]},
+ {max_corner[0], max_corner[1], max_corner[2]}});
+ }
+} // namespace ArborX
+
+#endif
+
+#endif
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_arborx_bvh_h
+#define dealii_arborx_bvh_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_ARBORX
+# include <deal.II/arborx/access_traits.h>
+
+# include <ArborX_LinearBVH.hpp>
+# include <Kokkos_Core.hpp>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * This namespace contains wrappers for the ArborX library.
+ */
+namespace ArborXWrappers
+{
+ /**
+ * This class implements a wrapper around ArborX::BVH. BVH stands for Bounding
+ * Volume Hierarchy.
+ *
+ * From [Wikipedia](https://en.wikipedia.org/wiki/Bounding_Volume_Hierarchy):
+ * <blockquote>
+ * A bounding volume hierarchy (BVH) is a tree structure on a set of geometric
+ * objects. All geometric objects are wrapped in bounding volumes that form
+ * the leaf nodes of the tree. These nodes are then grouped as small sets and
+ * enclosed within larger bounding volumes. These, in turn, are also grouped
+ * and enclosed within other larger bounding volumes in a recursive fashion,
+ * eventually resulting in a tree structure with a single bounding volume at
+ * the top of the tree. Bounding volume hierarchies are used to support
+ * several operations on sets of geometric objects efficiently, such as in
+ * collision detection and ray tracing.
+ * </blockquote>
+ *
+ * Because ArborX uses Kokkos, Kokkos needs to be initialized and finalized
+ * before using this class.
+ */
+ class BVH
+ {
+ public:
+ /**
+ * Constructor. Use a vector of BoundingBox @p bounding_boxes as primitives.
+ */
+ template <int dim, typename Number>
+ BVH(const std::vector<BoundingBox<dim, Number>> &bounding_boxes);
+
+ /**
+ * Return the indices of those BoundingBox objects that satisfy the @p queries.
+ * Because @p queries can contain multiple queries, the function returns a pair
+ * of indices and offsets.
+ *
+ * Below is an example piece of code that does the following: Let us assume
+ * that we have a set of bounding boxes for objects -- say, the bounding
+ * boxes of each of the cells in a triangulation, or the bounding boxes for
+ * each of the parts of a triangulation in a parallel triangulation. We will
+ * store those in the `bvh_bounding_boxes` array below.
+ *
+ * Let us then also assume that we have a set of other bounding boxes, let's
+ * say for small objects that are moving around in our domain. We will put
+ * these bounding boxes into the `bb_intersect` array. The question we would
+ * then like to answer is the following: with which of the BVH bounding
+ * box(es) do each of the bb_intersect bounding boxes intersect? In other
+ * words, in which cell(s) or partition(s) are the particles?
+ *
+ * This query is answered by the following piece of code:
+ *
+ * @code
+ * const std::vector<BoundingBox<dim>> query_bounding_boxes = ...
+ * ArborXWrappers::BoundingBoxIntersect bb_intersect(query_bounding_boxes);
+ *
+ * const std::vector<BoundingBox<dim>> 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<int> bvh_bounding_box_indices;
+ * for (int i = offset[j]; i < offset[j+1]; ++i)
+ * bvh_bounding_box_indices.push_back(indices[i]);
+ * @endcode
+ *
+ * In many other applications, we are interested not only in finding which
+ * bounding boxes another bounding box lies in, but in fact which bounding
+ * boxes individual points lie in -- say, if instead of objects we have
+ * point-like particles moving around. In that case, we would need to answer
+ * a query for points, and this can be done as follows:
+ *
+ * @code
+ * const std::vector<Point<dim>> query_points = ...
+ * ArborXWrappers::PointIntersect 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
+ */
+ template <typename QueryType>
+ std::pair<std::vector<int>, std::vector<int>>
+ query(const QueryType &queries);
+
+ private:
+ /**
+ * Underlying ArborX object.
+ */
+ ArborX::BVH<Kokkos::HostSpace> bvh;
+ };
+
+
+
+ template <int dim, typename Number>
+ BVH::BVH(const std::vector<BoundingBox<dim, Number>> &bounding_boxes)
+ : bvh(Kokkos::DefaultHostExecutionSpace{}, bounding_boxes)
+ {}
+
+
+
+ template <typename QueryType>
+ std::pair<std::vector<int>, std::vector<int>>
+ BVH::query(const QueryType &queries)
+ {
+ Kokkos::View<int *, Kokkos::HostSpace> indices("indices", 0);
+
+ Kokkos::View<int *, Kokkos::HostSpace> offset("offset", 0);
+ ArborX::query(
+ bvh, Kokkos::DefaultHostExecutionSpace{}, queries, indices, offset);
+ std::vector<int> indices_vector;
+ indices_vector.insert(indices_vector.begin(),
+ indices.data(),
+ indices.data() + indices.extent(0));
+ std::vector<int> offset_vector;
+ offset_vector.insert(offset_vector.begin(),
+ offset.data(),
+ offset.data() + offset.extent(0));
+
+ return {indices_vector, offset_vector};
+ }
+} // namespace ArborXWrappers
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
+#endif
ADD_SUBDIRECTORY(non_matching)
ADD_SUBDIRECTORY(simplex)
ADD_SUBDIRECTORY(sundials)
+ADD_SUBDIRECTORY(arborx)
FOREACH(build ${DEAL_II_BUILD_TYPES})
STRING(TOLOWER ${build} build_lowercase)
--- /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.
+##
+## ---------------------------------------------------------------------
+
+INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
+
+SET(_src
+ access_traits.cc
+ )
+
+SET(_inst
+ access_traits.inst.in
+ )
+
+FILE(GLOB _header
+ ${CMAKE_SOURCE_DIR}/include/deal.II/arborx/*.h
+ )
+
+DEAL_II_ADD_LIBRARY(obj_arborx OBJECT ${_src} ${_header} ${_inst})
+EXPAND_INSTANTIATIONS(obj_arborx "${_inst}")
--- /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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/arborx/access_traits.h>
+
+#ifdef DEAL_II_WITH_ARBORX
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace ArborXWrappers
+{
+ // ------------------- PointIntersect ------------------- //
+ template <int dim, typename Number>
+ PointIntersect::PointIntersect(
+ const std::vector<dealii::Point<dim, Number>> &dim_points)
+ {
+ static_assert(dim != 1, "dim equal to one is not supported.");
+
+ const unsigned int size = dim_points.size();
+ points.reserve(size);
+ for (unsigned int i = 0; i < size; ++i)
+ {
+ points.emplace_back(static_cast<float>(dim_points[i][0]),
+ static_cast<float>(dim_points[i][1]),
+ dim == 2 ? 0.f :
+ static_cast<float>(dim_points[i][2]));
+ }
+ }
+
+
+
+ std::size_t
+ PointIntersect::size() const
+ {
+ return points.size();
+ }
+
+
+
+ const dealii::Point<3, float> &
+ PointIntersect::get(unsigned int i) const
+ {
+ return points[i];
+ }
+
+
+
+ // ------------------- BoundingBoxIntersect ------------------- //
+ template <int dim, typename Number>
+ BoundingBoxIntersect::BoundingBoxIntersect(
+ const std::vector<dealii::BoundingBox<dim, Number>> &bb)
+ {
+ const unsigned int size = bb.size();
+ bounding_boxes.reserve(size);
+ dealii::Point<3, float> min_corner_arborx(0., 0., 0.);
+ dealii::Point<3, float> max_corner_arborx(0., 0., 0.);
+ for (unsigned int i = 0; i < size; ++i)
+ {
+ auto boundary_points = bb[i].get_boundary_points();
+ dealii::Point<dim, Number> min_corner = boundary_points.first;
+ dealii::Point<dim, Number> max_corner = boundary_points.second;
+ for (int d = 0; d < dim; ++d)
+ {
+ min_corner_arborx[d] = static_cast<float>(min_corner[d]);
+ max_corner_arborx[d] = static_cast<float>(max_corner[d]);
+ }
+ bounding_boxes.emplace_back(
+ std::make_pair(min_corner_arborx, max_corner_arborx));
+ }
+ }
+
+
+
+ std::size_t
+ BoundingBoxIntersect::size() const
+ {
+ return bounding_boxes.size();
+ }
+
+
+
+ const dealii::BoundingBox<3, float> &
+ BoundingBoxIntersect::get(unsigned int i) const
+ {
+ return bounding_boxes[i];
+ }
+} // namespace ArborXWrappers
+
+DEAL_II_NAMESPACE_CLOSE
+
+namespace ArborX
+{
+ // ------------------- Point Primitives AccessTraits ------------------- //
+ template <int dim, typename Number>
+ std::size_t
+ AccessTraits<std::vector<dealii::Point<dim, Number>>, PrimitivesTag>::size(
+ const std::vector<dealii::Point<dim, Number>> &v)
+ {
+ return v.size();
+ }
+
+
+
+ template <int dim, typename Number>
+ Point
+ AccessTraits<std::vector<dealii::Point<dim, Number>>, PrimitivesTag>::get(
+ const std::vector<dealii::Point<dim, Number>> &v,
+ std::size_t i)
+ {
+ // ArborX assumes that the point coordinates use float and that the point
+ // is 3D
+ return {{static_cast<float>(v[i][0]),
+ static_cast<float>(v[i][1]),
+ dim == 2 ? 0 : static_cast<float>(v[i][2])}};
+ }
+
+
+
+ // ----------------- BoundingBox Primitives AccessTraits ----------------- //
+ template <int dim, typename Number>
+ std::size_t
+ AccessTraits<std::vector<dealii::BoundingBox<dim, Number>>, PrimitivesTag>::
+ size(const std::vector<dealii::BoundingBox<dim, Number>> &v)
+ {
+ return v.size();
+ }
+
+
+
+ template <int dim, typename Number>
+ Box
+ AccessTraits<std::vector<dealii::BoundingBox<dim, Number>>, PrimitivesTag>::
+ get(const std::vector<dealii::BoundingBox<dim, Number>> &v, std::size_t i)
+ {
+ const auto boundary_points = v[i].get_boundary_points();
+ const dealii::Point<dim, Number> min_corner = boundary_points.first;
+ const dealii::Point<dim, Number> max_corner = boundary_points.second;
+ // ArborX assumes that the bounding box coordinates use float and that the
+ // bounding box is 3D
+ return {{static_cast<float>(min_corner[0]),
+ static_cast<float>(min_corner[1]),
+ dim == 2 ? 0.f : static_cast<float>(min_corner[2])},
+ {static_cast<float>(max_corner[0]),
+ static_cast<float>(max_corner[1]),
+ dim == 2 ? 0.f : static_cast<float>(max_corner[2])}};
+ }
+} // namespace ArborX
+
+// ----------------------- Instantiations --------------------//
+# include "access_traits.inst"
+
+#endif
--- /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.
+//
+// ---------------------------------------------------------------------
+
+
+for (DIM : DIMENSIONS; SCALAR : REAL_SCALARS)
+ {
+#if DIM > 1
+ DEAL_II_NAMESPACE_OPEN
+ namespace ArborXWrappers
+ {
+ template PointIntersect::PointIntersect(
+ const std::vector<dealii::Point<DIM, SCALAR>> &dim_points);
+ template BoundingBoxIntersect::BoundingBoxIntersect(
+ const std::vector<dealii::BoundingBox<DIM, SCALAR>> &bb);
+ \}
+ DEAL_II_NAMESPACE_CLOSE
+
+ namespace ArborX
+ {
+ template struct AccessTraits<std::vector<dealii::Point<DIM, SCALAR>>,
+ PrimitivesTag>;
+ template struct AccessTraits<
+ std::vector<dealii::BoundingBox<DIM, SCALAR>>,
+ PrimitivesTag>;
+ \}
+#endif
+ }
projects / dealii.git / commitdiff