version 3.10-dev
boundingboxtree.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3//
4// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5// SPDX-License-Identifier: GPL-3.0-or-later
6//
19#ifndef DUMUX_GEOMETRY_BOUNDINGBOXTREE_HH
20#define DUMUX_GEOMETRY_BOUNDINGBOXTREE_HH
21
22#include <vector>
23#include <array>
24#include <algorithm>
25#include <memory>
26#include <numeric>
27#include <type_traits>
28#include <iostream>
29
30#include <dune/common/promotiontraits.hh>
31#include <dune/common/timer.hh>
32#include <dune/common/fvector.hh>
33
34namespace Dumux {
35
54template <class GeometricEntitySet>
56{
57 enum { dimworld = GeometricEntitySet::dimensionworld };
58 using ctype = typename GeometricEntitySet::ctype;
59
65 struct BoundingBoxNode
66 {
67 std::size_t child0;
68 std::size_t child1;
69 };
70
71public:
73 using EntitySet = GeometricEntitySet;
74
76 BoundingBoxTree() = default;
77
79 BoundingBoxTree(std::shared_ptr<const GeometricEntitySet> set)
80 { build(set); }
81
83 void build(std::shared_ptr<const GeometricEntitySet> set)
84 {
85 // set the pointer to the entity set
86 entitySet_ = set;
87
88 // clear all internal data
89 boundingBoxNodes_.clear();
90 boundingBoxCoordinates_.clear();
91
92 // start the timer
93 Dune::Timer timer;
94
95 // Create bounding boxes for all elements
96 const auto numLeaves = set->size();
97
98 // reserve enough space for the nodes and the coordinates
99 const auto numNodes = 2*numLeaves - 1;
100 boundingBoxNodes_.reserve(numNodes);
101 boundingBoxCoordinates_.reserve(numNodes*2*dimworld);
102
103 // create a vector for leaf boxes (min and max for all dims)
104 std::vector<ctype> leafBoxes(2*dimworld*numLeaves);
105
106 for (const auto& geometricEntity : *set)
107 computeEntityBoundingBox_(leafBoxes.data() + 2*dimworld*set->index(geometricEntity), geometricEntity);
108
109 // create the leaf partition, the set of available indices (to be sorted)
110 std::vector<std::size_t> leafPartition(numLeaves);
111 std::iota(leafPartition.begin(), leafPartition.end(), 0);
112
113 // Recursively build the bounding box tree
114 build_(leafBoxes, leafPartition.begin(), leafPartition.end());
115
116 // We are done, log output
117 std::cout << "Computed bounding box tree with " << numBoundingBoxes()
118 << " nodes for " << numLeaves << " grid entities in "
119 << timer.stop() << " seconds." << std::endl;
120 }
121
123 const EntitySet& entitySet() const
124 { return *entitySet_; }
125
129
131 const BoundingBoxNode& getBoundingBoxNode(std::size_t nodeIdx) const
132 { return boundingBoxNodes_[nodeIdx]; }
133
135 const ctype* getBoundingBoxCoordinates(std::size_t nodeIdx) const
136 { return boundingBoxCoordinates_.data() + 2*dimworld*nodeIdx; }
137
139 std::size_t numBoundingBoxes() const
140 { return boundingBoxNodes_.size(); }
141
144 bool isLeaf(const BoundingBoxNode& node, std::size_t nodeIdx) const
145 { return node.child0 == nodeIdx; }
146
147private:
148
150 std::vector<BoundingBoxNode> boundingBoxNodes_;
151
153 std::vector<ctype> boundingBoxCoordinates_;
154
156 std::shared_ptr<const EntitySet> entitySet_;
157
159 template <class Entity>
160 void computeEntityBoundingBox_(ctype* b, const Entity& entity) const
161 {
162 // get the bounding box coordinates
163 ctype* xMin = b;
164 ctype* xMax = b + dimworld;
165
166 // get mesh entity data
167 auto geometry = entity.geometry();
168
169 // Get coordinates of first vertex
170 auto corner = geometry.corner(0);
171 for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx)
172 xMin[dimIdx] = xMax[dimIdx] = corner[dimIdx];
173
174 // Compute the min and max over the remaining vertices
175 for (std::size_t cornerIdx = 1; cornerIdx < geometry.corners(); ++cornerIdx)
176 {
177 corner = geometry.corner(cornerIdx);
178 for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx)
179 {
180 using std::max;
181 using std::min;
182 xMin[dimIdx] = min(xMin[dimIdx], corner[dimIdx]);
183 xMax[dimIdx] = max(xMax[dimIdx], corner[dimIdx]);
184 }
185 }
186 }
187
189 std::size_t build_(const std::vector<ctype>& leafBoxes,
190 const std::vector<std::size_t>::iterator& begin,
191 const std::vector<std::size_t>::iterator& end)
192 {
193 assert(begin < end);
194
195 // If we reached the end of the recursion, i.e. only a leaf box is left
196 if (end - begin == 1)
197 {
198 // Get the bounding box coordinates for the leaf
199 const std::size_t leafNodeIdx = *begin;
200 const auto beginCoords = leafBoxes.begin() + 2*dimworld*leafNodeIdx;
201 const auto endCoords = beginCoords + 2*dimworld;
202
203 // Store the data in the bounding box
204 // leaf nodes are indicated by setting child0 to
205 // the node itself and child1 to the index of the entity in the bounding box.
206 return addBoundingBox_(BoundingBoxNode{numBoundingBoxes(), leafNodeIdx}, beginCoords, endCoords);
207 }
208
209 // Compute the bounding box of all bounding boxes in the range [begin, end]
210 const auto bCoords = computeBBoxOfBBoxes_(leafBoxes, begin, end);
211
212 // sort bounding boxes along the longest axis
213 const auto axis = computeLongestAxis_(bCoords);
214
215 // nth_element sorts the range to make sure that middle points to the coordinate median in axis direction
216 // this is the most expensive part of the algorithm
217 auto middle = begin + (end - begin)/2;
218 std::nth_element(begin, middle, end, [&leafBoxes, axis](std::size_t i, std::size_t j)
219 {
220 const ctype* bi = leafBoxes.data() + 2*dimworld*i;
221 const ctype* bj = leafBoxes.data() + 2*dimworld*j;
222 return bi[axis] + bi[axis + dimworld] < bj[axis] + bj[axis + dimworld];
223 });
224
225 // split the bounding boxes into two at the middle iterator and call build recursively, each
226 // call resulting in a new node of this bounding box, i.e. the root will be added at the end of the process.
227 return addBoundingBox_(BoundingBoxNode{build_(leafBoxes, begin, middle), build_(leafBoxes, middle, end)},
228 bCoords.begin(), bCoords.end());
229 }
230
232 template <class Iterator>
233 std::size_t addBoundingBox_(BoundingBoxNode&& node,
234 const Iterator& coordBegin,
235 const Iterator& coordEnd)
236 {
237 // Add the bounding box
238 boundingBoxNodes_.emplace_back(node);
239
240 // Add the bounding box's coordinates
241 boundingBoxCoordinates_.insert(boundingBoxCoordinates_.end(), coordBegin, coordEnd);
242
243 // return the index of the new node
244 return boundingBoxNodes_.size() - 1;
245 }
246
248 std::array<ctype, 2*dimworld>
249 computeBBoxOfBBoxes_(const std::vector<ctype>& leafBoxes,
250 const std::vector<std::size_t>::iterator& begin,
251 const std::vector<std::size_t>::iterator& end)
252 {
253 std::array<ctype, 2*dimworld> bBoxCoords;
254
255 // copy the iterator and get coordinates of first box
256 auto it = begin;
257 const auto* bFirst = leafBoxes.data() + 2*dimworld*(*it);
258
259 for (int coordIdx = 0; coordIdx < 2*dimworld; ++coordIdx)
260 bBoxCoords[coordIdx] = bFirst[coordIdx];
261
262 // Compute min and max with the remaining boxes
263 for (; it != end; ++it)
264 {
265 const auto* b = leafBoxes.data() + 2*dimworld*(*it);
266 for (int coordIdx = 0; coordIdx < dimworld; ++coordIdx)
267 if (b[coordIdx] < bBoxCoords[coordIdx])
268 bBoxCoords[coordIdx] = b[coordIdx];
269 for (int coordIdx = dimworld; coordIdx < 2*dimworld; ++coordIdx)
270 if (b[coordIdx] > bBoxCoords[coordIdx])
271 bBoxCoords[coordIdx] = b[coordIdx];
272 }
273
274 return bBoxCoords;
275 }
276
278 std::size_t computeLongestAxis_(const std::array<ctype, 2*dimworld>& bCoords)
279 {
280 std::array<ctype, dimworld> axisLength;
281 for (int coordIdx = 0; coordIdx < dimworld; ++coordIdx)
282 axisLength[coordIdx] = bCoords[dimworld + coordIdx] - bCoords[coordIdx];
283
284 return std::distance(axisLength.begin(), std::max_element(axisLength.begin(), axisLength.end()));
285 }
286};
287
293template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 3, int> = 0>
294inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b)
295{
296 static constexpr ctype eps_ = 1.0e-7;
297
298 using std::max;
299 const auto dx = b[3] - b[0];
300 const auto dy = b[4] - b[1];
301 const auto dz = b[5] - b[2];
302 const ctype eps = max({dx, dy, dz})*eps_;
303 return (b[0] - eps <= point[0] && point[0] <= b[3] + eps &&
304 b[1] - eps <= point[1] && point[1] <= b[4] + eps &&
305 b[2] - eps <= point[2] && point[2] <= b[5] + eps);
306}
307
313template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 2, int> = 0>
314inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b)
315{
316 static constexpr ctype eps_ = 1.0e-7;
317
318 using std::max;
319 const auto dx = b[2] - b[0];
320 const auto dy = b[3] - b[1];
321 const ctype eps = max(dx, dy)*eps_;
322 return (b[0] - eps <= point[0] && point[0] <= b[2] + eps &&
323 b[1] - eps <= point[1] && point[1] <= b[3] + eps);
324}
325
331template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 1, int> = 0>
332inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b)
333{
334 static constexpr ctype eps_ = 1.0e-7;
335 const ctype eps0 = eps_*(b[1] - b[0]);
336 return b[0] - eps0 <= point[0] && point[0] <= b[1] + eps0;
337}
338
344template<class ctype, int dimworld>
345inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point,
346 const Dune::FieldVector<ctype, dimworld>& min,
347 const Dune::FieldVector<ctype, dimworld>& max)
348{
349 std::array<ctype, 2*dimworld> bBox;
350 std::copy(min.begin(), min.end(), bBox.begin());
351 std::copy(max.begin(), max.end(), bBox.begin()+dimworld);
352 return intersectsPointBoundingBox(point, bBox.data());
353}
354
360template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 3, int> = 0>
361inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
362{
363 using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType;
364 static constexpr ctype eps_ = 1.0e-7;
365 const ctype eps0 = eps_*std::max(b[3]-b[0], a[3]-a[0]);
366 const ctype eps1 = eps_*std::max(b[4]-b[1], a[4]-a[1]);
367 const ctype eps2 = eps_*std::max(b[5]-b[2], a[5]-a[2]);
368 return (b[0] - eps0 <= a[3] && a[0] <= b[3] + eps0 &&
369 b[1] - eps1 <= a[4] && a[1] <= b[4] + eps1 &&
370 b[2] - eps2 <= a[5] && a[2] <= b[5] + eps2);
371
372}
373
379template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 2, int> = 0>
380inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
381{
382 using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType;
383 static constexpr ctype eps_ = 1.0e-7;
384 const ctype eps0 = eps_*std::max(b[2]-b[0], a[2]-a[0]);
385 const ctype eps1 = eps_*std::max(b[3]-b[1], a[3]-a[1]);
386 return (b[0] - eps0 <= a[2] && a[0] <= b[2] + eps0 &&
387 b[1] - eps1 <= a[3] && a[1] <= b[3] + eps1);
388}
389
395template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 1, int> = 0>
396inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
397{
398 using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType;
399 static constexpr ctype eps_ = 1.0e-7;
400 const ctype eps0 = eps_*std::max(b[1]-b[0], a[1]-a[0]);
401 return b[0] - eps0 <= a[1] && a[0] <= b[1] + eps0;
402}
403
409template<int dimworld, class ctype>
410inline ctype squaredDistancePointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b)
411{
412 ctype squaredDistance = 0.0;
413 for (int d = 0; d < dimworld; ++d)
414 {
415 if (point[d] < b[d])
416 squaredDistance += (point[d] - b[d])*(point[d] - b[d]);
417 if (point[d] > b[d+dimworld])
418 squaredDistance += (point[d] - b[d+dimworld])*(point[d] - b[d+dimworld]);
419 }
420 return squaredDistance;
421}
422
423} // end namespace Dumux
424
425#endif
An axis-aligned bounding box volume tree implementation.
Definition: boundingboxtree.hh:56
BoundingBoxTree()=default
Default Constructor.
const ctype * getBoundingBoxCoordinates(std::size_t nodeIdx) const
Get an existing bounding box for a given node.
Definition: boundingboxtree.hh:135
bool isLeaf(const BoundingBoxNode &node, std::size_t nodeIdx) const
Definition: boundingboxtree.hh:144
const EntitySet & entitySet() const
the entity set this tree was built with
Definition: boundingboxtree.hh:123
GeometricEntitySet EntitySet
the type of entity set this tree was built with
Definition: boundingboxtree.hh:73
void build(std::shared_ptr< const GeometricEntitySet > set)
Build up bounding box tree for a grid with leafGridView.
Definition: boundingboxtree.hh:83
BoundingBoxTree(std::shared_ptr< const GeometricEntitySet > set)
Constructor with gridView.
Definition: boundingboxtree.hh:79
std::size_t numBoundingBoxes() const
Get the number of bounding boxes currently in the tree.
Definition: boundingboxtree.hh:139
const BoundingBoxNode & getBoundingBoxNode(std::size_t nodeIdx) const
Interface to be used by other bounding box trees.
Definition: boundingboxtree.hh:131
static ctype distance(const Dune::FieldVector< ctype, dimWorld > &a, const Dune::FieldVector< ctype, dimWorld > &b)
Compute the shortest distance between two points.
Definition: distance.hh:282
static ctype squaredDistance(const Dune::FieldVector< ctype, dimWorld > &a, const Dune::FieldVector< ctype, dimWorld > &b)
Compute the shortest squared distance between two points.
Definition: distance.hh:291
constexpr BJ bj
Tag for the Beavers-Joseph slip condition.
Definition: slipcondition.hh:42
Definition: adapt.hh:17
bool intersectsBoundingBoxBoundingBox(const ctypea *a, const ctypeb *b)
Check whether a bounding box is intersecting another bounding box (dimworld == 3)
Definition: boundingboxtree.hh:361
ctype squaredDistancePointBoundingBox(const Dune::FieldVector< ctype, dimworld > &point, const ctype *b)
Compute squared distance between point and bounding box.
Definition: boundingboxtree.hh:410
bool intersectsPointBoundingBox(const Dune::FieldVector< ctype, dimworld > &point, const ctype *b)
Check whether a point is intersectin a bounding box (dimworld == 3)
Definition: boundingboxtree.hh:294