From c5fbbbfaa73349ec5a01d9f40a2fb94da795b379 Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Fri, 8 Jan 2021 20:33:27 +0000 Subject: [PATCH] Add wrapper for ArborX --- include/deal.II/arborx/access_traits.h | 250 +++++++++++++++++++++++++ include/deal.II/arborx/bvh.h | 163 ++++++++++++++++ source/CMakeLists.txt | 1 + source/arborx/CMakeLists.txt | 31 +++ source/arborx/access_traits.cc | 164 ++++++++++++++++ source/arborx/access_traits.inst.in | 39 ++++ 6 files changed, 648 insertions(+) create mode 100644 include/deal.II/arborx/access_traits.h create mode 100644 include/deal.II/arborx/bvh.h create mode 100644 source/arborx/CMakeLists.txt create mode 100644 source/arborx/access_traits.cc create mode 100644 source/arborx/access_traits.inst.in diff --git a/include/deal.II/arborx/access_traits.h b/include/deal.II/arborx/access_traits.h new file mode 100644 index 0000000000..ce9204c1f8 --- /dev/null +++ b/include/deal.II/arborx/access_traits.h @@ -0,0 +1,250 @@ +// --------------------------------------------------------------------- +// +// 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 + +#ifdef DEAL_II_WITH_ARBORX +# include + +# include + + +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 + PointIntersect(const std::vector> &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> 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 + BoundingBoxIntersect( + const std::vector> &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> 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 as + * primitive. + */ + template + struct AccessTraits>, PrimitivesTag> + { + using memory_space = Kokkos::HostSpace; + + /** + * Return the size of the vector @p v. + */ + static std::size_t + size(const std::vector> &v); + + /** + * Return an ArborX::Point from the dealii::Point `v[i]`. + */ + static Point + get(const std::vector> &v, std::size_t i); + }; + + + + /** + * This struct allows ArborX to use std::vector as + * primitive. + */ + template + struct AccessTraits>, + PrimitivesTag> + { + using memory_space = Kokkos::HostSpace; + + /** + * Return the size of the vector @p v. + */ + static std::size_t + size(const std::vector> &v); + + /** + * Return an ArborX::Box from the dealii::BoundingBox `v[i]`. + */ + static Box + get(const std::vector> &v, std::size_t i); + }; + + + + /** + * This struct allows ArborX to use PointIntersect in a predicate. + */ + template <> + struct AccessTraits + { + 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 + { + 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::size( + const dealii::ArborXWrappers::PointIntersect &pt_intersect) + { + return pt_intersect.size(); + } + + + + inline auto + AccessTraits::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:: + size(const dealii::ArborXWrappers::BoundingBoxIntersect &bb_intersect) + { + return bb_intersect.size(); + } + + + + inline auto + AccessTraits:: + 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 diff --git a/include/deal.II/arborx/bvh.h b/include/deal.II/arborx/bvh.h new file mode 100644 index 0000000000..98fec0d674 --- /dev/null +++ b/include/deal.II/arborx/bvh.h @@ -0,0 +1,163 @@ +// --------------------------------------------------------------------- +// +// 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 + +#ifdef DEAL_II_WITH_ARBORX +# include + +# include +# include + +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): + *
+ * 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. + *
+ * + * 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 + BVH(const std::vector> &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> query_bounding_boxes = ... + * ArborXWrappers::BoundingBoxIntersect bb_intersect(query_bounding_boxes); + * + * const std::vector> bvh_bounding_boxes = ... + * ArborxWrappers::BVH bvh(bvh_bounding_boxes); + * + * auto [indices, offset] = bvh.query(bb_intersect); + * @endcode + * + * The elements of `bvh_bounding_boxes` that intersect the `j`th BoundingBox + * of `query_bounding_boxes` are given by: + * + * @code + * std::vector bvh_bounding_box_indices; + * for (int i = offset[j]; i < offset[j+1]; ++i) + * bvh_bounding_box_indices.push_back(indices[i]); + * @endcode + * + * In many other applications, we are interested not only in finding which + * bounding boxes another bounding box lies in, but in fact which bounding + * boxes individual points lie in -- say, if instead of objects we have + * point-like particles moving around. In that case, we would need to answer + * a query for points, and this can be done as follows: + * + * @code + * const std::vector> query_points = ... + * ArborXWrappers::PointIntersect pt_intersect(query_points); + * + * const std::vector> bvh_bounding_boxes = ... + * ArborxWrappers::BVH bvh(bvh_bounding_boxes); + * + * auto [indices, offset] = bvh.query(pt_intersect); + * @endcode + */ + template + std::pair, std::vector> + query(const QueryType &queries); + + private: + /** + * Underlying ArborX object. + */ + ArborX::BVH bvh; + }; + + + + template + BVH::BVH(const std::vector> &bounding_boxes) + : bvh(Kokkos::DefaultHostExecutionSpace{}, bounding_boxes) + {} + + + + template + std::pair, std::vector> + BVH::query(const QueryType &queries) + { + Kokkos::View indices("indices", 0); + + Kokkos::View offset("offset", 0); + ArborX::query( + bvh, Kokkos::DefaultHostExecutionSpace{}, queries, indices, offset); + std::vector indices_vector; + indices_vector.insert(indices_vector.begin(), + indices.data(), + indices.data() + indices.extent(0)); + std::vector offset_vector; + offset_vector.insert(offset_vector.begin(), + offset.data(), + offset.data() + offset.extent(0)); + + return {indices_vector, offset_vector}; + } +} // namespace ArborXWrappers + +DEAL_II_NAMESPACE_CLOSE + +#endif +#endif diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index a6d85d8b68..460d78fa34 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -73,6 +73,7 @@ ADD_SUBDIRECTORY(optimization/rol) ADD_SUBDIRECTORY(non_matching) ADD_SUBDIRECTORY(simplex) ADD_SUBDIRECTORY(sundials) +ADD_SUBDIRECTORY(arborx) FOREACH(build ${DEAL_II_BUILD_TYPES}) STRING(TOLOWER ${build} build_lowercase) diff --git a/source/arborx/CMakeLists.txt b/source/arborx/CMakeLists.txt new file mode 100644 index 0000000000..fa08009cde --- /dev/null +++ b/source/arborx/CMakeLists.txt @@ -0,0 +1,31 @@ +## --------------------------------------------------------------------- +## +## 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}") diff --git a/source/arborx/access_traits.cc b/source/arborx/access_traits.cc new file mode 100644 index 0000000000..b6612bd9a8 --- /dev/null +++ b/source/arborx/access_traits.cc @@ -0,0 +1,164 @@ +// --------------------------------------------------------------------- +// +// 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 + +#ifdef DEAL_II_WITH_ARBORX + +DEAL_II_NAMESPACE_OPEN + +namespace ArborXWrappers +{ + // ------------------- PointIntersect ------------------- // + template + PointIntersect::PointIntersect( + const std::vector> &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(dim_points[i][0]), + static_cast(dim_points[i][1]), + dim == 2 ? 0.f : + static_cast(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 + BoundingBoxIntersect::BoundingBoxIntersect( + const std::vector> &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 min_corner = boundary_points.first; + dealii::Point max_corner = boundary_points.second; + for (int d = 0; d < dim; ++d) + { + min_corner_arborx[d] = static_cast(min_corner[d]); + max_corner_arborx[d] = static_cast(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 + std::size_t + AccessTraits>, PrimitivesTag>::size( + const std::vector> &v) + { + return v.size(); + } + + + + template + Point + AccessTraits>, PrimitivesTag>::get( + const std::vector> &v, + std::size_t i) + { + // ArborX assumes that the point coordinates use float and that the point + // is 3D + return {{static_cast(v[i][0]), + static_cast(v[i][1]), + dim == 2 ? 0 : static_cast(v[i][2])}}; + } + + + + // ----------------- BoundingBox Primitives AccessTraits ----------------- // + template + std::size_t + AccessTraits>, PrimitivesTag>:: + size(const std::vector> &v) + { + return v.size(); + } + + + + template + Box + AccessTraits>, PrimitivesTag>:: + get(const std::vector> &v, std::size_t i) + { + const auto boundary_points = v[i].get_boundary_points(); + const dealii::Point min_corner = boundary_points.first; + const dealii::Point max_corner = boundary_points.second; + // ArborX assumes that the bounding box coordinates use float and that the + // bounding box is 3D + return {{static_cast(min_corner[0]), + static_cast(min_corner[1]), + dim == 2 ? 0.f : static_cast(min_corner[2])}, + {static_cast(max_corner[0]), + static_cast(max_corner[1]), + dim == 2 ? 0.f : static_cast(max_corner[2])}}; + } +} // namespace ArborX + +// ----------------------- Instantiations --------------------// +# include "access_traits.inst" + +#endif diff --git a/source/arborx/access_traits.inst.in b/source/arborx/access_traits.inst.in new file mode 100644 index 0000000000..c1d094aed3 --- /dev/null +++ b/source/arborx/access_traits.inst.in @@ -0,0 +1,39 @@ +// --------------------------------------------------------------------- +// +// 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> &dim_points); + template BoundingBoxIntersect::BoundingBoxIntersect( + const std::vector> &bb); + \} + DEAL_II_NAMESPACE_CLOSE + + namespace ArborX + { + template struct AccessTraits>, + PrimitivesTag>; + template struct AccessTraits< + std::vector>, + PrimitivesTag>; + \} +#endif + } -- 2.39.5