GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/discretization/pq1bubble/geometryhelper.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 93 114 81.6%
Functions: 49 117 41.9%
Branches: 67 147 45.6%

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 PQ1BubbleDiscretization
10 * \brief Helper class constructing the dual grid finite volume geometries
11 * for the cvfe discretizazion method
12 */
13 #ifndef DUMUX_DISCRETIZATION_PQ1BUBBLE_GEOMETRY_HELPER_HH
14 #define DUMUX_DISCRETIZATION_PQ1BUBBLE_GEOMETRY_HELPER_HH
15
16 #include <array>
17
18 #include <dune/common/exceptions.hh>
19
20 #include <dune/geometry/type.hh>
21 #include <dune/geometry/referenceelements.hh>
22 #include <dune/geometry/multilineargeometry.hh>
23 #include <dune/common/reservedvector.hh>
24
25 #include <dumux/common/math.hh>
26 #include <dumux/geometry/volume.hh>
27 #include <dumux/discretization/box/boxgeometryhelper.hh>
28
29 namespace Dumux {
30
31 //! Traits for an efficient corner storage for the PQ1Bubble method
32 template <class ct>
33 struct PQ1BubbleMLGeometryTraits : public Dune::MultiLinearGeometryTraits<ct>
34 {
35 // we use static vectors to store the corners as we know
36 // the maximum number of corners in advance (2^dim)
37 template< int mydim, int cdim >
38 struct CornerStorage
39 {
40 using Type = Dune::ReservedVector< Dune::FieldVector< ct, cdim >, (1<<mydim)+1>;
41 };
42 };
43
44 namespace Detail::PQ1Bubble {
45
46 template<Dune::GeometryType::Id gt>
47 struct OverlappingScvCorners;
48
49 template<>
50 struct OverlappingScvCorners<Dune::GeometryTypes::line>
51 {
52 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
53 static constexpr std::array<std::array<Key, 2>, 1> keys = {{
54 { Key{0, 1}, Key{1, 1} }
55 }};
56 };
57
58 template<>
59 struct OverlappingScvCorners<Dune::GeometryTypes::triangle>
60 {
61 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
62 static constexpr std::array<std::array<Key, 3>, 1> keys = {{
63 { Key{0, 1}, Key{1, 1}, Key{2, 1} }
64 }};
65 };
66
67 template<>
68 struct OverlappingScvCorners<Dune::GeometryTypes::quadrilateral>
69 {
70 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
71 static constexpr std::array<std::array<Key, 4>, 1> keys = {{
72 { Key{2, 1}, Key{1, 1}, Key{0, 1}, Key{3, 1} }
73 }};
74 };
75
76 template<>
77 struct OverlappingScvCorners<Dune::GeometryTypes::tetrahedron>
78 {
79 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
80 static constexpr std::array<std::array<Key, 4>, 1> keys = {{
81 { Key{0, 1}, Key{1, 1}, Key{2, 1}, Key{3, 1} }
82 }};
83 };
84
85 template<>
86 struct OverlappingScvCorners<Dune::GeometryTypes::hexahedron>
87 {
88 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
89 static constexpr std::array<std::array<Key, 6>, 1> keys = {{
90 { Key{0, 1}, Key{2, 1}, Key{3, 1}, Key{1, 1}, Key{4, 1}, Key{5, 1} }
91 }};
92 };
93
94
95 template<Dune::GeometryType::Id gt>
96 struct OverlappingScvfCorners;
97
98 template<>
99 struct OverlappingScvfCorners<Dune::GeometryTypes::line>
100 {
101 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
102 static constexpr std::array<std::array<Key, 1>, 1> keys = {{
103 { Key{0, 0} }
104 }};
105 };
106
107 template<>
108 struct OverlappingScvfCorners<Dune::GeometryTypes::triangle>
109 {
110 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
111 static constexpr std::array<std::array<Key, 2>, 3> keys = {{
112 { Key{0, 1}, Key{1, 1} },
113 { Key{0, 1}, Key{2, 1} },
114 { Key{1, 1}, Key{2, 1} }
115 }};
116 };
117
118 template<>
119 struct OverlappingScvfCorners<Dune::GeometryTypes::quadrilateral>
120 {
121 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
122 static constexpr std::array<std::array<Key, 2>, 4> keys = {{
123 { Key{0, 1}, Key{2, 1} },
124 { Key{2, 1}, Key{1, 1} },
125 { Key{0, 1}, Key{3, 1} },
126 { Key{1, 1}, Key{3, 1} }
127 }};
128 };
129
130 template<>
131 struct OverlappingScvfCorners<Dune::GeometryTypes::tetrahedron>
132 {
133 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
134 static constexpr std::array<std::array<Key, 3>, 10> keys = {{
135 { Key{0, 1}, Key{1, 1}, Key{2, 1} },
136 { Key{0, 1}, Key{1, 1}, Key{3, 1} },
137 { Key{0, 1}, Key{2, 1}, Key{3, 1} },
138 { Key{1, 1}, Key{2, 1}, Key{3, 1} }
139 }};
140 };
141
142 template<>
143 struct OverlappingScvfCorners<Dune::GeometryTypes::hexahedron>
144 {
145 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
146 static constexpr std::array<std::array<Key, 3>, 8> keys = {{
147 { Key{4, 1}, Key{0, 1}, Key{2, 1} },
148 { Key{4, 1}, Key{2, 1}, Key{1, 1} },
149 { Key{4, 1}, Key{0, 1}, Key{3, 1} },
150 { Key{4, 1}, Key{1, 1}, Key{3, 1} },
151 { Key{5, 1}, Key{0, 1}, Key{2, 1} },
152 { Key{5, 1}, Key{2, 1}, Key{1, 1} },
153 { Key{5, 1}, Key{0, 1}, Key{3, 1} },
154 { Key{5, 1}, Key{1, 1}, Key{3, 1} }
155 }};
156 };
157
158 } // end namespace Detail::PQ1Bubble
159
160 /*!
161 * \ingroup PQ1BubbleDiscretization
162 * \brief A class to create sub control volume and sub control volume face geometries per element
163 */
164 template <class GridView, class ScvType, class ScvfType>
165 class PQ1BubbleGeometryHelper
166 {
167 using Scalar = typename GridView::ctype;
168 using GlobalPosition = typename Dune::FieldVector<Scalar, GridView::dimensionworld>;
169 using ScvCornerStorage = typename ScvType::Traits::CornerStorage;
170 using ScvfCornerStorage = typename ScvfType::Traits::CornerStorage;
171 using LocalIndexType = typename ScvType::Traits::LocalIndexType;
172
173 using Element = typename GridView::template Codim<0>::Entity;
174 using Intersection = typename GridView::Intersection;
175
176 static constexpr auto dim = GridView::dimension;
177 static constexpr auto dimWorld = GridView::dimensionworld;
178 public:
179
180 PQ1BubbleGeometryHelper(const typename Element::Geometry& geometry)
181 : geo_(geometry)
182
2/4
✓ Branch 1 taken 12000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12000 times.
✗ Branch 5 not taken.
190010 , boxHelper_(geometry)
183 {}
184
185 //! Create a vector with the scv corners
186 432029 ScvCornerStorage getScvCorners(unsigned int localScvIdx) const
187 {
188 // proceed according to number of corners of the element
189 432029 const auto type = geo_.type();
190 432029 const auto numBoxScv = boxHelper_.numScv();
191 // reuse box geometry helper for the corner scvs
192
2/2
✓ Branch 0 taken 337716 times.
✓ Branch 1 taken 94313 times.
432029 if (localScvIdx < numBoxScv)
193 337716 return boxHelper_.getScvCorners(localScvIdx);
194
195 94313 const auto localOverlappingScvIdx = localScvIdx-numBoxScv;
196
1/2
✓ Branch 0 taken 94313 times.
✗ Branch 1 not taken.
94313 if (type == Dune::GeometryTypes::triangle)
197 {
198 using Corners = Detail::PQ1Bubble::OverlappingScvCorners<Dune::GeometryTypes::triangle>;
199 79096 return Detail::Box::keyToCornerStorage<ScvCornerStorage>(geo_, Corners::keys[localOverlappingScvIdx]);
200 }
201
1/2
✓ Branch 0 taken 54765 times.
✗ Branch 1 not taken.
54765 else if (type == Dune::GeometryTypes::quadrilateral)
202 {
203 using Corners = Detail::PQ1Bubble::OverlappingScvCorners<Dune::GeometryTypes::quadrilateral>;
204 100704 return Detail::Box::keyToCornerStorage<ScvCornerStorage>(geo_, Corners::keys[localOverlappingScvIdx]);
205 }
206
1/2
✓ Branch 0 taken 4413 times.
✗ Branch 1 not taken.
4413 else if (type == Dune::GeometryTypes::tetrahedron)
207 {
208 using Corners = Detail::PQ1Bubble::OverlappingScvCorners<Dune::GeometryTypes::tetrahedron>;
209 8824 return Detail::Box::keyToCornerStorage<ScvCornerStorage>(geo_, Corners::keys[localOverlappingScvIdx]);
210 }
211
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
1 else if (type == Dune::GeometryTypes::hexahedron)
212 {
213 using Corners = Detail::PQ1Bubble::OverlappingScvCorners<Dune::GeometryTypes::hexahedron>;
214 2 return Detail::Box::keyToCornerStorage<ScvCornerStorage>(geo_, Corners::keys[localOverlappingScvIdx]);
215 }
216 else
217 DUNE_THROW(Dune::NotImplemented, "PQ1Bubble scv geometries for dim=" << dim
218 << " dimWorld=" << dimWorld
219 << " type=" << type);
220 }
221
222 432029 Dune::GeometryType getScvGeometryType(unsigned int localScvIdx) const
223 {
224 // proceed according to number of corners of the element
225 432029 const auto type = geo_.type();
226 432029 const auto numBoxScv = boxHelper_.numScv();
227
2/2
✓ Branch 0 taken 337716 times.
✓ Branch 1 taken 94313 times.
432029 if (localScvIdx < numBoxScv)
228 337716 return Dune::GeometryTypes::cube(dim);
229
1/2
✓ Branch 0 taken 94313 times.
✗ Branch 1 not taken.
94313 else if (type == Dune::GeometryTypes::simplex(dim))
230 43960 return Dune::GeometryTypes::simplex(dim);
231
1/2
✓ Branch 0 taken 50353 times.
✗ Branch 1 not taken.
50353 else if (type == Dune::GeometryTypes::quadrilateral)
232 50352 return Dune::GeometryTypes::quadrilateral;
233
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
1 else if (type == Dune::GeometryTypes::hexahedron)
234 1 return Dune::GeometryTypes::none(dim); // octahedron
235 else
236 DUNE_THROW(Dune::NotImplemented, "PQ1Bubble scv geometries for dim=" << dim
237 << " dimWorld=" << dimWorld
238 << " type=" << type);
239 }
240
241 //! Create a vector with the corners of sub control volume faces
242 684264 ScvfCornerStorage getScvfCorners(unsigned int localScvfIdx) const
243 {
244 // proceed according to number of corners
245 684264 const auto type = geo_.type();
246 684264 const auto numBoxScvf = boxHelper_.numInteriorScvf();
247 // reuse box geometry helper for the corner scvs
248
2/2
✓ Branch 0 taken 346548 times.
✓ Branch 1 taken 337716 times.
684264 if (localScvfIdx < numBoxScvf)
249 346548 return boxHelper_.getScvfCorners(localScvfIdx);
250
251 337716 const auto localOverlappingScvfIdx = localScvfIdx-numBoxScvf;
252
1/2
✓ Branch 0 taken 337716 times.
✗ Branch 1 not taken.
337716 if (type == Dune::GeometryTypes::triangle)
253 {
254 using Corners = Detail::PQ1Bubble::OverlappingScvfCorners<Dune::GeometryTypes::triangle>;
255 237288 return Detail::Box::keyToCornerStorage<ScvfCornerStorage>(geo_, Corners::keys[localOverlappingScvfIdx]);
256 }
257
1/2
✓ Branch 0 taken 219072 times.
✗ Branch 1 not taken.
219072 else if (type == Dune::GeometryTypes::quadrilateral)
258 {
259 using Corners = Detail::PQ1Bubble::OverlappingScvfCorners<Dune::GeometryTypes::quadrilateral>;
260 402816 return Detail::Box::keyToCornerStorage<ScvfCornerStorage>(geo_, Corners::keys[localOverlappingScvfIdx]);
261 }
262
1/2
✓ Branch 0 taken 17664 times.
✗ Branch 1 not taken.
17664 else if (type == Dune::GeometryTypes::tetrahedron)
263 {
264 using Corners = Detail::PQ1Bubble::OverlappingScvfCorners<Dune::GeometryTypes::tetrahedron>;
265 35296 return Detail::Box::keyToCornerStorage<ScvfCornerStorage>(geo_, Corners::keys[localOverlappingScvfIdx]);
266 }
267
1/2
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
16 else if (type == Dune::GeometryTypes::hexahedron)
268 {
269 using Corners = Detail::PQ1Bubble::OverlappingScvfCorners<Dune::GeometryTypes::hexahedron>;
270 32 return Detail::Box::keyToCornerStorage<ScvfCornerStorage>(geo_, Corners::keys[localOverlappingScvfIdx]);
271 }
272 else
273 DUNE_THROW(Dune::NotImplemented, "PQ1Bubble scvf geometries for dim=" << dim
274 << " dimWorld=" << dimWorld
275 << " type=" << type);
276 }
277
278 440224 Dune::GeometryType getInteriorScvfGeometryType(unsigned int localScvfIdx) const
279 {
280
1/2
✓ Branch 2 taken 75200 times.
✗ Branch 3 not taken.
684264 const auto numBoxScvf = boxHelper_.numInteriorScvf();
281
6/6
✓ Branch 0 taken 224524 times.
✓ Branch 1 taken 215700 times.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 12 times.
✓ Branch 4 taken 8 times.
✓ Branch 5 taken 12 times.
440264 if (localScvfIdx < numBoxScvf)
282 224524 return Dune::GeometryTypes::cube(dim-1);
283 else
284 215700 return Dune::GeometryTypes::simplex(dim-1);
285 }
286
287 //! Create the sub control volume face geometries on the boundary
288 ScvfCornerStorage getBoundaryScvfCorners(unsigned int localFacetIndex,
289 unsigned int indexInFacet) const
290 {
291
1/2
✓ Branch 1 taken 12882 times.
✗ Branch 2 not taken.
16730 return boxHelper_.getBoundaryScvfCorners(localFacetIndex, indexInFacet);
292 }
293
294 Dune::GeometryType getBoundaryScvfGeometryType(unsigned int localScvfIdx) const
295 {
296 return Dune::GeometryTypes::cube(dim-1);
297 }
298
299 template<int d = dimWorld, std::enable_if_t<(d==3), int> = 0>
300 44140 GlobalPosition normal(const ScvfCornerStorage& p, const std::array<LocalIndexType, 2>& scvPair)
301 {
302 308980 auto normal = Dumux::crossProduct(p[1]-p[0], p[2]-p[0]);
303 88280 normal /= normal.two_norm();
304
305 88280 GlobalPosition v = dofPosition(scvPair[1]) - dofPosition(scvPair[0]);
306
307 44140 const auto s = v*normal;
308
4/4
✓ Branch 0 taken 8828 times.
✓ Branch 1 taken 35312 times.
✓ Branch 2 taken 8828 times.
✓ Branch 3 taken 35312 times.
88280 if (std::signbit(s))
309 normal *= -1;
310
311 44140 return normal;
312 }
313
314 template<int d = dimWorld, std::enable_if_t<(d==2), int> = 0>
315 640104 GlobalPosition normal(const ScvfCornerStorage& p, const std::array<LocalIndexType, 2>& scvPair)
316 {
317 //! obtain normal vector by 90° counter-clockwise rotation of t
318 1920312 const auto t = p[1] - p[0];
319 1920312 GlobalPosition normal({-t[1], t[0]});
320 1280208 normal /= normal.two_norm();
321
322 1280208 GlobalPosition v = dofPosition(scvPair[1]) - dofPosition(scvPair[0]);
323
324 640104 const auto s = v*normal;
325
4/4
✓ Branch 0 taken 232204 times.
✓ Branch 1 taken 407900 times.
✓ Branch 2 taken 232204 times.
✓ Branch 3 taken 407900 times.
1280208 if (std::signbit(s))
326 normal *= -1;
327
328 640104 return normal;
329 }
330
331 //! the wrapped element geometry
332 const typename Element::Geometry& elementGeometry() const
333 { return geo_; }
334
335 //! number of interior sub control volume faces
336 624500 std::size_t numInteriorScvf() const
337 {
338 624500 return boxHelper_.numInteriorScvf() + referenceElement(geo_).size(dim);
339 }
340
341 //! number of boundary sub control volume faces for face localFacetIndex
342 std::size_t numBoundaryScvf(unsigned int localFacetIndex) const
343 {
344 return referenceElement(geo_).size(localFacetIndex, 1, dim);
345 }
346
347 //! number of sub control volumes (number of codim-1 entities)
348 std::size_t numScv() const
349 {
350
9/13
✓ Branch 2 taken 12000 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 288912 times.
✓ Branch 6 taken 166721 times.
✓ Branch 7 taken 21101 times.
✓ Branch 8 taken 12000 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 12000 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 61600 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 49600 times.
✓ Branch 17 taken 12000 times.
526334 return boxHelper_.numScv() + 1;
351 }
352
353 //! get scv volume
354 22069 Scalar scvVolume(unsigned int localScvIdx, const ScvCornerStorage& p) const
355 {
356
1/2
✓ Branch 1 taken 49600 times.
✗ Branch 2 not taken.
432021 const auto scvType = getScvGeometryType(localScvIdx);
357 if constexpr (dim == 3)
358
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 22068 times.
22069 if (scvType == Dune::GeometryTypes::none(dim))
359 1 return octahedronVolume_(p);
360
361 432020 return Dumux::convexPolytopeVolume<dim>(
362 scvType,
363 3826928 [&](unsigned int i){ return p[i]; }
364 );
365 }
366
367 template<class DofMapper>
368 276912 auto dofIndex(const DofMapper& dofMapper, const Element& element, unsigned int localScvIdx) const
369 {
370
2/2
✓ Branch 1 taken 215700 times.
✓ Branch 2 taken 61212 times.
276912 if (localScvIdx < numScv()-1)
371 215700 return dofMapper.subIndex(element, localScvIdx, dim);
372 else
373 61212 return dofMapper.index(element);
374 }
375
376 1600509 GlobalPosition dofPosition(unsigned int localScvIdx) const
377 {
378
2/2
✓ Branch 1 taken 1218088 times.
✓ Branch 2 taken 382421 times.
1600509 if (localScvIdx < numScv()-1)
379 1218088 return geo_.corner(localScvIdx);
380 else
381 382421 return geo_.center();
382 }
383
384 684244 std::array<LocalIndexType, 2> getScvPairForScvf(unsigned int localScvfIndex) const
385 {
386
2/2
✓ Branch 1 taken 346536 times.
✓ Branch 2 taken 337708 times.
684244 const auto numEdges = referenceElement(geo_).size(dim-1);
387
2/2
✓ Branch 0 taken 346536 times.
✓ Branch 1 taken 337708 times.
684244 if (localScvfIndex < numEdges)
388 return {
389 346536 static_cast<LocalIndexType>(referenceElement(geo_).subEntity(localScvfIndex, dim-1, 0, dim)),
390 346536 static_cast<LocalIndexType>(referenceElement(geo_).subEntity(localScvfIndex, dim-1, 1, dim))
391 693072 };
392 else
393 return {
394 337708 static_cast<LocalIndexType>(numScv()-1),
395 337708 static_cast<LocalIndexType>(localScvfIndex-numEdges)
396 337708 };
397 }
398
399 std::array<LocalIndexType, 2> getScvPairForBoundaryScvf(unsigned int localFacetIndex, unsigned int localIsScvfIndex) const
400 {
401 const LocalIndexType insideScvIdx
402 = static_cast<LocalIndexType>(referenceElement(geo_).subEntity(localFacetIndex, 1, localIsScvfIndex, dim));
403 return { insideScvIdx, insideScvIdx };
404 }
405
406 bool isOverlappingScvf(unsigned int localScvfIndex) const
407 {
408
6/6
✓ Branch 1 taken 215700 times.
✓ Branch 2 taken 224524 times.
✓ Branch 4 taken 159608 times.
✓ Branch 5 taken 84412 times.
✓ Branch 6 taken 37600 times.
✓ Branch 7 taken 37600 times.
684244 if (localScvfIndex < boxHelper_.numInteriorScvf())
409 return false;
410 else
411 return true;
412 }
413
414 bool isOverlappingBoundaryScvf(unsigned int localFacetIndex) const
415 {
416 return false;
417 }
418
419 bool isOverlappingScv(unsigned int localScvIndex) const
420 {
421
6/6
✓ Branch 1 taken 61212 times.
✓ Branch 2 taken 215700 times.
✓ Branch 4 taken 70701 times.
✓ Branch 5 taken 84408 times.
✓ Branch 6 taken 12000 times.
✓ Branch 7 taken 37600 times.
432021 if (localScvIndex < boxHelper_.numScv())
422 return false;
423 else
424 return true;
425 }
426
427 private:
428 1 Scalar octahedronVolume_(const ScvCornerStorage& p) const
429 {
430 using std::abs;
431 return 1.0/6.0 * (
432 11 abs(Dumux::tripleProduct(p[4]-p[0], p[1]-p[0], p[2]-p[0]))
433 10 + abs(Dumux::tripleProduct(p[5]-p[0], p[1]-p[0], p[2]-p[0]))
434 2 );
435 }
436
437 const typename Element::Geometry& geo_; //!< Reference to the element geometry
438 Dumux::BoxGeometryHelper<GridView, dim, ScvType, ScvfType> boxHelper_;
439 };
440
441 } // end namespace Dumux
442
443 #endif
444