| 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-FileCopyrightText: 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 <limits> | ||
| 28 | #include <type_traits> | ||
| 29 | #include <iostream> | ||
| 30 | |||
| 31 | #include <dune/common/promotiontraits.hh> | ||
| 32 | #include <dune/common/timer.hh> | ||
| 33 | #include <dune/common/fvector.hh> | ||
| 34 | |||
| 35 | namespace Dumux { | ||
| 36 | |||
| 37 | |||
| 38 | #ifndef DOXYGEN | ||
| 39 | namespace Detail { | ||
| 40 | |||
| 41 | template<typename ctype> | ||
| 42 | inline constexpr ctype minimumBaseEpsilon = 10.0*std::numeric_limits<ctype>::epsilon(); | ||
| 43 | |||
| 44 | } // namespace Detail | ||
| 45 | #endif // DOXYGEN | ||
| 46 | |||
| 47 | /*! | ||
| 48 | * \ingroup Geometry | ||
| 49 | * \brief An axis-aligned bounding box volume tree implementation | ||
| 50 | * | ||
| 51 | * The class constructs a hierarchical structure of bounding box volumes around | ||
| 52 | * grid entities. This class can be used to efficiently compute intersections | ||
| 53 | * between a grid and other geometrical object. It only implements the intersection | ||
| 54 | * of two of such bounding box trees, so that two independent grids can be intersected. | ||
| 55 | * \tparam GeometricEntitySet has the following requirements | ||
| 56 | * * export dimensionworld, ctype | ||
| 57 | * * a size() member function returning the number of entities | ||
| 58 | * * begin() and end() member function returning at least forward iterators to entities | ||
| 59 | * * an index() method returning a consecutive index given an entity | ||
| 60 | * * an entity() method returning an entity given the consecutive index | ||
| 61 | * * entities have the following requirements: | ||
| 62 | * * a member function geometry() returning a geometry with the member functions | ||
| 63 | * * corner() and corners() returning global coordinates and number of corners | ||
| 64 | */ | ||
| 65 | template <class GeometricEntitySet> | ||
| 66 | class BoundingBoxTree | ||
| 67 | { | ||
| 68 | enum { dimworld = GeometricEntitySet::dimensionworld }; | ||
| 69 | using ctype = typename GeometricEntitySet::ctype; | ||
| 70 | |||
| 71 | /*! | ||
| 72 | * \brief Bounding box node data structure | ||
| 73 | * leaf nodes are indicated by setting child0 to | ||
| 74 | * the node itself and child1 to the index of the entity in the bounding box. | ||
| 75 | */ | ||
| 76 | struct BoundingBoxNode | ||
| 77 | { | ||
| 78 | std::size_t child0; | ||
| 79 | std::size_t child1; | ||
| 80 | }; | ||
| 81 | |||
| 82 | public: | ||
| 83 | //! the type of entity set this tree was built with | ||
| 84 | using EntitySet = GeometricEntitySet; | ||
| 85 | |||
| 86 | //! Default Constructor | ||
| 87 | 26 | BoundingBoxTree() = default; | |
| 88 | |||
| 89 | //! Constructor with gridView | ||
| 90 | 4750 | explicit BoundingBoxTree(std::shared_ptr<const GeometricEntitySet> set) | |
| 91 |
1/8✓ Branch 2 taken 2664 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
9500 | { build(set); } |
| 92 | |||
| 93 | //! Build up bounding box tree for a grid with leafGridView | ||
| 94 | 4822 | void build(std::shared_ptr<const GeometricEntitySet> set) | |
| 95 | { | ||
| 96 | // set the pointer to the entity set | ||
| 97 | 4822 | entitySet_ = set; | |
| 98 | |||
| 99 | // clear all internal data | ||
| 100 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2698 times.
|
4822 | boundingBoxNodes_.clear(); |
| 101 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2698 times.
|
4830 | boundingBoxCoordinates_.clear(); |
| 102 | |||
| 103 | // start the timer | ||
| 104 | 4822 | Dune::Timer timer; | |
| 105 | |||
| 106 | // Create bounding boxes for all elements | ||
| 107 | 4822 | const auto numLeaves = set->size(); | |
| 108 | |||
| 109 | // reserve enough space for the nodes and the coordinates | ||
| 110 | 4822 | const auto numNodes = 2*numLeaves - 1; | |
| 111 | 4822 | boundingBoxNodes_.reserve(numNodes); | |
| 112 | 4822 | boundingBoxCoordinates_.reserve(numNodes*2*dimworld); | |
| 113 | |||
| 114 | // create a vector for leaf boxes (min and max for all dims) | ||
| 115 |
1/2✓ Branch 2 taken 815 times.
✗ Branch 3 not taken.
|
4822 | std::vector<ctype> leafBoxes(2*dimworld*numLeaves); |
| 116 | |||
| 117 |
14/17✓ Branch 1 taken 55630 times.
✓ Branch 2 taken 964 times.
✓ Branch 4 taken 80584 times.
✓ Branch 5 taken 334179 times.
✓ Branch 7 taken 787522 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 768374 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 768433 times.
✓ Branch 13 taken 242 times.
✓ Branch 16 taken 141 times.
✗ Branch 17 not taken.
✓ Branch 6 taken 2864779 times.
✓ Branch 9 taken 19695 times.
✓ Branch 3 taken 5929 times.
✓ Branch 0 taken 1744 times.
✓ Branch 15 taken 19122 times.
|
17656070 | for (const auto& geometricEntity : *set) |
| 118 |
3/5✓ Branch 1 taken 774303 times.
✓ Branch 2 taken 98442 times.
✓ Branch 4 taken 768374 times.
✗ Branch 5 not taken.
✗ Branch 3 not taken.
|
6103006 | computeEntityBoundingBox_(leafBoxes.data() + 2*dimworld*set->index(geometricEntity), geometricEntity); |
| 119 | |||
| 120 | // create the leaf partition, the set of available indices (to be sorted) | ||
| 121 |
1/4✓ Branch 1 taken 2702 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
4822 | std::vector<std::size_t> leafPartition(numLeaves); |
| 122 |
1/2✓ Branch 1 taken 2702 times.
✗ Branch 2 not taken.
|
9644 | std::iota(leafPartition.begin(), leafPartition.end(), 0); |
| 123 | |||
| 124 | // Recursively build the bounding box tree | ||
| 125 |
2/6✓ Branch 1 taken 2702 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2702 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
4822 | build_(leafBoxes, leafPartition.begin(), leafPartition.end()); |
| 126 | |||
| 127 | // We are done, log output | ||
| 128 |
1/2✓ Branch 1 taken 2702 times.
✗ Branch 2 not taken.
|
4822 | std::cout << "Computed bounding box tree with " << numBoundingBoxes() |
| 129 |
2/4✓ Branch 1 taken 2702 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2702 times.
✗ Branch 5 not taken.
|
4822 | << " nodes for " << numLeaves << " grid entities in " |
| 130 |
4/8✓ Branch 2 taken 2702 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2702 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 2702 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 1617 times.
✗ Branch 11 not taken.
|
4822 | << timer.stop() << " seconds." << std::endl; |
| 131 | 9643 | } | |
| 132 | |||
| 133 | //! the entity set this tree was built with | ||
| 134 | 4884490 | const EntitySet& entitySet() const | |
| 135 |
11/24✗ Branch 0 not taken.
✓ Branch 1 taken 2254 times.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 36 times.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✓ Branch 4 taken 83 times.
✓ Branch 5 taken 11340 times.
✓ Branch 3 taken 24 times.
✓ Branch 2 taken 11829 times.
✓ Branch 8 taken 4640 times.
✓ Branch 9 taken 1 times.
✓ Branch 12 taken 143776 times.
✗ Branch 6 not taken.
|
4640048 | { return *entitySet_; } |
| 136 | |||
| 137 | ///////////////////////////////////////////////////// | ||
| 138 | //! Interface to be used by other bounding box trees | ||
| 139 | ///////////////////////////////////////////////////// | ||
| 140 | |||
| 141 | //! Get an existing bounding box for a given node | ||
| 142 | 30195427 | const BoundingBoxNode& getBoundingBoxNode(std::size_t nodeIdx) const | |
| 143 |
24/24✓ Branch 0 taken 2783711 times.
✓ Branch 1 taken 2171654 times.
✓ Branch 2 taken 5439337 times.
✓ Branch 3 taken 5140836 times.
✓ Branch 4 taken 271887 times.
✓ Branch 5 taken 467418 times.
✓ Branch 6 taken 135151 times.
✓ Branch 7 taken 221909 times.
✓ Branch 8 taken 2202622 times.
✓ Branch 9 taken 4480198 times.
✓ Branch 10 taken 917311 times.
✓ Branch 11 taken 2385395 times.
✓ Branch 12 taken 2462 times.
✓ Branch 13 taken 1343 times.
✓ Branch 14 taken 1575 times.
✓ Branch 15 taken 1038 times.
✓ Branch 16 taken 86414 times.
✓ Branch 17 taken 60943 times.
✓ Branch 18 taken 65232 times.
✓ Branch 19 taken 48077 times.
✓ Branch 20 taken 1255532 times.
✓ Branch 21 taken 684532 times.
✓ Branch 22 taken 893421 times.
✓ Branch 23 taken 477429 times.
|
30195427 | { return boundingBoxNodes_[nodeIdx]; } |
| 144 | |||
| 145 | //! Get an existing bounding box for a given node | ||
| 146 | 27091171 | const ctype* getBoundingBoxCoordinates(std::size_t nodeIdx) const | |
| 147 |
24/24✓ Branch 0 taken 2787984 times.
✓ Branch 1 taken 2089395 times.
✓ Branch 2 taken 5441548 times.
✓ Branch 3 taken 5136945 times.
✓ Branch 4 taken 374715 times.
✓ Branch 5 taken 129184 times.
✓ Branch 6 taken 184215 times.
✓ Branch 7 taken 64062 times.
✓ Branch 8 taken 3479865 times.
✓ Branch 9 taken 1431601 times.
✓ Branch 10 taken 1652342 times.
✓ Branch 11 taken 741317 times.
✓ Branch 12 taken 2462 times.
✓ Branch 13 taken 1343 times.
✓ Branch 14 taken 1575 times.
✓ Branch 15 taken 1038 times.
✓ Branch 16 taken 86414 times.
✓ Branch 17 taken 60943 times.
✓ Branch 18 taken 65232 times.
✓ Branch 19 taken 48077 times.
✓ Branch 20 taken 1255532 times.
✓ Branch 21 taken 684532 times.
✓ Branch 22 taken 893421 times.
✓ Branch 23 taken 477429 times.
|
27091171 | { return boundingBoxCoordinates_.data() + 2*dimworld*nodeIdx; } |
| 148 | |||
| 149 | //! Get the number of bounding boxes currently in the tree | ||
| 150 | 4202239 | std::size_t numBoundingBoxes() const | |
| 151 |
33/47✓ Branch 1 taken 800 times.
✓ Branch 2 taken 120 times.
✓ Branch 4 taken 1098 times.
✓ Branch 5 taken 170 times.
✓ Branch 7 taken 59459 times.
✓ Branch 8 taken 50 times.
✓ Branch 11 taken 2 times.
✓ Branch 12 taken 120 times.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✓ Branch 18 taken 2 times.
✓ Branch 19 taken 140 times.
✓ Branch 21 taken 2 times.
✓ Branch 22 taken 9 times.
✓ Branch 25 taken 10 times.
✓ Branch 26 taken 120 times.
✓ Branch 28 taken 22 times.
✗ Branch 29 not taken.
✓ Branch 32 taken 6 times.
✓ Branch 33 taken 120 times.
✓ Branch 35 taken 6 times.
✗ Branch 36 not taken.
✓ Branch 39 taken 2 times.
✓ Branch 40 taken 120 times.
✓ Branch 42 taken 2 times.
✗ Branch 43 not taken.
✓ Branch 46 taken 2 times.
✗ Branch 47 not taken.
✓ Branch 49 taken 2 times.
✗ Branch 50 not taken.
✓ Branch 10 taken 282230 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 120 times.
✗ Branch 3 not taken.
✓ Branch 13 taken 14 times.
✓ Branch 16 taken 179 times.
✗ Branch 17 not taken.
✗ Branch 20 not taken.
✓ Branch 23 taken 120 times.
✗ Branch 24 not taken.
✗ Branch 27 not taken.
✓ Branch 30 taken 120 times.
✓ Branch 31 taken 16 times.
✓ Branch 34 taken 4 times.
✓ Branch 37 taken 128 times.
✗ Branch 38 not taken.
✗ Branch 41 not taken.
|
414893 | { return boundingBoxNodes_.size(); } |
| 152 | |||
| 153 | //! Check whether a bounding box node is a leaf node | ||
| 154 | //! Leaf nodes have itself as child0 | ||
| 155 | 21528582 | bool isLeaf(const BoundingBoxNode& node, std::size_t nodeIdx) const | |
| 156 |
24/24✓ Branch 0 taken 422263 times.
✓ Branch 1 taken 2455118 times.
✓ Branch 2 taken 297953 times.
✓ Branch 3 taken 5145849 times.
✓ Branch 4 taken 241967 times.
✓ Branch 5 taken 473955 times.
✓ Branch 6 taken 115965 times.
✓ Branch 7 taken 229070 times.
✓ Branch 8 taken 1915130 times.
✓ Branch 9 taken 4623970 times.
✓ Branch 10 taken 917311 times.
✓ Branch 11 taken 2385395 times.
✓ Branch 12 taken 700 times.
✓ Branch 13 taken 1762 times.
✓ Branch 14 taken 401 times.
✓ Branch 15 taken 1174 times.
✓ Branch 16 taken 14222 times.
✓ Branch 17 taken 72192 times.
✓ Branch 18 taken 10064 times.
✓ Branch 19 taken 55168 times.
✓ Branch 20 taken 290819 times.
✓ Branch 21 taken 964713 times.
✓ Branch 22 taken 213315 times.
✓ Branch 23 taken 680106 times.
|
21528582 | { return node.child0 == nodeIdx; } |
| 157 | |||
| 158 | private: | ||
| 159 | |||
| 160 | //! vector of bounding box nodes | ||
| 161 | std::vector<BoundingBoxNode> boundingBoxNodes_; | ||
| 162 | |||
| 163 | //! vector of bounding box coordinates | ||
| 164 | std::vector<ctype> boundingBoxCoordinates_; | ||
| 165 | |||
| 166 | //! a pointer to the entity set | ||
| 167 | std::shared_ptr<const EntitySet> entitySet_; | ||
| 168 | |||
| 169 | //! Compute the bounding box of a grid entity | ||
| 170 | template <class Entity> | ||
| 171 | 6101262 | void computeEntityBoundingBox_(ctype* b, const Entity& entity) const | |
| 172 | { | ||
| 173 | // get the bounding box coordinates | ||
| 174 | 6101262 | ctype* xMin = b; | |
| 175 | 6101262 | ctype* xMax = b + dimworld; | |
| 176 | |||
| 177 | // get mesh entity data | ||
| 178 | 6101262 | auto geometry = entity.geometry(); | |
| 179 | |||
| 180 | // Get coordinates of first vertex | ||
| 181 | 13373949 | auto corner = geometry.corner(0); | |
| 182 |
4/4✓ Branch 0 taken 8609346 times.
✓ Branch 1 taken 3983292 times.
✓ Branch 2 taken 139848 times.
✓ Branch 3 taken 528 times.
|
20969689 | for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx) |
| 183 | 14868427 | xMin[dimIdx] = xMax[dimIdx] = corner[dimIdx]; | |
| 184 | |||
| 185 | // Compute the min and max over the remaining vertices | ||
| 186 |
3/3✓ Branch 1 taken 3794364 times.
✓ Branch 2 taken 94698 times.
✓ Branch 0 taken 15323524 times.
|
32958162 | for (std::size_t cornerIdx = 1; cornerIdx < geometry.corners(); ++cornerIdx) |
| 187 | { | ||
| 188 | 26856900 | corner = geometry.corner(cornerIdx); | |
| 189 |
3/3✓ Branch 0 taken 38971314 times.
✓ Branch 1 taken 15700790 times.
✓ Branch 2 taken 245078 times.
|
97497515 | for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx) |
| 190 | { | ||
| 191 | using std::max; | ||
| 192 | using std::min; | ||
| 193 |
4/4✓ Branch 0 taken 1224912 times.
✓ Branch 1 taken 38239286 times.
✓ Branch 2 taken 8754432 times.
✓ Branch 3 taken 30709766 times.
|
72331197 | xMin[dimIdx] = min(xMin[dimIdx], corner[dimIdx]); |
| 194 |
2/2✓ Branch 0 taken 8754432 times.
✓ Branch 1 taken 30709766 times.
|
85299326 | xMax[dimIdx] = max(xMax[dimIdx], corner[dimIdx]); |
| 195 | } | ||
| 196 | } | ||
| 197 | 6099518 | } | |
| 198 | |||
| 199 | //! Build bounding box tree for all entities recursively | ||
| 200 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 7571990 times.
|
12201190 | std::size_t build_(const std::vector<ctype>& leafBoxes, |
| 201 | const std::vector<std::size_t>::iterator& begin, | ||
| 202 | const std::vector<std::size_t>::iterator& end) | ||
| 203 | { | ||
| 204 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 7571990 times.
|
12201190 | assert(begin < end); |
| 205 | |||
| 206 | // If we reached the end of the recursion, i.e. only a leaf box is left | ||
| 207 |
2/2✓ Branch 0 taken 3787346 times.
✓ Branch 1 taken 3784644 times.
|
12201190 | if (end - begin == 1) |
| 208 | { | ||
| 209 | // Get the bounding box coordinates for the leaf | ||
| 210 | 6103006 | const std::size_t leafNodeIdx = *begin; | |
| 211 | 6103006 | const auto beginCoords = leafBoxes.begin() + 2*dimworld*leafNodeIdx; | |
| 212 | 6103006 | const auto endCoords = beginCoords + 2*dimworld; | |
| 213 | |||
| 214 | // Store the data in the bounding box | ||
| 215 | // leaf nodes are indicated by setting child0 to | ||
| 216 | // the node itself and child1 to the index of the entity in the bounding box. | ||
| 217 | 6103006 | return addBoundingBox_(BoundingBoxNode{numBoundingBoxes(), leafNodeIdx}, beginCoords, endCoords); | |
| 218 | } | ||
| 219 | |||
| 220 | // Compute the bounding box of all bounding boxes in the range [begin, end] | ||
| 221 | 6098184 | const auto bCoords = computeBBoxOfBBoxes_(leafBoxes, begin, end); | |
| 222 | |||
| 223 | // sort bounding boxes along the longest axis | ||
| 224 | 6098184 | const auto axis = computeLongestAxis_(bCoords); | |
| 225 | |||
| 226 | // nth_element sorts the range to make sure that middle points to the coordinate median in axis direction | ||
| 227 | // this is the most expensive part of the algorithm | ||
| 228 | 6098184 | auto middle = begin + (end - begin)/2; | |
| 229 | 103430506 | std::nth_element(begin, middle, end, [&leafBoxes, axis](std::size_t i, std::size_t j) | |
| 230 | { | ||
| 231 | 151565169 | const ctype* bi = leafBoxes.data() + 2*dimworld*i; | |
| 232 | 151565169 | const ctype* bj = leafBoxes.data() + 2*dimworld*j; | |
| 233 | 151565169 | return bi[axis] + bi[axis + dimworld] < bj[axis] + bj[axis + dimworld]; | |
| 234 | }); | ||
| 235 | |||
| 236 | // split the bounding boxes into two at the middle iterator and call build recursively, each | ||
| 237 | // call resulting in a new node of this bounding box, i.e. the root will be added at the end of the process. | ||
| 238 | 6098184 | return addBoundingBox_(BoundingBoxNode{build_(leafBoxes, begin, middle), build_(leafBoxes, middle, end)}, | |
| 239 | 12196368 | bCoords.begin(), bCoords.end()); | |
| 240 | } | ||
| 241 | |||
| 242 | //! Add a new bounding box to the tree | ||
| 243 | template <class Iterator> | ||
| 244 | 15143980 | std::size_t addBoundingBox_(BoundingBoxNode&& node, | |
| 245 | const Iterator& coordBegin, | ||
| 246 | const Iterator& coordEnd) | ||
| 247 | { | ||
| 248 | // Add the bounding box | ||
| 249 | 15143980 | boundingBoxNodes_.emplace_back(node); | |
| 250 | |||
| 251 | // Add the bounding box's coordinates | ||
| 252 | 15143980 | boundingBoxCoordinates_.insert(boundingBoxCoordinates_.end(), coordBegin, coordEnd); | |
| 253 | |||
| 254 | // return the index of the new node | ||
| 255 | 15143980 | return boundingBoxNodes_.size() - 1; | |
| 256 | } | ||
| 257 | |||
| 258 | //! Compute the bounding box of a vector of bounding boxes | ||
| 259 | std::array<ctype, 2*dimworld> | ||
| 260 | 6098184 | computeBBoxOfBBoxes_(const std::vector<ctype>& leafBoxes, | |
| 261 | const std::vector<std::size_t>::iterator& begin, | ||
| 262 | const std::vector<std::size_t>::iterator& end) | ||
| 263 | { | ||
| 264 | std::array<ctype, 2*dimworld> bBoxCoords; | ||
| 265 | |||
| 266 | // copy the iterator and get coordinates of first box | ||
| 267 | 6098184 | auto it = begin; | |
| 268 | 6098184 | const auto* bFirst = leafBoxes.data() + 2*dimworld*(*it); | |
| 269 | |||
| 270 |
2/2✓ Branch 0 taken 17877982 times.
✓ Branch 1 taken 3784644 times.
|
35817670 | for (int coordIdx = 0; coordIdx < 2*dimworld; ++coordIdx) |
| 271 | 29719486 | bBoxCoords[coordIdx] = bFirst[coordIdx]; | |
| 272 | |||
| 273 | // Compute min and max with the remaining boxes | ||
| 274 |
2/2✓ Branch 0 taken 55408969 times.
✓ Branch 1 taken 3784644 times.
|
96848063 | for (; it != end; ++it) |
| 275 | { | ||
| 276 | 90749879 | const auto* b = leafBoxes.data() + 2*dimworld*(*it); | |
| 277 |
2/2✓ Branch 0 taken 132557576 times.
✓ Branch 1 taken 55408969 times.
|
314908250 | for (int coordIdx = 0; coordIdx < dimworld; ++coordIdx) |
| 278 |
2/2✓ Branch 0 taken 4408799 times.
✓ Branch 1 taken 128148777 times.
|
224021459 | if (b[coordIdx] < bBoxCoords[coordIdx]) |
| 279 | 6940676 | bBoxCoords[coordIdx] = b[coordIdx]; | |
| 280 |
2/2✓ Branch 0 taken 132557576 times.
✓ Branch 1 taken 55408969 times.
|
314771338 | for (int coordIdx = dimworld; coordIdx < 2*dimworld; ++coordIdx) |
| 281 |
2/2✓ Branch 0 taken 7069096 times.
✓ Branch 1 taken 125488480 times.
|
224021459 | if (b[coordIdx] > bBoxCoords[coordIdx]) |
| 282 | 11249945 | bBoxCoords[coordIdx] = b[coordIdx]; | |
| 283 | } | ||
| 284 | |||
| 285 | 6098184 | return bBoxCoords; | |
| 286 | } | ||
| 287 | |||
| 288 | //! Compute the bounding box of a vector of bounding boxes | ||
| 289 | 5638730 | std::size_t computeLongestAxis_(const std::array<ctype, 2*dimworld>& bCoords) | |
| 290 | { | ||
| 291 | std::array<ctype, dimworld> axisLength; | ||
| 292 |
2/2✓ Branch 0 taken 8002036 times.
✓ Branch 1 taken 3326006 times.
|
19559904 | for (int coordIdx = 0; coordIdx < dimworld; ++coordIdx) |
| 293 | 13921174 | axisLength[coordIdx] = bCoords[dimworld + coordIdx] - bCoords[coordIdx]; | |
| 294 | |||
| 295 | 5638730 | return std::distance(axisLength.begin(), std::max_element(axisLength.begin(), axisLength.end())); | |
| 296 | } | ||
| 297 | }; | ||
| 298 | |||
| 299 | /*! | ||
| 300 | * \brief Check whether a point is intersectin a bounding box (dimworld == 3) | ||
| 301 | * \param point The point | ||
| 302 | * \param b Pointer to bounding box coordinates | ||
| 303 | */ | ||
| 304 | template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 3, int> = 0> | ||
| 305 | 25049136 | inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b) | |
| 306 | { | ||
| 307 | static constexpr ctype eps_ = 1.0e-7; | ||
| 308 | |||
| 309 | using std::max; | ||
| 310 | 25049136 | const auto dx = b[3] - b[0]; | |
| 311 | 25049136 | const auto dy = b[4] - b[1]; | |
| 312 | 25049136 | const auto dz = b[5] - b[2]; | |
| 313 | 25049136 | const ctype eps = max({dx, dy, dz})*eps_; | |
| 314 |
6/6✓ Branch 0 taken 23899873 times.
✓ Branch 1 taken 1149263 times.
✓ Branch 2 taken 1158152 times.
✓ Branch 3 taken 22741721 times.
✓ Branch 4 taken 953988 times.
✓ Branch 5 taken 21787733 times.
|
25049136 | return (b[0] - eps <= point[0] && point[0] <= b[3] + eps && |
| 315 |
6/6✓ Branch 0 taken 953988 times.
✓ Branch 1 taken 21787733 times.
✓ Branch 2 taken 982602 times.
✓ Branch 3 taken 20805131 times.
✓ Branch 4 taken 806265 times.
✓ Branch 5 taken 19998866 times.
|
22741721 | b[1] - eps <= point[1] && point[1] <= b[4] + eps && |
| 316 |
6/6✓ Branch 0 taken 23899873 times.
✓ Branch 1 taken 1149263 times.
✓ Branch 2 taken 806265 times.
✓ Branch 3 taken 19998866 times.
✓ Branch 4 taken 874072 times.
✓ Branch 5 taken 19124794 times.
|
45854267 | b[2] - eps <= point[2] && point[2] <= b[5] + eps); |
| 317 | } | ||
| 318 | |||
| 319 | /*! | ||
| 320 | * \brief Check whether a point is intersectin a bounding box (dimworld == 2) | ||
| 321 | * \param point The point | ||
| 322 | * \param b Pointer to bounding box coordinates | ||
| 323 | */ | ||
| 324 | template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 2, int> = 0> | ||
| 325 | 760440 | inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b) | |
| 326 | { | ||
| 327 | static constexpr ctype eps_ = 1.0e-7; | ||
| 328 | |||
| 329 | using std::max; | ||
| 330 | 760440 | const auto dx = b[2] - b[0]; | |
| 331 |
2/2✓ Branch 0 taken 371731 times.
✓ Branch 1 taken 388709 times.
|
760440 | const auto dy = b[3] - b[1]; |
| 332 | 760440 | const ctype eps = max(dx, dy)*eps_; | |
| 333 |
6/6✓ Branch 0 taken 661804 times.
✓ Branch 1 taken 98636 times.
✓ Branch 2 taken 99295 times.
✓ Branch 3 taken 562509 times.
✓ Branch 4 taken 98760 times.
✓ Branch 5 taken 463749 times.
|
760440 | return (b[0] - eps <= point[0] && point[0] <= b[2] + eps && |
| 334 |
6/6✓ Branch 0 taken 661804 times.
✓ Branch 1 taken 98636 times.
✓ Branch 2 taken 98760 times.
✓ Branch 3 taken 463749 times.
✓ Branch 4 taken 105961 times.
✓ Branch 5 taken 357788 times.
|
1322949 | b[1] - eps <= point[1] && point[1] <= b[3] + eps); |
| 335 | } | ||
| 336 | |||
| 337 | /*! | ||
| 338 | * \brief Check whether a point is intersectin a bounding box (dimworld == 1) | ||
| 339 | * \param point The point | ||
| 340 | * \param b Pointer to bounding box coordinates | ||
| 341 | */ | ||
| 342 | template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 1, int> = 0> | ||
| 343 | 214 | inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b) | |
| 344 | { | ||
| 345 | static constexpr ctype eps_ = 1.0e-7; | ||
| 346 | 214 | const ctype eps0 = eps_*(b[1] - b[0]); | |
| 347 |
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; |
| 348 | } | ||
| 349 | |||
| 350 | /*! | ||
| 351 | * \ingroup Geometry | ||
| 352 | * \brief Determine if a point intersects an axis-aligned bounding box | ||
| 353 | * The bounding box is given by the lower left corner (min) and the upper right corner (max) | ||
| 354 | */ | ||
| 355 | template<class ctype, int dimworld> | ||
| 356 | 12767952 | inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, | |
| 357 | const Dune::FieldVector<ctype, dimworld>& min, | ||
| 358 | const Dune::FieldVector<ctype, dimworld>& max) | ||
| 359 | { | ||
| 360 | std::array<ctype, 2*dimworld> bBox; | ||
| 361 | 12767952 | std::copy(min.begin(), min.end(), bBox.begin()); | |
| 362 | 12767952 | std::copy(max.begin(), max.end(), bBox.begin()+dimworld); | |
| 363 | 12767952 | return intersectsPointBoundingBox(point, bBox.data()); | |
| 364 | } | ||
| 365 | |||
| 366 | /*! | ||
| 367 | * \brief Check whether a bounding box is intersecting another bounding box (dimworld == 3) | ||
| 368 | * \param a Pointer to first bounding box coordinates | ||
| 369 | * \param b Pointer to second bounding box coordinates | ||
| 370 | */ | ||
| 371 | template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 3, int> = 0> | ||
| 372 | 1449824 | inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b) | |
| 373 | { | ||
| 374 | using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType; | ||
| 375 | static constexpr ctype eps_ = 1.0e-7; | ||
| 376 |
2/2✓ Branch 0 taken 650791 times.
✓ Branch 1 taken 799033 times.
|
1449824 | const ctype scale0 = std::max(b[3]-b[0], a[3]-a[0]); |
| 377 |
2/2✓ Branch 0 taken 648220 times.
✓ Branch 1 taken 801604 times.
|
1449824 | const ctype scale1 = std::max(b[4]-b[1], a[4]-a[1]); |
| 378 |
4/4✓ Branch 0 taken 773678 times.
✓ Branch 1 taken 676146 times.
✓ Branch 2 taken 487138 times.
✓ Branch 3 taken 962686 times.
|
2223502 | const ctype scale2 = std::max(b[5]-b[2], a[5]-a[2]); |
| 379 |
2/2✓ Branch 0 taken 679820 times.
✓ Branch 1 taken 770004 times.
|
1449824 | const ctype maxScale = std::max(scale0, std::max(scale1, scale2)); |
| 380 | 1449824 | const ctype minEps = maxScale*Detail::minimumBaseEpsilon<ctype>; | |
| 381 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1449824 times.
|
1449824 | const ctype eps0 = std::max(eps_*scale0, minEps); |
| 382 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1449824 times.
|
1449824 | const ctype eps1 = std::max(eps_*scale1, minEps); |
| 383 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1449824 times.
|
1449824 | const ctype eps2 = std::max(eps_*scale2, minEps); |
| 384 |
2/2✓ Branch 0 taken 1169361 times.
✓ Branch 1 taken 131704 times.
|
1301065 | return (b[0] - eps0 <= a[3] && a[0] <= b[3] + eps0 && |
| 385 |
4/4✓ Branch 0 taken 1107750 times.
✓ Branch 1 taken 61611 times.
✓ Branch 2 taken 1044745 times.
✓ Branch 3 taken 63005 times.
|
1169361 | b[1] - eps1 <= a[4] && a[1] <= b[4] + eps1 && |
| 386 |
6/6✓ Branch 0 taken 1301065 times.
✓ Branch 1 taken 148759 times.
✓ Branch 2 taken 991448 times.
✓ Branch 3 taken 53297 times.
✓ Branch 4 taken 58384 times.
✓ Branch 5 taken 933064 times.
|
2494569 | b[2] - eps2 <= a[5] && a[2] <= b[5] + eps2); |
| 387 | |||
| 388 | } | ||
| 389 | |||
| 390 | /*! | ||
| 391 | * \brief Check whether a bounding box is intersecting another bounding box (dimworld == 2) | ||
| 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 == 2, int> = 0> | ||
| 396 | 2138999 | 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 |
2/2✓ Branch 0 taken 2111540 times.
✓ Branch 1 taken 27459 times.
|
2138999 | const ctype scale0 = std::max(b[2]-b[0], a[2]-a[0]); |
| 401 |
4/4✓ Branch 0 taken 2116509 times.
✓ Branch 1 taken 22490 times.
✓ Branch 2 taken 950618 times.
✓ Branch 3 taken 1188381 times.
|
4255508 | const ctype scale1 = std::max(b[3]-b[1], a[3]-a[1]); |
| 402 | 2138999 | const ctype maxScale = std::max(scale0, scale1); | |
| 403 | 2138999 | const ctype minEps = maxScale*Detail::minimumBaseEpsilon<ctype>; | |
| 404 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2138999 times.
|
2138999 | const ctype eps0 = std::max(eps_*scale0, minEps); |
| 405 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2138999 times.
|
2138999 | const ctype eps1 = std::max(eps_*scale1, minEps); |
| 406 |
2/2✓ Branch 0 taken 1692032 times.
✓ Branch 1 taken 203997 times.
|
1896029 | return (b[0] - eps0 <= a[2] && a[0] <= b[2] + eps0 && |
| 407 |
6/6✓ Branch 0 taken 1896029 times.
✓ Branch 1 taken 242970 times.
✓ Branch 2 taken 1472810 times.
✓ Branch 3 taken 219222 times.
✓ Branch 4 taken 202478 times.
✓ Branch 5 taken 1270332 times.
|
3831031 | b[1] - eps1 <= a[3] && a[1] <= b[3] + eps1); |
| 408 | } | ||
| 409 | |||
| 410 | /*! | ||
| 411 | * \brief Check whether a bounding box is intersecting another bounding box (dimworld == 1) | ||
| 412 | * \param a Pointer to first bounding box coordinates | ||
| 413 | * \param b Pointer to second bounding box coordinates | ||
| 414 | */ | ||
| 415 | template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 1, int> = 0> | ||
| 416 | inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b) | ||
| 417 | { | ||
| 418 | using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType; | ||
| 419 | static constexpr ctype eps_ = 1.0e-7; | ||
| 420 | const ctype scale0 = std::max(b[1]-b[0], a[1]-a[0]); | ||
| 421 | const ctype eps0 = std::max(eps_*scale0, Detail::minimumBaseEpsilon<ctype>*scale0); | ||
| 422 | return b[0] - eps0 <= a[1] && a[0] <= b[1] + eps0; | ||
| 423 | } | ||
| 424 | |||
| 425 | /*! | ||
| 426 | * \brief Compute squared distance between point and bounding box | ||
| 427 | * \param point The point | ||
| 428 | * \param b Pointer to bounding box coordinates | ||
| 429 | */ | ||
| 430 | template<int dimworld, class ctype> | ||
| 431 | 11105602 | inline ctype squaredDistancePointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b) | |
| 432 | { | ||
| 433 | 11105602 | ctype squaredDistance = 0.0; | |
| 434 |
2/2✓ Branch 0 taken 308444 times.
✓ Branch 1 taken 154222 times.
|
44255994 | for (int d = 0; d < dimworld; ++d) |
| 435 | { | ||
| 436 |
2/2✓ Branch 0 taken 26776 times.
✓ Branch 1 taken 281668 times.
|
33150392 | if (point[d] < b[d]) |
| 437 | 7492563 | squaredDistance += (point[d] - b[d])*(point[d] - b[d]); | |
| 438 |
2/2✓ Branch 0 taken 112124 times.
✓ Branch 1 taken 196320 times.
|
33150392 | if (point[d] > b[d+dimworld]) |
| 439 | 18831684 | squaredDistance += (point[d] - b[d+dimworld])*(point[d] - b[d+dimworld]); | |
| 440 | } | ||
| 441 | 11105602 | return squaredDistance; | |
| 442 | } | ||
| 443 | |||
| 444 | } // end namespace Dumux | ||
| 445 | |||
| 446 | #endif | ||
| 447 |