Expose _bbox_coordinates of BoundingBoxTree to users

cpp/dolfinx/geometry/BoundingBoxTree.h

@@ -370,6 +370,20 @@ class BoundingBoxTree
   /// Return number of bounding boxes
   std::int32_t num_bboxes() const { return _bboxes.size() / 2; }
 
+  /// @brief Access coordinates of lower and upper corners of bounding boxes
+  /// (const version).
+  ///
+  /// @return The flattened row-major coordinate vector, where the shape is
+  /// `(2*num_bboxes, 3)`.
+  std::span<const T> bbox_coordinates() const { return _bbox_coordinates; }
+
+  /// @brief Access coordinates of lower and upper corners of bounding boxes
+  /// (non-const version)
+  ///
+  /// @return The flattened row-major coordinate vector, where the shape is
+  /// `(2*num_bboxes, 3)`.
+  std::span<T> bbox_coordinates() { return _bbox_coordinates; }
+
   /// Topological dimension of leaf entities
   int tdim() const { return _tdim; }
 

cpp/dolfinx/geometry/utils.h

@@ -197,15 +197,16 @@ constexpr bool is_leaf(std::array<int, 2> bbox)
 /// the bounds of the bounding box, b(0,i) <= x[i] <= b(1,i) for i = 0,
 /// 1, 2
 template <std::floating_point T>
-constexpr bool point_in_bbox(const std::array<T, 6>& b, std::span<const T, 3> x)
+constexpr bool point_in_bbox(std::span<const T, 6> b, std::span<const T, 3> x)
 {
   constexpr T rtol = 1e-14;
   bool in = true;
   for (std::size_t i = 0; i < 3; i++)
   {
     T eps = rtol * (b[i + 3] - b[i]);
-    in &= x[i] >= (b[i] - eps);
-    in &= x[i] <= (b[i + 3] + eps);
+    in &= (x[i] >= (b[i] - eps)) && (x[i] <= (b[i + 3] + eps));
+    if (not(in))
+      break;
   }
 
   return in;
@@ -310,10 +311,13 @@ void _compute_collisions_point(const geometry::BoundingBoxTree<T>& tree,
 {
   std::deque<std::int32_t> stack;
   std::int32_t next = tree.num_bboxes() - 1;
+  std::span<const T> bbox_coordinates = tree.bbox_coordinates();
+  auto view_bbox = [&bbox_coordinates](std::size_t node)
+  { return std::span<const T, 6>(bbox_coordinates.data() + 6 * node, 6); };
   while (next != -1)
   {
     const std::array<int, 2> bbox = tree.bbox(next);
-    if (is_leaf(bbox) and point_in_bbox(tree.get_bbox(next), p))
+    if (is_leaf(bbox) and point_in_bbox(view_bbox(next), p))
     {
       // If box is a leaf node then add it to the list of colliding
       // entities
@@ -324,8 +328,8 @@ void _compute_collisions_point(const geometry::BoundingBoxTree<T>& tree,
     {
       // Check whether the point collides with child nodes (left and
       // right)
-      bool left = point_in_bbox(tree.get_bbox(bbox[0]), p);
-      bool right = point_in_bbox(tree.get_bbox(bbox[1]), p);
+      bool left = point_in_bbox(view_bbox(bbox[0]), p);
+      bool right = point_in_bbox(view_bbox(bbox[1]), p);
       if (left and right)
       {
         // If the point collides with both child nodes, add the right

python/dolfinx/geometry.py

@@ -82,6 +82,16 @@ def num_bboxes(self) -> int:
         """Number of bounding boxes."""
         return self._cpp_object.num_bboxes
 
+    @property
+    def bbox_coordinates(self) -> typing.Union[npt.NDArray[np.float32], npt.NDArray[np.float64]]:
+        """Coordinates of lower and upper corners of bounding boxes.
+
+        Note:
+            Rows `2*ibbox` and `2*ibbox+1` correspond to the lower
+            and upper corners of bounding box `ibbox`, respectively.
+        """
+        return self._cpp_object.bbox_coordinates
+
     def get_bbox(self, i) -> npt.NDArray[np.floating]:
         """Get lower and upper corners of the ith bounding box.
 

python/dolfinx/wrappers/geometry.cpp

@@ -55,6 +55,18 @@ void declare_bbtree(nb::module_& m, std::string type)
           nb::arg("padding") = 0.0)
       .def_prop_ro("num_bboxes",
                    &dolfinx::geometry::BoundingBoxTree<T>::num_bboxes)
+      .def_prop_ro(
+          "bbox_coordinates",
+          [](dolfinx::geometry::BoundingBoxTree<T>& self)
+          {
+            std::span<T> bbox_coordinates = self.bbox_coordinates();
+            return nb::ndarray<T, nb::shape<-1, 3>, nb::numpy>(
+                bbox_coordinates.data(), {bbox_coordinates.size() / 3, 3});
+          },
+          nb::rv_policy::reference_internal,
+          "Return coordinates of bounding boxes."
+          "Row `2*ibbox` and `2*ibbox+1` correspond "
+          "to the lower and upper corners of bounding box `ibbox`.")
       .def(
           "get_bbox",
           [](const dolfinx::geometry::BoundingBoxTree<T>& self,

python/test/unit/geometry/test_bounding_box_tree.py

@@ -496,3 +496,23 @@ def test_surface_bbtree_collision(dtype):
 
     collisions = compute_collisions_trees(bbtree1, bbtree2)
     assert len(collisions) == 1
+
+
+@pytest.mark.parametrize("ct", [CellType.tetrahedron])
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_shift_bbtree(ct, dtype):
+    tdim = 3
+    mesh = create_unit_cube(MPI.COMM_WORLD, 3, 3, 3, ct, dtype=dtype)
+    bbtree = bb_tree(mesh, tdim, padding=0.0)
+    rng = np.random.default_rng(0)
+    points = rng.random((10, 3))
+
+    # Point-tree collisions pre-motion
+    collisions_pre = compute_collisions_points(bbtree, points)
+    # Shift everything
+    shift = np.array([1.0, 2.0, 3.0], dtype=dtype)
+    points[:] += shift
+    bbtree.bbox_coordinates[:] += shift
+
+    collisions_post = compute_collisions_points(bbtree, points)
+    assert (collisions_pre.array == collisions_post.array).all()