GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/io/grid/cakegridmanager.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 195 235 83.0%
Functions: 12 20 60.0%
Branches: 348 797 43.7%

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 InputOutput
10 * \brief Provides a grid manager for a piece of cake grid
11 */
12 #ifndef DUMUX_CAKE_GRID_MANAGER_HH
13 #define DUMUX_CAKE_GRID_MANAGER_HH
14
15 #include <vector>
16 #include <iostream>
17 #include <cmath>
18 #include <algorithm>
19
20 #include <dune/common/dynvector.hh>
21 #include <dune/common/float_cmp.hh>
22 #include <dune/common/math.hh>
23 #include <dune/grid/common/gridfactory.hh>
24 #include <dumux/common/parameters.hh>
25
26 namespace Dumux {
27
28 /*!
29 * \ingroup InputOutput
30 * \brief Provides a grid manager with a method for creating creating vectors
31 * with polar Coordinates and one for creating a Cartesian grid from
32 * these polar coordinates.
33 */
34 template <class Grid>
35
2/8
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 11 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
44 class CakeGridManager
36 {
37 using Scalar = typename Grid::ctype;
38
39 using GridFactory = Dune::GridFactory<Grid>;
40 using GridPointer = std::shared_ptr<Grid>;
41
42 enum { dim = Grid::dimension,
43 dimWorld = Grid::dimensionworld };
44
45 public:
46 /*!
47 * \brief Make the grid.
48 */
49 44 void init(const std::string& modelParamGroup = "")
50 {
51 static_assert(dim == 2 || dim == 3, "The CakeGridManager is only implemented for 2D and 3D.");
52
53 44 const bool verbose = getParamFromGroup<bool>(modelParamGroup, "Grid.Verbosity", false);
54
55 44 std::array<std::vector<Scalar>, dim> polarCoordinates;
56 // Indices specifying in which direction the piece of cake is oriented
57 44 Dune::FieldVector<int, dim> indices(-1);
58
1/2
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
44 createVectors(polarCoordinates, indices, modelParamGroup, verbose);
59
60
2/8
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
88 gridPtr() = createCakeGrid(polarCoordinates, indices, modelParamGroup, verbose);
61 44 loadBalance();
62 44 }
63
64 /*!
65 * \brief Create vectors containing polar coordinates of all points.
66 *
67 * All keys are expected to be in group GridParameterGroup.
68 * The following keys are recognized:
69 * - Radial : min/max value for radial coordinate
70 * - Angular : min/max value for angular coordinate
71 * - Axial : min/max value for axial coordinate
72 * Adding 0, 1 (or 2 in 3D) specifies in which direction (x, y and z, respectively)
73 * the radial, angular and axial direction are oriented
74 * - Cells : number of cells array for x-coordinate (Again, an added 0, 1 or 2 specifies x, y or z)
75 * - Grading : grading factor array for x-coordinate (Same here)
76 * - Verbosity : whether the grid construction should output to standard out
77 *
78 * The grading factor \f$ g \f$ specifies the ratio between the next and the current cell size:
79 * \f$ g = \frac{h_{i+1}}{h_i} \f$.
80 * Negative grading factors are converted to
81 * \f$ g = -\frac{1}{g_\textrm{negative}} \f$
82 * to avoid issues with imprecise fraction numbers.
83 */
84 44 static void createVectors(std::array<std::vector<Scalar>, dim> &polarCoordinates,
85 Dune::FieldVector<int, dim> &indices,
86 const std::string& modelParamGroup,
87 bool verbose = false)
88 {
89 // The positions
90 44 std::array<std::vector<Scalar>, dim> positions;
91
92
2/2
✓ Branch 0 taken 55 times.
✓ Branch 1 taken 22 times.
154 for (int i = 0; i < dim; ++i)
93 {
94
5/14
✓ Branch 1 taken 55 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 55 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 55 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 55 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 55 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
220 const bool hasRadial = hasParamInGroup(modelParamGroup, "Grid.Radial" + std::to_string(i));
95
5/14
✓ Branch 1 taken 55 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 55 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 55 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 55 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 55 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
220 const bool hasAngular = hasParamInGroup(modelParamGroup, "Grid.Angular" + std::to_string(i));
96
5/14
✓ Branch 1 taken 33 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 33 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 33 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 33 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 33 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
176 const bool hasAxial = (dim == 3) && hasParamInGroup(modelParamGroup, "Grid.Axial" + std::to_string(i));
97
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 55 times.
110 if (static_cast<int>(hasRadial) + static_cast<int>(hasAngular) + static_cast<int>(hasAxial) != 1)
98 DUNE_THROW(Dune::RangeError, "Multiple or no position vectors (radial, angular, axial) specified in coord direction: " << i << std::endl);
99
100
2/2
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 33 times.
110 if (hasRadial)
101 {
102
6/16
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 22 times.
✗ Branch 8 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 22 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 22 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 22 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
132 positions[i] = getParamFromGroup<std::vector<Scalar>>(modelParamGroup, "Grid.Radial" + std::to_string(i));
103 88 indices[0] = i; // Index specifying radial direction
104 }
105
2/2
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 11 times.
66 else if (hasAngular)
106 {
107
6/16
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 22 times.
✗ Branch 8 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 22 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 22 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 22 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
132 positions[i] = getParamFromGroup<std::vector<Scalar>>(modelParamGroup, "Grid.Angular" + std::to_string(i));
108 88 indices[1] = i; // Index specifying angular direction
109 }
110 else // hasAxial
111 {
112
6/20
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 11 times.
✗ Branch 8 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 11 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 11 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 11 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
66 positions[i] = getParamFromGroup<std::vector<Scalar>>(modelParamGroup, "Grid.Axial" + std::to_string(i));
113 44 indices[2] = i; // Index specifying axial direction
114 }
115
116
6/12
✓ Branch 0 taken 55 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 55 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 55 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 55 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 55 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 55 times.
660 if (!std::is_sorted(positions[i].begin(), positions[i].end()))
117 DUNE_THROW(Dune::GridError, "Make sure to specify a monotone increasing \"Positions\" array");
118
119
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 55 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 55 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 55 times.
330 if(positions[i].size() < 2)
120 DUNE_THROW(Dune::GridError, "Make sure to specify position arrays with at least two entries (min and max value).");
121 }
122
123 // check that all indices are specified i.e. > 0
124
2/2
✓ Branch 0 taken 55 times.
✓ Branch 1 taken 22 times.
154 for (int i = 0; i < dim; ++i)
125
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 55 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 55 times.
220 if (indices[i] < 0)
126 DUNE_THROW(Dune::RangeError, "Please specify Positions Angular and Radial and Axial correctly and unambiguously!" << std::endl);
127
128 // the number of cells (has a default)
129 44 std::array<std::vector<int>, dim> cells;
130
2/2
✓ Branch 0 taken 55 times.
✓ Branch 1 taken 22 times.
154 for (int i = 0; i < dim; ++i)
131 {
132
4/8
✓ Branch 1 taken 55 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 55 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 55 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 55 times.
✗ Branch 11 not taken.
440 cells[i].resize(positions[i].size()-1, 1.0);
133
7/18
✓ Branch 1 taken 55 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 55 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 55 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 55 times.
✗ Branch 11 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 55 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 55 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 55 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
440 cells[i] = getParamFromGroup<std::vector<int>>(modelParamGroup, "Grid.Cells" + std::to_string(i), cells[i]);
134
5/10
✗ Branch 0 not taken.
✓ Branch 1 taken 55 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 55 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 55 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 55 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 55 times.
550 if (cells[i].size() + 1 != positions[i].size())
135 DUNE_THROW(Dune::RangeError, "Make sure to specify the correct length of \"Cells\" and \"Positions\" arrays");
136 }
137
138 // grading factor (has a default)
139 44 std::array<std::vector<Scalar>, dim> grading;
140
2/2
✓ Branch 0 taken 55 times.
✓ Branch 1 taken 22 times.
154 for (int i = 0; i < dim; ++i)
141 {
142
4/8
✓ Branch 1 taken 55 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 55 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 55 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 55 times.
✗ Branch 11 not taken.
440 grading[i].resize(positions[i].size()-1, 1.0);
143
7/18
✓ Branch 1 taken 55 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 55 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 55 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 55 times.
✗ Branch 11 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 55 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 55 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 55 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
440 grading[i] = getParamFromGroup<std::vector<Scalar>>(modelParamGroup, "Grid.Grading" + std::to_string(i), grading[i]);
144
5/10
✗ Branch 0 not taken.
✓ Branch 1 taken 55 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 55 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 55 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 55 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 55 times.
550 if (grading[i].size() + 1 != positions[i].size())
145 DUNE_THROW(Dune::RangeError, "Make sure to specify the correct length of \"Grading\" and \"Positions\" arrays");
146 }
147
148 // make the grid
149 44 std::array<std::vector<Scalar>, dim> globalPositions;
150 using std::pow;
151 using Dune::power;
152
2/2
✓ Branch 0 taken 55 times.
✓ Branch 1 taken 22 times.
154 for (int dimIdx = 0; dimIdx < dim; dimIdx++)
153 {
154 // Each grid direction is subdivided into (numCells + 1) points
155 std::size_t numGlobalPositions = 1;
156
6/6
✓ Branch 0 taken 99 times.
✓ Branch 1 taken 55 times.
✓ Branch 2 taken 99 times.
✓ Branch 3 taken 55 times.
✓ Branch 4 taken 99 times.
✓ Branch 5 taken 55 times.
924 for (int zoneIdx = 0; zoneIdx < cells[dimIdx].size(); ++zoneIdx)
157 594 numGlobalPositions += cells[dimIdx][zoneIdx];
158
159
2/4
✓ Branch 1 taken 55 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 55 times.
✗ Branch 5 not taken.
220 globalPositions[dimIdx].resize(numGlobalPositions);
160 std::size_t posIdx = 0;
161
6/6
✓ Branch 0 taken 99 times.
✓ Branch 1 taken 55 times.
✓ Branch 2 taken 99 times.
✓ Branch 3 taken 55 times.
✓ Branch 4 taken 99 times.
✓ Branch 5 taken 55 times.
924 for (int zoneIdx = 0; zoneIdx < cells[dimIdx].size(); ++zoneIdx)
162 {
163
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 99 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 99 times.
396 const Scalar lower = positions[dimIdx][zoneIdx];
164
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 99 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 99 times.
396 const Scalar upper = positions[dimIdx][zoneIdx+1];
165
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 99 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 99 times.
396 const int numCells = cells[dimIdx][zoneIdx];
166
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 99 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 99 times.
396 Scalar gradingFactor = grading[dimIdx][zoneIdx];
167 198 const Scalar length = upper - lower;
168 198 Scalar height = 1.0;
169 198 bool increasingCellSize = false;
170
171
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 99 times.
198 if (verbose)
172 {
173 std::cout << "dim " << dimIdx
174 << " lower " << lower
175 << " upper " << upper
176 << " numCells " << numCells
177 << " grading " << gradingFactor;
178 }
179
180
2/2
✓ Branch 0 taken 33 times.
✓ Branch 1 taken 66 times.
198 if (gradingFactor > 1.0)
181 66 increasingCellSize = true;
182
183 // take absolute values and reverse cell size increment to achieve
184 // reverse behavior for negative values
185
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 88 times.
198 if (gradingFactor < 0.0)
186 {
187 using std::abs;
188
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
22 gradingFactor = abs(gradingFactor);
189
190
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
22 if (gradingFactor < 1.0)
191 increasingCellSize = true;
192 }
193
194
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 99 times.
✓ Branch 2 taken 55 times.
✓ Branch 3 taken 44 times.
198 const bool useGrading = Dune::FloatCmp::eq(gradingFactor, 1.0) ? false : true;
195
196 // if the grading factor is exactly 1.0 do equal spacing
197 if (!useGrading)
198 {
199 110 height = 1.0 / numCells;
200
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 55 times.
110 if (verbose)
201 std::cout << " -> h " << height * length << std::endl;
202 }
203 // if grading factor is not 1.0, do power law spacing
204 else
205 {
206 88 height = (1.0 - gradingFactor) / (1.0 - power(gradingFactor, numCells));
207
208
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 44 times.
88 if (verbose)
209 {
210 std::cout << " -> grading_eff " << gradingFactor
211 << " h_min " << height * power(gradingFactor, 0) * length
212 << " h_max " << height * power(gradingFactor, numCells-1) * length
213 << std::endl;
214 }
215 }
216
217 // the positions inside a global segment
218
3/8
✓ Branch 1 taken 99 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 99 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 99 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
792 Dune::DynamicVector<Scalar> localPositions(numCells, 0.0);
219
2/2
✓ Branch 0 taken 286 times.
✓ Branch 1 taken 99 times.
770 for (int i = 1; i < numCells; i++)
220 {
221 572 Scalar hI = height;
222
2/2
✓ Branch 0 taken 88 times.
✓ Branch 1 taken 198 times.
572 if (useGrading)
223 {
224
2/2
✓ Branch 0 taken 66 times.
✓ Branch 1 taken 22 times.
176 if (increasingCellSize)
225 264 hI *= power(gradingFactor, i-1);
226
227 else
228 88 hI *= power(gradingFactor, numCells-i);
229 }
230 1716 localPositions[i] = localPositions[i-1] + hI;
231 }
232
233 localPositions *= length;
234 localPositions += lower;
235
236
2/2
✓ Branch 0 taken 385 times.
✓ Branch 1 taken 99 times.
968 for (int i = 0; i < numCells; ++i)
237 3080 globalPositions[dimIdx][posIdx++] = localPositions[i];
238 }
239
240 550 globalPositions[dimIdx][posIdx] = positions[dimIdx].back();
241 }
242
243
4/8
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 22 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 22 times.
✗ Branch 11 not taken.
176 polarCoordinates[0] = globalPositions[indices[0]];
244
4/8
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 22 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 22 times.
✗ Branch 11 not taken.
176 polarCoordinates[1] = globalPositions[indices[1]];
245 if (dim == 3)
246
5/10
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 11 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 11 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 11 times.
✗ Branch 14 not taken.
110 polarCoordinates[2] = globalPositions[dim - indices[0] - indices[1]];
247
248 // convert angular coordinates into radians
249 220 std::transform(polarCoordinates[1].begin(), polarCoordinates[1].end(),
250 88 polarCoordinates[1].begin(),
251 396 [](Scalar s){ return s*M_PI/180; });
252 44 }
253
254 /*!
255 * \brief Creates Cartesian grid from polar coordinates.
256 *
257 * \param polarCoordinates Vector containing radial, angular and axial coordinates (in this order)
258 * \param indices Indices specifying the radial, angular and axial direction (in this order)
259 * \param modelParamGroup name of the model parameter group
260 * \param verbose if the output should be verbose
261 */
262 44 std::unique_ptr<Grid> createCakeGrid(std::array<std::vector<Scalar>, dim> &polarCoordinates,
263 Dune::FieldVector<int, dim> &indices,
264 const std::string& modelParamGroup,
265 bool verbose = false)
266 {
267
1/2
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
108 GridFactory gridFactory;
268
269 // create grid on rank 0, return empty grid on all other ranks
270
7/7
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 18 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 12 times.
✓ Branch 4 taken 10 times.
✓ Branch 5 taken 12 times.
✓ Branch 6 taken 6 times.
92 if (gridFactory.comm().rank() != 0)
271
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
16 return std::unique_ptr<Grid>(gridFactory.createGrid());
272
273
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 10 times.
28 const auto& dR = polarCoordinates[0];
274
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 10 times.
28 const auto& dA = polarCoordinates[1];
275
276 // For the special case of 360°, the last element is connected with the first one.
277 // Thus, one needs to distinguish here between all other cases, were the normal procedure
278 // is applied for all elements until the last one (= dA.size() - 1) and the 360° case,
279 // where it is only applied until the second to last element (dA.size() - 2).
280
8/8
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 10 times.
✓ Branch 4 taken 4 times.
✓ Branch 5 taken 10 times.
✓ Branch 6 taken 10 times.
✓ Branch 7 taken 4 times.
92 const int dAMax = Dune::FloatCmp::eq(dA[dA.size()-1], 2*M_PI) ? dA.size() - 2 : dA.size() - 1;
281 28 constexpr auto type = Dune::GeometryTypes::cube(dim);
282
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
28 const bool hasHole = dR[0] >= 1.0e-8*dR.back();
283
284 // create nodes
285 if constexpr (dim == 3)
286 {
287 14 constexpr auto prismType = Dune::GeometryTypes::prism;
288
3/6
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 7 times.
✗ Branch 7 not taken.
42 std::vector<Scalar> dZ = polarCoordinates[2];
289
2/2
✓ Branch 0 taken 58 times.
✓ Branch 1 taken 7 times.
130 for (int j = 0; j <= dAMax; ++j)
290 {
291
4/4
✓ Branch 0 taken 696 times.
✓ Branch 1 taken 58 times.
✓ Branch 2 taken 696 times.
✓ Branch 3 taken 58 times.
1624 for (int l = 0; l <= dZ.size() - 1; ++l)
292 {
293
6/6
✓ Branch 0 taken 96 times.
✓ Branch 1 taken 600 times.
✓ Branch 2 taken 3384 times.
✓ Branch 3 taken 696 times.
✓ Branch 4 taken 3384 times.
✓ Branch 5 taken 696 times.
13944 for (int i = hasHole ? 0 : 1; i <= dR.size()- 1; ++i)
294 {
295 // transform into Cartesian coordinates
296 using std::cos;
297 using std::sin;
298 6768 Dune::FieldVector <double, dim> v(0.0);
299
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 3384 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3384 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 3384 times.
20304 v[indices[2]] = dZ[l];
300
4/8
✗ Branch 0 not taken.
✓ Branch 1 taken 3384 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3384 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 3384 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 3384 times.
27072 v[indices[0]] = cos(dA[j])*dR[i];
301
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 3384 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3384 times.
13536 v[indices[1]] = sin(dA[j])*dR[i];
302
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3384 times.
6768 if (verbose)
303 printCoordinate(v);
304
1/2
✓ Branch 1 taken 3384 times.
✗ Branch 2 not taken.
6768 gridFactory.insertVertex(v);
305 }
306 }
307 }
308
309 // if the well is not considered, we have to add the points lying on the center-line
310 // this has to be done separately, otherwise these points would occur multiple times
311
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 6 times.
14 if(!hasHole)
312 {
313
4/4
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 1 times.
28 for (int l = 0; l <= dZ.size() - 1; ++l)
314 {
315 24 Dune::FieldVector <double, dim> v(0.0);
316
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 12 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 12 times.
72 v[indices[2]] = dZ[l];
317
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 12 times.
48 v[indices[0]] = 0.0;
318
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 12 times.
48 v[indices[1]] = 0.0;
319
320
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
24 if (verbose)
321 printCoordinate(v);
322
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
24 gridFactory.insertVertex(v);
323 }
324 }
325
326
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
14 if (verbose)
327 std::cout << "Filled node vector" << std::endl;
328
329 // assign nodes
330 14 unsigned int z = 0;
331 14 unsigned int t = 0;
332
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 1 times.
14 unsigned int rSize = hasHole ? dR.size() : dR.size()-1;
333 14 unsigned int zSize = dZ.size();
334
4/4
✓ Branch 0 taken 56 times.
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 56 times.
✓ Branch 3 taken 7 times.
238 for (int j = 0; j < dA.size() - 1; ++j)
335 {
336
4/4
✓ Branch 0 taken 616 times.
✓ Branch 1 taken 56 times.
✓ Branch 2 taken 616 times.
✓ Branch 3 taken 56 times.
1456 for (int l = 0; l < dZ.size() - 1; ++l)
337 {
338
2/2
✓ Branch 0 taken 561 times.
✓ Branch 1 taken 55 times.
1232 if (j < dAMax)
339 {
340
4/4
✓ Branch 0 taken 2244 times.
✓ Branch 1 taken 561 times.
✓ Branch 2 taken 2244 times.
✓ Branch 3 taken 561 times.
6732 for (int i = 0; i < dR.size() - 1; ++i)
341 {
342
2/2
✓ Branch 0 taken 77 times.
✓ Branch 1 taken 2167 times.
4488 if(!hasHole && i==0)
343 {
344
3/8
✓ Branch 1 taken 77 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 77 times.
✓ Branch 5 taken 77 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
462 std::vector<unsigned int> vid({rSize*zSize*(dAMax+1) + l, z, z+rSize*zSize,
345
1/2
✓ Branch 1 taken 77 times.
✗ Branch 2 not taken.
154 rSize*zSize*(dAMax+1) + l+1, z+rSize, z+rSize*zSize+rSize});
346
347
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 77 times.
154 if (verbose)
348 printIndices(vid);
349
350
1/2
✓ Branch 1 taken 77 times.
✗ Branch 2 not taken.
154 gridFactory.insertElement(prismType, vid);
351 }
352 else
353 {
354
2/6
✓ Branch 1 taken 2167 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2167 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
8668 std::vector<unsigned int> vid({z, z+1,
355 8668 z+rSize*zSize, z+rSize*zSize+1,
356 8668 z+rSize, z+rSize+1,
357
1/2
✓ Branch 1 taken 2167 times.
✗ Branch 2 not taken.
4334 z+rSize*zSize+rSize, z+rSize*zSize+rSize+1});
358
359
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2167 times.
4334 if (verbose)
360 printIndices(vid);
361
362
1/2
✓ Branch 1 taken 2167 times.
✗ Branch 2 not taken.
4334 gridFactory.insertElement(type, vid);
363
364
1/2
✓ Branch 0 taken 2167 times.
✗ Branch 1 not taken.
4334 z++;
365 }
366 }
367 1122 z++;
368 }
369 else
370 {
371 // assign nodes for 360°-cake
372
4/4
✓ Branch 0 taken 220 times.
✓ Branch 1 taken 55 times.
✓ Branch 2 taken 220 times.
✓ Branch 3 taken 55 times.
660 for (int i = 0; i < dR.size() - 1; ++i)
373 {
374
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 209 times.
440 if(!hasHole && i==0)
375 {
376
3/8
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 11 times.
✓ Branch 5 taken 11 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
66 std::vector<unsigned int> vid({rSize*zSize*(dAMax+1) + l, z, t,
377
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
22 rSize*zSize*(dAMax+1) + l+1, z+rSize, t+rSize});
378
379
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
22 if (verbose)
380 printIndices(vid);
381
382
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
22 gridFactory.insertElement(prismType, vid);
383 }
384 else
385 {
386
2/8
✓ Branch 1 taken 209 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 209 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
836 std::vector<unsigned int> vid({z, z+1,
387 418 t, t+1,
388 836 z+rSize, z+rSize+1,
389
1/2
✓ Branch 1 taken 209 times.
✗ Branch 2 not taken.
418 t+rSize, t+rSize+1});
390
391
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 209 times.
418 if (verbose)
392 printIndices(vid);
393
394
1/2
✓ Branch 1 taken 209 times.
✗ Branch 2 not taken.
418 gridFactory.insertElement(type, vid);
395 418 t++;
396
1/2
✓ Branch 0 taken 209 times.
✗ Branch 1 not taken.
418 z++;
397 }
398 }
399 110 t++;
400 110 z++;
401
402
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 55 times.
110 if (verbose)
403 std::cout << "assign nodes 360° ends..." << std::endl;
404 }
405 }
406 112 z += rSize;
407 }
408 }
409
410 // for dim = 2
411 else
412 {
413 14 constexpr auto triangleType = Dune::GeometryTypes::simplex(dim);
414
2/2
✓ Branch 0 taken 58 times.
✓ Branch 1 taken 7 times.
130 for (int j = 0; j <= dAMax; ++j)
415 {
416
6/6
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 50 times.
✓ Branch 2 taken 282 times.
✓ Branch 3 taken 58 times.
✓ Branch 4 taken 282 times.
✓ Branch 5 taken 58 times.
1162 for (int i = hasHole ? 0 : 1; i <= dR.size()- 1; ++i)
417 {
418 // transform into Cartesian coordinates
419 564 Dune::FieldVector <double, dim> v(0.0);
420
421
4/8
✗ Branch 0 not taken.
✓ Branch 1 taken 282 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 282 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 282 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 282 times.
2256 v[indices[0]] = cos(dA[j])*dR[i];
422
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 282 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 282 times.
1128 v[indices[1]] = sin(dA[j])*dR[i];
423
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 282 times.
564 if(verbose)
424 printCoordinate(v);
425
1/2
✓ Branch 1 taken 282 times.
✗ Branch 2 not taken.
564 gridFactory.insertVertex(v);
426 }
427 }
428
429 // if the well is not considered, we have to add the point located at the center
430 // this has to be done separately, otherwise this point would occur multiple times
431
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 6 times.
14 if(!hasHole)
432 {
433 2 Dune::FieldVector <double, dim> v(0.0);
434
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
4 v[indices[0]] = 0.0;
435
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
4 v[indices[1]] = 0.0;
436
437
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
2 if(verbose)
438 printCoordinate(v);
439
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
2 gridFactory.insertVertex(v);
440 }
441
442
1/2
✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
28 std::cout << "Filled node vector" << std::endl;
443
444 // assign nodes
445 14 unsigned int z = 0;
446 14 unsigned int t = 0;
447
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 1 times.
14 unsigned int rSize = hasHole ? dR.size() : dR.size()-1;
448
4/4
✓ Branch 0 taken 56 times.
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 56 times.
✓ Branch 3 taken 7 times.
126 for (int j = 0; j < dA.size() - 1; ++j)
449 {
450
2/2
✓ Branch 0 taken 51 times.
✓ Branch 1 taken 5 times.
112 if (j < dAMax)
451 {
452
4/4
✓ Branch 0 taken 204 times.
✓ Branch 1 taken 51 times.
✓ Branch 2 taken 204 times.
✓ Branch 3 taken 51 times.
612 for (int i = 0; i < dR.size() - 1; ++i)
453 {
454
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 197 times.
408 if(!hasHole && i==0)
455 {
456
4/10
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 7 times.
✓ Branch 8 taken 7 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
56 std::vector<unsigned int> vid({rSize*(dAMax+1), z, z+rSize});
457
458
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
14 if (verbose)
459 printIndices(vid);
460
461
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
14 gridFactory.insertElement(triangleType, vid);
462 }
463 else
464 {
465
3/8
✓ Branch 1 taken 197 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 197 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 197 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
1182 std::vector<unsigned int> vid({z, z+1, z+rSize, z+rSize+1});
466
467
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 197 times.
394 if (verbose)
468 printIndices(vid);
469
470
1/2
✓ Branch 1 taken 197 times.
✗ Branch 2 not taken.
394 gridFactory.insertElement(type, vid);
471
1/2
✓ Branch 0 taken 197 times.
✗ Branch 1 not taken.
394 z++;
472 }
473 }
474 102 z++;
475 }
476 else
477 {
478 // assign nodes for 360°-cake
479
4/4
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 20 times.
✓ Branch 3 taken 5 times.
60 for (int i = 0; i < dR.size() - 1; ++i)
480 {
481
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 19 times.
40 if(!hasHole && i==0)
482 {
483
4/10
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
8 std::vector<unsigned int> vid({rSize*(dAMax+1), z, t});
484
485
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
2 if (verbose)
486 printIndices(vid);
487
488
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
2 gridFactory.insertElement(triangleType, vid);
489 }
490 else
491 {
492
3/8
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 19 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 19 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
114 std::vector<unsigned int> vid({z, z+1, t, t+1});
493
494
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 19 times.
38 if (verbose)
495 printIndices(vid);
496
497
1/2
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
38 gridFactory.insertElement(type, vid);
498 38 t++;
499
1/2
✓ Branch 0 taken 19 times.
✗ Branch 1 not taken.
38 z++;
500 }
501 }
502 10 t++;
503 10 z++;
504
505
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
10 if (verbose)
506 std::cout << "assign nodes 360 ends..." << std::endl;
507 }
508 }
509 }
510
511 // create the grid
512
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
28 return std::unique_ptr<Grid>(gridFactory.createGrid());
513 }
514
515 /*!
516 * \brief Returns a reference to the grid.
517 */
518 Grid& grid()
519 {
520
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6 times.
22 return *gridPtr();
521 }
522
523 /*!
524 * \brief Distributes the grid on all processes of a parallel
525 * computation.
526 */
527 void loadBalance()
528 {
529
2/4
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11 times.
✗ Branch 5 not taken.
22 gridPtr()->loadBalance();
530 }
531
532 protected:
533 static void printCoordinate(const Dune::FieldVector <double, dim>& v)
534 { std::cout << "Coordinates of: " << v << std::endl; }
535
536 static void printIndices(const std::vector<unsigned int>& vid)
537 {
538 std::cout << "element vertex indices: ";
539 for (int k = 0; k < vid.size(); ++k)
540 std::cout << vid[k] << " ";
541 std::cout << std::endl;
542 }
543
544 /*!
545 * \brief Returns a reference to the shared pointer to the grid.
546 */
547 GridPointer& gridPtr()
548 {
549
6/12
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 11 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 11 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 6 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 6 times.
66 return cakeGrid_;
550 }
551
552 private:
553 GridPointer cakeGrid_;
554 };
555
556 } // end namespace Dumux
557
558 #endif
559