GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/discretization/porenetwork/fvelementgeometry.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 48 48 100.0%
Functions: 11 11 100.0%
Branches: 77 128 60.2%

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 * \file
9 * \ingroup PoreNetworkDiscretization
10 * \brief Base class for the local geometry for porenetworks
11 */
12 #ifndef DUMUX_DISCRETIZATION_PNM_FV_ELEMENT_GEOMETRY_HH
13 #define DUMUX_DISCRETIZATION_PNM_FV_ELEMENT_GEOMETRY_HH
14
15 #include <optional>
16 #include <array>
17 #include <utility>
18
19 #include <dune/geometry/type.hh>
20
21 #include <dumux/common/indextraits.hh>
22 #include <dumux/discretization/scvandscvfiterators.hh>
23
24 namespace Dumux::PoreNetwork {
25
26 /*!
27 * \ingroup PoreNetworkDiscretization
28 * \brief Base class for the local geometry for porenetworks
29 * \tparam GG the finite volume grid geometry type
30 * \tparam enableFVGridGeometryCache if the grid geometry is cached or not
31 */
32 template<class GG, bool enableFVGridGeometryCache>
33 class PNMFVElementGeometry;
34
35 //! specialization in case the FVElementGeometries are stored
36 template<class GG>
37 class PNMFVElementGeometry<GG, true>
38 {
39 using GridView = typename GG::GridView;
40 static constexpr int dim = GridView::dimension;
41 static constexpr int dimWorld = GridView::dimensionworld;
42 using GridIndexType = typename IndexTraits<GridView>::GridIndex;
43 using LocalIndexType = typename IndexTraits<GridView>::LocalIndex;
44 using CoordScalar = typename GridView::ctype;
45 using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType;
46 public:
47 //! export type of subcontrol volume
48 using SubControlVolume = typename GG::SubControlVolume;
49 //! export type of subcontrol volume face
50 using SubControlVolumeFace = typename GG::SubControlVolumeFace;
51 //! export type of finite volume grid geometry
52 using GridGeometry = GG;
53 //! export element type
54 using Element = typename GridView::template Codim<0>::Entity;
55 //! the maximum number of scvs per element
56 static constexpr std::size_t maxNumElementScvs = 2;
57
58 //! Constructor
59 PNMFVElementGeometry(const GridGeometry& gridGeometry)
60 : gridGeometryPtr_(&gridGeometry) {}
61
62 //! Get a sub control volume with a local scv index
63 const SubControlVolume& scv(LocalIndexType scvIdx) const
64 {
65 return gridGeometry().scvs(eIdx_)[scvIdx];
66 }
67
68 //! Get a sub control volume face with a local scvf index
69 const SubControlVolumeFace& scvf(LocalIndexType scvfIdx) const
70 {
71 return gridGeometry().scvfs(eIdx_)[scvfIdx];
72 }
73
74 //! iterator range for sub control volumes. Iterates over
75 //! all scvs of the bound element.
76 //! This is a free function found by means of ADL
77 //! To iterate over all sub control volumes of this FVElementGeometry use
78 //! for (auto&& scv : scvs(fvGeometry))
79 friend inline auto scvs(const PNMFVElementGeometry& fvGeometry)
80 {
81 const auto& g = fvGeometry.gridGeometry();
82 using Iter = typename std::decay_t<decltype(g.scvs(fvGeometry.eIdx_))>::const_iterator;
83 return Dune::IteratorRange<Iter>(g.scvs(fvGeometry.eIdx_).begin(), g.scvs(fvGeometry.eIdx_).end());
84 }
85
86 //! iterator range for sub control volumes faces. Iterates over
87 //! all scvfs of the bound element.
88 //! This is a free function found by means of ADL
89 //! To iterate over all sub control volume faces of this FVElementGeometry use
90 //! for (auto&& scvf : scvfs(fvGeometry))
91 friend inline auto scvfs(const PNMFVElementGeometry& fvGeometry)
92 {
93 const auto& g = fvGeometry.gridGeometry();
94 using Iter = typename std::decay_t<decltype(g.scvfs(fvGeometry.eIdx_))>::const_iterator;
95 return Dune::IteratorRange<Iter>(g.scvfs(fvGeometry.eIdx_).begin(), g.scvfs(fvGeometry.eIdx_).end());
96 }
97
98 //! Get a local finite element basis
99 const FeLocalBasis& feLocalBasis() const
100 {
101 return gridGeometry().feCache().get(element_->type()).localBasis();
102 }
103
104 //! The total number of sub control volumes
105 std::size_t numScv() const
106 {
107 return gridGeometry().scvs(eIdx_).size();
108 }
109
110 //! The total number of sub control volume faces
111 std::size_t numScvf() const
112 {
113 return gridGeometry().scvfs(eIdx_).size();
114 }
115
116 /*!
117 * \brief bind the local view (r-value overload)
118 * This overload is called when an instance of this class is a temporary in the usage context
119 * This allows a usage like this: `const auto view = localView(...).bind(element);`
120 */
121 PNMFVElementGeometry bind(const Element& element) &&
122 {
123 this->bindElement(element);
124 return std::move(*this);
125 }
126
127 //! this function is for compatibility reasons with cc methods
128 //! The box stencil is always element-local so bind and bindElement
129 //! are identical.
130 void bind(const Element& element) &
131 { this->bindElement(element); }
132
133 /*!
134 * \brief bind the local view (r-value overload)
135 * This overload is called when an instance of this class is a temporary in the usage context
136 * This allows a usage like this: `const auto view = localView(...).bindElement(element);`
137 */
138 PNMFVElementGeometry bindElement(const Element& element) &&
139 {
140 this->bindElement(element);
141 return std::move(*this);
142 }
143
144 //! Binding of an element, has to be called before using the fvgeometries
145 //! Prepares all the volume variables within the element
146 //! For compatibility reasons with the FVGeometry cache being disabled
147 void bindElement(const Element& element) &
148 {
149 element_ = element;
150 eIdx_ = gridGeometry().elementMapper().index(element);
151 }
152
153 //! Returns true if bind/bindElement has already been called
154 bool isBound() const
155 { return static_cast<bool>(element_); }
156
157 //! The bound element
158 const Element& element() const
159 { return *element_; }
160
161 //! The global finite volume geometry we are a restriction of
162 const GridGeometry& gridGeometry() const
163 { return *gridGeometryPtr_; }
164
165 //! Returns whether one of the geometry's scvfs lies on a boundary
166 bool hasBoundaryScvf() const
167 { return gridGeometry().hasBoundaryScvf(eIdx_); }
168
169 //! Create the geometry of a given sub control volume
170 typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const
171 {
172 const auto geo = element().geometry();
173 return {
174 Dune::GeometryTypes::simplex(dim),
175 std::array{{ geo.corner(scv.indexInElement()), geo.center() }}
176 };
177 }
178
179 //! Create the geometry of a given sub control volume face
180 typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const
181 { return { scvf.center() }; }
182
183 private:
184 std::optional<Element> element_;
185 const GridGeometry* gridGeometryPtr_;
186
187 GridIndexType eIdx_;
188 };
189
190 //! specialization in case the FVElementGeometries are not stored
191 template<class GG>
192
6/9
✓ Branch 0 taken 892 times.
✓ Branch 1 taken 1011839 times.
✓ Branch 3 taken 892 times.
✓ Branch 4 taken 316416 times.
✓ Branch 6 taken 189 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1754916 times.
✗ Branch 9 not taken.
✗ Branch 5 not taken.
3085144 class PNMFVElementGeometry<GG, false>
193 {
194 using GridView = typename GG::GridView;
195 static constexpr int dim = GridView::dimension;
196 static constexpr int dimWorld = GridView::dimensionworld;
197 using GridIndexType = typename IndexTraits<GridView>::GridIndex;
198 using LocalIndexType = typename IndexTraits<GridView>::LocalIndex;
199 using CoordScalar = typename GridView::ctype;
200 using FeLocalBasis = typename GG::FeCache::FiniteElementType::Traits::LocalBasisType;
201
202 public:
203 //! export type of subcontrol volume
204 using SubControlVolume = typename GG::SubControlVolume;
205 //! export type of subcontrol volume face
206 using SubControlVolumeFace = typename GG::SubControlVolumeFace;
207 //! export type of finite volume grid geometry
208 using GridGeometry = GG;
209 //! export element type
210 using Element = typename GridView::template Codim<0>::Entity;
211 //! the maximum number of scvs per element
212 static constexpr std::size_t maxNumElementScvs = 2;
213
214 //! Constructor
215 2807832 PNMFVElementGeometry(const GridGeometry& gridGeometry)
216
15/33
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✓ Branch 24 taken 1 times.
✗ Branch 25 not taken.
✓ Branch 29 taken 338742 times.
✗ Branch 30 not taken.
✓ Branch 32 taken 338742 times.
✗ Branch 33 not taken.
✓ Branch 35 taken 182830 times.
✗ Branch 36 not taken.
✓ Branch 38 taken 182830 times.
✗ Branch 39 not taken.
✓ Branch 43 taken 2 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 2 times.
✗ Branch 47 not taken.
✓ Branch 49 taken 1 times.
✗ Branch 50 not taken.
✓ Branch 52 taken 1 times.
✗ Branch 53 not taken.
✓ Branch 6 taken 952 times.
5615664 : gridGeometryPtr_(&gridGeometry) {}
217
218 //! Get a sub control volume with a local scv index
219 47295213 const SubControlVolume& scv(LocalIndexType scvIdx) const
220 {
221
12/14
✓ Branch 0 taken 2522980 times.
✓ Branch 1 taken 2123811 times.
✓ Branch 2 taken 1092264 times.
✓ Branch 3 taken 2443 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 2522974 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1092264 times.
✓ Branch 8 taken 22032 times.
✓ Branch 9 taken 16808 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 32058562 times.
✓ Branch 12 taken 18206 times.
✓ Branch 13 taken 25653 times.
43561546 return scvs_[scvIdx];
222 }
223
224 //! Get a sub control volume face with a local scvf index
225 32064942 const SubControlVolumeFace& scvf(LocalIndexType scvfIdx) const
226 {
227
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 31950154 times.
✓ Branch 2 taken 93628 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 21160 times.
32064942 return scvfs_[0];
228 }
229
230 //! iterator range for sub control volumes. Iterates over
231 //! all scvs of the bound element.
232 //! This is a free function found by means of ADL
233 //! To iterate over all sub control volumes of this FVElementGeometry use
234 //! for (auto&& scv : scvs(fvGeometry))
235 122490180 friend inline auto scvs(const PNMFVElementGeometry& fvGeometry)
236 {
237 using Iter = typename std::decay_t<decltype(fvGeometry.scvs_)>::const_iterator;
238
2/4
✓ Branch 0 taken 2522976 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1092264 times.
✗ Branch 3 not taken.
122490180 return Dune::IteratorRange<Iter>(fvGeometry.scvs_.begin(), fvGeometry.scvs_.end());
239 }
240
241 //! iterator range for sub control volumes faces. Iterates over
242 //! all scvfs of the bound element.
243 //! This is a free function found by means of ADL
244 //! To iterate over all sub control volume faces of this FVElementGeometry use
245 //! for (auto&& scvf : scvfs(fvGeometry))
246 45501314 friend inline auto scvfs(const PNMFVElementGeometry& fvGeometry)
247 {
248 using Iter = typename std::decay_t<decltype(fvGeometry.scvfs_)>::const_iterator;
249
1/2
✓ Branch 0 taken 38840 times.
✗ Branch 1 not taken.
45501314 return Dune::IteratorRange<Iter>(fvGeometry.scvfs_.begin(), fvGeometry.scvfs_.end());
250 }
251
252 //! Get a local finite element basis
253 const FeLocalBasis& feLocalBasis() const
254 {
255 return gridGeometry().feCache().get(element_->type()).localBasis();
256 }
257
258 //! The total number of sub control volumes
259 6169863 std::size_t numScv() const
260 {
261
2/4
✓ Branch 0 taken 2522976 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1092264 times.
✗ Branch 3 not taken.
19567574 return scvs_.size();
262 }
263
264 //! The total number of sub control volume faces
265 std::size_t numScvf() const
266 {
267 5421023 return scvfs_.size();
268 }
269
270 /*!
271 * \brief bind the local view (r-value overload)
272 * This overload is called when an instance of this class is a temporary in the usage context
273 * This allows a usage like this: `const auto view = localView(...).bind(element);`
274 */
275 318007 PNMFVElementGeometry bind(const Element& element) &&
276 {
277
1/2
✓ Branch 1 taken 318007 times.
✗ Branch 2 not taken.
318007 this->bindElement(element);
278 318007 return std::move(*this);
279 }
280
281 //! this function is for compatibility reasons with cc methods
282 //! The box stencil is always element-local so bind and bindElement
283 //! are identical.
284 1424371 void bind(const Element& element) &
285
7/22
✓ Branch 1 taken 36822 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 51716 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 15 taken 1007 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 2184 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 2184 times.
✗ Branch 22 not taken.
✓ Branch 24 taken 1007 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 1007 times.
✗ Branch 28 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
1424371 { this->bindElement(element); }
286
287 /*!
288 * \brief bind the local view (r-value overload)
289 * This overload is called when an instance of this class is a temporary in the usage context
290 * This allows a usage like this: `const auto view = localView(...).bindElement(element);`
291 */
292 1330 PNMFVElementGeometry bindElement(const Element& element) &&
293 {
294
1/2
✓ Branch 1 taken 1330 times.
✗ Branch 2 not taken.
1330 this->bindElement(element);
295 1330 return std::move(*this);
296 }
297
298 //! Binding of an element, has to be called before using the fvgeometries
299 //! Prepares all the volume variables within the element
300 //! For compatibility reasons with the FVGeometry cache being disabled
301 3955347 void bindElement(const Element& element) &
302 {
303
2/2
✓ Branch 0 taken 1166603 times.
✓ Branch 1 taken 2788744 times.
3955347 element_ = element;
304 7910694 eIdx_ = gridGeometry().elementMapper().index(element);
305 3955347 makeElementGeometries(element);
306 3955347 }
307
308 //! Returns true if bind/bindElement has already been called
309 bool isBound() const
310 { return static_cast<bool>(element_); }
311
312 //! The bound element
313 3343377 const Element& element() const
314
2/4
✓ Branch 8 taken 88498 times.
✗ Branch 9 not taken.
✓ Branch 5 taken 5130 times.
✗ Branch 6 not taken.
3346532 { return *element_; }
315
316 //! The global finite volume geometry we are a restriction of
317 109306837 const GridGeometry& gridGeometry() const
318
13/14
✓ Branch 4 taken 155312 times.
✓ Branch 5 taken 1722274 times.
✓ Branch 6 taken 76020 times.
✓ Branch 7 taken 847020 times.
✓ Branch 9 taken 13988 times.
✓ Branch 10 taken 6790 times.
✓ Branch 8 taken 16990 times.
✓ Branch 11 taken 7308 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 35427936 times.
✓ Branch 2 taken 350280 times.
✓ Branch 3 taken 4224 times.
✓ Branch 1 taken 65462 times.
✓ Branch 0 taken 29440 times.
54089414 { return *gridGeometryPtr_; }
319
320 //! Returns whether one of the geometry's scvfs lies on a boundary
321 4622 bool hasBoundaryScvf() const
322
4/4
✓ Branch 0 taken 2206 times.
✓ Branch 1 taken 1464 times.
✓ Branch 2 taken 476 times.
✓ Branch 3 taken 476 times.
4622 { return hasBoundaryScvf_; }
323
324 //! Create the geometry of a given sub control volume
325 typename SubControlVolume::Traits::Geometry geometry(const SubControlVolume& scv) const
326 {
327 const auto geo = element().geometry();
328 return {
329 Dune::GeometryTypes::simplex(dim),
330 std::array{{ geo.corner(scv.indexInElement()), geo.center() }}
331 };
332 }
333
334 //! Create the geometry of a given sub control volume face
335 typename SubControlVolumeFace::Traits::Geometry geometry(const SubControlVolumeFace& scvf) const
336 { return { scvf.center() }; }
337
338 private:
339
340 3955347 void makeElementGeometries(const Element& element)
341 {
342 3955347 hasBoundaryScvf_ = false;
343
344 // get the element geometry
345 11866041 auto elementGeometry = element.geometry();
346
347 // construct the sub control volumes
348
2/2
✓ Branch 1 taken 7910694 times.
✓ Branch 2 taken 3955347 times.
11866041 for (LocalIndexType scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx)
349 {
350 // get associated dof index
351 7910694 const auto dofIdxGlobal = gridGeometry().vertexMapper().subIndex(element, scvLocalIdx, dim);
352
353 // get the corners
354 15816828 auto corners = std::array{elementGeometry.corner(scvLocalIdx), elementGeometry.center()};
355
356 // get the fractional volume associated with the scv
357 7910694 const auto volume = gridGeometry().poreVolume(dofIdxGlobal) / gridGeometry().coordinationNumber(dofIdxGlobal);
358
359 // add scv to the local container
360
2/2
✓ Branch 1 taken 2104898 times.
✓ Branch 2 taken 5805796 times.
7910694 scvs_[scvLocalIdx] = SubControlVolume(dofIdxGlobal,
361 scvLocalIdx,
362 eIdx_,
363 std::move(corners),
364 volume);
365
366
2/2
✓ Branch 0 taken 2104898 times.
✓ Branch 1 taken 5805796 times.
7910694 if (gridGeometry().poreLabel(dofIdxGlobal) > 0)
367 2104898 hasBoundaryScvf_ = true;
368 }
369
370 // construct the inner sub control volume face
371 11863761 auto unitOuterNormal = elementGeometry.corner(1) - elementGeometry.corner(0);
372 3955347 unitOuterNormal /= unitOuterNormal.two_norm();
373 LocalIndexType scvfLocalIdx = 0;
374 7910694 scvfs_[0] = SubControlVolumeFace(elementGeometry.center(),
375 std::move(unitOuterNormal),
376 3955347 gridGeometry().throatCrossSectionalArea(gridGeometry().elementMapper().index(element)),
377 scvfLocalIdx++,
378 std::array<LocalIndexType, 2>({0, 1}));
379 3955347 }
380
381 //! The bound element
382 std::optional<Element> element_;
383 GridIndexType eIdx_;
384
385 //! The global geometry this is a restriction of
386 const GridGeometry* gridGeometryPtr_;
387
388 //! vectors to store the geometries locally after binding an element
389 std::array<SubControlVolume, 2> scvs_;
390 std::array<SubControlVolumeFace, 1> scvfs_;
391
392 bool hasBoundaryScvf_ = false;
393 };
394
395 } // end namespace Dumux::PoreNetwork
396
397 #endif
398