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