| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // | ||
| 3 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 4 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 5 | // | ||
| 6 | /*! | ||
| 7 | * \file | ||
| 8 | * \ingroup Core | ||
| 9 | * \brief defines intersection mappers. | ||
| 10 | */ | ||
| 11 | #ifndef DUMUX_INTERSECTIONITERATOR_HH | ||
| 12 | #define DUMUX_INTERSECTIONITERATOR_HH | ||
| 13 | |||
| 14 | #include <vector> | ||
| 15 | #include <unordered_map> | ||
| 16 | |||
| 17 | #include <dune/grid/common/mcmgmapper.hh> | ||
| 18 | #include <dune/grid/common/rangegenerators.hh> | ||
| 19 | |||
| 20 | namespace Dumux { | ||
| 21 | |||
| 22 | /*! | ||
| 23 | * \ingroup Core | ||
| 24 | * \brief defines a standard intersection mapper for mapping of global DOFs assigned | ||
| 25 | * to faces. It only works for conforming grids, without hanging nodes. | ||
| 26 | */ | ||
| 27 | template<class GridView> | ||
| 28 | ✗ | class ConformingGridIntersectionMapper | |
| 29 | { | ||
| 30 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 31 | using GridIndexType = typename GridView::IndexSet::IndexType; | ||
| 32 | |||
| 33 | static constexpr int codimIntersection = 1; | ||
| 34 | public: | ||
| 35 | |||
| 36 | 116 | ConformingGridIntersectionMapper(const GridView& gridView) | |
| 37 | : gridView_(gridView) | ||
| 38 |
4/7✗ Branch 0 not taken.
✓ Branch 1 taken 114 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 108 times.
✓ Branch 4 taken 6 times.
✓ Branch 5 taken 7 times.
✗ Branch 6 not taken.
|
230 | , indexSet_(&gridView.indexSet()) |
| 39 | {} | ||
| 40 | |||
| 41 | void update (const GridView& gridView) | ||
| 42 | { | ||
| 43 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
|
62 | gridView_ = gridView; |
| 44 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
|
64 | indexSet_ = &gridView_.indexSet(); |
| 45 | } | ||
| 46 | |||
| 47 | void update (GridView&& gridView) | ||
| 48 | { | ||
| 49 | gridView_ = std::move(gridView); | ||
| 50 | indexSet_ = &gridView_.indexSet(); | ||
| 51 | } | ||
| 52 | |||
| 53 | //! The total number of intersections | ||
| 54 | std::size_t numIntersections() const | ||
| 55 | { | ||
| 56 | return gridView_.size(1); | ||
| 57 | } | ||
| 58 | |||
| 59 | /*! | ||
| 60 | * \brief The number of faces for a given element | ||
| 61 | * | ||
| 62 | * \param element The element | ||
| 63 | */ | ||
| 64 | ✗ | std::size_t numFaces(const Element& element) const | |
| 65 | { | ||
| 66 |
2/4✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 25 times.
✗ Branch 5 not taken.
|
179500 | return element.subEntities(1); |
| 67 | } | ||
| 68 | |||
| 69 | ✗ | GridIndexType globalIntersectionIndex(const Element& element, const std::size_t localFaceIdx) const | |
| 70 | { | ||
| 71 |
8/15✓ Branch 1 taken 18 times.
✓ Branch 2 taken 672630 times.
✓ Branch 3 taken 656012 times.
✓ Branch 4 taken 6452 times.
✓ Branch 5 taken 7268 times.
✓ Branch 6 taken 848 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 14 taken 839328 times.
✗ Branch 15 not taken.
✓ Branch 18 taken 68160 times.
✗ Branch 19 not taken.
|
9586512 | return indexSet_->subIndex(element, localFaceIdx, codimIntersection); |
| 72 | } | ||
| 73 | |||
| 74 | private: | ||
| 75 | GridView gridView_; | ||
| 76 | const typename GridView::IndexSet* indexSet_; | ||
| 77 | }; | ||
| 78 | |||
| 79 | /*! | ||
| 80 | * \ingroup Core | ||
| 81 | * \brief defines an intersection mapper for mapping of global DOFs assigned | ||
| 82 | * to faces which also works for non-conforming grids and corner-point grids. | ||
| 83 | */ | ||
| 84 | template<class GridView> | ||
| 85 | class NonConformingGridIntersectionMapper | ||
| 86 | { | ||
| 87 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 88 | using Intersection = typename GridView::Intersection; | ||
| 89 | using GridIndexType = typename GridView::IndexSet::IndexType; | ||
| 90 | |||
| 91 | public: | ||
| 92 | NonConformingGridIntersectionMapper(const GridView& gridview) | ||
| 93 | : gridView_(gridview), | ||
| 94 | indexSet_(&gridView_.indexSet()), | ||
| 95 | numIntersections_(gridView_.size(1)), | ||
| 96 | intersectionMapGlobal_(gridView_.size(0)) | ||
| 97 | { | ||
| 98 | |||
| 99 | } | ||
| 100 | |||
| 101 | //! The total number of intersections | ||
| 102 | std::size_t numIntersections() const | ||
| 103 | { | ||
| 104 | return numIntersections_; | ||
| 105 | } | ||
| 106 | |||
| 107 | GridIndexType globalIntersectionIndex(const Element& element, const std::size_t localFaceIdx) const | ||
| 108 | { | ||
| 109 | return (intersectionMapGlobal_[index(element)].find(localFaceIdx))->second; //use find() for const function! | ||
| 110 | } | ||
| 111 | |||
| 112 | std::size_t numFaces(const Element& element) | ||
| 113 | { | ||
| 114 | return intersectionMapGlobal_[index(element)].size(); | ||
| 115 | } | ||
| 116 | |||
| 117 | void update (const GridView& gridView) | ||
| 118 | { | ||
| 119 | gridView_ = gridView; | ||
| 120 | indexSet_ = &gridView_.indexSet(); | ||
| 121 | update_(); | ||
| 122 | } | ||
| 123 | |||
| 124 | void update (GridView&& gridView) | ||
| 125 | { | ||
| 126 | gridView_ = std::move(gridView); | ||
| 127 | indexSet_ = &gridView_.indexSet(); | ||
| 128 | update_(); | ||
| 129 | } | ||
| 130 | |||
| 131 | private: | ||
| 132 | GridIndexType index(const Element& element) const | ||
| 133 | { | ||
| 134 | return indexSet_->index(element); | ||
| 135 | } | ||
| 136 | |||
| 137 | void update_() | ||
| 138 | { | ||
| 139 | intersectionMapGlobal_.clear(); | ||
| 140 | intersectionMapGlobal_.resize(gridView_.size(0)); | ||
| 141 | |||
| 142 | int globalIntersectionIdx = 0; | ||
| 143 | for (const auto& element : elements(gridView_)) | ||
| 144 | { | ||
| 145 | int eIdx = index(element); | ||
| 146 | int fIdx = 0; | ||
| 147 | |||
| 148 | // run through all intersections with neighbors | ||
| 149 | for (const auto& intersection : intersections(gridView_, element)) | ||
| 150 | { | ||
| 151 | if (intersection.neighbor()) | ||
| 152 | { | ||
| 153 | auto neighbor = intersection.outside(); | ||
| 154 | int eIdxN = index(neighbor); | ||
| 155 | |||
| 156 | if (element.level() > neighbor.level() || (element.level() == neighbor.level() && eIdx < eIdxN)) | ||
| 157 | { | ||
| 158 | int fIdxN = 0; | ||
| 159 | for (const auto& intersectionNeighbor : intersections(gridView_, neighbor)) | ||
| 160 | { | ||
| 161 | if (intersectionNeighbor.neighbor()) | ||
| 162 | { | ||
| 163 | if (intersectionNeighbor.outside() == element) | ||
| 164 | { | ||
| 165 | break; | ||
| 166 | } | ||
| 167 | } | ||
| 168 | fIdxN++; | ||
| 169 | } | ||
| 170 | intersectionMapGlobal_[eIdx][fIdx] = globalIntersectionIdx; | ||
| 171 | intersectionMapGlobal_[eIdxN][fIdxN] = globalIntersectionIdx; | ||
| 172 | globalIntersectionIdx++; | ||
| 173 | } | ||
| 174 | } | ||
| 175 | else | ||
| 176 | { | ||
| 177 | intersectionMapGlobal_[eIdx][fIdx] = globalIntersectionIdx; | ||
| 178 | globalIntersectionIdx++; | ||
| 179 | } | ||
| 180 | |||
| 181 | fIdx++; | ||
| 182 | } | ||
| 183 | } | ||
| 184 | numIntersections_ = globalIntersectionIdx; | ||
| 185 | } | ||
| 186 | |||
| 187 | GridView gridView_; | ||
| 188 | const typename GridView::IndexSet* indexSet_; | ||
| 189 | unsigned int numIntersections_; | ||
| 190 | std::vector<std::unordered_map<int, int> > intersectionMapGlobal_; | ||
| 191 | }; | ||
| 192 | |||
| 193 | /*! | ||
| 194 | * \ingroup Core | ||
| 195 | * \brief defines an intersection mapper for mapping of global DOFs assigned | ||
| 196 | * to faces which also works for adaptive grids. | ||
| 197 | */ | ||
| 198 | template<class GridView> | ||
| 199 | class IntersectionMapper | ||
| 200 | { | ||
| 201 | using Grid = typename GridView::Grid; | ||
| 202 | enum {dim=Grid::dimension}; | ||
| 203 | using Element = typename Grid::template Codim<0>::Entity; | ||
| 204 | using Intersection = typename GridView::Intersection; | ||
| 205 | using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>; | ||
| 206 | |||
| 207 | public: | ||
| 208 | IntersectionMapper(const GridView& gridview) | ||
| 209 | : gridView_(gridview) | ||
| 210 | , elementMapper_(gridView_, Dune::mcmgElementLayout()) | ||
| 211 | , size_(gridView_.size(1)) | ||
| 212 | , intersectionMapGlobal_(gridView_.size(0)) | ||
| 213 | , intersectionMapLocal_(gridView_.size(0)) | ||
| 214 | { | ||
| 215 | const auto element = *gridView_.template begin<0>(); | ||
| 216 | |||
| 217 | int fIdx = 0; | ||
| 218 | for (const auto& intersection : intersections(gridView_, element)) | ||
| 219 | { | ||
| 220 | int idxInInside = intersection.indexInInside(); | ||
| 221 | |||
| 222 | standardLocalIdxMap_[idxInInside] = fIdx; | ||
| 223 | |||
| 224 | fIdx++; | ||
| 225 | } | ||
| 226 | } | ||
| 227 | |||
| 228 | const ElementMapper& elementMapper() const | ||
| 229 | { | ||
| 230 | return elementMapper_; | ||
| 231 | } | ||
| 232 | |||
| 233 | /*! | ||
| 234 | * \brief Map element to array index. | ||
| 235 | * | ||
| 236 | * \param element Reference to element which should be mapped. | ||
| 237 | * \return An index in the range 0 ... Max number of elements in set - 1. | ||
| 238 | */ | ||
| 239 | int index(const Element& element) const | ||
| 240 | { | ||
| 241 | return elementMapper_.index(element); | ||
| 242 | } | ||
| 243 | |||
| 244 | /*! | ||
| 245 | * \brief Map interface fIdx'th interface of element index to array index. | ||
| 246 | * | ||
| 247 | * \param elemIdx Index of element. | ||
| 248 | * \param fIdx Index of interface to map. | ||
| 249 | * \return An index in the range 0 ... Max number of intersections in set - 1. | ||
| 250 | */ | ||
| 251 | int subIndex(int elemIdx, int fIdx) | ||
| 252 | { | ||
| 253 | return intersectionMapGlobal_[elemIdx][fIdx]; | ||
| 254 | } | ||
| 255 | |||
| 256 | /*! | ||
| 257 | * \brief Map interface fIdx'th interface of element index to array index. | ||
| 258 | * | ||
| 259 | * \param elemIdx Index of element. | ||
| 260 | * \param fIdx Index of interface to map. | ||
| 261 | * \return An index in the range 0 ... Max number of intersections in set - 1. | ||
| 262 | */ | ||
| 263 | int subIndex(int elemIdx, int fIdx) const | ||
| 264 | { | ||
| 265 | return (intersectionMapGlobal_[elemIdx].find(fIdx))->second;//use find() for const function! | ||
| 266 | } | ||
| 267 | |||
| 268 | /*! | ||
| 269 | * \brief Map interface fIdx'th interface of element to array index. | ||
| 270 | * | ||
| 271 | * \param element Reference to element. | ||
| 272 | * \param fIdx Index of interface to map. | ||
| 273 | * \return An index in the range 0 ... Max number of intersections in set - 1. | ||
| 274 | */ | ||
| 275 | int subIndex(const Element& element, int fIdx) | ||
| 276 | { | ||
| 277 | return intersectionMapGlobal_[index(element)][fIdx]; | ||
| 278 | } | ||
| 279 | |||
| 280 | int subIndex(const Element& element, int fIdx) const | ||
| 281 | { | ||
| 282 | return intersectionMapGlobal_[index(element)].find(fIdx)->second;//use find() for const function! | ||
| 283 | } | ||
| 284 | |||
| 285 | int maplocal(int elemIdx, int fIdx) | ||
| 286 | { | ||
| 287 | return intersectionMapLocal_[elemIdx][fIdx]; | ||
| 288 | } | ||
| 289 | |||
| 290 | int maplocal(int elemIdx, int fIdx) const | ||
| 291 | { | ||
| 292 | return (intersectionMapLocal_[elemIdx].find(fIdx))->second;//use find() for const function! | ||
| 293 | } | ||
| 294 | |||
| 295 | |||
| 296 | int maplocal(const Element& element, int fIdx) | ||
| 297 | { | ||
| 298 | return intersectionMapLocal_[index(element)][fIdx]; | ||
| 299 | } | ||
| 300 | |||
| 301 | int maplocal(const Element& element, int fIdx) const | ||
| 302 | { | ||
| 303 | return (intersectionMapLocal_[index(element)].find(fIdx))->second;//use find() for const function! | ||
| 304 | } | ||
| 305 | |||
| 306 | // return number intersections | ||
| 307 | unsigned int size () const | ||
| 308 | { | ||
| 309 | return size_; | ||
| 310 | } | ||
| 311 | |||
| 312 | unsigned int size (int elemIdx) const | ||
| 313 | { | ||
| 314 | return intersectionMapLocal_[elemIdx].size(); | ||
| 315 | } | ||
| 316 | |||
| 317 | unsigned int size (const Element& element) const | ||
| 318 | { | ||
| 319 | return intersectionMapLocal_[index(element)].size(); | ||
| 320 | } | ||
| 321 | |||
| 322 | void update (const GridView &gridView) | ||
| 323 | { | ||
| 324 | gridView_ = gridView; | ||
| 325 | update_(); | ||
| 326 | } | ||
| 327 | |||
| 328 | void update (GridView&& gridView) | ||
| 329 | { | ||
| 330 | gridView_ = std::move(gridView); | ||
| 331 | update_(); | ||
| 332 | } | ||
| 333 | |||
| 334 | protected: | ||
| 335 | void update_() | ||
| 336 | { | ||
| 337 | elementMapper_.update(gridView_); | ||
| 338 | |||
| 339 | intersectionMapGlobal_.clear(); | ||
| 340 | intersectionMapGlobal_.resize(elementMapper_.size()); | ||
| 341 | intersectionMapLocal_.clear(); | ||
| 342 | intersectionMapLocal_.resize(elementMapper_.size()); | ||
| 343 | |||
| 344 | for (const auto& element : elements(gridView_)) | ||
| 345 | { | ||
| 346 | int eIdxGlobal = index(element); | ||
| 347 | |||
| 348 | int fIdx = 0; | ||
| 349 | // run through all intersections with neighbors | ||
| 350 | for (const auto& intersection : intersections(gridView_, element)) | ||
| 351 | { | ||
| 352 | int indexInInside = intersection.indexInInside(); | ||
| 353 | intersectionMapLocal_[eIdxGlobal][fIdx] = indexInInside; | ||
| 354 | |||
| 355 | fIdx++; | ||
| 356 | } | ||
| 357 | } | ||
| 358 | int globalIntersectionIdx = 0; | ||
| 359 | for (const auto& element : elements(gridView_)) | ||
| 360 | { | ||
| 361 | int eIdxGlobal = index(element); | ||
| 362 | |||
| 363 | int fIdx = 0; | ||
| 364 | // run through all intersections with neighbors | ||
| 365 | for (const auto& intersection : intersections(gridView_, element)) | ||
| 366 | { | ||
| 367 | if (intersection.neighbor()) | ||
| 368 | { | ||
| 369 | auto neighbor = intersection.outside(); | ||
| 370 | int globalIdxNeighbor = index(neighbor); | ||
| 371 | |||
| 372 | if (element.level() > neighbor.level() || (element.level() == neighbor.level() && eIdxGlobal < globalIdxNeighbor)) | ||
| 373 | { | ||
| 374 | |||
| 375 | int faceIdxNeighbor = 0; | ||
| 376 | if (size(globalIdxNeighbor) == 2 * dim) | ||
| 377 | { | ||
| 378 | faceIdxNeighbor = standardLocalIdxMap_[intersection.indexInOutside()]; | ||
| 379 | } | ||
| 380 | else | ||
| 381 | { | ||
| 382 | for (const auto& intersectionNeighbor : intersections(gridView_, neighbor)) | ||
| 383 | { | ||
| 384 | if (intersectionNeighbor.neighbor()) | ||
| 385 | { | ||
| 386 | if (intersectionNeighbor.outside() == element) | ||
| 387 | { | ||
| 388 | break; | ||
| 389 | } | ||
| 390 | } | ||
| 391 | faceIdxNeighbor++; | ||
| 392 | } | ||
| 393 | } | ||
| 394 | |||
| 395 | intersectionMapGlobal_[eIdxGlobal][fIdx] = globalIntersectionIdx; | ||
| 396 | intersectionMapGlobal_[globalIdxNeighbor][faceIdxNeighbor] = globalIntersectionIdx; | ||
| 397 | globalIntersectionIdx ++; | ||
| 398 | } | ||
| 399 | } | ||
| 400 | else | ||
| 401 | { | ||
| 402 | intersectionMapGlobal_[eIdxGlobal][fIdx] = globalIntersectionIdx; | ||
| 403 | globalIntersectionIdx ++; | ||
| 404 | } | ||
| 405 | fIdx++; | ||
| 406 | } | ||
| 407 | } | ||
| 408 | size_ = globalIntersectionIdx; | ||
| 409 | } | ||
| 410 | GridView gridView_; | ||
| 411 | ElementMapper elementMapper_; | ||
| 412 | unsigned int size_; | ||
| 413 | std::vector<std::unordered_map<int, int> > intersectionMapGlobal_; | ||
| 414 | std::vector<std::unordered_map<int, int> > intersectionMapLocal_; | ||
| 415 | std::unordered_map<int, int> standardLocalIdxMap_; | ||
| 416 | }; | ||
| 417 | |||
| 418 | } // end namespace Dumux | ||
| 419 | |||
| 420 | #endif | ||
| 421 |