3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
geometry/boundingboxtree.hh
Go to the documentation of this file.
1/*****************************************************************************
2 * See the file COPYING for full copying permissions. *
3 * *
4 * This program is free software: you can redistribute it and/or modify *
5 * it under the terms of the GNU General Public License as published by *
6 * the Free Software Foundation, either version 3 of the License, or *
7 * (at your option) any later version. *
8 * *
9 * This program is distributed in the hope that it will be useful, *
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
12 * GNU General Public License for more details. *
13 * *
14 * You should have received a copy of the GNU General Public License *
15 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
16 *****************************************************************************/
29#ifndef DUMUX_GEOMETRY_BOUNDINGBOXTREE_HH
30#define DUMUX_GEOMETRY_BOUNDINGBOXTREE_HH
31
32#include <vector>
33#include <array>
34#include <algorithm>
35#include <memory>
36#include <numeric>
37#include <type_traits>
38#include <iostream>
39
40#include <dune/common/promotiontraits.hh>
41#include <dune/common/timer.hh>
42#include <dune/common/fvector.hh>
43
44namespace Dumux {
45
64template <class GeometricEntitySet>
66{
67 enum { dimworld = GeometricEntitySet::dimensionworld };
68 using ctype = typename GeometricEntitySet::ctype;
69
75 struct BoundingBoxNode
76 {
77 std::size_t child0;
78 std::size_t child1;
79 };
80
81public:
83 using EntitySet = GeometricEntitySet;
84
86 BoundingBoxTree() = default;
87
89 BoundingBoxTree(std::shared_ptr<const GeometricEntitySet> set)
90 { build(set); }
91
93 void build(std::shared_ptr<const GeometricEntitySet> set)
94 {
95 // set the pointer to the entity set
96 entitySet_ = set;
97
98 // clear all internal data
99 boundingBoxNodes_.clear();
100 boundingBoxCoordinates_.clear();
101
102 // start the timer
103 Dune::Timer timer;
104
105 // Create bounding boxes for all elements
106 const auto numLeaves = set->size();
107
108 // reserve enough space for the nodes and the coordinates
109 const auto numNodes = 2*numLeaves - 1;
110 boundingBoxNodes_.reserve(numNodes);
111 boundingBoxCoordinates_.reserve(numNodes*2*dimworld);
112
113 // create a vector for leaf boxes (min and max for all dims)
114 std::vector<ctype> leafBoxes(2*dimworld*numLeaves);
115
116 for (const auto& geometricEntity : *set)
117 computeEntityBoundingBox_(leafBoxes.data() + 2*dimworld*set->index(geometricEntity), geometricEntity);
118
119 // create the leaf partition, the set of available indices (to be sorted)
120 std::vector<std::size_t> leafPartition(numLeaves);
121 std::iota(leafPartition.begin(), leafPartition.end(), 0);
122
123 // Recursively build the bounding box tree
124 build_(leafBoxes, leafPartition.begin(), leafPartition.end());
125
126 // We are done, log output
127 std::cout << "Computed bounding box tree with " << numBoundingBoxes()
128 << " nodes for " << numLeaves << " grid entites in "
129 << timer.stop() << " seconds." << std::endl;
130 }
131
133 const EntitySet& entitySet() const
134 { return *entitySet_; }
135
139
141 const BoundingBoxNode& getBoundingBoxNode(std::size_t nodeIdx) const
142 { return boundingBoxNodes_[nodeIdx]; }
143
145 const ctype* getBoundingBoxCoordinates(std::size_t nodeIdx) const
146 { return boundingBoxCoordinates_.data() + 2*dimworld*nodeIdx; }
147
149 std::size_t numBoundingBoxes() const
150 { return boundingBoxNodes_.size(); }
151
154 bool isLeaf(const BoundingBoxNode& node, std::size_t nodeIdx) const
155 { return node.child0 == nodeIdx; }
156
157private:
158
160 std::vector<BoundingBoxNode> boundingBoxNodes_;
161
163 std::vector<ctype> boundingBoxCoordinates_;
164
166 std::shared_ptr<const EntitySet> entitySet_;
167
169 template <class Entity>
170 void computeEntityBoundingBox_(ctype* b, const Entity& entity) const
171 {
172 // get the bounding box coordinates
173 ctype* xMin = b;
174 ctype* xMax = b + dimworld;
175
176 // get mesh entity data
177 auto geometry = entity.geometry();
178
179 // Get coordinates of first vertex
180 auto corner = geometry.corner(0);
181 for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx)
182 xMin[dimIdx] = xMax[dimIdx] = corner[dimIdx];
183
184 // Compute the min and max over the remaining vertices
185 for (std::size_t cornerIdx = 1; cornerIdx < geometry.corners(); ++cornerIdx)
186 {
187 corner = geometry.corner(cornerIdx);
188 for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx)
189 {
190 using std::max;
191 using std::min;
192 xMin[dimIdx] = min(xMin[dimIdx], corner[dimIdx]);
193 xMax[dimIdx] = max(xMax[dimIdx], corner[dimIdx]);
194 }
195 }
196 }
197
199 std::size_t build_(const std::vector<ctype>& leafBoxes,
200 const std::vector<std::size_t>::iterator& begin,
201 const std::vector<std::size_t>::iterator& end)
202 {
203 assert(begin < end);
204
205 // If we reached the end of the recursion, i.e. only a leaf box is left
206 if (end - begin == 1)
207 {
208 // Get the bounding box coordinates for the leaf
209 const std::size_t leafNodeIdx = *begin;
210 const auto beginCoords = leafBoxes.begin() + 2*dimworld*leafNodeIdx;
211 const auto endCoords = beginCoords + 2*dimworld;
212
213 // Store the data in the bounding box
214 // leaf nodes are indicated by setting child0 to
215 // the node itself and child1 to the index of the entity in the bounding box.
216 return addBoundingBox_(BoundingBoxNode{numBoundingBoxes(), leafNodeIdx}, beginCoords, endCoords);
217 }
218
219 // Compute the bounding box of all bounding boxes in the range [begin, end]
220 const auto bCoords = computeBBoxOfBBoxes_(leafBoxes, begin, end);
221
222 // sort bounding boxes along the longest axis
223 const auto axis = computeLongestAxis_(bCoords);
224
225 // nth_element sorts the range to make sure that middle points to the coordinate median in axis direction
226 // this is the most expensive part of the algorithm
227 auto middle = begin + (end - begin)/2;
228 std::nth_element(begin, middle, end, [&leafBoxes, axis](std::size_t i, std::size_t j)
229 {
230 const ctype* bi = leafBoxes.data() + 2*dimworld*i;
231 const ctype* bj = leafBoxes.data() + 2*dimworld*j;
232 return bi[axis] + bi[axis + dimworld] < bj[axis] + bj[axis + dimworld];
233 });
234
235 // split the bounding boxes into two at the middle iterator and call build recursively, each
236 // call resulting in a new node of this bounding box, i.e. the root will be added at the end of the process.
237 return addBoundingBox_(BoundingBoxNode{build_(leafBoxes, begin, middle), build_(leafBoxes, middle, end)},
238 bCoords.begin(), bCoords.end());
239 }
240
242 template <class Iterator>
243 std::size_t addBoundingBox_(BoundingBoxNode&& node,
244 const Iterator& coordBegin,
245 const Iterator& coordEnd)
246 {
247 // Add the bounding box
248 boundingBoxNodes_.emplace_back(node);
249
250 // Add the bounding box's coordinates
251 boundingBoxCoordinates_.insert(boundingBoxCoordinates_.end(), coordBegin, coordEnd);
252
253 // return the index of the new node
254 return boundingBoxNodes_.size() - 1;
255 }
256
258 std::array<ctype, 2*dimworld>
259 computeBBoxOfBBoxes_(const std::vector<ctype>& leafBoxes,
260 const std::vector<std::size_t>::iterator& begin,
261 const std::vector<std::size_t>::iterator& end)
262 {
263 std::array<ctype, 2*dimworld> bBoxCoords;
264
265 // copy the iterator and get coordinates of first box
266 auto it = begin;
267 const auto* bFirst = leafBoxes.data() + 2*dimworld*(*it);
268
269 for (int coordIdx = 0; coordIdx < 2*dimworld; ++coordIdx)
270 bBoxCoords[coordIdx] = bFirst[coordIdx];
271
272 // Compute min and max with the remaining boxes
273 for (; it != end; ++it)
274 {
275 const auto* b = leafBoxes.data() + 2*dimworld*(*it);
276 for (int coordIdx = 0; coordIdx < dimworld; ++coordIdx)
277 if (b[coordIdx] < bBoxCoords[coordIdx])
278 bBoxCoords[coordIdx] = b[coordIdx];
279 for (int coordIdx = dimworld; coordIdx < 2*dimworld; ++coordIdx)
280 if (b[coordIdx] > bBoxCoords[coordIdx])
281 bBoxCoords[coordIdx] = b[coordIdx];
282 }
283
284 return bBoxCoords;
285 }
286
288 std::size_t computeLongestAxis_(const std::array<ctype, 2*dimworld>& bCoords)
289 {
290 std::array<ctype, dimworld> axisLength;
291 for (int coordIdx = 0; coordIdx < dimworld; ++coordIdx)
292 axisLength[coordIdx] = bCoords[dimworld + coordIdx] - bCoords[coordIdx];
293
294 return std::distance(axisLength.begin(), std::max_element(axisLength.begin(), axisLength.end()));
295 }
296};
297
303template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 3, int> = 0>
304inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b)
305{
306 static constexpr ctype eps_ = 1.0e-7;
307
308 using std::max;
309 const auto dx = b[3] - b[0];
310 const auto dy = b[4] - b[1];
311 const auto dz = b[5] - b[2];
312 const ctype eps = max({dx, dy, dz})*eps_;
313 return (b[0] - eps <= point[0] && point[0] <= b[3] + eps &&
314 b[1] - eps <= point[1] && point[1] <= b[4] + eps &&
315 b[2] - eps <= point[2] && point[2] <= b[5] + eps);
316}
317
323template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 2, int> = 0>
324inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b)
325{
326 static constexpr ctype eps_ = 1.0e-7;
327
328 using std::max;
329 const auto dx = b[2] - b[0];
330 const auto dy = b[3] - b[1];
331 const ctype eps = max(dx, dy)*eps_;
332 return (b[0] - eps <= point[0] && point[0] <= b[2] + eps &&
333 b[1] - eps <= point[1] && point[1] <= b[3] + eps);
334}
335
341template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 1, int> = 0>
342inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b)
343{
344 static constexpr ctype eps_ = 1.0e-7;
345 const ctype eps0 = eps_*(b[1] - b[0]);
346 return b[0] - eps0 <= point[0] && point[0] <= b[1] + eps0;
347}
348
354template<class ctype, int dimworld>
355inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point,
356 const Dune::FieldVector<ctype, dimworld>& min,
357 const Dune::FieldVector<ctype, dimworld>& max)
358{
359 std::array<ctype, 2*dimworld> bBox;
360 std::copy(min.begin(), min.end(), bBox.begin());
361 std::copy(max.begin(), max.end(), bBox.begin()+dimworld);
362 return intersectsPointBoundingBox(point, bBox.data());
363}
364
370template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 3, int> = 0>
371inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
372{
373 using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType;
374 static constexpr ctype eps_ = 1.0e-7;
375 const ctype eps0 = eps_*std::max(b[3]-b[0], a[3]-a[0]);
376 const ctype eps1 = eps_*std::max(b[4]-b[1], a[4]-a[1]);
377 const ctype eps2 = eps_*std::max(b[5]-b[2], a[5]-a[2]);
378 return (b[0] - eps0 <= a[3] && a[0] <= b[3] + eps0 &&
379 b[1] - eps1 <= a[4] && a[1] <= b[4] + eps1 &&
380 b[2] - eps2 <= a[5] && a[2] <= b[5] + eps2);
381
382}
383
389template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 2, int> = 0>
390inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
391{
392 using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType;
393 static constexpr ctype eps_ = 1.0e-7;
394 const ctype eps0 = eps_*std::max(b[2]-b[0], a[2]-a[0]);
395 const ctype eps1 = eps_*std::max(b[3]-b[1], a[3]-a[1]);
396 return (b[0] - eps0 <= a[2] && a[0] <= b[2] + eps0 &&
397 b[1] - eps1 <= a[3] && a[1] <= b[3] + eps1);
398}
399
405template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 1, int> = 0>
406inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
407{
408 using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType;
409 static constexpr ctype eps_ = 1.0e-7;
410 const ctype eps0 = eps_*std::max(b[1]-b[0], a[1]-a[0]);
411 return b[0] - eps0 <= a[1] && a[0] <= b[1] + eps0;
412}
413
414} // end namespace Dumux
415
416#endif
ctype distance(const Dune::FieldVector< ctype, dimWorld > &a, const Dune::FieldVector< ctype, dimWorld > &b)
Compute the shortest distance between two points.
Definition: geometry/distance.hh:138
Definition: adapt.hh:29
bool intersectsBoundingBoxBoundingBox(const ctypea *a, const ctypeb *b)
Check whether a bounding box is intersecting another bounding box (dimworld == 3)
Definition: geometry/boundingboxtree.hh:371
bool intersectsPointBoundingBox(const Dune::FieldVector< ctype, dimworld > &point, const ctype *b)
Check whether a point is intersectin a bounding box (dimworld == 3)
Definition: geometry/boundingboxtree.hh:304
An axis-aligned bounding box volume tree implementation.
Definition: geometry/boundingboxtree.hh:66
BoundingBoxTree()=default
Default Constructor.
const ctype * getBoundingBoxCoordinates(std::size_t nodeIdx) const
Get an existing bounding box for a given node.
Definition: geometry/boundingboxtree.hh:145
bool isLeaf(const BoundingBoxNode &node, std::size_t nodeIdx) const
Definition: geometry/boundingboxtree.hh:154
const EntitySet & entitySet() const
the entity set this tree was built with
Definition: geometry/boundingboxtree.hh:133
GeometricEntitySet EntitySet
the type of entity set this tree was built with
Definition: geometry/boundingboxtree.hh:83
void build(std::shared_ptr< const GeometricEntitySet > set)
Build up bounding box tree for a grid with leafGridView.
Definition: geometry/boundingboxtree.hh:93
BoundingBoxTree(std::shared_ptr< const GeometricEntitySet > set)
Constructor with gridView.
Definition: geometry/boundingboxtree.hh:89
std::size_t numBoundingBoxes() const
Get the number of bounding boxes currently in the tree.
Definition: geometry/boundingboxtree.hh:149
const BoundingBoxNode & getBoundingBoxNode(std::size_t nodeIdx) const
Interface to be used by other bounding box trees.
Definition: geometry/boundingboxtree.hh:141