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 if constexpr (dim == 3) {
843 numFaces =
e.subEntities(1);
849 numBoundarySegments_ = 0;
852 const Geometry& geometry =
e.geometry();
853 geometryType_ = geometry.type();
854 const auto&
referenceElement = Dune::ReferenceElements<CoordScalar,dim>::general(geometryType_);
861 void updatePrimaryTopology(
const Element& element)
869 void update(
const Element&
e)
873 const Geometry& geometry =
e.geometry();
874 geometryType_ = geometry.type();
876 const auto&
referenceElement = Dune::ReferenceElements<CoordScalar,dim>::general(geometryType_);
878 elementVolume = geometry.volume();
880 elementGlobal = geometry.global(elementLocal);
885 subContVol[
vert].
global = geometry.global(subContVol[
vert].local);
903 fillSubContVolData_();
906 for (
unsigned k = 0;
k < numEdges;
k++) {
907 unsigned short i =
static_cast<unsigned short>(
referenceElement.subEntity(
static_cast<int>(
k), dim-1, 0, dim));
908 unsigned short j =
static_cast<unsigned short>(
referenceElement.subEntity(
static_cast<int>(
k), dim-1, 1, dim));
909 if (numEdges == 4 && (i == 2 || j == 2))
911 subContVolFace[
k].
i = i;
912 subContVolFace[
k].j = j;
919 LocalPosition ipLocal_;
921 if constexpr (dim == 1) {
924 subContVolFace[
k].
area_ = 1.0;
927 else if constexpr (dim == 2) {
931 for (
unsigned m = 0;
m < dimWorld; ++
m)
936 for (
unsigned m = 0;
m < dimWorld; ++
m)
937 diffVec[
m] = subContVol[j].global[
m] - subContVol[i].global[
m];
939 if (subContVolFace[
k].normal_ *
diffVec < 0)
945 else if constexpr (dim == 3) {
954 normalOfQuadrilateral3D(subContVolFace[
k].normal_,
956 elementGlobal, faceCoord[
leftFace]);
962 subContVolFace[
k].
ipGlobal_ = geometry.global(ipLocal_);
966 IntersectionIterator
endit = gridView_.iend(
e);
967 for (IntersectionIterator it = gridView_.ibegin(
e); it !=
endit; ++it) {
968 const typename IntersectionIterator::Intersection& intersection = *it ;
970 if ( ! intersection.boundary())
973 unsigned face =
static_cast<unsigned>(intersection.indexInInside());
978 unsigned bfIdx = numBoundarySegments_;
979 ++numBoundarySegments_;
981 if constexpr (dim == 1) {
985 else if constexpr (dim == 2) {
989 boundaryFace_[
bfIdx].
area_ = 0.5 * intersection.geometry().volume();
991 else if constexpr (dim == 3) {
1007 throw std::logic_error(
"Not implemented:VcfvStencil for dim = "+std::to_string(dim));
1014 boundaryFace_[
bfIdx].
normal_ = intersection.centerUnitOuterNormal();
1018 updateScvGeometry(
e);
1021 void updateScvGeometry(
const Element& element)
1023 auto geomType = element.geometry().type();
1030 &VcfvScvGeometries<Scalar, dim, ElementType::simplex>::get(
vertIdx);
1037 &VcfvScvGeometries<Scalar, dim, ElementType::cube>::get(
vertIdx);
1041 throw std::logic_error(
"Not implemented: SCV geometries for non hexahedron elements");
1044#if HAVE_DUNE_LOCALFUNCTIONS
1045 void updateCenterGradients()
1048 const auto&
geom = element_.geometry();
1050 std::vector<ShapeJacobian>
localJac;
1055 const auto& globalPos = subContVol[
scvIdx].geometry().center();
1057 const auto jacInvT =
geom.jacobianInverseTransposed(globalPos);
1065 unsigned numDof()
const
1066 {
return numVertices; }
1068 unsigned numPrimaryDof()
const
1069 {
return numDof(); }
1071 Dune::PartitionType partitionType(
unsigned scvIdx)
const
1074 const SubControlVolume& subControlVolume(
unsigned scvIdx)
const
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(dofIdx < numDof());
1100 return static_cast<unsigned>(vertexMapper_.subIndex(element_,
static_cast<int>(dofIdx), dim));
1109 assert(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_()
1121 if constexpr (dim == 1) {
1123 subContVol[0].
volume_ = 0.5*elementVolume;
1124 subContVol[1].
volume_ = 0.5*elementVolume;
1126 else if constexpr (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 throw std::logic_error(
"Not implemented:VcfvStencil dim = "+std::to_string(dim)
1157 +
", numVertices = "+std::to_string(numVertices));
1160 else if constexpr (dim == 3) {
1161 switch (numVertices) {
1163 for (
unsigned k = 0;
k < numVertices;
k++)
1164 subContVol[
k].volume_ = elementVolume / 4.0;
1168 hexahedronVolume(subContVol[0].global,
1177 hexahedronVolume(subContVol[1].global,
1186 hexahedronVolume(subContVol[2].global,
1195 hexahedronVolume(subContVol[3].global,
1203 subContVol[4].
volume_ = elementVolume -
1211 hexahedronVolume(subContVol[0].global,
1220 hexahedronVolume(subContVol[1].global,
1229 hexahedronVolume(subContVol[2].global,
1238 hexahedronVolume(edgeCoord[0],
1242 subContVol[3].global,
1247 hexahedronVolume(edgeCoord[1],
1251 subContVol[4].global,
1256 hexahedronVolume(edgeCoord[2],
1260 subContVol[5].global,
1267 hexahedronVolume(subContVol[0].global,
1276 hexahedronVolume(subContVol[1].global,
1285 hexahedronVolume(subContVol[2].global,
1294 hexahedronVolume(subContVol[3].global,
1303 hexahedronVolume(edgeCoord[0],
1307 subContVol[4].global,
1312 hexahedronVolume(edgeCoord[1],
1316 subContVol[5].global,
1321 hexahedronVolume(edgeCoord[2],
1325 subContVol[6].global,
1330 hexahedronVolume(edgeCoord[3],
1334 subContVol[7].global,
1340 throw std::logic_error(
"Not implemented:VcfvStencil for dim = "+std::to_string(dim)
1341 +
", numVertices = "+std::to_string(numVertices));
1345 throw std::logic_error(
"Not implemented:VcfvStencil for dim = "+std::to_string(dim));
1347#if __GNUC__ || __clang__
1348#pragma GCC diagnostic pop
1351 const GridView& gridView_;
1352 const Mapper& vertexMapper_;
1356#if HAVE_DUNE_LOCALFUNCTIONS
1361 LocalPosition elementLocal;
1363 GlobalPosition elementGlobal;
1365 Scalar elementVolume;
1367 SubControlVolume subContVol[maxNC];
1369 SubControlVolumeFace subContVolFace[maxNE];
1372 unsigned numBoundarySegments_;
1374 GlobalPosition edgeCoord[maxNE];
1376 GlobalPosition faceCoord[maxNF];
1378 unsigned numVertices;
1383 Dune::GeometryType geometryType_;