DOLFINx 0.9.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
BoundingBoxTree.h
1// Copyright (C) 2013-2022 Chris N. Richardson, Anders Logg, Garth N. Wells,
2// Jørgen S. Dokken, Sarah Roggendorf
3//
4// This file is part of DOLFINx (https://www.fenicsproject.org)
5//
6// SPDX-License-Identifier: LGPL-3.0-or-later
7
8#pragma once
9
10#include <algorithm>
11#include <array>
12#include <cassert>
13#include <cstdint>
14#include <dolfinx/mesh/utils.h>
15#include <mpi.h>
16#include <span>
17#include <string>
18#include <vector>
19
21{
22namespace impl_bb
23{
24//-----------------------------------------------------------------------------
25// Compute bounding box of mesh entity. The bounding box is defined by (lower
26// left corner, top right corner). Storage flattened row-major
27template <std::floating_point T>
28std::array<T, 6> compute_bbox_of_entity(const mesh::Mesh<T>& mesh, int dim,
29 std::int32_t index)
30{
31 // Get the geometrical indices for the mesh entity
32 std::span<const T> xg = mesh.geometry().x();
33
34 // FIXME: return of small dynamic array is expensive
35 std::span<const std::int32_t> entity(&index, 1);
36 const std::vector<std::int32_t> vertex_indices
37 = mesh::entities_to_geometry(mesh, dim, entity, false);
38
39 std::array<T, 6> b;
40 auto b0 = std::span(b).template subspan<0, 3>();
41 auto b1 = std::span(b).template subspan<3, 3>();
42 std::copy_n(std::next(xg.begin(), 3 * vertex_indices.front()), 3, b0.begin());
43 std::copy_n(std::next(xg.begin(), 3 * vertex_indices.front()), 3, b1.begin());
44
45 // Compute min and max over vertices
46 for (std::int32_t local_vertex : vertex_indices)
47 {
48 for (std::size_t j = 0; j < 3; ++j)
49 {
50 b0[j] = std::min(b0[j], xg[3 * local_vertex + j]);
51 b1[j] = std::max(b1[j], xg[3 * local_vertex + j]);
52 }
53 }
54
55 return b;
56}
57//-----------------------------------------------------------------------------
58// Compute bounding box of bounding boxes. Each bounding box is defined as a
59// tuple (corners, entity_index). The corners of the bounding box is flattened
60// row-major as (lower left corner, top right corner).
61template <std::floating_point T>
62std::array<T, 6> compute_bbox_of_bboxes(
63 std::span<const std::pair<std::array<T, 6>, std::int32_t>> leaf_bboxes)
64{
65 // Compute min and max over remaining boxes
66 std::array<T, 6> b = leaf_bboxes.front().first;
67 for (auto [box, _] : leaf_bboxes)
68 {
69 std::transform(box.cbegin(), std::next(box.cbegin(), 3), b.cbegin(),
70 b.begin(), [](auto a, auto b) { return std::min(a, b); });
71 std::transform(std::next(box.cbegin(), 3), box.cend(),
72 std::next(b.cbegin(), 3), std::next(b.begin(), 3),
73 [](auto a, auto b) { return std::max(a, b); });
74 }
75
76 return b;
77}
78//------------------------------------------------------------------------------
79template <std::floating_point T>
80std::int32_t _build_from_leaf(
81 std::span<std::pair<std::array<T, 6>, std::int32_t>> leaf_bboxes,
82 std::vector<int>& bboxes, std::vector<T>& bbox_coordinates)
83{
84 if (leaf_bboxes.size() == 1)
85 {
86 // Reached leaf
87
88 // Get bounding box coordinates for leaf
89 const auto [b, entity_index] = leaf_bboxes.front();
90
91 // Store bounding box data
92 bboxes.push_back(entity_index);
93 bboxes.push_back(entity_index);
94 std::copy_n(b.begin(), 6, std::back_inserter(bbox_coordinates));
95 return bboxes.size() / 2 - 1;
96 }
97 else
98 {
99 // Compute bounding box of all bounding boxes
100 std::array b = compute_bbox_of_bboxes<T>(leaf_bboxes);
101
102 // Sort bounding boxes along longest axis
103 std::array<T, 3> b_diff;
104 std::transform(std::next(b.cbegin(), 3), b.cend(), b.cbegin(),
105 b_diff.begin(), std::minus<T>());
106 const std::size_t axis = std::distance(
107 b_diff.begin(), std::max_element(b_diff.begin(), b_diff.end()));
108
109 auto middle = std::next(leaf_bboxes.begin(), leaf_bboxes.size() / 2);
110 std::nth_element(leaf_bboxes.begin(), middle, leaf_bboxes.end(),
111 [axis](auto& p0, auto& p1) -> bool
112 {
113 auto x0 = p0.first[axis] + p0.first[3 + axis];
114 auto x1 = p1.first[axis] + p1.first[3 + axis];
115 return x0 < x1;
116 });
117
118 // Split bounding boxes into two groups and call recursively
119 assert(!leaf_bboxes.empty());
120 std::size_t part = leaf_bboxes.size() / 2;
121 std::int32_t bbox0
122 = _build_from_leaf(leaf_bboxes.first(part), bboxes, bbox_coordinates);
123 std::int32_t bbox1 = _build_from_leaf(
124 leaf_bboxes.last(leaf_bboxes.size() - part), bboxes, bbox_coordinates);
125
126 // Store bounding box data. Note that root box will be added last.
127 bboxes.push_back(bbox0);
128 bboxes.push_back(bbox1);
129 std::copy_n(b.begin(), 6, std::back_inserter(bbox_coordinates));
130 return bboxes.size() / 2 - 1;
131 }
132}
133//-----------------------------------------------------------------------------
134template <std::floating_point T>
135std::pair<std::vector<std::int32_t>, std::vector<T>> build_from_leaf(
136 std::vector<std::pair<std::array<T, 6>, std::int32_t>>& leaf_bboxes)
137{
138 std::vector<std::int32_t> bboxes;
139 std::vector<T> bbox_coordinates;
140 impl_bb::_build_from_leaf<T>(leaf_bboxes, bboxes, bbox_coordinates);
141 return {std::move(bboxes), std::move(bbox_coordinates)};
142}
143//-----------------------------------------------------------------------------
144template <std::floating_point T>
145std::int32_t
146_build_from_point(std::span<std::pair<std::array<T, 3>, std::int32_t>> points,
147 std::vector<std::int32_t>& bboxes,
148 std::vector<T>& bbox_coordinates)
149{
150 // Reached leaf
151 if (points.size() == 1)
152 {
153 // Store bounding box data
154
155 // Index of entity contained in leaf
156 const std::int32_t c1 = points[0].second;
157 bboxes.push_back(c1);
158 bboxes.push_back(c1);
159 bbox_coordinates.insert(bbox_coordinates.end(), points[0].first.begin(),
160 points[0].first.end());
161 bbox_coordinates.insert(bbox_coordinates.end(), points[0].first.begin(),
162 points[0].first.end());
163 return bboxes.size() / 2 - 1;
164 }
165
166 // Compute bounding box of all points
167 auto [min, max] = std::ranges::minmax_element(points);
168 std::array<T, 3> b0 = min->first;
169 std::array<T, 3> b1 = max->first;
170
171 // Sort bounding boxes along longest axis
172 std::array<T, 3> b_diff;
173 std::ranges::transform(b1, b0, b_diff.begin(), std::minus<T>());
174 const std::size_t axis
175 = std::distance(b_diff.begin(), std::ranges::max_element(b_diff));
176
177 auto middle = std::next(points.begin(), points.size() / 2);
178 std::nth_element(points.begin(), middle, points.end(),
179 [axis](auto& p0, auto&& p1) -> bool
180 { return p0.first[axis] < p1.first[axis]; });
181
182 // Split bounding boxes into two groups and call recursively
183 assert(!points.empty());
184 std::size_t part = points.size() / 2;
185 std::int32_t bbox0
186 = _build_from_point(points.first(part), bboxes, bbox_coordinates);
187 std::int32_t bbox1 = _build_from_point(points.last(points.size() - part),
188 bboxes, bbox_coordinates);
189
190 // Store bounding box data. Note that root box will be added last.
191 bboxes.push_back(bbox0);
192 bboxes.push_back(bbox1);
193 bbox_coordinates.insert(bbox_coordinates.end(), b0.begin(), b0.end());
194 bbox_coordinates.insert(bbox_coordinates.end(), b1.begin(), b1.end());
195 return bboxes.size() / 2 - 1;
196}
197//-----------------------------------------------------------------------------
198} // namespace impl_bb
199
202template <std::floating_point T>
204{
205private:
206 static std::vector<std::int32_t> range(mesh::Topology& topology, int tdim)
207 {
208 topology.create_entities(tdim);
209 auto map = topology.index_map(tdim);
210 assert(map);
211 const std::int32_t num_entities = map->size_local() + map->num_ghosts();
212 std::vector<std::int32_t> r(num_entities);
213 std::iota(r.begin(), r.end(), 0);
214 return r;
215 }
216
217public:
227 std::span<const std::int32_t> entities, double padding = 0)
228 : _tdim(tdim)
229 {
230 if (tdim < 0 or tdim > mesh.topology()->dim())
231 {
232 throw std::runtime_error(
233 "Dimension must be non-negative and less than or "
234 "equal to the topological dimension of the mesh");
235 }
236
237 // Initialize entities of given dimension if they don't exist
238 mesh.topology_mutable()->create_entities(tdim);
239 mesh.topology_mutable()->create_connectivity(tdim, mesh.topology()->dim());
240
241 // Create bounding boxes for all mesh entities (leaves)
242 std::vector<std::pair<std::array<T, 6>, std::int32_t>> leaf_bboxes;
243 leaf_bboxes.reserve(entities.size());
244 for (std::int32_t e : entities)
245 {
246 std::array<T, 6> b = impl_bb::compute_bbox_of_entity(mesh, tdim, e);
247 std::transform(b.cbegin(), std::next(b.cbegin(), 3), b.begin(),
248 [padding](auto x) { return x - padding; });
249 std::transform(std::next(b.begin(), 3), b.end(), std::next(b.begin(), 3),
250 [padding](auto x) { return x + padding; });
251 leaf_bboxes.emplace_back(b, e);
252 }
253
254 // Recursively build the bounding box tree from the leaves
255 if (!leaf_bboxes.empty())
256 std::tie(_bboxes, _bbox_coordinates)
257 = impl_bb::build_from_leaf(leaf_bboxes);
258
259 spdlog::info("Computed bounding box tree with {} nodes for {} entities",
260 num_bboxes(), entities.size());
261 }
262
269 BoundingBoxTree(const mesh::Mesh<T>& mesh, int tdim, T padding = 0)
271 mesh, tdim, range(mesh.topology_mutable(), tdim), padding)
272 {
273 // Do nothing
274 }
275
278 BoundingBoxTree(std::vector<std::pair<std::array<T, 3>, std::int32_t>> points)
279 : _tdim(0)
280 {
281 // Recursively build the bounding box tree from the leaves
282 if (!points.empty())
283 {
284 _bboxes.clear();
285 impl_bb::_build_from_point(std::span(points), _bboxes, _bbox_coordinates);
286 }
287
288 spdlog::info("Computed bounding box tree with {} nodes for {} points.",
289 num_bboxes(), points.size());
290 }
291
294
296 BoundingBoxTree(const BoundingBoxTree& tree) = delete;
297
300
302 BoundingBoxTree& operator=(const BoundingBoxTree& other) = default;
303
305 ~BoundingBoxTree() = default;
306
312 std::array<T, 6> get_bbox(std::size_t node) const
313 {
314 std::array<T, 6> x;
315 std::copy_n(_bbox_coordinates.data() + 6 * node, 6, x.begin());
316 return x;
317 }
318
325 {
326 // Build tree for each rank
327 const int mpi_size = dolfinx::MPI::size(comm);
328
329 // Send root node coordinates to all processes
330 // This is to counteract the fact that a process might have 0 bounding box
331 // causing false positives on process collisions around (0,0,0)
332 constexpr T max_val = std::numeric_limits<T>::max();
333 std::array<T, 6> send_bbox
334 = {max_val, max_val, max_val, max_val, max_val, max_val};
335 if (num_bboxes() > 0)
336 std::copy_n(std::prev(_bbox_coordinates.end(), 6), 6, send_bbox.begin());
337 std::vector<T> recv_bbox(mpi_size * 6);
338 MPI_Allgather(send_bbox.data(), 6, dolfinx::MPI::mpi_type<T>(),
339 recv_bbox.data(), 6, dolfinx::MPI::mpi_type<T>(), comm);
340
341 std::vector<std::pair<std::array<T, 6>, std::int32_t>> _recv_bbox(mpi_size);
342 for (std::size_t i = 0; i < _recv_bbox.size(); ++i)
343 {
344 std::copy_n(std::next(recv_bbox.begin(), 6 * i), 6,
345 _recv_bbox[i].first.begin());
346 _recv_bbox[i].second = i;
347 }
348
349 auto [global_bboxes, global_coords] = impl_bb::build_from_leaf(_recv_bbox);
350 BoundingBoxTree global_tree(std::move(global_bboxes),
351 std::move(global_coords));
352
353 spdlog::info("Computed global bounding box tree with {} boxes.",
354 global_tree.num_bboxes());
355
356 return global_tree;
357 }
358
360 std::int32_t num_bboxes() const { return _bboxes.size() / 2; }
361
363 int tdim() const { return _tdim; }
364
366 std::string str() const
367 {
368 std::stringstream s;
369 tree_print(s, _bboxes.size() / 2 - 1);
370 return s.str();
371 }
372
380 std::array<std::int32_t, 2> bbox(std::size_t node) const
381 {
382 assert(2 * node + 1 < _bboxes.size());
383 return {_bboxes[2 * node], _bboxes[2 * node + 1]};
384 }
385
386private:
387 // Constructor
388 BoundingBoxTree(std::vector<std::int32_t>&& bboxes,
389 std::vector<T>&& bbox_coords)
390 : _tdim(0), _bboxes(bboxes), _bbox_coordinates(bbox_coords)
391 {
392 // Do nothing
393 }
394
395 // Topological dimension of leaf entities
396 int _tdim;
397
398 // Print out recursively, for debugging
399 void tree_print(std::stringstream& s, std::int32_t i) const
400 {
401 s << "[";
402 for (std::size_t j = 0; j < 2; ++j)
403 {
404 for (std::size_t k = 0; k < 3; ++k)
405 s << _bbox_coordinates[6 * i + j * 3 + k] << " ";
406 if (j == 0)
407 s << "]->"
408 << "[";
409 }
410 s << "]\n";
411
412 if (_bboxes[2 * i] == _bboxes[2 * i + 1])
413 s << "leaf containing entity (" << _bboxes[2 * i + 1] << ")";
414 else
415 {
416 s << "{";
417 tree_print(s, _bboxes[2 * i]);
418 s << ", \n";
419 tree_print(s, _bboxes[2 * i + 1]);
420 s << "}\n";
421 }
422 }
423
424 // List of bounding boxes (parent-child-entity relations)
425 std::vector<std::int32_t> _bboxes;
426
427 // List of bounding box coordinates
428 std::vector<T> _bbox_coordinates;
429};
430} // namespace dolfinx::geometry
Definition BoundingBoxTree.h:204
BoundingBoxTree(const mesh::Mesh< T > &mesh, int tdim, T padding=0)
Definition BoundingBoxTree.h:269
BoundingBoxTree(std::vector< std::pair< std::array< T, 3 >, std::int32_t > > points)
Definition BoundingBoxTree.h:278
BoundingBoxTree(const mesh::Mesh< T > &mesh, int tdim, std::span< const std::int32_t > entities, double padding=0)
Definition BoundingBoxTree.h:226
BoundingBoxTree & operator=(const BoundingBoxTree &other)=default
Copy assignment.
BoundingBoxTree & operator=(BoundingBoxTree &&other)=default
Move assignment.
BoundingBoxTree create_global_tree(MPI_Comm comm) const
Definition BoundingBoxTree.h:324
int tdim() const
Topological dimension of leaf entities.
Definition BoundingBoxTree.h:363
BoundingBoxTree(BoundingBoxTree &&tree)=default
Move constructor.
std::int32_t num_bboxes() const
Return number of bounding boxes.
Definition BoundingBoxTree.h:360
std::array< T, 6 > get_bbox(std::size_t node) const
Return bounding box coordinates for a given node in the tree,.
Definition BoundingBoxTree.h:312
std::array< std::int32_t, 2 > bbox(std::size_t node) const
Definition BoundingBoxTree.h:380
BoundingBoxTree(const BoundingBoxTree &tree)=delete
Copy constructor.
~BoundingBoxTree()=default
Destructor.
std::string str() const
Print out for debugging.
Definition BoundingBoxTree.h:366
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition Mesh.h:23
Topology stores the topology of a mesh, consisting of mesh entities and connectivity (incidence relat...
Definition Topology.h:44
std::shared_ptr< const common::IndexMap > index_map(int dim) const
Get the IndexMap that described the parallel distribution of the mesh entities.
Definition Topology.cpp:815
std::int32_t create_entities(int dim)
Create entities of given topological dimension.
Definition Topology.cpp:831
Functions supporting mesh operations.
int size(MPI_Comm comm)
Definition MPI.cpp:72
constexpr MPI_Datatype mpi_type()
MPI Type.
Definition MPI.h:273
Geometry data structures and algorithms.
Definition BoundingBoxTree.h:21
std::vector< std::int32_t > entities_to_geometry(const Mesh< T > &mesh, int dim, std::span< const std::int32_t > entities, bool permute=false)
Compute the geometry degrees of freedom associated with the closure of a given set of cell entities.
Definition utils.h:633