GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/discretization/staggered/fvelementgeometry.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 4 113 3.5%
Functions: 0 10 0.0%
Branches: 70 4885 1.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 StaggeredDiscretization
10 * \copydoc Dumux::StaggeredFVElementGeometry
11 */
12 #ifndef DUMUX_DISCRETIZATION_STAGGERED_FV_ELEMENT_GEOMETRY_HH
13 #define DUMUX_DISCRETIZATION_STAGGERED_FV_ELEMENT_GEOMETRY_HH
14
15 #include <optional>
16 #include <bitset>
17
18 #include <dumux/common/indextraits.hh>
19 #include <dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh>
20 #include <dumux/discretization/facecentered/staggered/normalaxis.hh>
21
22 namespace Dumux {
23
24 /*!
25 * \ingroup StaggeredDiscretization
26 * \brief Stencil-local finite volume geometry (scvs and scvfs) for staggered models
27 * This builds up the sub control volumes and sub control volume faces
28 * for each element in the local scope we are restricting to, e.g. stencil or element.
29 * \tparam GG the finite volume grid geometry type
30 * \tparam enableGridGeometryCache if the grid geometry is cached or not
31 * \note This class is specialized for versions with and without caching the fv geometries on the grid view
32 */
33 template<class GG, bool enableGridGeometryCache>
34 class StaggeredFVElementGeometry;
35
36 /*!
37 * \ingroup StaggeredDiscretization
38 * \brief Base class for the finite volume geometry vector for staggered models
39 * This locally builds up the sub control volumes and sub control volume faces
40 * for each element.
41 * Specialization for grid caching enabled
42 * \tparam GG the finite volume grid geometry type
43 */
44 template<class GG>
45
6/8
✓ Branch 0 taken 14152 times.
✓ Branch 1 taken 896032 times.
✓ Branch 3 taken 14152 times.
✓ Branch 4 taken 896032 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 412912 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 412912 times.
2739192 class StaggeredFVElementGeometry<GG, true> : public CCTpfaFVElementGeometry<GG, true>
46 {
47 using ParentType = CCTpfaFVElementGeometry<GG, true>;
48 using GridView = typename GG::GridView;
49 using GridIndexType = typename IndexTraits<GridView>::GridIndex;
50 using LocalIndexType = typename IndexTraits<GridView>::LocalIndex;
51 public:
52 //! export type of the element
53 using Element = typename GridView::template Codim<0>::Entity;
54 //! export type of subcontrol volume face
55 using SubControlVolumeFace = typename GG::SubControlVolumeFace;
56
57
22/38
✓ Branch 1 taken 1322 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1322 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 44 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 6 times.
✓ Branch 12 taken 44 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 6 times.
✓ Branch 15 taken 12 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✓ Branch 18 taken 12 times.
✓ Branch 19 taken 79 times.
✓ Branch 20 taken 1 times.
✓ Branch 21 taken 17 times.
✓ Branch 22 taken 79 times.
✓ Branch 23 taken 1 times.
✓ Branch 24 taken 17 times.
✓ Branch 25 taken 79 times.
✓ Branch 26 taken 1 times.
✗ Branch 27 not taken.
✓ Branch 28 taken 79 times.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✓ Branch 37 taken 63 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 63 times.
✗ Branch 41 not taken.
132028 using ParentType::ParentType;
58
59 //! Constructor getting a auxiliary cell center of face specific FvGridGeometry type.
60 //! Needed for the multi-domain framework.
61 template<class CellCenterOrFaceFVGridGeometry>
62 StaggeredFVElementGeometry(const CellCenterOrFaceFVGridGeometry& gridGeometry)
63
12/24
✓ Branch 1 taken 1355444 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1355444 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1355444 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1355444 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 50 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 50 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 50 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 50 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 50 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 50 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 17 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 17 times.
✗ Branch 35 not taken.
5422110 : ParentType(gridGeometry.actualGridGeometry()) {}
64
65 //! Get a sub control volume face with an element index and a local scvf index
66 using ParentType::scvf;
67 const SubControlVolumeFace& scvf(GridIndexType eIdx, LocalIndexType localScvfIdx) const
68 {
69
30/39
✓ Branch 0 taken 2400 times.
✓ Branch 1 taken 600 times.
✓ Branch 2 taken 2400 times.
✓ Branch 3 taken 600 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✓ Branch 23 taken 53878 times.
✓ Branch 24 taken 346528 times.
✓ Branch 25 taken 278986698 times.
✓ Branch 26 taken 400406 times.
✓ Branch 27 taken 278986698 times.
✓ Branch 30 taken 48 times.
✓ Branch 31 taken 284208837 times.
✓ Branch 32 taken 48 times.
✓ Branch 33 taken 288104029 times.
✗ Branch 34 not taken.
✓ Branch 35 taken 3895192 times.
✓ Branch 36 taken 10782296 times.
✓ Branch 37 taken 202085552 times.
✓ Branch 38 taken 14485338 times.
✓ Branch 39 taken 274943710 times.
✓ Branch 40 taken 4358938 times.
✓ Branch 41 taken 72867958 times.
✓ Branch 42 taken 668404 times.
✓ Branch 43 taken 10732 times.
✓ Branch 44 taken 180432 times.
✓ Branch 45 taken 5156616 times.
✓ Branch 46 taken 202632 times.
✓ Branch 47 taken 5637904 times.
✓ Branch 48 taken 69908 times.
✓ Branch 49 taken 483020 times.
✓ Branch 50 taken 35200 times.
✓ Branch 51 taken 800 times.
2857529270 return this->gridGeometry().scvf(eIdx, localScvfIdx);
70 }
71 };
72
73 /*!
74 * \ingroup StaggeredDiscretization
75 * \brief Base class for the finite volume geometry vector for staggered models
76 * This locally builds up the sub control volumes and sub control volume faces
77 * for each element.
78 * Specialization for grid caching disabled
79 * \tparam GG the finite volume grid geometry type
80 */
81 template<class GG>
82 class StaggeredFVElementGeometry<GG, false>
83 {
84 using ThisType = StaggeredFVElementGeometry<GG, false>;
85 using GridView = typename GG::GridView;
86 using GridIndexType = typename IndexTraits<GridView>::GridIndex;
87 using LocalIndexType = typename IndexTraits<GridView>::LocalIndex;
88 public:
89 //! export type of the element
90 using Element = typename GridView::template Codim<0>::Entity;
91 //! export type of subcontrol volume
92 using SubControlVolume = typename GG::SubControlVolume;
93 //! export type of subcontrol volume face
94 using SubControlVolumeFace = typename GG::SubControlVolumeFace;
95 //! export type of finite volume grid geometry
96 using GridGeometry = GG;
97
98 //! Constructor getting a auxiliary cell center of face specific FvGridGeometry type.
99 //! Needed for the multi-domain framework.
100 template<class CellCenterOrFaceFVGridGeometry>
101 StaggeredFVElementGeometry(const CellCenterOrFaceFVGridGeometry& gridGeometry)
102 : gridGeometryPtr_(&gridGeometry.actualGridGeometry()) {}
103
104 //! Constructor
105 StaggeredFVElementGeometry(const GridGeometry& gridGeometry)
106 : gridGeometryPtr_(&gridGeometry) {}
107
108 //! Get a sub control volume face with an element index and a local scvf index
109 const SubControlVolumeFace& scvf(GridIndexType eIdx, LocalIndexType localScvfIdx) const
110 {
111 return scvf(this->gridGeometry().localToGlobalScvfIndex(eIdx, localScvfIdx));
112 }
113
114 //! Get an element sub control volume with a global scv index
115 //! We separate element and neighbor scvs to speed up mapping
116 const SubControlVolume& scv(GridIndexType scvIdx) const
117 {
118 if (scvIdx == scvIndices_[0])
119 return scvs_[0];
120 else
121 return neighborScvs_[findLocalIndex_(scvIdx, neighborScvIndices_)];
122 }
123
124 //! Get an element sub control volume face with a global scvf index
125 //! We separate element and neighbor scvfs to speed up mapping
126 const SubControlVolumeFace& scvf(GridIndexType scvfIdx) const
127 {
128 auto it = std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx);
129 if (it != scvfIndices_.end())
130 return scvfs_[std::distance(scvfIndices_.begin(), it)];
131 else
132 return neighborScvfs_[findLocalIndex_(scvfIdx, neighborScvfIndices_)];
133 }
134
135 //! iterator range for sub control volumes. Iterates over
136 //! all scvs of the bound element (not including neighbor scvs)
137 //! This is a free function found by means of ADL
138 //! To iterate over all sub control volumes of this FVElementGeometry use
139 //! for (auto&& scv : scvs(fvGeometry))
140 friend inline Dune::IteratorRange<typename std::array<SubControlVolume, 1>::const_iterator>
141 scvs(const ThisType& g)
142 {
143 using IteratorType = typename std::array<SubControlVolume, 1>::const_iterator;
144 return Dune::IteratorRange<IteratorType>(g.scvs_.begin(), g.scvs_.end());
145 }
146
147 //! iterator range for sub control volumes faces. Iterates over
148 //! all scvfs of the bound element (not including neighbor scvfs)
149 //! This is a free function found by means of ADL
150 //! To iterate over all sub control volume faces of this FVElementGeometry use
151 //! for (auto&& scvf : scvfs(fvGeometry))
152 friend inline Dune::IteratorRange<typename std::vector<SubControlVolumeFace>::const_iterator>
153 scvfs(const ThisType& g)
154 {
155 using IteratorType = typename std::vector<SubControlVolumeFace>::const_iterator;
156 return Dune::IteratorRange<IteratorType>(g.scvfs_.begin(), g.scvfs_.end());
157 }
158
159 //! number of sub control volumes in this fv element geometry
160 std::size_t numScv() const
161 { return scvs_.size(); }
162
163 //! number of sub control volumes in this fv element geometry
164 std::size_t numScvf() const
165 { return scvfs_.size(); }
166
167 /*!
168 * \brief bind the local view (r-value overload)
169 * This overload is called when an instance of this class is a temporary in the usage context
170 * This allows a usage like this: `const auto view = localView(...).bind(element);`
171 */
172 StaggeredFVElementGeometry bind(const Element& element) &&
173 {
174 this->bind_(element);
175 return std::move(*this);
176 }
177
178 //! bind this local view to a specific element (full stencil)
179 void bind(const Element& element) &
180 { this->bind_(element); }
181
182 /*!
183 * \brief bind the local view (r-value overload)
184 * This overload is called when an instance of this class is a temporary in the usage context
185 * This allows a usage like this: `const auto view = localView(...).bindElement(element);`
186 */
187 StaggeredFVElementGeometry bindElement(const Element& element) &&
188 {
189 this->bindElement_(element);
190 return std::move(*this);
191 }
192
193 //! bind this local view to a specific element
194 void bindElement(const Element& element) &
195 { this->bindElement_(element); }
196
197 //! Returns true if bind/bindElement has already been called
198 bool isBound() const
199 { return static_cast<bool>(element_); }
200
201 //! The bound element
202 const Element& element() const
203 { return *element_; }
204
205 //! The grid finite volume geometry we are a restriction of
206 const GridGeometry& gridGeometry() const
207 { return *gridGeometryPtr_; }
208
209 //! Returns whether one of the geometry's scvfs lies on a boundary
210 bool hasBoundaryScvf() const
211 { return hasBoundaryScvf_; }
212
213 //! Create the geometry of a given sub control volume
214 typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const
215 { return gridGeometryPtr_->element(scv.dofIndex()).geometry(); }
216
217 //! Create the geometry of a given sub control volume face
218 typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const
219 {
220 const auto insideElementIndex = scvf.insideScvIdx();
221 const auto elementGeometry = (insideElementIndex != scvIndices_[0]) ?
222 element_->geometry() :
223 gridGeometryPtr_->element(insideElementIndex).geometry();
224 const auto center = elementGeometry.center();
225 const auto normalAxis = Dumux::normalAxis(scvf.unitOuterNormal());
226
227 auto lowerLeft = elementGeometry.corner(0);
228 auto upperRight = elementGeometry.corner(elementGeometry.corners()-1);
229
230 // shift corners to scvf plane and halve lateral faces
231 lowerLeft[normalAxis] = center[normalAxis];
232 upperRight[normalAxis] = center[normalAxis];
233
234 auto inPlaneAxes = std::move(std::bitset<SubControlVolumeFace::Traits::dimWorld>{}.set());
235 inPlaneAxes.set(normalAxis, false);
236
237 return {lowerLeft, upperRight, inPlaneAxes};
238 }
239
240 private:
241 //! Binding of an element preparing the geometries only inside the element
242 void bindElement_(const Element& element)
243 {
244 clear_();
245 element_ = element;
246 scvfs_.reserve(element.subEntities(1));
247 scvfIndices_.reserve(element.subEntities(1));
248 makeElementGeometries_();
249 }
250
251 //! Binding of an element preparing the geometries of the whole stencil
252 //! called by the local jacobian to prepare element assembly
253 void bind_(const Element& element)
254 {
255 bindElement_(element);
256
257 neighborScvs_.reserve(element.subEntities(1));
258 neighborScvfIndices_.reserve(element.subEntities(1));
259 neighborScvfs_.reserve(element.subEntities(1));
260
261 std::vector<GridIndexType> handledNeighbors;
262 handledNeighbors.reserve(element.subEntities(1));
263 for (const auto& intersection : intersections(gridGeometry().gridView(), element))
264 {
265 if (intersection.neighbor())
266 {
267 const auto outside = intersection.outside();
268 const auto outsideIdx = gridGeometry().elementMapper().index(outside);
269
270 // make outside geometries only if not done yet (could happen on non-conforming grids)
271 if ( std::find(handledNeighbors.begin(), handledNeighbors.end(), outsideIdx) == handledNeighbors.end() )
272 {
273 makeNeighborGeometries_(outside, outsideIdx);
274 handledNeighbors.push_back(outsideIdx);
275 }
276 }
277 }
278 }
279
280 //! create scvs and scvfs of the bound element
281 void makeElementGeometries_()
282 {
283 const auto& element = *element_;
284 const auto eIdx = gridGeometry().elementMapper().index(element);
285 scvs_[0] = SubControlVolume(element.geometry(), eIdx);
286 scvIndices_[0] = eIdx;
287
288 const auto& scvFaceIndices = gridGeometry().scvfIndicesOfScv(eIdx);
289 const auto& neighborVolVarIndices = gridGeometry().neighborVolVarIndices(eIdx);
290
291 typename GridGeometry::GeometryHelper geometryHelper(element, gridGeometry().gridView());
292
293 int scvfCounter = 0;
294 for (const auto& intersection : intersections(gridGeometry().gridView(), element))
295 {
296 const auto& scvfNeighborVolVarIndex = neighborVolVarIndices[scvfCounter];
297
298 if (intersection.neighbor() || intersection.boundary())
299 {
300 geometryHelper.updateLocalFace(gridGeometry().intersectionMapper(), intersection);
301 std::vector<GridIndexType> scvIndices{eIdx, scvfNeighborVolVarIndex};
302 scvfs_.emplace_back(intersection,
303 intersection.geometry(),
304 scvFaceIndices[scvfCounter],
305 scvIndices,
306 geometryHelper);
307 scvfIndices_.emplace_back(scvFaceIndices[scvfCounter]);
308 scvfCounter++;
309
310 if (intersection.boundary())
311 hasBoundaryScvf_ = true;
312 }
313 }
314 }
315
316 //! create the necessary scvs and scvfs of the neighbor elements to the bound elements
317 void makeNeighborGeometries_(const Element& element, const GridIndexType eIdx)
318 {
319 // using ScvfGridIndexStorage = typename SubControlVolumeFace::Traits::GridIndexStorage;
320
321 // create the neighbor scv
322 neighborScvs_.emplace_back(element.geometry(), eIdx);
323 neighborScvIndices_.push_back(eIdx);
324
325 const auto& scvFaceIndices = gridGeometry().scvfIndicesOfScv(eIdx);
326 const auto& neighborVolVarIndices = gridGeometry().neighborVolVarIndices(eIdx);
327
328 typename GridGeometry::GeometryHelper geometryHelper(element, gridGeometry().gridView());
329
330 int scvfCounter = 0;
331 for (const auto& intersection : intersections(gridGeometry().gridView(), element))
332 {
333 const auto& scvfNeighborVolVarIndex = neighborVolVarIndices[scvfCounter];
334 geometryHelper.updateLocalFace(gridGeometry().intersectionMapper(), intersection);
335
336 if (intersection.neighbor())
337 {
338 // only create subcontrol faces where the outside element is the bound element
339 if (intersection.outside() == *element_)
340 {
341 std::vector<GridIndexType> scvIndices{eIdx, scvfNeighborVolVarIndex};
342 neighborScvfs_.emplace_back(intersection,
343 intersection.geometry(),
344 scvFaceIndices[scvfCounter],
345 scvIndices,
346 geometryHelper);
347
348 neighborScvfIndices_.push_back(scvFaceIndices[scvfCounter]);
349 }
350 scvfCounter++;
351 }
352 else if (intersection.boundary())
353 scvfCounter++;
354 }
355 }
356
357 const LocalIndexType findLocalIndex_(const GridIndexType idx,
358 const std::vector<GridIndexType>& indices) const
359 {
360 auto it = std::find(indices.begin(), indices.end(), idx);
361 assert(it != indices.end() && "Could not find the scv/scvf! Make sure to properly bind this class!");
362 return std::distance(indices.begin(), it);
363 }
364
365 //! Clear all local data
366 void clear_()
367 {
368 scvfIndices_.clear();
369 scvfs_.clear();
370
371 neighborScvIndices_.clear();
372 neighborScvfIndices_.clear();
373 neighborScvs_.clear();
374 neighborScvfs_.clear();
375
376 hasBoundaryScvf_ = false;
377 }
378
379 std::optional<Element> element_; //!< the element to which this fvgeometry is bound
380 const GridGeometry* gridGeometryPtr_; //!< the grid fvgeometry
381
382 // local storage after binding an element
383 std::array<GridIndexType, 1> scvIndices_;
384 std::array<SubControlVolume, 1> scvs_;
385
386 std::vector<GridIndexType> scvfIndices_;
387 std::vector<SubControlVolumeFace> scvfs_;
388
389 std::vector<GridIndexType> neighborScvIndices_;
390 std::vector<SubControlVolume> neighborScvs_;
391
392 std::vector<GridIndexType> neighborScvfIndices_;
393 std::vector<SubControlVolumeFace> neighborScvfs_;
394
395 bool hasBoundaryScvf_ = false;
396 };
397
398
399 } // end namespace
400
401 #endif
402