| 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-FileCopyrightText: 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 PoreNetworkDiscretization | ||
| 10 | * \brief Base class for the finite volume geometry for porenetwork models | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_DISCRETIZATION_PNM_GRID_GEOMETRY_HH | ||
| 13 | #define DUMUX_DISCRETIZATION_PNM_GRID_GEOMETRY_HH | ||
| 14 | |||
| 15 | #include <string> | ||
| 16 | #include <utility> | ||
| 17 | #include <unordered_map> | ||
| 18 | #include <functional> | ||
| 19 | |||
| 20 | #include <dune/common/exceptions.hh> | ||
| 21 | #include <dune/localfunctions/lagrange/lagrangelfecache.hh> | ||
| 22 | |||
| 23 | #include <dumux/discretization/method.hh> | ||
| 24 | #include <dumux/common/indextraits.hh> | ||
| 25 | #include <dumux/common/defaultmappertraits.hh> | ||
| 26 | #include <dumux/common/parameters.hh> | ||
| 27 | |||
| 28 | #include <dumux/discretization/basegridgeometry.hh> | ||
| 29 | #include <dumux/discretization/porenetwork/fvelementgeometry.hh> | ||
| 30 | #include <dumux/discretization/porenetwork/subcontrolvolume.hh> | ||
| 31 | #include <dumux/discretization/porenetwork/subcontrolvolumeface.hh> | ||
| 32 | #include <dumux/porenetwork/common/throatproperties.hh> | ||
| 33 | #include <dumux/porenetwork/common/poreproperties.hh> | ||
| 34 | #include <dumux/discretization/extrusion.hh> | ||
| 35 | |||
| 36 | namespace Dumux::PoreNetwork { | ||
| 37 | |||
| 38 | /*! | ||
| 39 | * \ingroup PoreNetworkDiscretization | ||
| 40 | * \brief Base class for geometry data extraction from the grid data format | ||
| 41 | */ | ||
| 42 | template<class Scalar, class GridView> | ||
| 43 | class DefaultPNMData | ||
| 44 | { | ||
| 45 | using GridIndex = typename IndexTraits<GridView>::GridIndex; | ||
| 46 | using SmallLocalIndex = typename IndexTraits<GridView>::SmallLocalIndex; | ||
| 47 | using Label = std::int_least8_t; | ||
| 48 | using Vertex = typename GridView::template Codim<GridView::dimension>::Entity; | ||
| 49 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 50 | |||
| 51 | static const int dim = GridView::dimension; | ||
| 52 | |||
| 53 | public: | ||
| 54 | |||
| 55 | template<class GridData> | ||
| 56 | 28 | void update(const GridView& gridView, const GridData& gridData) | |
| 57 | { | ||
| 58 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 28 times.
|
28 | coordinationNumber_ = gridData.getCoordinationNumbers(); |
| 59 | |||
| 60 | 28 | const auto numThroats = gridView.size(0); | |
| 61 | 28 | throatInscribedRadius_.resize(numThroats); | |
| 62 | 28 | throatLength_.resize(numThroats); | |
| 63 | 28 | throatLabel_.resize(numThroats); | |
| 64 | 28 | throatCrossSectionalArea_.resize(numThroats); | |
| 65 | 28 | throatShapeFactor_.resize(numThroats); | |
| 66 | |||
| 67 | 28 | useSameGeometryForAllPores_ = true; | |
| 68 | 28 | useSameShapeForAllThroats_ = true; | |
| 69 | 28 | overwriteGridDataWithShapeSpecificValues_ = false; | |
| 70 | |||
| 71 | // first check if the same geometry shall be used for all entities ... | ||
| 72 |
3/4✓ Branch 2 taken 28 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 24 times.
✓ Branch 6 taken 4 times.
|
56 | if (hasParamInGroup(gridData.paramGroup(), "Grid.ThroatCrossSectionShape")) |
| 73 | { | ||
| 74 | 24 | const auto throatGeometryInput = getParamFromGroup<std::string>(gridData.paramGroup(), "Grid.ThroatCrossSectionShape"); | |
| 75 |
1/2✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
|
24 | const auto throatGeometry = Throat::shapeFromString(throatGeometryInput); |
| 76 |
1/2✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
|
24 | throatGeometry_.resize(1); |
| 77 |
1/2✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
|
24 | throatGeometry_[0] = throatGeometry; |
| 78 |
2/4✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
|
24 | overwriteGridDataWithShapeSpecificValues_ = getParamFromGroup<bool>(gridData.paramGroup(), "Grid.OverwriteGridDataWithShapeSpecificValues", true); |
| 79 | |||
| 80 |
3/6✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 24 times.
✗ Branch 8 not taken.
|
24 | std::cout << "Using '" << throatGeometryInput << "' as cross-sectional shape for all throats." << std::endl; |
| 81 | 24 | } | |
| 82 | else // .. otherwise, get the corresponding parameter index from the grid data and resize the respective vector | ||
| 83 | { | ||
| 84 | 4 | std::cout << "Reading shape factors for throats from grid data." << std::endl; | |
| 85 | 4 | useSameShapeForAllThroats_ = false; | |
| 86 | 4 | throatGeometry_.resize(numThroats); | |
| 87 | } | ||
| 88 | |||
| 89 | // get the vertex parameters | ||
| 90 | 28 | const auto numPores = gridView.size(dim); | |
| 91 | 28 | poreInscribedRadius_.resize(numPores); | |
| 92 | 28 | poreLabel_.resize(numPores); | |
| 93 | 28 | poreVolume_.resize(numPores); | |
| 94 | |||
| 95 | // first check if the same geometry shall be used for all entities ... | ||
| 96 |
2/4✓ Branch 2 taken 28 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 28 times.
✗ Branch 6 not taken.
|
56 | if (hasParamInGroup(gridData.paramGroup(), "Grid.PoreGeometry")) |
| 97 | { | ||
| 98 | 28 | const auto poreGeometryInput = getParamFromGroup<std::string>(gridData.paramGroup(), "Grid.PoreGeometry"); | |
| 99 |
1/2✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
|
28 | poreGeometry_.resize(1); |
| 100 |
2/4✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
|
28 | poreGeometry_[0] = Pore::shapeFromString(poreGeometryInput);; |
| 101 | |||
| 102 |
3/6✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 28 times.
✗ Branch 8 not taken.
|
28 | std::cout << "Using '" << poreGeometryInput << "' as geometry for all pores." << std::endl; |
| 103 | 28 | } | |
| 104 | else // .. otherwise, get the corresponding parameter index from the grid data and resize the respective vector | ||
| 105 | { | ||
| 106 | ✗ | std::cout << "Reading pore shapes from grid data." << std::endl; | |
| 107 | ✗ | useSameGeometryForAllPores_ = false; | |
| 108 | ✗ | poreGeometry_.resize(numPores); | |
| 109 | } | ||
| 110 | |||
| 111 | |||
| 112 |
2/2✓ Branch 1 taken 18855 times.
✓ Branch 2 taken 28 times.
|
37738 | for (const auto& vertex : vertices(gridView)) |
| 113 | { | ||
| 114 |
5/8✓ Branch 0 taken 25 times.
✓ Branch 1 taken 18830 times.
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 25 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 23 times.
✗ Branch 10 not taken.
|
18857 | static const auto poreInscribedRadiusIdx = gridData.parameterIndex("PoreInscribedRadius"); |
| 115 |
5/8✓ Branch 0 taken 25 times.
✓ Branch 1 taken 18830 times.
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 25 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 23 times.
✗ Branch 10 not taken.
|
18857 | static const auto poreLabelIdx = gridData.parameterIndex("PoreLabel"); |
| 116 | 18855 | const auto vIdx = gridView.indexSet().index(vertex); | |
| 117 | 18855 | const auto& params = gridData.parameters(vertex); | |
| 118 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 18855 times.
|
18855 | poreInscribedRadius_[vIdx] = params[poreInscribedRadiusIdx]; |
| 119 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 18855 times.
|
18855 | assert(poreInscribedRadius_[vIdx] > 0.0); |
| 120 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 18855 times.
|
18855 | poreLabel_[vIdx] = params[poreLabelIdx]; |
| 121 | |||
| 122 | 18855 | if (!useSameGeometryForAllPores()) | |
| 123 | ✗ | poreGeometry_[vIdx] = getPoreGeometry_(gridData, vertex); | |
| 124 | |||
| 125 | 18855 | poreVolume_[vIdx] = getPoreVolume_(gridData, vertex, vIdx); | |
| 126 | } | ||
| 127 | |||
| 128 |
3/4✗ Branch 0 not taken.
✓ Branch 1 taken 42157 times.
✓ Branch 3 taken 42157 times.
✓ Branch 4 taken 28 times.
|
84342 | for (const auto& element : elements(gridView)) |
| 129 | { | ||
| 130 | 42157 | const int eIdx = gridView.indexSet().index(element); | |
| 131 | 42157 | const auto& params = gridData.parameters(element); | |
| 132 |
5/8✓ Branch 0 taken 25 times.
✓ Branch 1 taken 42132 times.
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 25 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 23 times.
✗ Branch 10 not taken.
|
42159 | static const auto throatInscribedRadiusIdx = gridData.parameterIndex("ThroatInscribedRadius"); |
| 133 |
5/8✓ Branch 0 taken 25 times.
✓ Branch 1 taken 42132 times.
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 25 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 23 times.
✗ Branch 10 not taken.
|
42159 | static const auto throatLengthIdx = gridData.parameterIndex("ThroatLength"); |
| 134 |
2/2✓ Branch 0 taken 25 times.
✓ Branch 1 taken 42132 times.
|
42157 | throatInscribedRadius_[eIdx] = params[throatInscribedRadiusIdx]; |
| 135 | 42157 | throatLength_[eIdx] = params[throatLengthIdx]; | |
| 136 | |||
| 137 | // use a default value if no throat label is given by the grid | ||
| 138 |
4/6✓ Branch 0 taken 25 times.
✓ Branch 1 taken 42132 times.
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 25 times.
✗ Branch 7 not taken.
|
42157 | static const bool gridHasThroatLabel = gridData.gridHasElementParameter("ThroatLabel"); |
| 139 |
1/2✓ Branch 0 taken 42157 times.
✗ Branch 1 not taken.
|
42157 | if (gridHasThroatLabel) |
| 140 | { | ||
| 141 |
5/8✓ Branch 0 taken 25 times.
✓ Branch 1 taken 42132 times.
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 25 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 23 times.
✗ Branch 10 not taken.
|
42159 | static const auto throatLabelIdx = gridData.parameterIndex("ThroatLabel"); |
| 142 | 42157 | throatLabel_[eIdx] = params[throatLabelIdx]; | |
| 143 | } | ||
| 144 | else | ||
| 145 | { | ||
| 146 | ✗ | const auto vIdx0 = gridView.indexSet().subIndex(element, 0, dim); | |
| 147 | ✗ | const auto vIdx1 = gridView.indexSet().subIndex(element, 1, dim); | |
| 148 | |||
| 149 | ✗ | const auto poreLabel0 = poreLabel(vIdx0); | |
| 150 | ✗ | const auto poreLabel1 = poreLabel(vIdx1); | |
| 151 | |||
| 152 | ✗ | if (poreLabel0 >= 0 && poreLabel1 >= 0) | |
| 153 | { | ||
| 154 | ✗ | std::cout << "\n Warning: Throat " | |
| 155 | << eIdx << " connects two boundary pores with different pore labels. Using the greater pore label as throat label.\n" | ||
| 156 | ✗ | << "Set the throat labels explicitly in your grid file, if needed." << std::endl; | |
| 157 | } | ||
| 158 | |||
| 159 | using std::max; | ||
| 160 | ✗ | throatLabel_[eIdx] = max(poreLabel0, poreLabel1); | |
| 161 | } | ||
| 162 | |||
| 163 | 42157 | if (!useSameShapeForAllThroats()) | |
| 164 | { | ||
| 165 |
4/8✓ Branch 0 taken 1 times.
✓ Branch 1 taken 6381 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
6383 | static const auto throatShapeFactorIdx = gridData.parameterIndex("ThroatShapeFactor"); |
| 166 |
4/8✓ Branch 0 taken 1 times.
✓ Branch 1 taken 6381 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
6383 | static const auto throatAreaIdx = gridData.parameterIndex("ThroatCrossSectionalArea"); |
| 167 |
2/2✓ Branch 0 taken 4122 times.
✓ Branch 1 taken 2260 times.
|
6382 | throatShapeFactor_[eIdx] = params[throatShapeFactorIdx]; |
| 168 |
2/2✓ Branch 0 taken 4122 times.
✓ Branch 1 taken 2260 times.
|
10504 | throatGeometry_[eIdx] = Throat::shape(throatShapeFactor_[eIdx]); |
| 169 | 6382 | throatCrossSectionalArea_[eIdx] = params[throatAreaIdx]; | |
| 170 | } | ||
| 171 | else | ||
| 172 | { | ||
| 173 | 35775 | throatCrossSectionalArea_[eIdx] = getThroatCrossSectionalArea_(gridData, element, eIdx); | |
| 174 | 35775 | throatShapeFactor_[eIdx] = getThroatShapeFactor_(gridData, element, eIdx); | |
| 175 | } | ||
| 176 | |||
| 177 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 42157 times.
|
42157 | assert(throatInscribedRadius_[eIdx] > 0.0); |
| 178 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 42157 times.
|
42157 | assert(throatLength_[eIdx] > 0.0); |
| 179 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 42157 times.
|
42157 | assert(throatCrossSectionalArea_[eIdx] > 0.0); |
| 180 | |||
| 181 |
4/6✓ Branch 0 taken 25 times.
✓ Branch 1 taken 42132 times.
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 25 times.
✗ Branch 7 not taken.
|
42157 | static const bool addThroatVolumeToPoreVolume = getParamFromGroup<bool>(gridData.paramGroup(), "Grid.AddThroatVolumeToPoreVolume", false); |
| 182 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 42157 times.
|
42157 | if (addThroatVolumeToPoreVolume) |
| 183 | { | ||
| 184 | ✗ | for (int vIdxLocal = 0; vIdxLocal < 2; ++vIdxLocal) | |
| 185 | { | ||
| 186 | ✗ | const auto vIdx = gridView.indexSet().subIndex(element, vIdxLocal, dim); | |
| 187 | ✗ | poreVolume_[vIdx] += 0.5 * throatCrossSectionalArea_[eIdx] * throatLength_[eIdx]; | |
| 188 | } | ||
| 189 | } | ||
| 190 | } | ||
| 191 | |||
| 192 | 28 | maybeResizeContainers_(); | |
| 193 | 28 | } | |
| 194 | |||
| 195 | //! Returns the pore label (e.g. used for setting BCs) | ||
| 196 | 12411430 | Label poreLabel(const GridIndex dofIdxGlobal) const | |
| 197 |
33/36✓ Branch 0 taken 50536 times.
✓ Branch 1 taken 1461655 times.
✓ Branch 2 taken 279250 times.
✓ Branch 3 taken 2818774 times.
✓ Branch 4 taken 6054 times.
✓ Branch 5 taken 2731 times.
✓ Branch 6 taken 25729 times.
✓ Branch 7 taken 20531 times.
✓ Branch 8 taken 22905 times.
✓ Branch 9 taken 6949 times.
✓ Branch 10 taken 79295 times.
✓ Branch 11 taken 448251 times.
✓ Branch 12 taken 61067 times.
✓ Branch 13 taken 735053 times.
✓ Branch 14 taken 11353 times.
✓ Branch 15 taken 103368 times.
✓ Branch 16 taken 2422 times.
✓ Branch 17 taken 13242 times.
✓ Branch 18 taken 46412 times.
✓ Branch 19 taken 389227 times.
✓ Branch 20 taken 1408 times.
✓ Branch 21 taken 5387 times.
✓ Branch 22 taken 1365 times.
✓ Branch 23 taken 4018 times.
✓ Branch 24 taken 805 times.
✓ Branch 25 taken 3213 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 6790 times.
✓ Branch 28 taken 1732252 times.
✓ Branch 29 taken 2101022 times.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✓ Branch 32 taken 82760 times.
✓ Branch 33 taken 715720 times.
✓ Branch 34 taken 195 times.
✓ Branch 35 taken 1819 times.
|
11697899 | { return poreLabel_[dofIdxGlobal]; } |
| 198 | |||
| 199 | //! Returns the vector of pore labels | ||
| 200 | 24 | const std::vector<Label>& poreLabel() const | |
| 201 |
2/4✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
24 | { return poreLabel_; } |
| 202 | |||
| 203 | //! Returns the inscribed radius of the pore | ||
| 204 | 27527647 | Scalar poreInscribedRadius(const GridIndex dofIdxGlobal) const | |
| 205 |
3/4✓ Branch 0 taken 2522976 times.
✓ Branch 1 taken 41563 times.
✓ Branch 2 taken 1092264 times.
✗ Branch 3 not taken.
|
14347230 | { return poreInscribedRadius_[dofIdxGlobal]; } |
| 206 | |||
| 207 | //! Returns the vector of inscribed pore radii | ||
| 208 | const std::vector<Scalar>& poreInscribedRadius() const | ||
| 209 | { return poreInscribedRadius_; } | ||
| 210 | |||
| 211 | //! Returns the volume of the pore | ||
| 212 | 23505988 | Scalar poreVolume(const GridIndex dofIdxGlobal) const | |
| 213 | 7923358 | { return poreVolume_[dofIdxGlobal]; } | |
| 214 | |||
| 215 | //! Returns the vector of pore volumes | ||
| 216 | 3 | const std::vector<Scalar>& poreVolume() const | |
| 217 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
3 | { return poreVolume_; } |
| 218 | |||
| 219 | //! Returns the inscribed radius of the throat | ||
| 220 | 42408966 | Scalar throatInscribedRadius(const GridIndex eIdx) const | |
| 221 |
1/2✓ Branch 0 taken 722 times.
✗ Branch 1 not taken.
|
42389310 | { return throatInscribedRadius_[eIdx]; } |
| 222 | |||
| 223 | //! Returns the vector of inscribed throat radii | ||
| 224 | const std::vector<Scalar>& throatInscribedRadius() const | ||
| 225 | { return throatInscribedRadius_; } | ||
| 226 | |||
| 227 | //! Returns the length of the throat | ||
| 228 | 40488187 | Scalar throatLength(const GridIndex eIdx) const | |
| 229 | 40488187 | { return throatLength_[eIdx]; } | |
| 230 | |||
| 231 | //! Returns the vector of throat lengths | ||
| 232 | const std::vector<Scalar>& throatLength() const | ||
| 233 | { return throatLength_; } | ||
| 234 | |||
| 235 | //! Returns an index indicating if a throat is touching the domain boundary | ||
| 236 | 722 | Label throatLabel(const GridIndex eIdx) const | |
| 237 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
722 | { return throatLabel_[eIdx]; } |
| 238 | |||
| 239 | //! Returns the vector of throat labels | ||
| 240 | 24 | const std::vector<Label>& throatLabel() const | |
| 241 |
2/4✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
24 | { return throatLabel_; } |
| 242 | |||
| 243 | //! Returns the number of throats connected to a pore (coordination number) | ||
| 244 | 7910694 | SmallLocalIndex coordinationNumber(const GridIndex dofIdxGlobal) const | |
| 245 | 7910694 | { return coordinationNumber_[dofIdxGlobal]; } | |
| 246 | |||
| 247 | //! Returns the vector of coordination numbers | ||
| 248 | 24 | const std::vector<SmallLocalIndex>& coordinationNumber() const | |
| 249 |
2/4✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
1368803 | { return coordinationNumber_; } |
| 250 | |||
| 251 | //! the geometry of the pore | ||
| 252 | 3063616 | Pore::Shape poreGeometry(const GridIndex vIdx) const | |
| 253 | 3064436 | { return useSameGeometryForAllPores() ? poreGeometry_[0] : poreGeometry_[vIdx]; } | |
| 254 | |||
| 255 | //! Returns the vector of pore geometries | ||
| 256 | 3062796 | const std::vector<Pore::Shape>& poreGeometry() const | |
| 257 | { | ||
| 258 | 3062796 | if (useSameGeometryForAllPores()) | |
| 259 | { | ||
| 260 | // if a vector of pore geometries is requested (e.g., for vtk output), | ||
| 261 | // resize the container and fill it with the same value everywhere | ||
| 262 | 3062796 | const auto poreGeo = poreGeometry_[0]; | |
| 263 | 3062796 | poreGeometry_.resize(poreInscribedRadius_.size(), poreGeo); | |
| 264 | } | ||
| 265 | |||
| 266 | 3062796 | return poreGeometry_; | |
| 267 | } | ||
| 268 | |||
| 269 | //! Returns the throat's cross-sectional shape | ||
| 270 | 41538200 | Throat::Shape throatCrossSectionShape(const GridIndex eIdx) const | |
| 271 | 41538200 | { return useSameShapeForAllThroats() ? throatGeometry_[0] : throatGeometry_[eIdx]; } | |
| 272 | |||
| 273 | //! Returns the vector of cross-sectional shapes | ||
| 274 | const std::vector<Throat::Shape>& throatCrossSectionShape() const | ||
| 275 | { | ||
| 276 | if (useSameShapeForAllThroats()) | ||
| 277 | { | ||
| 278 | // if a vector of throat cross section shapes is requested (e.g., for vtk output), | ||
| 279 | // resize the container and fill it with the same value everywhere | ||
| 280 | const auto throatShape = throatGeometry_[0]; | ||
| 281 | throatGeometry_.resize(throatInscribedRadius_.size(), throatShape); | ||
| 282 | } | ||
| 283 | |||
| 284 | return throatGeometry_; | ||
| 285 | } | ||
| 286 | |||
| 287 | //! Returns the throat's cross-sectional area | ||
| 288 | 44566606 | Scalar throatCrossSectionalArea(const GridIndex eIdx) const | |
| 289 |
2/2✓ Branch 0 taken 153993 times.
✓ Branch 1 taken 867363 times.
|
44564326 | { return throatCrossSectionalArea_[eIdx]; } |
| 290 | |||
| 291 | //! Returns the vector of throat cross-sectional areas | ||
| 292 | 3 | const std::vector<Scalar>& throatCrossSectionalArea() const | |
| 293 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
3 | { return throatCrossSectionalArea_; } |
| 294 | |||
| 295 | //! Returns the throat's shape factor | ||
| 296 | 41436612 | Scalar throatShapeFactor(const GridIndex eIdx) const | |
| 297 |
3/4✓ Branch 0 taken 3894043 times.
✓ Branch 1 taken 35427936 times.
✓ Branch 2 taken 1092264 times.
✗ Branch 3 not taken.
|
41436612 | { return useSameShapeForAllThroats() ? throatShapeFactor_[0] : throatShapeFactor_[eIdx]; } |
| 298 | |||
| 299 | //! Returns the vector of throat shape factors | ||
| 300 | 2 | const std::vector<Scalar>& throatShapeFactor() const | |
| 301 | { | ||
| 302 | 2 | if (useSameShapeForAllThroats()) | |
| 303 | { | ||
| 304 | // if a vector of throat shape factors is requested (e.g., for vtk output), | ||
| 305 | // resize the container and fill it with the same value everywhere | ||
| 306 | 2 | const auto shapeFactor = throatShapeFactor_[0]; | |
| 307 | 2 | throatShapeFactor_.resize(throatInscribedRadius_.size(), shapeFactor); | |
| 308 | } | ||
| 309 | |||
| 310 | 2 | return throatShapeFactor_; | |
| 311 | } | ||
| 312 | |||
| 313 | //! Returns whether all pores feature the same shape | ||
| 314 | 6145299 | bool useSameGeometryForAllPores() const | |
| 315 |
7/10✗ Branch 0 not taken.
✓ Branch 1 taken 18855 times.
✓ Branch 2 taken 83320 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2979476 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 83320 times.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 2979476 times.
✗ Branch 9 not taken.
|
6144451 | { return useSameGeometryForAllPores_; } |
| 316 | |||
| 317 | //! Returns whether all throats feature the same cross-sectional shape | ||
| 318 | 45482220 | bool useSameShapeForAllThroats() const | |
| 319 |
11/18✓ Branch 0 taken 61536 times.
✓ Branch 1 taken 35432943 times.
✓ Branch 2 taken 418355 times.
✓ Branch 3 taken 30768 times.
✓ Branch 4 taken 2522980 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1092264 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 92459 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 2879460 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 959946 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 960111 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 959820 times.
✗ Branch 17 not taken.
|
45410642 | { return useSameShapeForAllThroats_; } |
| 320 | |||
| 321 | private: | ||
| 322 | |||
| 323 | //! determine the pore geometry provided as scalar value by the grid file | ||
| 324 | template<class GridData> | ||
| 325 | ✗ | Pore::Shape getPoreGeometry_(const GridData& gridData, const Vertex& vertex) const | |
| 326 | { | ||
| 327 | ✗ | static const auto poreGeometryIdx = gridData.parameterIndex("PoreGeometry"); | |
| 328 | using T = std::underlying_type_t<Pore::Shape>; | ||
| 329 | ✗ | const auto poreGeometryValue = static_cast<T>(gridData.parameters(vertex)[poreGeometryIdx]); | |
| 330 | ✗ | return static_cast<Pore::Shape>(poreGeometryValue); | |
| 331 | } | ||
| 332 | |||
| 333 | //! automatically determine the pore volume if not provided by the grid file | ||
| 334 | template<class GridData> | ||
| 335 | 18855 | Scalar getPoreVolume_(const GridData& gridData, const Vertex& vertex, const std::size_t vIdx) const | |
| 336 | { | ||
| 337 |
4/6✓ Branch 0 taken 25 times.
✓ Branch 1 taken 18830 times.
✓ Branch 3 taken 25 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 25 times.
✗ Branch 7 not taken.
|
18855 | static const bool gridHasPoreVolume = gridData.gridHasVertexParameter("PoreVolume"); |
| 338 | |||
| 339 |
2/2✓ Branch 0 taken 18035 times.
✓ Branch 1 taken 820 times.
|
18855 | if (gridHasPoreVolume) |
| 340 | { | ||
| 341 |
5/8✓ Branch 0 taken 14 times.
✓ Branch 1 taken 18021 times.
✓ Branch 3 taken 14 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 14 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 12 times.
✗ Branch 10 not taken.
|
18037 | static const auto poreVolumeIdx = gridData.parameterIndex("PoreVolume"); |
| 342 | 18035 | return gridData.parameters(vertex)[poreVolumeIdx]; | |
| 343 | } | ||
| 344 | else | ||
| 345 | { | ||
| 346 |
2/4✓ Branch 0 taken 820 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 820 times.
|
1640 | if (poreGeometry(vIdx) == Pore::Shape::cylinder) |
| 347 | { | ||
| 348 | ✗ | static const Scalar fixedHeight = getParamFromGroup<Scalar>(gridData.paramGroup(), "Grid.PoreHeight", -1.0); | |
| 349 | ✗ | const Scalar h = fixedHeight > 0.0 ? fixedHeight : gridData.getParameter(vertex, "PoreHeight"); | |
| 350 | ✗ | return Pore::volume(Pore::Shape::cylinder, poreInscribedRadius(vIdx), h); | |
| 351 | } | ||
| 352 | else | ||
| 353 | 820 | return Pore::volume(poreGeometry(vIdx), poreInscribedRadius(vIdx)); | |
| 354 | } | ||
| 355 | } | ||
| 356 | |||
| 357 | //! automatically determine throat cross-sectional area if not provided by the grid file | ||
| 358 | template<class GridData> | ||
| 359 | 35775 | Scalar getThroatCrossSectionalArea_(const GridData& gridData, const Element& element, const std::size_t eIdx) const | |
| 360 | { | ||
| 361 |
4/6✓ Branch 0 taken 24 times.
✓ Branch 1 taken 35751 times.
✓ Branch 3 taken 24 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 24 times.
✗ Branch 7 not taken.
|
35775 | static const bool gridHasThroatCrossSectionalArea = gridData.gridHasElementParameter("ThroatCrossSectionalArea"); |
| 362 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 35775 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
35775 | if (gridHasThroatCrossSectionalArea && !overwriteGridDataWithShapeSpecificValues_) |
| 363 | { | ||
| 364 | ✗ | static const auto throatAreaIdx = gridData.parameterIndex("ThroatCrossSectionalArea"); | |
| 365 | ✗ | return gridData.parameters(element)[throatAreaIdx]; | |
| 366 | } | ||
| 367 | else | ||
| 368 | { | ||
| 369 |
2/4✓ Branch 0 taken 35775 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 35775 times.
|
71550 | if (const auto shape = throatCrossSectionShape(eIdx); shape == Throat::Shape::rectangle) |
| 370 | { | ||
| 371 | ✗ | static const auto throatHeight = getParamFromGroup<Scalar>(gridData.paramGroup(), "Grid.ThroatHeight"); | |
| 372 | ✗ | return Throat::totalCrossSectionalAreaForRectangle(throatInscribedRadius_[eIdx], throatHeight); | |
| 373 | } | ||
| 374 | else | ||
| 375 | 35775 | return Throat::totalCrossSectionalArea(shape, throatInscribedRadius_[eIdx]); | |
| 376 | } | ||
| 377 | } | ||
| 378 | |||
| 379 | //! automatically determine throat shape factor if not provided by the grid file | ||
| 380 | template<class GridData> | ||
| 381 | 35775 | Scalar getThroatShapeFactor_(const GridData& gridData, const Element& element, const std::size_t eIdx) const | |
| 382 | { | ||
| 383 |
4/6✓ Branch 0 taken 24 times.
✓ Branch 1 taken 35751 times.
✓ Branch 3 taken 24 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 24 times.
✗ Branch 7 not taken.
|
35775 | static const bool gridHasThroatShapeFactor = gridData.gridHasElementParameter("ThroatShapeFactor"); |
| 384 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 35775 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
35775 | if (gridHasThroatShapeFactor && !overwriteGridDataWithShapeSpecificValues_) |
| 385 | { | ||
| 386 | ✗ | static const auto throatShapeFactorIdx = gridData.parameterIndex("ThroatShapeFactor"); | |
| 387 | ✗ | return gridData.parameters(element)[throatShapeFactorIdx]; | |
| 388 | } | ||
| 389 | else | ||
| 390 | { | ||
| 391 |
2/4✓ Branch 0 taken 35775 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 35775 times.
|
71550 | if (const auto shape = throatCrossSectionShape(eIdx); shape == Throat::Shape::rectangle) |
| 392 | { | ||
| 393 | ✗ | static const auto throatHeight = getParamFromGroup<Scalar>(gridData.paramGroup(), "Grid.ThroatHeight"); | |
| 394 | ✗ | return Throat::shapeFactorRectangle(throatInscribedRadius_[eIdx], throatHeight); | |
| 395 | } | ||
| 396 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 35775 times.
|
35775 | else if (shape == Throat::Shape::polygon || shape == Throat::Shape::scaleneTriangle) |
| 397 | { | ||
| 398 | ✗ | static const auto shapeFactor = getParamFromGroup<Scalar>(gridData.paramGroup(), "Grid.ThroatShapeFactor"); | |
| 399 | ✗ | return shapeFactor; | |
| 400 | } | ||
| 401 | else | ||
| 402 | 35775 | return Throat::shapeFactor<Scalar>(shape, throatInscribedRadius_[eIdx]); | |
| 403 | } | ||
| 404 | } | ||
| 405 | |||
| 406 | 28 | void maybeResizeContainers_() | |
| 407 | { | ||
| 408 | // check if all throat might have the same shape in order to save some memory | ||
| 409 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 24 times.
|
28 | if (!useSameShapeForAllThroats() && |
| 410 |
2/4✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
|
8 | std::adjacent_find(throatGeometry_.begin(), throatGeometry_.end(), std::not_equal_to<Throat::Shape>() ) == throatGeometry_.end()) |
| 411 | { | ||
| 412 | ✗ | std::cout << "All throats feature the same shape, resizing containers" << std::endl; | |
| 413 | ✗ | useSameShapeForAllThroats_ = true; | |
| 414 | ✗ | const Scalar shapeFactor = throatShapeFactor_[0]; | |
| 415 | ✗ | const auto throatGeometry = throatGeometry_[0]; | |
| 416 | ✗ | throatShapeFactor_.resize(1); | |
| 417 | ✗ | throatGeometry_.resize(1); | |
| 418 | ✗ | throatShapeFactor_[0] = shapeFactor; | |
| 419 | ✗ | throatGeometry_[0] = throatGeometry; | |
| 420 | } | ||
| 421 | |||
| 422 | // check if all throat might have the same shape in order to save some memory | ||
| 423 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 28 times.
|
28 | if (!useSameGeometryForAllPores() && |
| 424 | ✗ | std::adjacent_find(poreGeometry_.begin(), poreGeometry_.end(), std::not_equal_to<Pore::Shape>() ) == poreGeometry_.end()) | |
| 425 | { | ||
| 426 | ✗ | std::cout << "All pores feature the same shape, resizing containers" << std::endl; | |
| 427 | ✗ | useSameGeometryForAllPores_ = true; | |
| 428 | ✗ | const auto poreGeometry = poreGeometry_[0]; | |
| 429 | ✗ | poreGeometry_.resize(1); | |
| 430 | ✗ | poreGeometry_[0] = poreGeometry; | |
| 431 | } | ||
| 432 | 28 | } | |
| 433 | |||
| 434 | mutable std::vector<Pore::Shape> poreGeometry_; | ||
| 435 | std::vector<Scalar> poreInscribedRadius_; | ||
| 436 | std::vector<Scalar> poreVolume_; | ||
| 437 | std::vector<Label> poreLabel_; // 0:no, 1:general, 2:coupling1, 3:coupling2, 4:inlet, 5:outlet | ||
| 438 | std::vector<SmallLocalIndex> coordinationNumber_; | ||
| 439 | mutable std::vector<Throat::Shape> throatGeometry_; | ||
| 440 | mutable std::vector<Scalar> throatShapeFactor_; | ||
| 441 | std::vector<Scalar> throatInscribedRadius_; | ||
| 442 | std::vector<Scalar> throatLength_; | ||
| 443 | std::vector<Label> throatLabel_; // 0:no, 1:general, 2:coupling1, 3:coupling2, 4:inlet, 5:outlet | ||
| 444 | std::vector<Scalar> throatCrossSectionalArea_; | ||
| 445 | bool useSameGeometryForAllPores_; | ||
| 446 | bool useSameShapeForAllThroats_; | ||
| 447 | bool overwriteGridDataWithShapeSpecificValues_; | ||
| 448 | }; | ||
| 449 | |||
| 450 | /*! | ||
| 451 | * \ingroup PoreNetworkDiscretization | ||
| 452 | * \brief The default traits | ||
| 453 | * \tparam the grid view type | ||
| 454 | */ | ||
| 455 | template<class GridView, class MapperTraits = DefaultMapperTraits<GridView>> | ||
| 456 | struct PNMDefaultGridGeometryTraits | ||
| 457 | : public MapperTraits | ||
| 458 | { | ||
| 459 | using SubControlVolume = PNMSubControlVolume<GridView>; | ||
| 460 | using SubControlVolumeFace = PNMSubControlVolumeFace<GridView>; | ||
| 461 | |||
| 462 | template<class GridGeometry, bool enableCache> | ||
| 463 | using LocalView = PNMFVElementGeometry<GridGeometry, enableCache>; | ||
| 464 | |||
| 465 | using PNMData = DefaultPNMData<typename SubControlVolume::Traits::Scalar, GridView>; | ||
| 466 | }; | ||
| 467 | |||
| 468 | /*! | ||
| 469 | * \ingroup PoreNetworkDiscretization | ||
| 470 | * \brief Base class for the finite volume geometry for porenetwork models | ||
| 471 | * \note This class is specialized for versions with and without caching the fv geometries on the grid view | ||
| 472 | */ | ||
| 473 | template<class Scalar, | ||
| 474 | class GridView, | ||
| 475 | bool enableGridGeometryCache = false, | ||
| 476 | class Traits = PNMDefaultGridGeometryTraits<GridView> > | ||
| 477 | class GridGeometry; | ||
| 478 | |||
| 479 | /*! | ||
| 480 | * \ingroup PoreNetworkDiscretization | ||
| 481 | * \brief Base class for the finite volume geometry for porenetwork models | ||
| 482 | * \note For caching enabled we store the fv geometries for the whole grid view which is memory intensive but faster | ||
| 483 | */ | ||
| 484 | template<class Scalar, class GV, class Traits> | ||
| 485 | class GridGeometry<Scalar, GV, true, Traits> | ||
| 486 | : public BaseGridGeometry<GV, Traits> | ||
| 487 | , public Traits::PNMData | ||
| 488 | { | ||
| 489 | using ThisType = GridGeometry<Scalar, GV, true, Traits>; | ||
| 490 | using ParentType = BaseGridGeometry<GV, Traits>; | ||
| 491 | using GridIndexType = typename IndexTraits<GV>::GridIndex; | ||
| 492 | using LocalIndexType = typename IndexTraits<GV>::LocalIndex; | ||
| 493 | using PNMData = typename Traits::PNMData; | ||
| 494 | |||
| 495 | using Element = typename GV::template Codim<0>::Entity; | ||
| 496 | using CoordScalar = typename GV::ctype; | ||
| 497 | static const int dim = GV::dimension; | ||
| 498 | static const int dimWorld = GV::dimensionworld; | ||
| 499 | |||
| 500 | public: | ||
| 501 | //! export the discretization method this geometry belongs to | ||
| 502 | using DiscretizationMethod = DiscretizationMethods::Box; | ||
| 503 | static constexpr DiscretizationMethod discMethod{}; | ||
| 504 | |||
| 505 | //! export the type of the fv element geometry (the local view type) | ||
| 506 | using LocalView = typename Traits::template LocalView<ThisType, true>; | ||
| 507 | //! export the type of sub control volume | ||
| 508 | using SubControlVolume = typename Traits::SubControlVolume; | ||
| 509 | //! export the type of sub control volume | ||
| 510 | using SubControlVolumeFace = typename Traits::SubControlVolumeFace; | ||
| 511 | //! export the type of extrusion | ||
| 512 | using Extrusion = Extrusion_t<Traits>; | ||
| 513 | //! export dof mapper type | ||
| 514 | using DofMapper = typename Traits::VertexMapper; | ||
| 515 | //! export the finite element cache type | ||
| 516 | using FeCache = Dune::LagrangeLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>; | ||
| 517 | //! export the grid view type | ||
| 518 | using GridView = GV; | ||
| 519 | |||
| 520 | //! Constructor | ||
| 521 | template<class GridData> | ||
| 522 | GridGeometry(const GridView& gridView, const GridData& gridData) | ||
| 523 | : ParentType(gridView) | ||
| 524 | { | ||
| 525 | static_assert(GridView::dimension == 1, "Porenetwork model only allow GridView::dimension == 1!"); | ||
| 526 | update_(gridData); | ||
| 527 | } | ||
| 528 | |||
| 529 | //! the vertex mapper is the dofMapper | ||
| 530 | //! this is convenience to have better chance to have the same main files for box/tpfa/mpfa... | ||
| 531 | const DofMapper& dofMapper() const | ||
| 532 | { return this->vertexMapper(); } | ||
| 533 | |||
| 534 | //! The total number of sub control volumes | ||
| 535 | std::size_t numScv() const | ||
| 536 | { return numScv_; } | ||
| 537 | |||
| 538 | //! The total number of sun control volume faces | ||
| 539 | std::size_t numScvf() const | ||
| 540 | { return numScvf_; } | ||
| 541 | |||
| 542 | //! The total number of boundary sub control volume faces | ||
| 543 | //! For compatibility reasons with cc methods | ||
| 544 | std::size_t numBoundaryScvf() const | ||
| 545 | { return 0; } | ||
| 546 | |||
| 547 | //! The total number of degrees of freedom | ||
| 548 | std::size_t numDofs() const | ||
| 549 | { return this->vertexMapper().size(); } | ||
| 550 | |||
| 551 | //! update all fvElementGeometries (call this after grid adaption) | ||
| 552 | template<class GridData> | ||
| 553 | void update(const GridView& gridView, const GridData& gridData) | ||
| 554 | { | ||
| 555 | ParentType::update(gridView); | ||
| 556 | update_(gridData); | ||
| 557 | } | ||
| 558 | |||
| 559 | //! update all fvElementGeometries (call this after grid adaption) | ||
| 560 | template<class GridData> | ||
| 561 | void update(GridView&& gridView, const GridData& gridData) | ||
| 562 | { | ||
| 563 | ParentType::update(std::move(gridView)); | ||
| 564 | update_(gridData); | ||
| 565 | } | ||
| 566 | |||
| 567 | //! The finite element cache for creating local FE bases | ||
| 568 | const FeCache& feCache() const | ||
| 569 | { return feCache_; } | ||
| 570 | |||
| 571 | //! Get the local scvs for an element | ||
| 572 | const std::array<SubControlVolume, 2>& scvs(GridIndexType eIdx) const | ||
| 573 | { return scvs_[eIdx]; } | ||
| 574 | |||
| 575 | //! Get the local scvfs for an element | ||
| 576 | const std::array<SubControlVolumeFace, 1>& scvfs(GridIndexType eIdx) const | ||
| 577 | { return scvfs_[eIdx]; } | ||
| 578 | |||
| 579 | //! If a vertex / d.o.f. is on the boundary | ||
| 580 | bool dofOnBoundary(GridIndexType dofIdx) const | ||
| 581 | { return boundaryDofIndices_[dofIdx]; } | ||
| 582 | |||
| 583 | //! If a vertex / d.o.f. is on a periodic boundary (not implemented) | ||
| 584 | bool dofOnPeriodicBoundary(GridIndexType dofIdx) const | ||
| 585 | { return false; } | ||
| 586 | |||
| 587 | //! The index of the vertex / d.o.f. on the other side of the periodic boundary | ||
| 588 | GridIndexType periodicallyMappedDof(GridIndexType dofIdx) const | ||
| 589 | { DUNE_THROW(Dune::NotImplemented, "Periodic boundaries"); } | ||
| 590 | |||
| 591 | //! Returns whether one of the geometry's scvfs lies on a boundary | ||
| 592 | bool hasBoundaryScvf(GridIndexType eIdx) const | ||
| 593 | { return hasBoundaryScvf_[eIdx]; } | ||
| 594 | |||
| 595 | private: | ||
| 596 | |||
| 597 | template<class GridData> | ||
| 598 | void update_(const GridData& gridData) | ||
| 599 | { | ||
| 600 | PNMData::update(this->gridView(), gridData); | ||
| 601 | |||
| 602 | scvs_.clear(); | ||
| 603 | scvfs_.clear(); | ||
| 604 | |||
| 605 | auto numElements = this->gridView().size(0); | ||
| 606 | scvs_.resize(numElements); | ||
| 607 | scvfs_.resize(numElements); | ||
| 608 | hasBoundaryScvf_.resize(numElements, false); | ||
| 609 | |||
| 610 | boundaryDofIndices_.assign(numDofs(), false); | ||
| 611 | |||
| 612 | numScvf_ = numElements; | ||
| 613 | numScv_ = 2*numElements; | ||
| 614 | |||
| 615 | // Build the SCV and SCV faces | ||
| 616 | for (const auto& element : elements(this->gridView())) | ||
| 617 | { | ||
| 618 | // get the element geometry | ||
| 619 | auto eIdx = this->elementMapper().index(element); | ||
| 620 | auto elementGeometry = element.geometry(); | ||
| 621 | |||
| 622 | // construct the sub control volumes | ||
| 623 | for (LocalIndexType scvLocalIdx = 0; scvLocalIdx < elementGeometry.corners(); ++scvLocalIdx) | ||
| 624 | { | ||
| 625 | const auto dofIdxGlobal = this->vertexMapper().subIndex(element, scvLocalIdx, dim); | ||
| 626 | |||
| 627 | // get the corners | ||
| 628 | auto corners = std::array{elementGeometry.corner(scvLocalIdx), elementGeometry.center()}; | ||
| 629 | |||
| 630 | // get the fractional volume associated with the scv | ||
| 631 | const auto volume = this->poreVolume(dofIdxGlobal) / this->coordinationNumber(dofIdxGlobal); | ||
| 632 | |||
| 633 | scvs_[eIdx][scvLocalIdx] = SubControlVolume(dofIdxGlobal, | ||
| 634 | scvLocalIdx, | ||
| 635 | eIdx, | ||
| 636 | std::move(corners), | ||
| 637 | volume); | ||
| 638 | |||
| 639 | if (this->poreLabel(dofIdxGlobal) > 0) | ||
| 640 | { | ||
| 641 | if (boundaryDofIndices_[dofIdxGlobal]) | ||
| 642 | continue; | ||
| 643 | |||
| 644 | boundaryDofIndices_[dofIdxGlobal] = true; | ||
| 645 | hasBoundaryScvf_[eIdx] = true; | ||
| 646 | } | ||
| 647 | } | ||
| 648 | |||
| 649 | // construct the inner sub control volume face | ||
| 650 | auto unitOuterNormal = elementGeometry.corner(1)-elementGeometry.corner(0); | ||
| 651 | unitOuterNormal /= unitOuterNormal.two_norm(); | ||
| 652 | LocalIndexType scvfLocalIdx = 0; | ||
| 653 | scvfs_[eIdx][0] = SubControlVolumeFace(elementGeometry.center(), | ||
| 654 | std::move(unitOuterNormal), | ||
| 655 | this->throatCrossSectionalArea(this->elementMapper().index(element)), | ||
| 656 | scvfLocalIdx++, | ||
| 657 | std::array<LocalIndexType, 2>({0, 1})); | ||
| 658 | } | ||
| 659 | } | ||
| 660 | |||
| 661 | const FeCache feCache_; | ||
| 662 | |||
| 663 | std::vector<std::array<SubControlVolume, 2>> scvs_; | ||
| 664 | std::vector<std::array<SubControlVolumeFace, 1>> scvfs_; | ||
| 665 | |||
| 666 | std::size_t numScv_; | ||
| 667 | std::size_t numScvf_; | ||
| 668 | |||
| 669 | // vertices on the boundary | ||
| 670 | std::vector<bool> boundaryDofIndices_; | ||
| 671 | std::vector<bool> hasBoundaryScvf_; | ||
| 672 | }; | ||
| 673 | |||
| 674 | /*! | ||
| 675 | * \ingroup PoreNetworkDiscretization | ||
| 676 | * \brief Base class for the finite volume geometry for porenetwork models | ||
| 677 | * \note For caching disabled we store only some essential index maps to build up local systems on-demand in | ||
| 678 | * the corresponding FVElementGeometry | ||
| 679 | */ | ||
| 680 | template<class Scalar, class GV, class Traits> | ||
| 681 | class GridGeometry<Scalar, GV, false, Traits> | ||
| 682 | : public BaseGridGeometry<GV, Traits> | ||
| 683 | , public Traits::PNMData | ||
| 684 | { | ||
| 685 | using ThisType = GridGeometry<Scalar, GV, false, Traits>; | ||
| 686 | using ParentType = BaseGridGeometry<GV, Traits>; | ||
| 687 | using GridIndexType = typename IndexTraits<GV>::GridIndex; | ||
| 688 | using LocalIndexType = typename IndexTraits<GV>::LocalIndex; | ||
| 689 | using PNMData = typename Traits::PNMData; | ||
| 690 | |||
| 691 | static const int dim = GV::dimension; | ||
| 692 | static const int dimWorld = GV::dimensionworld; | ||
| 693 | |||
| 694 | using Element = typename GV::template Codim<0>::Entity; | ||
| 695 | using CoordScalar = typename GV::ctype; | ||
| 696 | |||
| 697 | public: | ||
| 698 | //! export the discretization method this geometry belongs to | ||
| 699 | using DiscretizationMethod = DiscretizationMethods::Box; | ||
| 700 | static constexpr DiscretizationMethod discMethod{}; | ||
| 701 | |||
| 702 | //! export the type of the fv element geometry (the local view type) | ||
| 703 | using LocalView = typename Traits::template LocalView<ThisType, false>; | ||
| 704 | //! export the type of sub control volume | ||
| 705 | using SubControlVolume = typename Traits::SubControlVolume; | ||
| 706 | //! export the type of sub control volume | ||
| 707 | using SubControlVolumeFace = typename Traits::SubControlVolumeFace; | ||
| 708 | //! export the type of extrusion | ||
| 709 | using Extrusion = Extrusion_t<Traits>; | ||
| 710 | //! export dof mapper type | ||
| 711 | using DofMapper = typename Traits::VertexMapper; | ||
| 712 | //! export the finite element cache type | ||
| 713 | using FeCache = Dune::LagrangeLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>; | ||
| 714 | //! export the grid view type | ||
| 715 | using GridView = GV; | ||
| 716 | |||
| 717 | //! Constructor | ||
| 718 | template<class GridData> | ||
| 719 | 26 | GridGeometry(const GridView& gridView, const GridData& gridData) | |
| 720 |
2/4✓ Branch 3 taken 26 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 26 times.
✗ Branch 7 not taken.
|
26 | : ParentType(gridView) |
| 721 | { | ||
| 722 | static_assert(GridView::dimension == 1, "Porenetwork model only allow GridView::dimension == 1!"); | ||
| 723 |
1/2✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
|
26 | update_(gridData); |
| 724 | 26 | } | |
| 725 | |||
| 726 | |||
| 727 | //! the vertex mapper is the dofMapper | ||
| 728 | //! this is convenience to have better chance to have the same main files for box/tpfa/mpfa... | ||
| 729 | 10914175 | const DofMapper& dofMapper() const | |
| 730 |
4/16✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 571 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 5405 times.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 5 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 6 times.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
|
10914749 | { return this->vertexMapper(); } |
| 731 | |||
| 732 | //! The total number of sub control volumes | ||
| 733 | std::size_t numScv() const | ||
| 734 | { return numScv_; } | ||
| 735 | |||
| 736 | //! The total number of sun control volume faces | ||
| 737 | std::size_t numScvf() const | ||
| 738 | { return numScvf_; } | ||
| 739 | |||
| 740 | //! The total number of boundary sub control volume faces | ||
| 741 | //! For compatibility reasons with cc methods | ||
| 742 | std::size_t numBoundaryScvf() const | ||
| 743 | { return 0; } | ||
| 744 | |||
| 745 | //! The total number of degrees of freedom | ||
| 746 | 434 | std::size_t numDofs() const | |
| 747 |
19/31✓ Branch 1 taken 3 times.
✓ Branch 2 taken 6 times.
✓ Branch 4 taken 5 times.
✓ Branch 5 taken 9 times.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✓ Branch 14 taken 1 times.
✓ Branch 16 taken 7 times.
✓ Branch 17 taken 1 times.
✓ Branch 19 taken 3 times.
✓ Branch 20 taken 1 times.
✓ Branch 22 taken 1 times.
✓ Branch 23 taken 1 times.
✓ Branch 25 taken 1 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 1 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 1 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 1 times.
✗ Branch 35 not taken.
✗ Branch 15 not taken.
✗ Branch 18 not taken.
✗ Branch 21 not taken.
✗ Branch 24 not taken.
✗ Branch 3 not taken.
✗ Branch 6 not taken.
✓ Branch 8 taken 6 times.
✗ Branch 9 not taken.
✓ Branch 12 taken 280 times.
|
403 | { return this->vertexMapper().size(); } |
| 748 | |||
| 749 | //! update all fvElementGeometries (call this after grid adaption) | ||
| 750 | template<class GridData> | ||
| 751 | 2 | void update(const GridView& gridView, const GridData& gridData) | |
| 752 | { | ||
| 753 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | ParentType::update(gridView); |
| 754 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | update_(gridData); |
| 755 | 2 | } | |
| 756 | |||
| 757 | //! update all fvElementGeometries (call this after grid adaption) | ||
| 758 | template<class GridData> | ||
| 759 | void update(GridView&& gridView, const GridData& gridData) | ||
| 760 | { | ||
| 761 | ParentType::update(std::move(gridView)); | ||
| 762 | update_(gridData); | ||
| 763 | } | ||
| 764 | |||
| 765 | //! The finite element cache for creating local FE bases | ||
| 766 | 5454381 | const FeCache& feCache() const | |
| 767 | 5454381 | { return feCache_; } | |
| 768 | |||
| 769 | //! If a vertex / d.o.f. is on the boundary | ||
| 770 | 3297598 | bool dofOnBoundary(GridIndexType dofIdx) const | |
| 771 |
6/6✓ Branch 0 taken 222918 times.
✓ Branch 1 taken 2132846 times.
✓ Branch 2 taken 85041 times.
✓ Branch 3 taken 855817 times.
✓ Branch 4 taken 493 times.
✓ Branch 5 taken 483 times.
|
3297598 | { return boundaryDofIndices_[dofIdx]; } |
| 772 | |||
| 773 | //! If a vertex / d.o.f. is on a periodic boundary (not implemented) | ||
| 774 | bool dofOnPeriodicBoundary(GridIndexType dofIdx) const | ||
| 775 | { return false; } | ||
| 776 | |||
| 777 | //! The index of the vertex / d.o.f. on the other side of the periodic boundary | ||
| 778 | GridIndexType periodicallyMappedDof(GridIndexType dofIdx) const | ||
| 779 | { DUNE_THROW(Dune::NotImplemented, "Periodic boundaries"); } | ||
| 780 | |||
| 781 | private: | ||
| 782 | |||
| 783 | template<class GridData> | ||
| 784 | 28 | void update_(const GridData& gridData) | |
| 785 | { | ||
| 786 | 28 | PNMData::update(this->gridView(), gridData); | |
| 787 | |||
| 788 | 28 | boundaryDofIndices_.assign(numDofs(), false); | |
| 789 | |||
| 790 | // save global data on the grid's scvs and scvfs | ||
| 791 | 28 | numScvf_ = this->gridView().size(0); | |
| 792 | 28 | numScv_ = 2*numScvf_; | |
| 793 | |||
| 794 |
2/2✓ Branch 2 taken 42157 times.
✓ Branch 3 taken 28 times.
|
84342 | for (const auto& element : elements(this->gridView())) |
| 795 | { | ||
| 796 | // treat boundaries | ||
| 797 |
2/2✓ Branch 0 taken 84314 times.
✓ Branch 1 taken 42157 times.
|
126471 | for (LocalIndexType vIdxLocal = 0; vIdxLocal < 2; ++vIdxLocal) |
| 798 | { | ||
| 799 | 84314 | const auto vIdxGlobal = this->vertexMapper().subIndex(element, vIdxLocal, dim); | |
| 800 |
2/2✓ Branch 0 taken 14159 times.
✓ Branch 1 taken 70155 times.
|
84314 | if (this->poreLabel(vIdxGlobal) > 0) |
| 801 | { | ||
| 802 |
2/2✓ Branch 0 taken 9390 times.
✓ Branch 1 taken 4769 times.
|
14159 | if (boundaryDofIndices_[vIdxGlobal]) |
| 803 | 9390 | continue; | |
| 804 | |||
| 805 | 4769 | boundaryDofIndices_[vIdxGlobal] = true; | |
| 806 | } | ||
| 807 | } | ||
| 808 | } | ||
| 809 | 28 | } | |
| 810 | |||
| 811 | const FeCache feCache_; | ||
| 812 | |||
| 813 | // Information on the global number of geometries | ||
| 814 | std::size_t numScv_; | ||
| 815 | std::size_t numScvf_; | ||
| 816 | |||
| 817 | // vertices on the boundary | ||
| 818 | std::vector<bool> boundaryDofIndices_; | ||
| 819 | }; | ||
| 820 | |||
| 821 | } // end namespace Dumux::PoreNetwork | ||
| 822 | |||
| 823 | #endif | ||
| 824 |