GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/multidomain/facet/codimonegridadapter.hh
Date: 2025-06-14 19:21:29
Exec Total Coverage
Lines: 68 69 98.6%
Functions: 63 63 100.0%
Branches: 67 108 62.0%
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::CodimOneGridAdapter
12 */
13 #ifndef DUMUX_FACETCOUPLING_CODIM_ONE_GRID_ADAPTER_HH
14 #define DUMUX_FACETCOUPLING_CODIM_ONE_GRID_ADAPTER_HH
15
16 #include <cassert>
17 #include <vector>
18 #include <memory>
19
20 #include <dune/grid/common/mcmgmapper.hh>
21
22 #include <dumux/common/indextraits.hh>
23
24 namespace Dumux {
25
26 /*!
27 * \ingroup FacetCoupling
28 * \brief Adapter that allows retrieving information on a d-dimensional grid for
29 * entities of a (d-1)-dimensional grid. This lower-dimensional grid is
30 * assumed to be facet-conforming to the d-dimensional grid. This class
31 * can be used in the context of models where a sub-domain lives on the
32 * facets of a bulk grid.
33 *
34 * \tparam Embeddings Class containing the embedments of entities between grids of codim 1
35 * \tparam bulkGridId The grid id of the d-dimensional grid within the hierarchy
36 * \tparam facetGridId The grid id of the (d-1)-dimensional grid within the hierarchy
37 */
38 template<class Embeddings, int bulkGridId = 0, int facetGridId = 1>
39 class CodimOneGridAdapter
40 {
41 // Extract some types of the facet-conforming grid of codimension one
42 using FacetGridView = typename Embeddings::template GridView<facetGridId>;
43 using FacetGridVertex = typename FacetGridView::template Codim<FacetGridView::dimension>::Entity;
44 using FacetGridElement = typename FacetGridView::template Codim<0>::Entity;
45 using FacetGridIndexType = typename IndexTraits<FacetGridView>::GridIndex;
46
47 // Extract some types of the bulk grid
48 using BulkGridView = typename Embeddings::template GridView<bulkGridId>;
49 using BulkMapper = Dune::MultipleCodimMultipleGeomTypeMapper<BulkGridView>;
50 using BulkGridElement = typename BulkGridView::template Codim<0>::Entity;
51 using BulkGridIntersection = typename BulkGridView::Intersection;
52 using BulkGridVertex = typename BulkGridView::template Codim<BulkGridView::dimension>::Entity;
53 using BulkIndexType = typename IndexTraits<BulkGridView>::GridIndex;
54
55 // check if provided id combination makes sense
56 static_assert( int(FacetGridView::dimension) == int(BulkGridView::dimension) - 1,
57 "Grid dimension mismatch! Please check the provided domain ids!" );
58 static_assert( int(FacetGridView::dimensionworld) == int(BulkGridView::dimensionworld),
59 "Grid world dimension mismatch! All grids must have the same world dimension" );
60
61 public:
62
63 //! The constructor
64 83 CodimOneGridAdapter(std::shared_ptr<const Embeddings> embeddings)
65 83 : embeddingsPtr_(embeddings)
66
2/4
✓ Branch 1 taken 79 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 77 times.
✗ Branch 5 not taken.
 166 , bulkVertexMapper_(embeddings->template gridView<bulkGridId>(), Dune::mcmgVertexLayout())
67 {
68 // bulk insertion to grid index map
69
1/2
✓ Branch 1 taken 77 times.
✗ Branch 2 not taken.
 83 const auto& bulkGridView = embeddings->template gridView<bulkGridId>();
70
3/5
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 1 taken 77 times.
✓ Branch 4 taken 77 times.
✗ Branch 5 not taken.
 83 bulkInsertionToGridVIdx_.resize(bulkGridView.size(BulkGridView::dimension));
71
7/9
✓ Branch 1 taken 3467 times.
✓ Branch 2 taken 2 times.
✓ Branch 4 taken 77 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 129033 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 68 times.
✓ Branch 10 taken 129033 times.
✓ Branch 3 taken 3285 times.
 267489 for (const auto& v : vertices(bulkGridView))
72
5/8
✓ Branch 1 taken 132318 times.
✓ Branch 2 taken 114 times.
✓ Branch 4 taken 3285 times.
✓ Branch 5 taken 129033 times.
✓ Branch 7 taken 3285 times.
✗ Branch 8 not taken.
✗ Branch 3 not taken.
✗ Branch 6 not taken.
 133598 bulkInsertionToGridVIdx_[embeddings->template insertionIndex<bulkGridId>(v)] = bulkVertexMapper_.index(v);
73
74 // Determine map from the hierarchy's vertex idx to bulk insertion idx
75 // There is one unique set of vertex indices within the hierarchy.
76 // Obtain the hierarchy indices that make up the bulk grid. These
77 // are ordered corresponding to their insertion (thus loopIdx = insertionIdx)
78
1/2
✓ Branch 1 taken 79 times.
✗ Branch 2 not taken.
 83 hierarchyToBulkInsertionIdx_.resize(embeddingsPtr_->numVerticesInHierarchy());
79
1/2
✓ Branch 1 taken 79 times.
✗ Branch 2 not taken.
 83 bulkGridHasHierarchyVertex_.resize(embeddingsPtr_->numVerticesInHierarchy(), false);
80 83 const auto& bulkHierarchyIndices = embeddingsPtr_->gridHierarchyIndices(bulkGridId);
81
2/2
✓ Branch 0 taken 132432 times.
✓ Branch 1 taken 79 times.
 133681 for (std::size_t insIdx = 0; insIdx < bulkHierarchyIndices.size(); ++insIdx)
82 {
83 133598 hierarchyToBulkInsertionIdx_[ bulkHierarchyIndices[insIdx] ] = insIdx;
84 133598 bulkGridHasHierarchyVertex_[ bulkHierarchyIndices[insIdx] ] = true;
85 }
86
87 // determine which bulk vertices lie on facet elements
88
3/5
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 1 taken 77 times.
✓ Branch 4 taken 77 times.
✗ Branch 5 not taken.
 83 bulkVertexIsOnFacetGrid_.resize(bulkGridView.size(BulkGridView::dimension), false);
89 83 const auto& facetGridView = embeddings->template gridView<facetGridId>();
90
3/4
✓ Branch 1 taken 2919 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2919 times.
✓ Branch 5 taken 79 times.
 6169 for (const auto& v : vertices(facetGridView))
91 {
92
2/4
✓ Branch 1 taken 2919 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2919 times.
✗ Branch 4 not taken.
 3043 const auto insIdx = embeddings->template insertionIndex<facetGridId>(v);
93
1/2
✓ Branch 0 taken 2919 times.
✗ Branch 1 not taken.
 3043 const auto hierarchyInsIdx = embeddings->gridHierarchyIndices(facetGridId)[insIdx];
94
95
1/2
✓ Branch 0 taken 2919 times.
✗ Branch 1 not taken.
 3043 if (bulkGridHasHierarchyVertex_[hierarchyInsIdx])
96 3043 bulkVertexIsOnFacetGrid_[ getBulkGridVertexIndex_(hierarchyInsIdx) ] = true;
97 }
98
99 // determine the bulk vertex indices that make up facet elements & connectivity
100
1/2
✓ Branch 2 taken 79 times.
✗ Branch 3 not taken.
 83 facetElementCorners_.resize(facetGridView.size(0));
101
3/5
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 1 taken 77 times.
✓ Branch 4 taken 77 times.
✗ Branch 5 not taken.
 83 facetElementsAtBulkVertex_.resize(bulkGridView.size(BulkGridView::dimension));
102
103 83 std::size_t facetElementCounter = 0;
104
2/2
✓ Branch 0 taken 3140 times.
✓ Branch 1 taken 79 times.
 3407 for (const auto& element : elements(facetGridView))
105 {
106
1/2
✓ Branch 0 taken 3140 times.
✗ Branch 1 not taken.
 3324 if (isEmbedded(element))
107 {
108 // obtain the bulk vertex indices of the corners of this element
109
1/2
✓ Branch 1 taken 3140 times.
✗ Branch 2 not taken.
 3324 const auto numCorners = element.subEntities(FacetGridView::dimension);
110
1/2
✓ Branch 1 taken 3140 times.
✗ Branch 2 not taken.
 3324 std::vector<BulkIndexType> cornerIndices(numCorners);
111
2/2
✓ Branch 0 taken 7138 times.
✓ Branch 1 taken 3140 times.
 11006 for (int i = 0; i < numCorners; ++i)
112
1/2
✓ Branch 1 taken 7138 times.
✗ Branch 2 not taken.
 7682 cornerIndices[i] = bulkGridVertexIndex(element.template subEntity<FacetGridView::dimension>(i));
113
114 // update connectivity map facetVertex -> facetElements
115
2/2
✓ Branch 0 taken 7138 times.
✓ Branch 1 taken 3140 times.
 11006 for (auto bulkVIdx : cornerIndices)
116
1/4
✓ Branch 1 taken 7138 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 7682 facetElementsAtBulkVertex_[bulkVIdx].push_back(facetElementCounter);
117
118 // update facet elements (identified by corners - store them sorted!)
119 3324 std::sort(cornerIndices.begin(), cornerIndices.end());
120
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 320 times.
 3324 facetElementCorners_[facetElementCounter] = std::move(cornerIndices);
121 3324 }
122
123 3324 facetElementCounter++;
124 }
125 83 }
126
127 /*!
128 * \brief Returns the index within the d-dimensional grid of a vertex
129 * of the (d-1)-dimensional grid.
130 * \note Leads to undefined behaviour if called for a vertex which doesn't
131 * exist on the d-dimensional grid
132 */
133 13954 BulkIndexType bulkGridVertexIndex(const FacetGridVertex& v) const
134 {
135 13954 const auto insIdx = embeddingsPtr_->template insertionIndex<facetGridId>(v);
136 13954 const auto hierarchyInsIdx = embeddingsPtr_->gridHierarchyIndices(facetGridId)[insIdx];
137 13954 return getBulkGridVertexIndex_(hierarchyInsIdx);
138 }
139
140 /*!
141 * \brief Returns true if the vertex of the d-dimensional grid with
142 * the given vertex index also exists on the (d-1)-dimensional grid
143 */
144 4460 bool isOnFacetGrid(const BulkGridVertex& v) const
145 {
146 4460 const auto bulkInsIdx = embeddingsPtr_->template insertionIndex<bulkGridId>(v);
147 4460 const auto bulkVIdx = bulkInsertionToGridVIdx_[bulkInsIdx];
148 4460 return bulkVertexIsOnFacetGrid_[bulkVIdx];
149 }
150
151 /*!
152 * \brief Returns true if the given intersection coincides with a facet grid
153 *
154 * \param element An element of the bulk grid
155 * \param intersection An bulk grid intersection within the given element
156 */
157 bool isOnFacetGrid(const BulkGridElement& element, const BulkGridIntersection& intersection) const
158 {
159 // Intersection lies on facet grid, if the corners of the intersection make up a facet element
160 const auto refElement = referenceElement(element);
161 const auto numCorners = intersection.geometry().corners();
162 const auto facetIdx = intersection.indexInInside();
163
164 std::vector<BulkIndexType> cornerIndices(numCorners);
165 for (int i = 0; i < numCorners; ++i)
166 cornerIndices[i] = bulkVertexMapper_.subIndex( element,
167 refElement.subEntity(facetIdx, 1, i, BulkGridView::dimension),
168 BulkGridView::dimension );
169
170 return composeFacetElement(cornerIndices);
171 }
172
173 /*!
174 * \brief Returns true if a given set of bulk vertex indices make up a facet grid element
175 * \note The vertex indices are NOT the dof indices in the context of models where there
176 * are multiple dofs at one vertex (enriched nodal dofs). Here, we expect the vertex
177 * indices (unique index per vertex).
178 */
179 template<class IndexStorage>
180 513970 bool composeFacetElement(const IndexStorage& bulkVertexIndices) const
181 {
182 // set up a vector containing all element indices the vertices are connected to
183 513970 std::vector<std::size_t> facetElemIndices;
184
2/2
✓ Branch 0 taken 1073424 times.
✓ Branch 1 taken 496774 times.
 1638442 for (auto bulkVIdx : bulkVertexIndices)
185 1124472 facetElemIndices.insert( facetElemIndices.end(),
186
1/2
✓ Branch 1 taken 1073424 times.
✗ Branch 2 not taken.
 1124472 facetElementsAtBulkVertex_[bulkVIdx].begin(),
187
1/2
✓ Branch 1 taken 1073424 times.
✗ Branch 2 not taken.
 1124472 facetElementsAtBulkVertex_[bulkVIdx].end() );
188
189 // if no facet elements are connected to the vertices this is not on facet grid
190
2/2
✓ Branch 0 taken 44488 times.
✓ Branch 1 taken 452286 times.
 513970 if (facetElemIndices.size() == 0)
191 return false;
192
193 // make the container unique
194 51084 std::sort(facetElemIndices.begin(), facetElemIndices.end());
195 51084 facetElemIndices.erase(std::unique(facetElemIndices.begin(), facetElemIndices.end()), facetElemIndices.end());
196
197 // check if given indices make up a facet element
198
1/2
✓ Branch 1 taken 44488 times.
✗ Branch 2 not taken.
 51084 auto cornerIndexCopy = bulkVertexIndices;
199 51084 std::sort(cornerIndexCopy.begin(), cornerIndexCopy.end());
200
2/2
✓ Branch 0 taken 210916 times.
✓ Branch 1 taken 36250 times.
 294585 for (const auto& facetElemIdx : facetElemIndices)
201 {
202
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 210916 times.
 252549 const auto& facetElemCorners = facetElementCorners_[facetElemIdx];
203
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 210916 times.
 252549 if (facetElemCorners.size() != cornerIndexCopy.size())
204 ✗ continue;
205
206
3/4
✓ Branch 0 taken 210916 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 202678 times.
✓ Branch 3 taken 8238 times.
 252549 if ( std::equal(cornerIndexCopy.begin(), cornerIndexCopy.end(),
207 facetElemCorners.begin(), facetElemCorners.end()) )
208
1/2
✓ Branch 0 taken 5024 times.
✗ Branch 1 not taken.
 51084 return true;
209 }
210
211 // no corresponding facet element found
212 return false;
213 513970 }
214
215 /*!
216 * \brief Returns true if a (d-1)-dimensional element is embedded in
217 * the d-dimensional domain
218 */
219 3140 bool isEmbedded(const FacetGridElement& e) const
220
2/4
✓ Branch 1 taken 2964 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 176 times.
✗ Branch 5 not taken.
 3140 { return numEmbedments(e) > 0; }
221
222 /*!
223 * \brief Returns the number of d-dimensional elements in which the
224 * given (d-1)-dimensional element is embedded in
225 */
226 4408 std::size_t numEmbedments(const FacetGridElement& e) const
227
1/2
✓ Branch 1 taken 384 times.
✗ Branch 2 not taken.
 4408 { return embeddingsPtr_->template adjoinedEntityIndices<facetGridId>(e).size(); }
228
229 private:
230 //! Obtains the bulk grid vertex index from a given insertion index on the hierarchy
231 15719 BulkIndexType getBulkGridVertexIndex_(BulkIndexType hierarchyInsertionIdx) const
232 {
233
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14441 times.
 15719 assert(bulkGridHasHierarchyVertex_[hierarchyInsertionIdx]);
234 15719 return bulkInsertionToGridVIdx_[ hierarchyToBulkInsertionIdx_[hierarchyInsertionIdx] ];
235 }
236
237 // shared pointer to the embedment data
238 std::shared_ptr<const Embeddings> embeddingsPtr_;
239
240 // vertex mapper of the bulk grid
241 BulkMapper bulkVertexMapper_;
242
243 // data stored on grid vertices
244 std::vector<bool> bulkVertexIsOnFacetGrid_;
245 std::vector<BulkIndexType> bulkInsertionToGridVIdx_;
246 std::vector<BulkIndexType> hierarchyToBulkInsertionIdx_;
247 std::vector<bool> bulkGridHasHierarchyVertex_;
248
249 // data stored for elements on the codim one grid
250 std::vector< std::vector<BulkIndexType> > facetElementsAtBulkVertex_;
251 std::vector< std::vector<BulkIndexType> > facetElementCorners_;
252 };
253
254 } // end namespace Dumux
255
256 #endif
257