3.4
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
helper.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_CC_MPFA_HELPER_HH
25#define DUMUX_DISCRETIZATION_CC_MPFA_HELPER_HH
26
27#include <dune/common/exceptions.hh>
28#include <dune/common/fvector.hh>
29#include <dune/common/fmatrix.hh>
30#include <dune/common/parallel/mpihelper.hh>
31#include <dune/geometry/type.hh>
32
33#include <dumux/common/math.hh>
34
35namespace Dumux {
36
41template<class GridGeometry, int dim, int dimWorld>
43
48template<class GridGeometry>
49class MpfaDimensionHelper<GridGeometry, /*dim*/2, /*dimWorld*/2>
50{
51 using GridView = typename GridGeometry::GridView;
52 using CoordScalar = typename GridView::ctype;
53 using Element = typename GridView::template Codim<0>::Entity;
54
55 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
56
57 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
58 using ScvfCornerVector = typename SubControlVolumeFace::Traits::CornerStorage;
59
60 // Container to store the positions of intersections required for scvf
61 // corner computation. In 2d, these are the center plus the two corners
62 using ScvfPositionsOnIntersection = std::array<GlobalPosition, 3>;
63
64public:
69 template< class ScvBasis >
70 static ScvBasis calculateInnerNormals(const ScvBasis& scvBasis)
71 {
72 static const Dune::FieldMatrix<CoordScalar, 2, 2> R = {{0.0, 1.0}, {-1.0, 0.0}};
73
74 ScvBasis innerNormals;
75 R.mv(scvBasis[1], innerNormals[0]);
76 R.mv(scvBasis[0], innerNormals[1]);
77
78 // adjust sign depending on basis being a RHS
79 if (!isRightHandSystem(scvBasis))
80 innerNormals[0] *= -1.0;
81 else
82 innerNormals[1] *= -1.0;
83
84 return innerNormals;
85 }
86
92 template< class ScvBasis >
93 static CoordScalar calculateDetX(const ScvBasis& scvBasis)
94 {
95 using std::abs;
96 return abs(crossProduct<CoordScalar>(scvBasis[0], scvBasis[1]));
97 }
98
104 static std::size_t getGlobalNumScvf(const GridView& gridView)
105 {
106 assert(gridView.size(Dune::GeometryTypes::triangle)
107 + gridView.size(Dune::GeometryTypes::quadrilateral) == gridView.size(0));
108
109 return gridView.size(Dune::GeometryTypes::triangle)
110 * getNumLocalScvfs(Dune::GeometryTypes::triangle)
111 + gridView.size(Dune::GeometryTypes::quadrilateral)
112 * getNumLocalScvfs(Dune::GeometryTypes::quadrilateral);
113 }
114
119 template< class ScvBasis >
120 static bool isRightHandSystem(const ScvBasis& scvBasis)
121 { return !std::signbit(crossProduct<CoordScalar>(scvBasis[0], scvBasis[1])); }
122
133 template<class ElementGeometry, class ReferenceElement>
134 static ScvfPositionsOnIntersection computeScvfCornersOnIntersection(const ElementGeometry& eg,
135 const ReferenceElement& refElement,
136 unsigned int indexInElement,
137 unsigned int numCorners)
138 {
139 ScvfPositionsOnIntersection p;
140
141 // compute facet center and corners
142 p[0] = 0.0;
143 for (unsigned int c = 0; c < numCorners; ++c)
144 {
145 // codim = 1, dim = 2
146 p[c+1] = eg.global(refElement.position(refElement.subEntity(indexInElement, 1, c, 2), 2));
147 p[0] += p[c+1];
148 }
149 p[0] /= numCorners;
150
151 return p;
152 }
153
162 static ScvfCornerVector getScvfCorners(const ScvfPositionsOnIntersection& p,
163 unsigned int numIntersectionCorners,
164 unsigned int cornerIdx)
165 {
166 // make sure the given input is admissible
167 assert(cornerIdx < 2 && "provided index exceeds the number of corners of facets in 2d");
168
169 // create & return the scvf corner vector
170 return cornerIdx == 0 ? ScvfCornerVector({p[0], p[1]})
171 : ScvfCornerVector({p[0], p[2]});
172 }
173
178 static CoordScalar computeScvfArea(const ScvfCornerVector& scvfCorners)
179 { return (scvfCorners[1]-scvfCorners[0]).two_norm(); }
180
185 static std::size_t getNumLocalScvfs(const Dune::GeometryType& gt)
186 {
187 if (gt == Dune::GeometryTypes::triangle)
188 return 6;
189 else if (gt == Dune::GeometryTypes::quadrilateral)
190 return 8;
191 else
192 DUNE_THROW(Dune::NotImplemented, "Mpfa for 2d geometry type " << gt);
193 }
194};
195
201template<class GridGeometry>
202class MpfaDimensionHelper<GridGeometry, /*dim*/2, /*dimWorld*/3>
203: public MpfaDimensionHelper<GridGeometry, 2, 2>
204{
205 using GridView = typename GridGeometry::GridView;
206 using CoordScalar = typename GridView::ctype;
207public:
208
213 template< class ScvBasis >
214 static ScvBasis calculateInnerNormals(const ScvBasis& scvBasis)
215 {
216 // compute vector normal to the basis plane
217 const auto normal = [&scvBasis] ()
218 {
219 auto n = crossProduct<CoordScalar>(scvBasis[0], scvBasis[1]);
220 n /= n.two_norm();
221 return n;
222 } ();
223
224 // compute inner normals using the normal vector
225 ScvBasis innerNormals;
226 innerNormals[0] = crossProduct<CoordScalar>(scvBasis[1], normal);
227 innerNormals[1] = crossProduct<CoordScalar>(normal, scvBasis[0]);
228
229 return innerNormals;
230 }
231
240 template< class ScvBasis >
241 static CoordScalar calculateDetX(const ScvBasis& scvBasis)
242 {
243 using std::abs;
244 return abs(crossProduct<CoordScalar>(scvBasis[0], scvBasis[1]).two_norm());
245 }
246
254 template< class ScvBasis >
255 static constexpr bool isRightHandSystem(const ScvBasis& scvBasis) { return true; }
256};
262template<class GridGeometry>
263class MpfaDimensionHelper<GridGeometry, /*dim*/3, /*dimWorld*/3>
264{
265 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
266 using ScvfCornerVector = typename SubControlVolumeFace::Traits::CornerStorage;
267
268 // Be picky about the dimensions
269 using GridView = typename GridGeometry::GridView;
270 using Element = typename GridView::template Codim<0>::Entity;
271 using CoordScalar = typename GridView::ctype;
272 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
273
274 // container to store the positions of intersections required for
275 // scvf corner computation. Maximum number of points needed is 9
276 // for the supported geometry types (quadrilateral facet)
277 using ScvfPositionsOnIntersection = std::array<GlobalPosition, 9>;
278
279public:
280
286 template< class ScvBasis >
287 static ScvBasis calculateInnerNormals(const ScvBasis& scvBasis)
288 {
289 ScvBasis innerNormals;
290
291 innerNormals[0] = crossProduct<CoordScalar>(scvBasis[1], scvBasis[2]);
292 innerNormals[1] = crossProduct<CoordScalar>(scvBasis[2], scvBasis[0]);
293 innerNormals[2] = crossProduct<CoordScalar>(scvBasis[0], scvBasis[1]);
294
295 if (!isRightHandSystem(scvBasis))
296 std::for_each(innerNormals.begin(), innerNormals.end(), [] (auto& n) { n *= -1.0; });
297
298 return innerNormals;
299 }
300
307 template< class ScvBasis >
308 static CoordScalar calculateDetX(const ScvBasis& scvBasis)
309 {
310 using std::abs;
311 return abs(tripleProduct<CoordScalar>(scvBasis[0], scvBasis[1], scvBasis[2]));
312 }
313
320 static std::size_t getGlobalNumScvf(const GridView& gridView)
321 {
322 assert(gridView.size(Dune::GeometryTypes::tetrahedron)
323 + gridView.size(Dune::GeometryTypes::pyramid)
324 + gridView.size(Dune::GeometryTypes::prism)
325 + gridView.size(Dune::GeometryTypes::hexahedron) == gridView.size(0));
326
327 return gridView.size(Dune::GeometryTypes::tetrahedron)
328 * getNumLocalScvfs(Dune::GeometryTypes::tetrahedron)
329 + gridView.size(Dune::GeometryTypes::pyramid)
330 * getNumLocalScvfs(Dune::GeometryTypes::pyramid)
331 + gridView.size(Dune::GeometryTypes::prism)
332 * getNumLocalScvfs(Dune::GeometryTypes::prism)
333 + gridView.size(Dune::GeometryTypes::hexahedron)
334 * getNumLocalScvfs(Dune::GeometryTypes::hexahedron);
335 }
336
341 template< class ScvBasis >
342 static bool isRightHandSystem(const ScvBasis& scvBasis)
343 { return !std::signbit(tripleProduct<CoordScalar>(scvBasis[0], scvBasis[1], scvBasis[2])); }
344
355 template<class ElementGeometry, class ReferenceElement>
356 static ScvfPositionsOnIntersection computeScvfCornersOnIntersection(const ElementGeometry& eg,
357 const ReferenceElement& refElement,
358 unsigned int indexInElement,
359 unsigned int numCorners)
360 {
361 // The size of ScvfPositionsOnIntersection doesn't allow for faces with more than four corners!
362 ScvfPositionsOnIntersection p;
363 if (numCorners > 4)
364 DUNE_THROW(Dune::InvalidStateException, "Mpfa implementation cannot handle faces with more than 4 corners");
365
366 // compute facet center and corners
367 p[0] = 0.0;
368 for (unsigned int c = 0; c < numCorners; ++c)
369 {
370 // codim = 1, dim = 3
371 p[c+1] = eg.global(refElement.position(refElement.subEntity(indexInElement, 1, c, 3), 3));
372 p[0] += p[c+1];
373 }
374 p[0] /= numCorners;
375
376 // proceed according to number of corners
377 switch (numCorners)
378 {
379 case 3: // triangle
380 {
381 // add edge midpoints, there are 3 for triangles
382 p[numCorners+1] = p[2] + p[1];
383 p[numCorners+1] /= 2;
384 p[numCorners+2] = p[3] + p[1];
385 p[numCorners+2] /= 2;
386 p[numCorners+3] = p[3] + p[2];
387 p[numCorners+3] /= 2;
388 return p;
389 }
390 case 4: // quadrilateral
391 {
392 // add edge midpoints, there are 4 for quadrilaterals
393 p[numCorners+1] = p[3] + p[1];
394 p[numCorners+1] /= 2;
395 p[numCorners+2] = p[4] + p[2];
396 p[numCorners+2] /= 2;
397 p[numCorners+3] = p[2] + p[1];
398 p[numCorners+3] /= 2;
399 p[numCorners+4] = p[4] + p[3];
400 p[numCorners+4] /= 2;
401 return p;
402 }
403 default:
404 DUNE_THROW(Dune::NotImplemented,
405 "Mpfa scvf corners for dim = 3, dimWorld = 3, corners = " << numCorners);
406 }
407 }
408
417 static ScvfCornerVector getScvfCorners(const ScvfPositionsOnIntersection& p,
418 unsigned int numIntersectionCorners,
419 unsigned int cornerIdx)
420 {
421 // proceed according to number of corners
422 // we assume the ordering according to the above method computeScvfCornersOnIntersection()
423 switch (numIntersectionCorners)
424 {
425 case 3: // triangle
426 {
427 // Only build the maps the first time we encounter a triangle
428 static const std::uint8_t vo = 1; // vertex offset in point vector p
429 static const std::uint8_t eo = 4; // edge offset in point vector p
430 static const std::uint8_t map[3][4] =
431 {
432 {0, eo+1, eo+0, vo+0},
433 {0, eo+0, eo+2, vo+1},
434 {0, eo+2, eo+1, vo+2}
435 };
436
437 return ScvfCornerVector( {p[map[cornerIdx][0]],
438 p[map[cornerIdx][1]],
439 p[map[cornerIdx][2]],
440 p[map[cornerIdx][3]]} );
441 }
442 case 4: // quadrilateral
443 {
444 // Only build the maps the first time we encounter a quadrilateral
445 static const std::uint8_t vo = 1; // vertex offset in point vector p
446 static const std::uint8_t eo = 5; // face offset in point vector p
447 static const std::uint8_t map[4][4] =
448 {
449 {0, eo+0, eo+2, vo+0},
450 {0, eo+2, eo+1, vo+1},
451 {0, eo+3, eo+0, vo+2},
452 {0, eo+1, eo+3, vo+3}
453 };
454
455 return ScvfCornerVector( {p[map[cornerIdx][0]],
456 p[map[cornerIdx][1]],
457 p[map[cornerIdx][2]],
458 p[map[cornerIdx][3]]} );
459 }
460 default:
461 DUNE_THROW(Dune::NotImplemented,
462 "Mpfa scvf corners for dim = 3, dimWorld = 3, corners = " << numIntersectionCorners);
463 }
464 }
465
470 static CoordScalar computeScvfArea(const ScvfCornerVector& scvfCorners)
471 {
472 // after Wolfram alpha quadrilateral area
473 return 0.5*Dumux::crossProduct(scvfCorners[3]-scvfCorners[0], scvfCorners[2]-scvfCorners[1]).two_norm();
474 }
475
481 static std::size_t getNumLocalScvfs(const Dune::GeometryType& gt)
482 {
483 if (gt == Dune::GeometryTypes::tetrahedron)
484 return 12;
485 else if (gt == Dune::GeometryTypes::pyramid)
486 return 16;
487 else if (gt == Dune::GeometryTypes::prism)
488 return 18;
489 else if (gt == Dune::GeometryTypes::hexahedron)
490 return 24;
491 else
492 DUNE_THROW(Dune::NotImplemented, "Mpfa for 3d geometry type " << gt);
493 }
494};
495
502template<class GridGeometry>
503class CCMpfaHelper : public MpfaDimensionHelper<GridGeometry,
504 GridGeometry::GridView::dimension,
505 GridGeometry::GridView::dimensionworld>
506{
507 using PrimaryInteractionVolume = typename GridGeometry::GridIVIndexSets::PrimaryInteractionVolume;
508 using SecondaryInteractionVolume = typename GridGeometry::GridIVIndexSets::SecondaryInteractionVolume;
509
510 using VertexMapper = typename GridGeometry::VertexMapper;
511 using FVElementGeometry = typename GridGeometry::LocalView;
512 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
513 using ScvfCornerVector = typename SubControlVolumeFace::Traits::CornerStorage;
514
515 using GridView = typename GridGeometry::GridView;
516 static constexpr int dim = GridView::dimension;
517
518 using Element = typename GridView::template Codim<0>::Entity;
519 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
520 using CoordScalar = typename GridView::ctype;
521
522public:
529 template< class Scalar >
530 static GlobalPosition getScvfIntegrationPoint(const ScvfCornerVector& scvfCorners, Scalar q)
531 {
532 // ordering -> first corner: facet center, last corner: vertex
533 if (q == 0.0)
534 return scvfCorners[0];
535
536 auto ip = scvfCorners.back() - scvfCorners.front();
537 ip *= q;
538 ip += scvfCorners[0];
539 return ip;
540 }
541
548 static std::vector<bool> findGhostVertices(const GridView& gridView, const VertexMapper& vertexMapper)
549 {
550 std::vector<bool> ghostVertices(gridView.size(dim), false);
551
552 // if not run in parallel, skip the rest
553 if (Dune::MPIHelper::getCollectiveCommunication().size() == 1)
554 return ghostVertices;
555
556 // mpfa methods cannot yet handle ghost cells
557 if (gridView.ghostSize(0) > 0)
558 DUNE_THROW(Dune::InvalidStateException, "Mpfa methods in parallel do not work with ghost cells. Use overlap cells instead.");
559
560 // mpfa methods have to have overlapping cells
561 if (gridView.overlapSize(0) == 0)
562 DUNE_THROW(Dune::InvalidStateException, "Grid no overlaping cells. This is required by mpfa methods in parallel.");
563
564 for (const auto& element : elements(gridView))
565 {
566 for (const auto& is : intersections(gridView, element))
567 {
568 if (!is.neighbor() && !is.boundary())
569 {
570 const auto refElement = referenceElement(element);
571
572 for (int isVertex = 0; isVertex < is.geometry().corners(); ++isVertex)
573 {
574 const auto vIdxLocal = refElement.subEntity(is.indexInInside(), 1, isVertex, dim);
575 const auto vIdxGlobal = vertexMapper.subIndex(element, vIdxLocal, dim);
576 ghostVertices[vIdxGlobal] = true;
577 }
578 }
579 }
580 }
581
582 return ghostVertices;
583 }
584
587 static constexpr bool considerSecondaryIVs()
588 { return !std::is_same<PrimaryInteractionVolume, SecondaryInteractionVolume>::value; }
589
591 template<typename V, typename T>
592 static bool vectorContainsValue(const V& vector, const T value)
593 { return std::find(vector.begin(), vector.end(), value) != vector.end(); }
594};
595
596} // end namespace Dumux
597
598#endif
Define some often used mathematical functions.
Vector normal(const Vector &v)
Create a vector normal to the given one (v is expected to be non-zero)
Definition: normal.hh:36
Dune::FieldVector< Scalar, 3 > crossProduct(const Dune::FieldVector< Scalar, 3 > &vec1, const Dune::FieldVector< Scalar, 3 > &vec2)
Cross product of two vectors in three-dimensional Euclidean space.
Definition: math.hh:654
Definition: adapt.hh:29
std::size_t numCorners(Shape shape)
Returns the number of corners of a given geometry.
Definition: throatproperties.hh:215
Dimension-specific helper class to get data required for mpfa scheme.
Definition: helper.hh:42
static std::size_t getNumLocalScvfs(const Dune::GeometryType &gt)
Calculates the number of scvfs in a given element geometry type.
Definition: helper.hh:185
static ScvfPositionsOnIntersection computeScvfCornersOnIntersection(const ElementGeometry &eg, const ReferenceElement &refElement, unsigned int indexInElement, unsigned int numCorners)
Returns a vector containing the positions on a given intersection that are relevant for scvf corner c...
Definition: helper.hh:134
static ScvBasis calculateInnerNormals(const ScvBasis &scvBasis)
Calculates the inner normal vectors to a given scv basis.
Definition: helper.hh:70
static CoordScalar computeScvfArea(const ScvfCornerVector &scvfCorners)
Calculates the area of an scvf.
Definition: helper.hh:178
static ScvfCornerVector getScvfCorners(const ScvfPositionsOnIntersection &p, unsigned int numIntersectionCorners, unsigned int cornerIdx)
Returns the corners of the sub control volume face constructed in a corner (vertex) of an intersectio...
Definition: helper.hh:162
static CoordScalar calculateDetX(const ScvBasis &scvBasis)
Calculates the determinant of an scv basis. This is equal to the cross product for dim = dimWorld = 2...
Definition: helper.hh:93
static std::size_t getGlobalNumScvf(const GridView &gridView)
Returns the global number of scvfs in the grid. Assumes the grid to be made up of only basic geometry...
Definition: helper.hh:104
static bool isRightHandSystem(const ScvBasis &scvBasis)
Checks whether or not a given scv basis forms a right hand system.
Definition: helper.hh:120
static constexpr bool isRightHandSystem(const ScvBasis &scvBasis)
Checks whether or not a given scv basis forms a right hand system. Note that for dim = 2 < dimWorld =...
Definition: helper.hh:255
static ScvBasis calculateInnerNormals(const ScvBasis &scvBasis)
Calculates the inner normal vectors to a given scv basis.
Definition: helper.hh:214
static CoordScalar calculateDetX(const ScvBasis &scvBasis)
Calculates the determinant of an scv basis. For dim = 2 < dimWorld = 3 this is actually not the deter...
Definition: helper.hh:241
static CoordScalar calculateDetX(const ScvBasis &scvBasis)
Calculates the determinant of an scv basis. This is equal to the cross product for dim = dimWorld = 2...
Definition: helper.hh:308
static ScvfPositionsOnIntersection computeScvfCornersOnIntersection(const ElementGeometry &eg, const ReferenceElement &refElement, unsigned int indexInElement, unsigned int numCorners)
Returns a vector containing the positions on a given intersection that are relevant for scvf corner c...
Definition: helper.hh:356
static ScvBasis calculateInnerNormals(const ScvBasis &scvBasis)
Calculates the inner normal vectors to a given scv basis.
Definition: helper.hh:287
static bool isRightHandSystem(const ScvBasis &scvBasis)
Checks whether or not a given scv basis forms a right hand system.
Definition: helper.hh:342
static std::size_t getGlobalNumScvf(const GridView &gridView)
Returns the global number of scvfs in the grid. Assumes the grid to be made up of only basic geometry...
Definition: helper.hh:320
static ScvfCornerVector getScvfCorners(const ScvfPositionsOnIntersection &p, unsigned int numIntersectionCorners, unsigned int cornerIdx)
Returns the corners of the sub control volume face constructed in a corner (vertex) of an intersectio...
Definition: helper.hh:417
static std::size_t getNumLocalScvfs(const Dune::GeometryType &gt)
Calculates the number of scvfs in a given element geometry type.
Definition: helper.hh:481
static CoordScalar computeScvfArea(const ScvfCornerVector &scvfCorners)
Calculates the area of an scvf.
Definition: helper.hh:470
Helper class to get the required information on an interaction volume.
Definition: helper.hh:506
static constexpr bool considerSecondaryIVs()
Definition: helper.hh:587
static bool vectorContainsValue(const V &vector, const T value)
returns whether or not a value exists in a vector
Definition: helper.hh:592
static std::vector< bool > findGhostVertices(const GridView &gridView, const VertexMapper &vertexMapper)
Returns a vector which maps true to each vertex on processor boundaries and false otherwise.
Definition: helper.hh:548
static GlobalPosition getScvfIntegrationPoint(const ScvfCornerVector &scvfCorners, Scalar q)
Calculates the integration point on an scvf.
Definition: helper.hh:530