GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/io/grid/porenetwork/structuredlatticegridcreator.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 172 221 77.8%
Functions: 33 61 54.1%
Branches: 236 613 38.5%
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 PoreNetworkModels
10 * \brief Creates a network grid from a structured lattice. Connections can be randomly deleted.
11 */
12 #ifndef DUMUX_IO_STRUCTURED_LATTICE_GRID_CREATOR_HH
13 #define DUMUX_IO_STRUCTURED_LATTICE_GRID_CREATOR_HH
14
15 #if HAVE_DUNE_FOAMGRID
16
17 #include <vector>
18 #include <memory>
19 #include <type_traits>
20 #include <random>
21
22 #include <dune/common/concept.hh>
23 #include <dune/common/exceptions.hh>
24 #include <dune/grid/common/gridfactory.hh>
25 #include <dune/grid/yaspgrid.hh>
26 #include <dune/geometry/referenceelements.hh>
27
28 // FoamGrid specific includes
29 #include <dune/foamgrid/foamgrid.hh>
30 #include <dune/foamgrid/dgffoam.hh>
31
32 #include <dumux/io/grid/gridmanager_yasp.hh>
33 #include <dumux/common/parameters.hh>
34 #include <dumux/common/exceptions.hh>
35
36 namespace Dumux::PoreNetwork {
37
38 namespace Concept {
39
40 /*!
41 * \ingroup PoreNetworkModels
42 * \brief The element selector concept
43 */
44 template<class GlobalPosition>
45 struct LowDimElementSelector
46 {
47 template<class F>
48 auto require(F&& f) -> decltype(
49 bool(f(std::declval<const GlobalPosition&>(), std::declval<const GlobalPosition&>()))
50 );
51 };
52 } // end namespace Concept
53
54 /*!
55 * \ingroup PoreNetworkModels
56 * \brief Creates a network grid from a structured lattice. Connections can be randomly deleted.
57 */
58 template<int dimWorld>
59 class StructuredLatticeGridCreator
60 {
61 using GridType = Dune::FoamGrid<1, dimWorld>;
62 using Element = typename GridType::template Codim<0>::Entity;
63 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
64 using CoordScalar = typename GridType::ctype;
65 using HostGrid = Dune::YaspGrid<dimWorld, Dune::TensorProductCoordinates<CoordScalar, dimWorld>>;
66 using HostGridView = typename HostGrid::LeafGridView;
67 using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dimWorld>;
68 // grid factory expects unsigned int, will yield compiler warning for std::size_t
69 using VertexIndexPair = std::pair<unsigned int, unsigned int>;
70
71 struct Diagonal
72 {
73 VertexIndexPair localVertexIndices;
74 unsigned int directionNumber;
75 };
76 public:
77
78 using Grid = GridType;
79
80 void init(const std::string& paramGroup = "")
81 {
82 auto useAllElements = [](const GlobalPosition& a, const GlobalPosition& b){ return true; };
83
3/6
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
 18 init(useAllElements, paramGroup);
84 }
85
86 template<class LowDimElementSelector, // cppcheck-suppress syntaxError
87 typename std::enable_if_t<Dune::models<Concept::LowDimElementSelector<GlobalPosition>, LowDimElementSelector>(), int> = 0>
88 28 void init(const LowDimElementSelector& lowDimElementSelector,
89 const std::string& paramGroup = "")
90 {
91 28 paramGroup_ = paramGroup;
92
4/10
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 18 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 18 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 18 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 56 removeThroatsOnBoundary_ = getParamFromGroup<std::vector<std::size_t>>(paramGroup, "Grid.RemoveThroatsOnBoundary",
93 std::vector<std::size_t>());
94
95 28 setElementSelector_(lowDimElementSelector);
96 28 initRandomNumberGenerator_();
97
98 // create the host grid
99 28 const auto hostGridInputData = getHostGridInputData_();
100 using HostGrid = Dune::YaspGrid<dimWorld, Dune::TensorProductCoordinates<CoordScalar, dimWorld>>;
101 using HastGridManager = GridManager<HostGrid>;
102
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 56 HastGridManager hostGridManager;
103
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 28 hostGridManager.init(hostGridInputData.positions, hostGridInputData.cells, hostGridInputData.grading, paramGroup_);
104
5/10
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 18 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 18 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 18 times.
✓ Branch 12 taken 18 times.
✗ Branch 13 not taken.
 28 hostGridView_ = std::make_unique<HostGridView>(hostGridManager.grid().leafGridView());
105
106 28 hostGridLowerLeft_ = hostGridInputData.lowerLeft;
107 28 hostGridUpperRight_ = hostGridInputData.upperRight;
108
109
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 28 convertHostGridToLowDimGrid_();
110 28 }
111
112 /*!
113 * \brief Call loadBalance() function of the grid.
114 */
115 void loadBalance()
116 {
117 if (Dune::MPIHelper::getCommunication().size() > 1)
118 gridPtr_->loadBalance();
119 }
120
121 /*!
122 * \brief Returns a reference to the grid.
123 */
124 Grid& grid()
125 { return *gridPtr(); }
126
127 /*!
128 * \brief Returns a reference to the grid pointer (std::shared_ptr<Grid>)
129 */
130 std::shared_ptr<Grid>& gridPtr()
131 { return gridPtr_; }
132
133 private:
134
135 template<class LowDimElementSelector>
136 28 void setElementSelector_(const LowDimElementSelector& lowDimElementSelector)
137 {
138
4/4
✓ Branch 0 taken 13 times.
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 13 times.
✓ Branch 3 taken 5 times.
 56 if (!removeThroatsOnBoundary_.empty())
139 {
140 37384 auto finalSelector = [&](const GlobalPosition& a, const GlobalPosition& b)
141 {
142
9/12
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 36819 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 232 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 7 times.
✓ Branch 11 taken 285 times.
✓ Branch 13 taken 7 times.
✗ Branch 14 not taken.
 37349 static const auto lambdaToRemoveThroatsOnBoundary = getLambdaToRemoveThroatsOnBoundary_();
143 using std::min;
144 37349 return min(lambdaToRemoveThroatsOnBoundary(a, b), lowDimElementSelector(a, b));
145 };
146 22 considerLowDimElement_ = finalSelector;
147 }
148 else
149 6 considerLowDimElement_ = lowDimElementSelector;
150 28 }
151
152 28 void initRandomNumberGenerator_()
153 {
154 28 directionProbability_ = getDirectionProbability_();
155
156
8/14
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 18 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 18 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 18 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 7 times.
✓ Branch 12 taken 11 times.
✓ Branch 13 taken 7 times.
✓ Branch 14 taken 11 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
 140 if (hasParamInGroup(paramGroup_, "Grid.DeletionRandomNumberSeed"))
157 {
158 11 const auto seed = getParamFromGroup<std::size_t>(paramGroup_, "Grid.DeletionRandomNumberSeed");
159 11 generator_.seed(seed);
160 }
161 else
162 {
163
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 9 not taken.
✓ Branch 10 taken 11 times.
✓ Branch 12 taken 11 times.
✗ Branch 13 not taken.
 51 std::random_device rd;
164
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
 34 generator_.seed(rd());
165 }
166 28 }
167
168 28 void convertHostGridToLowDimGrid_()
169 {
170 28 resizeContainers_();
171
172
1/2
✓ Branch 2 taken 15535 times.
✗ Branch 3 not taken.
 31602 for (const auto& element : elements(*hostGridView_))
173 {
174 15773 const auto geometry = element.geometry();
175 31546 const auto refElement = ReferenceElements::general(geometry.type());
176 15773 treatEdges_(element, refElement, geometry);
177 16033 treatDiagonals_(element, refElement, geometry);
178 }
179
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 18 times.
 56 if (elementSet_.empty())
180 ✗ DUNE_THROW(Dune::GridError, "Trying to create pore network with zero elements!");
181
182 // make the actual network grid
183 28 Dune::GridFactory<Grid> factory;
184
4/4
✓ Branch 0 taken 17063 times.
✓ Branch 1 taken 18 times.
✓ Branch 2 taken 17063 times.
✓ Branch 3 taken 18 times.
 17479 for (auto&& vertex : vertexSet_)
185
1/2
✓ Branch 1 taken 17063 times.
✗ Branch 2 not taken.
 17395 factory.insertVertex(vertex);
186
4/4
✓ Branch 0 taken 43907 times.
✓ Branch 1 taken 18 times.
✓ Branch 2 taken 43907 times.
✓ Branch 3 taken 18 times.
 88730 for (auto&& element : elementSet_)
187
4/10
✓ Branch 1 taken 43907 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43907 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 43907 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 43907 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 177292 factory.insertElement(Dune::GeometryTypes::cube(1), {element.first, element.second});
188
189
5/12
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 18 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 18 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 18 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 18 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 28 gridPtr_ = std::shared_ptr<Grid>(factory.createGrid());
190 28 }
191
192 28 void resizeContainers_()
193 {
194 56 vertexInserted_.resize(hostGridView_->size(dimWorld), false);
195 56 hostVertexIdxToVertexIdx_.resize(hostGridView_->size(dimWorld));
196 56 edgeInserted_.resize(hostGridView_->size(dimWorld-1), false);
197 56 faceInserted_.resize(hostGridView_->size(dimWorld-2), false);
198 28 }
199
200 template<class HostGridElement>
201 214784 bool maybeInsertElementAndVertices_(const HostGridElement& hostGridElement,
202 const typename HostGridElement::Geometry& hostGridElementGeometry,
203 std::size_t localVertexIdx0,
204 std::size_t localVertexIdx1,
205 double directionProbability)
206 {
207
2/2
✓ Branch 1 taken 43907 times.
✓ Branch 2 taken 168229 times.
 214784 if (neglectElement_(hostGridElementGeometry, localVertexIdx0, localVertexIdx1, directionProbability))
208 return false;
209 else
210 {
211 44323 insertElementAndVertices_(hostGridElement, hostGridElementGeometry, localVertexIdx0, localVertexIdx1, directionProbability);
212 44323 return true;
213 }
214 }
215
216 template<class HostGridElement>
217 44323 void insertElementAndVertices_(const HostGridElement& hostGridElement,
218 const typename HostGridElement::Geometry& hostGridElementGeometry,
219 std::size_t localVertexIdx0,
220 std::size_t localVertexIdx1,
221 double directionProbability)
222 {
223 // get global vertex indices w.r.t. host grid
224
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 43907 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 43907 times.
 88646 const auto vIdxGlobal0 = hostGridView_->indexSet().subIndex(hostGridElement, localVertexIdx0, dimWorld);
225
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 43907 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 43907 times.
 88646 const auto vIdxGlobal1 = hostGridView_->indexSet().subIndex(hostGridElement, localVertexIdx1, dimWorld);
226
227 219951 auto getGobalVertexIdx = [&](auto globalIdx, auto localIdx)
228 {
229
18/18
✓ Branch 0 taken 16731 times.
✓ Branch 1 taken 70251 times.
✓ Branch 2 taken 16731 times.
✓ Branch 3 taken 70251 times.
✓ Branch 4 taken 16731 times.
✓ Branch 5 taken 70251 times.
✓ Branch 6 taken 138 times.
✓ Branch 7 taken 138 times.
✓ Branch 8 taken 138 times.
✓ Branch 9 taken 138 times.
✓ Branch 10 taken 138 times.
✓ Branch 11 taken 138 times.
✓ Branch 12 taken 194 times.
✓ Branch 13 taken 362 times.
✓ Branch 14 taken 194 times.
✓ Branch 15 taken 362 times.
✓ Branch 16 taken 194 times.
✓ Branch 17 taken 362 times.
 263442 if (!vertexInserted_[globalIdx])
230 {
231
0/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 34126 vertexInserted_[globalIdx] = true;
232
0/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 34126 hostVertexIdxToVertexIdx_[globalIdx] = vertexSet_.size();
233
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 34126 vertexSet_.push_back(hostGridElementGeometry.corner(localIdx));
234 }
235
236 175628 return hostVertexIdxToVertexIdx_[globalIdx];
237 };
238
239 44323 const auto newVertexIdx0 = getGobalVertexIdx(vIdxGlobal0, localVertexIdx0);
240 44323 const auto newVertexIdx1 = getGobalVertexIdx(vIdxGlobal1, localVertexIdx1);
241 44323 elementSet_.emplace_back(newVertexIdx0, newVertexIdx1);
242 44323 }
243
244 28 auto getDirectionProbability_() const
245 {
246 std::array<double, numDirections_()> directionProbability;
247 28 directionProbability.fill(0.0); // do not delete any throats
248
249 // convenience option to delete all diagonal throats
250
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 18 times.
 28 if (getParamFromGroup<bool>(paramGroup_, "Grid.RegularLattice", false))
251 {
252 directionProbability.fill(1.0); // delete all throats ...
253 ✗ for (int i = 0; i < dimWorld; ++i)
254 ✗ directionProbability[i] = 0.0; // ... but not the ones parallel to the main axes
255
256 return directionProbability;
257 }
258
259 // get user input or print out help message explaining correct usage
260
8/14
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 18 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 18 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 18 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 14 times.
✓ Branch 12 taken 4 times.
✓ Branch 13 taken 14 times.
✓ Branch 14 taken 4 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
 140 if (hasParamInGroup(paramGroup_, "Grid.DeletionProbability"))
261 {
262 try
263 {
264
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 24 directionProbability = getParamFromGroup<decltype(directionProbability)>(paramGroup_, "Grid.DeletionProbability");
265 }
266 ✗ catch(Dumux::ParameterException &e)
267 {
268 ✗ throwDirectionError_();
269 }
270 }
271
272 return directionProbability;
273 }
274
275 ✗ void throwDirectionError_() const
276 {
277 // define directions (to be used for user specified anisotropy)
278 ✗ Dune::FieldVector<std::string, numDirections_()> directions;
279 if (dimWorld < 3) // 2D
280 {
281 // x, y
282 ✗ directions[0] = "1: (1, 0)\n";
283 ✗ directions[1] = "2: (0, 1)\n";
284 // diagonals through cell midpoint
285 ✗ directions[2] = "3: (1, 1)\n";
286 ✗ directions[3] = "4: (1, -1)\n";
287 }
288 else // 3D
289 {
290 // x, y, z
291 ✗ directions[0] = " 1: (1, 0, 0)\n";
292 ✗ directions[1] = "2: (0, 1, 0)\n";
293 ✗ directions[2] = "3: (0, 0, 1)\n";
294 //face diagonals
295 ✗ directions[3] = "4: (1, 1, 0)\n";
296 ✗ directions[4] = "5: (1, -1, 0)\n";
297 ✗ directions[5] = "6: (1, 0, 1)\n";
298 ✗ directions[6] = "7: (1, 0, -1)\n";
299 ✗ directions[7] = "8: (0, 1, 1)\n";
300 ✗ directions[8] = "9: (0, 1, -1)\n";
301 // diagonals through cell midpoint
302 ✗ directions[9] = "10: (1, 1, 1)\n";
303 ✗ directions[10] = "11: (1, 1, -1)\n";
304 ✗ directions[11] = "12: (-1, 1, 1)\n";
305 ✗ directions[12] = "13: (-1, -1, 1)\n";
306 }
307 ✗ DUNE_THROW(ParameterException, "You must specify probabilities for all directions (" << numDirections_() << ") \n" << directions << "\nExample (3D):\n\n"
308 << "DeletionProbability = 0.5 0.5 0 0 0 0 0 0 0 0 0 0 0 \n\n"
309 << "deletes approximately 50% of all throats in x and y direction, while no deletion in any other direction takes place.\n" );
310 }
311
312 static constexpr std::size_t numDirections_()
313 { return (dimWorld < 3) ? 4 : 13; }
314
315 template<class Geometry>
316 214784 bool neglectElement_(Geometry& hostGridElementGeometry,
317 std::size_t localVertexIdx0,
318 std::size_t localVertexIdx1,
319 double directionProbability)
320 {
321 429568 if (randomNumer_() < directionProbability)
322 return true;
323
324 // TODO Change order of execution: check considerLowDimElement_ before checking the random number
325 // TODO This change will alter the reference solution because randomNumer_() gets called less, therefore the random numbers are different
326
3/8
✗ Branch 0 not taken.
✓ Branch 1 taken 49269 times.
✓ Branch 2 taken 5362 times.
✓ Branch 3 taken 43907 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 199340 if (!considerLowDimElement_(hostGridElementGeometry.corner(localVertexIdx0), hostGridElementGeometry.corner(localVertexIdx1)))
327 5512 return true;
328
329 return false;
330 }
331
332 auto randomNumer_()
333
6/6
✓ Branch 1 taken 48703 times.
✓ Branch 2 taken 160785 times.
✓ Branch 4 taken 234 times.
✓ Branch 5 taken 1838 times.
✓ Branch 7 taken 332 times.
✓ Branch 8 taken 244 times.
 212136 { return uniformdist_(generator_); }
334
335 28 auto getHostGridInputData_() const
336 {
337 ✗ struct InputData
338 {
339 std::array<std::vector<int>, dimWorld> cells;
340 std::array<std::vector<CoordScalar>, dimWorld> positions;
341 std::array<std::vector<CoordScalar>, dimWorld> grading;
342 GlobalPosition lowerLeft;
343 GlobalPosition upperRight;
344 28 } inputData;
345
346 // convenience references
347 28 auto& cells = inputData.cells;
348 28 auto& positions = inputData.positions;
349 28 auto& grading = inputData.grading;
350 28 auto& lowerLeft = inputData.lowerLeft;
351 28 auto& upperRight = inputData.upperRight;
352
353 // try to get the pore positions explicitly ...
354
2/2
✓ Branch 0 taken 45 times.
✓ Branch 1 taken 18 times.
 95 for (int i = 0; i < dimWorld; ++i)
355
8/22
✓ Branch 1 taken 45 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 45 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 45 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 45 times.
✗ Branch 11 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 45 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 45 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 45 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 45 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
 268 positions[i] = getParamFromGroup<std::vector<CoordScalar>>(paramGroup_, "Grid.Positions" + std::to_string(i), std::vector<CoordScalar>{});
356
16/84
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 10 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 10 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 7 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 7 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 1 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✓ Branch 29 taken 8 times.
✗ Branch 30 not taken.
✓ Branch 31 taken 8 times.
✗ Branch 32 not taken.
✓ Branch 33 taken 8 times.
✗ Branch 34 not taken.
✓ Branch 35 taken 8 times.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✓ Branch 45 taken 2 times.
✗ Branch 46 not taken.
✓ Branch 47 taken 2 times.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✗ Branch 69 not taken.
✗ Branch 70 not taken.
✗ Branch 71 not taken.
✗ Branch 72 not taken.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✗ Branch 75 not taken.
✗ Branch 76 not taken.
✓ Branch 77 taken 8 times.
✗ Branch 78 not taken.
✓ Branch 79 taken 8 times.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✗ Branch 83 not taken.
 148 if (std::none_of(positions.begin(), positions.end(), [&](auto& v){ return v.empty(); }))
357 {
358 ✗ for (int i = 0; i < dimWorld; ++i)
359 {
360 ✗ cells[i].resize(positions[i].size()-1, 1.0);
361 ✗ grading[i].resize(positions[i].size()-1, 1.0);
362 }
363 }
364 else // .. or calculate positions from input data
365 {
366
1/4
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 56 const auto lowerLeft = getParamFromGroup<GlobalPosition>(paramGroup_, "Grid.LowerLeft", GlobalPosition(0.0));
367
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 28 const auto upperRight = getParamFromGroup<GlobalPosition>(paramGroup_, "Grid.UpperRight");
368
2/4
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 18 times.
✗ Branch 4 not taken.
 84 const auto numPores = getParamFromGroup<std::vector<unsigned int>>(paramGroup_, "Grid.NumPores");
369
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 18 times.
 56 if (numPores.size() != dimWorld)
370 ✗ DUNE_THROW(ParameterException, "Grid.NumPores has to be a space-separated list of " << dimWorld << " integers!");
371
372
2/2
✓ Branch 0 taken 45 times.
✓ Branch 1 taken 18 times.
 95 for (int i = 0; i < dimWorld; ++i)
373 {
374
3/6
✓ Branch 1 taken 45 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 45 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 45 times.
✗ Branch 8 not taken.
 201 positions[i].push_back(lowerLeft[i]);
375
3/6
✓ Branch 1 taken 45 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 45 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 45 times.
✗ Branch 8 not taken.
 201 positions[i].push_back(upperRight[i]);
376
3/6
✓ Branch 1 taken 45 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 45 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 45 times.
✗ Branch 8 not taken.
 201 cells[i].push_back(numPores[i] - 1);
377
4/8
✓ Branch 1 taken 45 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 45 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 45 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 45 times.
✗ Branch 11 not taken.
 268 grading[i].resize(positions[i].size()-1, 1.0);
378
7/20
✓ Branch 1 taken 45 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 45 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 45 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 45 times.
✗ Branch 11 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 45 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 45 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 45 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
 268 grading[i] = getParamFromGroup<std::vector<CoordScalar>>(paramGroup_, "Grid.Grading" + std::to_string(i), grading[i]);
379 }
380 }
381
382 // get the lower left position
383 ✗ lowerLeft = [&]()
384 {
385 28 GlobalPosition result;
386
2/4
✓ Branch 0 taken 45 times.
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 95 for (int i = 0; i < dimWorld; ++i)
387 201 result[i] = positions[i][0];
388 28 return result;
389 28 }();
390
391 // get the upper right position
392 ✗ upperRight = [&]()
393 {
394 28 GlobalPosition result;
395
2/4
✓ Branch 0 taken 45 times.
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 95 for (int i = 0; i < dimWorld; ++i)
396 268 result[i] = positions[i].back();
397 28 return result;
398 28 }();
399
400 28 return inputData;
401 }
402
403 template<class HostGridElement, class ReferenceElement, class Geometry>
404 15773 void treatEdges_(const HostGridElement& element, const ReferenceElement& refElement, const Geometry& geometry)
405 {
406 // go over all edges and add them as elements if they passed all the tests
407
4/4
✓ Branch 0 taken 185148 times.
✓ Branch 1 taken 15517 times.
✓ Branch 2 taken 185148 times.
✓ Branch 3 taken 15517 times.
 202969 for (unsigned int edgeIdx = 0; edgeIdx < element.subEntities(dimWorld-1); ++edgeIdx)
408 {
409 187196 const auto vIdxLocal0 = refElement.subEntity(edgeIdx, dimWorld-1, 0, dimWorld);
410 187196 const auto vIdxLocal1 = refElement.subEntity(edgeIdx, dimWorld-1, 1, dimWorld);
411
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 185148 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 185148 times.
 374392 const auto edgeIdxGlobal = hostGridView_->indexSet().subIndex(element, edgeIdx, dimWorld-1);
412
413
4/4
✓ Branch 0 taken 52364 times.
✓ Branch 1 taken 132784 times.
✓ Branch 2 taken 52364 times.
✓ Branch 3 taken 132784 times.
 374392 if (edgeInserted_[edgeIdxGlobal])
414 continue;
415 else
416 53284 edgeInserted_[edgeIdxGlobal] = true;
417
418 53284 std::size_t directionNumber = 0;
419
420 if(dimWorld == 2 ) // 2D
421 {
422
2/2
✓ Branch 0 taken 166 times.
✓ Branch 1 taken 166 times.
 652 if(edgeIdx < 2) // y-direction
423 directionNumber = 0;
424 else // x-direction
425 326 directionNumber = 1;
426 }
427 else // 3D
428 {
429
2/2
✓ Branch 0 taken 34688 times.
✓ Branch 1 taken 17344 times.
 52632 if(edgeIdx < 4) // z-direction
430 directionNumber = 2;
431
4/4
✓ Branch 0 taken 32729 times.
✓ Branch 1 taken 1959 times.
✓ Branch 2 taken 17344 times.
✓ Branch 3 taken 15385 times.
 35088 else if(edgeIdx == 4 || edgeIdx == 5 || edgeIdx == 8 || edgeIdx == 9) // y-direction
432 directionNumber = 1;
433 else if(edgeIdx == 6 || edgeIdx == 7 || edgeIdx == 10 || edgeIdx == 11) // x-direction
434 directionNumber = 0;
435 }
436 106568 maybeInsertElementAndVertices_(element, geometry, vIdxLocal0, vIdxLocal1, directionProbability_[directionNumber]);
437 }
438 15773 }
439
440 template<class HostGridElement, class ReferenceElement, class Geometry>
441 ✗ void treatDiagonals_(const HostGridElement& element, const ReferenceElement& refElement, const Geometry& geometry)
442 {
443 if constexpr (dimWorld < 3)
444 132 treatDiagonalConnections2D_(element, geometry);
445 else
446 15385 treatDiagonalConnections3D_(element, refElement, geometry);
447 ✗ }
448
449 template<class HostGridElement, class Geometry>
450 void treatDiagonalConnections2D_(const HostGridElement& element,
451 const Geometry& geometry)
452 {
453 static constexpr std::array<Diagonal, 2> diagonals{ Diagonal{std::make_pair(0, 3), 2},
454 Diagonal{std::make_pair(1, 2), 3} };
455
456 132 treatIntersectingDiagonals_(element, geometry, diagonals);
457 }
458
459 template<class HostGridElement, class ReferenceElement, class Geometry>
460 15385 void treatDiagonalConnections3D_(const HostGridElement& element,
461 const ReferenceElement& refElement,
462 const Geometry& geometry)
463 {
464 // set diagonals on host grid element faces
465
4/4
✓ Branch 0 taken 92310 times.
✓ Branch 1 taken 15385 times.
✓ Branch 2 taken 92310 times.
✓ Branch 3 taken 15385 times.
 107695 for (auto faceIdx = 0; faceIdx < element.subEntities(dimWorld-2); ++faceIdx)
466 {
467
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 92310 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 92310 times.
 184620 const auto faceIdxGlobal = hostGridView_->indexSet().subIndex(element, faceIdx, dimWorld-2);
468 // face already checked?
469
4/4
✓ Branch 0 taken 43326 times.
✓ Branch 1 taken 48984 times.
✓ Branch 2 taken 43326 times.
✓ Branch 3 taken 48984 times.
 184620 if (faceInserted_[faceIdxGlobal])
470 43326 continue;
471 else
472 48984 faceInserted_[faceIdxGlobal] = true;
473
474 // get local vertex indices
475 std::array<unsigned int, 4> vIdxLocal;
476
2/2
✓ Branch 0 taken 195936 times.
✓ Branch 1 taken 48984 times.
 244920 for (int i = 0; i < 4; i++)
477 195936 vIdxLocal[i] = refElement.subEntity(faceIdx, dimWorld-2, i, dimWorld);
478
479 ✗ const auto directionNumbers = [&]()
480 {
481 48984 if (faceIdx < 2)
482 ✗ return std::make_pair<unsigned int, unsigned int>(8,7);
483
2/4
✓ Branch 0 taken 16328 times.
✓ Branch 1 taken 16328 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 32656 else if (faceIdx < 4)
484 ✗ return std::make_pair<unsigned int, unsigned int>(6,5);
485 else
486 ✗ return std::make_pair<unsigned int, unsigned int>(4,3);
487
2/4
✓ Branch 0 taken 32656 times.
✓ Branch 1 taken 16328 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 48984 }();
488
489 146952 const std::array<Diagonal, 2> diagonals{ Diagonal{std::make_pair(vIdxLocal[1], vIdxLocal[2]), directionNumbers.first},
490 146952 Diagonal{std::make_pair(vIdxLocal[0], vIdxLocal[3]), directionNumbers.second} };
491
492 48984 treatIntersectingDiagonals_(element, geometry, diagonals);
493 }
494
495 static constexpr std::array<Diagonal, 4> diagonals{ Diagonal{std::make_pair(0, 7), 9},
496 Diagonal{std::make_pair(3, 4), 10},
497 Diagonal{std::make_pair(1, 6), 11},
498 Diagonal{std::make_pair(2, 5), 12}, };
499
500 15385 treatIntersectingDiagonals_(element, geometry, diagonals);
501 15385 }
502
503 template<class HostGridElement, class Geometry, class Diagonals>
504 128998 void treatIntersectingDiagonals_(const HostGridElement& element,
505 const Geometry& geometry,
506 const Diagonals& diagonals)
507 {
508
4/6
✓ Branch 0 taken 27 times.
✓ Branch 1 taken 64474 times.
✓ Branch 3 taken 27 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 27 times.
✗ Branch 7 not taken.
 128998 static const bool allowIntersectingDiagonals = getParamFromGroup<bool>(paramGroup_, "Grid.AllowIntersectingDiagonals", true);
509
510 288770 auto treat = [&](const Diagonal& diagonal)
511 {
512 159772 return maybeInsertElementAndVertices_(element, geometry,
513 319544 diagonal.localVertexIndices.first, diagonal.localVertexIndices.second,
514 319544 directionProbability_[diagonal.directionNumber]);
515 };
516
517
1/2
✓ Branch 0 taken 64501 times.
✗ Branch 1 not taken.
 128998 if (allowIntersectingDiagonals)
518 {
519 // insert all diagonals
520
2/2
✓ Branch 0 taken 159772 times.
✓ Branch 1 taken 64501 times.
 706530 for (const auto& diagonal : diagonals)
521 319536 treat(diagonal);
522 }
523 else
524 {
525 ✗ auto order = createOrderedList_(diagonals.size());
526 ✗ std::shuffle(order.begin(), order.end(), generator_);
527 ✗ for (auto i : order)
528 {
529 ✗ const auto& diagonal = diagonals[i];
530 ✗ const bool inserted = treat(diagonal);
531 ✗ if (inserted)
532 ✗ return;
533 }
534 }
535 }
536
537 ✗ std::vector<std::size_t> createOrderedList_(const std::size_t size) const
538 {
539 ✗ std::vector<std::size_t> result(size);
540 ✗ std::iota(result.begin(), result.end(), 0);
541 ✗ return result;
542 }
543
544 auto getLambdaToRemoveThroatsOnBoundary_() const
545 {
546 2117 auto negletElementsOnBoundary = [&](const GlobalPosition& a, const GlobalPosition& b)
547 {
548 74698 const auto center = 0.5 * (a + b);
549 37653 const auto eps = (a-b).two_norm() * 1e-5;
550
551 37349 bool neglectElement = false;
552
12/12
✓ Branch 0 taken 210310 times.
✓ Branch 1 taken 31611 times.
✓ Branch 2 taken 210310 times.
✓ Branch 3 taken 31611 times.
✓ Branch 4 taken 1170 times.
✓ Branch 5 taken 138 times.
✓ Branch 6 taken 1170 times.
✓ Branch 7 taken 138 times.
✓ Branch 8 taken 638 times.
✓ Branch 9 taken 238 times.
✓ Branch 10 taken 638 times.
✓ Branch 11 taken 238 times.
 318803 for (auto i : removeThroatsOnBoundary_)
553 {
554
13/15
✓ Branch 0 taken 36811 times.
✓ Branch 1 taken 36290 times.
✓ Branch 2 taken 35770 times.
✓ Branch 3 taken 35052 times.
✓ Branch 4 taken 34091 times.
✓ Branch 5 taken 32958 times.
✓ Branch 6 taken 184 times.
✓ Branch 7 taken 168 times.
✓ Branch 8 taken 156 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 172 times.
✓ Branch 11 taken 164 times.
✓ Branch 12 taken 156 times.
✓ Branch 13 taken 146 times.
✗ Branch 14 not taken.
 212118 switch(i)
555 {
556 111651 case 0: neglectElement = center[0] < hostGridLowerLeft_[0] + eps; break;
557 110028 case 1: neglectElement = center[0] > hostGridUpperRight_[0] - eps; break;
558 36132 case 2: if constexpr (dimWorld > 1)
559 {
560 72264 neglectElement = center[1] < hostGridLowerLeft_[1] + eps;
561 36132 break;
562 }
563 35148 case 3: if constexpr (dimWorld > 1)
564 {
565 70296 neglectElement = center[1] > hostGridUpperRight_[1] - eps;
566 35148 break;
567 }
568 34687 case 4: if constexpr (dimWorld > 2)
569 {
570 68074 neglectElement = center[2] < hostGridLowerLeft_[2] + eps;
571 34037 break;
572 }
573 32908 case 5: if constexpr (dimWorld > 2)
574 {
575 65816 neglectElement = center[2] > hostGridUpperRight_[2] - eps;
576 32908 break;
577 }
578 }
579
580
6/6
✓ Branch 0 taken 205098 times.
✓ Branch 1 taken 5212 times.
✓ Branch 2 taken 1074 times.
✓ Branch 3 taken 96 times.
✓ Branch 4 taken 584 times.
✓ Branch 5 taken 54 times.
 212118 if (neglectElement)
581 return false;
582 }
583
584 return true;
585 };
586
587 return negletElementsOnBoundary;
588 }
589
590 std::string paramGroup_;
591 GlobalPosition hostGridLowerLeft_;
592 GlobalPosition hostGridUpperRight_;
593 std::vector<std::size_t> removeThroatsOnBoundary_;
594 std::unique_ptr<const HostGridView> hostGridView_;
595 std::function<bool(const GlobalPosition&, const GlobalPosition&)> considerLowDimElement_;
596
597 std::vector<GlobalPosition> vertexSet_;
598 std::vector<VertexIndexPair> elementSet_;
599
600 std::vector<bool> vertexInserted_;
601 std::vector<std::size_t> hostVertexIdxToVertexIdx_;
602 std::vector<std::size_t> edgeInserted_;
603 std::vector<std::size_t> faceInserted_;
604
605 mutable std::mt19937 generator_;
606 std::uniform_real_distribution<> uniformdist_{0, 1};
607 std::array<double, numDirections_()> directionProbability_;
608
609 std::shared_ptr<Grid> gridPtr_;
610 };
611
612 } // end namespace Dumux::PoreNetwork
613
614 #endif // HAVE_DUNE_FOAMGRID
615
616 #endif
617