3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
staggeredgeometryhelper.hh
Go to the documentation of this file.
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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_DISCRETIZATION_STAGGERED_GEOMETRY_HELPER_HH
25#define DUMUX_DISCRETIZATION_STAGGERED_GEOMETRY_HELPER_HH
26
27#include <dune/geometry/multilineargeometry.hh>
28
29#include <dumux/common/math.hh>
31#include <type_traits>
32#include <algorithm>
33#include <array>
34#include <bitset>
35
36namespace Dumux {
37namespace Detail {
38
70template<class GridView, int upwindSchemeOrder>
72{
73 using Scalar = typename GridView::ctype;
74 using GlobalPosition = typename GridView::template Codim<0>::Entity::Geometry::GlobalCoordinate;
77
78 std::bitset<upwindSchemeOrder> hasParallelNeighbor;
81 std::array<GridIndexType, upwindSchemeOrder> parallelDofs;
82 std::array<Scalar, upwindSchemeOrder> parallelCellWidths;
83 bool hasOuterLateral = false;
84 std::pair<GridIndexType, GridIndexType> lateralPair;
88};
89
94template<class GridView, int upwindSchemeOrder>
96{
97 using Scalar = typename GridView::ctype;
99
102 std::bitset<upwindSchemeOrder-1> hasForwardNeighbor;
103 std::bitset<upwindSchemeOrder-1> hasBackwardNeighbor;
104 std::array<GridIndexType, upwindSchemeOrder-1> inAxisForwardDofs;
105 std::array<GridIndexType, upwindSchemeOrder-1> inAxisBackwardDofs;
107 std::array<Scalar, upwindSchemeOrder-1> inAxisForwardDistances;
108 std::array<Scalar, upwindSchemeOrder-1> inAxisBackwardDistances;
109};
110
115template<class GridView>
116struct AxisData<GridView, 1>
117{
118 using Scalar = typename GridView::ctype;
120
124};
125
126} // namespace Detail
127
132template<class Vector>
133inline static unsigned int directionIndex(Vector&& vector)
134{
135 const auto eps = 1e-8;
136 return std::find_if(vector.begin(), vector.end(), [eps](const auto& x) { return std::abs(x) > eps; } ) - vector.begin();
137}
138
144template<class GridView, int upwindSchemeOrder>
146{
147 using Scalar = typename GridView::ctype;
148 static constexpr auto dim = GridView::dimension;
149
150 using Element = typename GridView::template Codim<0>::Entity;
151 using Intersection = typename GridView::Intersection;
152 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
153
154 using GridIndexType = typename IndexTraits<GridView>::GridIndex;
155 using SmallLocalIndexType = typename IndexTraits<GridView>::SmallLocalIndex;
156
157 //TODO include assert that checks for quad geometry
158 static constexpr auto codimIntersection = 1;
159 static constexpr auto numFacets = dim * 2;
160 static constexpr auto numPairs = 2 * (dim - 1);
161
162 static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
163
164public:
167
168 FreeFlowStaggeredGeometryHelper(const Element& element, const GridView& gridView)
169 : element_(element)
170 , gridView_(gridView)
171 { }
172
174 template<class IntersectionMapper>
175 void updateLocalFace(const IntersectionMapper&, const Intersection& intersection)
176 {
177 intersection_ = intersection;
178 fillAxisData_();
179 fillPairData_();
180 }
181
185 SmallLocalIndexType localFaceIndex() const
186 {
187 return intersection_.indexInInside();
188 }
189
194 {
195 return axisData_;
196 }
197
201 std::array<PairData, numPairs> pairData() const
202 {
203 return pairData_;
204 }
205
209 unsigned int directionIndex() const
210 {
211 return Dumux::directionIndex(intersection_.centerUnitOuterNormal());
212 }
213
217 unsigned int directionIndex(const Intersection& intersection) const
218 {
219 return Dumux::directionIndex(std::move(intersection.centerUnitOuterNormal()));
220 }
221
222private:
223
228 void fillAxisData_()
229 {
230 if constexpr (useHigherOrder)
231 fillOuterAxisData_();
232
233 const auto inIdx = intersection_.indexInInside();
234 const auto oppIdx = localOppositeIdx_(inIdx);
235
236 // Set the self Dof
237 axisData_.selfDof = gridView_.indexSet().subIndex(intersection_.inside(), inIdx, codimIntersection);
238
239 // Set the opposite Dof
240 axisData_.oppositeDof = gridView_.indexSet().subIndex(intersection_.inside(), oppIdx, codimIntersection);
241
242 // Set the Self to Opposite Distance
243 const auto self = getFacet_(inIdx, element_);
244 const auto opposite = getFacet_(oppIdx, element_);
245 axisData_.selfToOppositeDistance = (self.geometry().center() - opposite.geometry().center()).two_norm();
246
247 }
248
254 void fillOuterAxisData_()
255 {
256 // reset the axis data struct
257 axisData_ = {};
258
259 // Set the forward dofs
260 const auto thisFaceLocalIdx = intersection_.indexInInside();
261 const auto& faceCenterPos = getFacet_(thisFaceLocalIdx, element_).geometry().center();
262 if (intersection_.neighbor())
263 addForwardNeighborAxisData_(intersection_.outside(), 0, thisFaceLocalIdx, faceCenterPos);
264
265 // Set the backward dofs
266 const auto oppositeFaceLocalIdx = localOppositeIdx_(thisFaceLocalIdx);
267 const auto& oppositeFaceCenterPos = getFacet_(oppositeFaceLocalIdx, element_).geometry().center();
268 for (const auto& intersection : intersections(gridView_, element_))
269 {
270 if (intersection.indexInInside() == oppositeFaceLocalIdx && intersection.neighbor())
271 {
272 addBackwardNeighborAxisData_(intersection.outside(), 0, oppositeFaceLocalIdx, oppositeFaceCenterPos);
273 break;
274 }
275 }
276 }
277
281 void addForwardNeighborAxisData_(const Element& neighbor,
282 const std::size_t distance,
283 const unsigned int localIntersectionIndex,
284 const GlobalPosition& faceCenterPos)
285 {
286 // fill the forward axis data for this element in the stencil
287 axisData_.hasForwardNeighbor.set(distance, true);
288 axisData_.inAxisForwardDofs[distance] = gridView_.indexSet().subIndex(neighbor, localIntersectionIndex, codimIntersection);
289 const auto& forwardFacePos = getFacet_(localIntersectionIndex, neighbor).geometry().center();
290 axisData_.inAxisForwardDistances[distance] = (forwardFacePos - faceCenterPos).two_norm();
291
292 // recursion end if we added all axis data for the given stencil size (numParallelFaces)
293 const auto numForwardBackwardAxisDofs = axisData_.inAxisForwardDofs.size();
294 if (distance >= numForwardBackwardAxisDofs-1)
295 return;
296
297 for (const auto& intersection : intersections(gridView_, neighbor))
298 {
299 if (intersection.indexInInside() == localIntersectionIndex && intersection.neighbor())
300 {
301 addForwardNeighborAxisData_(intersection.outside(), distance+1, localIntersectionIndex, forwardFacePos);
302 break;
303 }
304 }
305 }
306
310 void addBackwardNeighborAxisData_(const Element& neighbor,
311 const std::size_t distance,
312 const unsigned int localIntersectionIndex,
313 const GlobalPosition& faceCenterPos)
314 {
315 // fill the forward axis data for this element in the stencil
316 axisData_.hasBackwardNeighbor.set(distance, true);
317 axisData_.inAxisBackwardDofs[distance] = gridView_.indexSet().subIndex(neighbor, localIntersectionIndex, codimIntersection);
318 const auto& backwardFacePos = getFacet_(localIntersectionIndex, neighbor).geometry().center();
319 axisData_.inAxisBackwardDistances[distance] = (backwardFacePos - faceCenterPos).two_norm();
320
321 // recursion end if we added all axis data for the given stencil size (numParallelFaces)
322 const auto numForwardBackwardAxisDofs = axisData_.inAxisForwardDofs.size();
323 if (distance >= numForwardBackwardAxisDofs-1)
324 return;
325
326 for (const auto& intersection : intersections(gridView_, neighbor))
327 {
328 if (intersection.indexInInside() == localIntersectionIndex && intersection.neighbor())
329 {
330 addBackwardNeighborAxisData_(intersection.outside(), distance+1, localIntersectionIndex, backwardFacePos);
331 break;
332 }
333 }
334 }
335
340 void fillPairData_()
341 {
342 // reset the pair data structs
343 pairData_ = {};
344
345 // set basic global positions
346 const auto& selfElementCenter = element_.geometry().center();
347 const auto& selfFacetCenter = intersection_.geometry().center();
348
349 // get the inner lateral Dof Index
350 SmallLocalIndexType numPairInnerLateralIdx = 0;
351 for (const auto& innerElementIntersection : intersections(gridView_, element_))
352 {
353 if (facetIsNormal_(innerElementIntersection.indexInInside(), intersection_.indexInInside()))
354 {
355 const auto innerElementIntersectionIdx = innerElementIntersection.indexInInside();
356 setLateralPairFirstInfo_(innerElementIntersectionIdx, element_, numPairInnerLateralIdx);
357
358 const auto distance = innerElementIntersection.geometry().center() - selfElementCenter;
359 pairData_[numPairInnerLateralIdx].lateralStaggeredFaceCenter = selfFacetCenter+distance;
360
361 numPairInnerLateralIdx++;
362 }
363 }
364
365 // get the outer lateral Dof Index
366 SmallLocalIndexType numPairOuterLateralIdx = 0;
367 if (intersection_.neighbor())
368 {
369 // the direct neighbor element and the respective intersection index
370 const auto& directNeighborElement = intersection_.outside();
371
372 for (const auto& directNeighborElementIntersection : intersections(gridView_, directNeighborElement))
373 {
374 // skip the directly neighboring face itself and its opposing one
375 if (facetIsNormal_(directNeighborElementIntersection.indexInInside(), intersection_.indexInOutside()))
376 {
377 const auto directNeighborElemIsIdx = directNeighborElementIntersection.indexInInside();
378 setLateralPairSecondInfo_(directNeighborElemIsIdx, directNeighborElement, numPairOuterLateralIdx);
379 numPairOuterLateralIdx++;
380 }
381 }
382 }
383 else // intersection is on boundary
384 {
385 for (const auto& intersection : intersections(gridView_, element_))
386 {
387 if (facetIsNormal_(intersection.indexInInside(), intersection_.indexInInside()))
388 {
389 assert(!pairData_[numPairOuterLateralIdx].hasOuterLateral);
390
391 const auto normalDistanceoffset = selfFacetCenter - selfElementCenter;
392 pairData_[numPairOuterLateralIdx].lateralDistance = normalDistanceoffset.two_norm();
393 numPairOuterLateralIdx++;
394 }
395 }
396 }
397
398 // fill the pair data with the parallel dofs and distances
399 const auto parallelLocalIdx = intersection_.indexInInside();
400 SmallLocalIndexType numPairParallelIdx = 0;
401 for (const auto& intersection : intersections(gridView_, element_))
402 {
403 if (facetIsNormal_(intersection.indexInInside(), parallelLocalIdx))
404 {
405 if (intersection.neighbor())
406 {
407 // recursively insert parallel neighbor faces into pair data
408 const auto parallelAxisIdx = directionIndex(intersection);
409 const auto localLateralIntersectionIndex = intersection.indexInInside();
410
411 /*
412 * ------------
413 * | |
414 * | |
415 * | |
416 * iiiiiiiiiii*bbbbbbbbbbb
417 * | o zzzzzzzz |
418 * | o zzzzzzzz |
419 * | o zzzzzzzz |
420 * -----------------------
421 *
422 * i:intersection,o:intersection_, b: outerIntersection, z: intersection_.outside()
423 */
424 if (intersection_.neighbor())
425 for (const auto& outerIntersection : intersections(gridView_, intersection_.outside()))
426 if (intersection.indexInInside() == outerIntersection.indexInInside())
427 if (!outerIntersection.neighbor())
428 pairData_[numPairParallelIdx].hasHalfParallelNeighbor = true;
429
430 /* ------------
431 * | o
432 * | o
433 * | o
434 * iiiiiiiiiii------------
435 * | zzzzzzzz b |
436 * | zzzzzzzz b |
437 * | zzzzzzzz b |
438 * -----------------------
439 *
440 * i:intersection,o:intersection_, b: outerIntersection, z: intersection.outside()
441 */
442 if (!intersection_.neighbor())
443 for (const auto& outerIntersection : intersections(gridView_, intersection.outside()))
444 if (intersection_.indexInInside() == outerIntersection.indexInInside())
445 if (outerIntersection.neighbor())
446 pairData_[numPairParallelIdx].hasCornerParallelNeighbor = true;
447
448
449 addParallelNeighborPairData_(intersection.outside(), 0, localLateralIntersectionIndex, parallelLocalIdx, parallelAxisIdx, numPairParallelIdx);
450 }
451
452 ++numPairParallelIdx;
453 }
454 }
455 }
456
457 // zero distance means direct neighbor to element_
458 void addParallelNeighborPairData_(const Element& neighbor,
459 const std::size_t distance,
460 const unsigned int localLateralIntersectionIndex,
461 const unsigned int parallelLocalIdx,
462 const unsigned int parallelAxisIdx,
463 const SmallLocalIndexType dataIdx)
464 {
465 // fill the pair data for this element in the stencil
466 pairData_[dataIdx].hasParallelNeighbor.set(distance, true);
467 pairData_[dataIdx].parallelDofs[distance] = gridView_.indexSet().subIndex(neighbor, parallelLocalIdx, codimIntersection);
468 pairData_[dataIdx].parallelCellWidths[distance] = setParallelPairCellWidths_(neighbor, parallelAxisIdx);
469
470 // recursion end if we added all parallel data for the given stencil size (numParallelFaces)
471 const auto numParallelFaces = pairData_[0].parallelCellWidths.size();
472 if (distance >= numParallelFaces-1)
473 return;
474
475 // continue with neighbor's neighbor element in the same direction, if we find one
476 for (const auto& lateralIntersection : intersections(gridView_, neighbor))
477 {
478 // we only expect to find one intersection matching this condition
479 if (lateralIntersection.indexInInside() == localLateralIntersectionIndex)
480 {
481 // if there is no neighbor recursion ends here
482 if (lateralIntersection.neighbor())
483 addParallelNeighborPairData_(lateralIntersection.outside(), distance+1,
484 localLateralIntersectionIndex, parallelLocalIdx, parallelAxisIdx, dataIdx);
485 break;
486 }
487 }
488 }
489
495 int localOppositeIdx_(const int idx) const
496 {
497 return (idx % 2) ? (idx - 1) : (idx + 1);
498 }
499
506 bool facetIsNormal_(const int selfIdx, const int otherIdx) const
507 {
508 return !(selfIdx == otherIdx || localOppositeIdx_(selfIdx) == otherIdx);
509 }
510
511 auto getFacet_(const int localFacetIdx, const Element& element) const
512 {
513 return element.template subEntity <1> (localFacetIdx);
514 }
515
517 void setLateralPairFirstInfo_(const int isIdx, const Element& element, const int numPairsIdx)
518 {
519 // store the inner lateral dofIdx
520 pairData_[numPairsIdx].lateralPair.first = gridView_.indexSet().subIndex(element, isIdx, codimIntersection);
521
522 // store the local lateral facet index
523 pairData_[numPairsIdx].localLateralFaceIdx = isIdx;
524 }
525
527 void setLateralPairSecondInfo_(const int isIdx, const Element& element, const int numPairsIdx)
528 {
529 // store the dofIdx
530 pairData_[numPairsIdx].hasOuterLateral = true;
531 pairData_[numPairsIdx].lateralPair.second = gridView_.indexSet().subIndex(element, isIdx, codimIntersection);
532
533 // set basic global positions
534 const auto& selfFacetCenter = intersection_.geometry().center();
535 const auto& selfElementCenter = element_.geometry().center();
536
537 const auto& neighborElement = intersection_.outside();
538 const auto& neighborElementCenter = neighborElement.geometry().center();
539 const auto& neighborFacetCenter = getFacet_(intersection_.indexInOutside(), neighborElement).geometry().center();
540
541 const Scalar insideLateralDistance = (selfFacetCenter - selfElementCenter).two_norm();
542 const Scalar outsideLateralDistance = (neighborFacetCenter - neighborElementCenter).two_norm();
543
544 pairData_[numPairsIdx].lateralDistance = insideLateralDistance + outsideLateralDistance;
545 }
546
548 Scalar setParallelPairCellWidths_(const Element& element, const int parallelAxisIdx) const
549 {
550 std::vector<GlobalPosition> faces;
551 faces.reserve(numFacets);
552 for (const auto& intersection : intersections(gridView_, element))
553 faces.push_back(intersection.geometry().center());
554
555 switch (parallelAxisIdx)
556 {
557 case 0:
558 return (faces[1] - faces[0]).two_norm();
559 case 1:
560 return (faces[3] - faces[2]).two_norm();
561 case 2:
562 {
563 assert(dim == 3);
564 return (faces[5] - faces[4]).two_norm();
565 }
566 default:
567 DUNE_THROW(Dune::InvalidStateException, "Something went terribly wrong");
568 }
569 }
570
571 Intersection intersection_;
572 const Element& element_;
573 const GridView gridView_;
574 AxisData axisData_;
575 std::array<PairData, numPairs> pairData_;
576};
577
578} // end namespace Dumux
579
580#endif
Defines the index types used for grid and local indices.
Define some often used mathematical functions.
ctype distance(const Dune::FieldVector< ctype, dimWorld > &a, const Dune::FieldVector< ctype, dimWorld > &b)
Compute the shortest distance between two points.
Definition: geometry/distance.hh:138
static unsigned int directionIndex(Vector &&vector)
Returns the dirction index of the facet (0 = x, 1 = y, 2 = z)
Definition: staggeredgeometryhelper.hh:133
Definition: adapt.hh:29
typename GridView::IndexSet::IndexType GridIndex
Definition: indextraits.hh:39
std::uint_least8_t SmallLocalIndex
Definition: indextraits.hh:41
defines an intersection mapper for mapping of global DOFs assigned to faces which also works for adap...
Definition: intersectionmapper.hh:185
Parallel Data stored per sub face.
Definition: staggeredgeometryhelper.hh:72
std::pair< GridIndexType, GridIndexType > lateralPair
Definition: staggeredgeometryhelper.hh:84
typename IndexTraits< GridView >::GridIndex GridIndexType
Definition: staggeredgeometryhelper.hh:75
bool hasOuterLateral
Definition: staggeredgeometryhelper.hh:83
typename GridView::ctype Scalar
Definition: staggeredgeometryhelper.hh:73
GlobalPosition lateralStaggeredFaceCenter
Definition: staggeredgeometryhelper.hh:87
std::bitset< upwindSchemeOrder > hasParallelNeighbor
Definition: staggeredgeometryhelper.hh:78
bool hasHalfParallelNeighbor
Definition: staggeredgeometryhelper.hh:79
bool hasCornerParallelNeighbor
Definition: staggeredgeometryhelper.hh:80
SmallLocalIndexType localLateralFaceIdx
Definition: staggeredgeometryhelper.hh:85
typename IndexTraits< GridView >::SmallLocalIndex SmallLocalIndexType
Definition: staggeredgeometryhelper.hh:76
std::array< Scalar, upwindSchemeOrder > parallelCellWidths
Definition: staggeredgeometryhelper.hh:82
typename GridView::template Codim< 0 >::Entity::Geometry::GlobalCoordinate GlobalPosition
Definition: staggeredgeometryhelper.hh:74
std::array< GridIndexType, upwindSchemeOrder > parallelDofs
Definition: staggeredgeometryhelper.hh:81
Scalar lateralDistance
Definition: staggeredgeometryhelper.hh:86
In Axis Data stored per sub face.
Definition: staggeredgeometryhelper.hh:96
Scalar selfToOppositeDistance
Definition: staggeredgeometryhelper.hh:106
GridIndexType selfDof
Definition: staggeredgeometryhelper.hh:100
std::array< Scalar, upwindSchemeOrder-1 > inAxisForwardDistances
Definition: staggeredgeometryhelper.hh:107
std::bitset< upwindSchemeOrder-1 > hasBackwardNeighbor
Definition: staggeredgeometryhelper.hh:103
std::array< GridIndexType, upwindSchemeOrder-1 > inAxisForwardDofs
Definition: staggeredgeometryhelper.hh:104
std::array< Scalar, upwindSchemeOrder-1 > inAxisBackwardDistances
Definition: staggeredgeometryhelper.hh:108
std::bitset< upwindSchemeOrder-1 > hasForwardNeighbor
Definition: staggeredgeometryhelper.hh:102
typename IndexTraits< GridView >::GridIndex GridIndexType
Definition: staggeredgeometryhelper.hh:98
typename GridView::ctype Scalar
Definition: staggeredgeometryhelper.hh:97
GridIndexType oppositeDof
Definition: staggeredgeometryhelper.hh:101
std::array< GridIndexType, upwindSchemeOrder-1 > inAxisBackwardDofs
Definition: staggeredgeometryhelper.hh:105
GridIndexType oppositeDof
Definition: staggeredgeometryhelper.hh:122
typename IndexTraits< GridView >::GridIndex GridIndexType
Definition: staggeredgeometryhelper.hh:119
Scalar selfToOppositeDistance
Definition: staggeredgeometryhelper.hh:123
typename GridView::ctype Scalar
Definition: staggeredgeometryhelper.hh:118
GridIndexType selfDof
Definition: staggeredgeometryhelper.hh:121
Helper class constructing the dual grid finite volume geometries for the free flow staggered discreti...
Definition: staggeredgeometryhelper.hh:146
Detail::AxisData< GridView, upwindSchemeOrder > AxisData
Definition: staggeredgeometryhelper.hh:166
SmallLocalIndexType localFaceIndex() const
Returns the local index of the face (i.e. the intersection)
Definition: staggeredgeometryhelper.hh:185
FreeFlowStaggeredGeometryHelper(const Element &element, const GridView &gridView)
Definition: staggeredgeometryhelper.hh:168
AxisData axisData() const
Returns a copy of the axis data.
Definition: staggeredgeometryhelper.hh:193
unsigned int directionIndex() const
Returns the direction index of the primary facet (0 = x, 1 = y, 2 = z)
Definition: staggeredgeometryhelper.hh:209
void updateLocalFace(const IntersectionMapper &, const Intersection &intersection)
update the local face
Definition: staggeredgeometryhelper.hh:175
std::array< PairData, numPairs > pairData() const
Returns a copy of the pair data.
Definition: staggeredgeometryhelper.hh:201
unsigned int directionIndex(const Intersection &intersection) const
Returns the direction index of the facet passed as an argument (0 = x, 1 = y, 2 = z)
Definition: staggeredgeometryhelper.hh:217