| 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 InputOutput | ||
| 10 | * \brief Grid manager specialization for SubGrid | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_IO_GRID_MANAGER_SUB_HH | ||
| 13 | #define DUMUX_IO_GRID_MANAGER_SUB_HH | ||
| 14 | |||
| 15 | #include <memory> | ||
| 16 | #include <utility> | ||
| 17 | #include <algorithm> | ||
| 18 | #include <stack> | ||
| 19 | #include <vector> | ||
| 20 | #include <numeric> | ||
| 21 | |||
| 22 | #include <dune/common/shared_ptr.hh> | ||
| 23 | #include <dune/common/concept.hh> | ||
| 24 | #include <dune/grid/yaspgrid.hh> | ||
| 25 | |||
| 26 | // SubGrid specific includes | ||
| 27 | #if HAVE_DUNE_SUBGRID | ||
| 28 | #include <dune/subgrid/subgrid.hh> | ||
| 29 | #include <dumux/io/rasterimagereader.hh> | ||
| 30 | #include <dumux/io/rasterimagewriter.hh> | ||
| 31 | #endif | ||
| 32 | |||
| 33 | #include <dumux/io/grid/gridmanager_base.hh> | ||
| 34 | |||
| 35 | #include <dumux/common/parameters.hh> | ||
| 36 | #include <dumux/common/boundaryflag.hh> | ||
| 37 | |||
| 38 | #if HAVE_DUNE_SUBGRID | ||
| 39 | namespace Dumux { | ||
| 40 | namespace Concept { | ||
| 41 | |||
| 42 | /*! | ||
| 43 | * \ingroup InputOutput | ||
| 44 | * \brief The element selector concept | ||
| 45 | */ | ||
| 46 | template<class Element> | ||
| 47 | struct ElementSelector | ||
| 48 | { | ||
| 49 | template<class F> | ||
| 50 | auto require(F&& f) -> decltype( | ||
| 51 | bool(f(std::declval<const Element&>())) | ||
| 52 | ); | ||
| 53 | }; | ||
| 54 | } // end namespace Concept | ||
| 55 | |||
| 56 | /*! | ||
| 57 | * \ingroup InputOutput | ||
| 58 | * \brief The base class for grid managers for dune-subgrid. | ||
| 59 | */ | ||
| 60 | template <class HostGrid, class HostGridManager = GridManager<HostGrid>> | ||
| 61 |
8/16✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
|
21 | class SubGridManagerBase |
| 62 | : public GridManagerBase<Dune::SubGrid<HostGrid::dimension, HostGrid>> | ||
| 63 | { | ||
| 64 | static constexpr int dim = HostGrid::dimension; | ||
| 65 | using HostElement = typename HostGrid::template Codim<0>::Entity; | ||
| 66 | using GlobalPosition = typename HostElement::Geometry::GlobalCoordinate; | ||
| 67 | |||
| 68 | public: | ||
| 69 | using Grid = Dune::SubGrid<dim, HostGrid>; | ||
| 70 | |||
| 71 | /*! | ||
| 72 | * \brief Make the grid using an externally created host grid. | ||
| 73 | */ | ||
| 74 | template<class ES, | ||
| 75 | typename std::enable_if_t<Dune::models<Concept::ElementSelector<HostElement>, ES>(), int> = 0> | ||
| 76 | ✗ | void init(HostGrid& hostGrid, | |
| 77 | const ES& selector, | ||
| 78 | const std::string& paramGroup = "") | ||
| 79 | { | ||
| 80 | ✗ | this->gridPtr() = std::make_unique<Grid>(hostGrid); | |
| 81 | ✗ | updateSubGrid_(selector); | |
| 82 | ✗ | loadBalance(); | |
| 83 | ✗ | } | |
| 84 | |||
| 85 | /*! | ||
| 86 | * \brief Make the grid and create the host grid internally. | ||
| 87 | */ | ||
| 88 | template<class ES, | ||
| 89 | typename std::enable_if_t<Dune::models<Concept::ElementSelector<HostElement>, ES>(), int> = 0> | ||
| 90 | 10 | void init(const ES& selector, | |
| 91 | const std::string& paramGroup = "") | ||
| 92 | { | ||
| 93 | 10 | initHostGrid_(paramGroup); | |
| 94 |
3/8✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 8 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
20 | this->gridPtr() = std::make_unique<Grid>(hostGridManager_->grid()); |
| 95 | 10 | updateSubGrid_(selector); | |
| 96 | 10 | loadBalance(); | |
| 97 | 10 | } | |
| 98 | |||
| 99 | /*! | ||
| 100 | * \brief Call loadBalance() function of the grid. | ||
| 101 | */ | ||
| 102 | 31 | void loadBalance() | |
| 103 | { | ||
| 104 |
2/4✗ Branch 1 not taken.
✓ Branch 2 taken 23 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 23 times.
|
31 | if (Dune::MPIHelper::getCommunication().size() > 1) |
| 105 | ✗ | this->grid().loadBalance(); | |
| 106 | 31 | } | |
| 107 | |||
| 108 | /*! | ||
| 109 | * \brief Update the existing subgrid. | ||
| 110 | */ | ||
| 111 | template<class ES, | ||
| 112 | typename std::enable_if_t<Dune::models<Concept::ElementSelector<HostElement>, ES>(), int> = 0> | ||
| 113 | void update(const ES& selector) | ||
| 114 | { | ||
| 115 | updateSubGrid_(selector); | ||
| 116 | loadBalance(); | ||
| 117 | } | ||
| 118 | |||
| 119 | protected: | ||
| 120 | |||
| 121 | /*! | ||
| 122 | * \brief Update the subgrid. | ||
| 123 | */ | ||
| 124 | template<class ES, | ||
| 125 | typename std::enable_if_t<Dune::models<Concept::ElementSelector<HostElement>, ES>(), int> = 0> | ||
| 126 | 37 | void updateSubGrid_(const ES& selector) | |
| 127 | { | ||
| 128 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
37 | auto& subgridPtr = this->gridPtr(); |
| 129 | |||
| 130 | // A container to store the host grid elements' ids. | ||
| 131 |
1/2✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
|
74 | std::set<typename HostGrid::Traits::GlobalIdSet::IdType> elementsForSubgrid; |
| 132 |
2/4✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
|
74 | const auto& globalIDset = subgridPtr->getHostGrid().globalIdSet(); |
| 133 | |||
| 134 | // Construct the subgrid. | ||
| 135 |
2/4✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
|
74 | subgridPtr->createBegin(); |
| 136 | |||
| 137 | // Loop over all elements of the host grid and use the selector to | ||
| 138 | // choose which elements to add to the subgrid. | ||
| 139 |
2/4✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 18 times.
✗ Branch 5 not taken.
|
74 | auto hostGridView = subgridPtr->getHostGrid().leafGridView(); |
| 140 |
12/16✓ Branch 1 taken 30849 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 86114 times.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 78331 times.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 768 times.
✓ Branch 10 taken 3 times.
✓ Branch 11 taken 768 times.
✓ Branch 12 taken 3 times.
✓ Branch 14 taken 768 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 768 times.
✗ Branch 18 not taken.
|
255231 | for (const auto& e : elements(hostGridView)) |
| 141 |
4/4✓ Branch 1 taken 82439 times.
✓ Branch 2 taken 4428 times.
✓ Branch 3 taken 50590 times.
✓ Branch 4 taken 28509 times.
|
153958 | if (selector(e)) |
| 142 |
2/4✓ Branch 1 taken 53930 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 53930 times.
✗ Branch 5 not taken.
|
113498 | elementsForSubgrid.insert(globalIDset.template id<0>(e)); |
| 143 | |||
| 144 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 22 times.
|
74 | if (elementsForSubgrid.empty()) |
| 145 | ✗ | DUNE_THROW(Dune::GridError, "No elements in subgrid"); | |
| 146 | |||
| 147 |
2/4✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
|
74 | subgridPtr->insertSetPartial(elementsForSubgrid); |
| 148 |
2/4✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
|
74 | subgridPtr->createEnd(); |
| 149 | 37 | } | |
| 150 | |||
| 151 | 11 | void initHostGrid_(const std::string& paramGroup) | |
| 152 | { | ||
| 153 |
2/4✗ Branch 1 not taken.
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 9 times.
|
11 | hostGridManager_ = std::make_unique<HostGridManager>(); |
| 154 | 22 | hostGridManager_->init(paramGroup); | |
| 155 | 11 | } | |
| 156 | |||
| 157 | 2 | void initHostGrid_(const GlobalPosition& lowerLeft, | |
| 158 | const GlobalPosition& upperRight, | ||
| 159 | const std::array<int, dim>& cells, | ||
| 160 | const std::string& paramGroup, | ||
| 161 | const int overlap = 1) | ||
| 162 | { | ||
| 163 |
2/4✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
|
2 | hostGridManager_ = std::make_unique<HostGridManager>(); |
| 164 | 6 | hostGridManager_->init(lowerLeft, upperRight, cells, paramGroup, overlap); | |
| 165 | 2 | } | |
| 166 | |||
| 167 | /*! | ||
| 168 | * \brief Returns a reference to the host grid. | ||
| 169 | */ | ||
| 170 | HostGrid& hostGrid_() | ||
| 171 | { | ||
| 172 |
0/2✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
40002 | return hostGridManager_->grid(); |
| 173 | } | ||
| 174 | |||
| 175 | std::unique_ptr<HostGridManager> hostGridManager_; | ||
| 176 | }; | ||
| 177 | |||
| 178 | /*! | ||
| 179 | * \ingroup InputOutput | ||
| 180 | * \brief Provides a grid manager for SubGrids | ||
| 181 | * from information in the input file | ||
| 182 | * | ||
| 183 | * The following keys are recognized: | ||
| 184 | * - All parameters that the host grid knows | ||
| 185 | */ | ||
| 186 | template<int dim, class HostGrid> | ||
| 187 |
11/19✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 1 times.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✓ Branch 23 taken 1 times.
✗ Branch 24 not taken.
|
12 | class GridManager<Dune::SubGrid<dim, HostGrid>> |
| 188 | : public SubGridManagerBase<HostGrid, GridManager<HostGrid>> | ||
| 189 | {}; | ||
| 190 | |||
| 191 | /*! | ||
| 192 | * \ingroup InputOutput | ||
| 193 | * \brief Provides a grid manager for SubGrids | ||
| 194 | * from information in the input file | ||
| 195 | * | ||
| 196 | * The following keys are recognized: | ||
| 197 | * - All parameters that the host grid knows | ||
| 198 | * - Image: the image file if the sub grid is constructed from a raster image | ||
| 199 | * (in that case the host grid has to be any 2D YaspGrid) | ||
| 200 | */ | ||
| 201 | template<int dim, class Coordinates> | ||
| 202 |
11/19✓ Branch 1 taken 3 times.
✓ Branch 2 taken 7 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 2 times.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✓ Branch 26 taken 1 times.
✗ Branch 27 not taken.
✓ Branch 31 taken 1 times.
✗ Branch 32 not taken.
|
28 | class GridManager<Dune::SubGrid<dim, Dune::YaspGrid<dim, Coordinates>>> |
| 203 | : public SubGridManagerBase<Dune::YaspGrid<dim, Coordinates>, GridManager<Dune::YaspGrid<dim, Coordinates>>> | ||
| 204 | { | ||
| 205 | using ParentType = SubGridManagerBase<Dune::YaspGrid<dim, Coordinates>, | ||
| 206 | GridManager<Dune::YaspGrid<dim, Coordinates>>>; | ||
| 207 | public: | ||
| 208 | using typename ParentType::Grid; | ||
| 209 | using ParentType::init; | ||
| 210 | |||
| 211 | /*! | ||
| 212 | * \brief Make the subgrid without host grid and element selector | ||
| 213 | * This means we try to construct the element selector from the input file | ||
| 214 | * \param paramGroup the parameter file group to check | ||
| 215 | */ | ||
| 216 | 2 | void init(const std::string& paramGroup = "") | |
| 217 | { | ||
| 218 | // check if there is an image file we can construct the element selector from | ||
| 219 |
6/14✓ 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 9 not taken.
✓ Branch 10 taken 2 times.
✓ Branch 11 taken 2 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
|
6 | if (hasParamInGroup(paramGroup, "Grid.Image")) |
| 220 | { | ||
| 221 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
4 | const auto imgFileName = getParamFromGroup<std::string>(paramGroup, "Grid.Image"); |
| 222 |
2/6✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
4 | const auto ext = this->getFileExtension(imgFileName); |
| 223 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | if (ext == "pbm") |
| 224 | { | ||
| 225 | if (dim != 2) | ||
| 226 | DUNE_THROW(Dune::GridError, "Portable Bitmap Format only supports dim == 2"); | ||
| 227 | |||
| 228 | // read image | ||
| 229 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
6 | const auto img = NetPBMReader::readPBM(imgFileName); |
| 230 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | createGridFromImage_(img, paramGroup); |
| 231 | } | ||
| 232 | else | ||
| 233 | ✗ | DUNE_THROW(Dune::IOError, "The SubGridManager doesn't support image files with extension: *." << ext); | |
| 234 | |||
| 235 | } | ||
| 236 | ✗ | else if (hasParamInGroup(paramGroup, "Grid.BinaryMask")) | |
| 237 | { | ||
| 238 | ✗ | const auto maskFileName = getParamFromGroup<std::string>(paramGroup, "Grid.BinaryMask"); | |
| 239 | ✗ | std::ifstream mask(maskFileName, std::ios_base::binary); | |
| 240 | ✗ | std::vector<char> buffer(std::istreambuf_iterator<char>(mask), std::istreambuf_iterator<char>{}); | |
| 241 | ✗ | const auto cells = getParamFromGroup<std::array<int, dim>>(paramGroup, "Grid.Cells"); | |
| 242 | ✗ | if (const auto c = std::accumulate(cells.begin(), cells.end(), 1, std::multiplies<int>{}); c != buffer.size()) | |
| 243 | ✗ | DUNE_THROW(Dune::IOError, "Grid dimensions doesn't match number of cells specified " << c << ":" << buffer.size()); | |
| 244 | |||
| 245 | ✗ | maybePostProcessBinaryMask_(buffer, cells, paramGroup); | |
| 246 | |||
| 247 | using GlobalPosition = typename ParentType::Grid::template Codim<0>::Geometry::GlobalCoordinate; | ||
| 248 | ✗ | const auto lowerLeft = GlobalPosition(0.0); | |
| 249 | ✗ | const auto upperRight = [&]{ | |
| 250 | ✗ | if (hasParamInGroup(paramGroup, "Grid.PixelDimensions")) | |
| 251 | { | ||
| 252 | ✗ | auto upperRight = getParamFromGroup<GlobalPosition>(paramGroup, "Grid.PixelDimensions"); | |
| 253 | ✗ | for (int i = 0; i < upperRight.size(); ++i) | |
| 254 | ✗ | upperRight[i] *= cells[i]; | |
| 255 | ✗ | upperRight += lowerLeft; | |
| 256 | ✗ | return upperRight; | |
| 257 | } | ||
| 258 | else | ||
| 259 | ✗ | return getParamFromGroup<GlobalPosition>(paramGroup, "Grid.UpperRight"); | |
| 260 | ✗ | }(); | |
| 261 | |||
| 262 | // construct the host grid | ||
| 263 | ✗ | this->initHostGrid_(lowerLeft, upperRight, cells, paramGroup); | |
| 264 | |||
| 265 | // check if the marker is customized, per default | ||
| 266 | // we use all cells that are encoded as 0 | ||
| 267 | ✗ | const char marker = getParamFromGroup<char>(paramGroup, "Grid.Marker", 0); | |
| 268 | ✗ | const auto elementSelector = [&](const auto& element) | |
| 269 | { | ||
| 270 | ✗ | auto level0 = element; | |
| 271 | ✗ | while(level0.hasFather()) | |
| 272 | ✗ | level0 = level0.father(); | |
| 273 | ✗ | const auto eIdx = this->hostGrid_().levelGridView(0).indexSet().index(level0); | |
| 274 | ✗ | return buffer[eIdx] != marker; | |
| 275 | }; | ||
| 276 | |||
| 277 | // create the grid | ||
| 278 | ✗ | this->gridPtr() = std::make_unique<Grid>(this->hostGrid_()); | |
| 279 | ✗ | this->updateSubGrid_(elementSelector); | |
| 280 | ✗ | this->loadBalance(); | |
| 281 | } | ||
| 282 | else | ||
| 283 | { | ||
| 284 | ✗ | std::cerr << "No element selector provided and Grid.Image or Grid.BinaryMask not found" << std::endl; | |
| 285 | ✗ | std::cerr << "Constructing sub-grid with all elements of the host grid" << std::endl; | |
| 286 | |||
| 287 | // create host grid | ||
| 288 | using GlobalPosition = typename ParentType::Grid::template Codim<0>::Geometry::GlobalCoordinate; | ||
| 289 | ✗ | const auto lowerLeft = getParamFromGroup<GlobalPosition>(paramGroup, "Grid.LowerLeft", GlobalPosition(0.0)); | |
| 290 | ✗ | const auto upperRight = getParamFromGroup<GlobalPosition>(paramGroup, "Grid.UpperRight"); | |
| 291 | ✗ | const auto cells = getParamFromGroup<std::array<int, dim>>(paramGroup, "Grid.Cells"); | |
| 292 | ✗ | this->initHostGrid_(lowerLeft, upperRight, cells, paramGroup); | |
| 293 | ✗ | const auto elementSelector = [](const auto& element){ return true; }; | |
| 294 | // create the grid | ||
| 295 | ✗ | this->gridPtr() = std::make_unique<Grid>(this->hostGrid_()); | |
| 296 | ✗ | this->updateSubGrid_(elementSelector); | |
| 297 | ✗ | this->loadBalance(); | |
| 298 | } | ||
| 299 | 2 | } | |
| 300 | |||
| 301 | private: | ||
| 302 | template<class Img> | ||
| 303 | 2 | void createGridFromImage_(const Img& img, const std::string& paramGroup) | |
| 304 | { | ||
| 305 | using GlobalPosition = typename ParentType::Grid::template Codim<0>::Geometry::GlobalCoordinate; | ||
| 306 |
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 9 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
6 | const bool repeated = hasParamInGroup(paramGroup,"Grid.Repeat"); |
| 307 | |||
| 308 | // get the number of cells | ||
| 309 | 2 | const std::array<int, dim> imageDimensions{static_cast<int>(img.header().nCols), static_cast<int>(img.header().nRows)}; | |
| 310 | 4 | std::array<int, dim> cells{imageDimensions[0], imageDimensions[1]}; | |
| 311 | |||
| 312 | 2 | std::array<int, dim> repeatsDefault; repeatsDefault.fill(1); | |
| 313 | 2 | const auto numRepeats = getParamFromGroup<std::array<int, dim>>(paramGroup, "Grid.Repeat", repeatsDefault); | |
| 314 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
|
6 | for (int i = 0; i < dim; i++) |
| 315 | 16 | cells[i] = cells[i] * numRepeats[i]; | |
| 316 | |||
| 317 | // get the corner coordinates | ||
| 318 | 2 | const auto [lowerLeft, upperRight] = [&]() | |
| 319 | { | ||
| 320 | 2 | const auto lowerLeft = getParamFromGroup<GlobalPosition>(paramGroup, "Grid.LowerLeft", GlobalPosition(0.0)); | |
| 321 |
7/12✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 times.
✓ Branch 9 taken 1 times.
✓ Branch 10 taken 1 times.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
10 | if (hasParamInGroup(paramGroup, "Grid.PixelDimensions")) |
| 322 | { | ||
| 323 | 1 | auto upperRight = getParamFromGroup<GlobalPosition>(paramGroup, "Grid.PixelDimensions"); | |
| 324 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
|
3 | for (int i = 0; i < upperRight.size(); ++i) |
| 325 | 6 | upperRight[i] *= cells[i]; | |
| 326 | 1 | upperRight += lowerLeft; | |
| 327 | 1 | return std::make_pair(lowerLeft, upperRight); | |
| 328 | } | ||
| 329 | else | ||
| 330 | 2 | return std::make_pair(lowerLeft, getParamFromGroup<GlobalPosition>(paramGroup, "Grid.UpperRight")); | |
| 331 |
1/2✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
8 | }(); |
| 332 | |||
| 333 | // construct the host grid | ||
| 334 | 2 | this->initHostGrid_(lowerLeft, upperRight, cells, paramGroup); | |
| 335 | |||
| 336 | // check if the marker is customized, per default | ||
| 337 | // we mark all cells that are encoded as 0 | ||
| 338 | 2 | const bool marked = getParamFromGroup<bool>(paramGroup, "Grid.Marker", false); | |
| 339 | |||
| 340 | // Create the element selector for a single image | ||
| 341 | 5002 | const auto elementSelector = [&](const auto& element) | |
| 342 | { | ||
| 343 |
3/6✗ Branch 1 not taken.
✓ Branch 2 taken 2500 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2500 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2500 times.
|
2500 | auto eIdx = this->hostGrid_().leafGridView().indexSet().index(element); |
| 344 | |||
| 345 | // if the hostgrid was refined, get the index of the original, un-refined | ||
| 346 | // host grid element which fits with the image's indices | ||
| 347 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 2500 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2500 times.
|
5000 | if (element.hasFather()) |
| 348 | { | ||
| 349 | ✗ | auto level0Element = element.father(); | |
| 350 | ✗ | while (level0Element.hasFather()) | |
| 351 | ✗ | level0Element = level0Element.father(); | |
| 352 | |||
| 353 | ✗ | assert(level0Element.level() == 0); | |
| 354 | ✗ | eIdx = this->hostGrid_().levelGridView(0).indexSet().index(level0Element); | |
| 355 | } | ||
| 356 | 7500 | return img[eIdx] == marked; | |
| 357 | }; | ||
| 358 | |||
| 359 | // Create the element selector for a repeated image | ||
| 360 | 75002 | const auto repeatedElementSelector = [&](const auto& element) | |
| 361 | { | ||
| 362 |
3/6✗ Branch 1 not taken.
✓ Branch 2 taken 37500 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 37500 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 37500 times.
|
37500 | auto eIdx = this->hostGrid_().leafGridView().indexSet().index(element); |
| 363 | |||
| 364 | // if the hostgrid was refined, get the index of the original, un-refined | ||
| 365 | // host grid element which fits with the image's indices | ||
| 366 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 37500 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 37500 times.
|
75000 | if (element.hasFather()) |
| 367 | { | ||
| 368 | ✗ | auto level0Element = element.father(); | |
| 369 | ✗ | while (level0Element.hasFather()) | |
| 370 | ✗ | level0Element = level0Element.father(); | |
| 371 | |||
| 372 | ✗ | assert(level0Element.level() == 0); | |
| 373 | ✗ | eIdx = this->hostGrid_().levelGridView(0).indexSet().index(level0Element); | |
| 374 | } | ||
| 375 | |||
| 376 | // figure out the size of the repeat, and the size of the target repeated grid | ||
| 377 | 75000 | const int numCols = imageDimensions[0]; | |
| 378 | 75000 | const int numRows = imageDimensions[1]; | |
| 379 | |||
| 380 | // map the eIdx to the original img index | ||
| 381 | 112500 | const int repeatUnitIndex = eIdx % (numCols * numRepeats[0] * numRows); | |
| 382 | 37500 | const int imgI = repeatUnitIndex % numCols; | |
| 383 | 75000 | const int imgJ = repeatUnitIndex / (numCols * numRepeats[0]); | |
| 384 | 112500 | return img[ (imgJ * numCols + imgI) ] == marked; | |
| 385 | }; | ||
| 386 | |||
| 387 | // create the grid | ||
| 388 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
|
2 | if (repeated) |
| 389 | { | ||
| 390 |
3/8✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
1 | this->gridPtr() = std::make_unique<Grid>(this->hostGrid_()); |
| 391 | 1 | this->updateSubGrid_(repeatedElementSelector); | |
| 392 | } | ||
| 393 | else | ||
| 394 | { | ||
| 395 |
3/8✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
1 | this->gridPtr() = std::make_unique<Grid>(this->hostGrid_()); |
| 396 | 1 | this->updateSubGrid_(elementSelector); | |
| 397 | } | ||
| 398 | |||
| 399 | 2 | this->loadBalance(); | |
| 400 | 2 | } | |
| 401 | |||
| 402 | private: | ||
| 403 | ✗ | void maybePostProcessBinaryMask_(std::vector<char>& mask, const std::array<int, dim>& cells, const std::string& paramGroup) const | |
| 404 | { | ||
| 405 | // implements pruning for directional connectivity scenario (e.g. flow from left to right) | ||
| 406 | // all cells that don't connect the boundaries in X (or Y or Z) direction are removed | ||
| 407 | ✗ | if (hasParamInGroup(paramGroup, "Grid.KeepOnlyConnected")) | |
| 408 | { | ||
| 409 | ✗ | std::vector<int> marker(mask.size(), 0); | |
| 410 | ✗ | const std::string direction = getParamFromGroup<std::string>(paramGroup, "Grid.KeepOnlyConnected"); | |
| 411 | if constexpr (dim == 3) | ||
| 412 | { | ||
| 413 | if (direction == "Y") | ||
| 414 | { | ||
| 415 | std::cout << "Keeping only cells connected to both boundaries in y-direction" << std::endl; | ||
| 416 | flood_(mask, marker, cells, 0, 0, 0, cells[0], 1, cells[2], 1); | ||
| 417 | flood_(mask, marker, cells, 0, cells[1]-1, 0, cells[0], cells[1], cells[2], 2); | ||
| 418 | } | ||
| 419 | else if (direction == "Z") | ||
| 420 | { | ||
| 421 | std::cout << "Keeping only cells connected to both boundaries in z-direction" << std::endl; | ||
| 422 | flood_(mask, marker, cells, 0, 0, 0, cells[0], cells[1], 1, 1); | ||
| 423 | flood_(mask, marker, cells, 0, 0, cells[2]-1, cells[0], cells[1], cells[2], 2); | ||
| 424 | } | ||
| 425 | else | ||
| 426 | { | ||
| 427 | std::cout << "Keeping only cells connected to both boundaries in x-direction" << std::endl; | ||
| 428 | flood_(mask, marker, cells, 0, 0, 0, 1, cells[1], cells[2], 1); | ||
| 429 | flood_(mask, marker, cells, cells[0]-1, 0, 0, cells[0], cells[1], cells[2], 2); | ||
| 430 | } | ||
| 431 | |||
| 432 | for (unsigned int i = 0; i < cells[0]; ++i) | ||
| 433 | for (unsigned int j = 0; j < cells[1]; ++j) | ||
| 434 | for (unsigned int k = 0; k < cells[2]; ++k) | ||
| 435 | if (const auto ijk = getIJK_(i, j, k, cells); marker[ijk] < 2) | ||
| 436 | mask[ijk] = false; | ||
| 437 | } | ||
| 438 | |||
| 439 | else if constexpr (dim == 2) | ||
| 440 | { | ||
| 441 | ✗ | if (direction == "Y") | |
| 442 | { | ||
| 443 | ✗ | std::cout << "Keeping only cells connected to both boundaries in y-direction" << std::endl; | |
| 444 | ✗ | flood_(mask, marker, cells, 0, 0, cells[0], 1, 1); | |
| 445 | ✗ | flood_(mask, marker, cells, 0, cells[1]-1, cells[0], cells[1], 2); | |
| 446 | } | ||
| 447 | else | ||
| 448 | { | ||
| 449 | ✗ | std::cout << "Keeping only cells connected to both boundaries in x-direction" << std::endl; | |
| 450 | ✗ | flood_(mask, marker, cells, 0, 0, 1, cells[1], 1); | |
| 451 | ✗ | flood_(mask, marker, cells, cells[0]-1, 0, cells[0], cells[1], 2); | |
| 452 | } | ||
| 453 | |||
| 454 | ✗ | for (unsigned int i = 0; i < cells[0]; ++i) | |
| 455 | ✗ | for (unsigned int j = 0; j < cells[1]; ++j) | |
| 456 | ✗ | if (const auto ij = getIJ_(i, j, cells); marker[ij] < 2) | |
| 457 | ✗ | mask[ij] = false; | |
| 458 | } | ||
| 459 | else | ||
| 460 | DUNE_THROW(Dune::NotImplemented, "KeepOnlyConnected only implemented for 2D and 3D"); | ||
| 461 | } | ||
| 462 | ✗ | } | |
| 463 | |||
| 464 | void flood_(std::vector<char>& mask, std::vector<int>& markers, const std::array<int, 3>& cells, | ||
| 465 | unsigned int iMin, unsigned int jMin, unsigned int kMin, | ||
| 466 | unsigned int iMax, unsigned int jMax, unsigned int kMax, | ||
| 467 | int marker) const | ||
| 468 | { | ||
| 469 | // flood-fill with marker starting from all seed cells in the range | ||
| 470 | for (unsigned int i = iMin; i < iMax; ++i) | ||
| 471 | for (unsigned int j = jMin; j < jMax; ++j) | ||
| 472 | for (unsigned int k = kMin; k < kMax; ++k) | ||
| 473 | fill_(mask, markers, cells, i, j, k, marker); | ||
| 474 | } | ||
| 475 | |||
| 476 | void flood_(std::vector<char>& mask, std::vector<int>& markers, const std::array<int, 2>& cells, | ||
| 477 | unsigned int iMin, unsigned int jMin, | ||
| 478 | unsigned int iMax, unsigned int jMax, | ||
| 479 | int marker) const | ||
| 480 | { | ||
| 481 | // flood-fill with marker starting from all seed cells in the range | ||
| 482 | ✗ | for (unsigned int i = iMin; i < iMax; ++i) | |
| 483 | ✗ | for (unsigned int j = jMin; j < jMax; ++j) | |
| 484 | ✗ | fill_(mask, markers, cells, i, j, marker); | |
| 485 | } | ||
| 486 | |||
| 487 | void fill_(std::vector<char>& mask, std::vector<int>& markers, const std::array<int, 3>& cells, | ||
| 488 | unsigned int i, unsigned int j, unsigned int k, int marker) const | ||
| 489 | { | ||
| 490 | std::stack<std::tuple<unsigned int, unsigned int, unsigned int>> st; | ||
| 491 | st.push(std::make_tuple(i, j, k)); | ||
| 492 | while(!st.empty()) | ||
| 493 | { | ||
| 494 | std::tie(i, j, k) = st.top(); | ||
| 495 | st.pop(); | ||
| 496 | |||
| 497 | if (i >= cells[0] || j >= cells[1] || k >= cells[2]) | ||
| 498 | continue; | ||
| 499 | |||
| 500 | const auto ijk = getIJK_(i, j, k, cells); | ||
| 501 | if (!mask[ijk] || markers[ijk] >= marker) | ||
| 502 | continue; | ||
| 503 | |||
| 504 | markers[ijk] += 1; | ||
| 505 | |||
| 506 | // we rely on -1 wrapping over for unsigned int = 0 | ||
| 507 | st.push(std::make_tuple(i+1, j, k)); | ||
| 508 | st.push(std::make_tuple(i-1, j, k)); | ||
| 509 | st.push(std::make_tuple(i, j+1, k)); | ||
| 510 | st.push(std::make_tuple(i, j-1, k)); | ||
| 511 | st.push(std::make_tuple(i, j, k+1)); | ||
| 512 | st.push(std::make_tuple(i, j, k-1)); | ||
| 513 | } | ||
| 514 | } | ||
| 515 | |||
| 516 | ✗ | void fill_(std::vector<char>& mask, std::vector<int>& markers, const std::array<int, 2>& cells, | |
| 517 | unsigned int i, unsigned int j, int marker) const | ||
| 518 | { | ||
| 519 | ✗ | std::stack<std::tuple<unsigned int, unsigned int>> st; | |
| 520 | ✗ | st.push(std::make_tuple(i, j)); | |
| 521 | ✗ | while(!st.empty()) | |
| 522 | { | ||
| 523 | ✗ | std::tie(i, j) = st.top(); | |
| 524 | ✗ | st.pop(); | |
| 525 | |||
| 526 | ✗ | if (i >= cells[0] || j >= cells[1]) | |
| 527 | continue; | ||
| 528 | |||
| 529 | ✗ | const auto ij = getIJ_(i, j, cells); | |
| 530 | ✗ | if (!mask[ij] || markers[ij] >= marker) | |
| 531 | continue; | ||
| 532 | |||
| 533 | ✗ | markers[ij] += 1; | |
| 534 | |||
| 535 | // we rely on -1 wrapping over for unsigned int = 0 | ||
| 536 | ✗ | st.push(std::make_tuple(i+1, j)); | |
| 537 | ✗ | st.push(std::make_tuple(i-1, j)); | |
| 538 | ✗ | st.push(std::make_tuple(i, j+1)); | |
| 539 | ✗ | st.push(std::make_tuple(i, j-1)); | |
| 540 | } | ||
| 541 | ✗ | } | |
| 542 | |||
| 543 | int getIJK_(int i, int j, int k, const std::array<int, 3>& cells) const | ||
| 544 | { return i + j*cells[0] + k*cells[0]*cells[1]; } | ||
| 545 | |||
| 546 | ✗ | int getIJ_(int i, int j, const std::array<int, 2>& cells) const | |
| 547 | ✗ | { return i + j*cells[0]; } | |
| 548 | }; | ||
| 549 | |||
| 550 | //! dune-subgrid doesn't have this implemented | ||
| 551 | template<int dim, class HostGrid> | ||
| 552 | class BoundaryFlag<Dune::SubGrid<dim, HostGrid>> | ||
| 553 | { | ||
| 554 | public: | ||
| 555 | BoundaryFlag() : flag_(-1) {} | ||
| 556 | |||
| 557 | template<class Intersection> | ||
| 558 | 1831792 | BoundaryFlag(const Intersection& i) : flag_(-1) {} | |
| 559 | |||
| 560 | using value_type = int; | ||
| 561 | |||
| 562 | value_type get() const | ||
| 563 | { DUNE_THROW(Dune::NotImplemented, "Sub-grid doesn't implement boundary segment indices!"); } | ||
| 564 | |||
| 565 | private: | ||
| 566 | int flag_; | ||
| 567 | }; | ||
| 568 | |||
| 569 | } // end namespace Dumux | ||
| 570 | #endif // HAVE_DUNE_SUBGRID | ||
| 571 | #endif | ||
| 572 |