GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/discretization/staggered/freeflow/elementvolumevariables.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 32 99 32.3%
Functions: 137 327 41.9%
Branches: 78 1046 7.5%

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 StaggeredDiscretization
10 * \copydoc Dumux::StaggeredElementVolumeVariables
11 */
12 #ifndef DUMUX_DISCRETIZATION_STAGGERED_ELEMENT_VOLUMEVARIABLES_HH
13 #define DUMUX_DISCRETIZATION_STAGGERED_ELEMENT_VOLUMEVARIABLES_HH
14
15 #include <algorithm>
16 #include <cassert>
17 #include <vector>
18 #include <utility>
19
20 #include <dune/common/exceptions.hh>
21 #include <dumux/discretization/staggered/elementsolution.hh>
22 #include <dumux/common/typetraits/vector.hh>
23
24 namespace Dumux {
25
26 /*!
27 * \ingroup StaggeredDiscretization
28 * \brief Base class for the element volume variables vector for the staggered model
29 */
30 template<class GVV, bool cachingEnabled>
31 class StaggeredElementVolumeVariables
32 {};
33
34 /*!
35 * \ingroup StaggeredDiscretization
36 * \brief Class for the element volume variables vector for the staggered model.
37 Specialization in case the volume variables are stored globally.
38 */
39 template<class GVV>
40 class StaggeredElementVolumeVariables<GVV, /*cachingEnabled*/true>
41 {
42 public:
43 //! export type of the grid volume variables
44 using GridVolumeVariables = GVV;
45
46 //! export type of the volume variables
47 using VolumeVariables = typename GridVolumeVariables::VolumeVariables;
48
49 //! Constructor
50 StaggeredElementVolumeVariables(const GridVolumeVariables& gridVolVars)
51 : gridVolVarsPtr_(&gridVolVars)
52
40/80
✓ Branch 1 taken 1322 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1322 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1322 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1322 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1322 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1355444 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1355444 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 1355444 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 1355444 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 1355444 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 1355444 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 1355444 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 1355444 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 1355444 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 1355444 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 1355444 times.
✗ Branch 47 not taken.
✓ Branch 49 taken 1355444 times.
✗ Branch 50 not taken.
✓ Branch 52 taken 1355444 times.
✗ Branch 53 not taken.
✓ Branch 55 taken 1355444 times.
✗ Branch 56 not taken.
✓ Branch 58 taken 1355444 times.
✗ Branch 59 not taken.
✓ Branch 61 taken 1355444 times.
✗ Branch 62 not taken.
✓ Branch 64 taken 1355444 times.
✗ Branch 65 not taken.
✓ Branch 67 taken 1355444 times.
✗ Branch 68 not taken.
✓ Branch 70 taken 1355444 times.
✗ Branch 71 not taken.
✓ Branch 73 taken 1355444 times.
✗ Branch 74 not taken.
✓ Branch 76 taken 50 times.
✗ Branch 77 not taken.
✓ Branch 79 taken 50 times.
✗ Branch 80 not taken.
✓ Branch 82 taken 50 times.
✗ Branch 83 not taken.
✓ Branch 85 taken 50 times.
✗ Branch 86 not taken.
✓ Branch 88 taken 50 times.
✗ Branch 89 not taken.
✓ Branch 91 taken 80 times.
✗ Branch 92 not taken.
✓ Branch 94 taken 80 times.
✗ Branch 95 not taken.
✓ Branch 97 taken 80 times.
✗ Branch 98 not taken.
✓ Branch 100 taken 80 times.
✗ Branch 101 not taken.
✓ Branch 103 taken 80 times.
✗ Branch 104 not taken.
✓ Branch 106 taken 80 times.
✗ Branch 107 not taken.
✓ Branch 109 taken 80 times.
✗ Branch 110 not taken.
✓ Branch 112 taken 80 times.
✗ Branch 113 not taken.
✓ Branch 115 taken 80 times.
✗ Branch 116 not taken.
✓ Branch 118 taken 80 times.
✗ Branch 119 not taken.
27116540 , numScv_(gridVolVars.problem().gridGeometry().numScv())
53 {}
54
55 //! operator for the access with an scv
56 template<class SubControlVolume, typename std::enable_if_t<!std::is_integral<SubControlVolume>::value, int> = 0>
57 710334810 const VolumeVariables& operator [](const SubControlVolume& scv) const
58 {
59
2/4
✓ Branch 0 taken 710334810 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 710334810 times.
✗ Branch 3 not taken.
1420669620 if (scv.dofIndex() < numScv_)
60 2131004430 return gridVolVars().volVars(scv.dofIndex());
61 else
62 return boundaryVolumeVariables_[getLocalIdx_(scv.dofIndex())];
63 }
64
65 //! operator for the access with an index
66 //! needed for Staggered methods for the access to the boundary volume variables
67 5516132696 const VolumeVariables& operator [](const std::size_t scvIdx) const
68 {
69
2/2
✓ Branch 0 taken 5387597557 times.
✓ Branch 1 taken 128535139 times.
5516132696 if (scvIdx < numScv_)
70 10775195114 return gridVolVars().volVars(scvIdx);
71 else
72 128535139 return boundaryVolumeVariables_[getLocalIdx_(scvIdx)];
73 }
74
75 /*!
76 * \brief bind the local view (r-value overload)
77 * This overload is called when an instance of this class is a temporary in the usage context
78 * This allows a usage like this: `const auto view = localView(...).bind(element);`
79 */
80 template<class FVElementGeometry, class SolutionVector>
81 StaggeredElementVolumeVariables bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
82 const FVElementGeometry& fvGeometry,
83 const SolutionVector& sol) &&
84 {
85 this->bind_(element, fvGeometry, sol);
86 return std::move(*this);
87 }
88
89 template<class FVElementGeometry, class SolutionVector>
90 void bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
91 const FVElementGeometry& fvGeometry,
92 const SolutionVector& sol) &
93
1/2
✓ Branch 1 taken 89392 times.
✗ Branch 2 not taken.
119616 { this->bind_(element, fvGeometry, sol); }
94
95 /*!
96 * \brief bind the local view (r-value overload)
97 * This overload is called when an instance of this class is a temporary in the usage context
98 * This allows a usage like this: `const auto view = localView(...).bind(element);`
99 */
100 template<class FVElementGeometry, class SolutionVector>
101 StaggeredElementVolumeVariables bindElement(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
102 const FVElementGeometry& fvGeometry,
103 const SolutionVector& sol) &&
104 {
105 this->bindElement_(element, fvGeometry, sol);
106 return std::move(*this);
107 }
108
109 template<class FVElementGeometry, class SolutionVector>
110 void bindElement(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
111 const FVElementGeometry& fvGeometry,
112 const SolutionVector& sol) &
113 { this->bindElement_(element, fvGeometry, sol); }
114
115 //! The global volume variables object we are a restriction of
116 const GridVolumeVariables& gridVolVars() const
117 { return *gridVolVarsPtr_; }
118
119 private:
120
121 //! Clear all local storage
122 void clear_()
123 {
124 1318046 boundaryVolVarIndices_.clear();
125
2/2
✓ Branch 0 taken 95973 times.
✓ Branch 1 taken 563050 times.
659023 boundaryVolumeVariables_.clear();
126 }
127
128 //! Binding of an element, prepares the volume variables within the element stencil
129 //! called by the local jacobian to prepare element assembly. Specialization callable with MultiTypeBlockVector.
130 template<class FVElementGeometry, class SolutionVector, typename std::enable_if_t<isMultiTypeBlockVector<SolutionVector>::value, int> = 0>
131 void bind_(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
132 const FVElementGeometry& fvGeometry,
133 const SolutionVector& sol)
134 {
135 // forward to the actual method
136
6/12
✓ Branch 1 taken 469136 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 469136 times.
✗ Branch 5 not taken.
✓ Branch 11 taken 15112 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 15112 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 15112 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 15112 times.
✗ Branch 21 not taken.
6420496 bind_(element, fvGeometry, sol[FVElementGeometry::GridGeometry::cellCenterIdx()]);
137 }
138
139 //! Binding of an element, prepares the volume variables within the element stencil
140 //! called by the local jacobian to prepare element assembly
141 template<class FVElementGeometry, class SolutionVector, typename std::enable_if_t<!isMultiTypeBlockVector<SolutionVector>::value, int> = 0>
142 3299640 void bind_(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
143 const FVElementGeometry& fvGeometry,
144 const SolutionVector& sol)
145 {
146
4/4
✓ Branch 0 taken 659023 times.
✓ Branch 1 taken 2640617 times.
✓ Branch 2 taken 659023 times.
✓ Branch 3 taken 2640617 times.
6599280 if (!fvGeometry.hasBoundaryScvf())
147 return;
148
149
2/2
✓ Branch 0 taken 95973 times.
✓ Branch 1 taken 563050 times.
659023 clear_();
150 1318046 boundaryVolVarIndices_.reserve(fvGeometry.numScvf());
151 1318046 boundaryVolumeVariables_.reserve(fvGeometry.numScvf());
152
153 // handle the boundary volume variables
154
6/6
✓ Branch 0 taken 2636092 times.
✓ Branch 1 taken 659023 times.
✓ Branch 2 taken 2636092 times.
✓ Branch 3 taken 659023 times.
✓ Branch 4 taken 1934724 times.
✓ Branch 5 taken 701368 times.
3954138 for (auto&& scvf : scvfs(fvGeometry))
155 {
156 // if we are not on a boundary, skip the rest
157
2/2
✓ Branch 0 taken 1934724 times.
✓ Branch 1 taken 701368 times.
2636092 if (!scvf.boundary())
158 1934724 continue;
159
160 701368 const auto& problem = gridVolVars().problem();
161 701368 auto boundaryPriVars = gridVolVars().getBoundaryPriVars(problem, sol, element, scvf);
162 1402736 const auto elemSol = elementSolution<FVElementGeometry>(std::move(boundaryPriVars));
163
2/4
✓ Branch 1 taken 134252 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 134252 times.
✗ Branch 5 not taken.
1402736 auto&& scvI = fvGeometry.scv(scvf.insideScvIdx());
164
165
1/2
✓ Branch 1 taken 134252 times.
✗ Branch 2 not taken.
701368 VolumeVariables volVars;
166
1/2
✓ Branch 1 taken 701368 times.
✗ Branch 2 not taken.
701368 volVars.update(elemSol,
167 problem,
168 element,
169 scvI);
170
171
1/2
✓ Branch 1 taken 701368 times.
✗ Branch 2 not taken.
701368 boundaryVolumeVariables_.emplace_back(std::move(volVars));
172
3/8
✓ Branch 1 taken 701368 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 701368 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 701368 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
1402736 boundaryVolVarIndices_.push_back(scvf.outsideScvIdx());
173 }
174 }
175
176 //! function to prepare the vol vars within the element
177 template<class FVElementGeometry, class SolutionVector>
178 void bindElement_(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
179 const FVElementGeometry& fvGeometry,
180 const SolutionVector& sol)
181 {}
182
183 const GridVolumeVariables* gridVolVarsPtr_;
184
185 //! map a global scv index to the local storage index
186 128535139 int getLocalIdx_(const int volVarIdx) const
187 {
188 385605417 auto it = std::find(boundaryVolVarIndices_.begin(), boundaryVolVarIndices_.end(), volVarIdx);
189
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 128535139 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 128535139 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 128535139 times.
385605417 assert(it != boundaryVolVarIndices_.end() && "Could not find the current volume variables for volVarIdx!");
190 385605417 return std::distance(boundaryVolVarIndices_.begin(), it);
191 }
192
193 std::vector<std::size_t> boundaryVolVarIndices_;
194 std::vector<VolumeVariables> boundaryVolumeVariables_;
195 const std::size_t numScv_;
196 };
197
198
199 /*!
200 * \ingroup StaggeredDiscretization
201 * \brief Class for the element volume variables vector for the staggered model.
202 Specialization in case the volume variables are not stored globally.
203 */
204 template<class GVV>
205 class StaggeredElementVolumeVariables<GVV, /*cachingEnabled*/false>
206 {
207 using PrimaryVariables = typename GVV::VolumeVariables::PrimaryVariables;
208
209 public:
210 //! export type of the grid volume variables
211 using GridVolumeVariables = GVV;
212
213 //! export type of the volume variables
214 using VolumeVariables = typename GridVolumeVariables::VolumeVariables;
215
216 //! Constructor
217 StaggeredElementVolumeVariables(const GridVolumeVariables& gridVolVars)
218 : gridVolVarsPtr_(&gridVolVars) {}
219
220 /*!
221 * \brief bind the local view (r-value overload)
222 * This overload is called when an instance of this class is a temporary in the usage context
223 * This allows a usage like this: `const auto view = localView(...).bind(element);`
224 */
225 template<class FVElementGeometry, class SolutionVector>
226 StaggeredElementVolumeVariables bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
227 const FVElementGeometry& fvGeometry,
228 const SolutionVector& sol) &&
229 {
230 this->bind_(element, fvGeometry, sol);
231 return std::move(*this);
232 }
233
234 template<class FVElementGeometry, class SolutionVector>
235 void bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
236 const FVElementGeometry& fvGeometry,
237 const SolutionVector& sol) &
238 { this->bind_(element, fvGeometry, sol); }
239
240 /*!
241 * \brief bind the local view (r-value overload)
242 * This overload is called when an instance of this class is a temporary in the usage context
243 * This allows a usage like this: `const auto view = localView(...).bind(element);`
244 */
245 template<class FVElementGeometry, class SolutionVector>
246 StaggeredElementVolumeVariables bindElement(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
247 const FVElementGeometry& fvGeometry,
248 const SolutionVector& sol) &&
249 {
250 this->bindElement_(element, fvGeometry, sol);
251 return std::move(*this);
252 }
253
254 template<class FVElementGeometry, class SolutionVector>
255 void bindElement(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
256 const FVElementGeometry& fvGeometry,
257 const SolutionVector& sol) &
258 { this->bindElement_(element, fvGeometry, sol); }
259
260 //! const operator for the access with an scv
261 template<class SubControlVolume, typename std::enable_if_t<!std::is_integral<SubControlVolume>::value, int> = 0>
262 const VolumeVariables& operator [](const SubControlVolume& scv) const
263 { return volumeVariables_[getLocalIdx_(scv.dofIndex())]; }
264
265 //! operator for the access with an scv
266 template<class SubControlVolume, typename std::enable_if_t<!std::is_integral<SubControlVolume>::value, int> = 0>
267 VolumeVariables& operator [](const SubControlVolume& scv)
268 { return volumeVariables_[getLocalIdx_(scv.dofIndex())]; }
269
270 //! const operator for the access with an index
271 const VolumeVariables& operator [](std::size_t scvIdx) const
272 { return volumeVariables_[getLocalIdx_(scvIdx)]; }
273
274 //! operator for the access with an index
275 VolumeVariables& operator [](std::size_t scvIdx)
276 { return volumeVariables_[getLocalIdx_(scvIdx)]; }
277
278 //! The global volume variables object we are a restriction of
279 const GridVolumeVariables& gridVolVars() const
280 { return *gridVolVarsPtr_; }
281
282 private:
283 //! Binding of an element, prepares the volume variables within the element stencil
284 //! called by the local jacobian to prepare element assembly. Specialization callable with MultiTypeBlockVector.
285 template<class FVElementGeometry, class SolutionVector, typename std::enable_if_t<isMultiTypeBlockVector<SolutionVector>::value, int> = 0>
286 void bind_(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
287 const FVElementGeometry& fvGeometry,
288 const SolutionVector& sol)
289 {
290 // forward to the actual method
291 bind_(element, fvGeometry, sol[FVElementGeometry::GridGeometry::cellCenterIdx()]);
292 }
293
294 //! Binding of an element, prepares the volume variables within the element stencil
295 //! called by the local jacobian to prepare element assembly
296 template<class FVElementGeometry, class SolutionVector, typename std::enable_if_t<!isMultiTypeBlockVector<SolutionVector>::value, int> = 0>
297 void bind_(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
298 const FVElementGeometry& fvGeometry,
299 const SolutionVector& sol)
300 {
301 clear_();
302
303 const auto& problem = gridVolVars().problem();
304 const auto& gridGeometry = fvGeometry.gridGeometry();
305 const auto globalI = gridGeometry.elementMapper().index(element);
306 const auto& map = gridGeometry.connectivityMap();
307 constexpr auto cellCenterIdx = FVElementGeometry::GridGeometry::cellCenterIdx();
308 const auto& connectivityMapI = map(cellCenterIdx, cellCenterIdx, globalI);
309 const auto numDofs = connectivityMapI.size();
310
311 auto&& scvI = fvGeometry.scv(globalI);
312
313 // resize local containers to the required size (for internal elements)
314 volumeVariables_.resize(numDofs+1);
315 volVarIndices_.resize(numDofs+1);
316 int localIdx = 0;
317
318 // Lambda to update the volume variables of the given index
319 auto doVolVarUpdate = [&](int globalJ)
320 {
321 const auto& elementJ = gridGeometry.element(globalJ);
322 auto&& scvJ = fvGeometry.scv(globalJ);
323 const auto elemSol = makeElementSolutionFromCellCenterPrivars<PrimaryVariables>(sol[globalJ]);
324 volumeVariables_[localIdx].update(elemSol,
325 problem,
326 elementJ,
327 scvJ);
328 volVarIndices_[localIdx] = scvJ.dofIndex();
329 ++localIdx;
330 };
331
332 // Update the volume variables of the element at hand
333 doVolVarUpdate(globalI);
334
335 // Update the volume variables of the neighboring elements
336 for (const auto& globalJ : connectivityMapI)
337 doVolVarUpdate(globalJ);
338
339 if (fvGeometry.hasBoundaryScvf())
340 {
341 // Update boundary volume variables
342 for (auto&& scvf : scvfs(fvGeometry))
343 {
344 // if we are not on a boundary, skip to the next scvf
345 if (!scvf.boundary())
346 continue;
347
348 volumeVariables_.resize(localIdx+1);
349 volVarIndices_.resize(localIdx+1);
350
351 auto boundaryPriVars = gridVolVars().getBoundaryPriVars(problem, sol, element, scvf);
352 auto elemSol = elementSolution<FVElementGeometry>(std::move(boundaryPriVars));
353 volumeVariables_[localIdx].update(elemSol,
354 problem,
355 element,
356 scvI);
357 volVarIndices_[localIdx] = scvf.outsideScvIdx();
358 ++localIdx;
359 }
360 }
361 }
362
363 //! Binding of an element, prepares the volume variables within the element stencil
364 //! called by the local jacobian to prepare element assembly. Specialization callable with MultiTypeBlockVector.
365 template<class FVElementGeometry, class SolutionVector, typename std::enable_if_t<isMultiTypeBlockVector<SolutionVector>::value, int> = 0>
366 void bindElement_(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
367 const FVElementGeometry& fvGeometry,
368 const SolutionVector& sol)
369 {
370 // forward to the actual method
371 bindElement_(element, fvGeometry, sol[FVElementGeometry::GridGeometry::cellCenterIdx()]);
372 }
373
374 //! Binding of an element, prepares only the volume variables of the element.
375 //! Specialization for Staggered models
376 template<class FVElementGeometry, class SolutionVector, typename std::enable_if_t<!isMultiTypeBlockVector<SolutionVector>::value, int> = 0>
377 void bindElement_(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
378 const FVElementGeometry& fvGeometry,
379 const SolutionVector& sol)
380 {
381 clear_();
382
383 const auto globalI = fvGeometry.gridGeometry().elementMapper().index(element);
384 volumeVariables_.resize(1);
385 volVarIndices_.resize(1);
386
387 // update the volume variables of the element
388 auto&& scv = fvGeometry.scv(globalI);
389
390 const auto elemSol = makeElementSolutionFromCellCenterPrivars<PrimaryVariables>(sol[globalI]);
391 volumeVariables_[0].update(elemSol,
392 gridVolVars().problem(),
393 element,
394 scv);
395 volVarIndices_[0] = scv.dofIndex();
396 }
397
398 //! Clear all local storage
399 void clear_()
400 {
401 volVarIndices_.clear();
402 volumeVariables_.clear();
403 }
404
405 const GridVolumeVariables* gridVolVarsPtr_;
406
407 int getLocalIdx_(const int volVarIdx) const
408 {
409 auto it = std::find(volVarIndices_.begin(), volVarIndices_.end(), volVarIdx);
410 assert(it != volVarIndices_.end() && "Could not find the current volume variables for volVarIdx!");
411 return std::distance(volVarIndices_.begin(), it);
412 }
413
414 std::vector<std::size_t> volVarIndices_;
415 std::vector<VolumeVariables> volumeVariables_;
416 };
417
418 } // end namespace Dumux
419
420 #endif
421