| 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 FreeFlowPoreNetworkCoupling | ||
| 10 | * \brief A grid creator that matches a free-flow grid to a PNM grid. | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_MULTIDOMAIN_BOUNDARY_FREEFLOW_PORENETWORK_SNAPPY_GRID_MANAGER_HH | ||
| 13 | #define DUMUX_MULTIDOMAIN_BOUNDARY_FREEFLOW_PORENETWORK_SNAPPY_GRID_MANAGER_HH | ||
| 14 | |||
| 15 | #include <bitset> | ||
| 16 | #include <optional> | ||
| 17 | #include <dune/common/float_cmp.hh> | ||
| 18 | #include <dune/geometry/axisalignedcubegeometry.hh> | ||
| 19 | #include <dumux/geometry/intersectspointgeometry.hh> | ||
| 20 | #include <dumux/io/grid/gridmanager.hh> | ||
| 21 | |||
| 22 | namespace Dumux::PoreNetwork { | ||
| 23 | |||
| 24 | namespace Detail { | ||
| 25 | /*! | ||
| 26 | * \ingroup FreeFlowPoreNetworkCoupling | ||
| 27 | * \brief A helper for the grid creator that matches a free-flow grid to a PNM grid. | ||
| 28 | */ | ||
| 29 | template<class GridView> | ||
| 30 | class SnappyGridManagerHelper | ||
| 31 | { | ||
| 32 | using Scalar = typename GridView::ctype; | ||
| 33 | static constexpr auto dim = GridView::dimension; | ||
| 34 | static constexpr auto dimWorld = GridView::dimensionworld; | ||
| 35 | using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>; | ||
| 36 | using Line = Dune::AxisAlignedCubeGeometry< Scalar, 1, dimWorld>; | ||
| 37 | using Plane = Dune::AxisAlignedCubeGeometry< Scalar, dimWorld-1, dimWorld>; | ||
| 38 | |||
| 39 | public: | ||
| 40 | |||
| 41 | /*! | ||
| 42 | * \brief Returns the direction index of a unit vector (0 = x, 1 = y, 2 = z) | ||
| 43 | */ | ||
| 44 | 12 | static unsigned int directionIndex(const GlobalPosition& vector) | |
| 45 | { | ||
| 46 | 12 | const auto eps = 1e-8*vector.two_norm(); | |
| 47 |
7/30✗ 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.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✓ Branch 20 taken 2 times.
✓ Branch 21 taken 10 times.
✓ Branch 22 taken 2 times.
✓ Branch 23 taken 10 times.
✓ Branch 24 taken 10 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 10 times.
✗ Branch 27 not taken.
✗ Branch 32 not taken.
✓ Branch 33 taken 12 times.
|
92 | return std::find_if(vector.begin(), vector.end(), [eps](const auto& x) { return std::abs(x) > eps; } ) - vector.begin(); |
| 48 | } | ||
| 49 | |||
| 50 | 4 | static auto couplingPlaneBoundingBox(const GlobalPosition& gridLowerLeft, | |
| 51 | const GlobalPosition& gridUpperRight, | ||
| 52 | const GlobalPosition& couplingPlaneNormal, | ||
| 53 | const std::string& modelParamGroup) | ||
| 54 | { | ||
| 55 | 4 | const auto couplingPlaneNormalDirectionIndex = directionIndex(couplingPlaneNormal); | |
| 56 | |||
| 57 | // get the spatial extent of the coupling plane | ||
| 58 | 4 | const auto couplingPlaneLowerLeft = getParamFromGroup<GlobalPosition>(modelParamGroup, "Grid.CouplingPlaneLowerLeft", gridLowerLeft); | |
| 59 | 4 | const auto couplingPlaneUpperRight = getParamFromGroup<GlobalPosition>(modelParamGroup, "Grid.CouplingPlaneUpperRight", | |
| 60 | [&]() | ||
| 61 | { | ||
| 62 | 4 | GlobalPosition tmp(gridUpperRight); | |
| 63 | 4 | tmp[couplingPlaneNormalDirectionIndex] = gridLowerLeft[couplingPlaneNormalDirectionIndex]; | |
| 64 | 4 | return tmp; | |
| 65 | 8 | }()); | |
| 66 | 8 | return std::make_pair(couplingPlaneLowerLeft, couplingPlaneUpperRight); | |
| 67 | } | ||
| 68 | |||
| 69 | 2 | static Plane makeCouplingPlane(const GlobalPosition& gridLowerLeft, | |
| 70 | const GlobalPosition& gridUpperRight, | ||
| 71 | const GlobalPosition& couplingPlaneNormal, | ||
| 72 | const std::string& modelParamGroup) | ||
| 73 | { | ||
| 74 | 2 | auto [couplingPlaneLowerLeft, couplingPlaneUpperRight] = couplingPlaneBoundingBox(gridLowerLeft, gridUpperRight, couplingPlaneNormal, modelParamGroup); | |
| 75 | 2 | const auto couplingPlaneNormalDirectionIndex = directionIndex(couplingPlaneNormal); | |
| 76 | 4 | auto inPlaneAxes = std::move(std::bitset<dimWorld>{}.set()); | |
| 77 | 2 | inPlaneAxes.set(couplingPlaneNormalDirectionIndex, false); | |
| 78 | 2 | return Plane(std::move(couplingPlaneLowerLeft), std::move(couplingPlaneUpperRight), inPlaneAxes); | |
| 79 | } | ||
| 80 | |||
| 81 | /*! | ||
| 82 | * \brief Creates lines parallel to the bounding box axis of the coupling plane. | ||
| 83 | * Creates one line for dim == 2 and two lines for dim == 3. The array size is dim, therefore the entry | ||
| 84 | * for the normal to the coupling plane is empty. | ||
| 85 | * | ||
| 86 | * \param gridLowerLeft The lower left position of the bulk grid | ||
| 87 | * \param gridUpperRight The upper right position of the bulk grid | ||
| 88 | * \param couplingPlaneNormal The normal vector of the coupling plane | ||
| 89 | * \param modelParamGroup The bulk model parameter group | ||
| 90 | */ | ||
| 91 | 2 | static std::array<std::optional<Line>, dim> makeAxisParallelLinesFromCouplingPlane(const GlobalPosition& gridLowerLeft, | |
| 92 | const GlobalPosition& gridUpperRight, | ||
| 93 | const GlobalPosition& couplingPlaneNormal, | ||
| 94 | const std::string& modelParamGroup) | ||
| 95 | { | ||
| 96 | 2 | const auto couplingPlaneNormalDirectionIndex = directionIndex(couplingPlaneNormal); | |
| 97 | |||
| 98 | // get the spatial extent of the coupling plane | ||
| 99 | 2 | auto [a, b] = couplingPlaneBoundingBox(gridLowerLeft, gridUpperRight, couplingPlaneNormal, modelParamGroup); | |
| 100 | // structured bindings cannot be captured by a lambda, therefore we introduce the helper variables a and b here | ||
| 101 | 2 | const auto& couplingPlaneLowerLeft = a; | |
| 102 | 2 | const auto& couplingPlaneUpperRight = b; | |
| 103 | |||
| 104 | // Create an array of lines. The entry for dim == couplingPlaneNormalDirectionIndex remains empty | ||
| 105 | 2 | std::array<std::optional<Line>, dim> lines; | |
| 106 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
|
6 | for (int i = 0; i < dim; ++i) |
| 107 | { | ||
| 108 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
|
4 | if (i == couplingPlaneNormalDirectionIndex) |
| 109 | continue; | ||
| 110 | |||
| 111 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
|
4 | lines[i] = [&]() |
| 112 | { | ||
| 113 | 4 | GlobalPosition tmp(couplingPlaneLowerLeft); | |
| 114 | 2 | tmp[i] = couplingPlaneUpperRight[i]; | |
| 115 | 2 | auto axis = std::bitset<dimWorld>(); | |
| 116 | 2 | axis.set(i, true); | |
| 117 | 2 | return Line(couplingPlaneLowerLeft, tmp, axis); | |
| 118 | 6 | }(); | |
| 119 | } | ||
| 120 | |||
| 121 | 2 | return lines; | |
| 122 | } | ||
| 123 | |||
| 124 | /*! | ||
| 125 | * \brief Returns the lowDim positions intersecting with a given line. | ||
| 126 | * | ||
| 127 | * \param bulkGridLowerLeft The lower left position of the bulk grid | ||
| 128 | * \param bulkGridUpperRight The upper right position of the bulk grid | ||
| 129 | * \param couplingPlaneNormal The normal vector of the coupling plane | ||
| 130 | * \param lowDimGridView The lowDim gridView | ||
| 131 | * \param lowDimGridData The lowDim grid data | ||
| 132 | * \param modelParamGroup The bulk model parameter group | ||
| 133 | */ | ||
| 134 | template<class LowDimGridView, class LowDimGridData> | ||
| 135 | static auto getPointsOnLine(const Dune::FieldVector<Scalar, 3>& bulkGridLowerLeft, | ||
| 136 | const Dune::FieldVector<Scalar, 3>& bulkGridUpperRight, | ||
| 137 | const Dune::FieldVector<Scalar, 3>& couplingPlaneNormal, | ||
| 138 | const LowDimGridView& lowDimGridView, | ||
| 139 | const LowDimGridData& lowDimGridData, | ||
| 140 | const std::string& modelParamGroup) | ||
| 141 | { | ||
| 142 | const auto couplingPlane = makeCouplingPlane(bulkGridLowerLeft, bulkGridUpperRight, couplingPlaneNormal, modelParamGroup); | ||
| 143 | const auto couplingPlaneNormalDirectionIndex = directionIndex(couplingPlaneNormal); | ||
| 144 | using ScalarVector = std::vector<Scalar>; | ||
| 145 | |||
| 146 | std::array<ScalarVector, 3> coordinates; | ||
| 147 | for (auto& c : coordinates) | ||
| 148 | c.reserve(lowDimGridView.size(1)); | ||
| 149 | |||
| 150 | ScalarVector poreRadius; | ||
| 151 | poreRadius.reserve(lowDimGridView.size(1)); | ||
| 152 | |||
| 153 | auto makeUnique = [](ScalarVector& v) | ||
| 154 | { | ||
| 155 | std::sort(v.begin(), v.end()); | ||
| 156 | v.erase(std::unique(v.begin(), v.end()), v.end()); | ||
| 157 | }; | ||
| 158 | |||
| 159 | for (const auto& vertex : vertices(lowDimGridView)) | ||
| 160 | { | ||
| 161 | const auto& pos = vertex.geometry().center(); | ||
| 162 | if (intersectsPointGeometry(pos, couplingPlane)) | ||
| 163 | { | ||
| 164 | poreRadius.push_back(lowDimGridData.parameters(vertex)[lowDimGridData.parameterIndex("PoreInscribedRadius")]); | ||
| 165 | for (int i = 0; i < 3; ++i) | ||
| 166 | coordinates[i].push_back(pos[i]); | ||
| 167 | } | ||
| 168 | } | ||
| 169 | |||
| 170 | if (std::any_of(poreRadius.begin(), poreRadius.end(), [&](auto& r) { return Dune::FloatCmp::eq(r, poreRadius[0]); })) | ||
| 171 | DUNE_THROW(Dune::GridError, "SnappyGridCreator in 3D currently only supports equal pore radii at the interface"); | ||
| 172 | |||
| 173 | for (auto& c : coordinates) | ||
| 174 | makeUnique(c); | ||
| 175 | |||
| 176 | struct Data { Scalar pos; Scalar radius; }; | ||
| 177 | std::array<std::optional<std::vector<Data>>, 3> result; | ||
| 178 | |||
| 179 | for (int i = 0; i < 3; ++i) | ||
| 180 | { | ||
| 181 | if (i != couplingPlaneNormalDirectionIndex) | ||
| 182 | { | ||
| 183 | std::vector<Data> tmp(coordinates[i].size()); | ||
| 184 | for (int poreIdx = 0; poreIdx < tmp.size(); ++poreIdx) | ||
| 185 | tmp[poreIdx] = Data{coordinates[i][poreIdx], poreRadius[0]}; | ||
| 186 | |||
| 187 | result[i] = std::move(tmp); | ||
| 188 | } | ||
| 189 | } | ||
| 190 | |||
| 191 | return result; | ||
| 192 | } | ||
| 193 | |||
| 194 | /*! | ||
| 195 | * \brief Returns the lowDim positions intersecting with a given line. | ||
| 196 | * | ||
| 197 | * \param bulkGridLowerLeft The lower left position of the bulk grid | ||
| 198 | * \param bulkGridUpperRight The upper right position of the bulk grid | ||
| 199 | * \param couplingPlaneNormal The normal vector of the coupling plane | ||
| 200 | * \param lowDimGridView The lowDim gridView | ||
| 201 | * \param lowDimGridData The lowDim grid data | ||
| 202 | * \param modelParamGroup The bulk model parameter group | ||
| 203 | */ | ||
| 204 | template<class LowDimGridView, class LowDimGridData> | ||
| 205 | 2 | static auto getPointsOnLine(const Dune::FieldVector<Scalar, 2>& bulkGridLowerLeft, | |
| 206 | const Dune::FieldVector<Scalar, 2>& bulkGridUpperRight, | ||
| 207 | const Dune::FieldVector<Scalar, 2>& couplingPlaneNormal, | ||
| 208 | const LowDimGridView& lowDimGridView, | ||
| 209 | const LowDimGridData& lowDimGridData, | ||
| 210 | const std::string& modelParamGroup) | ||
| 211 | { | ||
| 212 |
1/2✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
4 | std::vector<bool> vertexVisited(lowDimGridView.size(LowDimGridView::dimension), false); |
| 213 | |||
| 214 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | const auto axisParallelLines = makeAxisParallelLinesFromCouplingPlane(bulkGridLowerLeft, bulkGridUpperRight, couplingPlaneNormal, modelParamGroup); |
| 215 | |||
| 216 | //The line to check for intersections | ||
| 217 | 2 | const auto line = [&]() | |
| 218 | { | ||
| 219 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
6 | for (const auto& l : axisParallelLines) |
| 220 |
2/4✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
|
4 | if (l.has_value()) |
| 221 | 4 | return *l; | |
| 222 | ✗ | DUNE_THROW(Dune::InvalidStateException, "No line found"); | |
| 223 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | }(); |
| 224 | |||
| 225 | 6 | const auto orientation = line.corner(1) - line.corner(0); | |
| 226 | 2 | const auto dirIdx = directionIndex(orientation); | |
| 227 | |||
| 228 | struct Data { Scalar pos; Scalar radius;}; | ||
| 229 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
6 | std::vector<Data> data; |
| 230 | |||
| 231 |
4/4✓ Branch 1 taken 11 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 11 times.
✓ Branch 4 taken 2 times.
|
13 | for (const auto& element : elements(lowDimGridView)) |
| 232 | { | ||
| 233 |
2/2✓ Branch 0 taken 22 times.
✓ Branch 1 taken 11 times.
|
33 | for (std::size_t localVertexIdx = 0; localVertexIdx < 2; ++localVertexIdx) |
| 234 | { | ||
| 235 |
3/4✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 4 times.
✓ Branch 6 taken 18 times.
|
22 | if (const auto& pos = element.geometry().corner(localVertexIdx); intersectsPointGeometry(pos, line)) |
| 236 | { | ||
| 237 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | const auto& vertex = element.template subEntity<LowDimGridView::dimension>(localVertexIdx); |
| 238 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
|
4 | const auto vIdx = lowDimGridView.indexSet().index(vertex); |
| 239 |
7/16✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 4 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 4 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 4 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
|
8 | const auto poreRadius = lowDimGridData.parameters(vertex)[lowDimGridData.parameterIndex("PoreInscribedRadius")]; |
| 240 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | assert(poreRadius > 0.0); |
| 241 | |||
| 242 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
|
12 | if (vertexVisited[vIdx]) |
| 243 | ✗ | continue; | |
| 244 | else | ||
| 245 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
|
8 | vertexVisited[vIdx] = true; |
| 246 | |||
| 247 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
|
8 | data.push_back({pos[dirIdx], poreRadius}); |
| 248 | } | ||
| 249 | } | ||
| 250 | } | ||
| 251 | |||
| 252 | 6 | std::sort(data.begin(), data.end(), | |
| 253 | ✗ | [](const auto& pore0, const auto& pore1) | |
| 254 | { | ||
| 255 | ✗ | return pore0.pos < pore1.pos; | |
| 256 | }); | ||
| 257 | |||
| 258 | // check if any intersections were found | ||
| 259 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
|
4 | if (data.empty()) |
| 260 | { | ||
| 261 | ✗ | std::cout << "line boundaries: " << std::endl; | |
| 262 | ✗ | for (int k = 0; k < line.corners(); ++k) | |
| 263 | ✗ | std::cout << "(" << line.corner(k) << ") "; | |
| 264 | |||
| 265 | ✗ | std::cout << std::endl; | |
| 266 | ✗ | DUNE_THROW(Dune::GridError, "No coupled pores found"); | |
| 267 | } | ||
| 268 | |||
| 269 | 2 | std::array<std::optional<std::vector<Data>>, 2> result; | |
| 270 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
4 | result[dirIdx] = data; |
| 271 | 2 | return result; | |
| 272 | } | ||
| 273 | }; | ||
| 274 | |||
| 275 | } // end namespace Detail | ||
| 276 | |||
| 277 | /*! | ||
| 278 | * \ingroup FreeFlowPoreNetworkCoupling | ||
| 279 | * \brief A grid creator that matches a free-flow grid to a PNM grid. | ||
| 280 | */ | ||
| 281 | template<int dim, class OtherGridCreator, class DiscMethod = DiscretizationMethods::None> | ||
| 282 | class SnappyGridManager : public Dumux::GridManager<Dune::YaspGrid<dim, Dune::TensorProductCoordinates<typename OtherGridCreator::Grid::ctype, OtherGridCreator::Grid::LeafGridView::dimensionworld>>> | ||
| 283 | { | ||
| 284 | using Scalar = typename OtherGridCreator::Grid::ctype; | ||
| 285 | using OtherGrid = typename OtherGridCreator::Grid; | ||
| 286 | using IntVector = std::vector<int>; | ||
| 287 | using ScalarVector = std::vector<Scalar>; | ||
| 288 | |||
| 289 | static constexpr auto dimWorld = OtherGrid::LeafGridView::dimensionworld; | ||
| 290 | using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>; | ||
| 291 | |||
| 292 | using ParentType = Dumux::GridManager<Dune::YaspGrid<dim, Dune::TensorProductCoordinates<typename OtherGridCreator::Grid::ctype, OtherGridCreator::Grid::LeafGridView::dimensionworld>>>; | ||
| 293 | using LowDimGridData = typename OtherGridCreator::GridData; | ||
| 294 | |||
| 295 | struct GridConstructionData | ||
| 296 | { | ||
| 297 | std::array<ScalarVector, dim> positions; | ||
| 298 | std::array<IntVector, dim> cells; | ||
| 299 | std::array<ScalarVector, dim> grading; | ||
| 300 | std::array<std::optional<ScalarVector>, dim> interfacePositions; | ||
| 301 | }; | ||
| 302 | |||
| 303 | public: | ||
| 304 | using Grid = Dune::YaspGrid<dim, Dune::TensorProductCoordinates<Scalar, dim>>; | ||
| 305 | using ParentType::ParentType; | ||
| 306 | |||
| 307 | //! Make the grid | ||
| 308 | 2 | void init(const OtherGrid& otherGrid, const LowDimGridData& otherData, const std::string& modelParamGroup = "") | |
| 309 | { | ||
| 310 | 2 | modelParamGroup_ = modelParamGroup; | |
| 311 | 4 | std::array<ScalarVector, dim> positions; | |
| 312 | 2 | std::array<IntVector, dim> cells; | |
| 313 | 2 | std::array<ScalarVector, dim> grading; | |
| 314 | 2 | std::array<std::optional<ScalarVector>, dim> interFacePositions; | |
| 315 | using SnappyGridManagerHelper = Detail::SnappyGridManagerHelper<typename Grid::LeafGridView>; | ||
| 316 | |||
| 317 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | const auto lowerLeft = getParamFromGroup<GlobalPosition>(modelParamGroup, "Grid.LowerLeft"); |
| 318 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | const auto upperRight = getParamFromGroup<GlobalPosition>(modelParamGroup, "Grid.UpperRight"); |
| 319 | |||
| 320 |
2/4✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
|
4 | static const auto couplingPlaneNormal = getParamFromGroup<GlobalPosition>(modelParamGroup, |
| 321 | "Grid.CouplinglineNormal", | ||
| 322 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
8 | [](){ GlobalPosition tmp(0.0); tmp[tmp.size()-1] = 1.0; return tmp; }()); |
| 323 | |||
| 324 | 2 | const auto couplingPlaneNormalDirectionIndex = SnappyGridManagerHelper::directionIndex(couplingPlaneNormal); | |
| 325 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | const auto couplingPlane = SnappyGridManagerHelper::makeCouplingPlane(lowerLeft, upperRight, couplingPlaneNormal, modelParamGroup); |
| 326 | |||
| 327 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | std::cout << "plane: \n"; |
| 328 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
|
6 | for (int i = 0; i < couplingPlane.corners(); ++i) |
| 329 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
|
8 | std::cout << couplingPlane.corner(i) << std::endl; |
| 330 | |||
| 331 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
4 | const auto pointsOnAxisParallelLines = SnappyGridManagerHelper::getPointsOnLine(lowerLeft, upperRight, couplingPlaneNormal, otherGrid.leafGridView(), otherData, modelParamGroup_); |
| 332 | |||
| 333 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
|
6 | for (int i = 0; i < dim; ++i) |
| 334 | { | ||
| 335 | // set the lower left positions | ||
| 336 |
3/6✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
|
12 | positions[i].push_back(lowerLeft[i]); |
| 337 | |||
| 338 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
|
4 | if (i != couplingPlaneNormalDirectionIndex) // treat cells parallel to the coupling plane |
| 339 | { | ||
| 340 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
6 | if (!pointsOnAxisParallelLines[i].has_value()) |
| 341 | ✗ | DUNE_THROW(Dune::GridError, "Something went wrong with the coupling plane normal"); | |
| 342 | |||
| 343 | 4 | const auto& pointsOnLine = (*pointsOnAxisParallelLines[i]); | |
| 344 | |||
| 345 |
2/4✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
|
4 | std::cout << "point " << i << std::endl; |
| 346 |
4/4✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 2 times.
|
14 | for (auto x : pointsOnLine) |
| 347 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
|
4 | std::cout << x.pos << std::endl; |
| 348 | |||
| 349 | // check for user-defined additional points in the upstream area | ||
| 350 |
7/18✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
|
10 | const ScalarVector upstreamPositions = getParamFromGroup<ScalarVector>(modelParamGroup, "Grid.UpstreamPositions" + std::to_string(i), ScalarVector{}); |
| 351 | |||
| 352 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | addUpstreamPositions_(i, upstreamPositions, positions, pointsOnLine); |
| 353 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | addUpstreamCells_(i, upstreamPositions, cells); |
| 354 | |||
| 355 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | addCouplingPositions_(i, positions, pointsOnLine, interFacePositions, lowerLeft, upperRight); |
| 356 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | addCouplingCells_(i, cells, pointsOnLine); |
| 357 | |||
| 358 | // check for user-defined additional points in the downstream area | ||
| 359 |
9/24✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 2 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.
|
10 | const ScalarVector downstreamPositions = getParamFromGroup<ScalarVector>(modelParamGroup, "Grid.DownstreamPositions" + std::to_string(i), ScalarVector{}); |
| 360 | |||
| 361 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | addDownstreamPositions_(i, downstreamPositions, positions, upperRight); |
| 362 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | addDownstreamCells_(i, downstreamPositions, cells); |
| 363 | |||
| 364 | // set the upper right positions | ||
| 365 |
3/6✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
|
6 | positions[i].push_back(upperRight[i]); |
| 366 | |||
| 367 | // handle grading of the cells | ||
| 368 | // use 1 by default | ||
| 369 |
4/10✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
8 | grading[i].resize(positions[i].size()-1, 1.0); |
| 370 | |||
| 371 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | addUpstreamGrading_(i, upstreamPositions, grading); |
| 372 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | addDownstreamGrading_(i, downstreamPositions, positions, grading); |
| 373 | } | ||
| 374 | else // i == couplingPlaneNormalDirectionIndex | ||
| 375 | { | ||
| 376 | // treat the cells normal to the coupling plane | ||
| 377 | if (i != couplingPlaneNormalDirectionIndex) | ||
| 378 | DUNE_THROW(Dune::GridError, "Something went wrong with the coupling plane normal"); | ||
| 379 | |||
| 380 |
8/20✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
|
8 | const ScalarVector normalPositions = getParamFromGroup<ScalarVector>(modelParamGroup, "Grid.Positions" + std::to_string(i), ScalarVector{}); |
| 381 | |||
| 382 | // add positions, cells and grading | ||
| 383 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | addNormalPositions_(i, normalPositions, positions, upperRight); |
| 384 |
3/6✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
|
6 | positions[i].push_back(upperRight[i]); |
| 385 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | addNormalCells_(i, normalPositions, cells); |
| 386 |
4/10✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
8 | grading[i].resize(positions[i].size()-1, 1.0); |
| 387 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | addNormalGrading_(i, normalPositions, grading); |
| 388 | } | ||
| 389 | |||
| 390 |
5/10✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 4 times.
|
20 | if (positions[i].size() != cells[i].size() + 1) |
| 391 | ✗ | DUNE_THROW(Dune::GridError, "Wrong number of cells in " << i); | |
| 392 | } | ||
| 393 | |||
| 394 | // forward to the actual grid creator | ||
| 395 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | ParentType::init(positions, cells, grading, modelParamGroup); |
| 396 | |||
| 397 | 2 | gridConstructionData_ = GridConstructionData{std::move(positions), std::move(cells), std::move(grading), std::move(interFacePositions)}; | |
| 398 | 2 | } | |
| 399 | |||
| 400 | //! Return data used to create the snappy grid (needed for Dune::TensorProductCoordinates) and the locations of pores intersecting with the interface | ||
| 401 | const GridConstructionData& getGridConstructionData() const | ||
| 402 | { return gridConstructionData_; } | ||
| 403 | |||
| 404 | private: | ||
| 405 | |||
| 406 | ///////////////////////////////////////////////////// | ||
| 407 | //// Inlet ////////////////////////////////////////// | ||
| 408 | ///////////////////////////////////////////////////// | ||
| 409 | |||
| 410 | template<class PointsOnLine> | ||
| 411 | 2 | void addUpstreamPositions_(const int directionIndex, | |
| 412 | const ScalarVector& upstreamPositions, | ||
| 413 | std::array<ScalarVector, dim>& positions, | ||
| 414 | const PointsOnLine& points) | ||
| 415 | { | ||
| 416 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
|
4 | if (!upstreamPositions.empty()) |
| 417 | { | ||
| 418 | ✗ | const Scalar gridLowerLeft = positions[directionIndex][0]; | |
| 419 | ✗ | for (const Scalar pos : upstreamPositions) | |
| 420 | { | ||
| 421 | ✗ | if ((pos < points[0].pos - points[0].radius) && (pos > gridLowerLeft)) | |
| 422 | ✗ | positions[directionIndex].push_back(pos); | |
| 423 | else | ||
| 424 | ✗ | DUNE_THROW(Dune::RangeError, "Make sure to set positions only in the inlet"); | |
| 425 | } | ||
| 426 | } | ||
| 427 | 2 | } | |
| 428 | |||
| 429 | 2 | void addUpstreamCells_(const int directionIndex, | |
| 430 | const ScalarVector& upstreamPositions, | ||
| 431 | std::array<IntVector, dim>& cells) | ||
| 432 | { | ||
| 433 |
8/22✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✓ Branch 18 taken 2 times.
✗ 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.
|
12 | const IntVector cellsUpstream = getParamFromGroup<IntVector>(modelParamGroup_, |
| 434 | "Grid.UpstreamCells" + std::to_string(directionIndex), | ||
| 435 | IntVector{}); | ||
| 436 | |||
| 437 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
|
4 | if (cellsUpstream.empty()) |
| 438 | ✗ | return; | |
| 439 | |||
| 440 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
6 | if (cellsUpstream.size() != upstreamPositions.size() + 1) |
| 441 | ✗ | DUNE_THROW(Dumux::ParameterException, "UpstreamCells" << std::to_string(directionIndex) << " must equal UpstreamPositions" << std::to_string(directionIndex) << " + 1"); | |
| 442 | |||
| 443 |
4/4✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 2 times.
|
6 | for (int c : cellsUpstream) |
| 444 | { | ||
| 445 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if (c < 1) |
| 446 | ✗ | DUNE_THROW(Dumux::ParameterException, "UpstreamCells" << std::to_string(directionIndex) << " must have values greater than zero."); | |
| 447 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
4 | cells[directionIndex].push_back(c); |
| 448 | } | ||
| 449 | } | ||
| 450 | |||
| 451 | 2 | void addUpstreamGrading_(const int directionIndex, | |
| 452 | const ScalarVector& upstreamPositions, | ||
| 453 | std::array<ScalarVector, dim>& grading) | ||
| 454 | { | ||
| 455 |
2/4✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
|
4 | if (upstreamPositions.empty()) |
| 456 | 2 | return; | |
| 457 | |||
| 458 | ✗ | const ScalarVector upstreamGrading = getParamFromGroup<ScalarVector>(modelParamGroup_, "Grid.UpstreamGrading" + std::to_string(directionIndex), ScalarVector{}); | |
| 459 | |||
| 460 | ✗ | if (!upstreamGrading.empty()) | |
| 461 | { | ||
| 462 | ✗ | if (upstreamGrading.size() != upstreamPositions.size() + 1) | |
| 463 | ✗ | DUNE_THROW(Dune::RangeError, "UpstreamGrading" << std::to_string(directionIndex) << " must equal UpstreamPositions" << std::to_string(directionIndex) << " + 1"); | |
| 464 | |||
| 465 | ✗ | for (int i = 0; i < upstreamPositions.size() + 1; ++i) | |
| 466 | ✗ | grading[directionIndex][i] = upstreamGrading[i]; | |
| 467 | } | ||
| 468 | } | ||
| 469 | |||
| 470 | ///////////////////////////////////////////////////// | ||
| 471 | //// Outlet ///////////////////////////////////////// | ||
| 472 | ///////////////////////////////////////////////////// | ||
| 473 | |||
| 474 | 2 | void addDownstreamPositions_(const int directionIndex, | |
| 475 | const ScalarVector& downstreamPositions, | ||
| 476 | std::array<ScalarVector, dim>& gridPositions, | ||
| 477 | const GlobalPosition& gridUpperRight) | ||
| 478 | { | ||
| 479 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
|
4 | if (!downstreamPositions.empty()) |
| 480 | { | ||
| 481 | ✗ | for (const Scalar pos : downstreamPositions) | |
| 482 | { | ||
| 483 | ✗ | if ((pos > gridPositions[directionIndex].back()) && (pos < gridUpperRight[directionIndex])) | |
| 484 | ✗ | gridPositions[directionIndex].push_back(pos); | |
| 485 | else | ||
| 486 | ✗ | DUNE_THROW(Dune::RangeError, "Make sure to set ascending positions only in the outlet"); | |
| 487 | } | ||
| 488 | } | ||
| 489 | 2 | } | |
| 490 | |||
| 491 | 2 | void addDownstreamCells_(const int directionIndex, | |
| 492 | const ScalarVector& downstreamPositions, | ||
| 493 | std::array<IntVector, dim>& cells) | ||
| 494 | { | ||
| 495 |
8/22✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✓ Branch 18 taken 2 times.
✗ 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.
|
12 | const IntVector downstreamcells = getParamFromGroup<IntVector>(modelParamGroup_, |
| 496 | "Grid.DownstreamCells" + std::to_string(directionIndex), | ||
| 497 | IntVector{}); | ||
| 498 | |||
| 499 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
|
4 | if (downstreamcells.empty()) |
| 500 | ✗ | return; | |
| 501 | |||
| 502 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
6 | if (downstreamcells.size() != downstreamPositions.size() + 1) |
| 503 | ✗ | DUNE_THROW(Dumux::ParameterException, "DownstreamCells" << std::to_string(directionIndex) << " must equal DownstreamPositions" << std::to_string(directionIndex) << " + 1"); | |
| 504 | |||
| 505 |
4/4✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 2 times.
|
6 | for (int c : downstreamcells) |
| 506 | { | ||
| 507 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if (c < 1) |
| 508 | ✗ | DUNE_THROW(Dumux::ParameterException, "DownstreamCells" << std::to_string(directionIndex) << " must have values greater than zero."); | |
| 509 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
4 | cells[directionIndex].push_back(c); |
| 510 | } | ||
| 511 | } | ||
| 512 | |||
| 513 | 2 | void addDownstreamGrading_(const int directionIndex, | |
| 514 | const ScalarVector& downstreamPositions, | ||
| 515 | std::array<ScalarVector, dim>& gridPositions, | ||
| 516 | std::array<ScalarVector, dim>& grading) | ||
| 517 | { | ||
| 518 |
2/4✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
|
4 | if (downstreamPositions.empty()) |
| 519 | 2 | return; | |
| 520 | |||
| 521 | ✗ | const ScalarVector downstreamGrading = getParamFromGroup<ScalarVector>(modelParamGroup_, "Grid.DownstreamGrading" + std::to_string(directionIndex), ScalarVector{}); | |
| 522 | |||
| 523 | ✗ | if (!downstreamGrading.empty()) | |
| 524 | { | ||
| 525 | ✗ | if (downstreamGrading.size() != downstreamPositions.size() + 1) | |
| 526 | ✗ | DUNE_THROW(Dune::RangeError, "DownstreamGrading" << std::to_string(directionIndex) << " must equal DownstreamPositions" << std::to_string(directionIndex) << " + 1"); | |
| 527 | |||
| 528 | ✗ | const int offSet = gridPositions[directionIndex].size() - downstreamPositions.size() - 2; | |
| 529 | ✗ | for (int i = 0; i < downstreamPositions.size() + 1; ++i) | |
| 530 | ✗ | grading[directionIndex][offSet + i] = downstreamGrading[i]; | |
| 531 | } | ||
| 532 | } | ||
| 533 | |||
| 534 | ///////////////////////////////////////////////////// | ||
| 535 | //// Coupling interface////////////////////////////// | ||
| 536 | ///////////////////////////////////////////////////// | ||
| 537 | |||
| 538 | template<class PointsOnLine> | ||
| 539 | 2 | void addCouplingPositions_(const int directionIndex, | |
| 540 | std::array<ScalarVector, dim>& positions, | ||
| 541 | const PointsOnLine& points, | ||
| 542 | std::array<std::optional<ScalarVector>, dim>& interFacePositions, | ||
| 543 | const GlobalPosition& gridLowerLeft, | ||
| 544 | const GlobalPosition& gridUpperRight) | ||
| 545 | { | ||
| 546 | // create an empty vector for the given directionIndex | ||
| 547 |
1/2✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
|
6 | interFacePositions[directionIndex] = ScalarVector{}; |
| 548 | |||
| 549 | // set the points for the pore body positions | ||
| 550 |
4/4✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 2 times.
|
10 | for (const auto& point : points) |
| 551 | { | ||
| 552 | 4 | const auto left = point.pos - point.radius; | |
| 553 | 4 | const auto right = point.pos + point.radius; | |
| 554 | |||
| 555 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
|
12 | if (left < positions[directionIndex].back()) |
| 556 | ✗ | DUNE_THROW(Dune::RangeError, "Pore body radii are too large, they intersect!"); | |
| 557 | |||
| 558 |
2/4✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
|
8 | if (left > gridLowerLeft[directionIndex]) |
| 559 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
|
8 | positions[directionIndex].push_back(left); |
| 560 | |||
| 561 |
2/4✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
|
8 | if (right < gridUpperRight[directionIndex]) |
| 562 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
|
8 | positions[directionIndex].push_back(right); |
| 563 | |||
| 564 |
3/6✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
|
12 | interFacePositions[directionIndex]->push_back(left); |
| 565 |
3/6✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
|
12 | interFacePositions[directionIndex]->push_back(right); |
| 566 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
|
12 | assert(interFacePositions[directionIndex].has_value()); |
| 567 | } | ||
| 568 | 2 | } | |
| 569 | |||
| 570 | template<class PointsOnLine> | ||
| 571 | 2 | void addCouplingCells_(const int directionIndex, | |
| 572 | std::array<IntVector, dim>& cells, | ||
| 573 | const PointsOnLine& points) | ||
| 574 | { | ||
| 575 | // set the number of cells above the porous medium | ||
| 576 | 2 | const int cellsPerPore = getParamFromGroup<int>(modelParamGroup_, "Grid.CellsPerPore"); | |
| 577 | 2 | const int fixedCellsBetweenPores = getParamFromGroup<int>(modelParamGroup_, "Grid.FixedCellsBetweenPores", -1); | |
| 578 | |||
| 579 | // set the number of cells between the pore bodies | ||
| 580 |
4/4✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 2 times.
|
6 | for (int i = 0; i < points.size(); ++i) |
| 581 | { | ||
| 582 | 8 | cells[directionIndex].push_back(cellsPerPore); | |
| 583 |
4/4✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 2 times.
|
8 | if (i < points.size() -1) |
| 584 | { | ||
| 585 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if (fixedCellsBetweenPores > 0) |
| 586 | ✗ | cells[directionIndex].push_back(fixedCellsBetweenPores); | |
| 587 | else | ||
| 588 | { | ||
| 589 | // set number of cells such that the spacing of the cells between the left and right pore body | ||
| 590 | // is the average of the spacings for the cells coupled to those two pore bodies | ||
| 591 | 2 | const auto spacingLeft = points[i].radius*2.0 / cellsPerPore; | |
| 592 | 2 | const auto spacingRight = points[i+1].radius*2.0 / cellsPerPore; | |
| 593 | 2 | const auto avgSpacing = (spacingLeft + spacingRight) / 2; | |
| 594 | 4 | const auto lengthBetween = (points[i+1].pos - (points[i+1].radius)) | |
| 595 | 2 | - (points[i].pos + (points[i].radius)); | |
| 596 | 2 | const int cellsBetween = std::ceil(lengthBetween / avgSpacing); | |
| 597 | 4 | cells[directionIndex].push_back(cellsBetween); | |
| 598 | } | ||
| 599 | } | ||
| 600 | } | ||
| 601 | 2 | } | |
| 602 | |||
| 603 | ///////////////////////////////////////////////////// | ||
| 604 | //// Normal direction /////////////////////////////// | ||
| 605 | ///////////////////////////////////////////////////// | ||
| 606 | |||
| 607 | 2 | void addNormalPositions_(const int directionIndex, | |
| 608 | const ScalarVector& normalPositions, | ||
| 609 | std::array<ScalarVector, dim>& gridPositions, | ||
| 610 | const GlobalPosition& gridUpperRight) | ||
| 611 | { | ||
| 612 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
|
4 | if (!normalPositions.empty()) |
| 613 | { | ||
| 614 | ✗ | for (const Scalar pos : normalPositions) | |
| 615 | { | ||
| 616 | ✗ | if ((pos > gridPositions[directionIndex].back()) && (pos < gridUpperRight[directionIndex])) | |
| 617 | ✗ | gridPositions[directionIndex].push_back(pos); | |
| 618 | else | ||
| 619 | ✗ | DUNE_THROW(Dune::RangeError, "Make sure to set ascending normal positions"); | |
| 620 | } | ||
| 621 | } | ||
| 622 | 2 | } | |
| 623 | |||
| 624 | 2 | void addNormalCells_(const int directionIndex, | |
| 625 | const ScalarVector& normalPositions, | ||
| 626 | std::array<IntVector, dim>& cells) | ||
| 627 | { | ||
| 628 |
5/14✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 2 times.
✓ Branch 11 taken 2 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
|
6 | const IntVector cellsNormal = getParamFromGroup<IntVector>(modelParamGroup_, "Grid.Cells" + std::to_string(directionIndex)); |
| 629 | |||
| 630 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
6 | if (cellsNormal.size() != normalPositions.size() + 1) |
| 631 | ✗ | DUNE_THROW(Dumux::ParameterException, "Cells" << std::to_string(directionIndex) << " must be of size " << normalPositions.size() + 1); | |
| 632 | |||
| 633 |
4/4✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 2 times.
|
6 | for (int c : cellsNormal) |
| 634 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
4 | cells[directionIndex].push_back(c); |
| 635 | 2 | } | |
| 636 | |||
| 637 | 2 | void addNormalGrading_(const int directionIndex, | |
| 638 | const ScalarVector& normalPositions, | ||
| 639 | std::array<ScalarVector, dim>& grading) | ||
| 640 | { | ||
| 641 |
8/22✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
|
8 | const ScalarVector normalGrading = getParamFromGroup<ScalarVector>(modelParamGroup_, "Grid.Grading" + std::to_string(directionIndex), ScalarVector{}); |
| 642 | |||
| 643 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
|
4 | if (!normalGrading.empty()) |
| 644 | { | ||
| 645 | ✗ | if (normalGrading.size() != normalPositions.size() + 1) | |
| 646 | ✗ | DUNE_THROW(Dune::RangeError, "Grading" << std::to_string(directionIndex) << " must be of size " << normalPositions.size() + 1); | |
| 647 | |||
| 648 | ✗ | for (int i = 0; i < normalPositions.size() + 1; ++i) | |
| 649 | ✗ | grading[directionIndex][i] = normalGrading[i]; | |
| 650 | } | ||
| 651 | 2 | } | |
| 652 | |||
| 653 | |||
| 654 | std::string modelParamGroup_ = ""; | ||
| 655 | GridConstructionData gridConstructionData_; | ||
| 656 | }; | ||
| 657 | |||
| 658 | } // end namespace Dumux::PoreNetwork | ||
| 659 | |||
| 660 | #endif | ||
| 661 |