GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/multidomain/facet/enrichmenthelper.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 94 94 100.0%
Functions: 16 16 100.0%
Branches: 162 275 58.9%
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 /*!
9 * \file
10 * \ingroup FacetCoupling
11 * \copydoc Dumux::VertexEnrichmentHelper
12 */
13 #ifndef DUMUX_VERTEX_ENRICHMENT_HELPER_HH
14 #define DUMUX_VERTEX_ENRICHMENT_HELPER_HH
15
16 #include <cmath>
17 #include <vector>
18 #include <unordered_map>
19 #include <unordered_set>
20
21 #include <dune/common/exceptions.hh>
22 #include <dune/common/reservedvector.hh>
23
24 #include <dune/grid/common/mcmgmapper.hh>
25
26 #include <dumux/common/math.hh>
27 #include <dumux/common/indextraits.hh>
28
29 namespace Dumux {
30
31 /*!
32 * \ingroup FacetCoupling
33 * \brief Specialization of the enrichment helper class for 2d grids.
34 * In this case, we look for two-dimensional bulk grid elements that
35 * are enclosed by (lie in between) two 1-dimensional facet grid elements.
36 *
37 * \tparam GridView The grid view of the domain for which nodal dofs should be enriched.
38 * \tparam CodimOneGridView The grid view of a (dim-1)-dimensional grid conforming
39 * with the facets of this grid view, indicating on which facets
40 * nodal dofs should be enriched.
41 */
42 template< class GridView, class CodimOneGridView >
43 class VertexEnrichmentHelper
44 {
45 static constexpr int dim = GridView::dimension;
46 static constexpr int dimWorld = GridView::dimensionworld;
47 static_assert(dim == 2 || dim == 3, "Grid dimension must be two or three");
48 static_assert(dimWorld == int(CodimOneGridView::dimensionworld), "world dimension mismatch");
49
50 using Intersection = typename GridView::Intersection;
51 using MCMGMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
52 using GridIndexType = typename IndexTraits<GridView>::GridIndex;
53 using Element = typename GridView::template Codim<0>::Entity;
54 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
55
56 using ElementPath = std::vector< GridIndexType >;
57 using NodalElementPaths = std::vector< std::vector<ElementPath> >;
58
59 public:
60
61 /*!
62 * \brief Enriches the dof map subject to a (dim-1)-dimensional grid.
63 * \note This assumes conforming grids and assumes the index map to be
64 * initialized for the bulk grid already!
65 *
66 * \param indexMap The index map which is to be updated
67 * \param vertexMarkers Markers if vertices should be enriched
68 * \param gridView A view on the grid for which vertices should be enriched
69 * \param vertexMapper Maps vertices of the given grid view
70 * \param elementMapper Maps elements of the given grid view
71 * \param codimOneGridView The view on the facet-conforming grid of codim 1
72 * \param codimOneGridAdapter Adapter class that allows access to information on the d-
73 * dimensional grid for entities of the (d-1)-dimensional grid
74 * \return the number of dofs after enrichment
75 */
76 template< class IndexMap, class CodimOneGridAdapter >
77 28 static std::size_t enrich(IndexMap& indexMap,
78 const std::vector<bool>& vertexMarkers,
79 const GridView& gridView,
80 const MCMGMapper& vertexMapper,
81 const MCMGMapper& elementMapper,
82 const CodimOneGridView& codimOneGridView,
83 const CodimOneGridAdapter& codimOneGridAdapter)
84 {
85 // determine the bulk element paths around vertices marked for enrichment
86
1/2
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
 28 NodalElementPaths nodalPaths(gridView.size(dim));
87
8/11
✓ Branch 1 taken 603 times.
✓ Branch 2 taken 1 times.
✓ Branch 4 taken 25 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 53705 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 53705 times.
✓ Branch 10 taken 24 times.
✓ Branch 14 taken 24 times.
✗ Branch 15 not taken.
✓ Branch 3 taken 491 times.
 112726 for (const auto& e : elements(gridView))
88 {
89
1/2
✓ Branch 1 taken 54196 times.
✗ Branch 2 not taken.
 56299 const auto eIdx = elementMapper.index(e);
90
91
1/2
✓ Branch 1 taken 53793 times.
✗ Branch 2 not taken.
 56299 std::vector<unsigned int> handledFacets;
92
1/2
✓ Branch 1 taken 491 times.
✗ Branch 2 not taken.
 56299 const auto refElement = referenceElement(e);
93
13/16
✓ Branch 1 taken 54284 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 491 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2751 times.
✓ Branch 8 taken 204200 times.
✓ Branch 10 taken 259378 times.
✓ Branch 11 taken 32 times.
✓ Branch 12 taken 2719 times.
✗ Branch 13 not taken.
✓ Branch 17 taken 491 times.
✓ Branch 18 taken 435 times.
✓ Branch 19 taken 240 times.
✓ Branch 20 taken 1724 times.
✓ Branch 9 taken 2839 times.
✓ Branch 16 taken 56 times.
 1082287 for (const auto& is : intersections(gridView, e))
94 {
95 // do not start path search on boundary intersections
96
2/2
✓ Branch 0 taken 4824 times.
✓ Branch 1 taken 203600 times.
 216428 if (!is.neighbor())
97 5664 continue;
98
99 // if facet has been handled already, skip rest (necessary for e.g. Dune::FoamGrid)
100
4/5
✗ Branch 0 not taken.
✓ Branch 1 taken 203600 times.
✓ Branch 2 taken 32 times.
✓ Branch 3 taken 232 times.
✓ Branch 4 taken 201612 times.
 210828 if (std::find(handledFacets.begin(), handledFacets.end(), is.indexInInside()) != handledFacets.end())
101 64 continue;
102
103 // first, get indices of all vertices of this intersection
104
2/3
✓ Branch 2 taken 201844 times.
✓ Branch 3 taken 1724 times.
✗ Branch 4 not taken.
 210764 const auto numCorners = is.geometry().corners();
105
1/2
✓ Branch 1 taken 203568 times.
✗ Branch 2 not taken.
 210764 std::vector<GridIndexType> faceVertexIndices(numCorners);
106
2/2
✓ Branch 0 taken 419272 times.
✓ Branch 1 taken 203568 times.
 651392 for (int i = 0; i < numCorners; ++i)
107
1/2
✓ Branch 2 taken 419272 times.
✗ Branch 3 not taken.
 440628 faceVertexIndices[i] = vertexMapper.subIndex( e,
108 refElement.subEntity(is.indexInInside(), 1, i, dim),
109 dim );
110
111 // start path search on intersections that lie on a facet grid element
112
3/4
✓ Branch 1 taken 203568 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1896 times.
✓ Branch 4 taken 201672 times.
 210764 if (codimOneGridAdapter.composeFacetElement(faceVertexIndices))
113 {
114
1/4
✓ Branch 1 taken 1896 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 2088 handledFacets.push_back(is.indexInInside());
115
2/2
✓ Branch 0 taken 4160 times.
✓ Branch 1 taken 1896 times.
 6808 for (int i = 0; i < numCorners; ++i)
116 {
117 // set up path only around those vertices that are to be enriched
118
2/2
✓ Branch 0 taken 166 times.
✓ Branch 1 taken 3994 times.
 4720 const auto vIdxGlobal = faceVertexIndices[i];
119
2/2
✓ Branch 0 taken 166 times.
✓ Branch 1 taken 3994 times.
 4720 if (!vertexMarkers[vIdxGlobal])
120 198 continue;
121
122 // construct the path only if element is not yet contained in another path
123 4522 bool found = false;
124
4/4
✓ Branch 0 taken 2050 times.
✓ Branch 1 taken 2294 times.
✓ Branch 2 taken 4344 times.
✓ Branch 3 taken 1700 times.
 6864 for (const auto& path : nodalPaths[vIdxGlobal])
125
2/2
✓ Branch 0 taken 2050 times.
✓ Branch 1 taken 2294 times.
 5062 if (std::find(path.begin(), path.end(), eIdx) != path.end())
126 { found = true; break; }
127
128
2/2
✓ Branch 0 taken 1700 times.
✓ Branch 1 taken 2294 times.
 4522 if (!found)
129 {
130 1802 ElementPath path;
131
132 // Reserve enough memory so that reallocation during the recursive search
133 // does not happen and references are invalidated. Memory is not an issue
134 // here as after enrichment the container is thrown away again.
135 // TODO improve this!?
136
1/2
✓ Branch 1 taken 1700 times.
✗ Branch 2 not taken.
 1802 path.reserve(150);
137
1/2
✓ Branch 1 taken 1700 times.
✗ Branch 2 not taken.
 1802 path.push_back(eIdx);
138
1/2
✓ Branch 1 taken 1700 times.
✗ Branch 2 not taken.
 1802 continuePathSearch_(path, gridView, elementMapper, vertexMapper, codimOneGridAdapter, e, refElement, is, vIdxGlobal);
139
1/2
✓ Branch 1 taken 1700 times.
✗ Branch 2 not taken.
 1802 nodalPaths[vIdxGlobal].emplace_back(std::move(path));
140 1802 }
141 }
142 }
143 }
144 }
145
146 // determine the offsets for each bulk vertex index on the basis of the paths found per vertex
147
4/7
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 1 taken 25 times.
✓ Branch 4 taken 25 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 25 times.
✗ Branch 8 not taken.
 28 std::vector<std::size_t> bulkVertexIndexOffsets(gridView.size(dim), 0);
148
7/9
✓ Branch 1 taken 234 times.
✓ Branch 2 taken 1 times.
✓ Branch 4 taken 25 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 48985 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 24 times.
✓ Branch 10 taken 48985 times.
✓ Branch 3 taken 153 times.
 99659 for (const auto& v : vertices(gridView))
149 {
150
1/2
✓ Branch 1 taken 153 times.
✗ Branch 2 not taken.
 49778 const auto vIdx = vertexMapper.index(v);
151
2/2
✓ Branch 0 taken 855 times.
✓ Branch 1 taken 48340 times.
 49778 if (vertexMarkers[vIdx])
152 909 bulkVertexIndexOffsets[vIdx] = nodalPaths[vIdx].size()-1;
153 }
154
155 // ... and accumulate the offsets
156 28 std::size_t sumOffset = 0;
157 28 std::size_t size = 0;
158
2/2
✓ Branch 0 taken 49195 times.
✓ Branch 1 taken 26 times.
 49806 for (auto& nodalOffset : bulkVertexIndexOffsets)
159 {
160 49778 const auto os = nodalOffset;
161 49778 nodalOffset = sumOffset;
162 49778 sumOffset += os;
163
2/2
✓ Branch 0 taken 815 times.
✓ Branch 1 taken 48380 times.
 49778 size += (os == 0) ? 1 : os + 1;
164 }
165
166 // Now, finally set up the new index map
167
8/11
✓ Branch 1 taken 603 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 603 times.
✓ Branch 5 taken 1 times.
✓ Branch 7 taken 53706 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 53705 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 53705 times.
✓ Branch 13 taken 24 times.
✓ Branch 6 taken 491 times.
 168333 for (const auto& e : elements(gridView))
168 {
169 56299 const auto& eg = e.geometry();
170
1/2
✓ Branch 1 taken 54196 times.
✗ Branch 2 not taken.
 56299 const auto eIdx = elementMapper.index(e);
171
3/3
✓ Branch 1 taken 54460 times.
✓ Branch 2 taken 88 times.
✓ Branch 0 taken 208128 times.
 272663 for (int i = 0; i < eg.corners(); ++i)
172 {
173
1/2
✓ Branch 1 taken 208392 times.
✗ Branch 2 not taken.
 216364 const auto origVIdx = vertexMapper.subIndex(e, i, dim);
174
175 // it the node itself is not enriched, simply add offset
176
2/2
✓ Branch 0 taken 203124 times.
✓ Branch 1 taken 5268 times.
 216364 if (!vertexMarkers[origVIdx])
177 210052 indexMap[eIdx][i] += bulkVertexIndexOffsets[origVIdx];
178
179 // find the local index of the path the element belongs to
180 else
181 {
182 6312 bool found = false;
183 6312 const auto& paths = nodalPaths[ origVIdx ];
184
1/2
✓ Branch 0 taken 7617 times.
✗ Branch 1 not taken.
 9041 for (int pathIdx = 0; pathIdx < paths.size(); ++pathIdx)
185 {
186
2/2
✓ Branch 0 taken 5268 times.
✓ Branch 1 taken 2349 times.
 9041 const auto& curPath = paths[pathIdx];
187
2/2
✓ Branch 0 taken 5268 times.
✓ Branch 1 taken 2349 times.
 9041 if ( std::find(curPath.begin(), curPath.end(), eIdx) != curPath.end() )
188 {
189 6312 indexMap[eIdx][i] += bulkVertexIndexOffsets[origVIdx] + pathIdx;
190 found = true;
191 break;
192 }
193 }
194
195 if (!found)
196
0/20
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
 216364 DUNE_THROW(Dune::InvalidStateException, "Element not found in any path");
197 }
198 }
199 }
200
201 // return the number of dofs after enrichment
202 28 return size;
203 28 }
204
205 private:
206 template< class ReferenceElement, class CodimOneGridAdapter >
207 6312 static void continuePathSearch_(ElementPath& path,
208 const GridView& gridView,
209 const MCMGMapper& elementMapper,
210 const MCMGMapper& vertexMapper,
211 const CodimOneGridAdapter& codimOneGridAdapter,
212 const Element& element,
213 const ReferenceElement& refElement,
214 const Intersection& prevIntersection,
215 GridIndexType vIdxGlobal)
216 {
217 // on 2d/3d grids, we need to find 1/2 more intersections at vertex
218 static constexpr int numIsToFind = dim == 3 ? 2 : 1;
219
220 // keep track of facets handled already while searching
221 6312 unsigned int foundCounter = 0;
222 6312 std::vector<unsigned int> handledFacets;
223
13/19
✓ Branch 1 taken 5268 times.
✗ Branch 2 not taken.
✓ Branch 9 taken 148 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 302 times.
✓ Branch 12 taken 13363 times.
✓ Branch 16 taken 1463 times.
✓ Branch 17 taken 1595 times.
✓ Branch 6 taken 16677 times.
✓ Branch 7 taken 324 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 18952 times.
✓ Branch 15 taken 1091 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 280 times.
✓ Branch 4 taken 324 times.
✗ Branch 5 not taken.
✗ Branch 8 not taken.
✓ Branch 20 taken 929 times.
 89003 for (const auto& is : intersections(gridView, element))
224 {
225 // skip if on previous intersection again
226
5/5
✓ Branch 0 taken 280 times.
✓ Branch 1 taken 17786 times.
✓ Branch 2 taken 72 times.
✓ Branch 3 taken 4516 times.
✓ Branch 4 taken 12269 times.
 21919 if (is.indexInInside() == prevIntersection.indexInInside())
227 13742 continue;
228
229 // skip intersection if handled already (necessary for e.g. Dune::FoamGrid)
230
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 13282 times.
 32510 if (std::count(handledFacets.begin(), handledFacets.end(), is.indexInInside()))
231 48 continue;
232
233 // determine all vertex indices of this face
234
2/3
✓ Branch 2 taken 12353 times.
✓ Branch 3 taken 929 times.
✗ Branch 4 not taken.
 16207 const auto numCorners = is.geometry().corners();
235
1/2
✓ Branch 1 taken 13282 times.
✗ Branch 2 not taken.
 16207 std::vector<GridIndexType> faceVertexIndices(numCorners);
236
2/2
✓ Branch 0 taken 32218 times.
✓ Branch 1 taken 13282 times.
 57116 for (int i = 0; i < numCorners; ++i)
237
1/2
✓ Branch 2 taken 32218 times.
✗ Branch 3 not taken.
 40909 faceVertexIndices[i] = vertexMapper.subIndex( element,
238 refElement.subEntity(is.indexInInside(), 1, i, dim),
239 dim );
240
241 // we found another intersection at the vertex if it contains the vertex around which we rotate
242
2/2
✓ Branch 0 taken 7232 times.
✓ Branch 1 taken 6050 times.
 16207 if (std::find(faceVertexIndices.begin(), faceVertexIndices.end(), vIdxGlobal) != faceVertexIndices.end())
243 {
244
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 7232 times.
✗ Branch 2 not taken.
 9268 foundCounter++;
245
3/6
✓ Branch 1 taken 7232 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 692 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 9268 handledFacets.push_back(is.indexInInside());
246
247 // keep searching in the outside element only if ...
248 // ... this is a not (processor) boundary
249
2/2
✓ Branch 0 taken 62 times.
✓ Branch 1 taken 7170 times.
 9268 if (!is.neighbor())
250
1/2
✓ Branch 0 taken 268 times.
✗ Branch 1 not taken.
 2790 continue;
251
252 // ... this face does not lie on the facet grid
253
3/4
✓ Branch 1 taken 7170 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2294 times.
✓ Branch 4 taken 4876 times.
 9198 if (codimOneGridAdapter.composeFacetElement(faceVertexIndices))
254 2720 continue;
255
256 // ... outside element is not contained yet in path
257 6478 const auto outsideElement = is.outside();
258
3/4
✓ Branch 1 taken 4876 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3536 times.
✓ Branch 4 taken 1308 times.
 6478 const auto outsideElemIdx = elementMapper.index(outsideElement);
259
2/2
✓ Branch 0 taken 3568 times.
✓ Branch 1 taken 1308 times.
 6478 if (std::find(path.begin(), path.end(), outsideElemIdx) == path.end())
260 {
261
1/2
✓ Branch 1 taken 3270 times.
✗ Branch 2 not taken.
 4510 const auto idxInOutside = is.indexInOutside();
262
1/2
✓ Branch 1 taken 298 times.
✗ Branch 2 not taken.
 4510 const auto outsideRefElement = referenceElement(outsideElement);
263
1/2
✓ Branch 1 taken 3568 times.
✗ Branch 2 not taken.
 4510 path.push_back(outsideElemIdx);
264
265 // on 2d grids, there is only going to be one more
266
11/18
✓ Branch 1 taken 3568 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 298 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 298 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 202 times.
✓ Branch 12 taken 8646 times.
✓ Branch 14 taken 12877 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 1207 times.
✗ Branch 17 not taken.
✓ Branch 21 taken 298 times.
✓ Branch 22 taken 729 times.
✓ Branch 9 taken 52 times.
✗ Branch 10 not taken.
✗ Branch 13 not taken.
✓ Branch 6 taken 10939 times.
 57652 for (const auto& outsideIs : intersections(gridView, outsideElement))
267 {
268
2/2
✓ Branch 0 taken 3568 times.
✓ Branch 1 taken 8482 times.
 15245 if (outsideIs.indexInInside() == idxInOutside)
269 {
270 // if this is the last intersection to find, return
271
2/2
✓ Branch 0 taken 2443 times.
✓ Branch 1 taken 1125 times.
 4510 if (foundCounter == numIsToFind)
272
1/2
✓ Branch 1 taken 2443 times.
✗ Branch 2 not taken.
 2846 return continuePathSearch_(path,
273 gridView,
274 elementMapper,
275 vertexMapper,
276 codimOneGridAdapter,
277 outsideElement,
278 outsideRefElement,
279 outsideIs,
280 2846 vIdxGlobal);
281
282 // otherwise, put one more search on stack
283 else
284
1/2
✓ Branch 1 taken 1125 times.
✗ Branch 2 not taken.
 1664 continuePathSearch_(path,
285 gridView,
286 elementMapper,
287 vertexMapper,
288 codimOneGridAdapter,
289 outsideElement,
290 outsideRefElement,
291 outsideIs,
292 vIdxGlobal);
293 }
294 }
295 }
296 5900 }
297 }
298
299
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2825 times.
 3466 if (foundCounter != numIsToFind)
300
1/30
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✓ Branch 45 taken 428 times.
✗ Branch 46 not taken.
 428 DUNE_THROW(Dune::InvalidStateException, "Found " << foundCounter << " instead of " << numIsToFind << " intersections around vertex");
301 6312 }
302 };
303
304 } // end namespace Dumux
305
306 #endif
307