| Line | Branch | Exec | Source |
|---|---|---|---|
| 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 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup Geometry | ||
| 10 | * \brief An axis-aligned bounding box volume hierarchy for dune grids | ||
| 11 | * | ||
| 12 | * Dumux implementation of an AABB tree | ||
| 13 | * Inspired by the AABB tree implementation in DOLFIN by Anders Logg which has the | ||
| 14 | * following license info: DOLFIN is free software: you can redistribute it and/or modify | ||
| 15 | * it under the terms of the GNU Lesser General Public License as published by | ||
| 16 | * the Free Software Foundation, either version 3 of the License, or | ||
| 17 | * (at your option) any later version. | ||
| 18 | */ | ||
| 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 | |||
| 34 | namespace Dumux { | ||
| 35 | |||
| 36 | /*! | ||
| 37 | * \ingroup Geometry | ||
| 38 | * \brief An axis-aligned bounding box volume tree implementation | ||
| 39 | * | ||
| 40 | * The class constructs a hierarchical structure of bounding box volumes around | ||
| 41 | * grid entities. This class can be used to efficiently compute intersections | ||
| 42 | * between a grid and other geometrical object. It only implements the intersection | ||
| 43 | * of two of such bounding box trees, so that two independent grids can be intersected. | ||
| 44 | * \tparam GeometricEntitySet has the following requirements | ||
| 45 | * * export dimensionworld, ctype | ||
| 46 | * * a size() member function returning the number of entities | ||
| 47 | * * begin() and end() member function returning at least forward iterators to entities | ||
| 48 | * * an index() method returning a consecutive index given an entity | ||
| 49 | * * an entity() method returning an entity given the consecutive index | ||
| 50 | * * entities have the following requirements: | ||
| 51 | * * a member function geometry() returning a geometry with the member functions | ||
| 52 | * * corner() and corners() returning global coordinates and number of corners | ||
| 53 | */ | ||
| 54 | template <class GeometricEntitySet> | ||
| 55 | class BoundingBoxTree | ||
| 56 | { | ||
| 57 | enum { dimworld = GeometricEntitySet::dimensionworld }; | ||
| 58 | using ctype = typename GeometricEntitySet::ctype; | ||
| 59 | |||
| 60 | /*! | ||
| 61 | * \brief Bounding box node data structure | ||
| 62 | * leaf nodes are indicated by setting child0 to | ||
| 63 | * the node itself and child1 to the index of the entity in the bounding box. | ||
| 64 | */ | ||
| 65 | struct BoundingBoxNode | ||
| 66 | { | ||
| 67 | std::size_t child0; | ||
| 68 | std::size_t child1; | ||
| 69 | }; | ||
| 70 | |||
| 71 | public: | ||
| 72 | //! the type of entity set this tree was built with | ||
| 73 | using EntitySet = GeometricEntitySet; | ||
| 74 | |||
| 75 | //! Default Constructor | ||
| 76 |
8/16✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 4 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 4 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 4 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 4 times.
✗ Branch 23 not taken.
|
72 | BoundingBoxTree() = default; |
| 77 | |||
| 78 | //! Constructor with gridView | ||
| 79 | 4615 | BoundingBoxTree(std::shared_ptr<const GeometricEntitySet> set) | |
| 80 |
2/12✓ Branch 5 taken 2560 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2560 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
|
23075 | { build(set); } |
| 81 | |||
| 82 | //! Build up bounding box tree for a grid with leafGridView | ||
| 83 | 4687 | void build(std::shared_ptr<const GeometricEntitySet> set) | |
| 84 | { | ||
| 85 | // set the pointer to the entity set | ||
| 86 | 4687 | entitySet_ = set; | |
| 87 | |||
| 88 | // clear all internal data | ||
| 89 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2594 times.
|
4687 | boundingBoxNodes_.clear(); |
| 90 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2594 times.
|
4687 | boundingBoxCoordinates_.clear(); |
| 91 | |||
| 92 | // start the timer | ||
| 93 | 4687 | Dune::Timer timer; | |
| 94 | |||
| 95 | // Create bounding boxes for all elements | ||
| 96 | 9374 | const auto numLeaves = set->size(); | |
| 97 | |||
| 98 | // reserve enough space for the nodes and the coordinates | ||
| 99 | 4687 | const auto numNodes = 2*numLeaves - 1; | |
| 100 | 4687 | boundingBoxNodes_.reserve(numNodes); | |
| 101 | 4687 | boundingBoxCoordinates_.reserve(numNodes*2*dimworld); | |
| 102 | |||
| 103 | // create a vector for leaf boxes (min and max for all dims) | ||
| 104 | 9374 | std::vector<ctype> leafBoxes(2*dimworld*numLeaves); | |
| 105 | |||
| 106 |
16/18✓ Branch 0 taken 56259 times.
✓ Branch 1 taken 2432 times.
✓ Branch 2 taken 56259 times.
✓ Branch 3 taken 125396 times.
✓ Branch 4 taken 887 times.
✓ Branch 5 taken 123685 times.
✓ Branch 6 taken 166 times.
✓ Branch 7 taken 730 times.
✓ Branch 8 taken 3751 times.
✓ Branch 9 taken 3271465 times.
✓ Branch 10 taken 224 times.
✓ Branch 11 taken 764451 times.
✓ Branch 12 taken 224 times.
✓ Branch 13 taken 11160 times.
✓ Branch 14 taken 764451 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 764451 times.
✗ Branch 18 not taken.
|
9823576 | for (const auto& geometricEntity : *set) |
| 107 |
4/8✓ Branch 1 taken 781540 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 860856 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 781540 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 770380 times.
✗ Branch 11 not taken.
|
16481261 | computeEntityBoundingBox_(leafBoxes.data() + 2*dimworld*set->index(geometricEntity), geometricEntity); |
| 108 | |||
| 109 | // create the leaf partition, the set of available indices (to be sorted) | ||
| 110 |
3/8✓ Branch 1 taken 2598 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2598 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2598 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
14061 | std::vector<std::size_t> leafPartition(numLeaves); |
| 111 | 14061 | std::iota(leafPartition.begin(), leafPartition.end(), 0); | |
| 112 | |||
| 113 | // Recursively build the bounding box tree | ||
| 114 |
3/8✓ Branch 1 taken 2598 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2598 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2598 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
14061 | build_(leafBoxes, leafPartition.begin(), leafPartition.end()); |
| 115 | |||
| 116 | // We are done, log output | ||
| 117 |
1/2✓ Branch 2 taken 2598 times.
✗ Branch 3 not taken.
|
9374 | std::cout << "Computed bounding box tree with " << numBoundingBoxes() |
| 118 |
1/2✓ Branch 2 taken 2598 times.
✗ Branch 3 not taken.
|
9374 | << " nodes for " << numLeaves << " grid entities in " |
| 119 |
3/6✓ Branch 2 taken 2598 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 2598 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2598 times.
✗ Branch 9 not taken.
|
9374 | << timer.stop() << " seconds." << std::endl; |
| 120 | 4687 | } | |
| 121 | |||
| 122 | //! the entity set this tree was built with | ||
| 123 | const EntitySet& entitySet() const | ||
| 124 |
30/43✓ Branch 1 taken 16404 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 453 times.
✓ Branch 4 taken 16404 times.
✓ Branch 5 taken 314256 times.
✓ Branch 6 taken 453 times.
✓ Branch 7 taken 11396 times.
✓ Branch 8 taken 314256 times.
✓ Branch 9 taken 25831 times.
✓ Branch 10 taken 11396 times.
✓ Branch 11 taken 316422 times.
✓ Branch 12 taken 25831 times.
✓ Branch 13 taken 11408 times.
✓ Branch 14 taken 316422 times.
✓ Branch 15 taken 4640 times.
✓ Branch 16 taken 11408 times.
✓ Branch 17 taken 314662 times.
✓ Branch 18 taken 4640 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 314662 times.
✓ Branch 21 taken 1968 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 779344 times.
✓ Branch 24 taken 1968 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 779344 times.
✗ Branch 27 not taken.
✓ Branch 29 taken 7056 times.
✗ Branch 30 not taken.
✓ Branch 32 taken 7056 times.
✗ Branch 33 not taken.
✓ Branch 35 taken 7056 times.
✗ Branch 36 not taken.
✓ Branch 38 taken 7056 times.
✗ Branch 39 not taken.
✓ Branch 41 taken 7056 times.
✗ Branch 42 not taken.
✓ Branch 44 taken 7056 times.
✗ Branch 45 not taken.
✓ Branch 47 taken 28224 times.
✗ Branch 48 not taken.
✓ Branch 50 taken 28224 times.
✗ Branch 51 not taken.
|
13190792 | { return *entitySet_; } |
| 125 | |||
| 126 | ///////////////////////////////////////////////////// | ||
| 127 | //! Interface to be used by other bounding box trees | ||
| 128 | ///////////////////////////////////////////////////// | ||
| 129 | |||
| 130 | //! Get an existing bounding box for a given node | ||
| 131 | const BoundingBoxNode& getBoundingBoxNode(std::size_t nodeIdx) const | ||
| 132 |
48/48✓ Branch 0 taken 2782898 times.
✓ Branch 1 taken 2170638 times.
✓ Branch 2 taken 2782898 times.
✓ Branch 3 taken 2170638 times.
✓ Branch 4 taken 1706105 times.
✓ Branch 5 taken 1229640 times.
✓ Branch 6 taken 1706105 times.
✓ Branch 7 taken 1229640 times.
✓ Branch 8 taken 5640330 times.
✓ Branch 9 taken 5541430 times.
✓ Branch 10 taken 5640330 times.
✓ Branch 11 taken 5541430 times.
✓ Branch 12 taken 198637 times.
✓ Branch 13 taken 265226 times.
✓ Branch 14 taken 198637 times.
✓ Branch 15 taken 265226 times.
✓ Branch 16 taken 1807792 times.
✓ Branch 17 taken 4359807 times.
✓ Branch 18 taken 1807792 times.
✓ Branch 19 taken 4359807 times.
✓ Branch 20 taken 926872 times.
✓ Branch 21 taken 2381898 times.
✓ Branch 22 taken 926872 times.
✓ Branch 23 taken 2381898 times.
✓ Branch 24 taken 2888 times.
✓ Branch 25 taken 1331 times.
✓ Branch 26 taken 2888 times.
✓ Branch 27 taken 1331 times.
✓ Branch 28 taken 2136 times.
✓ Branch 29 taken 1089 times.
✓ Branch 30 taken 2136 times.
✓ Branch 31 taken 1089 times.
✓ Branch 32 taken 517680 times.
✓ Branch 33 taken 204669 times.
✓ Branch 34 taken 517680 times.
✓ Branch 35 taken 204669 times.
✓ Branch 36 taken 496500 times.
✓ Branch 37 taken 191801 times.
✓ Branch 38 taken 496500 times.
✓ Branch 39 taken 191801 times.
✓ Branch 40 taken 1259898 times.
✓ Branch 41 taken 692794 times.
✓ Branch 42 taken 1259898 times.
✓ Branch 43 taken 692794 times.
✓ Branch 44 taken 893742 times.
✓ Branch 45 taken 477682 times.
✓ Branch 46 taken 893742 times.
✓ Branch 47 taken 477682 times.
|
67506966 | { return boundingBoxNodes_[nodeIdx]; } |
| 133 | |||
| 134 | //! Get an existing bounding box for a given node | ||
| 135 | const ctype* getBoundingBoxCoordinates(std::size_t nodeIdx) const | ||
| 136 |
48/48✓ Branch 0 taken 2786613 times.
✓ Branch 1 taken 2088517 times.
✓ Branch 2 taken 2786613 times.
✓ Branch 3 taken 2088517 times.
✓ Branch 4 taken 1706359 times.
✓ Branch 5 taken 1227652 times.
✓ Branch 6 taken 1706359 times.
✓ Branch 7 taken 1227652 times.
✓ Branch 8 taken 5743235 times.
✓ Branch 9 taken 5203181 times.
✓ Branch 10 taken 5743235 times.
✓ Branch 11 taken 5203181 times.
✓ Branch 12 taken 247699 times.
✓ Branch 13 taken 107377 times.
✓ Branch 14 taken 247699 times.
✓ Branch 15 taken 107377 times.
✓ Branch 16 taken 3085033 times.
✓ Branch 17 taken 1311208 times.
✓ Branch 18 taken 3085033 times.
✓ Branch 19 taken 1311208 times.
✓ Branch 20 taken 1653642 times.
✓ Branch 21 taken 742171 times.
✓ Branch 22 taken 1653642 times.
✓ Branch 23 taken 742171 times.
✓ Branch 24 taken 2888 times.
✓ Branch 25 taken 1331 times.
✓ Branch 26 taken 2888 times.
✓ Branch 27 taken 1331 times.
✓ Branch 28 taken 2136 times.
✓ Branch 29 taken 1089 times.
✓ Branch 30 taken 2136 times.
✓ Branch 31 taken 1089 times.
✓ Branch 32 taken 517680 times.
✓ Branch 33 taken 204669 times.
✓ Branch 34 taken 517680 times.
✓ Branch 35 taken 204669 times.
✓ Branch 36 taken 496500 times.
✓ Branch 37 taken 191801 times.
✓ Branch 38 taken 496500 times.
✓ Branch 39 taken 191801 times.
✓ Branch 40 taken 1268158 times.
✓ Branch 41 taken 688442 times.
✓ Branch 42 taken 1268158 times.
✓ Branch 43 taken 688442 times.
✓ Branch 44 taken 893742 times.
✓ Branch 45 taken 477682 times.
✓ Branch 46 taken 893742 times.
✓ Branch 47 taken 477682 times.
|
61297610 | { return boundingBoxCoordinates_.data() + 2*dimworld*nodeIdx; } |
| 137 | |||
| 138 | //! Get the number of bounding boxes currently in the tree | ||
| 139 | std::size_t numBoundingBoxes() const | ||
| 140 |
35/50✓ Branch 1 taken 774 times.
✓ Branch 2 taken 120 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 911 times.
✓ Branch 5 taken 170 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 59308 times.
✓ Branch 8 taken 50 times.
✓ Branch 9 taken 120 times.
✓ Branch 10 taken 4789 times.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 120 times.
✓ Branch 13 taken 277359 times.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 124 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 2 times.
✓ Branch 19 taken 142 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 2 times.
✓ Branch 22 taken 4 times.
✓ Branch 23 taken 120 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 10 times.
✓ Branch 26 taken 120 times.
✗ Branch 27 not taken.
✓ Branch 28 taken 22 times.
✗ Branch 29 not taken.
✓ Branch 30 taken 120 times.
✓ Branch 31 taken 8 times.
✓ Branch 32 taken 6 times.
✓ Branch 33 taken 120 times.
✓ Branch 34 taken 4 times.
✓ Branch 35 taken 6 times.
✗ Branch 36 not taken.
✓ Branch 37 taken 124 times.
✗ Branch 38 not taken.
✓ Branch 39 taken 2 times.
✓ Branch 40 taken 124 times.
✗ Branch 41 not taken.
✓ Branch 42 taken 2 times.
✓ Branch 43 taken 16 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 18 times.
✗ Branch 47 not taken.
✓ Branch 49 taken 6 times.
✗ Branch 50 not taken.
✓ Branch 52 taken 8 times.
✗ Branch 53 not taken.
|
414411 | { return boundingBoxNodes_.size(); } |
| 141 | |||
| 142 | //! Check whether a bounding box node is a leaf node | ||
| 143 | //! Leaf nodes have itself as child0 | ||
| 144 | ✗ | bool isLeaf(const BoundingBoxNode& node, std::size_t nodeIdx) const | |
| 145 |
20/20✓ Branch 0 taken 236281 times.
✓ Branch 1 taken 1413134 times.
✓ Branch 2 taken 236281 times.
✓ Branch 3 taken 1413134 times.
✓ Branch 4 taken 17267 times.
✓ Branch 5 taken 72536 times.
✓ Branch 6 taken 17267 times.
✓ Branch 7 taken 72536 times.
✓ Branch 8 taken 4580 times.
✓ Branch 9 taken 9335 times.
✓ Branch 10 taken 4580 times.
✓ Branch 11 taken 9335 times.
✓ Branch 12 taken 256 times.
✓ Branch 13 taken 308 times.
✓ Branch 14 taken 256 times.
✓ Branch 15 taken 308 times.
✓ Branch 16 taken 143776 times.
✓ Branch 17 taken 287492 times.
✓ Branch 18 taken 143776 times.
✓ Branch 19 taken 287492 times.
|
4369930 | { return node.child0 == nodeIdx; } |
| 146 | |||
| 147 | private: | ||
| 148 | |||
| 149 | //! vector of bounding box nodes | ||
| 150 | std::vector<BoundingBoxNode> boundingBoxNodes_; | ||
| 151 | |||
| 152 | //! vector of bounding box coordinates | ||
| 153 | std::vector<ctype> boundingBoxCoordinates_; | ||
| 154 | |||
| 155 | //! a pointer to the entity set | ||
| 156 | std::shared_ptr<const EntitySet> entitySet_; | ||
| 157 | |||
| 158 | //! Compute the bounding box of a grid entity | ||
| 159 | template <class Entity> | ||
| 160 | 5417427 | void computeEntityBoundingBox_(ctype* b, const Entity& entity) const | |
| 161 | { | ||
| 162 | // get the bounding box coordinates | ||
| 163 | 5445171 | ctype* xMin = b; | |
| 164 | 5445171 | ctype* xMax = b + dimworld; | |
| 165 | |||
| 166 | // get mesh entity data | ||
| 167 |
1/2✓ Branch 2 taken 2178 times.
✗ Branch 3 not taken.
|
5988533 | auto geometry = entity.geometry(); |
| 168 | |||
| 169 | // Get coordinates of first vertex | ||
| 170 | 5445171 | auto corner = geometry.corner(0); | |
| 171 |
12/12✓ Branch 0 taken 3629827 times.
✓ Branch 1 taken 7768082 times.
✓ Branch 2 taken 216546 times.
✓ Branch 3 taken 528 times.
✓ Branch 4 taken 24000 times.
✓ Branch 5 taken 8160 times.
✓ Branch 6 taken 7200 times.
✓ Branch 7 taken 2400 times.
✓ Branch 8 taken 34020 times.
✓ Branch 9 taken 11340 times.
✓ Branch 10 taken 11820 times.
✓ Branch 11 taken 4100 times.
|
18964092 | for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx) |
| 172 | 27023150 | xMin[dimIdx] = xMax[dimIdx] = corner[dimIdx]; | |
| 173 | |||
| 174 | // Compute the min and max over the remaining vertices | ||
| 175 |
4/4✓ Branch 0 taken 14422958 times.
✓ Branch 1 taken 3423159 times.
✓ Branch 2 taken 14422891 times.
✓ Branch 3 taken 3423126 times.
|
35800534 | for (std::size_t cornerIdx = 1; cornerIdx < geometry.corners(); ++cornerIdx) |
| 176 | { | ||
| 177 | 24923310 | corner = geometry.corner(cornerIdx); | |
| 178 |
3/3✓ Branch 0 taken 37110668 times.
✓ Branch 1 taken 14551260 times.
✓ Branch 2 taken 125558 times.
|
91617509 | for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx) |
| 179 | { | ||
| 180 | using std::max; | ||
| 181 | using std::min; | ||
| 182 |
4/4✓ Branch 0 taken 1219247 times.
✓ Branch 1 taken 36145281 times.
✓ Branch 2 taken 1219247 times.
✓ Branch 3 taken 36130609 times.
|
133373706 | xMin[dimIdx] = min(xMin[dimIdx], corner[dimIdx]); |
| 183 |
4/4✓ Branch 0 taken 8075293 times.
✓ Branch 1 taken 29289235 times.
✓ Branch 2 taken 8060621 times.
✓ Branch 3 taken 29289235 times.
|
146759075 | xMax[dimIdx] = max(xMax[dimIdx], corner[dimIdx]); |
| 184 | } | ||
| 185 | } | ||
| 186 | 5417427 | } | |
| 187 | |||
| 188 | //! Build bounding box tree for all entities recursively | ||
| 189 | 10941143 | 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 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 6899208 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6899208 times.
|
21882286 | assert(begin < end); |
| 194 | |||
| 195 | // If we reached the end of the recursion, i.e. only a leaf box is left | ||
| 196 |
4/4✓ Branch 0 taken 3450903 times.
✓ Branch 1 taken 3448305 times.
✓ Branch 2 taken 3450903 times.
✓ Branch 3 taken 3448305 times.
|
21882286 | if (end - begin == 1) |
| 197 | { | ||
| 198 | // Get the bounding box coordinates for the leaf | ||
| 199 | 5472915 | const std::size_t leafNodeIdx = *begin; | |
| 200 | 10945830 | const auto beginCoords = leafBoxes.begin() + 2*dimworld*leafNodeIdx; | |
| 201 | 5472915 | 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 | 10945830 | return addBoundingBox_(BoundingBoxNode{numBoundingBoxes(), leafNodeIdx}, beginCoords, endCoords); | |
| 207 | } | ||
| 208 | |||
| 209 | // Compute the bounding box of all bounding boxes in the range [begin, end] | ||
| 210 | 5468228 | const auto bCoords = computeBBoxOfBBoxes_(leafBoxes, begin, end); | |
| 211 | |||
| 212 | // sort bounding boxes along the longest axis | ||
| 213 | 5482873 | 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 | 5468228 | auto middle = begin + (end - begin)/2; | |
| 218 | 174470186 | std::nth_element(begin, middle, end, [&leafBoxes, axis](std::size_t i, std::size_t j) | |
| 219 | { | ||
| 220 | 274777664 | const ctype* bi = leafBoxes.data() + 2*dimworld*i; | |
| 221 | 274777664 | const ctype* bj = leafBoxes.data() + 2*dimworld*j; | |
| 222 | 137388832 | 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 |
5/10✓ Branch 1 taken 3448305 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3448305 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3448305 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3448305 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3448305 times.
✗ Branch 14 not taken.
|
16404684 | bCoords.begin(), bCoords.end()); |
| 229 | } | ||
| 230 | |||
| 231 | //! Add a new bounding box to the tree | ||
| 232 | template <class Iterator> | ||
| 233 | 6896610 | std::size_t addBoundingBox_(BoundingBoxNode&& node, | |
| 234 | const Iterator& coordBegin, | ||
| 235 | const Iterator& coordEnd) | ||
| 236 | { | ||
| 237 | // Add the bounding box | ||
| 238 | 6896610 | boundingBoxNodes_.emplace_back(node); | |
| 239 | |||
| 240 | // Add the bounding box's coordinates | ||
| 241 |
0/6✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
20689830 | boundingBoxCoordinates_.insert(boundingBoxCoordinates_.end(), coordBegin, coordEnd); |
| 242 | |||
| 243 | // return the index of the new node | ||
| 244 | 13793220 | return boundingBoxNodes_.size() - 1; | |
| 245 | } | ||
| 246 | |||
| 247 | //! Compute the bounding box of a vector of bounding boxes | ||
| 248 | std::array<ctype, 2*dimworld> | ||
| 249 | 5468228 | 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 | 5468228 | auto it = begin; | |
| 257 | 5468228 | const auto* bFirst = leafBoxes.data() + 2*dimworld*(*it); | |
| 258 | |||
| 259 |
2/2✓ Branch 0 taken 3448305 times.
✓ Branch 1 taken 16521322 times.
|
32645282 | for (int coordIdx = 0; coordIdx < 2*dimworld; ++coordIdx) |
| 260 | 54354108 | bBoxCoords[coordIdx] = bFirst[coordIdx]; | |
| 261 | |||
| 262 | // Compute min and max with the remaining boxes | ||
| 263 |
4/4✓ Branch 0 taken 50188143 times.
✓ Branch 1 taken 3448305 times.
✓ Branch 2 taken 50188143 times.
✓ Branch 3 taken 3448305 times.
|
91744439 | for (; it != end; ++it) |
| 264 | { | ||
| 265 | 80807983 | const auto* b = leafBoxes.data() + 2*dimworld*(*it); | |
| 266 |
2/2✓ Branch 0 taken 50188143 times.
✓ Branch 1 taken 122059830 times.
|
284970374 | for (int coordIdx = 0; coordIdx < dimworld; ++coordIdx) |
| 267 |
4/4✓ Branch 0 taken 4061816 times.
✓ Branch 1 taken 117998014 times.
✓ Branch 2 taken 4061816 times.
✓ Branch 3 taken 117998014 times.
|
408050958 | if (b[coordIdx] < bBoxCoords[coordIdx]) |
| 268 | 12573600 | bBoxCoords[coordIdx] = b[coordIdx]; | |
| 269 |
2/2✓ Branch 0 taken 122059830 times.
✓ Branch 1 taken 50188143 times.
|
284833462 | for (int coordIdx = dimworld; coordIdx < 2*dimworld; ++coordIdx) |
| 270 |
4/4✓ Branch 0 taken 6466075 times.
✓ Branch 1 taken 115593755 times.
✓ Branch 2 taken 6466075 times.
✓ Branch 3 taken 115593755 times.
|
408050958 | if (b[coordIdx] > bBoxCoords[coordIdx]) |
| 271 | 20237312 | bBoxCoords[coordIdx] = b[coordIdx]; | |
| 272 | } | ||
| 273 | |||
| 274 | 5468228 | return bBoxCoords; | |
| 275 | } | ||
| 276 | |||
| 277 | //! Compute the bounding box of a vector of bounding boxes | ||
| 278 | ✗ | std::size_t computeLongestAxis_(const std::array<ctype, 2*dimworld>& bCoords) | |
| 279 | { | ||
| 280 | std::array<ctype, dimworld> axisLength; | ||
| 281 |
30/30✓ Branch 0 taken 423600 times.
✓ Branch 1 taken 905465 times.
✓ Branch 2 taken 2777384 times.
✓ Branch 3 taken 6649153 times.
✓ Branch 4 taken 20870 times.
✓ Branch 5 taken 62016 times.
✓ Branch 6 taken 160538 times.
✓ Branch 7 taken 477168 times.
✓ Branch 8 taken 498 times.
✓ Branch 9 taken 1040 times.
✓ Branch 10 taken 454 times.
✓ Branch 11 taken 956 times.
✓ Branch 12 taken 406 times.
✓ Branch 13 taken 1172 times.
✓ Branch 14 taken 8034 times.
✓ Branch 15 taken 24102 times.
✓ Branch 16 taken 8034 times.
✓ Branch 17 taken 24102 times.
✓ Branch 18 taken 2274 times.
✓ Branch 19 taken 6822 times.
✓ Branch 20 taken 2274 times.
✓ Branch 21 taken 6462 times.
✓ Branch 22 taken 11218 times.
✓ Branch 23 taken 33294 times.
✓ Branch 24 taken 10858 times.
✓ Branch 25 taken 32574 times.
✓ Branch 26 taken 3618 times.
✓ Branch 27 taken 10854 times.
✓ Branch 28 taken 3618 times.
✓ Branch 29 taken 10854 times.
|
11694339 | for (int coordIdx = 0; coordIdx < dimworld; ++coordIdx) |
| 282 | 32984136 | axisLength[coordIdx] = bCoords[dimworld + coordIdx] - bCoords[coordIdx]; | |
| 283 | |||
| 284 | 10301034 | return std::distance(axisLength.begin(), std::max_element(axisLength.begin(), axisLength.end())); | |
| 285 | } | ||
| 286 | }; | ||
| 287 | |||
| 288 | /*! | ||
| 289 | * \brief Check whether a point is intersectin a bounding box (dimworld == 3) | ||
| 290 | * \param point The point | ||
| 291 | * \param b Pointer to bounding box coordinates | ||
| 292 | */ | ||
| 293 | template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 3, int> = 0> | ||
| 294 | 25049136 | inline 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 | 25049136 | const auto dx = b[3] - b[0]; | |
| 300 | 25049136 | const auto dy = b[4] - b[1]; | |
| 301 | 25049136 | const auto dz = b[5] - b[2]; | |
| 302 | 25049136 | const ctype eps = max({dx, dy, dz})*eps_; | |
| 303 |
6/6✓ Branch 0 taken 23899873 times.
✓ Branch 1 taken 1149263 times.
✓ Branch 2 taken 22741721 times.
✓ Branch 3 taken 1158152 times.
✓ Branch 4 taken 22741721 times.
✓ Branch 5 taken 1158152 times.
|
25049136 | return (b[0] - eps <= point[0] && point[0] <= b[3] + eps && |
| 304 |
8/8✓ Branch 0 taken 21787733 times.
✓ Branch 1 taken 953988 times.
✓ Branch 2 taken 21787733 times.
✓ Branch 3 taken 953988 times.
✓ Branch 4 taken 20805131 times.
✓ Branch 5 taken 982602 times.
✓ Branch 6 taken 20805131 times.
✓ Branch 7 taken 982602 times.
|
45483442 | b[1] - eps <= point[1] && point[1] <= b[4] + eps && |
| 305 |
10/10✓ Branch 0 taken 23899873 times.
✓ Branch 1 taken 1149263 times.
✓ Branch 2 taken 19998866 times.
✓ Branch 3 taken 806265 times.
✓ Branch 4 taken 19998866 times.
✓ Branch 5 taken 806265 times.
✓ Branch 6 taken 874072 times.
✓ Branch 7 taken 19124794 times.
✓ Branch 8 taken 874072 times.
✓ Branch 9 taken 19124794 times.
|
66659398 | b[2] - eps <= point[2] && point[2] <= b[5] + eps); |
| 306 | } | ||
| 307 | |||
| 308 | /*! | ||
| 309 | * \brief Check whether a point is intersectin a bounding box (dimworld == 2) | ||
| 310 | * \param point The point | ||
| 311 | * \param b Pointer to bounding box coordinates | ||
| 312 | */ | ||
| 313 | template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 2, int> = 0> | ||
| 314 | 757997 | inline 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 | 757997 | const auto dx = b[2] - b[0]; | |
| 320 | 757997 | const auto dy = b[3] - b[1]; | |
| 321 |
2/2✓ Branch 0 taken 370536 times.
✓ Branch 1 taken 387461 times.
|
757997 | const ctype eps = max(dx, dy)*eps_; |
| 322 |
6/6✓ Branch 0 taken 659671 times.
✓ Branch 1 taken 98326 times.
✓ Branch 2 taken 560653 times.
✓ Branch 3 taken 99018 times.
✓ Branch 4 taken 560653 times.
✓ Branch 5 taken 99018 times.
|
757997 | return (b[0] - eps <= point[0] && point[0] <= b[2] + eps && |
| 323 |
10/10✓ Branch 0 taken 659671 times.
✓ Branch 1 taken 98326 times.
✓ Branch 2 taken 462236 times.
✓ Branch 3 taken 98417 times.
✓ Branch 4 taken 462236 times.
✓ Branch 5 taken 98417 times.
✓ Branch 6 taken 105636 times.
✓ Branch 7 taken 356600 times.
✓ Branch 8 taken 105636 times.
✓ Branch 9 taken 356600 times.
|
1879303 | b[1] - eps <= point[1] && point[1] <= b[3] + eps); |
| 324 | } | ||
| 325 | |||
| 326 | /*! | ||
| 327 | * \brief Check whether a point is intersectin a bounding box (dimworld == 1) | ||
| 328 | * \param point The point | ||
| 329 | * \param b Pointer to bounding box coordinates | ||
| 330 | */ | ||
| 331 | template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 1, int> = 0> | ||
| 332 | inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b) | ||
| 333 | { | ||
| 334 | static constexpr ctype eps_ = 1.0e-7; | ||
| 335 | 214 | const ctype eps0 = eps_*(b[1] - b[0]); | |
| 336 |
4/4✓ Branch 0 taken 128 times.
✓ Branch 1 taken 86 times.
✓ Branch 2 taken 104 times.
✓ Branch 3 taken 24 times.
|
214 | return b[0] - eps0 <= point[0] && point[0] <= b[1] + eps0; |
| 337 | } | ||
| 338 | |||
| 339 | /*! | ||
| 340 | * \ingroup Geometry | ||
| 341 | * \brief Determine if a point intersects an axis-aligned bounding box | ||
| 342 | * The bounding box is given by the lower left corner (min) and the upper right corner (max) | ||
| 343 | */ | ||
| 344 | template<class ctype, int dimworld> | ||
| 345 | 12767952 | inline 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 | 51071808 | std::copy(min.begin(), min.end(), bBox.begin()); | |
| 351 | 51071808 | std::copy(max.begin(), max.end(), bBox.begin()+dimworld); | |
| 352 | 25535904 | return intersectsPointBoundingBox(point, bBox.data()); | |
| 353 | } | ||
| 354 | |||
| 355 | /*! | ||
| 356 | * \brief Check whether a bounding box is intersecting another bounding box (dimworld == 3) | ||
| 357 | * \param a Pointer to first bounding box coordinates | ||
| 358 | * \param b Pointer to second bounding box coordinates | ||
| 359 | */ | ||
| 360 | template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 3, int> = 0> | ||
| 361 | 1444586 | inline 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 |
2/2✓ Branch 0 taken 645559 times.
✓ Branch 1 taken 799027 times.
|
1444586 | const ctype eps0 = eps_*std::max(b[3]-b[0], a[3]-a[0]); |
| 366 |
2/2✓ Branch 0 taken 644464 times.
✓ Branch 1 taken 800122 times.
|
1444586 | const ctype eps1 = eps_*std::max(b[4]-b[1], a[4]-a[1]); |
| 367 |
2/2✓ Branch 0 taken 769910 times.
✓ Branch 1 taken 674676 times.
|
1444586 | const ctype eps2 = eps_*std::max(b[5]-b[2], a[5]-a[2]); |
| 368 |
2/2✓ Branch 0 taken 1164379 times.
✓ Branch 1 taken 131588 times.
|
1295967 | return (b[0] - eps0 <= a[3] && a[0] <= b[3] + eps0 && |
| 369 |
4/4✓ Branch 0 taken 1102768 times.
✓ Branch 1 taken 61611 times.
✓ Branch 2 taken 1039763 times.
✓ Branch 3 taken 63005 times.
|
1164379 | b[1] - eps1 <= a[4] && a[1] <= b[4] + eps1 && |
| 370 |
6/6✓ Branch 0 taken 1295967 times.
✓ Branch 1 taken 148619 times.
✓ Branch 2 taken 986466 times.
✓ Branch 3 taken 53297 times.
✓ Branch 4 taken 58384 times.
✓ Branch 5 taken 928082 times.
|
2484349 | b[2] - eps2 <= a[5] && a[2] <= b[5] + eps2); |
| 371 | |||
| 372 | } | ||
| 373 | |||
| 374 | /*! | ||
| 375 | * \brief Check whether a bounding box is intersecting another bounding box (dimworld == 2) | ||
| 376 | * \param a Pointer to first bounding box coordinates | ||
| 377 | * \param b Pointer to second bounding box coordinates | ||
| 378 | */ | ||
| 379 | template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 2, int> = 0> | ||
| 380 | 2137844 | inline 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 |
2/2✓ Branch 0 taken 2111163 times.
✓ Branch 1 taken 26681 times.
|
2137844 | const ctype eps0 = eps_*std::max(b[2]-b[0], a[2]-a[0]); |
| 385 |
2/2✓ Branch 0 taken 2116201 times.
✓ Branch 1 taken 21643 times.
|
2137844 | const ctype eps1 = eps_*std::max(b[3]-b[1], a[3]-a[1]); |
| 386 |
2/2✓ Branch 0 taken 1691074 times.
✓ Branch 1 taken 203897 times.
|
1894971 | return (b[0] - eps0 <= a[2] && a[0] <= b[2] + eps0 && |
| 387 |
6/6✓ Branch 0 taken 1894971 times.
✓ Branch 1 taken 242873 times.
✓ Branch 2 taken 1471988 times.
✓ Branch 3 taken 219086 times.
✓ Branch 4 taken 202393 times.
✓ Branch 5 taken 1269595 times.
|
3828918 | b[1] - eps1 <= a[3] && a[1] <= b[3] + eps1); |
| 388 | } | ||
| 389 | |||
| 390 | /*! | ||
| 391 | * \brief Check whether a bounding box is intersecting another bounding box (dimworld == 1) | ||
| 392 | * \param a Pointer to first bounding box coordinates | ||
| 393 | * \param b Pointer to second bounding box coordinates | ||
| 394 | */ | ||
| 395 | template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 1, int> = 0> | ||
| 396 | inline 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 | |||
| 404 | /*! | ||
| 405 | * \brief Compute squared distance between point and bounding box | ||
| 406 | * \param point The point | ||
| 407 | * \param b Pointer to bounding box coordinates | ||
| 408 | */ | ||
| 409 | template<int dimworld, class ctype> | ||
| 410 | 11105630 | inline ctype squaredDistancePointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b) | |
| 411 | { | ||
| 412 | 11105630 | ctype squaredDistance = 0.0; | |
| 413 |
2/2✓ Branch 0 taken 308444 times.
✓ Branch 1 taken 154222 times.
|
44256088 | for (int d = 0; d < dimworld; ++d) |
| 414 | { | ||
| 415 |
4/4✓ Branch 0 taken 26776 times.
✓ Branch 1 taken 281668 times.
✓ Branch 2 taken 26776 times.
✓ Branch 3 taken 281668 times.
|
66300916 | if (point[d] < b[d]) |
| 416 | 22477392 | squaredDistance += (point[d] - b[d])*(point[d] - b[d]); | |
| 417 |
4/4✓ Branch 0 taken 112124 times.
✓ Branch 1 taken 196320 times.
✓ Branch 2 taken 112124 times.
✓ Branch 3 taken 196320 times.
|
66300916 | if (point[d] > b[d+dimworld]) |
| 418 | 56495163 | squaredDistance += (point[d] - b[d+dimworld])*(point[d] - b[d+dimworld]); | |
| 419 | } | ||
| 420 | 11105630 | return squaredDistance; | |
| 421 | } | ||
| 422 | |||
| 423 | } // end namespace Dumux | ||
| 424 | |||
| 425 | #endif | ||
| 426 |