Go to the documentation of this file.
27 #ifndef EWOMS_VCFV_STENCIL_HH
28 #define EWOMS_VCFV_STENCIL_HH
32 #include <opm/common/Exceptions.hpp>
33 #include <opm/common/ErrorMacros.hpp>
35 #include <dune/grid/common/intersectioniterator.hh>
36 #include <dune/grid/common/mcmgmapper.hh>
37 #include <dune/geometry/referenceelements.hh>
38 #include <dune/localfunctions/lagrange/pqkfactory.hh>
39 #include <dune/common/version.hh>
48 template < class Scalar, int dim, int basicGeomType>
49 class VcfvScvGeometries;
54 template < class Scalar>
55 class VcfvScvGeometries<Scalar, 1, Dune::GeometryType::cube>
59 static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::simplex;
67 Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] = {
78 for ( unsigned scvIdx = 0; scvIdx < numScv; ++scvIdx)
79 scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
84 static const ScvLocalGeometry & get( int scvIdx)
85 { return scvGeoms_[scvIdx]; }
88 static ScvLocalGeometry scvGeoms_[numScv];
91 template < class Scalar>
92 typename VcfvScvGeometries<Scalar, 1, Dune::GeometryType::cube>::ScvLocalGeometry
93 VcfvScvGeometries<Scalar, 1, Dune::GeometryType::cube>::scvGeoms_[
94 VcfvScvGeometries<Scalar, 1, Dune::GeometryType::cube>::numScv];
96 template < class Scalar>
97 class VcfvScvGeometries<Scalar, 1, Dune::GeometryType::simplex>
101 static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::simplex;
106 static const ScvLocalGeometry & get( int scvIdx)
107 { OPM_THROW(std::logic_error,
109 "VcfvScvGeometries<Scalar, 1, Dune::GeometryType::simplex>"); }
115 template < class Scalar>
116 class VcfvScvGeometries<Scalar, 2, Dune::GeometryType::simplex>
120 static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::simplex;
125 static const ScvLocalGeometry & get( int scvIdx)
126 { return scvGeoms_[scvIdx]; }
131 Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
155 for ( int scvIdx = 0; scvIdx < numScv; ++scvIdx)
156 scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
162 static ScvLocalGeometry scvGeoms_[numScv];
165 template < class Scalar>
166 typename VcfvScvGeometries<Scalar, 2, Dune::GeometryType::simplex>::ScvLocalGeometry
167 VcfvScvGeometries<Scalar, 2, Dune::GeometryType::simplex>::scvGeoms_[
168 VcfvScvGeometries<Scalar, 2, Dune::GeometryType::simplex>::numScv];
170 template < class Scalar>
171 class VcfvScvGeometries<Scalar, 2, Dune::GeometryType::cube>
175 static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::cube;
180 static const ScvLocalGeometry & get( int scvIdx)
181 { return scvGeoms_[scvIdx]; }
186 Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
217 for ( int scvIdx = 0; scvIdx < numScv; ++scvIdx)
218 scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
223 static ScvLocalGeometry scvGeoms_[numScv];
226 template < class Scalar>
227 typename VcfvScvGeometries<Scalar, 2, Dune::GeometryType::cube>::ScvLocalGeometry
228 VcfvScvGeometries<Scalar, 2, Dune::GeometryType::cube>::scvGeoms_[
229 VcfvScvGeometries<Scalar, 2, Dune::GeometryType::cube>::numScv];
234 template < class Scalar>
235 class VcfvScvGeometries<Scalar, 3, Dune::GeometryType::simplex>
239 static const Dune::GeometryType::BasicType basicType = Dune::GeometryType::simplex;
244 static const ScvLocalGeometry & get( int scvIdx)
245 { return scvGeoms_[scvIdx]; }
250 Scalar scvCorners[numScv][ScvLocalGeometry::numCorners][dim] =
256 { 1.0/3, 1.0/3, 0.0 },
259 { 1.0/3, 0.0, 1.0/3 },
260 { 0.0, 1.0/3, 1.0/3 },
261 { 1.0/4, 1.0/4, 1.0/4 },
267 { 1.0/3, 1.0/3, 0.0 },
268 { 1.0/2, 1.0/2, 0.0 },
270 { 1.0/3, 0.0, 1.0/3 },
271 { 1.0/2, 0.0, 1.0/2 },
272 { 1.0/4, 1.0/4, 1.0/4 },
273 { 1.0/3, 1.0/3, 1.0/3 },
278 { 1.0/3, 1.0/3, 0.0 },
280 { 1.0/2, 1.0/2, 0.0 },
282 { 0.0, 1.0/3, 1.0/3 },
283 { 1.0/4, 1.0/4, 1.0/4 },
284 { 0.0, 1.0/2, 1.0/2 },
285 { 1.0/3, 1.0/3, 1.0/3 },
290 { 1.0/3, 0.0, 1.0/3 },
291 { 0.0, 1.0/3, 1.0/3 },
292 { 1.0/4, 1.0/4, 1.0/4 },
295 { 1.0/2, 0.0, 1.0/2 },
296 { 0.0, 1.0/2, 1.0/2 },
297 { 1.0/3, 1.0/3, 1.0/3 },
301 for ( int scvIdx = 0; scvIdx < numScv; ++scvIdx)
302 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( int 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 ( int scvIdx = 0; scvIdx < numScv; ++scvIdx)
432 scvGeoms_[scvIdx].setCorners(scvCorners[scvIdx], ScvLocalGeometry::numCorners);
437 static ScvLocalGeometry scvGeoms_[numScv];
440 template < class Scalar>
441 typename VcfvScvGeometries<Scalar, 3, Dune::GeometryType::cube>::ScvLocalGeometry
442 VcfvScvGeometries<Scalar, 3, Dune::GeometryType::cube>::scvGeoms_[
443 VcfvScvGeometries<Scalar, 3, Dune::GeometryType::cube>::numScv];
468 template < class Scalar, class Gr idView>
471 enum{dim = GridView::dimension};
472 enum{dimWorld = GridView::dimensionworld};
473 enum{maxNC = (dim < 3 ? 4 : 8)};
474 enum{maxNE = (dim < 3 ? 4 : 12)};
475 enum{maxNF = (dim < 3 ? 1 : 6)};
476 enum{maxBF = (dim < 3 ? 8 : 24)};
477 typedef typename GridView::ctype CoordScalar;
478 typedef typename GridView::Traits::template Codim<0>::Entity Element;
479 typedef typename GridView::Traits::template Codim<0>::EntityPointer ElementPointer;
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 typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> LocalFiniteElementCache;
489 typedef typename LocalFiniteElementCache::FiniteElementType LocalFiniteElement;
490 typedef typename LocalFiniteElement::Traits::LocalBasisType::Traits LocalBasisTraits;
491 typedef typename LocalBasisTraits::JacobianType ShapeJacobian;
493 Scalar quadrilateralArea( const GlobalPosition& p0,
494 const GlobalPosition& p1,
495 const GlobalPosition& p2,
496 const GlobalPosition& p3)
498 return 0.5*std::abs((p3[0] - p1[0])*(p2[1] - p0[1]) - (p3[1] - p1[1])*(p2[0] - p0[0]));
501 void crossProduct(DimVector& c, const DimVector& a, const DimVector& b)
503 c[0] = a[1]*b[2] - a[2]*b[1];
504 c[1] = a[2]*b[0] - a[0]*b[2];
505 c[2] = a[0]*b[1] - a[1]*b[0];
508 Scalar pyramidVolume( const GlobalPosition& p0,
509 const GlobalPosition& p1,
510 const GlobalPosition& p2,
511 const GlobalPosition& p3,
512 const GlobalPosition& p4)
514 DimVector a(p2); a -= p0;
515 DimVector b(p3); b -= p1;
518 crossProduct(n, a, b);
522 return 1.0/6.0*(n*a);
525 Scalar prismVolume( const GlobalPosition& p0,
526 const GlobalPosition& p1,
527 const GlobalPosition& p2,
528 const GlobalPosition& p3,
529 const GlobalPosition& p4,
530 const GlobalPosition& p5)
533 for ( int k = 0; k < dimWorld; ++k)
536 for ( int k = 0; k < dimWorld; ++k)
539 crossProduct(m, a, b);
541 for ( int k = 0; k < dimWorld; ++k)
542 a[k] = p1[k] - p0[k];
543 for ( int k = 0; k < dimWorld; ++k)
544 b[k] = p2[k] - p0[k];
546 crossProduct(n, a, b);
549 for ( int k = 0; k < dimWorld; ++k)
550 a[k] = p5[k] - p0[k];
552 return std::abs(1.0/6.0*(n*a));
555 Scalar hexahedronVolume( const GlobalPosition& p0,
556 const GlobalPosition& p1,
557 const GlobalPosition& p2,
558 const GlobalPosition& p3,
559 const GlobalPosition& p4,
560 const GlobalPosition& p5,
561 const GlobalPosition& p6,
562 const GlobalPosition& p7)
565 prismVolume(p0,p1,p2,p4,p5,p6)
566 + prismVolume(p0,p2,p3,p4,p6,p7);
569 void normalOfQuadrilateral3D(DimVector &normal,
570 const GlobalPosition& p0,
571 const GlobalPosition& p1,
572 const GlobalPosition& p2,
573 const GlobalPosition& p3)
576 for ( int k = 0; k < dimWorld; ++k)
579 for ( int k = 0; k < dimWorld; ++k)
582 crossProduct(normal, a, b);
586 Scalar quadrilateralArea3D( const GlobalPosition& p0,
587 const GlobalPosition& p1,
588 const GlobalPosition& p2,
589 const GlobalPosition& p3)
592 normalOfQuadrilateral3D(normal, p0, p1, p2, p3);
593 return normal.two_norm();
596 void getFaceIndices( int numElemVertices, int k, int& leftFace, int& rightFace)
598 static const int edgeToFaceTet[2][6] = {
602 static const int edgeToFacePyramid[2][8] = {
603 {1, 2, 3, 4, 1, 3, 4, 2},
604 {0, 0, 0, 0, 3, 2, 1, 4}
606 static const int edgeToFacePrism[2][9] = {
607 {1, 0, 2, 0, 1, 2, 4, 4, 4},
608 {0, 2, 1, 3, 3, 3, 0, 1, 2}
610 static const int edgeToFaceHex[2][12] = {
611 {0, 2, 3, 1, 4, 1, 2, 4, 0, 5, 5, 3},
612 {2, 1, 0, 3, 0, 4, 4, 3, 5, 1, 2, 5}
615 switch (numElemVertices) {
617 leftFace = edgeToFaceTet[0][k];
618 rightFace = edgeToFaceTet[1][k];
621 leftFace = edgeToFacePyramid[0][k];
622 rightFace = edgeToFacePyramid[1][k];
625 leftFace = edgeToFacePrism[0][k];
626 rightFace = edgeToFacePrism[1][k];
629 leftFace = edgeToFaceHex[0][k];
630 rightFace = edgeToFaceHex[1][k];
633 OPM_THROW(std::logic_error,
634 "Not implemented: VcfvStencil::getFaceIndices for "
635 << numElemVertices << " vertices");
640 void getEdgeIndices( int numElemVertices, int face, int vert, int& leftEdge, int& rightEdge)
642 static const int faceAndVertexToLeftEdgeTet[4][4] = {
648 static const int faceAndVertexToRightEdgeTet[4][4] = {
654 static const int faceAndVertexToLeftEdgePyramid[5][5] = {
661 static const int faceAndVertexToRightEdgePyramid[5][5] = {
668 static const int faceAndVertexToLeftEdgePrism[5][6] = {
669 { 3, 3, -1, 0, 1, -1},
670 { 4, -1, 4, 0, -1, 2},
671 {-1, 5, 5, -1, 1, 2},
672 { 3, 3, 5, -1, -1, -1},
673 {-1, -1, -1, 6, 6, 8}
675 static const int faceAndVertexToRightEdgePrism[5][6] = {
676 { 0, 1, -1, 6, 6, -1},
677 { 0, -1, 2, 7, -1, 7},
678 {-1, 1, 2, -1, 8, 8},
679 { 4, 5, 4, -1, -1, -1},
680 {-1, -1, -1, 7, 8, 7}
682 static const int faceAndVertexToLeftEdgeHex[6][8] = {
683 { 0, -1, 4, -1, 8, -1, 2, -1},
684 {-1, 5, -1, 3, -1, 1, -1, 9},
685 { 6, 1, -1, -1, 0, 10, -1, -1},
686 {-1, -1, 2, 7, -1, -1, 11, 3},
687 { 4, 6, 7, 5, -1, -1, -1, -1},
688 {-1, -1, -1, -1, 10, 9, 8, 11}
690 static const int faceAndVertexToRightEdgeHex[6][8] = {
691 { 4, -1, 2, -1, 0, -1, 8, -1},
692 {-1, 1, -1, 5, -1, 9, -1, 3},
693 { 0, 6, -1, -1, 10, 1, -1, -1},
694 {-1, -1, 7, 3, -1, -1, 2, 11},
695 { 6, 5, 4, 7, -1, -1, -1, -1},
696 {-1, -1, -1, -1, 8, 10, 11, 9}
699 switch (numElemVertices) {
701 leftEdge = faceAndVertexToLeftEdgeTet[face][vert];
702 rightEdge = faceAndVertexToRightEdgeTet[face][vert];
705 leftEdge = faceAndVertexToLeftEdgePyramid[face][vert];
706 rightEdge = faceAndVertexToRightEdgePyramid[face][vert];
709 leftEdge = faceAndVertexToLeftEdgePrism[face][vert];
710 rightEdge = faceAndVertexToRightEdgePrism[face][vert];
713 leftEdge = faceAndVertexToLeftEdgeHex[face][vert];
714 rightEdge = faceAndVertexToRightEdgeHex[face][vert];
717 OPM_THROW(std::logic_error,
718 "Not implemented: VcfvStencil::getFaceIndices for "
719 << numElemVertices << " vertices");
725 typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView,
734 const GlobalPosition corner( int cornerIdx) const
737 const GlobalPosition global( const LocalPosition &localPos) const
738 { return element_->geometry().global(localPos); }
813 : gridView_(gridView)
814 , vertexMapper_(gridView)
815 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
816 , element_(*gridView.template begin<0>())
818 , elementPtr_(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();
840 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
843 numVertices = e.subEntities(dim);
844 numEdges = e.subEntities(dim-1);
845 numFaces = (dim<3)?0:e.subEntities(1);
847 elementPtr_ = ElementPointer(e);
849 numVertices = e.template count<dim>();
850 numEdges = e.template count<dim-1>();
851 numFaces = (dim<3)?0:e.template count</*codim=*/1>();
854 numBoundarySegments_ = 0;
861 const Geometry& geometry = e.geometry();
862 geometryType_ = geometry.type();
864 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3)
865 const typename Dune::ReferenceElementContainer<CoordScalar,dim>::value_type&
866 referenceElement = Dune::ReferenceElements<CoordScalar,dim>::general(geometryType_);
868 const typename Dune::GenericReferenceElementContainer<CoordScalar,dim>::value_type&
869 referenceElement = Dune::GenericReferenceElements<CoordScalar,
870 dim>::general(geometryType_);
873 elementVolume = geometry.volume();
874 elementLocal = referenceElement.position(0,0);
875 elementGlobal = geometry.global(elementLocal);
878 for ( int vert = 0; vert < numVertices; vert++) {
879 subContVol[vert]. local = referenceElement.position(vert, dim);
880 subContVol[vert]. global = geometry.global(subContVol[vert].local);
884 for ( int edge = 0; edge < numEdges; edge++) {
885 edgeCoord[edge] = geometry.global(referenceElement.position(edge, dim-1));
889 for ( int face = 0; face < numFaces; face++) {
890 faceCoord[face] = geometry.global(referenceElement.position(face, 1));
898 fillSubContVolData_();
901 for ( int k = 0; k < numEdges; k++) {
902 int i = referenceElement.subEntity(k, dim-1, 0, dim);
903 int j = referenceElement.subEntity(k, dim-1, 1, dim);
904 if (numEdges == 4 && (i == 2 || j == 2))
906 subContVolFace[k]. i = i;
907 subContVolFace[k]. j = j;
914 LocalPosition ipLocal_;
918 subContVolFace[k]. normal_ = 1.0;
919 subContVolFace[k]. area_ = 1.0;
920 ipLocal_ = subContVolFace[k]. ipLocal_;
923 ipLocal_ = referenceElement.position(k, dim-1) + elementLocal;
925 subContVolFace[k]. ipLocal_ = ipLocal_;
926 for ( int m = 0; m < dimWorld; ++m)
927 diffVec[m] = elementGlobal[m] - edgeCoord[k][m];
928 subContVolFace[k]. normal_[0] = diffVec[1];
929 subContVolFace[k]. normal_[1] = -diffVec[0];
931 for ( int m = 0; m < dimWorld; ++m)
932 diffVec[m] = subContVol[j].global[m] - subContVol[i].global[m];
934 if (subContVolFace[k].normal_ * diffVec < 0)
935 subContVolFace[k]. normal_ *= -1;
937 subContVolFace[k]. area_ = subContVolFace[k]. normal_.two_norm();
943 getFaceIndices(numVertices, k, leftFace, rightFace);
944 ipLocal_ = referenceElement.position(k, dim-1) + elementLocal
945 + referenceElement.position(leftFace, 1)
946 + referenceElement.position(rightFace, 1);
948 subContVolFace[k]. ipLocal_ = ipLocal_;
949 normalOfQuadrilateral3D(subContVolFace[k].normal_,
950 edgeCoord[k], faceCoord[rightFace],
951 elementGlobal, faceCoord[leftFace]);
952 subContVolFace[k]. area_ = subContVolFace[k]. normal_.two_norm();
957 subContVolFace[k]. ipGlobal_ = geometry.global(ipLocal_);
961 IntersectionIterator endit = gridView_.iend(e);
962 for (IntersectionIterator it = gridView_.ibegin(e); it != endit; ++it) {
963 const typename IntersectionIterator::Intersection& intersection = *it ;
965 if ( ! intersection.boundary())
968 int face = intersection.indexInInside();
969 int numVerticesOfFace = referenceElement.size(face, 1, dim);
970 for ( int vertInFace = 0; vertInFace < numVerticesOfFace; vertInFace++)
972 int vertInElement = referenceElement.subEntity(face, 1, vertInFace, dim);
973 int bfIdx = numBoundarySegments_;
974 ++numBoundarySegments_;
977 boundaryFace_[bfIdx]. ipLocal_ = referenceElement.position(vertInElement, dim);
978 boundaryFace_[bfIdx]. area_ = 1.0;
981 boundaryFace_[bfIdx]. ipLocal_ = referenceElement.position(vertInElement, dim)
982 + referenceElement.position(face, 1);
983 boundaryFace_[bfIdx]. ipLocal_ *= 0.5;
984 boundaryFace_[bfIdx]. area_ = 0.5 * intersection.geometry().volume();
989 getEdgeIndices(numVertices, face, vertInElement, leftEdge, rightEdge);
990 boundaryFace_[bfIdx]. ipLocal_ = referenceElement.position(vertInElement, dim)
991 + referenceElement.position(face, 1)
992 + referenceElement.position(leftEdge, dim-1)
993 + referenceElement.position(rightEdge, dim-1);
994 boundaryFace_[bfIdx]. ipLocal_ *= 0.25;
995 boundaryFace_[bfIdx]. area_ =
996 quadrilateralArea3D(subContVol[vertInElement].global,
997 edgeCoord[rightEdge],
999 edgeCoord[leftEdge]);
1002 OPM_THROW(std::logic_error, "Not implemented:VcfvStencil for dim = " << dim);
1004 boundaryFace_[bfIdx]. ipGlobal_ = geometry.global(boundaryFace_[bfIdx].ipLocal_);
1005 boundaryFace_[bfIdx]. i = vertInElement;
1006 boundaryFace_[bfIdx]. j = vertInElement;
1009 boundaryFace_[bfIdx]. normal_ = intersection.centerUnitOuterNormal();
1018 auto geomType = element.geometry().type();
1021 if (geomType.isTriangle() || geomType.isTetrahedron()) {
1022 for ( int vertIdx = 0; vertIdx < numVertices; ++vertIdx) {
1025 &VcfvScvGeometries<Scalar, dim, Dune::GeometryType::simplex>::get(vertIdx);
1028 else if (geomType.isLine() || geomType.isQuadrilateral() || geomType.isHexahedron()) {
1029 for ( int vertIdx = 0; vertIdx < numVertices; ++vertIdx) {
1032 &VcfvScvGeometries<Scalar, dim, Dune::GeometryType::cube>::get(vertIdx);
1036 OPM_THROW(std::logic_error,
1037 "Not implemented: SCV geometries for non hexahedron elements");
1042 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
1043 const auto &localFiniteElement = feCache_.get(element_.type());
1044 const auto &geom = element_.geometry();
1046 const auto &localFiniteElement = feCache_.get(elementPtr_->type());
1047 const auto &geom = elementPtr_->geometry();
1049 std::vector<ShapeJacobian> localJac;
1051 for ( int scvIdx = 0; scvIdx < numVertices; ++ scvIdx) {
1053 localFiniteElement.localBasis().evaluateJacobian(localCenter, localJac);
1056 const auto jacInvT = geom.jacobianInverseTransposed(globalPos);
1057 for ( int vert = 0; vert < numVertices; vert++) {
1058 jacInvT.mv(localJac[vert][0], subContVol[scvIdx].gradCenter[vert]);
1064 { return numVertices; }
1070 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
1071 { return element_.template subEntity<dim>(scvIdx)-> partitionType(); }
1073 { return elementPtr_->template subEntity<dim>(scvIdx)-> partitionType(); }
1078 assert(0 <= scvIdx && scvIdx < numDof());
1079 return subContVol[scvIdx];
1083 { return numEdges; }
1086 { return numBoundarySegments_; }
1089 { return subContVolFace[faceIdx]; }
1092 { return boundaryFace_[bfIdx]; }
1100 assert(0 <= dofIdx && dofIdx < numDof());
1102 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
1103 return vertexMapper_.subIndex(element_, dofIdx, dim);
1105 return vertexMapper_.map(*elementPtr_, dofIdx, dim);
1110 void fillSubContVolData_()
1114 subContVol[0]. volume_ = 0.5*elementVolume;
1115 subContVol[1]. volume_ = 0.5*elementVolume;
1117 else if (dim == 2) {
1118 switch (numVertices) {
1120 subContVol[0]. volume_ = elementVolume/3;
1121 subContVol[1]. volume_ = elementVolume/3;
1122 subContVol[2]. volume_ = elementVolume/3;
1126 quadrilateralArea(subContVol[0].global,
1131 quadrilateralArea(subContVol[1].global,
1136 quadrilateralArea(subContVol[2].global,
1141 quadrilateralArea(subContVol[3].global,
1147 OPM_THROW(std::logic_error,
1148 "Not implemented:VcfvStencil dim = " << dim
1149 << ", numVertices = " << numVertices);
1152 else if (dim == 3) {
1153 switch (numVertices) {
1155 for ( int k = 0; k < numVertices; k++)
1156 subContVol[k].volume_ = elementVolume / 4.0;
1159 subContVol[0]. volume_ = hexahedronVolume(subContVol[0].global,
1167 subContVol[1]. volume_ = hexahedronVolume(subContVol[1].global,
1175 subContVol[2]. volume_ = hexahedronVolume(subContVol[2].global,
1183 subContVol[3]. volume_ = hexahedronVolume(subContVol[3].global,
1191 subContVol[4]. volume_ = elementVolume -
1198 subContVol[0]. volume_ = hexahedronVolume(subContVol[0].global,
1206 subContVol[1]. volume_ = hexahedronVolume(subContVol[1].global,
1214 subContVol[2]. volume_ = hexahedronVolume(subContVol[2].global,
1222 subContVol[3]. volume_ = hexahedronVolume(edgeCoord[0],
1226 subContVol[3].global,
1230 subContVol[4]. volume_ = hexahedronVolume(edgeCoord[1],
1234 subContVol[4].global,
1238 subContVol[5]. volume_ = hexahedronVolume(edgeCoord[2],
1242 subContVol[5].global,
1248 subContVol[0]. volume_ = hexahedronVolume(subContVol[0].global,
1256 subContVol[1]. volume_ = hexahedronVolume(subContVol[1].global,
1264 subContVol[2]. volume_ = hexahedronVolume(subContVol[2].global,
1272 subContVol[3]. volume_ = hexahedronVolume(subContVol[3].global,
1280 subContVol[4]. volume_ = hexahedronVolume(edgeCoord[0],
1284 subContVol[4].global,
1288 subContVol[5]. volume_ = hexahedronVolume(edgeCoord[1],
1292 subContVol[5].global,
1296 subContVol[6]. volume_ = hexahedronVolume(edgeCoord[2],
1300 subContVol[6].global,
1304 subContVol[7]. volume_ = hexahedronVolume(edgeCoord[3],
1308 subContVol[7].global,
1314 OPM_THROW(std::logic_error,
1315 "Not implemented:VcfvStencil for dim = " << dim
1316 << ", numVertices = " << numVertices);
1320 OPM_THROW(std::logic_error, "Not implemented:VcfvStencil for dim = " << dim);
1324 VertexMapper vertexMapper_;
1326 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
1329 ElementPointer elementPtr_;
1332 static LocalFiniteElementCache feCache_;
1335 LocalPosition elementLocal;
1337 GlobalPosition elementGlobal;
1339 Scalar elementVolume;
1341 SubControlVolume subContVol[maxNC];
1343 SubControlVolumeFace subContVolFace[maxNE];
1345 BoundaryFace boundaryFace_[maxBF];
1346 int numBoundarySegments_;
1348 GlobalPosition edgeCoord[maxNE];
1350 GlobalPosition faceCoord[maxNF];
1357 Dune::GeometryType geometryType_;
1360 template< class Scalar, class Gr idView>
1361 typename VcfvStencil<Scalar, GridView>::LocalFiniteElementCache
1362 VcfvStencil<Scalar, GridView>::feCache_;
int interiorIndex() const Definition: vcfvstencil.hh:782
const ScvLocalGeometry & localGeometry() const Definition: vcfvstencil.hh:740
Scalar area_ area of face Definition: vcfvstencil.hh:806
const SubControlVolume & subControlVolume(int scvIdx) const Definition: vcfvstencil.hh:1076
Quadrature geometry for quadrilaterals.
void updateCenterGradients() Definition: vcfvstencil.hh:1040
Scalar area() const Definition: vcfvstencil.hh:788
const GlobalPosition & center() const Returns the center of weight of the polyhedron. Definition: quadraturegeometries.hh:67
const GlobalPosition global(const LocalPosition &localPos) const Definition: vcfvstencil.hh:737
int numInteriorFaces() const Definition: vcfvstencil.hh:1082
SubControlVolumeFace BoundaryFace compatibility typedef Definition: vcfvstencil.hh:810
int exteriorIndex() const Definition: vcfvstencil.hh:785
GlobalPosition global global vertex position Definition: vcfvstencil.hh:767
LocalPosition local local vertex position Definition: vcfvstencil.hh:765
DimVector normal_ normal on face pointing to CV j or outward of the domain with length equal to |scvf| ... Definition: vcfvstencil.hh:804
int i scvf seperates corner i and j of elem Definition: vcfvstencil.hh:798
const ScvLocalGeometry * localGeometry_ Definition: vcfvstencil.hh:743
finite volume intersected with element Definition: vcfvstencil.hh:747
int globalSpaceIndex(int dofIdx) const Return the global space index given the index of a degree of freedom. Definition: vcfvstencil.hh:1098
int j Definition: vcfvstencil.hh:798
Quadrature geometry for quadrilaterals. Definition: quadraturegeometries.hh:37
Dune::FieldVector< DimVector, maxNC > gradCenter derivative of shape function at the center of the sub control volume Definition: vcfvstencil.hh:774
GlobalPosition ipGlobal_ integration point in global coords Definition: vcfvstencil.hh:802
const GlobalPosition corner(int cornerIdx) const Definition: vcfvstencil.hh:734
Definition: cartesianindexmapper.hh:31
Scalar volume_ space occupied by the sub-control volume Definition: vcfvstencil.hh:769
ScvGeometry geometry_ The geometry of the sub-control volume in local coordinates. Definition: vcfvstencil.hh:771
Dune::PartitionType partitionType(int scvIdx) const Definition: vcfvstencil.hh:1069
interior face of a sub control volume Definition: vcfvstencil.hh:777
void update(const Element &e) Definition: vcfvstencil.hh:857
Represents the finite volume geometry of a single element in the VCFV discretization. Definition: vcfvstencil.hh:469
VcfvStencil(const GridView &gridView) Definition: vcfvstencil.hh:812
const GlobalPosition & integrationPos() const Definition: vcfvstencil.hh:794
int numDof() const Definition: vcfvstencil.hh:1063
Definition: baseauxiliarymodule.hh:35
const DimVector & normal() const Definition: vcfvstencil.hh:779
const Element * element_ Definition: vcfvstencil.hh:744
const ScvLocalGeometry & localGeometry() const Definition: vcfvstencil.hh:758
const GlobalPosition & globalPos() const Definition: vcfvstencil.hh:749
void updateTopology(const Element &e) Update the non-geometric part of the stencil. Definition: vcfvstencil.hh:838
const GlobalPosition center() const Definition: vcfvstencil.hh:731
Dune::MultipleCodimMultipleGeomTypeMapper< GridView, Dune::MCMGVertexLayout > VertexMapper Definition: vcfvstencil.hh:726
int numPrimaryDof() const Definition: vcfvstencil.hh:1066
int numBoundaryFaces() const Definition: vcfvstencil.hh:1085
const ScvGeometry & geometry() const Definition: vcfvstencil.hh:761
const GlobalPosition & corner(int cornerIdx) const Return the position of the corner with a given index. Definition: quadraturegeometries.hh:125
const BoundaryFace & boundaryFace(int bfIdx) const Definition: vcfvstencil.hh:1091
LocalPosition ipLocal_ integration point in local coords Definition: vcfvstencil.hh:800
Scalar volume() const Definition: vcfvstencil.hh:755
const LocalPosition & localPos() const Definition: vcfvstencil.hh:791
const GlobalPosition center() const Definition: vcfvstencil.hh:752
const SubControlVolumeFace & interiorFace(int faceIdx) const Definition: vcfvstencil.hh:1088
Definition: vcfvstencil.hh:728
void updateScvGeometry(const Element &element) Definition: vcfvstencil.hh:1016
|