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 Helper class to assign parameters to a generated grid | ||
11 | */ | ||
12 | #ifndef DUMUX_IO_PARAMETERS_FOR_GENERATED_GRID | ||
13 | #define DUMUX_IO_PARAMETERS_FOR_GENERATED_GRID | ||
14 | |||
15 | #include <algorithm> | ||
16 | #include <array> | ||
17 | #include <numeric> | ||
18 | #include <random> | ||
19 | #include <string> | ||
20 | #include <vector> | ||
21 | |||
22 | #include <dune/common/exceptions.hh> | ||
23 | #include <dune/grid/common/exceptions.hh> | ||
24 | #include <dune/geometry/axisalignedcubegeometry.hh> | ||
25 | |||
26 | #include <dumux/common/parameters.hh> | ||
27 | #include <dumux/common/random.hh> | ||
28 | #include <dumux/common/stringutilities.hh> | ||
29 | #include <dumux/geometry/intersectspointgeometry.hh> | ||
30 | #include <dumux/porenetwork/common/throatproperties.hh> | ||
31 | #include <dumux/porenetwork/common/poreproperties.hh> | ||
32 | |||
33 | namespace Dumux::PoreNetwork { | ||
34 | |||
35 | /*! | ||
36 | * \ingroup PoreNetworkModels | ||
37 | * \brief Helper class to assign parameters to a generated grid | ||
38 | */ | ||
39 | template <class Grid, class Scalar> | ||
40 | ✗ | class ParametersForGeneratedGrid | |
41 | { | ||
42 | using GridView = typename Grid::LeafGridView; | ||
43 | using Element = typename Grid::template Codim<0>::Entity; | ||
44 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
45 | using Vertex = typename Grid::template Codim<Grid::dimension>::Entity; | ||
46 | |||
47 | static constexpr auto dim = Grid::dimension; | ||
48 | static constexpr auto dimWorld = Grid::dimensionworld; | ||
49 | using BoundaryList = std::array<int, 2*dimWorld>; | ||
50 | |||
51 | public: | ||
52 | |||
53 | 33 | ParametersForGeneratedGrid(const GridView& gridView, const std::string& paramGroup) | |
54 | : gridView_(gridView) | ||
55 | , paramGroup_(paramGroup) | ||
56 |
1/5✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
99 | , priorityList_(getPriorityList_()) |
57 | { | ||
58 |
1/2✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
|
33 | computeBoundingBox_(); |
59 |
1/2✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
|
33 | boundaryFaceIndex_ = getBoundaryFacemarkerInput_(); |
60 | 33 | } | |
61 | |||
62 | /*! | ||
63 | * \brief Returns the boundary face marker index at given position | ||
64 | * | ||
65 | * \param pos The current position | ||
66 | */ | ||
67 | int boundaryFaceMarkerAtPos(const GlobalPosition& pos) const | ||
68 | { | ||
69 | // set the priority which decides the order the vertices on the boundary are indexed | ||
70 | // by default, vertices on min/max faces in x direction have the highest priority, followed by y and z | ||
71 |
6/12✓ Branch 0 taken 78 times.
✓ Branch 1 taken 36 times.
✓ Branch 2 taken 88616 times.
✓ Branch 3 taken 12659 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 568 times.
✓ Branch 9 taken 72 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
118980 | for (auto boundaryIdx : priorityList_) |
72 | { | ||
73 |
6/12✓ Branch 0 taken 4 times.
✓ Branch 1 taken 74 times.
✓ Branch 2 taken 4058 times.
✓ Branch 3 taken 84558 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 122 times.
✓ Branch 9 taken 446 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
89262 | if (onBoundary_(pos, boundaryIdx)) |
74 | 8368 | return boundaryFaceIndex_[boundaryIdx]; | |
75 | } | ||
76 | return -1; | ||
77 | } | ||
78 | |||
79 | /*! | ||
80 | * \brief Computes and returns the label of a given throat | ||
81 | * | ||
82 | * \param poreLabels The labels of the pores adjacent to the throat. | ||
83 | */ | ||
84 | 44289 | int throatLabel(const std::array<int, 2>& poreLabels) const | |
85 | { | ||
86 |
6/6✓ Branch 0 taken 9795 times.
✓ Branch 1 taken 34040 times.
✓ Branch 2 taken 9795 times.
✓ Branch 3 taken 34040 times.
✓ Branch 4 taken 9795 times.
✓ Branch 5 taken 34040 times.
|
132867 | if (poreLabels[0] == poreLabels[1]) // both vertices are inside the domain or on the same boundary face |
87 | return poreLabels[0]; | ||
88 |
4/4✓ Branch 0 taken 5385 times.
✓ Branch 1 taken 4410 times.
✓ Branch 2 taken 5385 times.
✓ Branch 3 taken 4410 times.
|
20042 | if (poreLabels[0] == -1) // vertex1 is inside the domain, vertex2 is on a boundary face |
89 | return poreLabels[1]; | ||
90 |
4/4✓ Branch 0 taken 898 times.
✓ Branch 1 taken 4487 times.
✓ Branch 2 taken 898 times.
✓ Branch 3 taken 4487 times.
|
11054 | if (poreLabels[1] == -1) // vertex2 is inside the domain, vertex1 is on a boundary face |
91 | return poreLabels[0]; | ||
92 | |||
93 | // use the priority list to find out which pore label is favored | ||
94 |
1/2✓ Branch 0 taken 1182 times.
✗ Branch 1 not taken.
|
3104 | for(const auto i : priorityList_) |
95 | { | ||
96 |
6/6✓ Branch 0 taken 700 times.
✓ Branch 1 taken 482 times.
✓ Branch 2 taken 700 times.
✓ Branch 3 taken 482 times.
✓ Branch 4 taken 700 times.
✓ Branch 5 taken 482 times.
|
3708 | if (poreLabels[0] == boundaryFaceIndex_[i]) |
97 | return poreLabels[0]; | ||
98 |
6/6✓ Branch 0 taken 284 times.
✓ Branch 1 taken 416 times.
✓ Branch 2 taken 284 times.
✓ Branch 3 taken 416 times.
✓ Branch 4 taken 284 times.
✓ Branch 5 taken 416 times.
|
2196 | if (poreLabels[1] == boundaryFaceIndex_[i]) |
99 | return poreLabels[1]; | ||
100 | } | ||
101 | |||
102 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Something went wrong with the throat labels"); | |
103 | } | ||
104 | |||
105 | /*! | ||
106 | * \brief Assign parameters for generically created grids | ||
107 | */ | ||
108 | template <class SetParameter, class GetParameter> | ||
109 | 33 | void assignParameters(const SetParameter& setParameter, | |
110 | const GetParameter& getParameter, | ||
111 | const std::size_t numSubregions) | ||
112 | { | ||
113 | using cytpe = typename GlobalPosition::value_type; | ||
114 | using InternalBoundingBox = Dune::AxisAlignedCubeGeometry<cytpe, dimWorld, dimWorld>; | ||
115 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 20 times.
|
33 | std::vector<InternalBoundingBox> internalBoundingBoxes; |
116 | |||
117 | // divide the network into subregions, if specified | ||
118 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 20 times.
|
33 | if (numSubregions > 0) |
119 | { | ||
120 | // get bounding boxes of subregions | ||
121 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
|
2 | for (int i = 0; i < numSubregions; ++i) |
122 | { | ||
123 |
7/20✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
|
3 | auto lowerLeft = getParamFromGroup<GlobalPosition>(paramGroup_, "Grid.Subregion" + std::to_string(i) + ".LowerLeft"); |
124 |
7/22✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✗ 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.
✗ Branch 25 not taken.
|
3 | auto upperRight = getParamFromGroup<GlobalPosition>(paramGroup_, "Grid.Subregion" + std::to_string(i) + ".UpperRight"); |
125 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | internalBoundingBoxes.emplace_back(std::move(lowerLeft), std::move(upperRight)); |
126 | } | ||
127 | } | ||
128 | |||
129 | // get the maximum possible pore body radii such that pore bodies do not intersect | ||
130 | // (requires ThroatRegionId, if specified) | ||
131 |
3/4✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 20 times.
|
66 | const std::vector<Scalar> maxPoreRadius = getMaxPoreRadii_(numSubregions, getParameter); |
132 |
3/8✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 21 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 21 times.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
99 | std::vector<bool> poreRadiusLimited(gridView_.size(dim), false); |
133 | |||
134 | // get a helper function for getting the pore radius of a pore body not belonging to a subregion | ||
135 |
1/4✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
66 | auto defaultPoreRadius = poreRadiusGenerator_(-1); |
136 | |||
137 | // get helper functions for pore body radii on subregions | ||
138 |
1/2✓ Branch 0 taken 21 times.
✗ Branch 1 not taken.
|
66 | std::vector<decltype(defaultPoreRadius)> subregionPoreRadius; |
139 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 21 times.
|
34 | for (int i = 0; i < numSubregions; ++i) |
140 |
3/8✓ 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 not taken.
✗ Branch 9 not taken.
|
1 | subregionPoreRadius.emplace_back(poreRadiusGenerator_(i)); |
141 | |||
142 | // get helper function for pore volume | ||
143 |
1/2✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
|
33 | const auto poreVolume = poreVolumeGenerator_(getParameter); |
144 | |||
145 | // treat the pore body parameters (label, radius and maybe regionId) | ||
146 |
5/6✓ Branch 1 taken 16951 times.
✓ Branch 2 taken 21 times.
✓ Branch 3 taken 16951 times.
✓ Branch 4 taken 21 times.
✓ Branch 6 taken 16951 times.
✗ Branch 7 not taken.
|
34679 | for (const auto& vertex : vertices(gridView_)) |
147 | { | ||
148 |
2/4✓ Branch 1 taken 16951 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 16951 times.
|
17323 | const auto& pos = vertex.geometry().center(); |
149 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 16951 times.
|
17323 | const auto vIdxGlobal = gridView_.indexSet().index(vertex); |
150 | 17323 | const auto poreLabel = boundaryFaceMarkerAtPos(pos); | |
151 |
6/12✓ Branch 1 taken 16951 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16951 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16951 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 16951 times.
✓ Branch 11 taken 15951 times.
✓ Branch 12 taken 1000 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
|
34646 | setParameter(vertex, "PoreLabel", poreLabel); |
152 | |||
153 | // sets the minimum of the given value and the maximum possible pore body radius | ||
154 | // and keeps track of capped radii | ||
155 | 17820 | auto setRadiusAndLogIfCapped = [&](const Scalar value) | |
156 | { | ||
157 |
8/12✓ Branch 0 taken 5188 times.
✓ Branch 1 taken 11556 times.
✓ Branch 2 taken 5188 times.
✓ Branch 3 taken 11556 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 138 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 138 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 194 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 194 times.
|
34152 | if (value > maxPoreRadius[vIdxGlobal]) |
158 | { | ||
159 |
2/12✓ Branch 1 taken 5188 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5188 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
|
10376 | poreRadiusLimited[vIdxGlobal] = true; |
160 |
4/30✓ Branch 1 taken 5188 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5188 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5188 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 5188 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 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.
|
20752 | setParameter(vertex, "PoreInscribedRadius", maxPoreRadius[vIdxGlobal]); |
161 | } | ||
162 | else | ||
163 |
12/30✓ Branch 1 taken 11556 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11556 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 11556 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 11556 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 14 taken 138 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 138 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 138 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 138 times.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✓ Branch 27 taken 194 times.
✗ Branch 28 not taken.
✓ Branch 30 taken 194 times.
✗ Branch 31 not taken.
✓ Branch 33 taken 194 times.
✗ Branch 34 not taken.
✓ Branch 35 taken 194 times.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
|
47552 | setParameter(vertex, "PoreInscribedRadius", value); |
164 | }; | ||
165 | |||
166 |
2/2✓ Branch 0 taken 15951 times.
✓ Branch 1 taken 1000 times.
|
17323 | if (numSubregions == 0) // assign radius if no subregions are specified |
167 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 15951 times.
✓ Branch 3 taken 15951 times.
✗ Branch 4 not taken.
|
32646 | setRadiusAndLogIfCapped(defaultPoreRadius(vertex, poreLabel)); |
168 | else // assign region ids and radii to vertices if they are within a subregion | ||
169 | { | ||
170 | // default value for vertices not belonging to a subregion | ||
171 |
5/12✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1000 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1000 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 1000 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1000 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
|
2000 | setParameter(vertex, "PoreRegionId", -1); |
172 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 1000 times.
✓ Branch 3 taken 1000 times.
✗ Branch 4 not taken.
|
2000 | setRadiusAndLogIfCapped(defaultPoreRadius(vertex, poreLabel)); |
173 | |||
174 |
2/2✓ Branch 0 taken 1000 times.
✓ Branch 1 taken 1000 times.
|
2000 | for (int id = 0; id < numSubregions; ++id) |
175 | { | ||
176 |
1/2✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
|
1000 | const auto& subregion = internalBoundingBoxes[id]; |
177 |
3/8✓ Branch 1 taken 1000 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 125 times.
✓ Branch 6 taken 875 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
1000 | if (intersectsPointGeometry(vertex.geometry().center(), subregion)) |
178 | { | ||
179 |
5/12✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 125 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 125 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 125 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 125 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
|
250 | setParameter(vertex, "PoreRegionId", id); |
180 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 125 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 125 times.
✓ Branch 5 taken 125 times.
✗ Branch 6 not taken.
|
375 | setRadiusAndLogIfCapped(subregionPoreRadius[id](vertex, poreLabel)); |
181 | } | ||
182 | } | ||
183 | } | ||
184 | |||
185 |
5/12✓ Branch 1 taken 16951 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16951 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16951 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 16951 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 16951 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
|
17323 | setParameter(vertex, "PoreVolume", poreVolume(vertex, poreLabel)); |
186 | } | ||
187 | |||
188 | // treat throat parameters | ||
189 |
1/2✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
|
33 | auto defaultThroatInscribedRadius = throatInscribedRadiusGenerator_(-1, getParameter); |
190 |
1/2✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
|
33 | auto defaultThroatLength = throatLengthGenerator_(-1, getParameter); |
191 | |||
192 | // get helper functions for throat radii and lengths on subregions | ||
193 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
33 | std::vector<decltype(defaultThroatInscribedRadius)> subregionThroatInscribedRadius; |
194 |
4/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 20 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 20 times.
|
67 | std::vector<decltype(defaultThroatLength)> subregionThroatLength; |
195 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 21 times.
|
34 | for (int i = 0; i < numSubregions; ++i) |
196 | { | ||
197 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
1 | subregionThroatInscribedRadius.emplace_back(throatInscribedRadiusGenerator_(i, getParameter)); |
198 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1 | subregionThroatLength.emplace_back(throatLengthGenerator_(i, getParameter)); |
199 | } | ||
200 | |||
201 | // set throat parameters | ||
202 |
6/6✓ Branch 1 taken 43835 times.
✓ Branch 2 taken 21 times.
✓ Branch 3 taken 43835 times.
✓ Branch 4 taken 21 times.
✓ Branch 5 taken 33359 times.
✓ Branch 6 taken 10476 times.
|
88611 | for (const auto& element : elements(gridView_)) |
203 | { | ||
204 |
2/2✓ Branch 0 taken 33359 times.
✓ Branch 1 taken 10476 times.
|
44289 | if (numSubregions == 0) // assign values if no subregions are specified |
205 | { | ||
206 |
6/14✓ Branch 1 taken 33359 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 33359 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 33359 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 33359 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 33359 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 33359 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
|
67626 | setParameter(element, "ThroatInscribedRadius", defaultThroatInscribedRadius(element)); |
207 |
5/12✓ Branch 1 taken 33359 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 33359 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 33359 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 33359 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 33359 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
33813 | setParameter(element, "ThroatLength", defaultThroatLength(element)); |
208 | } | ||
209 | else // assign values to throats if they are within a subregion | ||
210 | { | ||
211 | // default value for elements not belonging to a subregion | ||
212 |
5/12✓ Branch 1 taken 10476 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10476 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 10476 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 10476 times.
✓ Branch 12 taken 10476 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
20952 | setParameter(element, "ThroatRegionId", -1); |
213 |
6/14✓ Branch 1 taken 10476 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10476 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 10476 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 10476 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 10476 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 10476 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
|
20952 | setParameter(element, "ThroatInscribedRadius", defaultThroatInscribedRadius(element)); |
214 |
5/12✓ Branch 1 taken 10476 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10476 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 10476 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 10476 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 10476 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
10476 | setParameter(element, "ThroatLength", defaultThroatLength(element)); |
215 | |||
216 |
2/2✓ Branch 0 taken 10476 times.
✓ Branch 1 taken 10476 times.
|
20952 | for (int id = 0; id < numSubregions; ++id) |
217 | { | ||
218 |
1/2✓ Branch 1 taken 10476 times.
✗ Branch 2 not taken.
|
10476 | const auto& subregion = internalBoundingBoxes[id]; |
219 |
3/8✓ Branch 1 taken 10476 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1625 times.
✓ Branch 5 taken 8851 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
20952 | if (intersectsPointGeometry(element.geometry().center(), subregion)) |
220 | { | ||
221 |
5/12✓ Branch 1 taken 1625 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1625 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1625 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 1625 times.
✓ Branch 12 taken 1625 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
3250 | setParameter(element, "ThroatRegionId", id); |
222 |
7/16✓ Branch 1 taken 1625 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1625 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1625 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1625 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1625 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 1625 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 1625 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
|
4875 | setParameter(element, "ThroatInscribedRadius", subregionThroatInscribedRadius[id](element)); |
223 |
6/14✓ Branch 1 taken 1625 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1625 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1625 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1625 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1625 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 1625 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
|
3250 | setParameter(element, "ThroatLength", subregionThroatLength[id](element)); |
224 | } | ||
225 | } | ||
226 | } | ||
227 | |||
228 | // set the throat label | ||
229 |
1/2✓ Branch 1 taken 43835 times.
✗ Branch 2 not taken.
|
44289 | const auto vertex0 = element.template subEntity<dim>(0); |
230 |
1/2✓ Branch 1 taken 43835 times.
✗ Branch 2 not taken.
|
44289 | const auto vertex1 = element.template subEntity<dim>(1); |
231 |
5/12✓ Branch 1 taken 43835 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43835 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 43835 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 43835 times.
✓ Branch 12 taken 43835 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
88578 | const std::array poreLabels{static_cast<int>(getParameter(vertex0, "PoreLabel")), |
232 |
5/12✓ Branch 1 taken 43835 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43835 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 43835 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 43835 times.
✓ Branch 12 taken 43835 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
88578 | static_cast<int>(getParameter(vertex1, "PoreLabel"))}; |
233 |
5/12✓ Branch 1 taken 43835 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43835 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 43835 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 43835 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 43835 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
|
44289 | setParameter(element, "ThroatLabel", throatLabel(poreLabels)); |
234 | } | ||
235 | |||
236 | 99 | const auto numPoreRadiusLimited = std::count(poreRadiusLimited.begin(), poreRadiusLimited.end(), true); | |
237 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 19 times.
|
33 | if (numPoreRadiusLimited > 0) |
238 |
3/6✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
|
8 | std::cout << "*******\nWarning! " << numPoreRadiusLimited << " out of " << poreRadiusLimited.size() |
239 |
1/2✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
|
4 | << " pore body radii have been capped automatically in order to prevent intersecting pores\n*******" << std::endl; |
240 | 33 | } | |
241 | |||
242 | private: | ||
243 | |||
244 | /*! | ||
245 | * \brief Returns a list of boundary face priorities from user specified input or default values if no input is given | ||
246 | * | ||
247 | * This essentially determines the index of a node on an edge or corner. For instance in a 2D case, a list of {0,1,2,3} | ||
248 | * will give highest priority to the "x-direction" and lowest to the "diagonal-direction". | ||
249 | */ | ||
250 | BoundaryList getPriorityList_() const | ||
251 | { | ||
252 | 33 | const auto list = [&]() | |
253 | { | ||
254 | BoundaryList priorityList; | ||
255 | 99 | std::iota(priorityList.begin(), priorityList.end(), 0); | |
256 | |||
257 |
7/12✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 21 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 21 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 3 times.
✓ Branch 9 taken 18 times.
✓ Branch 10 taken 3 times.
✓ Branch 11 taken 18 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
165 | if (hasParamInGroup(paramGroup_, "Grid.PriorityList")) |
258 | { | ||
259 | try { | ||
260 | // priorities can also be set in the input file | ||
261 | 3 | priorityList = getParamFromGroup<BoundaryList>(paramGroup_, "Grid.PriorityList"); | |
262 | } | ||
263 | // make sure that a priority for each direction is set | ||
264 | ✗ | catch(Dune::RangeError& e) { | |
265 | ✗ | DUNE_THROW(Dumux::ParameterException, "You must specify priorities for all directions (" << dimWorld << ") \n" << e.what()); | |
266 | } | ||
267 | // make sure each direction is only set once | ||
268 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
6 | if (!isUnique_(priorityList)) |
269 | ✗ | DUNE_THROW(Dumux::ParameterException, "You must specify priorities for all directions (duplicate directions)"); | |
270 | |||
271 | //make sure that the directions are correct (ranging from 0 to dimWorld-1) | ||
272 |
10/42✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 3 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 3 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 3 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 3 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 3 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 3 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.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
|
21 | if (std::any_of(priorityList.begin(), priorityList.end(), []( const int i ){ return (i < 0 || i >= 2*dimWorld); })) |
273 | ✗ | DUNE_THROW(Dumux::ParameterException, "You must specify priorities for correct directions (0-" << 2*dimWorld-1 << ")"); | |
274 | } | ||
275 | 33 | return priorityList; | |
276 |
4/7✓ Branch 1 taken 23 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 15 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
|
57 | }(); |
277 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
12 | return list; |
278 | } | ||
279 | |||
280 | 33 | void computeBoundingBox_() | |
281 | { | ||
282 | // calculate the bounding box of the local partition of the grid view | ||
283 |
4/4✓ Branch 1 taken 17105 times.
✓ Branch 2 taken 21 times.
✓ Branch 3 taken 17105 times.
✓ Branch 4 taken 21 times.
|
17510 | for (const auto& vertex : vertices(gridView_)) |
284 | { | ||
285 |
2/2✓ Branch 0 taken 51030 times.
✓ Branch 1 taken 17105 times.
|
69349 | for (int i = 0; i < dimWorld; i++) |
286 | { | ||
287 | using std::min; | ||
288 | using std::max; | ||
289 |
8/8✓ Branch 1 taken 63 times.
✓ Branch 2 taken 50967 times.
✓ Branch 3 taken 63 times.
✓ Branch 4 taken 50967 times.
✓ Branch 5 taken 63 times.
✓ Branch 6 taken 50967 times.
✓ Branch 7 taken 61 times.
✓ Branch 8 taken 50929 times.
|
51963 | bBoxMin_[i] = min(bBoxMin_[i], vertex.geometry().corner(0)[i]); |
290 |
8/8✓ Branch 1 taken 326 times.
✓ Branch 2 taken 50704 times.
✓ Branch 3 taken 326 times.
✓ Branch 4 taken 50704 times.
✓ Branch 5 taken 326 times.
✓ Branch 6 taken 50704 times.
✓ Branch 7 taken 288 times.
✓ Branch 8 taken 50702 times.
|
52700 | bBoxMax_[i] = max(bBoxMax_[i], vertex.geometry().corner(0)[i]); |
291 | } | ||
292 | } | ||
293 | 33 | } | |
294 | |||
295 | /*! | ||
296 | * \brief Returns a list of boundary face indices from user specified input or default values if no input is given | ||
297 | */ | ||
298 | 33 | BoundaryList getBoundaryFacemarkerInput_() const | |
299 | { | ||
300 | BoundaryList boundaryFaceMarker; | ||
301 | 99 | std::fill(boundaryFaceMarker.begin(), boundaryFaceMarker.end(), 0); | |
302 |
1/2✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
|
33 | boundaryFaceMarker[0] = 1; |
303 |
1/2✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
|
33 | boundaryFaceMarker[1] = 1; |
304 | |||
305 |
8/14✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 21 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 21 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 21 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 19 times.
✓ Branch 12 taken 2 times.
✓ Branch 13 taken 19 times.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
|
165 | if (hasParamInGroup(paramGroup_, "Grid.BoundaryPoreLabels")) |
306 | { | ||
307 | 62 | const auto input = getParamFromGroup<std::vector<std::string>>(paramGroup_, "Grid.BoundaryPoreLabels"); | |
308 |
4/4✓ Branch 0 taken 88 times.
✓ Branch 1 taken 19 times.
✓ Branch 2 taken 88 times.
✓ Branch 3 taken 19 times.
|
365 | for (const auto& entry : input) |
309 | { | ||
310 |
2/4✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 88 times.
✗ Branch 5 not taken.
|
408 | const std::string errorMessage = "You must specify BoundaryPoreLabels in the format pos:num, where pos can be xMin, xMax, yMin, yMax, zMin, zMax and num is the corresponding label.\n" |
311 | "Example (2D, defaults are used for the remaining boundaries): xMin:2 yMax:3\n"; | ||
312 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 88 times.
|
136 | if (entry.find(':') == std::string::npos) |
313 | ✗ | DUNE_THROW(Dumux::ParameterException, errorMessage); | |
314 | |||
315 |
14/28✓ Branch 0 taken 19 times.
✓ Branch 1 taken 69 times.
✓ Branch 3 taken 19 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 19 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 19 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 19 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 19 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 19 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 19 times.
✗ Branch 22 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 19 times.
✗ Branch 27 not taken.
✓ Branch 28 taken 19 times.
✓ Branch 31 taken 114 times.
✓ Branch 32 taken 19 times.
✗ Branch 33 not taken.
✓ Branch 34 taken 114 times.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
|
322 | static const std::map<std::string, int> labels = {{"xMin", 0}, {"xMax", 1}, {"yMin", 2}, {"yMax", 3}, {"zMin", 4}, {"zMax", 5}}; |
316 |
3/10✓ Branch 2 taken 88 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 88 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 88 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
272 | const auto splitEntry = split(entry, ":"); |
317 |
5/12✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 88 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 88 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 88 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 88 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
544 | const std::string location = std::string(splitEntry[0].begin(), splitEntry[0].end()); |
318 | 136 | int value = 0; | |
319 | try { | ||
320 |
7/16✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 88 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 88 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 88 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 88 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 88 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 88 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
|
544 | value = std::stoi(std::string(splitEntry[1].begin(), splitEntry[1].end())); |
321 | } | ||
322 | ✗ | catch(...) { | |
323 | ✗ | DUNE_THROW(Dumux::ParameterException, errorMessage); | |
324 | } | ||
325 | |||
326 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 88 times.
|
272 | if (splitEntry.size() != 2) |
327 | ✗ | DUNE_THROW(Dumux::ParameterException, errorMessage); | |
328 | 136 | if (!labels.count(location)) | |
329 | ✗ | DUNE_THROW(Dumux::ParameterException, errorMessage); | |
330 | else | ||
331 |
3/6✓ Branch 1 taken 88 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 88 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 88 times.
|
136 | boundaryFaceMarker[labels.at(location)] = value; |
332 | } | ||
333 | } | ||
334 | 33 | return boundaryFaceMarker; | |
335 | } | ||
336 | |||
337 | // returns the maximum possible pore body radii such that pore bodies do not intersect | ||
338 | template <class GetParameter> | ||
339 | 33 | std::vector<Scalar> getMaxPoreRadii_(std::size_t numSubregions, const GetParameter& getParameter) const | |
340 | { | ||
341 | 33 | const auto numVertices = gridView_.size(dim); | |
342 |
1/2✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
|
66 | std::vector<Scalar> maxPoreRadius(numVertices, std::numeric_limits<Scalar>::max()); |
343 | |||
344 |
2/4✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
|
33 | if (!getParamFromGroup<bool>(paramGroup_, "Grid.CapPoreRadii", true)) |
345 | return maxPoreRadius; | ||
346 | |||
347 | // check for a user-specified fixed throat length | ||
348 |
1/4✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
33 | const Scalar inputThroatLength = getParamFromGroup<Scalar>(paramGroup_, "Grid.ThroatLength", -1.0); |
349 |
4/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 20 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 20 times.
|
66 | std::vector<Scalar> subregionInputThroatLengths; |
350 | |||
351 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 20 times.
|
33 | if (numSubregions > 0) |
352 | { | ||
353 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
|
2 | for (int i = 0; i < numSubregions; ++i) |
354 | { | ||
355 | // adapt the parameter group if there are subregions | ||
356 |
5/14✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
|
3 | const std::string paramGroup = paramGroup_ + ".SubRegion" + std::to_string(i); |
357 |
1/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
1 | const Scalar subregionInputThroatLength = getParamFromGroup<Scalar>(paramGroup, "Grid.ThroatLength", -1.0); |
358 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | subregionInputThroatLengths.push_back(subregionInputThroatLength); |
359 | } | ||
360 | } | ||
361 | |||
362 |
6/6✓ Branch 1 taken 43835 times.
✓ Branch 2 taken 21 times.
✓ Branch 3 taken 43835 times.
✓ Branch 4 taken 21 times.
✓ Branch 5 taken 10476 times.
✓ Branch 6 taken 33359 times.
|
44322 | for (const auto& element : elements(gridView_)) |
363 | { | ||
364 | // do not cap the pore radius if a user-specified throat length if given | ||
365 |
2/2✓ Branch 0 taken 10476 times.
✓ Branch 1 taken 33359 times.
|
44289 | if (numSubregions > 0) |
366 | { | ||
367 | // check subregions | ||
368 |
5/12✓ Branch 1 taken 10476 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10476 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 10476 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 10476 times.
✓ Branch 11 taken 10476 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
|
20952 | const auto subregionId = getParameter(element, "ThroatRegionId"); |
369 |
1/2✓ Branch 0 taken 10476 times.
✗ Branch 1 not taken.
|
10476 | if (subregionId >= 0) |
370 | { | ||
371 | // throat lies within a subregion | ||
372 |
2/4✓ Branch 0 taken 10476 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10476 times.
✗ Branch 3 not taken.
|
20952 | if (subregionInputThroatLengths[subregionId] > 0.0) |
373 | 246 | continue; | |
374 | } | ||
375 | ✗ | else if (inputThroatLength > 0.0) // throat does not lie within a subregion | |
376 | continue; | ||
377 | } | ||
378 |
2/2✓ Branch 0 taken 33236 times.
✓ Branch 1 taken 123 times.
|
33813 | else if (inputThroatLength > 0.0) // no subregions |
379 | continue; | ||
380 | |||
381 | // No fixed throat lengths given, check for max. pore radius. | ||
382 | // We define this as 1/2 of the length (minus a user specified value) of the shortest pore throat attached to the pore body. | ||
383 |
3/4✓ Branch 1 taken 43712 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 18 times.
✓ Branch 4 taken 43694 times.
|
44043 | const Scalar delta = element.geometry().volume(); |
384 |
4/6✓ Branch 0 taken 18 times.
✓ Branch 1 taken 43694 times.
✓ Branch 3 taken 18 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 18 times.
✗ Branch 7 not taken.
|
44043 | static const Scalar minThroatLength = getParamFromGroup<Scalar>(paramGroup_, "Grid.MinThroatLength", 1e-6); |
385 | 44043 | const Scalar maxRadius = (delta - minThroatLength)/2.0; | |
386 |
2/2✓ Branch 0 taken 87424 times.
✓ Branch 1 taken 43712 times.
|
132129 | for (int vIdxLocal = 0; vIdxLocal < 2 ; ++vIdxLocal) |
387 | { | ||
388 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 87424 times.
✓ Branch 3 taken 87424 times.
✗ Branch 4 not taken.
|
88086 | const int vIdxGlobal = gridView_.indexSet().subIndex(element, vIdxLocal, dim); |
389 |
4/4✓ Branch 0 taken 24192 times.
✓ Branch 1 taken 63232 times.
✓ Branch 2 taken 24192 times.
✓ Branch 3 taken 63232 times.
|
200633 | maxPoreRadius[vIdxGlobal] = std::min(maxPoreRadius[vIdxGlobal], maxRadius); |
390 | } | ||
391 | } | ||
392 | |||
393 | return maxPoreRadius; | ||
394 | } | ||
395 | |||
396 | // returns a function taking a vertex and a pore label and returning a radius | ||
397 | 34 | std::function<Scalar(const Vertex&, const int)> poreRadiusGenerator_(const int subregionId) const | |
398 | { | ||
399 | // adapt the parameter name if there are subregions | ||
400 |
11/20✓ Branch 0 taken 21 times.
✓ Branch 1 taken 1 times.
✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 21 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✓ Branch 18 taken 21 times.
✓ Branch 19 taken 1 times.
✓ Branch 20 taken 21 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
|
69 | const std::string prefix = subregionId < 0 ? "Grid." : "Grid.Subregion" + std::to_string(subregionId) + "."; |
401 | |||
402 | // prepare random number generation for lognormal parameter distribution | ||
403 | 34 | std::mt19937 generator; | |
404 | |||
405 | // check if pores for certain labels should be treated in a special way | ||
406 |
7/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 9 taken 22 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 22 times.
✓ Branch 13 taken 1 times.
✓ Branch 14 taken 21 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
|
170 | const auto poreLabelsToSetFixedRadius = getParamFromGroup<std::vector<int>>(paramGroup_, prefix + "PoreLabelsToSetFixedRadius", std::vector<int>{}); |
407 |
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 9 taken 22 times.
✗ Branch 10 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 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
|
170 | const auto poreLabelsToApplyFactorForRadius = getParamFromGroup<std::vector<int>>(paramGroup_, prefix + "PoreLabelsToApplyFactorForRadius", std::vector<int>{}); |
408 |
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 9 taken 22 times.
✗ Branch 10 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 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
|
170 | const auto poreRadiusForLabel = getParamFromGroup<std::vector<Scalar>>(paramGroup_, prefix + "FixedPoreRadiusForLabel", std::vector<Scalar>{}); |
409 |
7/18✓ 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 9 taken 22 times.
✗ Branch 10 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.
|
170 | const auto poreRadiusFactorForLabel = getParamFromGroup<std::vector<Scalar>>(paramGroup_, prefix + "PoreRadiusFactorForLabel", std::vector<Scalar>{}); |
410 | |||
411 | 130 | const auto generateFunction = [&](auto& poreRadiusDist) | |
412 | { | ||
413 | 33780 | return [=](const auto& vertex, const int poreLabel) mutable | |
414 | { | ||
415 | 33780 | const auto radius = poreRadiusDist(generator); | |
416 | |||
417 | // check if pores for certain labels should be treated in a special way | ||
418 |
4/8✓ Branch 0 taken 16704 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 16704 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 16704 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 16704 times.
|
67560 | if (poreLabelsToSetFixedRadius.empty() && poreLabelsToApplyFactorForRadius.empty()) |
419 | return radius; // nothing special to be done | ||
420 | |||
421 | // set a fixed radius for a given label | ||
422 |
0/8✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
22 | else if (!poreLabelsToSetFixedRadius.empty() || !poreRadiusForLabel.empty()) |
423 | { | ||
424 | ✗ | if (poreLabelsToSetFixedRadius.size() != poreRadiusForLabel.size()) | |
425 | ✗ | DUNE_THROW(Dumux::ParameterException, "PoreLabelsToSetFixedRadius must be of same size as FixedPoreRadiusForLabel"); | |
426 | |||
427 | ✗ | if (const auto it = std::find(poreLabelsToSetFixedRadius.begin(), poreLabelsToSetFixedRadius.end(), poreLabel); it != poreLabelsToSetFixedRadius.end()) | |
428 | ✗ | return poreRadiusForLabel[std::distance(poreLabelsToSetFixedRadius.begin(), it)]; | |
429 | } | ||
430 | |||
431 | // multiply the pore radius by a given value for a given label | ||
432 |
0/8✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
22 | else if (!poreLabelsToApplyFactorForRadius.empty() || !poreRadiusFactorForLabel.empty()) |
433 | { | ||
434 | ✗ | if (poreLabelsToApplyFactorForRadius.size() != poreRadiusFactorForLabel.size()) | |
435 | ✗ | DUNE_THROW(Dumux::ParameterException, "PoreLabelsToApplyFactorForRadius must be of same size as PoreRadiusFactorForLabel"); | |
436 | |||
437 | ✗ | if (const auto it = std::find(poreLabelsToApplyFactorForRadius.begin(), poreLabelsToApplyFactorForRadius.end(), poreLabel); it != poreLabelsToApplyFactorForRadius.end()) | |
438 | ✗ | return poreRadiusFactorForLabel[std::distance(poreLabelsToApplyFactorForRadius.begin(), it)] * radius; | |
439 | } | ||
440 | |||
441 | // default | ||
442 | return radius; | ||
443 | 22 | }; | |
444 | }; | ||
445 | |||
446 |
3/8✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 22 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
68 | const Scalar fixedPoreRadius = getParamFromGroup<Scalar>(paramGroup_, prefix + "PoreInscribedRadius", -1.0); |
447 | // return random radius according to a user-specified distribution | ||
448 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 16 times.
|
34 | if (fixedPoreRadius <= 0.0) |
449 | { | ||
450 | // allow to specify a seed to get reproducible results | ||
451 |
8/20✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 6 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 6 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 6 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 6 times.
✗ Branch 19 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 6 times.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
|
36 | const auto seed = getParamFromGroup<unsigned int>(paramGroup_, prefix + "ParameterRandomNumberSeed", std::random_device{}()); |
452 | 9 | generator.seed(seed); | |
453 | |||
454 |
3/6✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 6 times.
|
18 | const auto type = getParamFromGroup<std::string>(paramGroup_, prefix + "ParameterType", "lognormal"); |
455 |
1/2✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
|
9 | if (type == "lognormal") |
456 | { | ||
457 | // if we use a lognormal distribution, get the mean and standard deviation from input file | ||
458 |
15/36✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 6 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 6 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 6 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 6 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 6 times.
✓ Branch 23 taken 6 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 6 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 6 times.
✗ Branch 28 not taken.
✓ Branch 29 taken 6 times.
✗ Branch 30 not taken.
✓ Branch 32 taken 6 times.
✗ Branch 33 not taken.
✓ Branch 35 taken 6 times.
✗ Branch 36 not taken.
✓ Branch 38 taken 6 times.
✗ 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 not taken.
|
63 | const auto [meanPoreRadius, stddevPoreRadius] = getDistributionInputParams_("Lognormal", prefix, |
459 | "MeanPoreInscribedRadius", | ||
460 | "StandardDeviationPoreInscribedRadius"); | ||
461 | 9 | const Scalar variance = stddevPoreRadius*stddevPoreRadius; | |
462 | |||
463 | using std::log; | ||
464 | using std::sqrt; | ||
465 | 9 | const Scalar mu = log(meanPoreRadius/sqrt(1.0 + variance/(meanPoreRadius*meanPoreRadius))); | |
466 | 9 | const Scalar sigma = sqrt(log(1.0 + variance/(meanPoreRadius*meanPoreRadius))); | |
467 | |||
468 |
1/2✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
|
9 | Dumux::SimpleLogNormalDistribution<> poreRadiusDist(mu, sigma); |
469 |
2/4✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
|
9 | return generateFunction(poreRadiusDist); |
470 | } | ||
471 | ✗ | else if (type == "uniform") | |
472 | { | ||
473 | // if we use a uniform distribution, get the min and max from input file | ||
474 | ✗ | const auto [minPoreRadius, maxPoreRadius] = getDistributionInputParams_("Uniform", prefix, | |
475 | "MinPoreInscribedRadius", | ||
476 | "MaxPoreInscribedRadius"); | ||
477 | ✗ | Dumux::SimpleUniformDistribution<> poreRadiusDist(minPoreRadius, maxPoreRadius); | |
478 | ✗ | return generateFunction(poreRadiusDist); | |
479 | } | ||
480 | else | ||
481 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Unknown parameter type " << type); | |
482 | } | ||
483 | |||
484 | // always return a fixed constant radius | ||
485 | else | ||
486 | { | ||
487 | 25 | auto poreRadiusDist = [fixedPoreRadius](auto& gen){ return fixedPoreRadius; }; | |
488 |
2/6✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
25 | return generateFunction(poreRadiusDist); |
489 | } | ||
490 | } | ||
491 | |||
492 | // print helpful error message if params are not properly provided | ||
493 | 9 | std::array<Scalar, 2> getDistributionInputParams_(const std::string& distributionName, | |
494 | const std::string& prefix, | ||
495 | const std::string& paramName0, | ||
496 | const std::string& paramName1) const | ||
497 | { | ||
498 | try | ||
499 | { | ||
500 |
2/6✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
9 | return std::array{getParamFromGroup<Scalar>(paramGroup_, prefix + paramName0), |
501 |
4/12✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 6 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
36 | getParamFromGroup<Scalar>(paramGroup_, prefix + paramName1)}; |
502 | } | ||
503 | ✗ | catch (const Dumux::ParameterException& e) | |
504 | { | ||
505 | std::cout << "\n" << distributionName << " pore-size distribution needs input parameters " | ||
506 | << prefix + paramName0 << " and " << prefix + paramName1 << ".\n" | ||
507 | ✗ | << "Alternatively, use " << prefix << "PoreInscribedRadius to set a fixed inscribed pore radius." << std::endl; | |
508 | ✗ | DUNE_THROW(Dumux::ParameterException, e.what()); | |
509 | } | ||
510 | } | ||
511 | |||
512 | template <class GetParameter> | ||
513 | 33 | auto poreVolumeGenerator_(const GetParameter& getParameter) const | |
514 | { | ||
515 |
2/6✓ Branch 2 taken 21 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 21 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
33 | const auto geometry = Pore::shapeFromString(getParamFromGroup<std::string>(paramGroup_, "Grid.PoreGeometry")); |
516 |
2/6✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 21 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
66 | const auto capPoresOnBoundaries = getParamFromGroup<std::vector<int>>(paramGroup_, "Grid.CapPoresOnBoundaries", std::vector<int>{}); |
517 |
4/6✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 18 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 21 times.
|
37 | if (!isUnique_(capPoresOnBoundaries)) |
518 | ✗ | DUNE_THROW(Dune::InvalidStateException, "CapPoresOnBoundaries must not contain duplicates"); | |
519 | |||
520 | // automatically determine the pore volume if not provided by the grid file | ||
521 | 16951 | return [=] (const auto& vertex, const auto vIdx) | |
522 | { | ||
523 |
4/10✓ Branch 1 taken 16579 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16579 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16579 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 16579 times.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
50853 | const Scalar r = getParameter(vertex, "PoreInscribedRadius"); |
524 | 16951 | const Scalar volume = [&] | |
525 | { | ||
526 | 16951 | if (geometry == Pore::Shape::cylinder) | |
527 | { | ||
528 |
0/4✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
42 | static const Scalar fixedHeight = getParamFromGroup<Scalar>(paramGroup_, "Grid.PoreHeight", -1.0); |
529 |
0/10✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
40 | const Scalar h = fixedHeight > 0.0 ? fixedHeight : getParameter(vertex, "PoreHeight"); |
530 | 80 | return Pore::volume(geometry, r, h); | |
531 | } | ||
532 | else | ||
533 | 16911 | return Pore::volume(geometry, r); | |
534 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 16579 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
18109 | }(); |
535 | |||
536 |
4/4✓ Branch 0 taken 14778 times.
✓ Branch 1 taken 1801 times.
✓ Branch 2 taken 14778 times.
✓ Branch 3 taken 1801 times.
|
33902 | if (capPoresOnBoundaries.empty()) |
537 | return volume; | ||
538 | else | ||
539 | { | ||
540 | 14847 | std::size_t numCaps = 0; | |
541 | 14847 | Scalar factor = 1.0; | |
542 | 14847 | const auto& pos = vertex.geometry().center(); | |
543 |
4/4✓ Branch 0 taken 88668 times.
✓ Branch 1 taken 14778 times.
✓ Branch 2 taken 88668 times.
✓ Branch 3 taken 14778 times.
|
133554 | for (auto boundaryIdx : capPoresOnBoundaries) |
544 | { | ||
545 |
2/2✓ Branch 0 taken 3218 times.
✓ Branch 1 taken 85450 times.
|
89013 | if (onBoundary_(pos, boundaryIdx)) |
546 | { | ||
547 | 3266 | factor *= 0.5; | |
548 | 3266 | ++numCaps; | |
549 | } | ||
550 | } | ||
551 | |||
552 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 14778 times.
|
14847 | if (numCaps > dimWorld) |
553 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Pore " << vIdx << " at " << pos << " capped " << numCaps << " times. Capping should not happen more than " << dimWorld << " times"); | |
554 | |||
555 | 14847 | return volume * factor; | |
556 | } | ||
557 |
10/14✓ Branch 0 taken 2 times.
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 11 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✓ Branch 10 taken 7 times.
✓ Branch 11 taken 8 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
|
17012 | }; |
558 | } | ||
559 | |||
560 | // returns a lambda taking a element and returning a radius | ||
561 | template <class GetParameter> | ||
562 | 34 | auto throatInscribedRadiusGenerator_(const int subregionId, const GetParameter& getParameter) const | |
563 | { | ||
564 | // adapt the parameter name if there are subregions | ||
565 |
11/20✓ Branch 0 taken 21 times.
✓ Branch 1 taken 1 times.
✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 21 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✓ Branch 18 taken 21 times.
✓ Branch 19 taken 1 times.
✓ Branch 20 taken 21 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
|
35 | const std::string prefix = subregionId < 0 ? "Grid." : "Grid.Subregion" + std::to_string(subregionId) + "."; |
566 | |||
567 | // check for a user-specified fixed throat radius | ||
568 |
3/8✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 22 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
68 | const Scalar inputThroatInscribedRadius = getParamFromGroup<Scalar>(paramGroup_, prefix + "ThroatInscribedRadius", -1.0); |
569 | |||
570 | // shape parameter for calculation of throat radius | ||
571 |
3/10✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 22 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
68 | const Scalar throatN = getParamFromGroup<Scalar>(paramGroup_, prefix + "ThroatInscribedRadiusN", 0.1); |
572 | |||
573 | 45460 | return [=](const Element& element) | |
574 | { | ||
575 |
2/2✓ Branch 1 taken 43592 times.
✓ Branch 2 taken 1414 times.
|
45460 | const Scalar delta = element.geometry().volume(); |
576 |
4/4✓ Branch 0 taken 43592 times.
✓ Branch 1 taken 1414 times.
✓ Branch 2 taken 43592 times.
✓ Branch 3 taken 1414 times.
|
90920 | const std::array<Vertex, 2> vertices = {element.template subEntity<dim>(0), element.template subEntity<dim>(1)}; |
577 | |||
578 | // the element parameters (throat radius and length) | ||
579 |
2/2✓ Branch 0 taken 43592 times.
✓ Branch 1 taken 1414 times.
|
45460 | if (inputThroatInscribedRadius > 0.0) |
580 | return inputThroatInscribedRadius; | ||
581 | else | ||
582 | { | ||
583 |
6/14✓ Branch 1 taken 43592 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43592 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 43592 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 43592 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 43592 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 43592 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
|
131769 | const Scalar poreRadius0 = getParameter(vertices[0], "PoreInscribedRadius"); |
584 |
5/12✓ Branch 1 taken 43592 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43592 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 43592 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 43592 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 43592 times.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
|
131769 | const Scalar poreRadius1 = getParameter(vertices[1], "PoreInscribedRadius"); |
585 | 43923 | return Throat::averagedRadius(poreRadius0, poreRadius1, delta, throatN); | |
586 | } | ||
587 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 21 times.
|
35 | }; |
588 | } | ||
589 | |||
590 | // returns a lambda taking a element and returning a length | ||
591 | template <class GetParameter> | ||
592 | 34 | auto throatLengthGenerator_(int subregionId, const GetParameter& getParameter) const | |
593 | { | ||
594 | // adapt the parameter name if there are subregions | ||
595 |
11/20✓ Branch 0 taken 21 times.
✓ Branch 1 taken 1 times.
✓ Branch 3 taken 21 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 21 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✓ Branch 18 taken 21 times.
✓ Branch 19 taken 1 times.
✓ Branch 20 taken 21 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
|
35 | const std::string prefix = subregionId < 0 ? "Grid." : "Grid.Subregion" + std::to_string(subregionId) + "."; |
596 |
3/8✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 22 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
68 | const Scalar inputThroatLength = getParamFromGroup<Scalar>(paramGroup_, prefix + "ThroatLength", -1.0); |
597 | // decide whether to subtract the pore radii from the throat length or not | ||
598 |
3/10✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 22 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
68 | const bool subtractRadiiFromThroatLength = getParamFromGroup<bool>(paramGroup_, prefix + "SubtractPoreInscribedRadiiFromThroatLength", true); |
599 | |||
600 | 45460 | return [=](const Element& element) | |
601 | { | ||
602 |
1/2✓ Branch 0 taken 45006 times.
✗ Branch 1 not taken.
|
45460 | if (inputThroatLength > 0.0) |
603 | return inputThroatLength; | ||
604 | |||
605 |
1/2✓ Branch 1 taken 45006 times.
✗ Branch 2 not taken.
|
45337 | const Scalar delta = element.geometry().volume(); |
606 | 45337 | if (subtractRadiiFromThroatLength) | |
607 | { | ||
608 |
2/4✓ Branch 1 taken 45006 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 45006 times.
✗ Branch 5 not taken.
|
90674 | const std::array<Vertex, 2> vertices = {element.template subEntity<dim>(0), element.template subEntity<dim>(1)}; |
609 |
12/28✓ Branch 1 taken 45006 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 45006 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 45006 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 45006 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 45006 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 45006 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 45006 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 45006 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 45006 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 45006 times.
✗ Branch 27 not taken.
✓ Branch 28 taken 45006 times.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✓ Branch 31 taken 45006 times.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
|
226685 | const Scalar result = delta - getParameter(vertices[0], "PoreInscribedRadius") - getParameter(vertices[1], "PoreInscribedRadius"); |
610 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 45006 times.
|
45337 | if (result <= 0.0) |
611 | ✗ | DUNE_THROW(Dune::GridError, "Pore radii are so large they intersect! Something went wrong at element " << gridView_.indexSet().index(element)); | |
612 | else | ||
613 | return result; | ||
614 | } | ||
615 | else | ||
616 | return delta; | ||
617 |
3/6✓ Branch 0 taken 1 times.
✓ Branch 1 taken 21 times.
✓ Branch 2 taken 45006 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
45372 | }; |
618 | } | ||
619 | |||
620 | bool onBoundary_(const GlobalPosition& pos, std::size_t boundaryIdx) const | ||
621 | { | ||
622 | 178275 | constexpr auto eps = 1e-8; //TODO | |
623 | |||
624 | static constexpr std::array boundaryIdxToCoordinateIdx{0, 0, 1, 1, 2, 2}; | ||
625 |
10/18✓ Branch 0 taken 38 times.
✓ Branch 1 taken 40 times.
✓ Branch 2 taken 43034 times.
✓ Branch 3 taken 45022 times.
✓ Branch 4 taken 44592 times.
✓ Branch 5 taken 44636 times.
✓ Branch 6 taken 138 times.
✓ Branch 7 taken 207 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 254 times.
✓ Branch 13 taken 314 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
|
178275 | const auto coordinateIdx = boundaryIdxToCoordinateIdx[boundaryIdx]; |
626 | 178275 | auto isMaxBoundary = [] (int n) { return (n % 2 != 0); }; | |
627 | |||
628 |
10/18✓ Branch 0 taken 38 times.
✓ Branch 1 taken 40 times.
✓ Branch 2 taken 43034 times.
✓ Branch 3 taken 45022 times.
✓ Branch 4 taken 44592 times.
✓ Branch 5 taken 44636 times.
✓ Branch 6 taken 138 times.
✓ Branch 7 taken 207 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 254 times.
✓ Branch 13 taken 314 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
|
178275 | if (isMaxBoundary(boundaryIdx)) |
629 | 264092 | return pos[coordinateIdx] > bBoxMax_[coordinateIdx] - eps; | |
630 | else | ||
631 | 270577 | return pos[coordinateIdx] < bBoxMin_[coordinateIdx] + eps; | |
632 | } | ||
633 | |||
634 | // check if a container is unique | ||
635 | // copy the container such that is does not get altered by std::sort | ||
636 | template <class T> | ||
637 | ✗ | bool isUnique_(T t) const | |
638 | { | ||
639 | 9 | std::sort(t.begin(), t.end()); | |
640 |
2/4✗ Branch 3 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 3 times.
|
9 | return (std::unique(t.begin(), t.end()) == t.end()); |
641 | } | ||
642 | |||
643 | const GridView gridView_; | ||
644 | |||
645 | const std::string paramGroup_; | ||
646 | BoundaryList priorityList_; | ||
647 | BoundaryList boundaryFaceIndex_; | ||
648 | |||
649 | //! the bounding box of the whole domain | ||
650 | GlobalPosition bBoxMin_ = GlobalPosition(std::numeric_limits<double>::max()); | ||
651 | GlobalPosition bBoxMax_ = GlobalPosition(std::numeric_limits<double>::min()); | ||
652 | }; | ||
653 | |||
654 | } // namespace Dumux::PoreNetwork | ||
655 | |||
656 | #endif | ||
657 |