28 #ifndef EWOMS_VCFV_STENCIL_HH 29 #define EWOMS_VCFV_STENCIL_HH 33 #include <opm/common/Unused.hpp> 35 #include <opm/common/Exceptions.hpp> 36 #include <opm/common/ErrorMacros.hpp> 38 #include <dune/grid/common/intersectioniterator.hh> 39 #include <dune/grid/common/mcmgmapper.hh> 40 #include <dune/geometry/referenceelements.hh> 42 #if HAVE_DUNE_LOCALFUNCTIONS 43 #include <dune/localfunctions/lagrange/pqkfactory.hh> 44 #endif // HAVE_DUNE_LOCALFUNCTIONS 46 #include <dune/common/version.hh> 55 template <
class Scalar,
unsigned dim,
unsigned basicGeomType>
56 class VcfvScvGeometries;
61 template <
class Scalar>
62 class VcfvScvGeometries<Scalar, 1, Dune::GeometryType::cube>
66 static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::simplex;
74 Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] = {
85 for (
unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx)
86 scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
89 static const ScvLocalGeometry&
get(
unsigned scvIdx)
90 {
return scvGeoms_[scvIdx]; }
93 static ScvLocalGeometry scvGeoms_[numScv];
96 template <
class Scalar>
97 typename VcfvScvGeometries<Scalar, 1, Dune::GeometryType::cube>::ScvLocalGeometry
98 VcfvScvGeometries<Scalar, 1, Dune::GeometryType::cube>::scvGeoms_[
99 VcfvScvGeometries<Scalar, 1, Dune::GeometryType::cube>::numScv];
101 template <
class Scalar>
102 class VcfvScvGeometries<Scalar, 1, Dune::GeometryType::simplex>
106 static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::simplex;
111 static const ScvLocalGeometry&
get(
unsigned scvIdx OPM_UNUSED)
113 OPM_THROW(std::logic_error,
114 "Not implemented: VcfvScvGeometries<Scalar, 1, Dune::GeometryType::simplex>");
121 template <
class Scalar>
122 class VcfvScvGeometries<Scalar, 2, Dune::GeometryType::simplex>
126 static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::simplex;
131 static const ScvLocalGeometry&
get(
unsigned scvIdx)
132 {
return scvGeoms_[scvIdx]; }
137 Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
161 for (
unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx)
162 scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
166 static ScvLocalGeometry scvGeoms_[numScv];
169 template <
class Scalar>
170 typename VcfvScvGeometries<Scalar, 2, Dune::GeometryType::simplex>::ScvLocalGeometry
171 VcfvScvGeometries<Scalar, 2, Dune::GeometryType::simplex>::scvGeoms_[
172 VcfvScvGeometries<Scalar, 2, Dune::GeometryType::simplex>::numScv];
174 template <
class Scalar>
175 class VcfvScvGeometries<Scalar, 2, Dune::GeometryType::cube>
179 static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::cube;
184 static const ScvLocalGeometry&
get(
unsigned scvIdx)
185 {
return scvGeoms_[scvIdx]; }
190 Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
221 for (
unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx)
222 scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
225 static ScvLocalGeometry scvGeoms_[numScv];
228 template <
class Scalar>
229 typename VcfvScvGeometries<Scalar, 2, Dune::GeometryType::cube>::ScvLocalGeometry
230 VcfvScvGeometries<Scalar, 2, Dune::GeometryType::cube>::scvGeoms_[
231 VcfvScvGeometries<Scalar, 2, Dune::GeometryType::cube>::numScv];
236 template <
class Scalar>
237 class VcfvScvGeometries<Scalar, 3, Dune::GeometryType::simplex>
241 static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::simplex;
246 static const ScvLocalGeometry&
get(
unsigned scvIdx)
247 {
return scvGeoms_[scvIdx]; }
252 Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
258 { 1.0/3, 1.0/3, 0.0 },
261 { 1.0/3, 0.0, 1.0/3 },
262 { 0.0, 1.0/3, 1.0/3 },
263 { 1.0/4, 1.0/4, 1.0/4 },
269 { 1.0/3, 1.0/3, 0.0 },
270 { 1.0/2, 1.0/2, 0.0 },
272 { 1.0/3, 0.0, 1.0/3 },
273 { 1.0/2, 0.0, 1.0/2 },
274 { 1.0/4, 1.0/4, 1.0/4 },
275 { 1.0/3, 1.0/3, 1.0/3 },
280 { 1.0/3, 1.0/3, 0.0 },
282 { 1.0/2, 1.0/2, 0.0 },
284 { 0.0, 1.0/3, 1.0/3 },
285 { 1.0/4, 1.0/4, 1.0/4 },
286 { 0.0, 1.0/2, 1.0/2 },
287 { 1.0/3, 1.0/3, 1.0/3 },
292 { 1.0/3, 0.0, 1.0/3 },
293 { 0.0, 1.0/3, 1.0/3 },
294 { 1.0/4, 1.0/4, 1.0/4 },
297 { 1.0/2, 0.0, 1.0/2 },
298 { 0.0, 1.0/2, 1.0/2 },
299 { 1.0/3, 1.0/3, 1.0/3 },
303 for (
unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx)
304 scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
308 static ScvLocalGeometry scvGeoms_[numScv];
311 template <
class Scalar>
312 typename VcfvScvGeometries<Scalar, 3, Dune::GeometryType::simplex>::ScvLocalGeometry
313 VcfvScvGeometries<Scalar, 3, Dune::GeometryType::simplex>::scvGeoms_[
314 VcfvScvGeometries<Scalar, 3, Dune::GeometryType::simplex>::numScv];
316 template <
class Scalar>
317 class VcfvScvGeometries<Scalar, 3, Dune::GeometryType::cube>
321 static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::cube;
326 static const ScvLocalGeometry&
get(
unsigned scvIdx)
327 {
return scvGeoms_[scvIdx]; }
332 Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
338 { 1.0/2, 1.0/2, 0.0 },
341 { 1.0/2, 0.0, 1.0/2 },
342 { 0.0, 1.0/2, 1.0/2 },
343 { 1.0/2, 1.0/2, 1.0/2 },
349 { 1.0/2, 1.0/2, 0.0 },
352 { 1.0/2, 0.0, 1.0/2 },
354 { 1.0/2, 1.0/2, 1.0/2 },
355 { 1.0, 1.0/2, 1.0/2 },
360 { 1.0/2, 1.0/2, 0.0 },
364 { 0.0, 1.0/2, 1.0/2 },
365 { 1.0/2, 1.0/2, 1.0/2 },
367 { 1.0/2, 1.0, 1.0/2 },
371 { 1.0/2, 1.0/2, 0.0 },
376 { 1.0/2, 1.0/2, 1.0/2 },
377 { 1.0, 1.0/2, 1.0/2 },
378 { 1.0/2, 1.0, 1.0/2 },
384 { 1.0/2, 0.0, 1.0/2 },
385 { 0.0, 1.0/2, 1.0/2 },
386 { 1.0/2, 1.0/2, 1.0/2 },
391 { 1.0/2, 1.0/2, 1.0 },
395 { 1.0/2, 0.0, 1.0/2 },
397 { 1.0/2, 1.0/2, 1.0/2 },
398 { 1.0, 1.0/2, 1.0/2 },
402 { 1.0/2, 1.0/2, 1.0 },
407 { 0.0, 1.0/2, 1.0/2 },
408 { 1.0/2, 1.0/2, 1.0/2 },
410 { 1.0/2, 1.0, 1.0/2 },
413 { 1.0/2, 1.0/2, 1.0 },
419 { 1.0/2, 1.0/2, 1.0/2 },
420 { 1.0, 1.0/2, 1.0/2 },
421 { 1.0/2, 1.0, 1.0/2 },
424 { 1.0/2, 1.0/2, 1.0 },
431 for (
unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx)
432 scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
435 static ScvLocalGeometry scvGeoms_[numScv];
438 template <
class Scalar>
439 typename VcfvScvGeometries<Scalar, 3, Dune::GeometryType::cube>::ScvLocalGeometry
440 VcfvScvGeometries<Scalar, 3, Dune::GeometryType::cube>::scvGeoms_[
441 VcfvScvGeometries<Scalar, 3, Dune::GeometryType::cube>::numScv];
466 template <
class Scalar,
class Gr
idView>
469 enum{dim = GridView::dimension};
470 enum{dimWorld = GridView::dimensionworld};
471 enum{maxNC = (dim < 3 ? 4 : 8)};
472 enum{maxNE = (dim < 3 ? 4 : 12)};
473 enum{maxNF = (dim < 3 ? 1 : 6)};
474 enum{maxBF = (dim < 3 ? 8 : 24)};
475 typedef typename GridView::ctype CoordScalar;
476 typedef typename GridView::Traits::template Codim<0>::Entity Element;
478 typedef typename GridView::Traits::template Codim<dim>::Entity Entity;
480 typedef typename Element::Geometry Geometry;
481 typedef Dune::FieldVector<Scalar,dimWorld> DimVector;
482 typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
483 typedef Dune::FieldVector<CoordScalar,dim> LocalPosition;
484 typedef typename GridView::IntersectionIterator IntersectionIterator;
488 #if HAVE_DUNE_LOCALFUNCTIONS 489 typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> LocalFiniteElementCache;
490 typedef typename LocalFiniteElementCache::FiniteElementType LocalFiniteElement;
491 typedef typename LocalFiniteElement::Traits::LocalBasisType::Traits LocalBasisTraits;
492 typedef typename LocalBasisTraits::JacobianType ShapeJacobian;
493 #endif // HAVE_DUNE_LOCALFUNCTIONS 495 Scalar quadrilateralArea(
const GlobalPosition& p0,
496 const GlobalPosition& p1,
497 const GlobalPosition& p2,
498 const GlobalPosition& p3)
500 return 0.5*std::abs((p3[0] - p1[0])*(p2[1] - p0[1]) - (p3[1] - p1[1])*(p2[0] - p0[0]));
503 void crossProduct(DimVector& c,
const DimVector& a,
const DimVector& b)
505 c[0] = a[1]*b[2] - a[2]*b[1];
506 c[1] = a[2]*b[0] - a[0]*b[2];
507 c[2] = a[0]*b[1] - a[1]*b[0];
510 Scalar pyramidVolume(
const GlobalPosition& p0,
511 const GlobalPosition& p1,
512 const GlobalPosition& p2,
513 const GlobalPosition& p3,
514 const GlobalPosition& p4)
516 DimVector a(p2); a -= p0;
517 DimVector b(p3); b -= p1;
520 crossProduct(n, a, b);
524 return 1.0/6.0*(n*a);
527 Scalar prismVolume(
const GlobalPosition& p0,
528 const GlobalPosition& p1,
529 const GlobalPosition& p2,
530 const GlobalPosition& p3,
531 const GlobalPosition& p4,
532 const GlobalPosition& p5)
535 for (
unsigned k = 0; k < dimWorld; ++k)
538 for (
unsigned k = 0; k < dimWorld; ++k)
541 crossProduct(m, a, b);
543 for (
unsigned k = 0; k < dimWorld; ++k)
544 a[k] = p1[k] - p0[k];
545 for (
unsigned k = 0; k < dimWorld; ++k)
546 b[k] = p2[k] - p0[k];
548 crossProduct(n, a, b);
551 for (
unsigned k = 0; k < dimWorld; ++k)
552 a[k] = p5[k] - p0[k];
554 return std::abs(1.0/6.0*(n*a));
557 Scalar hexahedronVolume(
const GlobalPosition& p0,
558 const GlobalPosition& p1,
559 const GlobalPosition& p2,
560 const GlobalPosition& p3,
561 const GlobalPosition& p4,
562 const GlobalPosition& p5,
563 const GlobalPosition& p6,
564 const GlobalPosition& p7)
567 prismVolume(p0,p1,p2,p4,p5,p6)
568 + prismVolume(p0,p2,p3,p4,p6,p7);
571 void normalOfQuadrilateral3D(DimVector& normal,
572 const GlobalPosition& p0,
573 const GlobalPosition& p1,
574 const GlobalPosition& p2,
575 const GlobalPosition& p3)
578 for (
unsigned k = 0; k < dimWorld; ++k)
581 for (
unsigned k = 0; k < dimWorld; ++k)
584 crossProduct(normal, a, b);
588 Scalar quadrilateralArea3D(
const GlobalPosition& p0,
589 const GlobalPosition& p1,
590 const GlobalPosition& p2,
591 const GlobalPosition& p3)
594 normalOfQuadrilateral3D(normal, p0, p1, p2, p3);
595 return normal.two_norm();
598 void getFaceIndices(
unsigned numElemVertices,
unsigned k,
unsigned& leftFace,
unsigned& rightFace)
600 static const unsigned edgeToFaceTet[2][6] = {
604 static const unsigned edgeToFacePyramid[2][8] = {
605 {1, 2, 3, 4, 1, 3, 4, 2},
606 {0, 0, 0, 0, 3, 2, 1, 4}
608 static const unsigned edgeToFacePrism[2][9] = {
609 {1, 0, 2, 0, 1, 2, 4, 4, 4},
610 {0, 2, 1, 3, 3, 3, 0, 1, 2}
612 static const unsigned edgeToFaceHex[2][12] = {
613 {0, 2, 3, 1, 4, 1, 2, 4, 0, 5, 5, 3},
614 {2, 1, 0, 3, 0, 4, 4, 3, 5, 1, 2, 5}
617 switch (numElemVertices) {
619 leftFace = edgeToFaceTet[0][k];
620 rightFace = edgeToFaceTet[1][k];
623 leftFace = edgeToFacePyramid[0][k];
624 rightFace = edgeToFacePyramid[1][k];
627 leftFace = edgeToFacePrism[0][k];
628 rightFace = edgeToFacePrism[1][k];
631 leftFace = edgeToFaceHex[0][k];
632 rightFace = edgeToFaceHex[1][k];
635 OPM_THROW(std::logic_error,
636 "Not implemented: VcfvStencil::getFaceIndices for " 637 << numElemVertices <<
" vertices");
642 void getEdgeIndices(
unsigned numElemVertices,
unsigned face,
unsigned vert,
unsigned& leftEdge,
unsigned& rightEdge)
644 static const int faceAndVertexToLeftEdgeTet[4][4] = {
650 static const int faceAndVertexToRightEdgeTet[4][4] = {
656 static const int faceAndVertexToLeftEdgePyramid[5][5] = {
663 static const int faceAndVertexToRightEdgePyramid[5][5] = {
670 static const int faceAndVertexToLeftEdgePrism[5][6] = {
671 { 3, 3, -1, 0, 1, -1},
672 { 4, -1, 4, 0, -1, 2},
673 {-1, 5, 5, -1, 1, 2},
674 { 3, 3, 5, -1, -1, -1},
675 {-1, -1, -1, 6, 6, 8}
677 static const int faceAndVertexToRightEdgePrism[5][6] = {
678 { 0, 1, -1, 6, 6, -1},
679 { 0, -1, 2, 7, -1, 7},
680 {-1, 1, 2, -1, 8, 8},
681 { 4, 5, 4, -1, -1, -1},
682 {-1, -1, -1, 7, 8, 7}
684 static const int faceAndVertexToLeftEdgeHex[6][8] = {
685 { 0, -1, 4, -1, 8, -1, 2, -1},
686 {-1, 5, -1, 3, -1, 1, -1, 9},
687 { 6, 1, -1, -1, 0, 10, -1, -1},
688 {-1, -1, 2, 7, -1, -1, 11, 3},
689 { 4, 6, 7, 5, -1, -1, -1, -1},
690 {-1, -1, -1, -1, 10, 9, 8, 11}
692 static const int faceAndVertexToRightEdgeHex[6][8] = {
693 { 4, -1, 2, -1, 0, -1, 8, -1},
694 {-1, 1, -1, 5, -1, 9, -1, 3},
695 { 0, 6, -1, -1, 10, 1, -1, -1},
696 {-1, -1, 7, 3, -1, -1, 2, 11},
697 { 6, 5, 4, 7, -1, -1, -1, -1},
698 {-1, -1, -1, -1, 8, 10, 11, 9}
701 switch (numElemVertices) {
703 leftEdge =
static_cast<unsigned>(faceAndVertexToLeftEdgeTet[face][vert]);
704 rightEdge =
static_cast<unsigned>(faceAndVertexToRightEdgeTet[face][vert]);
707 leftEdge =
static_cast<unsigned>(faceAndVertexToLeftEdgePyramid[face][vert]);
708 rightEdge =
static_cast<unsigned>(faceAndVertexToRightEdgePyramid[face][vert]);
711 leftEdge =
static_cast<unsigned>(faceAndVertexToLeftEdgePrism[face][vert]);
712 rightEdge =
static_cast<unsigned>(faceAndVertexToRightEdgePrism[face][vert]);
715 leftEdge =
static_cast<unsigned>(faceAndVertexToLeftEdgeHex[face][vert]);
716 rightEdge =
static_cast<unsigned>(faceAndVertexToRightEdgeHex[face][vert]);
719 OPM_THROW(std::logic_error,
720 "Not implemented: VcfvStencil::getFaceIndices for " 721 << numElemVertices <<
" vertices");
727 typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView,
728 Dune::MCMGVertexLayout > VertexMapper;
735 const GlobalPosition center()
const 736 {
return global(localGeometry_->
center()); }
738 const GlobalPosition corner(
unsigned cornerIdx)
const 739 {
return global(localGeometry_->
corner(cornerIdx)); }
741 const GlobalPosition global(
const LocalPosition& localPos)
const 742 {
return element_->geometry().global(localPos); }
745 {
return *localGeometry_; }
748 const Element *element_;
753 const GlobalPosition& globalPos()
const 756 const GlobalPosition center()
const 759 Scalar volume()
const 783 const DimVector& normal()
const 786 unsigned short interiorIndex()
const 789 unsigned short exteriorIndex()
const 795 const LocalPosition& localPos()
const 798 const GlobalPosition& integrationPos()
const 816 VcfvStencil(
const GridView& gridView,
const VertexMapper& vertexMapper)
817 : gridView_(gridView)
818 , vertexMapper_( vertexMapper )
819 , element_(*gridView.template begin<0>())
821 static bool localGeometriesInitialized =
false;
822 if (!localGeometriesInitialized) {
823 localGeometriesInitialized =
true;
825 VcfvScvGeometries<Scalar, 1, Dune::GeometryType::cube>::init();
826 VcfvScvGeometries<Scalar, 2, Dune::GeometryType::cube>::init();
827 VcfvScvGeometries<Scalar, 2, Dune::GeometryType::simplex>::init();
828 VcfvScvGeometries<Scalar, 3, Dune::GeometryType::cube>::init();
829 VcfvScvGeometries<Scalar, 3, Dune::GeometryType::simplex>::init();
842 numVertices = e.subEntities(dim);
843 numEdges = e.subEntities(dim-1);
844 numFaces = (dim<3)?0:e.subEntities(1);
846 numBoundarySegments_ = 0;
849 const Geometry& geometry = e.geometry();
850 geometryType_ = geometry.type();
851 const typename Dune::ReferenceElementContainer<CoordScalar,dim>::value_type&
852 referenceElement = Dune::ReferenceElements<CoordScalar,dim>::general(geometryType_);
853 for (
unsigned vertexIdx = 0; vertexIdx < numVertices; vertexIdx++) {
854 subContVol[vertexIdx].
local = referenceElement.position(static_cast<int>(vertexIdx), dim);
855 subContVol[vertexIdx].
global = geometry.corner(static_cast<int>(vertexIdx));
859 void updatePrimaryTopology(
const Element& element)
867 void update(
const Element& e)
871 const Geometry& geometry = e.geometry();
872 geometryType_ = geometry.type();
874 const typename Dune::ReferenceElementContainer<CoordScalar,dim>::value_type&
875 referenceElement = Dune::ReferenceElements<CoordScalar,dim>::general(geometryType_);
877 elementVolume = geometry.volume();
878 elementLocal = referenceElement.position(0,0);
879 elementGlobal = geometry.global(elementLocal);
882 for (
unsigned vert = 0; vert < numVertices; vert++) {
883 subContVol[vert].
local = referenceElement.position(static_cast<int>(vert), dim);
884 subContVol[vert].
global = geometry.global(subContVol[vert].local);
888 for (
unsigned edge = 0; edge < numEdges; edge++) {
889 edgeCoord[edge] = geometry.global(referenceElement.position(static_cast<int>(edge), dim-1));
893 for (
unsigned face = 0; face < numFaces; face++) {
894 faceCoord[face] = geometry.global(referenceElement.position(static_cast<int>(face), 1));
902 fillSubContVolData_();
905 for (
unsigned k = 0; k < numEdges; k++) {
906 unsigned short i =
static_cast<unsigned short>(referenceElement.subEntity(static_cast<int>(k), dim-1, 0, dim));
907 unsigned short j =
static_cast<unsigned short>(referenceElement.subEntity(static_cast<int>(k), dim-1, 1, dim));
908 if (numEdges == 4 && (i == 2 || j == 2))
910 subContVolFace[k].
i = i;
911 subContVolFace[k].j = j;
918 LocalPosition ipLocal_;
922 subContVolFace[k].
normal_ = 1.0;
923 subContVolFace[k].
area_ = 1.0;
924 ipLocal_ = subContVolFace[k].
ipLocal_;
927 ipLocal_ = referenceElement.position(static_cast<int>(k), dim-1) + elementLocal;
929 subContVolFace[k].
ipLocal_ = ipLocal_;
930 for (
unsigned m = 0; m < dimWorld; ++m)
931 diffVec[m] = elementGlobal[m] - edgeCoord[k][m];
932 subContVolFace[k].
normal_[0] = diffVec[1];
933 subContVolFace[k].
normal_[1] = -diffVec[0];
935 for (
unsigned m = 0; m < dimWorld; ++m)
936 diffVec[m] = subContVol[j].global[m] - subContVol[i].global[m];
938 if (subContVolFace[k].normal_ * diffVec < 0)
939 subContVolFace[k].
normal_ *= -1;
941 subContVolFace[k].
area_ = subContVolFace[k].
normal_.two_norm();
947 getFaceIndices(numVertices, k, leftFace, rightFace);
948 ipLocal_ = referenceElement.position(static_cast<int>(k), dim-1) + elementLocal
949 + referenceElement.position(static_cast<int>(leftFace), 1)
950 + referenceElement.position(static_cast<int>(rightFace), 1);
952 subContVolFace[k].
ipLocal_ = ipLocal_;
953 normalOfQuadrilateral3D(subContVolFace[k].normal_,
954 edgeCoord[k], faceCoord[rightFace],
955 elementGlobal, faceCoord[leftFace]);
956 subContVolFace[k].
area_ = subContVolFace[k].
normal_.two_norm();
961 subContVolFace[k].
ipGlobal_ = geometry.global(ipLocal_);
965 IntersectionIterator endit = gridView_.iend(e);
966 for (IntersectionIterator it = gridView_.ibegin(e); it != endit; ++it) {
967 const typename IntersectionIterator::Intersection& intersection = *it ;
969 if ( ! intersection.boundary())
972 unsigned face =
static_cast<unsigned>(intersection.indexInInside());
973 unsigned numVerticesOfFace =
static_cast<unsigned>(referenceElement.size(static_cast<int>(face), 1, dim));
974 for (
unsigned vertInFace = 0; vertInFace < numVerticesOfFace; vertInFace++)
976 unsigned short vertInElement =
static_cast<unsigned short>(referenceElement.subEntity(static_cast<int>(face), 1, static_cast<int>(vertInFace), dim));
977 unsigned bfIdx = numBoundarySegments_;
978 ++numBoundarySegments_;
981 boundaryFace_[bfIdx].
ipLocal_ = referenceElement.position(static_cast<int>(vertInElement), dim);
982 boundaryFace_[bfIdx].
area_ = 1.0;
985 boundaryFace_[bfIdx].
ipLocal_ = referenceElement.position(static_cast<int>(vertInElement), dim)
986 + referenceElement.position(static_cast<int>(face), 1);
987 boundaryFace_[bfIdx].
ipLocal_ *= 0.5;
988 boundaryFace_[bfIdx].
area_ = 0.5 * intersection.geometry().volume();
993 getEdgeIndices(numVertices, face, vertInElement, leftEdge, rightEdge);
994 boundaryFace_[bfIdx].
ipLocal_ = referenceElement.position(static_cast<int>(vertInElement), dim)
995 + referenceElement.position(static_cast<int>(face), 1)
996 + referenceElement.position(static_cast<int>(leftEdge), dim-1)
997 + referenceElement.position(static_cast<int>(rightEdge), dim-1);
998 boundaryFace_[bfIdx].
ipLocal_ *= 0.25;
999 boundaryFace_[bfIdx].
area_ =
1000 quadrilateralArea3D(subContVol[vertInElement].global,
1001 edgeCoord[rightEdge],
1003 edgeCoord[leftEdge]);
1006 OPM_THROW(std::logic_error,
"Not implemented:VcfvStencil for dim = " << dim);
1008 boundaryFace_[bfIdx].
ipGlobal_ = geometry.global(boundaryFace_[bfIdx].ipLocal_);
1009 boundaryFace_[bfIdx].
i = vertInElement;
1010 boundaryFace_[bfIdx].j = vertInElement;
1013 boundaryFace_[bfIdx].
normal_ = intersection.centerUnitOuterNormal();
1017 updateScvGeometry(e);
1020 void updateScvGeometry(
const Element& element)
1022 auto geomType = element.geometry().type();
1025 if (geomType.isTriangle() || geomType.isTetrahedron()) {
1026 for (
unsigned vertIdx = 0; vertIdx < numVertices; ++vertIdx) {
1027 subContVol[vertIdx].
geometry_.element_ = &element;
1028 subContVol[vertIdx].
geometry_.localGeometry_ =
1029 &VcfvScvGeometries<Scalar, dim, Dune::GeometryType::simplex>::get(vertIdx);
1032 else if (geomType.isLine() || geomType.isQuadrilateral() || geomType.isHexahedron()) {
1033 for (
unsigned vertIdx = 0; vertIdx < numVertices; ++vertIdx) {
1034 subContVol[vertIdx].
geometry_.element_ = &element;
1035 subContVol[vertIdx].
geometry_.localGeometry_ =
1036 &VcfvScvGeometries<Scalar, dim, Dune::GeometryType::cube>::get(vertIdx);
1040 OPM_THROW(std::logic_error,
1041 "Not implemented: SCV geometries for non hexahedron elements");
1044 #if HAVE_DUNE_LOCALFUNCTIONS 1045 void updateCenterGradients()
1047 const auto& localFiniteElement = feCache_.get(element_.type());
1048 const auto& geom = element_.geometry();
1050 std::vector<ShapeJacobian> localJac;
1052 for (
unsigned scvIdx = 0; scvIdx < numVertices; ++ scvIdx) {
1053 const auto& localCenter = subContVol[scvIdx].localGeometry().
center();
1054 localFiniteElement.localBasis().evaluateJacobian(localCenter, localJac);
1055 const auto& globalPos = subContVol[scvIdx].geometry().center();
1057 const auto jacInvT = geom.jacobianInverseTransposed(globalPos);
1058 for (
unsigned vert = 0; vert < numVertices; vert++) {
1059 jacInvT.mv(localJac[vert][0], subContVol[scvIdx].gradCenter[vert]);
1065 unsigned numDof()
const 1066 {
return numVertices; }
1068 unsigned numPrimaryDof()
const 1069 {
return numDof(); }
1071 Dune::PartitionType partitionType(
unsigned scvIdx)
const 1072 {
return element_.template subEntity<dim>(scvIdx)->partitionType(); }
1074 const SubControlVolume& subControlVolume(
unsigned scvIdx)
const 1076 assert(0 <= scvIdx && scvIdx < numDof());
1077 return subContVol[scvIdx];
1080 unsigned numInteriorFaces()
const 1081 {
return numEdges; }
1083 unsigned numBoundaryFaces()
const 1084 {
return numBoundarySegments_; }
1086 const SubControlVolumeFace& interiorFace(
unsigned faceIdx)
const 1087 {
return subContVolFace[faceIdx]; }
1090 {
return boundaryFace_[bfIdx]; }
1098 assert(0 <= dofIdx && dofIdx < numDof());
1100 return static_cast<unsigned>(vertexMapper_.subIndex(element_, static_cast<int>(dofIdx), dim));
1109 assert(0 <= dofIdx && dofIdx < numDof());
1110 return element_.template subEntity<dim>(
static_cast<int>(dofIdx));
1114 #if __GNUC__ || __clang__ 1115 #pragma GCC diagnostic push 1116 #pragma GCC diagnostic ignored "-Wpragmas" 1117 #pragma GCC diagnostic ignored "-Warray-bounds" 1119 void fillSubContVolData_()
1123 subContVol[0].
volume_ = 0.5*elementVolume;
1124 subContVol[1].
volume_ = 0.5*elementVolume;
1126 else if (dim == 2) {
1127 switch (numVertices) {
1129 subContVol[0].
volume_ = elementVolume/3;
1130 subContVol[1].
volume_ = elementVolume/3;
1131 subContVol[2].
volume_ = elementVolume/3;
1135 quadrilateralArea(subContVol[0].global,
1140 quadrilateralArea(subContVol[1].global,
1145 quadrilateralArea(subContVol[2].global,
1150 quadrilateralArea(subContVol[3].global,
1156 OPM_THROW(std::logic_error,
1157 "Not implemented:VcfvStencil dim = " << dim
1158 <<
", numVertices = " << numVertices);
1161 else if (dim == 3) {
1162 switch (numVertices) {
1164 for (
unsigned k = 0; k < numVertices; k++)
1165 subContVol[k].volume_ = elementVolume / 4.0;
1169 hexahedronVolume(subContVol[0].global,
1178 hexahedronVolume(subContVol[1].global,
1187 hexahedronVolume(subContVol[2].global,
1196 hexahedronVolume(subContVol[3].global,
1204 subContVol[4].
volume_ = elementVolume -
1212 hexahedronVolume(subContVol[0].global,
1221 hexahedronVolume(subContVol[1].global,
1230 hexahedronVolume(subContVol[2].global,
1239 hexahedronVolume(edgeCoord[0],
1243 subContVol[3].global,
1248 hexahedronVolume(edgeCoord[1],
1252 subContVol[4].global,
1257 hexahedronVolume(edgeCoord[2],
1261 subContVol[5].global,
1268 hexahedronVolume(subContVol[0].global,
1277 hexahedronVolume(subContVol[1].global,
1286 hexahedronVolume(subContVol[2].global,
1295 hexahedronVolume(subContVol[3].global,
1304 hexahedronVolume(edgeCoord[0],
1308 subContVol[4].global,
1313 hexahedronVolume(edgeCoord[1],
1317 subContVol[5].global,
1322 hexahedronVolume(edgeCoord[2],
1326 subContVol[6].global,
1331 hexahedronVolume(edgeCoord[3],
1335 subContVol[7].global,
1341 OPM_THROW(std::logic_error,
1342 "Not implemented:VcfvStencil for dim = " << dim
1343 <<
", numVertices = " << numVertices);
1347 OPM_THROW(std::logic_error,
"Not implemented:VcfvStencil for dim = " << dim);
1349 #if __GNUC__ || __clang__ 1350 #pragma GCC diagnostic pop 1353 const GridView& gridView_;
1354 const VertexMapper& vertexMapper_;
1358 #if HAVE_DUNE_LOCALFUNCTIONS 1359 static LocalFiniteElementCache feCache_;
1360 #endif // HAVE_DUNE_LOCALFUNCTIONS 1363 LocalPosition elementLocal;
1365 GlobalPosition elementGlobal;
1367 Scalar elementVolume;
1369 SubControlVolume subContVol[maxNC];
1371 SubControlVolumeFace subContVolFace[maxNE];
1374 unsigned numBoundarySegments_;
1376 GlobalPosition edgeCoord[maxNE];
1378 GlobalPosition faceCoord[maxNF];
1380 unsigned numVertices;
1385 Dune::GeometryType geometryType_;
1388 #if HAVE_DUNE_LOCALFUNCTIONS 1389 template<
class Scalar,
class Gr
idView>
1390 typename VcfvStencil<Scalar, GridView>::LocalFiniteElementCache
1391 VcfvStencil<Scalar, GridView>::feCache_;
1392 #endif // HAVE_DUNE_LOCALFUNCTIONS DimVector normal_
normal on face pointing to CV j or outward of the domain with length equal to |scvf| ...
Definition: vcfvstencil.hh:808
Definition: baseauxiliarymodule.hh:37
const GlobalPosition & corner(unsigned cornerIdx) const
Return the position of the corner with a given index.
Definition: quadraturegeometries.hh:127
LocalPosition local
local vertex position
Definition: vcfvstencil.hh:769
GlobalPosition ipGlobal_
integration point in global coords
Definition: vcfvstencil.hh:806
Quadrature geometry for quadrilaterals.
Definition: quadraturegeometries.hh:39
finite volume intersected with element
Definition: vcfvstencil.hh:751
Scalar area_
area of face
Definition: vcfvstencil.hh:810
GlobalPosition global
global vertex position
Definition: vcfvstencil.hh:771
Entity entity(unsigned dofIdx) const
Return the global space index given the index of a degree of freedom.
Definition: vcfvstencil.hh:1107
unsigned short i
scvf seperates corner i and j of elem
Definition: vcfvstencil.hh:802
void updateTopology(const Element &e)
Update the non-geometric part of the stencil.
Definition: vcfvstencil.hh:838
interior face of a sub control volume
Definition: vcfvstencil.hh:781
ScvGeometry geometry_
The geometry of the sub-control volume in local coordinates.
Definition: vcfvstencil.hh:775
Definition: vcfvstencil.hh:732
const GlobalPosition & center() const
Returns the center of weight of the polyhedron.
Definition: quadraturegeometries.hh:69
LocalPosition ipLocal_
integration point in local coords
Definition: vcfvstencil.hh:804
Quadrature geometry for quadrilaterals.
VertexMapper Mapper
exported Mapper type
Definition: vcfvstencil.hh:730
SubControlVolumeFace BoundaryFace
compatibility typedef
Definition: vcfvstencil.hh:814
Represents the finite volume geometry of a single element in the VCFV discretization.
Definition: vcfvstencil.hh:467
Dune::FieldVector< DimVector, maxNC > gradCenter
derivative of shape function at the center of the sub control volume
Definition: vcfvstencil.hh:778
unsigned globalSpaceIndex(unsigned dofIdx) const
Return the global space index given the index of a degree of freedom.
Definition: vcfvstencil.hh:1096
Scalar volume_
space occupied by the sub-control volume
Definition: vcfvstencil.hh:773