]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add wrapper for ArborX
authorBruno Turcksin <bruno.turcksin@gmail.com>
Fri, 8 Jan 2021 20:33:27 +0000 (20:33 +0000)
committerBruno Turcksin <bruno.turcksin@gmail.com>
Fri, 12 Feb 2021 14:21:38 +0000 (14:21 +0000)
include/deal.II/arborx/access_traits.h [new file with mode: 0644]
include/deal.II/arborx/bvh.h [new file with mode: 0644]
source/CMakeLists.txt
source/arborx/CMakeLists.txt [new file with mode: 0644]
source/arborx/access_traits.cc [new file with mode: 0644]
source/arborx/access_traits.inst.in [new file with mode: 0644]

diff --git a/include/deal.II/arborx/access_traits.h b/include/deal.II/arborx/access_traits.h
new file mode 100644 (file)
index 0000000..ce9204c
--- /dev/null
@@ -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 <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
diff --git a/include/deal.II/arborx/bvh.h b/include/deal.II/arborx/bvh.h
new file mode 100644 (file)
index 0000000..98fec0d
--- /dev/null
@@ -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 <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
index a6d85d8b68025eb3a16929cf0e20b7eb5dab8c50..460d78fa34b75357f6f023678cfe92a278dc99d0 100644 (file)
@@ -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 (file)
index 0000000..fa08009
--- /dev/null
@@ -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 (file)
index 0000000..b6612bd
--- /dev/null
@@ -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 <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
diff --git a/source/arborx/access_traits.inst.in b/source/arborx/access_traits.inst.in
new file mode 100644 (file)
index 0000000..c1d094a
--- /dev/null
@@ -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<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
+  }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.