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