19#ifndef DUMUX_GEOMETRY_BOUNDINGBOXTREE_HH
20#define DUMUX_GEOMETRY_BOUNDINGBOXTREE_HH
30#include <dune/common/promotiontraits.hh>
31#include <dune/common/timer.hh>
32#include <dune/common/fvector.hh>
54template <
class GeometricEntitySet>
57 enum { dimworld = GeometricEntitySet::dimensionworld };
58 using ctype =
typename GeometricEntitySet::ctype;
65 struct BoundingBoxNode
83 void build(std::shared_ptr<const GeometricEntitySet> set)
89 boundingBoxNodes_.clear();
90 boundingBoxCoordinates_.clear();
96 const auto numLeaves = set->size();
99 const auto numNodes = 2*numLeaves - 1;
100 boundingBoxNodes_.reserve(numNodes);
101 boundingBoxCoordinates_.reserve(numNodes*2*dimworld);
104 std::vector<ctype> leafBoxes(2*dimworld*numLeaves);
106 for (
const auto& geometricEntity : *set)
107 computeEntityBoundingBox_(leafBoxes.data() + 2*dimworld*set->index(geometricEntity), geometricEntity);
110 std::vector<std::size_t> leafPartition(numLeaves);
111 std::iota(leafPartition.begin(), leafPartition.end(), 0);
114 build_(leafBoxes, leafPartition.begin(), leafPartition.end());
118 <<
" nodes for " << numLeaves <<
" grid entities in "
119 << timer.stop() <<
" seconds." << std::endl;
124 {
return *entitySet_; }
132 {
return boundingBoxNodes_[nodeIdx]; }
136 {
return boundingBoxCoordinates_.data() + 2*dimworld*nodeIdx; }
140 {
return boundingBoxNodes_.size(); }
144 bool isLeaf(
const BoundingBoxNode& node, std::size_t nodeIdx)
const
145 {
return node.child0 == nodeIdx; }
150 std::vector<BoundingBoxNode> boundingBoxNodes_;
153 std::vector<ctype> boundingBoxCoordinates_;
156 std::shared_ptr<const EntitySet> entitySet_;
159 template <
class Entity>
160 void computeEntityBoundingBox_(ctype* b,
const Entity& entity)
const
164 ctype* xMax = b + dimworld;
167 auto geometry = entity.geometry();
170 auto corner = geometry.corner(0);
171 for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx)
172 xMin[dimIdx] = xMax[dimIdx] = corner[dimIdx];
175 for (std::size_t cornerIdx = 1; cornerIdx < geometry.corners(); ++cornerIdx)
177 corner = geometry.corner(cornerIdx);
178 for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx)
182 xMin[dimIdx] = min(xMin[dimIdx], corner[dimIdx]);
183 xMax[dimIdx] = max(xMax[dimIdx], corner[dimIdx]);
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)
196 if (end - begin == 1)
199 const std::size_t leafNodeIdx = *begin;
200 const auto beginCoords = leafBoxes.begin() + 2*dimworld*leafNodeIdx;
201 const auto endCoords = beginCoords + 2*dimworld;
206 return addBoundingBox_(BoundingBoxNode{
numBoundingBoxes(), leafNodeIdx}, beginCoords, endCoords);
210 const auto bCoords = computeBBoxOfBBoxes_(leafBoxes, begin, end);
213 const auto axis = computeLongestAxis_(bCoords);
217 auto middle = begin + (end - begin)/2;
218 std::nth_element(begin, middle, end, [&leafBoxes, axis](std::size_t i, std::size_t j)
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];
227 return addBoundingBox_(BoundingBoxNode{build_(leafBoxes, begin, middle), build_(leafBoxes, middle, end)},
228 bCoords.begin(), bCoords.end());
232 template <
class Iterator>
233 std::size_t addBoundingBox_(BoundingBoxNode&& node,
234 const Iterator& coordBegin,
235 const Iterator& coordEnd)
238 boundingBoxNodes_.emplace_back(node);
241 boundingBoxCoordinates_.insert(boundingBoxCoordinates_.end(), coordBegin, coordEnd);
244 return boundingBoxNodes_.size() - 1;
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)
253 std::array<ctype, 2*dimworld> bBoxCoords;
257 const auto* bFirst = leafBoxes.data() + 2*dimworld*(*it);
259 for (
int coordIdx = 0; coordIdx < 2*dimworld; ++coordIdx)
260 bBoxCoords[coordIdx] = bFirst[coordIdx];
263 for (; it != end; ++it)
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];
278 std::size_t computeLongestAxis_(
const std::array<ctype, 2*dimworld>& bCoords)
280 std::array<ctype, dimworld> axisLength;
281 for (
int coordIdx = 0; coordIdx < dimworld; ++coordIdx)
282 axisLength[coordIdx] = bCoords[dimworld + coordIdx] - bCoords[coordIdx];
284 return std::distance(axisLength.begin(), std::max_element(axisLength.begin(), axisLength.end()));
293template<
class ctype,
int dimworld,
typename std::enable_if_t<dimworld == 3,
int> = 0>
296 static constexpr ctype eps_ = 1.0e-7;
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);
313template<
class ctype,
int dimworld,
typename std::enable_if_t<dimworld == 2,
int> = 0>
316 static constexpr ctype eps_ = 1.0e-7;
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);
331template<
class ctype,
int dimworld,
typename std::enable_if_t<dimworld == 1,
int> = 0>
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;
344template<
class ctype,
int dimworld>
346 const Dune::FieldVector<ctype, dimworld>& min,
347 const Dune::FieldVector<ctype, dimworld>& max)
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);
360template<
int dimworld,
class ctypea,
class ctypeb,
typename std::enable_if_t<dimworld == 3,
int> = 0>
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);
379template<
int dimworld,
class ctypea,
class ctypeb,
typename std::enable_if_t<dimworld == 2,
int> = 0>
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);
395template<
int dimworld,
class ctypea,
class ctypeb,
typename std::enable_if_t<dimworld == 1,
int> = 0>
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;
409template<
int dimworld,
class ctype>
413 for (
int d = 0; d < dimworld; ++d)
417 if (point[d] > b[d+dimworld])
418 squaredDistance += (point[d] - b[d+dimworld])*(point[d] - b[d+dimworld]);
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
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