Unify and fix BoundingBoxtree constructor (C++) (#3545)

* Try to unify and fix bb_tree constructors

* Add permissive nanobind

* Apply suggestions from code review

cpp/dolfinx/geometry/BoundingBoxTree.h

@@ -13,6 +13,7 @@
 #include <cstdint>
 #include <dolfinx/mesh/utils.h>
 #include <mpi.h>
+#include <optional>
 #include <span>
 #include <string>
 #include <vector>
@@ -203,6 +204,11 @@ template <std::floating_point T>
 class BoundingBoxTree
 {
 private:
+  /// @brief Return a range of entity indices of the given topological dimension
+  /// (including ghosts).
+  /// @param topology Topology to get entities from.
+  /// @param tdim Dimension of the entities.
+  /// @return Range of local indices, including ghosts.
   static std::vector<std::int32_t> range(mesh::Topology& topology, int tdim)
   {
     topology.create_entities(tdim);
@@ -220,28 +226,45 @@ class BoundingBoxTree
   /// @param[in] tdim Topological dimension of the mesh entities to
   /// build the bounding box tree for.
   /// @param[in] entities List of entity indices (local to process) to
-  /// compute the bounding box for (may be empty, if none).
+  /// compute the bounding box for. If `std::nullopt`, the bounding box tree is
+  /// computed for all local entities (including ghosts) of the given `tdim`.
   /// @param[in] padding Value to pad (extend) the the bounding box of
   /// each entity by.
   BoundingBoxTree(const mesh::Mesh<T>& mesh, int tdim,
-                  std::span<const std::int32_t> entities, double padding = 0)
+                  std::optional<std::span<const std::int32_t>> entities
+                  = std::nullopt,
+                  double padding = 0)
       : _tdim(tdim)
   {
+    // Initialize entities of given dimension if they don't exist
+    mesh.topology_mutable()->create_entities(tdim);
+
+    // Get input entities. If not provided, get all local entities of the given
+    // dimension (including ghosts)
+    std::span<const std::int32_t> entities_span;
+    std::optional<std::vector<std::int32_t>> local_range(std::nullopt);
+    if (entities)
+      entities_span = entities.value();
+    else
+    {
+      local_range.emplace(range(*mesh.topology_mutable(), tdim));
+      entities_span = std::span<const std::int32_t>(local_range->data(),
+                                                    local_range->size());
+    }
+
     if (tdim < 0 or tdim > mesh.topology()->dim())
     {
       throw std::runtime_error(
           "Dimension must be non-negative and less than or "
           "equal to the topological dimension of the mesh");
     }
 
-    // Initialize entities of given dimension if they don't exist
-    mesh.topology_mutable()->create_entities(tdim);
     mesh.topology_mutable()->create_connectivity(tdim, mesh.topology()->dim());
 
     // Create bounding boxes for all mesh entities (leaves)
     std::vector<std::pair<std::array<T, 6>, std::int32_t>> leaf_bboxes;
-    leaf_bboxes.reserve(entities.size());
-    for (std::int32_t e : entities)
+    leaf_bboxes.reserve(entities_span.size());
+    for (std::int32_t e : entities_span)
     {
       std::array<T, 6> b = impl_bb::compute_bbox_of_entity(mesh, tdim, e);
       std::transform(b.cbegin(), std::next(b.cbegin(), 3), b.begin(),
@@ -257,20 +280,7 @@ class BoundingBoxTree
           = impl_bb::build_from_leaf(leaf_bboxes);
 
     spdlog::info("Computed bounding box tree with {} nodes for {} entities",
-                 num_bboxes(), entities.size());
-  }
-
-  /// Constructor
-  /// @param[in] mesh The mesh for building the bounding box tree
-  /// @param[in] tdim The topological dimension of the mesh entities to
-  /// build the bounding box tree for
-  /// @param[in] padding Value to pad (extend) the the bounding box of
-  /// each entity by.
-  BoundingBoxTree(const mesh::Mesh<T>& mesh, int tdim, T padding = 0)
-      : BoundingBoxTree::BoundingBoxTree(
-            mesh, tdim, range(mesh.topology_mutable(), tdim), padding)
-  {
-    // Do nothing
+                 num_bboxes(), entities_span.size());
   }
 
   /// Constructor @param[in] points Cloud of points, with associated

python/dolfinx/geometry.py

@@ -121,8 +121,6 @@ def bb_tree(
     map = mesh.topology.index_map(dim)
     if map is None:
         raise RuntimeError(f"Mesh entities of dimension {dim} have not been created.")
-    if entities is None:
-        entities = np.arange(map.size_local + map.num_ghosts, dtype=np.int32)
 
     dtype = mesh.geometry.x.dtype
     if np.issubdtype(dtype, np.float32):

python/dolfinx/wrappers/geometry.cpp

@@ -16,8 +16,10 @@
 #include <nanobind/nanobind.h>
 #include <nanobind/ndarray.h>
 #include <nanobind/stl/array.h>
+#include <nanobind/stl/optional.h>
 #include <nanobind/stl/tuple.h>
 #include <nanobind/stl/vector.h>
+#include <optional>
 #include <span>
 
 namespace nb = nanobind;
@@ -34,16 +36,22 @@ void declare_bbtree(nb::module_& m, std::string type)
           "__init__",
           [](dolfinx::geometry::BoundingBoxTree<T>* bbt,
              const dolfinx::mesh::Mesh<T>& mesh, int dim,
-             nb::ndarray<const std::int32_t, nb::ndim<1>, nb::c_contig>
+             std::optional<
+                 nb::ndarray<const std::int32_t, nb::ndim<1>, nb::c_contig>>
                  entities,
              double padding)
           {
-            new (bbt) dolfinx::geometry::BoundingBoxTree<T>(
-                mesh, dim,
-                std::span<const std::int32_t>(entities.data(), entities.size()),
-                padding);
+            std::optional<std::span<const std::int32_t>> ents
+                = entities ? std::span<const std::int32_t>(
+                                 entities->data(),
+                                 entities->data() + entities->size())
+                           : std::optional<std::span<const std::int32_t>>(
+                                 std::nullopt);
+
+            new (bbt)
+                dolfinx::geometry::BoundingBoxTree<T>(mesh, dim, ents, padding);
           },
-          nb::arg("mesh"), nb::arg("dim"), nb::arg("entities"),
+          nb::arg("mesh"), nb::arg("dim"), nb::arg("entities").none(),
           nb::arg("padding") = 0.0)
       .def_prop_ro("num_bboxes",
                    &dolfinx::geometry::BoundingBoxTree<T>::num_bboxes)