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
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)
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)
544 for (
unsigned k = 0;
k < dimWorld; ++
k)
547 crossProduct(
n,
a,
b);
550 for (
unsigned k = 0;
k < dimWorld; ++
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)
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();
604 {1, 2, 3, 4, 1, 3, 4, 2},
605 {0, 0, 0, 0, 3, 2, 1, 4}
608 {1, 0, 2, 0, 1, 2, 4, 4, 4},
609 {0, 2, 1, 3, 3, 3, 0, 1, 2}
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}
634 throw std::logic_error(
"Not implemented: VcfvStencil::getFaceIndices for "
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}
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}
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}
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}
717 throw std::logic_error(
"Not implemented: VcfvStencil::getFaceIndices for "
725 using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
730 const GlobalPosition center()
const
731 {
return global(localGeometry_->
center()); }
733 const GlobalPosition corner(
unsigned cornerIdx)
const
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)
814 , element_(*gridView.
template begin<0>())
817 assert(
static_cast<int>(gridView.size(dimWorld)) ==
static_cast<int>(
mapper.size()));
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_);
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();
875 elementGlobal = geometry.global(elementLocal);
880 subContVol[
vert].
global = geometry.global(subContVol[
vert].local);
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_;
919 subContVolFace[
k].
area_ = 1.0;
926 for (
unsigned m = 0;
m < dimWorld; ++
m)
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)
949 normalOfQuadrilateral3D(subContVolFace[
k].normal_,
951 elementGlobal, faceCoord[
leftFace]);
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());
973 unsigned bfIdx = numBoundarySegments_;
974 ++numBoundarySegments_;
984 boundaryFace_[
bfIdx].
area_ = 0.5 * intersection.geometry().volume();
1002 throw std::logic_error(
"Not implemented:VcfvStencil for dim = "+std::to_string(dim));
1009 boundaryFace_[
bfIdx].
normal_ = intersection.centerUnitOuterNormal();
1013 updateScvGeometry(
e);
1016 void updateScvGeometry(
const Element& element)
1018 auto geomType = element.geometry().type();
1025 &VcfvScvGeometries<Scalar, dim, ElementType::simplex>::get(
vertIdx);
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()
1043 const auto&
geom = element_.geometry();
1045 std::vector<ShapeJacobian>
localJac;
1050 const auto& globalPos = subContVol[
scvIdx].geometry().center();
1052 const auto jacInvT =
geom.jacobianInverseTransposed(globalPos);
1060 unsigned numDof()
const
1061 {
return numVertices; }
1063 unsigned numPrimaryDof()
const
1064 {
return numDof(); }
1066 Dune::PartitionType partitionType(
unsigned scvIdx)
const
1069 const SubControlVolume& subControlVolume(
unsigned scvIdx)
const
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]; }
1095 return static_cast<unsigned>(vertexMapper_.subIndex(element_,
static_cast<int>(
dofIdx), dim));
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
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_;