GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/geometry/intersectionentityset.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 58 61 95.1%
Functions: 30 35 85.7%
Branches: 128 259 49.4%
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 A class representing the intersection entities of two geometric entity sets
11 */
12
13 #ifndef DUMUX_GEOMETRY_INTERSECTION_ENTITY_SET_HH
14 #define DUMUX_GEOMETRY_INTERSECTION_ENTITY_SET_HH
15
16 #include <algorithm>
17 #include <cassert>
18 #include <iostream>
19 #include <memory>
20 #include <type_traits>
21 #include <utility>
22 #include <vector>
23
24 #include <dune/common/indices.hh>
25 #include <dune/common/timer.hh>
26 #include <dune/common/iteratorrange.hh>
27 #include <dune/common/promotiontraits.hh>
28 #include <dune/common/reservedvector.hh>
29 #include <dune/geometry/affinegeometry.hh>
30 #include <dune/geometry/type.hh>
31
32 #include <dumux/geometry/boundingboxtree.hh>
33 #include <dumux/geometry/intersectingentities.hh>
34
35 namespace Dumux {
36
37 /*!
38 * \ingroup Geometry
39 * \brief A class representing the intersection entities two geometric entity sets
40 */
41 template<class DomainEntitySet, class TargetEntitySet>
42 class IntersectionEntitySet
43 {
44 using DomainTree = BoundingBoxTree<DomainEntitySet>;
45 using TargetTree = BoundingBoxTree<TargetEntitySet>;
46
47 using ctype = typename Dune::PromotionTraits<typename DomainEntitySet::ctype, typename TargetEntitySet::ctype>::PromotedType;
48
49 static constexpr int dimWorld = DomainEntitySet::dimensionworld;
50 static_assert(dimWorld == int(TargetEntitySet::dimensionworld), "Entity sets must have the same world dimension");
51
52 using GlobalPosition = Dune::FieldVector<ctype, dimWorld>;
53
54 static constexpr int dimDomain = DomainEntitySet::Entity::Geometry::mydimension;
55 static constexpr int dimTarget = TargetEntitySet::Entity::Geometry::mydimension;
56 static constexpr bool isMixedDimensional = dimDomain != dimTarget;
57
58 /*!
59 * \brief A class representing an intersection entity
60 */
61
3/19
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 321312 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 310 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 635568 times.
✗ Branch 18 not taken.
 1594109 class IntersectionEntity
62 {
63 static constexpr int dimIs = std::min(dimDomain, dimTarget);
64 using Geometry = Dune::AffineGeometry<ctype, dimIs, dimWorld>; // geometries are always simplices
65
66 // we can only have multiple neighbors in the mixeddimensional case and then only for the side with the largest dimension
67 using IndexStorage = std::pair<std::conditional_t<dimDomain <= dimTarget, Dune::ReservedVector<std::size_t, 1>, std::vector<std::size_t>>,
68 std::conditional_t<dimTarget <= dimDomain, Dune::ReservedVector<std::size_t, 1>, std::vector<std::size_t>>>;
69
70 static constexpr auto domainIdx = Dune::index_constant<0>{};
71 static constexpr auto targetIdx = Dune::index_constant<1>{};
72
73 public:
74 350430 IntersectionEntity(const DomainTree& domainTree, const TargetTree& targetTree)
75 : domainTree_(domainTree)
76 1051290 , targetTree_(targetTree)
77 {}
78
79 //! set the intersection geometry corners
80 355211 void setCorners(const std::vector<GlobalPosition>& corners)
81 {
82 355211 corners_ = corners;
83
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 350430 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 350430 times.
 710422 assert(corners.size() == dimIs + 1); // Only simplex intersections are allowed
84 355211 }
85
86 //! add a pair of neighbor elements
87 356208 void addNeighbors(std::size_t domain, std::size_t target)
88 {
89 1422838 if (numDomainNeighbors() == 0 && numTargetNeighbors() == 0)
90 {
91 710422 std::get<domainIdx>(neighbors_).push_back(domain);
92 710422 std::get<targetIdx>(neighbors_).push_back(target);
93 }
94 else if (dimDomain > dimTarget)
95 1994 std::get<domainIdx>(neighbors_).push_back(domain);
96
97 else if (dimTarget > dimDomain)
98 ✗ std::get<targetIdx>(neighbors_).push_back(target);
99
100 else
101 ✗ DUNE_THROW(Dune::InvalidStateException, "Cannot add more than one neighbor per side for equidimensional intersection!");
102 356208 }
103
104 //! get the intersection geometry
105 Geometry geometry() const
106
3/6
✓ Branch 1 taken 316020 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2873 times.
✓ Branch 5 taken 7056 times.
✗ Branch 6 not taken.
 732352 { return Geometry(Dune::GeometryTypes::simplex(dimIs), corners_); }
107
108 //! get the number of domain neighbors of this intersection
109 std::size_t numDomainNeighbors() const
110
19/24
✓ Branch 0 taken 348144 times.
✓ Branch 1 taken 997 times.
✓ Branch 2 taken 348144 times.
✓ Branch 3 taken 997 times.
✓ Branch 4 taken 349012 times.
✓ Branch 5 taken 345669 times.
✓ Branch 6 taken 349012 times.
✓ Branch 7 taken 345669 times.
✓ Branch 8 taken 9241 times.
✓ Branch 9 taken 34016 times.
✓ Branch 10 taken 9241 times.
✓ Branch 11 taken 9559 times.
✓ Branch 12 taken 24457 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 24457 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✓ Branch 18 taken 24457 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 1789 times.
✓ Branch 21 taken 1444 times.
✓ Branch 22 taken 1789 times.
✓ Branch 23 taken 1444 times.
 3511714 { return std::get<domainIdx>(neighbors_).size(); }
111
112 //! get the number of target neighbors of this intersection
113 std::size_t numTargetNeighbors() const
114
8/16
✓ Branch 0 taken 348144 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 348144 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2286 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2286 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1968 times.
✓ Branch 9 taken 2286 times.
✓ Branch 10 taken 1968 times.
✓ Branch 11 taken 2286 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
 709368 { return std::get<targetIdx>(neighbors_).size(); }
115
116 //! get the nth domain neighbor entity
117 typename DomainEntitySet::Entity domainEntity(unsigned int n = 0) const
118
16/32
✓ Branch 1 taken 339075 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 339075 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 339075 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 338757 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 635656 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 635656 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 635656 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 635656 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 7056 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 7056 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 7056 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 7056 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 28224 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 28224 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 28224 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 28224 times.
✗ Branch 47 not taken.
 4087498 { return domainTree_.entitySet().entity(std::get<domainIdx>(neighbors_)[n]); }
119
120 //! get the nth target neighbor entity
121 318 typename TargetEntitySet::Entity targetEntity(unsigned int n = 0) const
122
20/36
✓ Branch 1 taken 314256 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 314256 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 314256 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1968 times.
✓ Branch 10 taken 314256 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1968 times.
✓ Branch 13 taken 314256 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 1968 times.
✓ Branch 16 taken 314256 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 1968 times.
✓ Branch 19 taken 314256 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 314256 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 7056 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 7056 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 7056 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 7056 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 7056 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 7056 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 7056 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 7056 times.
✗ Branch 47 not taken.
 2757912 { return targetTree_.entitySet().entity(std::get<targetIdx>(neighbors_)[n]); }
123
124 private:
125 IndexStorage neighbors_;
126 std::vector<GlobalPosition> corners_;
127
128 const DomainTree& domainTree_;
129 const TargetTree& targetTree_;
130 };
131
132 using Intersections = std::vector<IntersectionEntity>;
133
134 public:
135 //! make intersection entity type available
136 using Entity = IntersectionEntity;
137 //! make entity iterator type available
138 using EntityIterator = typename Intersections::const_iterator;
139
140 /*!
141 * \brief Default constructor
142 */
143
4/8
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 15 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 15 times.
✗ Branch 11 not taken.
 312 IntersectionEntitySet() = default;
144
145 /*!
146 * \brief Build intersections
147 */
148 74 void build(std::shared_ptr<const DomainEntitySet> domainSet, std::shared_ptr<const TargetEntitySet> targetSet)
149 {
150
2/4
✗ Branch 1 not taken.
✓ Branch 2 taken 37 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 37 times.
 74 domainTree_ = std::make_shared<DomainTree>(domainSet);
151
2/4
✗ Branch 1 not taken.
✓ Branch 2 taken 37 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 37 times.
 74 targetTree_ = std::make_shared<TargetTree>(targetSet);
152 222 build(*domainTree_, *targetTree_);
153 74 }
154
155 /*!
156 * \brief Build intersections
157 */
158 void build(std::shared_ptr<const DomainTree> domainTree, std::shared_ptr<const TargetTree> targetTree)
159 {
160 // make sure the tree don't get out of scope
161 domainTree_ = domainTree;
162 targetTree_ = targetTree;
163 build(*domainTree_, *targetTree_);
164 }
165
166 /*!
167 * \brief Build intersections
168 * \note If you call this, make sure the bounding box tree stays alive for the life-time of this object
169 */
170 117 void build(const DomainTree& domainTree, const TargetTree& targetTree)
171 {
172 117 Dune::Timer timer;
173
174 // compute raw intersections
175 234 const auto rawIntersections = intersectingEntities(domainTree, targetTree);
176
177 // create a map to check whether intersection geometries were already inserted
178 // Note that this can only occur if the grids have different dimensionality.
179 // If this is the case, we keep track of the intersections using the indices of the lower-
180 // dimensional entities which is identical for all geometrically identical intersections.
181
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 234 std::vector<std::vector<std::vector<GlobalPosition>>> intersectionMap;
182
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 117 std::vector<std::vector<std::size_t>> intersectionIndex;
183 if constexpr (isMixedDimensional)
184 {
185
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 64 const auto numLowDimEntities = dimTarget < dimDomain ? targetTree.entitySet().size()
186 4 : domainTree.entitySet().size();
187
1/2
✓ Branch 1 taken 29 times.
✗ Branch 2 not taken.
 32 intersectionMap.resize(numLowDimEntities);
188
1/2
✓ Branch 1 taken 29 times.
✗ Branch 2 not taken.
 32 intersectionIndex.resize(numLowDimEntities);
189 }
190
191 // reserve memory for storing the intersections. In case of grids of
192 // different dimensionality this might be an overestimate. We get rid
193 // of the overhead memory at the end of this function though.
194
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 29 times.
 117 intersections_.clear();
195
2/4
✓ Branch 1 taken 78 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 78 times.
✗ Branch 5 not taken.
 234 intersections_.reserve(rawIntersections.size());
196
197
4/4
✓ Branch 0 taken 351427 times.
✓ Branch 1 taken 78 times.
✓ Branch 2 taken 351427 times.
✓ Branch 3 taken 78 times.
 356559 for (const auto& rawIntersection : rawIntersections)
198 {
199 356208 bool add = true;
200
201 // Check if intersection was already inserted.
202 // In this case we only add new neighbor information as the geometry is identical.
203 if constexpr (isMixedDimensional)
204 {
205 82664 const auto lowDimNeighborIdx = getLowDimNeighborIdx_(rawIntersection);
206
6/6
✓ Branch 0 taken 490324 times.
✓ Branch 1 taken 35864 times.
✓ Branch 2 taken 490324 times.
✓ Branch 3 taken 35864 times.
✓ Branch 4 taken 490324 times.
✓ Branch 5 taken 35864 times.
 614035 for (int i = 0; i < intersectionMap[lowDimNeighborIdx].size(); ++i)
207 {
208
2/2
✓ Branch 2 taken 997 times.
✓ Branch 3 taken 489327 times.
 982072 if (rawIntersection.cornersMatch(intersectionMap[lowDimNeighborIdx][i]))
209 {
210 997 add = false;
211 // only add the pair of neighbors using the insertionIndex
212
2/4
✓ Branch 1 taken 997 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 997 times.
✗ Branch 5 not taken.
 1994 auto idx = intersectionIndex[lowDimNeighborIdx][i];
213
2/4
✓ Branch 1 taken 997 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 997 times.
✗ Branch 5 not taken.
 1994 intersections_[idx].addNeighbors(rawIntersection.first(), rawIntersection.second());
214 break;
215 }
216 }
217 }
218
219 if(add)
220 {
221 // maybe add to the map
222 if constexpr (isMixedDimensional)
223 {
224
2/4
✓ Branch 1 taken 35864 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 35864 times.
✗ Branch 5 not taken.
 80670 intersectionMap[getLowDimNeighborIdx_(rawIntersection)].push_back(rawIntersection.corners());
225
2/4
✓ Branch 1 taken 35864 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 35864 times.
✗ Branch 5 not taken.
 80670 intersectionIndex[getLowDimNeighborIdx_(rawIntersection)].push_back(intersections_.size());
226 }
227
228 // add new intersection and add the neighbors
229
1/2
✓ Branch 1 taken 350430 times.
✗ Branch 2 not taken.
 355211 intersections_.emplace_back(domainTree, targetTree);
230
3/6
✓ Branch 1 taken 350430 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 350430 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 314566 times.
✗ Branch 8 not taken.
 1025298 intersections_.back().setCorners(rawIntersection.corners());
231
2/4
✓ Branch 1 taken 350430 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 350430 times.
✗ Branch 5 not taken.
 710422 intersections_.back().addNeighbors(rawIntersection.first(), rawIntersection.second());
232 }
233 }
234
235
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 60 times.
 117 intersections_.shrink_to_fit();
236
3/6
✓ Branch 2 taken 78 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 78 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 78 times.
✗ Branch 11 not taken.
 351 std::cout << "Computed " << size() << " intersection entities in " << timer.elapsed() << std::endl;
237 117 }
238
239 //! return begin iterator to intersection container
240 typename Intersections::const_iterator ibegin() const
241 82 { return intersections_.begin(); }
242
243 //! return end iterator to intersection container
244 typename Intersections::const_iterator iend() const
245 82 { return intersections_.end(); }
246
247 //! the number of intersections
248 std::size_t size() const
249
5/16
✗ Branch 0 not taken.
✓ Branch 1 taken 76 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 36 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
 124 { return intersections_.size(); }
250
251 /*!
252 * \brief Range generator to iterate with range-based for loops over all intersections
253 * as follows: for (const auto& is : intersections(glue)) { ... }
254 * \note free function
255 */
256 friend Dune::IteratorRange<EntityIterator> intersections(const IntersectionEntitySet& set)
257 246 { return {set.ibegin(), set.iend()}; }
258
259 private:
260 template<class RawIntersection,
261 bool enable = isMixedDimensional, std::enable_if_t<enable, int> = 0>
262 ✗ auto getLowDimNeighborIdx_(const RawIntersection& is)
263 {
264 if constexpr (dimTarget < dimDomain)
265
2/4
✓ Branch 1 taken 33896 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 33896 times.
✗ Branch 5 not taken.
 68789 return is.second();
266 else
267
2/4
✓ Branch 1 taken 1968 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1968 times.
✗ Branch 5 not taken.
 3936 return is.first();
268 }
269
270 Intersections intersections_;
271
272 std::shared_ptr<const DomainTree> domainTree_;
273 std::shared_ptr<const TargetTree> targetTree_;
274 };
275
276 } // end namespace Dumux
277
278 #endif
279