| 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 InputOutput | ||
| 10 | * \brief Grid manager specialization for YaspGrid | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_IO_GRID_MANAGER_YASP_HH | ||
| 13 | #define DUMUX_IO_GRID_MANAGER_YASP_HH | ||
| 14 | |||
| 15 | #include <dune/common/math.hh> | ||
| 16 | |||
| 17 | #include <dune/grid/yaspgrid.hh> | ||
| 18 | #include <dune/grid/io/file/dgfparser/dgfyasp.hh> | ||
| 19 | |||
| 20 | #ifndef DUMUX_IO_GRID_MANAGER_BASE_HH | ||
| 21 | #include <dumux/io/grid/gridmanager_base.hh> | ||
| 22 | #endif | ||
| 23 | |||
| 24 | #include <dumux/common/gridcapabilities.hh> | ||
| 25 | |||
| 26 | namespace Dumux { | ||
| 27 | |||
| 28 | /*! | ||
| 29 | * \ingroup InputOutput | ||
| 30 | * \brief Provides a grid manager for YaspGrids | ||
| 31 | * from information in the input file | ||
| 32 | * | ||
| 33 | * All keys are expected to be in group GridParameterGroup. | ||
| 34 | * The following keys are recognized: | ||
| 35 | * - File : a DGF file to load the coarse grid from | ||
| 36 | * - LowerLeft : lower left corner of the domain (only for EquidistantOffsetCoordinates, otherwise 0-origin) | ||
| 37 | * - UpperRight : upper right corner of the domain | ||
| 38 | * - Cells : the number of cells in each direction | ||
| 39 | * - Periodic : true or false for each direction | ||
| 40 | * - Overlap : overlap size in cells | ||
| 41 | * - Partitioning : a non-standard load-balancing, number of processors per direction | ||
| 42 | * - KeepPyhsicalOverlap : whether to keep the physical overlap | ||
| 43 | * in physical size or in number of cells upon refinement | ||
| 44 | * - Refinement : the number of global refines to apply initially. | ||
| 45 | * | ||
| 46 | */ | ||
| 47 | template<class Coordinates, int dim> | ||
| 48 |
4/6✓ Branch 1 taken 20 times.
✓ Branch 2 taken 32 times.
✓ Branch 3 taken 12 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
|
619 | class GridManager<Dune::YaspGrid<dim, Coordinates>> |
| 49 | : public GridManagerBase<Dune::YaspGrid<dim, Coordinates>> | ||
| 50 | { | ||
| 51 | using ct = typename Dune::YaspGrid<dim, Coordinates>::ctype; | ||
| 52 | using GlobalPosition = Dune::FieldVector<ct, dim>; | ||
| 53 | public: | ||
| 54 | using Grid = typename Dune::YaspGrid<dim, Coordinates>; | ||
| 55 | using ParentType = GridManagerBase<Grid>; | ||
| 56 | |||
| 57 | /*! | ||
| 58 | * \brief Make the grid. This is implemented by specializations of this method. | ||
| 59 | */ | ||
| 60 | 395 | void init(const std::string& modelParamGroup = "") | |
| 61 | { | ||
| 62 | // First try to create it from a DGF file in GridParameterGroup.File | ||
| 63 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 331 times.
|
790 | if (hasParamInGroup(modelParamGroup, "Grid.File")) |
| 64 | { | ||
| 65 | ✗ | const auto filename = getParamFromGroup<std::string>(modelParamGroup, "Grid.File"); | |
| 66 | ✗ | const std::string extension = ParentType::getFileExtension(filename); | |
| 67 | |||
| 68 | ✗ | if (extension != "dgf" && extension != "vti") | |
| 69 | ✗ | DUNE_THROW(Dune::IOError, "Grid type " << Dune::className<Grid>() << " doesn't support grid files with extension: *."<< extension); | |
| 70 | |||
| 71 | // VTK file formats for image grids | ||
| 72 | ✗ | if (extension == "vti") | |
| 73 | { | ||
| 74 | #if DUMUX_HAVE_GRIDFORMAT | ||
| 75 | VTKReader vtkReader(filename); | ||
| 76 | VTKReader::Data cellData, pointData; | ||
| 77 | |||
| 78 | const bool verbose = getParamFromGroup<bool>(modelParamGroup, "Grid.Verbosity", false); | ||
| 79 | std::shared_ptr<Detail::VTKReader::ImageGrid<ct, dim>> gridInfo | ||
| 80 | = vtkReader.readStructuredGrid<ct, dim>(cellData, pointData, verbose); | ||
| 81 | const std::bitset<dim> periodic; | ||
| 82 | const int overlap = getParamFromGroup<int>(modelParamGroup, "Grid.Overlap", 1); | ||
| 83 | const auto upperRight = gridInfo->upperRight(); | ||
| 84 | const auto cells = gridInfo->cells(); | ||
| 85 | const auto lowerLeft = gridInfo->lowerLeft(); | ||
| 86 | |||
| 87 | if constexpr (std::is_same_v<Dune::EquidistantCoordinates<ct, dim>, Coordinates>) | ||
| 88 | init_(upperRight, cells, modelParamGroup, overlap, periodic); | ||
| 89 | else | ||
| 90 | init_(lowerLeft, upperRight, cells, modelParamGroup, overlap, periodic); | ||
| 91 | |||
| 92 | ParentType::gridData_ = std::make_shared<GridData<Grid>>( | ||
| 93 | ParentType::gridPtr(), gridInfo, std::move(cellData), std::move(pointData) | ||
| 94 | ); | ||
| 95 | |||
| 96 | ParentType::enableVtkData_ = true; | ||
| 97 | #else | ||
| 98 | ✗ | DUNE_THROW(Dune::IOError, "Grid type " << Dune::className<Grid>() << " doesn't support grid files with extension: *."<< extension); | |
| 99 | #endif | ||
| 100 | } | ||
| 101 | else | ||
| 102 | { | ||
| 103 | ✗ | ParentType::makeGridFromDgfFile( | |
| 104 | getParamFromGroup<std::string>(modelParamGroup, "Grid.File") | ||
| 105 | ); | ||
| 106 | } | ||
| 107 | |||
| 108 | ✗ | postProcessing_(modelParamGroup); | |
| 109 | ✗ | return; | |
| 110 | ✗ | } | |
| 111 | |||
| 112 | // Then look for the necessary keys to construct from the input file | ||
| 113 |
1/2✓ Branch 2 taken 331 times.
✗ Branch 3 not taken.
|
790 | else if (hasParamInGroup(modelParamGroup, "Grid.UpperRight")) |
| 114 | { | ||
| 115 | // get the upper right corner coordinates | ||
| 116 | 395 | const auto upperRight = getParamFromGroup<GlobalPosition>(modelParamGroup, "Grid.UpperRight"); | |
| 117 | |||
| 118 | // number of cells in each direction | ||
| 119 | 395 | std::array<int, dim> cells; cells.fill(1); | |
| 120 | 395 | cells = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Cells", cells); | |
| 121 | |||
| 122 | // \todo TODO periodic boundaries with yasp (the periodicity concept of yasp grid is currently not supported, use dune-spgrid) | ||
| 123 | // const auto periodic = getParamFromGroup<std::bitset<dim>>(modelParamGroup, "Grid.Periodic", std::bitset<dim>{}); | ||
| 124 | 395 | const std::bitset<dim> periodic; | |
| 125 | |||
| 126 | // get the overlap | ||
| 127 | 395 | const int overlap = getParamFromGroup<int>(modelParamGroup, "Grid.Overlap", 1); | |
| 128 | |||
| 129 | if constexpr (std::is_same_v<Dune::EquidistantCoordinates<ct, dim>, Coordinates>) | ||
| 130 | 270 | init(upperRight, cells, modelParamGroup, overlap, periodic); | |
| 131 | else | ||
| 132 | { | ||
| 133 | 236 | const auto lowerLeft = getParamFromGroup<GlobalPosition>(modelParamGroup, "Grid.LowerLeft", GlobalPosition(0.0)); | |
| 134 | 125 | init(lowerLeft, upperRight, cells, modelParamGroup, overlap, periodic); | |
| 135 | } | ||
| 136 | } | ||
| 137 | |||
| 138 | // Didn't find a way to construct the grid | ||
| 139 | else | ||
| 140 | { | ||
| 141 | ✗ | const auto prefix = modelParamGroup.empty() ? modelParamGroup : modelParamGroup + "."; | |
| 142 |
0/30✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
|
395 | DUNE_THROW(ParameterException, "Please supply one of the parameters " |
| 143 | << prefix + "Grid.UpperRight" | ||
| 144 | << ", or a grid file in " << prefix + "Grid.File"); | ||
| 145 | |||
| 146 | ✗ | } | |
| 147 | } | ||
| 148 | |||
| 149 | /*! | ||
| 150 | * \brief Make the grid using input data not read from the input file. | ||
| 151 | * \note Use this function for EquidistantCoordinates where lowerLeft is always zero. | ||
| 152 | */ | ||
| 153 | 236 | void init(const GlobalPosition& upperRight, | |
| 154 | const std::array<int, dim>& cells, | ||
| 155 | const std::string& modelParamGroup = "", | ||
| 156 | const int overlap = 1, | ||
| 157 | const std::bitset<dim> periodic = std::bitset<dim>{}) | ||
| 158 | { | ||
| 159 | 236 | init_(upperRight, cells, modelParamGroup, overlap, periodic); | |
| 160 | 236 | postProcessing_(modelParamGroup); | |
| 161 | } | ||
| 162 | |||
| 163 | /*! | ||
| 164 | * \brief Make the grid using input data not read from the input file. | ||
| 165 | * \note Use this function for EquidistantOffsetCoordinates. | ||
| 166 | */ | ||
| 167 | 99 | void init(const GlobalPosition& lowerLeft, | |
| 168 | const GlobalPosition& upperRight, | ||
| 169 | const std::array<int, dim>& cells, | ||
| 170 | const std::string& modelParamGroup = "", | ||
| 171 | const int overlap = 1, | ||
| 172 | const std::bitset<dim> periodic = std::bitset<dim>{}) | ||
| 173 | { | ||
| 174 | 99 | init_(lowerLeft, upperRight, cells, modelParamGroup, overlap, periodic); | |
| 175 | 99 | postProcessing_(modelParamGroup); | |
| 176 | } | ||
| 177 | |||
| 178 | private: | ||
| 179 | /*! | ||
| 180 | * \brief Make the grid using input data not read from the input file. | ||
| 181 | * \note Use this function for EquidistantCoordinates where lowerLeft is always zero. | ||
| 182 | */ | ||
| 183 | 240 | void init_(const GlobalPosition& upperRight, | |
| 184 | const std::array<int, dim>& cells, | ||
| 185 | const std::string& modelParamGroup = "", | ||
| 186 | const int overlap = 1, | ||
| 187 | const std::bitset<dim> periodic = std::bitset<dim>{}) | ||
| 188 | { | ||
| 189 | static_assert(std::is_same_v<Dune::EquidistantCoordinates<ct, dim>, Coordinates>, | ||
| 190 | "Use init function taking lowerLeft as argument when working with EquidistantOffsetCoordinates"); | ||
| 191 | |||
| 192 |
1/2✓ Branch 2 taken 236 times.
✗ Branch 3 not taken.
|
480 | if (!hasParamInGroup(modelParamGroup, "Grid.Partitioning")) |
| 193 | { | ||
| 194 | // construct using default load balancing | ||
| 195 |
2/4✓ Branch 2 taken 236 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 236 times.
✗ Branch 6 not taken.
|
480 | ParentType::gridPtr() = std::make_unique<Grid>(upperRight, cells, periodic, overlap); |
| 196 | } | ||
| 197 | else | ||
| 198 | { | ||
| 199 | // construct using user defined partitioning | ||
| 200 | ✗ | const auto partitioning = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Partitioning"); | |
| 201 | ✗ | Dune::Yasp::FixedSizePartitioning<dim> lb(partitioning); | |
| 202 | ✗ | ParentType::gridPtr() = std::make_unique<Grid>(upperRight, cells, periodic, overlap, typename Grid::Communication(), &lb); | |
| 203 | ✗ | } | |
| 204 | 240 | } | |
| 205 | |||
| 206 | /*! | ||
| 207 | * \brief Make the grid using input data not read from the input file. | ||
| 208 | * \note Use this function for EquidistantOffsetCoordinates. | ||
| 209 | */ | ||
| 210 | 99 | void init_(const GlobalPosition& lowerLeft, | |
| 211 | const GlobalPosition& upperRight, | ||
| 212 | const std::array<int, dim>& cells, | ||
| 213 | const std::string& modelParamGroup = "", | ||
| 214 | const int overlap = 1, | ||
| 215 | const std::bitset<dim> periodic = std::bitset<dim>{}) | ||
| 216 | { | ||
| 217 | static_assert(std::is_same_v<Dune::EquidistantOffsetCoordinates<ct, dim>, Coordinates>, | ||
| 218 | "LowerLeft can only be specified with EquidistantOffsetCoordinates"); | ||
| 219 | |||
| 220 |
1/2✓ Branch 2 taken 99 times.
✗ Branch 3 not taken.
|
198 | if (!hasParamInGroup(modelParamGroup, "Grid.Partitioning")) |
| 221 | { | ||
| 222 | // construct using default load balancing | ||
| 223 |
2/4✓ Branch 2 taken 99 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 99 times.
✗ Branch 6 not taken.
|
198 | ParentType::gridPtr() = std::make_unique<Grid>(lowerLeft, upperRight, cells, periodic, overlap); |
| 224 | } | ||
| 225 | else | ||
| 226 | { | ||
| 227 | // construct using user defined partitioning | ||
| 228 | ✗ | const auto partitioning = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Partitioning"); | |
| 229 | ✗ | Dune::Yasp::FixedSizePartitioning<dim> lb(partitioning); | |
| 230 | ✗ | ParentType::gridPtr() = std::make_unique<Grid>(lowerLeft, upperRight, cells, periodic, overlap, typename Grid::Communication(), &lb); | |
| 231 | ✗ | } | |
| 232 | 99 | } | |
| 233 | |||
| 234 | /*! | ||
| 235 | * \brief Postprocessing for YaspGrid | ||
| 236 | */ | ||
| 237 | 399 | void postProcessing_(const std::string& modelParamGroup) | |
| 238 | { | ||
| 239 | // Check if should refine the grid | ||
| 240 | 399 | const bool keepPhysicalOverlap = getParamFromGroup<bool>(modelParamGroup, "Grid.KeepPhysicalOverlap", true); | |
| 241 | 399 | ParentType::grid().refineOptions(keepPhysicalOverlap); | |
| 242 | 399 | ParentType::maybeRefineGrid(modelParamGroup); | |
| 243 | 399 | ParentType::loadBalance(); | |
| 244 | 399 | } | |
| 245 | }; | ||
| 246 | |||
| 247 | /*! | ||
| 248 | * \ingroup InputOutput | ||
| 249 | * \brief Provides a grid manager for YaspGrids with different zones and grading | ||
| 250 | * | ||
| 251 | * All keys are expected to be in group GridParameterGroup. | ||
| 252 | * The following keys are recognized: | ||
| 253 | * - Positions0 : position array for x-coordinate | ||
| 254 | * - Positions1 : position array for y-coordinate | ||
| 255 | * - Positions2 : position array for z-coordinate | ||
| 256 | * - Cells0 : number of cells array for x-coordinate | ||
| 257 | * - Cells1 : number of cells array for y-coordinate | ||
| 258 | * - Cells2 : number of cells array for z-coordinate | ||
| 259 | * - Grading0 : grading factor array for x-coordinate | ||
| 260 | * - Grading1 : grading factor array for y-coordinate | ||
| 261 | * - Grading2 : grading factor array for z-coordinate | ||
| 262 | * - Verbosity : whether the grid construction should output to standard out | ||
| 263 | * - Periodic : true or false for each direction | ||
| 264 | * - Overlap : overlap size in cells | ||
| 265 | * - Partitioning : a non-standard load-balancing, number of processors per direction | ||
| 266 | * - KeepPyhsicalOverlap : whether to keep the physical overlap | ||
| 267 | * in physical size or in number of cells upon refinement | ||
| 268 | * - Refinement : the number of global refines to apply initially. | ||
| 269 | * | ||
| 270 | * The grading factor \f$ g \f$ specifies the ratio between the next and the current cell size: | ||
| 271 | * \f$ g = \frac{h_{i+1}}{h_i} \f$. | ||
| 272 | * Negative grading factors are converted to | ||
| 273 | * \f$ g = -\frac{1}{g_\textrm{negative}} \f$ | ||
| 274 | * to avoid issues with imprecise fraction numbers. | ||
| 275 | */ | ||
| 276 | template<class ctype, int dim> | ||
| 277 |
4/6✓ Branch 1 taken 11 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
|
168 | class GridManager<Dune::YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> >> |
| 278 | : public GridManagerBase<Dune::YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> > > | ||
| 279 | { | ||
| 280 | public: | ||
| 281 | using Grid = typename Dune::YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> >; | ||
| 282 | using ParentType = GridManagerBase<Grid>; | ||
| 283 | |||
| 284 | /*! | ||
| 285 | * \brief Make the grid. This is implemented by specializations of this method. | ||
| 286 | */ | ||
| 287 | 67 | void init(const std::string& modelParamGroup = "") | |
| 288 | { | ||
| 289 | // Only construction from the input file is possible | ||
| 290 | // Look for the necessary keys to construct from the input file | ||
| 291 | // The positions | ||
| 292 | 67 | std::array<std::vector<ctype>, dim> positions; | |
| 293 |
2/2✓ Branch 0 taken 123 times.
✓ Branch 1 taken 67 times.
|
190 | for (int i = 0; i < dim; ++i) |
| 294 |
3/6✓ Branch 2 taken 123 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 123 times.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 123 times.
|
369 | positions[i] = getParamFromGroup<std::vector<ctype>>(modelParamGroup, "Grid.Positions" + std::to_string(i)); |
| 295 | |||
| 296 | // the number of cells (has a default) | ||
| 297 | 67 | std::array<std::vector<int>, dim> cells; | |
| 298 |
2/2✓ Branch 0 taken 123 times.
✓ Branch 1 taken 67 times.
|
190 | for (int i = 0; i < dim; ++i) |
| 299 | { | ||
| 300 |
1/2✓ Branch 1 taken 123 times.
✗ Branch 2 not taken.
|
123 | cells[i].resize(positions[i].size()-1, 1.0); |
| 301 |
2/4✓ Branch 2 taken 123 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 123 times.
✗ Branch 6 not taken.
|
369 | cells[i] = getParamFromGroup<std::vector<int>>(modelParamGroup, "Grid.Cells" + std::to_string(i), cells[i]); |
| 302 | } | ||
| 303 | |||
| 304 | // grading factor (has a default) | ||
| 305 | 67 | std::array<std::vector<ctype>, dim> grading; | |
| 306 |
2/2✓ Branch 0 taken 123 times.
✓ Branch 1 taken 67 times.
|
190 | for (int i = 0; i < dim; ++i) |
| 307 | { | ||
| 308 |
1/2✓ Branch 1 taken 123 times.
✗ Branch 2 not taken.
|
123 | grading[i].resize(positions[i].size()-1, 1.0); |
| 309 |
3/6✓ Branch 2 taken 123 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 123 times.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 123 times.
|
369 | grading[i] = getParamFromGroup<std::vector<ctype>>(modelParamGroup, "Grid.Grading" + std::to_string(i), grading[i]); |
| 310 | } | ||
| 311 | |||
| 312 | // call the generic function | ||
| 313 |
1/2✓ Branch 1 taken 67 times.
✗ Branch 2 not taken.
|
67 | init(positions, cells, grading, modelParamGroup); |
| 314 | 134 | } | |
| 315 | |||
| 316 | /*! | ||
| 317 | * \brief Make the grid using input data not read from the input file. | ||
| 318 | */ | ||
| 319 | 101 | void init(const std::array<std::vector<ctype>, dim>& positions, | |
| 320 | const std::array<std::vector<int>, dim>& cells, | ||
| 321 | const std::array<std::vector<ctype>, dim>& grading, | ||
| 322 | const std::string& modelParamGroup = "") | ||
| 323 | { | ||
| 324 | // Additional parameters (they have a default) | ||
| 325 | 101 | const int overlap = getParamFromGroup<int>(modelParamGroup, "Grid.Overlap", 1); | |
| 326 | 101 | const bool verbose = getParamFromGroup<bool>(modelParamGroup, "Grid.Verbosity", false); | |
| 327 | // \todo TODO periodic boundaries with yasp (the periodicity concept of yasp grid is currently not supported, use dune-spgrid) | ||
| 328 | // const auto periodic = getParamFromGroup<std::bitset<dim>>(modelParamGroup, "Grid.Periodic", std::bitset<dim>{}); | ||
| 329 | 101 | const std::bitset<dim> periodic; | |
| 330 | |||
| 331 | // Some sanity checks | ||
| 332 |
2/2✓ Branch 0 taken 180 times.
✓ Branch 1 taken 91 times.
|
303 | for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx) |
| 333 | { | ||
| 334 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 180 times.
|
202 | if (cells[dimIdx].size() + 1 != positions[dimIdx].size()) |
| 335 | { | ||
| 336 | ✗ | DUNE_THROW(Dune::RangeError, "Make sure to specify correct \"Cells\" and \"Positions\" arrays"); | |
| 337 | } | ||
| 338 |
1/2✓ Branch 0 taken 180 times.
✗ Branch 1 not taken.
|
202 | if (grading[dimIdx].size() + 1 != positions[dimIdx].size()) |
| 339 | { | ||
| 340 | ✗ | DUNE_THROW(Dune::RangeError, "Make sure to specify correct \"Grading\" and \"Positions\" arrays"); | |
| 341 | } | ||
| 342 | ctype temp = std::numeric_limits<ctype>::lowest(); | ||
| 343 |
2/2✓ Branch 0 taken 394 times.
✓ Branch 1 taken 180 times.
|
640 | for (unsigned int posIdx = 0; posIdx < positions[dimIdx].size(); ++posIdx) |
| 344 | { | ||
| 345 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 394 times.
|
438 | if (temp > positions[dimIdx][posIdx]) |
| 346 | { | ||
| 347 | ✗ | DUNE_THROW(Dune::RangeError, "Make sure to specify a monotone increasing \"Positions\" array"); | |
| 348 | } | ||
| 349 | 438 | temp = positions[dimIdx][posIdx]; | |
| 350 | } | ||
| 351 | } | ||
| 352 | |||
| 353 | 101 | const auto globalPositions = computeGlobalPositions_(positions, cells, grading, verbose); | |
| 354 | |||
| 355 | // make the grid | ||
| 356 |
2/4✓ Branch 1 taken 91 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 91 times.
✗ Branch 5 not taken.
|
202 | if (!hasParamInGroup(modelParamGroup, "Grid.Partitioning")) |
| 357 | { | ||
| 358 | // construct using default load balancing | ||
| 359 |
2/4✓ Branch 1 taken 91 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 91 times.
|
101 | ParentType::gridPtr() = std::make_shared<Grid>(globalPositions, periodic, overlap); |
| 360 | } | ||
| 361 | else | ||
| 362 | { | ||
| 363 | // construct using user defined partitioning | ||
| 364 | ✗ | const auto partitioning = getParamFromGroup<std::array<int, dim>>(modelParamGroup, "Grid.Partitioning"); | |
| 365 | ✗ | Dune::Yasp::FixedSizePartitioning<dim> lb(partitioning); | |
| 366 | ✗ | ParentType::gridPtr() = std::make_shared<Grid>(globalPositions, periodic, overlap, typename Grid::Communication(), &lb); | |
| 367 | ✗ | } | |
| 368 | |||
| 369 |
1/2✓ Branch 1 taken 91 times.
✗ Branch 2 not taken.
|
101 | postProcessing_(modelParamGroup); |
| 370 | 101 | } | |
| 371 | |||
| 372 | private: | ||
| 373 | /*! | ||
| 374 | * \brief Postprocessing for YaspGrid | ||
| 375 | */ | ||
| 376 | 101 | void postProcessing_(const std::string& modelParamGroup) | |
| 377 | { | ||
| 378 | // Check if should refine the grid | ||
| 379 | 101 | const bool keepPhysicalOverlap = getParamFromGroup<bool>(modelParamGroup, "Grid.KeepPhysicalOverlap", true); | |
| 380 | 101 | ParentType::grid().refineOptions(keepPhysicalOverlap); | |
| 381 | 101 | ParentType::maybeRefineGrid(modelParamGroup); | |
| 382 | 101 | ParentType::loadBalance(); | |
| 383 | 101 | } | |
| 384 | |||
| 385 | //! Compute the global position tensor grid from the given positions, cells, and grading factors | ||
| 386 | std::array<std::vector<ctype>, dim> | ||
| 387 | 101 | computeGlobalPositions_(const std::array<std::vector<ctype>, dim>& positions, | |
| 388 | const std::array<std::vector<int>, dim>& cells, | ||
| 389 | const std::array<std::vector<ctype>, dim>& grading, | ||
| 390 | bool verbose = false) | ||
| 391 | { | ||
| 392 | 101 | std::array<std::vector<ctype>, dim> globalPositions; | |
| 393 | using std::pow; | ||
| 394 | using Dune::power; | ||
| 395 |
2/2✓ Branch 0 taken 180 times.
✓ Branch 1 taken 91 times.
|
303 | for (int dimIdx = 0; dimIdx < dim; dimIdx++) |
| 396 | { | ||
| 397 |
2/2✓ Branch 0 taken 214 times.
✓ Branch 1 taken 180 times.
|
438 | for (int zoneIdx = 0; zoneIdx < cells[dimIdx].size(); ++zoneIdx) |
| 398 | { | ||
| 399 |
2/2✓ Branch 0 taken 75 times.
✓ Branch 1 taken 139 times.
|
236 | ctype lower = positions[dimIdx][zoneIdx]; |
| 400 | 236 | ctype upper = positions[dimIdx][zoneIdx+1]; | |
| 401 | 236 | int numCells = cells[dimIdx][zoneIdx]; | |
| 402 | 236 | ctype gradingFactor = grading[dimIdx][zoneIdx]; | |
| 403 | 236 | ctype length = upper - lower; | |
| 404 | 236 | ctype height = 1.0; | |
| 405 | 236 | bool increasingCellSize = false; | |
| 406 | |||
| 407 |
2/2✓ Branch 0 taken 75 times.
✓ Branch 1 taken 139 times.
|
236 | if (verbose) |
| 408 | { | ||
| 409 | 75 | std::cout << "dim " << dimIdx | |
| 410 |
4/8✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 75 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 75 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 75 times.
✗ Branch 11 not taken.
|
75 | << " lower " << lower |
| 411 |
2/4✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 75 times.
✗ Branch 5 not taken.
|
75 | << " upper " << upper |
| 412 | 75 | << " numCells " << numCells | |
| 413 |
3/6✓ Branch 1 taken 75 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 75 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 75 times.
✗ Branch 8 not taken.
|
75 | << " grading " << gradingFactor; |
| 414 | } | ||
| 415 | |||
| 416 |
2/2✓ Branch 0 taken 186 times.
✓ Branch 1 taken 28 times.
|
236 | if (gradingFactor > 1.0) |
| 417 | { | ||
| 418 | increasingCellSize = true; | ||
| 419 | } | ||
| 420 | |||
| 421 | // take absolute values and reverse cell size increment to achieve | ||
| 422 | // reverse behavior for negative values | ||
| 423 |
2/2✓ Branch 0 taken 22 times.
✓ Branch 1 taken 164 times.
|
208 | if (gradingFactor < 0.0) |
| 424 | { | ||
| 425 | using std::abs; | ||
| 426 | 22 | gradingFactor = abs(gradingFactor); | |
| 427 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 17 times.
|
22 | if (gradingFactor < 1.0) |
| 428 | { | ||
| 429 | 5 | increasingCellSize = true; | |
| 430 | } | ||
| 431 | } | ||
| 432 | |||
| 433 | // if the grading factor is exactly 1.0 do equal spacing | ||
| 434 |
4/4✓ Branch 0 taken 204 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 159 times.
✓ Branch 3 taken 45 times.
|
236 | if (gradingFactor > 1.0 - 1e-7 && gradingFactor < 1.0 + 1e-7) |
| 435 | { | ||
| 436 | 181 | height = 1.0 / numCells; | |
| 437 |
2/2✓ Branch 0 taken 36 times.
✓ Branch 1 taken 123 times.
|
181 | if (verbose) |
| 438 | { | ||
| 439 |
2/4✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
|
36 | std::cout << " -> h " << height * length << std::endl; |
| 440 | } | ||
| 441 | } | ||
| 442 | // if grading factor is not 1.0, do power law spacing | ||
| 443 | else | ||
| 444 | { | ||
| 445 | 55 | height = (1.0 - gradingFactor) / (1.0 - power(gradingFactor, numCells)); | |
| 446 | |||
| 447 |
2/2✓ Branch 0 taken 16 times.
✓ Branch 1 taken 39 times.
|
55 | if (verbose) |
| 448 | { | ||
| 449 |
2/4✓ Branch 1 taken 39 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 39 times.
✗ Branch 5 not taken.
|
39 | std::cout << " -> grading_eff " << gradingFactor |
| 450 |
2/4✓ Branch 1 taken 39 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 39 times.
✗ Branch 5 not taken.
|
39 | << " h_min " << height * power(gradingFactor, 0) * length |
| 451 |
2/4✓ Branch 1 taken 39 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 39 times.
✗ Branch 5 not taken.
|
78 | << " h_max " << height * power(gradingFactor, numCells-1) * length |
| 452 | 236 | << std::endl; | |
| 453 | } | ||
| 454 | } | ||
| 455 | |||
| 456 | 236 | std::vector<ctype> localPositions; | |
| 457 |
1/2✓ Branch 1 taken 214 times.
✗ Branch 2 not taken.
|
236 | localPositions.push_back(0); |
| 458 |
2/2✓ Branch 0 taken 10113 times.
✓ Branch 1 taken 214 times.
|
10415 | for (int i = 0; i < numCells-1; i++) |
| 459 | { | ||
| 460 | 10179 | ctype hI = height; | |
| 461 |
4/4✓ Branch 0 taken 9564 times.
✓ Branch 1 taken 549 times.
✓ Branch 2 taken 5220 times.
✓ Branch 3 taken 4344 times.
|
10179 | if (!(gradingFactor < 1.0 + 1e-7 && gradingFactor > 1.0 - 1e-7)) |
| 462 | { | ||
| 463 |
2/2✓ Branch 0 taken 4414 times.
✓ Branch 1 taken 1355 times.
|
5769 | if (increasingCellSize) |
| 464 | { | ||
| 465 | 4414 | hI *= power(gradingFactor, i); | |
| 466 | } | ||
| 467 | else | ||
| 468 | { | ||
| 469 | 2710 | hI *= power(gradingFactor, numCells-i-1); | |
| 470 | } | ||
| 471 | } | ||
| 472 |
1/4✓ Branch 1 taken 10113 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
10179 | localPositions.push_back(localPositions[i] + hI); |
| 473 | } | ||
| 474 | |||
| 475 |
2/2✓ Branch 0 taken 10327 times.
✓ Branch 1 taken 214 times.
|
10651 | for (int i = 0; i < localPositions.size(); i++) |
| 476 | { | ||
| 477 | 10415 | localPositions[i] *= length; | |
| 478 | 10415 | localPositions[i] += lower; | |
| 479 | } | ||
| 480 | |||
| 481 |
2/2✓ Branch 0 taken 10327 times.
✓ Branch 1 taken 214 times.
|
10651 | for (unsigned int i = 0; i < localPositions.size(); ++i) |
| 482 | { | ||
| 483 |
1/2✓ Branch 1 taken 10327 times.
✗ Branch 2 not taken.
|
10415 | globalPositions[dimIdx].push_back(localPositions[i]); |
| 484 | } | ||
| 485 | } | ||
| 486 |
1/2✓ Branch 1 taken 180 times.
✗ Branch 2 not taken.
|
202 | globalPositions[dimIdx].push_back(positions[dimIdx].back()); |
| 487 | } | ||
| 488 | |||
| 489 | 101 | return globalPositions; | |
| 490 | ✗ | } | |
| 491 | }; | ||
| 492 | |||
| 493 | namespace Grid::Capabilities { | ||
| 494 | |||
| 495 | // To the best of our knowledge YaspGrid is view thread-safe | ||
| 496 | // This specialization can be removed after we depend on Dune release 2.9 and this is guaranteed by YaspGrid itself | ||
| 497 | template<class Coordinates, int dim> | ||
| 498 | struct MultithreadingSupported<Dune::YaspGrid<dim, Coordinates>> | ||
| 499 | { | ||
| 500 | template<class GV> | ||
| 501 | static bool eval(const GV& gv) // default is independent of the grid view | ||
| 502 | { return true; } | ||
| 503 | }; | ||
| 504 | |||
| 505 | } // end namespace Grid::Capabilities | ||
| 506 | |||
| 507 | } // end namespace Dumux | ||
| 508 | |||
| 509 | #endif | ||
| 510 |