GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/multidomain/facet/vertexmapper.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 67 72 93.1%
Functions: 47 81 58.0%
Branches: 130 271 48.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-FileCopyrightInfo: 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::EnrichedVertexDofMapper
12 */
13 #ifndef DUMUX_ENRICHED_VERTEX_DOF_MAPPER_HH
14 #define DUMUX_ENRICHED_VERTEX_DOF_MAPPER_HH
15
16 #include <vector>
17
18 #include <dune/common/exceptions.hh>
19 #include <dune/common/timer.hh>
20 #include <dune/grid/common/mcmgmapper.hh>
21
22 #include "enrichmenthelper.hh"
23
24 namespace Dumux {
25
26 /*!
27 * \ingroup FacetCoupling
28 * \brief An indicator class used to mark vertices for enrichment. This
29 * implementation marks all vertices of a given grid of codimension
30 * one for enrichment, except those that are connected to inmersed
31 * boundaries.
32 */
33 class EnrichmentIndicator
34 {
35
36 public:
37 /*!
38 * \brief Function that marks vertices for enrichment. This implementation
39 * works on the basis of a facet-conforming grid of codimension one
40 * which is used to determine the vertices that should be enriched.
41 * Effectively, all vertices that lie on the given (d-1)-dimensional
42 * grid are marked, except those that are connected to inmersed boundaries.
43 *
44 * \param vertexMarkers Stores for each vertex if it should be enriched
45 * \param gridView The d-dimensional grid for which vertices should be enriched
46 * \param vertexMapper Mapper that maps to the vertex indices of the given grid view
47 * \param codimOneGridView The view on the (d-1)-dimensional facet-conforming grid
48 * \param codimOneGridAdapter Adapter class that allows access to information on the d-
49 * dimensional grid for entities of the (d-1)-dimensional grid
50 */
51 template< class GridView,
52 class VertexMapper,
53 class CodimOneGridView,
54 class CodimOneGridAdapter >
55 28 static void markVerticesForEnrichment(std::vector<bool>& vertexMarkers,
56 const GridView& gridView,
57 const VertexMapper& vertexMapper,
58 const CodimOneGridView& codimOneGridView,
59 const CodimOneGridAdapter& codimOneGridAdapter)
60 {
61 static constexpr int dim = GridView::dimension;
62 static_assert(CodimOneGridView::dimension == dim-1, "Grid dimension mismatch");
63
64 // reset the markers
65 28 vertexMarkers.assign(gridView.size(dim), false);
66
67 // first find all bulk grid vertices on the boundary
68
4/8
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 25 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 24 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 24 times.
✗ Branch 11 not taken.
56 std::vector<bool> isOnBoundary(gridView.size(dim), false);
69
13/14
✓ Branch 1 taken 113 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 579 times.
✓ Branch 4 taken 26 times.
✓ Branch 5 taken 491 times.
✓ Branch 6 taken 89 times.
✓ Branch 7 taken 24 times.
✓ Branch 8 taken 491 times.
✓ Branch 9 taken 53705 times.
✓ Branch 10 taken 24 times.
✓ Branch 11 taken 53705 times.
✓ Branch 12 taken 24 times.
✓ Branch 14 taken 53705 times.
✗ Branch 15 not taken.
112626 for (const auto& e : elements(gridView))
70 {
71
1/2
✓ Branch 1 taken 54284 times.
✗ Branch 2 not taken.
56299 const auto refElem = referenceElement(e);
72
12/18
✓ Branch 1 taken 54284 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 491 times.
✓ Branch 5 taken 260253 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 491 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2455 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 206196 times.
✓ Branch 12 taken 264 times.
✓ Branch 14 taken 296 times.
✓ Branch 15 taken 491 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 54473 times.
✓ Branch 18 taken 1196 times.
✓ Branch 20 taken 1964 times.
✗ Branch 21 not taken.
600381 for (const auto& is : intersections(gridView, e))
73
4/4
✓ Branch 0 taken 800 times.
✓ Branch 1 taken 6012 times.
✓ Branch 2 taken 201644 times.
✓ Branch 3 taken 264 times.
217020 if (is.boundary())
74
7/8
✓ Branch 1 taken 15936 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1024 times.
✓ Branch 4 taken 32 times.
✓ Branch 5 taken 11352 times.
✓ Branch 6 taken 4584 times.
✓ Branch 7 taken 720 times.
✓ Branch 8 taken 240 times.
19968 for (int i = 0; i < is.geometry().corners(); ++i)
75
4/7
✗ Branch 0 not taken.
✓ Branch 1 taken 10392 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 720 times.
✓ Branch 4 taken 64 times.
✓ Branch 5 taken 10328 times.
✗ Branch 6 not taken.
14128 isOnBoundary[ vertexMapper.subIndex( e,
76 refElem.subEntity(is.indexInInside(), 1, i, dim),
77 26816 dim ) ] = true;
78 }
79
80 // mark all vertices on the lower-dimensional grid to be enriched (except immersed boundaries)
81
2/2
✓ Branch 1 taken 23 times.
✓ Branch 2 taken 3 times.
81 std::vector<typename GridView::IndexSet::IndexType> vertexIndicesStorage;
82
5/6
✓ Branch 1 taken 992 times.
✓ Branch 2 taken 26 times.
✓ Branch 3 taken 992 times.
✓ Branch 4 taken 26 times.
✓ Branch 6 taken 992 times.
✗ Branch 7 not taken.
1112 for (const auto& codimOneElement : elements(codimOneGridView))
83 {
84 // if a codimension one element has 2 or more embedments, we need to enrich
85
3/4
✓ Branch 1 taken 992 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 944 times.
✓ Branch 4 taken 48 times.
1084 if (codimOneGridAdapter.numEmbedments(codimOneElement) >= 2)
86 {
87
4/4
✓ Branch 0 taken 2072 times.
✓ Branch 1 taken 944 times.
✓ Branch 2 taken 2072 times.
✓ Branch 3 taken 944 times.
6760 for (int i = 0; i < codimOneElement.subEntities(dim-1); ++i)
88
2/6
✓ Branch 1 taken 2072 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2072 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
4688 vertexMarkers[ codimOneGridAdapter.bulkGridVertexIndex(codimOneElement.template subEntity<dim-1>(i)) ] = true;
89
90 // however, we have to exclude immersed boundaries
91
1/2
✓ Branch 1 taken 944 times.
✗ Branch 2 not taken.
1036 const auto refElem = referenceElement(codimOneElement);
92
5/8
✓ Branch 1 taken 944 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 3128 times.
✗ Branch 6 not taken.
✓ Branch 10 taken 126 times.
✓ Branch 11 taken 2058 times.
✓ Branch 13 taken 2184 times.
✗ Branch 14 not taken.
7828 for (const auto& intersection : intersections(codimOneGridView, codimOneElement))
93 {
94 // skip if intersection is not on boundary
95
4/4
✓ Branch 0 taken 126 times.
✓ Branch 1 taken 2058 times.
✓ Branch 2 taken 126 times.
✓ Branch 3 taken 2058 times.
4976 if (!intersection.boundary())
96 continue;
97
98 // obtain all bulk grid indices of the lower-dimensional intersection corners
99
2/4
✓ Branch 1 taken 126 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 126 times.
✗ Branch 5 not taken.
160 const auto numCorners = intersection.geometry().corners();
100
1/2
✓ Branch 1 taken 126 times.
✗ Branch 2 not taken.
160 vertexIndicesStorage.resize(numCorners);
101
2/2
✓ Branch 0 taken 206 times.
✓ Branch 1 taken 126 times.
432 for (int i = 0; i < numCorners; ++i)
102 {
103
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 206 times.
272 const auto vIdxLocal = refElem.subEntity(intersection.indexInInside(), 1, i, dim-1);
104
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 206 times.
✓ Branch 3 taken 206 times.
✗ Branch 4 not taken.
272 vertexIndicesStorage[i] = codimOneGridAdapter.bulkGridVertexIndex( codimOneElement.template subEntity<dim-1>(vIdxLocal) );
105 }
106
107 // if any of the vertices is on an immersed boundary, we must not enrich any of them
108
12/56
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 32 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ 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 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✓ Branch 20 taken 48 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 48 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 19 times.
✓ Branch 25 taken 27 times.
✓ Branch 26 taken 19 times.
✓ Branch 27 taken 27 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✓ Branch 31 taken 67 times.
✓ Branch 32 taken 25 times.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✓ Branch 48 taken 32 times.
✗ Branch 49 not taken.
✓ Branch 50 taken 32 times.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
732 if (std::any_of(vertexIndicesStorage.begin(), vertexIndicesStorage.end(), [&isOnBoundary] (auto idx) { return !isOnBoundary[idx]; }))
109 1122 std::for_each(vertexIndicesStorage.begin(), vertexIndicesStorage.end(), [&vertexMarkers] (auto idx) { vertexMarkers[idx] = false; });
110 }
111 }
112 }
113 28 }
114 };
115
116 /*!
117 * \ingroup FacetCoupling
118 * \brief A vertex mapper that allows for enrichment of nodes. Indication on where to
119 * enrich the nodes is done on the basis of a grid of codimension one living
120 * on the facets of the bulk grid.
121 *
122 * \tparam GV The Dune::GridView type
123 */
124 template<class GV>
125 class EnrichedVertexDofMapper
126 {
127 static constexpr int dim = GV::dimension;
128 static_assert(dim > 1, "Vertex dof enrichment mapper currently only works for dim > 1!");
129
130 using GIType = typename GV::IndexSet::IndexType;
131 using Vertex = typename GV::template Codim<dim>::Entity;
132 using Element = typename GV::template Codim<0>::Entity;
133 using MCMGMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GV>;
134
135 public:
136 //! export the underlying grid view type
137 using GridView = GV;
138 //! export the grid index type
139 using GridIndexType = GIType;
140
141 //! the constructor
142 28 EnrichedVertexDofMapper(const GV& gridView)
143 : gridView_(gridView)
144 , elementMapper_(gridView, Dune::mcmgElementLayout())
145
9/22
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 26 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 26 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 26 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 26 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 26 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 26 times.
✗ Branch 21 not taken.
✓ Branch 23 taken 26 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
196 , vertexMapper_(gridView, Dune::mcmgVertexLayout())
146 {
147
1/2
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
28 initialize_();
148 28 }
149
150 //! constructor taking a layout as additional argument (for compatibility)
151 26 EnrichedVertexDofMapper(const GV& gridView, Dune::MCMGLayout layout)
152 26 : EnrichedVertexDofMapper(gridView)
153 {
154
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 24 times.
52 if ( !( static_cast<bool>(layout(Dune::GeometryTypes::vertex, dim)) ) )
155 DUNE_THROW(Dune::InvalidStateException, "Vertex mapper only makes sense for vertex layout!");
156 26 }
157
158 //! map nodal subentity of codim 0 entity to the grid dof
159 20306633 GridIndexType subIndex(const Element& e, unsigned int i, unsigned int codim) const
160 {
161
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 19956491 times.
20306633 assert(codim == dim && "Only element corners can be mapped by this mapper");
162 20306633 return indexMap_[elementMapper_.index(e)][i];
163 }
164
165 //! map nodal subentity of codim 0 entity to the grid vertex index
166 GridIndexType vertexIndex(const Element& e, unsigned int i, unsigned int codim) const
167 {
168 assert(codim == dim && "Only element corners can be mapped by this mapper");
169
4/8
✓ Branch 1 taken 388576 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 37980 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 32 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 528 times.
✗ Branch 11 not taken.
427116 return vertexMapper_.subIndex(e, i, codim);
170 }
171
172 //! map nodal entity to the grid vertex index
173 GridIndexType vertexIndex(const Vertex& v) const
174 {
175 assert(Vertex::Geometry::mydimension == 0 && "Only vertices can be mapped by this mapper");
176
1/2
✓ Branch 1 taken 816 times.
✗ Branch 2 not taken.
7828 return vertexMapper_.index(v);
177 }
178
179 //! map vertex to the grid dof index
180 //! \note This is only valid if there are no enriched vertex dofs!
181 //! We therefore ask this in every call. This means quite some overhead, but
182 //! this mapper is not designed optimally for the case of no enriched nodes.
183 template< class EntityType >
184 2052 GridIndexType index(const EntityType& e) const
185 {
186
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1026 times.
2052 if (hasEnrichedVertices_)
187 DUNE_THROW(Dune::InvalidStateException, "Index map contains enriched vertex dofs. Direct mapping from vertex to index not possible.");
188
189 assert(EntityType::Geometry::mydimension == 0 && "Only vertices can be mapped by this mapper");
190 2052 return vertexMapper_.index(e);
191 }
192
193 //! returns the number of dofs managed by this mapper
194 std::size_t size() const
195 { return size_; }
196
197 //! returns true if a vertex dof had been enriched
198 2852 bool isEnriched(const Vertex& v)
199 2852 { return isEnriched_[ vertexMapper_.index(v) ]; }
200
201 //! the update here simply updates the non-enriched map
202 //! enrichment has to be done afterwards!
203 void update(const GV& gridView)
204 {
205 24 gridView_ = gridView;
206 24 initialize_();
207 }
208
209 void update(GV&& gridView)
210 {
211 gridView_ = std::move(gridView);
212 initialize_();
213 }
214
215 /*!
216 * \brief Enriches the dof map subject to a (dim-1)-dimensional grid.
217 * \note This assumes conforming grids!
218 *
219 * \param codimOneGridView The grid view of a (dim-1)-dimensional grid conforming
220 * with the facets of the grid view passed to the constructor,
221 * indicating on which facets nodal dofs should be enriched.
222 * \param codimOneGridAdapter Adapter class that allows access to information on the d-
223 * dimensional grid for entities of the (d-1)-dimensional grid
224 * \param verbose Enable/disable terminal output of the time necessary for enrichment
225 */
226 template<class CodimOneGridView, class CodimOneGridAdapter>
227 28 void enrich(const CodimOneGridView& codimOneGridView,
228 const CodimOneGridAdapter& codimOneGridAdapter,
229 bool verbose = false)
230 {
231 static const int codimOneDim = CodimOneGridView::dimension;
232 static_assert(codimOneDim == dim-1, "Grid dimension mismatch!");
233 static_assert(codimOneDim == 2 || codimOneDim == 1, "Inadmissible codimension one grid dimension");
234 static_assert(int(CodimOneGridView::dimensionworld) == int(GV::dimensionworld), "Grid world dimension mismatch");
235
236 // keep track of time
237 28 Dune::Timer watch;
238
239 // mark vertices for enrichment using the indicator
240 28 EnrichmentIndicator::markVerticesForEnrichment(isEnriched_, gridView_, vertexMapper_, codimOneGridView, codimOneGridAdapter);
241
242 // let the helper class do the enrichment of the index map
243 56 size_ = VertexEnrichmentHelper< GridView, CodimOneGridView >::enrich(indexMap_,
244 isEnriched_,
245 gridView_,
246 vertexMapper_,
247 28 elementMapper_,
248 codimOneGridView,
249 codimOneGridAdapter);
250
251 // check if new index map contains enriched dofs
252 84 hasEnrichedVertices_ = std::any_of(isEnriched_.begin(), isEnriched_.end(), [] (bool isEnriched) { return isEnriched; });
253
254
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 22 times.
28 if (verbose)
255 12 std::cout << "Vertex dof enrichment took " << watch.elapsed() << " seconds." << std::endl;
256 28 }
257
258 private:
259
260 //! initializes the mapper on the basis of the standard Dune mcmgmapper
261 54 void initialize_()
262 {
263 54 size_ = gridView_.size(dim);
264 54 hasEnrichedVertices_ = false;
265 54 indexMap_.resize(gridView_.size(0));
266 54 isEnriched_.resize(gridView_.size(dim), false);
267
11/14
✓ Branch 1 taken 1158 times.
✓ Branch 2 taken 50 times.
✓ Branch 3 taken 1158 times.
✓ Branch 4 taken 4 times.
✓ Branch 5 taken 46 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 106623 times.
✓ Branch 8 taken 46 times.
✓ Branch 9 taken 106623 times.
✓ Branch 10 taken 46 times.
✓ Branch 12 taken 106623 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 106623 times.
✗ Branch 16 not taken.
332915 for (const auto& e : elements(gridView_))
268 {
269
1/2
✓ Branch 1 taken 106623 times.
✗ Branch 2 not taken.
112793 const auto numCorners = e.geometry().corners();
270
1/2
✓ Branch 1 taken 106623 times.
✗ Branch 2 not taken.
111811 const auto eIdxGlobal = elementMapper_.index(e);
271
2/4
✓ Branch 1 taken 106623 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 106623 times.
✗ Branch 5 not taken.
223622 indexMap_[eIdxGlobal].resize(numCorners);
272
2/2
✓ Branch 0 taken 413932 times.
✓ Branch 1 taken 107781 times.
541687 for (unsigned int i = 0; i < numCorners; ++i)
273
1/2
✓ Branch 1 taken 409476 times.
✗ Branch 2 not taken.
429876 indexMap_[eIdxGlobal][i] = vertexMapper_.subIndex(e, i, dim);
274 }
275 54 }
276
277 // data members
278 std::size_t size_; //! number of dofs mapped to by this mapper
279 GV gridView_; //! the grid view
280 MCMGMapper elementMapper_; //! unmodified element mapper
281 MCMGMapper vertexMapper_; //! unmodified vertex mapper
282 bool hasEnrichedVertices_; //! keeps track of if vertices are enriched
283 std::vector<bool> isEnriched_; //! keeps track which vertices are enriched
284 std::vector< std::vector<GIType> > indexMap_; //! contains the new dof indices
285 };
286
287 } // end namespace Dumux
288
289 #endif
290