GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/discretization/box/boxgeometryhelper.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 89 107 83.2%
Functions: 165 269 61.3%
Branches: 213 411 51.8%
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 BoxDiscretization
10 * \brief Helper class constructing the dual grid finite volume geometries
11 * for the box discretizazion method
12 */
13 #ifndef DUMUX_DISCRETIZATION_BOX_GEOMETRY_HELPER_HH
14 #define DUMUX_DISCRETIZATION_BOX_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/typeindex.hh>
22 #include <dune/geometry/referenceelements.hh>
23 #include <dune/geometry/multilineargeometry.hh>
24
25 #include <dumux/common/math.hh>
26
27 namespace Dumux {
28
29 //! Traits for an efficient corner storage for box method sub control volumes
30 template <class ct>
31 struct BoxMLGeometryTraits : public Dune::MultiLinearGeometryTraits<ct>
32 {
33 // we use static vectors to store the corners as we know
34 // the number of corners in advance (2^(mydim) corners (1<<(mydim))
35 template< int mydim, int cdim >
36 struct CornerStorage
37 {
38 using Type = std::array< Dune::FieldVector< ct, cdim >, (1<<(mydim)) >;
39 };
40
41 // we know all scvfs will have the same geometry type
42 template< int mydim >
43 struct hasSingleGeometryType
44 {
45 static const bool v = true;
46 static const unsigned int topologyId = Dune::GeometryTypes::cube(mydim).id();
47 };
48 };
49
50 namespace Detail::Box {
51
52 template<Dune::GeometryType::Id gt>
53 struct ScvCorners;
54
55 template<>
56 struct ScvCorners<Dune::GeometryTypes::line>
57 {
58 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
59 static constexpr std::array<std::array<Key, 2>, 2> keys = {{
60 { Key{0, 1}, Key{0, 0} },
61 { Key{1, 1}, Key{0, 0} }
62 }};
63 };
64
65 template<>
66 struct ScvCorners<Dune::GeometryTypes::triangle>
67 {
68 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
69 static constexpr std::array<std::array<Key, 4>, 3> keys = {{
70 { Key{0, 2}, Key{0, 1}, Key{1, 1}, Key{0, 0} },
71 { Key{1, 2}, Key{2, 1}, Key{0, 1}, Key{0, 0} },
72 { Key{2, 2}, Key{1, 1}, Key{2, 1}, Key{0, 0} }
73 }};
74 };
75
76 template<>
77 struct ScvCorners<Dune::GeometryTypes::quadrilateral>
78 {
79 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
80 static constexpr std::array<std::array<Key, 4>, 4> keys = {{
81 { Key{0, 2}, Key{2, 1}, Key{0, 1}, Key{0, 0} },
82 { Key{1, 2}, Key{1, 1}, Key{2, 1}, Key{0, 0} },
83 { Key{2, 2}, Key{0, 1}, Key{3, 1}, Key{0, 0} },
84 { Key{3, 2}, Key{3, 1}, Key{1, 1}, Key{0, 0} }
85 }};
86 };
87
88 template<>
89 struct ScvCorners<Dune::GeometryTypes::tetrahedron>
90 {
91 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
92 static constexpr std::array<std::array<Key, 8>, 4> keys = {{
93 { Key{0, 3}, Key{0, 2}, Key{1, 2}, Key{0, 1}, Key{3, 2}, Key{1, 1}, Key{2, 1}, Key{0, 0} },
94 { Key{1, 3}, Key{2, 2}, Key{0, 2}, Key{0, 1}, Key{4, 2}, Key{3, 1}, Key{1, 1}, Key{0, 0} },
95 { Key{2, 3}, Key{1, 2}, Key{2, 2}, Key{0, 1}, Key{5, 2}, Key{2, 1}, Key{3, 1}, Key{0, 0} },
96 { Key{3, 3}, Key{3, 2}, Key{5, 2}, Key{2, 1}, Key{4, 2}, Key{1, 1}, Key{3, 1}, Key{0, 0} }
97 }};
98 };
99
100 template<>
101 struct ScvCorners<Dune::GeometryTypes::prism>
102 {
103 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
104 static constexpr std::array<std::array<Key, 8>, 6> keys = {{
105 { Key{0, 3}, Key{3, 2}, Key{4, 2}, Key{3, 1}, Key{0, 2}, Key{0, 1}, Key{1, 1}, Key{0, 0} },
106 { Key{1, 3}, Key{5, 2}, Key{3, 2}, Key{3, 1}, Key{1, 2}, Key{2, 1}, Key{0, 1}, Key{0, 0} },
107 { Key{2, 3}, Key{4, 2}, Key{5, 2}, Key{3, 1}, Key{2, 2}, Key{1, 1}, Key{2, 1}, Key{0, 0} },
108 { Key{3, 3}, Key{7, 2}, Key{6, 2}, Key{4, 1}, Key{0, 2}, Key{1, 1}, Key{0, 1}, Key{0, 0} },
109 { Key{4, 3}, Key{6, 2}, Key{8, 2}, Key{4, 1}, Key{1, 2}, Key{0, 1}, Key{2, 1}, Key{0, 0} },
110 { Key{5, 3}, Key{8, 2}, Key{7, 2}, Key{4, 1}, Key{2, 2}, Key{2, 1}, Key{1, 1}, Key{0, 0} }
111 }};
112 };
113
114 template<>
115 struct ScvCorners<Dune::GeometryTypes::hexahedron>
116 {
117 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
118 static constexpr std::array<std::array<Key, 8>, 8> keys = {{
119 { Key{0, 3}, Key{6, 2}, Key{4, 2}, Key{4, 1}, Key{0, 2}, Key{2, 1}, Key{0, 1}, Key{0, 0} },
120 { Key{1, 3}, Key{5, 2}, Key{6, 2}, Key{4, 1}, Key{1, 2}, Key{1, 1}, Key{2, 1}, Key{0, 0} },
121 { Key{2, 3}, Key{4, 2}, Key{7, 2}, Key{4, 1}, Key{2, 2}, Key{0, 1}, Key{3, 1}, Key{0, 0} },
122 { Key{3, 3}, Key{7, 2}, Key{5, 2}, Key{4, 1}, Key{3, 2}, Key{3, 1}, Key{1, 1}, Key{0, 0} },
123 { Key{4, 3}, Key{8, 2}, Key{10, 2}, Key{5, 1}, Key{0, 2}, Key{0, 1}, Key{2, 1}, Key{0, 0} },
124 { Key{5, 3}, Key{10, 2}, Key{9, 2}, Key{5, 1}, Key{1, 2}, Key{2, 1}, Key{1, 1}, Key{0, 0} },
125 { Key{6, 3}, Key{11, 2}, Key{8, 2}, Key{5, 1}, Key{2, 2}, Key{3, 1}, Key{0, 1}, Key{0, 0} },
126 { Key{7, 3}, Key{9, 2}, Key{11, 2}, Key{5, 1}, Key{3, 2}, Key{1, 1}, Key{3, 1}, Key{0, 0} }
127 }};
128 };
129
130 template<Dune::GeometryType::Id gt>
131 struct ScvfCorners;
132
133 template<>
134 struct ScvfCorners<Dune::GeometryTypes::line>
135 {
136 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
137 static constexpr std::array<std::array<Key, 1>, 1> keys = {{
138 { Key{0, 0} }
139 }};
140 };
141
142 template<>
143 struct ScvfCorners<Dune::GeometryTypes::triangle>
144 {
145 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
146 static constexpr std::array<std::array<Key, 2>, 3> keys = {{
147 { Key{0, 0}, Key{0, 1} },
148 { Key{1, 1}, Key{0, 0} },
149 { Key{0, 0}, Key{2, 1} }
150 }};
151 };
152
153 template<>
154 struct ScvfCorners<Dune::GeometryTypes::quadrilateral>
155 {
156 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
157 static constexpr std::array<std::array<Key, 2>, 4> keys = {{
158 { Key{0, 1}, Key{0, 0} },
159 { Key{0, 0}, Key{1, 1} },
160 { Key{0, 0}, Key{2, 1} },
161 { Key{3, 1}, Key{0, 0} }
162 }};
163 };
164
165 template<>
166 struct ScvfCorners<Dune::GeometryTypes::tetrahedron>
167 {
168 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
169 static constexpr std::array<std::array<Key, 4>, 6> keys = {{
170 { Key{0, 2}, Key{0, 1}, Key{1, 1}, Key{0, 0} },
171 { Key{0, 1}, Key{1, 2}, Key{0, 0}, Key{2, 1} },
172 { Key{2, 2}, Key{0, 1}, Key{3, 1}, Key{0, 0} },
173 { Key{2, 1}, Key{3, 2}, Key{0, 0}, Key{1, 1} },
174 { Key{3, 1}, Key{0, 0}, Key{4, 2}, Key{1, 1} },
175 { Key{5, 2}, Key{2, 1}, Key{3, 1}, Key{0, 0} }
176 }};
177 };
178
179 template<>
180 struct ScvfCorners<Dune::GeometryTypes::prism>
181 {
182 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
183 static constexpr std::array<std::array<Key, 4>, 9> keys = {{
184 { Key{0, 2}, Key{0, 1}, Key{1, 1}, Key{0, 0} },
185 { Key{1, 2}, Key{2, 1}, Key{0, 1}, Key{0, 0} },
186 { Key{2, 2}, Key{1, 1}, Key{2, 1}, Key{0, 0} },
187 { Key{3, 2}, Key{0, 1}, Key{3, 1}, Key{0, 0} },
188 { Key{4, 2}, Key{3, 1}, Key{1, 1}, Key{0, 0} },
189 { Key{5, 2}, Key{2, 1}, Key{3, 1}, Key{0, 0} },
190 { Key{6, 2}, Key{4, 1}, Key{0, 1}, Key{0, 0} },
191 { Key{7, 2}, Key{1, 1}, Key{4, 1}, Key{0, 0} },
192 { Key{8, 2}, Key{4, 1}, Key{2, 1}, Key{0, 0} }
193 }};
194 };
195
196 template<>
197 struct ScvfCorners<Dune::GeometryTypes::hexahedron>
198 {
199 using Key = std::pair<std::uint8_t, std::uint8_t>; // (i, codim)
200 static constexpr std::array<std::array<Key, 4>, 12> keys = {{
201 { Key{0, 1}, Key{0, 2}, Key{0, 0}, Key{2, 1} },
202 { Key{1, 1}, Key{0, 0}, Key{1, 2}, Key{2, 1} },
203 { Key{3, 1}, Key{2, 2}, Key{0, 0}, Key{0, 1} },
204 { Key{3, 2}, Key{3, 1}, Key{1, 1}, Key{0, 0} },
205 { Key{4, 1}, Key{4, 2}, Key{0, 0}, Key{0, 1} },
206 { Key{5, 2}, Key{4, 1}, Key{1, 1}, Key{0, 0} },
207 { Key{6, 2}, Key{4, 1}, Key{2, 1}, Key{0, 0} },
208 { Key{4, 1}, Key{7, 2}, Key{0, 0}, Key{3, 1} },
209 { Key{0, 0}, Key{0, 1}, Key{5, 1}, Key{8, 2} },
210 { Key{9, 2}, Key{1, 1}, Key{5, 1}, Key{0, 0} },
211 { Key{10, 2}, Key{2, 1}, Key{5, 1}, Key{0, 0} },
212 { Key{11, 2}, Key{5, 1}, Key{3, 1}, Key{0, 0} }
213 }};
214 };
215
216 // convert key array to global corner storage
217 template<class S, class Geo, class KeyArray, std::size_t... I>
218 176725542 S keyToCornerStorageImpl(const Geo& geo, const KeyArray& key, std::index_sequence<I...>)
219 {
220 using Dune::referenceElement;
221
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 76 times.
 176725542 const auto ref = referenceElement(geo);
222 // key is a pair of a local sub-entity index (first) and the sub-entity's codim (second)
223
36/54
✗ Branch 0 not taken.
✓ Branch 1 taken 88362771 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 88362771 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 88362771 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 24470613 times.
✓ Branch 8 taken 61224832 times.
✓ Branch 9 taken 24470613 times.
✓ Branch 10 taken 61224832 times.
✓ Branch 11 taken 24470613 times.
✓ Branch 12 taken 61224832 times.
✓ Branch 13 taken 14698043 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 42176424 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 42176424 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 42176408 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 9836471 times.
✓ Branch 22 taken 32253793 times.
✓ Branch 23 taken 9836471 times.
✓ Branch 24 taken 32253793 times.
✓ Branch 25 taken 251121 times.
✓ Branch 26 taken 32253793 times.
✓ Branch 27 taken 251121 times.
✗ Branch 28 not taken.
✓ Branch 29 taken 1719625 times.
✗ Branch 30 not taken.
✓ Branch 31 taken 1719625 times.
✗ Branch 32 not taken.
✓ Branch 33 taken 1719625 times.
✗ Branch 34 not taken.
✓ Branch 35 taken 251121 times.
✓ Branch 36 taken 1468504 times.
✓ Branch 37 taken 251120 times.
✓ Branch 38 taken 1468504 times.
✓ Branch 39 taken 251120 times.
✓ Branch 40 taken 1468504 times.
✓ Branch 41 taken 251120 times.
✗ Branch 42 not taken.
✓ Branch 43 taken 1719624 times.
✗ Branch 44 not taken.
✓ Branch 45 taken 1719624 times.
✗ Branch 46 not taken.
✓ Branch 47 taken 1719624 times.
✗ Branch 49 not taken.
✓ Branch 50 taken 1468504 times.
✗ Branch 51 not taken.
✓ Branch 52 taken 1468504 times.
✗ Branch 53 not taken.
✓ Branch 54 taken 1468504 times.
 772444322 return { geo.global(ref.position(key[I].first, key[I].second))... };
224 }
225
226 // convert key array to global corner storage
227 template<class S, class Geo, class T, std::size_t N, class Indices = std::make_index_sequence<N>>
228 S keyToCornerStorage(const Geo& geo, const std::array<T, N>& key)
229 {
230 87070123 return keyToCornerStorageImpl<S>(geo, key, Indices{});
231 }
232
233 // convert key array to global corner storage
234 // for the i-th sub-entity of codim c (e.g. the i-th facet/codim-1-entity for boundaries)
235 template<class S, class Geo, class KeyArray, std::size_t... I>
236 5771106 S subEntityKeyToCornerStorageImpl(const Geo& geo, unsigned int i, unsigned int c, const KeyArray& key, std::index_sequence<I...>)
237 {
238 using Dune::referenceElement;
239 5771106 const auto ref = referenceElement(geo);
240 // subEntity gives the subEntity number with respect to the codim-0 reference element
241 // key is a pair of a local sub-entity index (first) and the sub-entity's codim (second) but here w.r.t. the sub-entity i/c
242
7/13
✗ Branch 4 not taken.
✓ Branch 5 taken 5218024 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 2079240 times.
✓ Branch 12 taken 3138784 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 133424 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 490170 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 133424 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 490170 times.
 32583856 return { geo.global(ref.position(ref.subEntity(i, c, key[I].first, c+key[I].second), c+key[I].second))... };
243 }
244
245 // convert key array to global corner storage
246 // for the i-th sub-entity of codim c (e.g. the i-th facet/codim-1-entity for boundaries)
247 template<class S, class Geo, class T, std::size_t N, class Indices = std::make_index_sequence<N>>
248 S subEntityKeyToCornerStorage(const Geo& geo, unsigned int i, unsigned int c, const std::array<T, N>& key)
249 {
250 4828684 return subEntityKeyToCornerStorageImpl<S>(geo, i, c, key, Indices{});
251 }
252
253 } // end namespace Detail::Box
254
255 //! Create sub control volumes and sub control volume face geometries
256 template<class GridView, int dim, class ScvType, class ScvfType>
257 class BoxGeometryHelper;
258
259 //! A class to create sub control volume and sub control volume face geometries per element
260 template <class GridView, class ScvType, class ScvfType>
261 class BoxGeometryHelper<GridView, 1, ScvType, ScvfType>
262 {
263 private:
264 using Scalar = typename GridView::ctype;
265 using GlobalPosition = typename Dune::FieldVector<Scalar, GridView::dimensionworld>;
266 using ScvCornerStorage = typename ScvType::Traits::CornerStorage;
267 using ScvfCornerStorage = typename ScvfType::Traits::CornerStorage;
268 using ScvGeometry = typename ScvType::Traits::Geometry;
269 using ScvfGeometry = typename ScvfType::Traits::Geometry;
270
271 using Element = typename GridView::template Codim<0>::Entity;
272 using Intersection = typename GridView::Intersection;
273
274 static constexpr int dim = 1;
275 public:
276
277 2667326 BoxGeometryHelper(const typename Element::Geometry& geometry)
278 2667326 : geo_(geometry)
279 {}
280
281 //! Create a vector with the scv corners
282 ✗ ScvCornerStorage getScvCorners(unsigned int localScvIdx) const
283 {
284 using Corners = Detail::Box::ScvCorners<Dune::GeometryTypes::line>;
285
9/18
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 4 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 4 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 4 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 4 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 4 times.
✗ Branch 26 not taken.
 10670592 return Detail::Box::keyToCornerStorage<ScvCornerStorage>(geo_, Corners::keys[localScvIdx]);
286 }
287
288 ScvGeometry scvGeometry(unsigned int localScvIdx) const
289 {
290 return { Dune::GeometryTypes::line, getScvCorners(localScvIdx) };
291 }
292
293 //! Create a vector with the corners of sub control volume faces
294 ✗ ScvfCornerStorage getScvfCorners(unsigned int localScvfIdx) const
295 {
296 using Corners = Detail::Box::ScvfCorners<Dune::GeometryTypes::line>;
297
2/6
✓ Branch 3 taken 2667326 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2667326 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 5334652 return Detail::Box::keyToCornerStorage<ScvfCornerStorage>(geo_, Corners::keys[localScvfIdx]);
298 }
299
300 //! Create the sub control volume face geometries on the boundary
301 ✗ ScvfCornerStorage getBoundaryScvfCorners(unsigned int localFacetIndex,
302 unsigned int) const
303 {
304
4/7
✗ Branch 0 not taken.
✓ Branch 1 taken 957 times.
✓ Branch 2 taken 24488 times.
✓ Branch 3 taken 42264 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 49050 times.
✗ Branch 6 not taken.
 117716 return ScvfCornerStorage{{ geo_.corner(localFacetIndex) }};
305 }
306
307 //! get scvf normal vector
308 236548 GlobalPosition normal(const ScvfCornerStorage& scvfCorners,
309 const std::vector<unsigned int>& scvIndices) const
310 {
311 7562442 auto normal = geo_.corner(1) - geo_.corner(0);
312 7765430 normal /= normal.two_norm();
313 236548 return normal;
314 }
315
316 //! number of sub control volume faces (number of edges)
317 ✗ std::size_t numInteriorScvf() const
318 {
319 236508 return referenceElement(geo_).size(dim-1);
320 }
321
322 //! number of sub control volumes (number of vertices)
323 std::size_t numScv() const
324 {
325 return referenceElement(geo_).size(dim);
326 }
327
328 //! the wrapped element geometry
329 const typename Element::Geometry& elementGeometry() const
330 { return geo_; }
331
332 private:
333 const typename Element::Geometry& geo_; //!< Reference to the element geometry
334 };
335
336 //! A class to create sub control volume and sub control volume face geometries per element
337 template <class GridView, class ScvType, class ScvfType>
338 class BoxGeometryHelper<GridView, 2, ScvType, ScvfType>
339 {
340 using Scalar = typename GridView::ctype;
341 using GlobalPosition = typename Dune::FieldVector<Scalar, GridView::dimensionworld>;
342 using ScvCornerStorage = typename ScvType::Traits::CornerStorage;
343 using ScvfCornerStorage = typename ScvfType::Traits::CornerStorage;
344
345 using Element = typename GridView::template Codim<0>::Entity;
346 using Intersection = typename GridView::Intersection;
347
348 static constexpr auto dim = GridView::dimension;
349 static constexpr auto dimWorld = GridView::dimensionworld;
350 public:
351
352 10432264 BoxGeometryHelper(const typename Element::Geometry& geometry)
353
3/5
✓ Branch 0 taken 147328 times.
✓ Branch 1 taken 5858865 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7200 times.
✗ Branch 5 not taken.
 10432264 : geo_(geometry)
354 {}
355
356 //! Create a vector with the scv corners
357 7280640 ScvCornerStorage getScvCorners(unsigned int localScvIdx) const
358 {
359 // proceed according to number of corners of the element
360
1/2
✓ Branch 1 taken 17625753 times.
✗ Branch 2 not taken.
 38124421 const auto type = geo_.type();
361
2/3
✓ Branch 0 taken 6933220 times.
✓ Branch 1 taken 17625769 times.
✗ Branch 2 not taken.
 38124421 if (type == Dune::GeometryTypes::triangle)
362 {
363 using Corners = Detail::Box::ScvCorners<Dune::GeometryTypes::triangle>;
364
2/4
✓ Branch 1 taken 3936273 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3936273 times.
✗ Branch 5 not taken.
 24390306 return Detail::Box::keyToCornerStorage<ScvCornerStorage>(geo_, Corners::keys[localScvIdx]);
365 }
366
2/3
✓ Branch 0 taken 2731552 times.
✓ Branch 1 taken 13689496 times.
✗ Branch 2 not taken.
 25929268 else if (type == Dune::GeometryTypes::quadrilateral)
367 {
368 using Corners = Detail::Box::ScvCorners<Dune::GeometryTypes::quadrilateral>;
369
2/4
✓ Branch 1 taken 13689480 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13689480 times.
✗ Branch 5 not taken.
 51858536 return Detail::Box::keyToCornerStorage<ScvCornerStorage>(geo_, Corners::keys[localScvIdx]);
370 }
371 else
372 ✗ DUNE_THROW(Dune::NotImplemented, "Box scv geometries for dim=" << dim
373 << " dimWorld=" << dimWorld
374 << " type=" << type);
375 }
376
377 //! Create a vector with the corners of sub control volume faces
378 7280584 ScvfCornerStorage getScvfCorners(unsigned int localScvfIdx) const
379 {
380 // proceed according to number of corners
381
1/2
✓ Branch 1 taken 21576186 times.
✗ Branch 2 not taken.
 38124205 const auto type = geo_.type();
382
2/3
✓ Branch 0 taken 6933208 times.
✓ Branch 1 taken 21576186 times.
✗ Branch 2 not taken.
 38124205 if (type == Dune::GeometryTypes::triangle)
383 {
384 using Corners = Detail::Box::ScvfCorners<Dune::GeometryTypes::triangle>;
385
2/4
✓ Branch 1 taken 43086 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43086 times.
✗ Branch 5 not taken.
 24390258 return Detail::Box::keyToCornerStorage<ScvfCornerStorage>(geo_, Corners::keys[localScvfIdx]);
386 }
387
2/3
✓ Branch 0 taken 2731552 times.
✓ Branch 1 taken 21533100 times.
✗ Branch 2 not taken.
 25929076 else if (type == Dune::GeometryTypes::quadrilateral)
388 {
389 using Corners = Detail::Box::ScvfCorners<Dune::GeometryTypes::quadrilateral>;
390
2/4
✓ Branch 1 taken 21533100 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 21533100 times.
✗ Branch 5 not taken.
 51858152 return Detail::Box::keyToCornerStorage<ScvfCornerStorage>(geo_, Corners::keys[localScvfIdx]);
391 }
392 else
393 ✗ DUNE_THROW(Dune::NotImplemented, "Box scvf geometries for dim=" << dim
394 << " dimWorld=" << dimWorld
395 << " type=" << type);
396 }
397
398 //! Create the sub control volume face geometries on the boundary
399 ✗ ScvfCornerStorage getBoundaryScvfCorners(unsigned int localFacetIndex,
400 unsigned int indexInFacet) const
401 {
402 // we have to use the corresponding facet geometry as the intersection geometry
403 // might be rotated or flipped. This makes sure that the corners (dof location)
404 // and corresponding scvfs are sorted in the same way
405 using Corners = Detail::Box::ScvCorners<Dune::GeometryTypes::line>;
406 4511270 constexpr int facetCodim = 1;
407
8/13
✓ Branch 1 taken 6048 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2821532 times.
✓ Branch 4 taken 6048 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2821532 times.
✓ Branch 7 taken 3312 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2820812 times.
✓ Branch 10 taken 3248 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3248 times.
✗ Branch 14 not taken.
 9022540 return Detail::Box::subEntityKeyToCornerStorage<ScvfCornerStorage>(geo_, localFacetIndex, facetCodim, Corners::keys[indexInFacet]);
408 }
409
410 //! get scvf normal vector for dim == 2, dimworld == 3
411 template <int w = dimWorld>
412 typename std::enable_if<w == 3, GlobalPosition>::type
413 331566 normal(const ScvfCornerStorage& scvfCorners,
414 const std::vector<unsigned int>& scvIndices) const
415 {
416 385866 const auto v1 = geo_.corner(1) - geo_.corner(0);
417 385866 const auto v2 = geo_.corner(2) - geo_.corner(0);
418 331566 const auto v3 = Dumux::crossProduct(v1, v2);
419 994698 const auto t = scvfCorners[1] - scvfCorners[0];
420 331566 GlobalPosition normal = Dumux::crossProduct(v3, t);
421 663132 normal /= normal.two_norm();
422
423 //! ensure the right direction of the normal
424 385866 const auto v = geo_.corner(scvIndices[1]) - geo_.corner(scvIndices[0]);
425 331566 const auto s = v*normal;
426
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 331566 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 331566 times.
 663132 if (std::signbit(s))
427 normal *= -1;
428
429 331566 return normal;
430 }
431
432 //! get scvf normal vector for dim == 2, dimworld == 2
433 template <int w = dimWorld>
434 typename std::enable_if<w == 2, GlobalPosition>::type
435 37038367 normal(const ScvfCornerStorage& scvfCorners,
436 const std::vector<unsigned int>& scvIndices) const
437 {
438 //! obtain normal vector by 90° counter-clockwise rotation of t
439 111115101 const auto t = scvfCorners[1] - scvfCorners[0];
440 111115101 GlobalPosition normal({-t[1], t[0]});
441 74076734 normal /= normal.two_norm();
442
443 //! ensure the right direction of the normal
444 55244151 const auto v = geo_.corner(scvIndices[1]) - geo_.corner(scvIndices[0]);
445 37038367 const auto s = v*normal;
446
4/4
✓ Branch 0 taken 7614897 times.
✓ Branch 1 taken 29423470 times.
✓ Branch 2 taken 7614897 times.
✓ Branch 3 taken 29423470 times.
 74076734 if (std::signbit(s))
447 normal *= -1;
448
449 37038367 return normal;
450 }
451
452 //! number of sub control volume faces (number of edges)
453 ✗ std::size_t numInteriorScvf() const
454 {
455
16/16
✓ Branch 1 taken 284296 times.
✓ Branch 2 taken 284896 times.
✓ Branch 3 taken 284896 times.
✓ Branch 4 taken 284896 times.
✓ Branch 5 taken 600 times.
✓ Branch 7 taken 284896 times.
✓ Branch 8 taken 284896 times.
✓ Branch 9 taken 284896 times.
✓ Branch 10 taken 284896 times.
✓ Branch 12 taken 284296 times.
✓ Branch 13 taken 284296 times.
✓ Branch 14 taken 284296 times.
✓ Branch 15 taken 284896 times.
✓ Branch 16 taken 600 times.
✓ Branch 17 taken 600 times.
✓ Branch 18 taken 600 times.
 3189014 return referenceElement(geo_).size(dim-1);
456 }
457
458 //! number of sub control volumes (number of vertices)
459 ✗ std::size_t numScv() const
460 {
461
24/24
✓ Branch 1 taken 284896 times.
✓ Branch 2 taken 85748 times.
✓ Branch 3 taken 284896 times.
✓ Branch 4 taken 85748 times.
✓ Branch 8 taken 1139584 times.
✓ Branch 9 taken 370644 times.
✓ Branch 10 taken 1139584 times.
✓ Branch 11 taken 370644 times.
✓ Branch 13 taken 284896 times.
✓ Branch 14 taken 85748 times.
✓ Branch 15 taken 284896 times.
✓ Branch 16 taken 85748 times.
✓ Branch 24 taken 370644 times.
✓ Branch 25 taken 85748 times.
✓ Branch 26 taken 370644 times.
✓ Branch 27 taken 85748 times.
✓ Branch 29 taken 85748 times.
✓ Branch 30 taken 284896 times.
✓ Branch 31 taken 85748 times.
✓ Branch 32 taken 284896 times.
✓ Branch 34 taken 284896 times.
✓ Branch 35 taken 85748 times.
✓ Branch 36 taken 284896 times.
✓ Branch 37 taken 85748 times.
 3734092 return referenceElement(geo_).size(dim);
462 }
463
464 //! the wrapped element geometry
465 ✗ const typename Element::Geometry& elementGeometry() const
466 ✗ { return geo_; }
467
468 private:
469 const typename Element::Geometry& geo_; //!< Reference to the element geometry
470 };
471
472 //! A class to create sub control volume and sub control volume face geometries per element
473 template <class GridView, class ScvType, class ScvfType>
474 class BoxGeometryHelper<GridView, 3, ScvType, ScvfType>
475 {
476 using Scalar = typename GridView::ctype;
477 using GlobalPosition = typename Dune::FieldVector<Scalar, GridView::dimensionworld>;
478 using ScvCornerStorage = typename ScvType::Traits::CornerStorage;
479 using ScvfCornerStorage = typename ScvfType::Traits::CornerStorage;
480
481 using Element = typename GridView::template Codim<0>::Entity;
482 using Intersection = typename GridView::Intersection;
483
484 static constexpr auto dim = GridView::dimension;
485 static constexpr auto dimWorld = GridView::dimensionworld;
486
487 public:
488 453104 BoxGeometryHelper(const typename Element::Geometry& geometry)
489
5/8
✓ Branch 0 taken 67174 times.
✓ Branch 1 taken 157231 times.
✓ Branch 2 taken 8606 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 5120 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 10860 times.
✗ Branch 7 not taken.
 453104 : geo_(geometry)
490 {}
491
492 //! Create a vector with the scv corners
493 794824 ScvCornerStorage getScvCorners(unsigned int localScvIdx) const
494 {
495
3/5
✓ Branch 1 taken 117792 times.
✓ Branch 2 taken 34328 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 9008 times.
✗ Branch 5 not taken.
 1754920 const auto type = geo_.type();
496
4/6
✓ Branch 0 taken 759500 times.
✓ Branch 1 taken 117820 times.
✓ Branch 2 taken 34328 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 9008 times.
✗ Branch 5 not taken.
 1754920 if (type == Dune::GeometryTypes::tetrahedron)
497 {
498 using Corners = Detail::Box::ScvCorners<Dune::GeometryTypes::tetrahedron>;
499
4/8
✓ Branch 1 taken 152120 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 152120 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 9008 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 9008 times.
✗ Branch 11 not taken.
 2658776 return Detail::Box::keyToCornerStorage<ScvCornerStorage>(geo_, Corners::keys[localScvIdx]);
500 }
501
2/2
✓ Branch 0 taken 174396 times.
✓ Branch 1 taken 16 times.
 425532 else if (type == Dune::GeometryTypes::prism)
502 {
503 using Corners = Detail::Box::ScvCorners<Dune::GeometryTypes::prism>;
504 348792 return Detail::Box::keyToCornerStorage<ScvCornerStorage>(geo_, Corners::keys[localScvIdx]);
505 }
506
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
 251136 else if (type == Dune::GeometryTypes::hexahedron)
507 {
508 using Corners = Detail::Box::ScvCorners<Dune::GeometryTypes::hexahedron>;
509 502272 return Detail::Box::keyToCornerStorage<ScvCornerStorage>(geo_, Corners::keys[localScvIdx]);
510 }
511 else
512 ✗ DUNE_THROW(Dune::NotImplemented, "Box scv geometries for dim=" << dim
513 << " dimWorld=" << dimWorld
514 << " type=" << type);
515 }
516
517 //! Create a vector with the scvf corners
518 1192074 ScvfCornerStorage getScvfCorners(unsigned int localScvfIdx) const
519 {
520 // proceed according to number of corners
521
3/7
✓ Branch 1 taken 243564 times.
✓ Branch 2 taken 164616 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 51492 times.
✗ Branch 5 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 2591802 const auto type = geo_.type();
522
4/8
✓ Branch 0 taken 1139130 times.
✓ Branch 1 taken 243564 times.
✓ Branch 2 taken 164616 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 51492 times.
✗ Branch 5 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 2591802 if (type == Dune::GeometryTypes::tetrahedron)
523 {
524 using Corners = Detail::Box::ScvfCorners<Dune::GeometryTypes::tetrahedron>;
525
2/8
✓ Branch 1 taken 109944 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 109944 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 3961116 return Detail::Box::keyToCornerStorage<ScvfCornerStorage>(geo_, Corners::keys[localScvfIdx]);
526 }
527
3/4
✓ Branch 0 taken 261468 times.
✓ Branch 1 taken 185112 times.
✓ Branch 2 taken 164616 times.
✗ Branch 3 not taken.
 611244 else if (type == Dune::GeometryTypes::prism)
528 {
529 using Corners = Detail::Box::ScvfCorners<Dune::GeometryTypes::prism>;
530 522936 return Detail::Box::keyToCornerStorage<ScvfCornerStorage>(geo_, Corners::keys[localScvfIdx]);
531 }
532
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 185112 times.
✓ Branch 2 taken 164616 times.
✗ Branch 3 not taken.
 349776 else if (type == Dune::GeometryTypes::hexahedron)
533 {
534 using Corners = Detail::Box::ScvfCorners<Dune::GeometryTypes::hexahedron>;
535
4/7
✓ Branch 1 taken 185112 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 164616 times.
✓ Branch 4 taken 185112 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 164616 times.
✗ Branch 7 not taken.
 699552 return Detail::Box::keyToCornerStorage<ScvfCornerStorage>(geo_, Corners::keys[localScvfIdx]);
536 }
537 else
538 ✗ DUNE_THROW(Dune::NotImplemented, "Box scvf geometries for dim=" << dim
539 << " dimWorld=" << dimWorld
540 << " type=" << type);
541 }
542
543 //! Create the sub control volume face geometries on the boundary
544 663416 ScvfCornerStorage getBoundaryScvfCorners(unsigned localFacetIndex,
545 unsigned int indexInFacet) const
546 {
547 663416 constexpr int facetCodim = 1;
548
549 // we have to use the corresponding facet geometry as the intersection geometry
550 // might be rotated or flipped. This makes sure that the corners (dof location)
551 // and corresponding scvfs are sorted in the same way
552 using Dune::referenceElement;
553 663416 const auto type = referenceElement(geo_).type(localFacetIndex, facetCodim);
554
1/2
✓ Branch 0 taken 623594 times.
✗ Branch 1 not taken.
 663416 if (type == Dune::GeometryTypes::triangle)
555 {
556 using Corners = Detail::Box::ScvCorners<Dune::GeometryTypes::triangle>;
557 1023888 return Detail::Box::subEntityKeyToCornerStorage<ScvfCornerStorage>(geo_, localFacetIndex, facetCodim, Corners::keys[indexInFacet]);
558 }
559
1/2
✓ Branch 0 taken 143024 times.
✗ Branch 1 not taken.
 151472 else if (type == Dune::GeometryTypes::quadrilateral)
560 {
561 using Corners = Detail::Box::ScvCorners<Dune::GeometryTypes::quadrilateral>;
562 302944 return Detail::Box::subEntityKeyToCornerStorage<ScvfCornerStorage>(geo_, localFacetIndex, facetCodim, Corners::keys[indexInFacet]);
563 }
564 else
565 ✗ DUNE_THROW(Dune::NotImplemented, "Box boundary scvf geometries for dim=" << dim
566 << " dimWorld=" << dimWorld
567 << " type=" << type);
568 }
569
570 //! get scvf normal vector
571 2577654 GlobalPosition normal(const ScvfCornerStorage& p,
572 const std::vector<unsigned int>& scvIndices) const
573 {
574 18043578 auto normal = Dumux::crossProduct(p[1]-p[0], p[2]-p[0]);
575 5155308 normal /= normal.two_norm();
576
577 3304758 const auto v = geo_.corner(scvIndices[1]) - geo_.corner(scvIndices[0]);
578 2577654 const auto s = v*normal;
579
4/4
✓ Branch 0 taken 174312 times.
✓ Branch 1 taken 2338026 times.
✓ Branch 2 taken 174312 times.
✓ Branch 3 taken 2338026 times.
 5155308 if (std::signbit(s))
580 normal *= -1;
581
582 2577654 return normal;
583 }
584
585 //! number of sub control volume faces (number of edges)
586 ✗ std::size_t numInteriorScvf() const
587 {
588
12/12
✓ Branch 2 taken 26472 times.
✓ Branch 3 taken 17648 times.
✓ Branch 4 taken 26472 times.
✓ Branch 5 taken 17648 times.
✓ Branch 7 taken 26472 times.
✓ Branch 8 taken 17648 times.
✓ Branch 9 taken 26472 times.
✓ Branch 10 taken 17648 times.
✓ Branch 15 taken 17648 times.
✓ Branch 16 taken 26472 times.
✓ Branch 17 taken 17648 times.
✓ Branch 18 taken 26472 times.
 218768 return referenceElement(geo_).size(dim-1);
589 }
590
591 //! number of sub control volumes (number of vertices)
592 ✗ std::size_t numScv() const
593 {
594
24/24
✓ Branch 1 taken 17648 times.
✓ Branch 2 taken 4412 times.
✓ Branch 3 taken 17648 times.
✓ Branch 4 taken 4412 times.
✓ Branch 8 taken 88240 times.
✓ Branch 9 taken 22060 times.
✓ Branch 10 taken 88240 times.
✓ Branch 11 taken 22060 times.
✓ Branch 13 taken 17648 times.
✓ Branch 14 taken 4412 times.
✓ Branch 15 taken 17648 times.
✓ Branch 16 taken 4412 times.
✓ Branch 24 taken 22060 times.
✓ Branch 25 taken 4412 times.
✓ Branch 26 taken 22060 times.
✓ Branch 27 taken 4412 times.
✓ Branch 29 taken 4412 times.
✓ Branch 30 taken 17648 times.
✓ Branch 31 taken 4412 times.
✓ Branch 32 taken 17648 times.
✓ Branch 34 taken 17648 times.
✓ Branch 35 taken 4412 times.
✓ Branch 36 taken 17648 times.
✓ Branch 37 taken 4412 times.
 242660 return referenceElement(geo_).size(dim);
595 }
596
597 //! the wrapped element geometry
598 ✗ const typename Element::Geometry& elementGeometry() const
599 ✗ { return geo_; }
600
601 private:
602 const typename Element::Geometry& geo_; //!< Reference to the element geometry
603 };
604
605 } // end namespace Dumux
606
607 #endif
608