468 enum{dim = GridView::dimension};
469 enum{dimWorld = GridView::dimensionworld};
470 enum{maxNC = (dim < 3 ? 4 : 8)};
471 enum{maxNE = (dim < 3 ? 4 : 12)};
472 enum{maxNF = (dim < 3 ? 1 : 6)};
473 enum{maxBF = (dim < 3 ? 8 : 24)};
474 using CoordScalar =
typename GridView::ctype;
475 using Element =
typename GridView::Traits::template Codim<0>::Entity ;
477 using Entity =
typename GridView::Traits::template Codim<dim>::Entity ;
479 using Geometry =
typename Element::Geometry;
480 using DimVector = Dune::FieldVector<Scalar,dimWorld>;
481 using GlobalPosition = Dune::FieldVector<CoordScalar,dimWorld>;
482 using LocalPosition = Dune::FieldVector<CoordScalar,dim>;
483 using IntersectionIterator =
typename GridView::IntersectionIterator;
487#if HAVE_DUNE_LOCALFUNCTIONS
488 using LocalFiniteElementCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
489 using LocalFiniteElement =
typename LocalFiniteElementCache::FiniteElementType;
490 using LocalBasisTraits =
typename LocalFiniteElement::Traits::LocalBasisType::Traits;
491 using ShapeJacobian =
typename LocalBasisTraits::JacobianType;
494 Scalar quadrilateralArea(
const GlobalPosition& p0,
495 const GlobalPosition& p1,
496 const GlobalPosition& p2,
497 const GlobalPosition& p3)
499 return 0.5*std::abs((p3[0] - p1[0])*(p2[1] - p0[1]) - (p3[1] - p1[1])*(p2[0] - p0[0]));
502 void crossProduct(DimVector& c,
const DimVector& a,
const DimVector& b)
504 c[0] = a[1]*b[2] - a[2]*b[1];
505 c[1] = a[2]*b[0] - a[0]*b[2];
506 c[2] = a[0]*b[1] - a[1]*b[0];
509 Scalar pyramidVolume(
const GlobalPosition& p0,
510 const GlobalPosition& p1,
511 const GlobalPosition& p2,
512 const GlobalPosition& p3,
513 const GlobalPosition& p4)
515 DimVector a(p2); a -= p0;
516 DimVector b(p3); b -= p1;
519 crossProduct(n, a, b);
523 return 1.0/6.0*(n*a);
526 Scalar prismVolume(
const GlobalPosition& p0,
527 const GlobalPosition& p1,
528 const GlobalPosition& p2,
529 const GlobalPosition& p3,
530 const GlobalPosition& p4,
531 const GlobalPosition& p5)
534 for (
unsigned k = 0; k < dimWorld; ++k)
537 for (
unsigned k = 0; k < dimWorld; ++k)
540 crossProduct(m, a, b);
542 for (
unsigned k = 0; k < dimWorld; ++k)
543 a[k] = p1[k] - p0[k];
544 for (
unsigned k = 0; k < dimWorld; ++k)
545 b[k] = p2[k] - p0[k];
547 crossProduct(n, a, b);
550 for (
unsigned k = 0; k < dimWorld; ++k)
551 a[k] = p5[k] - p0[k];
553 return std::abs(1.0/6.0*(n*a));
556 Scalar hexahedronVolume(
const GlobalPosition& p0,
557 const GlobalPosition& p1,
558 const GlobalPosition& p2,
559 const GlobalPosition& p3,
560 const GlobalPosition& p4,
561 const GlobalPosition& p5,
562 const GlobalPosition& p6,
563 const GlobalPosition& p7)
566 prismVolume(p0,p1,p2,p4,p5,p6)
567 + prismVolume(p0,p2,p3,p4,p6,p7);
570 void normalOfQuadrilateral3D(DimVector& normal,
571 const GlobalPosition& p0,
572 const GlobalPosition& p1,
573 const GlobalPosition& p2,
574 const GlobalPosition& p3)
577 for (
unsigned k = 0; k < dimWorld; ++k)
580 for (
unsigned k = 0; k < dimWorld; ++k)
583 crossProduct(normal, a, b);
587 Scalar quadrilateralArea3D(
const GlobalPosition& p0,
588 const GlobalPosition& p1,
589 const GlobalPosition& p2,
590 const GlobalPosition& p3)
593 normalOfQuadrilateral3D(normal, p0, p1, p2, p3);
594 return normal.two_norm();
597 void getFaceIndices(
unsigned numElemVertices,
unsigned k,
unsigned& leftFace,
unsigned& rightFace)
599 static const unsigned edgeToFaceTet[2][6] = {
603 static const unsigned edgeToFacePyramid[2][8] = {
604 {1, 2, 3, 4, 1, 3, 4, 2},
605 {0, 0, 0, 0, 3, 2, 1, 4}
607 static const unsigned edgeToFacePrism[2][9] = {
608 {1, 0, 2, 0, 1, 2, 4, 4, 4},
609 {0, 2, 1, 3, 3, 3, 0, 1, 2}
611 static const unsigned edgeToFaceHex[2][12] = {
612 {0, 2, 3, 1, 4, 1, 2, 4, 0, 5, 5, 3},
613 {2, 1, 0, 3, 0, 4, 4, 3, 5, 1, 2, 5}
616 switch (numElemVertices) {
618 leftFace = edgeToFaceTet[0][k];
619 rightFace = edgeToFaceTet[1][k];
622 leftFace = edgeToFacePyramid[0][k];
623 rightFace = edgeToFacePyramid[1][k];
626 leftFace = edgeToFacePrism[0][k];
627 rightFace = edgeToFacePrism[1][k];
630 leftFace = edgeToFaceHex[0][k];
631 rightFace = edgeToFaceHex[1][k];
634 throw std::logic_error(
"Not implemented: VcfvStencil::getFaceIndices for "
635 +std::to_string(numElemVertices)+
" vertices");
640 void getEdgeIndices(
unsigned numElemVertices,
unsigned face,
unsigned vert,
unsigned& leftEdge,
unsigned& 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 =
static_cast<unsigned>(faceAndVertexToLeftEdgeTet[face][vert]);
702 rightEdge =
static_cast<unsigned>(faceAndVertexToRightEdgeTet[face][vert]);
705 leftEdge =
static_cast<unsigned>(faceAndVertexToLeftEdgePyramid[face][vert]);
706 rightEdge =
static_cast<unsigned>(faceAndVertexToRightEdgePyramid[face][vert]);
709 leftEdge =
static_cast<unsigned>(faceAndVertexToLeftEdgePrism[face][vert]);
710 rightEdge =
static_cast<unsigned>(faceAndVertexToRightEdgePrism[face][vert]);
713 leftEdge =
static_cast<unsigned>(faceAndVertexToLeftEdgeHex[face][vert]);
714 rightEdge =
static_cast<unsigned>(faceAndVertexToRightEdgeHex[face][vert]);
717 throw std::logic_error(
"Not implemented: VcfvStencil::getFaceIndices for "
718 +std::to_string(numElemVertices)+
" vertices");
725 using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
730 const GlobalPosition center()
const
731 {
return global(localGeometry_->
center()); }
733 const GlobalPosition corner(
unsigned cornerIdx)
const
734 {
return global(localGeometry_->
corner(cornerIdx)); }
736 const GlobalPosition global(
const LocalPosition& localPos)
const
737 {
return element_->geometry().global(localPos); }
740 {
return *localGeometry_; }
743 const Element *element_;
748 const GlobalPosition& globalPos()
const
751 const GlobalPosition center()
const
754 Scalar volume()
const
778 const DimVector& normal()
const
781 unsigned short interiorIndex()
const
784 unsigned short exteriorIndex()
const
790 const LocalPosition& localPos()
const
793 const GlobalPosition& integrationPos()
const
812 : gridView_(gridView)
813 , vertexMapper_(mapper )
814 , element_(*gridView.template begin<0>())
817 assert(
static_cast<int>(gridView.size(dimWorld)) ==
static_cast<int>(mapper.size()));
819 static bool localGeometriesInitialized =
false;
820 if (!localGeometriesInitialized) {
821 localGeometriesInitialized =
true;
823 VcfvScvGeometries<Scalar, 1, ElementType::cube>::init();
824 VcfvScvGeometries<Scalar, 2, ElementType::cube>::init();
825 VcfvScvGeometries<Scalar, 2, ElementType::simplex>::init();
826 VcfvScvGeometries<Scalar, 3, ElementType::cube>::init();
827 VcfvScvGeometries<Scalar, 3, ElementType::simplex>::init();
840 numVertices = e.subEntities(dim);
841 numEdges = e.subEntities(dim-1);
842 numFaces = (dim<3)?0:e.subEntities(1);
844 numBoundarySegments_ = 0;
847 const Geometry& geometry = e.geometry();
848 geometryType_ = geometry.type();
849 const auto& referenceElement = Dune::ReferenceElements<CoordScalar,dim>::general(geometryType_);
850 for (
unsigned vertexIdx = 0; vertexIdx < numVertices; vertexIdx++) {
851 subContVol[vertexIdx].
local = referenceElement.position(
static_cast<int>(vertexIdx), dim);
852 subContVol[vertexIdx].
global = geometry.corner(
static_cast<int>(vertexIdx));
856 void updatePrimaryTopology(
const Element& element)
864 void update(
const Element& e)
868 const Geometry& geometry = e.geometry();
869 geometryType_ = geometry.type();
871 const auto& referenceElement = Dune::ReferenceElements<CoordScalar,dim>::general(geometryType_);
873 elementVolume = geometry.volume();
874 elementLocal = referenceElement.position(0,0);
875 elementGlobal = geometry.global(elementLocal);
878 for (
unsigned vert = 0; vert < numVertices; vert++) {
879 subContVol[vert].
local = referenceElement.position(
static_cast<int>(vert), dim);
880 subContVol[vert].
global = geometry.global(subContVol[vert].local);
884 for (
unsigned edge = 0; edge < numEdges; edge++) {
885 edgeCoord[edge] = geometry.global(referenceElement.position(
static_cast<int>(edge), dim-1));
889 for (
unsigned face = 0; face < numFaces; face++) {
890 faceCoord[face] = geometry.global(referenceElement.position(
static_cast<int>(face), 1));
898 fillSubContVolData_();
901 for (
unsigned k = 0; k < numEdges; k++) {
902 unsigned short i =
static_cast<unsigned short>(referenceElement.subEntity(
static_cast<int>(k), dim-1, 0, dim));
903 unsigned short j =
static_cast<unsigned short>(referenceElement.subEntity(
static_cast<int>(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(
static_cast<int>(k), dim-1) + elementLocal;
925 subContVolFace[k].
ipLocal_ = ipLocal_;
926 for (
unsigned 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 (
unsigned 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(
static_cast<int>(k), dim-1) + elementLocal
945 + referenceElement.position(
static_cast<int>(leftFace), 1)
946 + referenceElement.position(
static_cast<int>(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 unsigned face =
static_cast<unsigned>(intersection.indexInInside());
969 unsigned numVerticesOfFace =
static_cast<unsigned>(referenceElement.size(
static_cast<int>(face), 1, dim));
970 for (
unsigned vertInFace = 0; vertInFace < numVerticesOfFace; vertInFace++)
972 unsigned short vertInElement =
static_cast<unsigned short>(referenceElement.subEntity(
static_cast<int>(face), 1,
static_cast<int>(vertInFace), dim));
973 unsigned bfIdx = numBoundarySegments_;
974 ++numBoundarySegments_;
977 boundaryFace_[bfIdx].
ipLocal_ = referenceElement.position(
static_cast<int>(vertInElement), dim);
978 boundaryFace_[bfIdx].
area_ = 1.0;
981 boundaryFace_[bfIdx].
ipLocal_ = referenceElement.position(
static_cast<int>(vertInElement), dim)
982 + referenceElement.position(
static_cast<int>(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(
static_cast<int>(vertInElement), dim)
991 + referenceElement.position(
static_cast<int>(face), 1)
992 + referenceElement.position(
static_cast<int>(leftEdge), dim-1)
993 + referenceElement.position(
static_cast<int>(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 throw std::logic_error(
"Not implemented:VcfvStencil for dim = "+std::to_string(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();
1013 updateScvGeometry(e);
1016 void updateScvGeometry(
const Element& element)
1018 auto geomType = element.geometry().type();
1021 if (geomType.isTriangle() || geomType.isTetrahedron()) {
1022 for (
unsigned vertIdx = 0; vertIdx < numVertices; ++vertIdx) {
1023 subContVol[vertIdx].
geometry_.element_ = &element;
1024 subContVol[vertIdx].
geometry_.localGeometry_ =
1025 &VcfvScvGeometries<Scalar, dim, ElementType::simplex>::get(vertIdx);
1028 else if (geomType.isLine() || geomType.isQuadrilateral() || geomType.isHexahedron()) {
1029 for (
unsigned vertIdx = 0; vertIdx < numVertices; ++vertIdx) {
1030 subContVol[vertIdx].
geometry_.element_ = &element;
1031 subContVol[vertIdx].
geometry_.localGeometry_ =
1032 &VcfvScvGeometries<Scalar, dim, ElementType::cube>::get(vertIdx);
1036 throw std::logic_error(
"Not implemented: SCV geometries for non hexahedron elements");
1039#if HAVE_DUNE_LOCALFUNCTIONS
1040 void updateCenterGradients()
1042 const auto& localFiniteElement = feCache_.get(element_.type());
1043 const auto& geom = element_.geometry();
1045 std::vector<ShapeJacobian> localJac;
1047 for (
unsigned scvIdx = 0; scvIdx < numVertices; ++ scvIdx) {
1048 const auto& localCenter = subContVol[scvIdx].localGeometry().
center();
1049 localFiniteElement.localBasis().evaluateJacobian(localCenter, localJac);
1050 const auto& globalPos = subContVol[scvIdx].geometry().center();
1052 const auto jacInvT = geom.jacobianInverseTransposed(globalPos);
1053 for (
unsigned vert = 0; vert < numVertices; vert++) {
1054 jacInvT.mv(localJac[vert][0], subContVol[scvIdx].gradCenter[vert]);
1060 unsigned numDof()
const
1061 {
return numVertices; }
1063 unsigned numPrimaryDof()
const
1064 {
return numDof(); }
1066 Dune::PartitionType partitionType(
unsigned scvIdx)
const
1067 {
return element_.template subEntity<dim>(scvIdx)->partitionType(); }
1069 const SubControlVolume& subControlVolume(
unsigned scvIdx)
const
1071 assert(scvIdx < numDof());
1072 return subContVol[scvIdx];
1075 unsigned numInteriorFaces()
const
1076 {
return numEdges; }
1078 unsigned numBoundaryFaces()
const
1079 {
return numBoundarySegments_; }
1081 const SubControlVolumeFace& interiorFace(
unsigned faceIdx)
const
1082 {
return subContVolFace[faceIdx]; }
1085 {
return boundaryFace_[bfIdx]; }
1093 assert(dofIdx < numDof());
1095 return static_cast<unsigned>(vertexMapper_.subIndex(element_,
static_cast<int>(dofIdx), dim));
1104 assert(dofIdx < numDof());
1105 return element_.template subEntity<dim>(
static_cast<int>(dofIdx));
1109#if __GNUC__ || __clang__
1110#pragma GCC diagnostic push
1111#pragma GCC diagnostic ignored "-Wpragmas"
1112#pragma GCC diagnostic ignored "-Warray-bounds"
1114 void fillSubContVolData_()
1118 subContVol[0].
volume_ = 0.5*elementVolume;
1119 subContVol[1].
volume_ = 0.5*elementVolume;
1121 else if (dim == 2) {
1122 switch (numVertices) {
1124 subContVol[0].
volume_ = elementVolume/3;
1125 subContVol[1].
volume_ = elementVolume/3;
1126 subContVol[2].
volume_ = elementVolume/3;
1130 quadrilateralArea(subContVol[0].global,
1135 quadrilateralArea(subContVol[1].global,
1140 quadrilateralArea(subContVol[2].global,
1145 quadrilateralArea(subContVol[3].global,
1151 throw std::logic_error(
"Not implemented:VcfvStencil dim = "+std::to_string(dim)
1152 +
", numVertices = "+std::to_string(numVertices));
1155 else if (dim == 3) {
1156 switch (numVertices) {
1158 for (
unsigned k = 0; k < numVertices; k++)
1159 subContVol[k].volume_ = elementVolume / 4.0;
1163 hexahedronVolume(subContVol[0].global,
1172 hexahedronVolume(subContVol[1].global,
1181 hexahedronVolume(subContVol[2].global,
1190 hexahedronVolume(subContVol[3].global,
1198 subContVol[4].
volume_ = elementVolume -
1206 hexahedronVolume(subContVol[0].global,
1215 hexahedronVolume(subContVol[1].global,
1224 hexahedronVolume(subContVol[2].global,
1233 hexahedronVolume(edgeCoord[0],
1237 subContVol[3].global,
1242 hexahedronVolume(edgeCoord[1],
1246 subContVol[4].global,
1251 hexahedronVolume(edgeCoord[2],
1255 subContVol[5].global,
1262 hexahedronVolume(subContVol[0].global,
1271 hexahedronVolume(subContVol[1].global,
1280 hexahedronVolume(subContVol[2].global,
1289 hexahedronVolume(subContVol[3].global,
1298 hexahedronVolume(edgeCoord[0],
1302 subContVol[4].global,
1307 hexahedronVolume(edgeCoord[1],
1311 subContVol[5].global,
1316 hexahedronVolume(edgeCoord[2],
1320 subContVol[6].global,
1325 hexahedronVolume(edgeCoord[3],
1329 subContVol[7].global,
1335 throw std::logic_error(
"Not implemented:VcfvStencil for dim = "+std::to_string(dim)
1336 +
", numVertices = "+std::to_string(numVertices));
1340 throw std::logic_error(
"Not implemented:VcfvStencil for dim = "+std::to_string(dim));
1342#if __GNUC__ || __clang__
1343#pragma GCC diagnostic pop
1346 const GridView& gridView_;
1347 const Mapper& vertexMapper_;
1351#if HAVE_DUNE_LOCALFUNCTIONS
1352 static LocalFiniteElementCache feCache_;
1356 LocalPosition elementLocal;
1358 GlobalPosition elementGlobal;
1360 Scalar elementVolume;
1362 SubControlVolume subContVol[maxNC];
1364 SubControlVolumeFace subContVolFace[maxNE];
1367 unsigned numBoundarySegments_;
1369 GlobalPosition edgeCoord[maxNE];
1371 GlobalPosition faceCoord[maxNF];
1373 unsigned numVertices;
1378 Dune::GeometryType geometryType_;