24#ifndef DUMUX_FVMPFAL3D_INTERACTIONVOLUMECONTAINER_ADAPTIVE_HH
25#define DUMUX_FVMPFAL3D_INTERACTIONVOLUMECONTAINER_ADAPTIVE_HH
45template<
class TypeTag>
54 dim = GridView::dimension, dimWorld = GridView::dimensionworld
60 using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
64 using Element =
typename GridView::Traits::template Codim<0>::Entity;
65 using Vertex =
typename GridView::Traits::template Codim<dim>::Entity;
66 using ElementGeometry =
typename Element::Geometry;
68 using GlobalPosition =
typename ElementGeometry::GlobalCoordinate;
69 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
71 using DimVector = Dune::FieldVector<Scalar, dim>;
72 using FaceVerticesVector = std::vector<Dune::FieldVector<std::set<int>, 2*dim> >;
76 pressureEqIdx = Indices::pressureEqIdx,
79 using IndexTranslator = IndexTranslatorAdaptive;
106 return faceVertices_[eIdxGlobal][fIdx];
118 DUNE_THROW(Dune::NotImplemented,
"Dimension not supported!");
126 Implementation &asImp_()
127 {
return *
static_cast<Implementation *
>(
this);}
130 const Implementation &asImp_()
const
131 {
return *
static_cast<const Implementation *
>(
this);}
134 FaceVerticesVector faceVertices_;
160template<
class TypeTag>
161void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeInnerInteractionVolume(InteractionVolume& interactionVolume,
162 const Vertex& vertex)
167 if (!interactionVolume.sameLevel())
170 std::vector < std::vector<int> > levelIdx(8, std::vector<int>(2));
171 for (
int i = 0; i < 8; i++)
174 levelIdx[i][1] = interactionVolume.getSubVolumeElement(i).level();
177 std::sort(levelIdx.begin(), levelIdx.end(), [](
const auto& a,
const auto& b) { return (a[1]<b[1]); });
184 for (
int i = 0; i < 8; i++)
186 int idx = levelIdx[i][0];
187 auto element = interactionVolume.getSubVolumeElement(idx);
189 const ElementGeometry& geometry = element.geometry();
191 const auto referenceElement = ReferenceElements::general(geometry.type());
197 DimVector edgeCoord(geometry.global(referenceElement.position(9, dim - 1)));
198 interactionVolume.setEdgePosition(edgeCoord, 2);
199 edgeCoord = geometry.global(referenceElement.position(3, dim - 1));
200 interactionVolume.setEdgePosition(edgeCoord, 0);
201 edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
202 interactionVolume.setEdgePosition(edgeCoord, 5);
208 DimVector edgeCoord(geometry.global(referenceElement.position(2, dim - 1)));
209 interactionVolume.setEdgePosition(edgeCoord, 0);
210 edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
211 interactionVolume.setEdgePosition(edgeCoord, 2);
212 edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
213 interactionVolume.setEdgePosition(edgeCoord, 3);
219 DimVector edgeCoord(geometry.global(referenceElement.position(1, dim - 1)));
220 interactionVolume.setEdgePosition(edgeCoord, 0);
221 edgeCoord = geometry.global(referenceElement.position(9, dim - 1));
222 interactionVolume.setEdgePosition(edgeCoord, 4);
223 edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
224 interactionVolume.setEdgePosition(edgeCoord, 5);
230 DimVector edgeCoord(geometry.global(referenceElement.position(0, dim - 1)));
231 interactionVolume.setEdgePosition(edgeCoord, 0);
232 edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
233 interactionVolume.setEdgePosition(edgeCoord, 4);
234 edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
235 interactionVolume.setEdgePosition(edgeCoord, 3);
241 DimVector edgeCoord(geometry.global(referenceElement.position(3, dim - 1)));
242 interactionVolume.setEdgePosition(edgeCoord, 1);
243 edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
244 interactionVolume.setEdgePosition(edgeCoord, 2);
245 edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
246 interactionVolume.setEdgePosition(edgeCoord, 5);
252 DimVector edgeCoord(geometry.global(referenceElement.position(2, dim - 1)));
253 interactionVolume.setEdgePosition(edgeCoord, 1);
254 edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
255 interactionVolume.setEdgePosition(edgeCoord, 2);
256 edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
257 interactionVolume.setEdgePosition(edgeCoord, 3);
263 DimVector edgeCoord(geometry.global(referenceElement.position(1, dim - 1)));
264 interactionVolume.setEdgePosition(edgeCoord, 1);
265 edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
266 interactionVolume.setEdgePosition(edgeCoord, 4);
267 edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
268 interactionVolume.setEdgePosition(edgeCoord, 5);
274 DimVector edgeCoord(geometry.global(referenceElement.position(0, dim - 1)));
275 interactionVolume.setEdgePosition(edgeCoord, 1);
276 edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
277 interactionVolume.setEdgePosition(edgeCoord, 4);
278 edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
279 interactionVolume.setEdgePosition(edgeCoord, 3);
285 ParentType::storeInnerInteractionVolume(interactionVolume, vertex,
false);
289 ParentType::storeInnerInteractionVolume(interactionVolume, vertex,
true);
315template<
class TypeTag>
316void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInteractionVolume(InteractionVolume& interactionVolume,
317 const Vertex& vertex)
319 const DimVector& centerPos = vertex.geometry().center();
321 interactionVolume.setCenterPosition(centerPos);
324 std::vector < std::vector<int> > levelIdx(8, std::vector<int>(2));
325 for (
int i = 0; i < 8; i++)
328 if (interactionVolume.hasSubVolumeElement(i))
329 levelIdx[i][1] = interactionVolume.getSubVolumeElement(i).level();
334 std::sort(levelIdx.begin(), levelIdx.end(), [](
const auto& a,
const auto& b) { return (a[1]<b[1]); });
341 for (
int i = 0; i < 8; i++)
343 if (levelIdx[i][1] < 0)
346 int idx = levelIdx[i][0];
348 auto element = interactionVolume.getSubVolumeElement(idx);
350 const ElementGeometry& geometry = element.geometry();
352 const auto referenceElement = ReferenceElements::general(geometry.type());
358 DimVector edgeCoord(geometry.global(referenceElement.position(9, dim - 1)));
359 interactionVolume.setEdgePosition(edgeCoord, 2);
360 edgeCoord = geometry.global(referenceElement.position(3, dim - 1));
361 interactionVolume.setEdgePosition(edgeCoord, 0);
362 edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
363 interactionVolume.setEdgePosition(edgeCoord, 5);
369 DimVector edgeCoord(geometry.global(referenceElement.position(2, dim - 1)));
370 interactionVolume.setEdgePosition(edgeCoord, 0);
371 edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
372 interactionVolume.setEdgePosition(edgeCoord, 2);
373 edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
374 interactionVolume.setEdgePosition(edgeCoord, 3);
380 DimVector edgeCoord(geometry.global(referenceElement.position(1, dim - 1)));
381 interactionVolume.setEdgePosition(edgeCoord, 0);
382 edgeCoord = geometry.global(referenceElement.position(9, dim - 1));
383 interactionVolume.setEdgePosition(edgeCoord, 4);
384 edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
385 interactionVolume.setEdgePosition(edgeCoord, 5);
391 DimVector edgeCoord(geometry.global(referenceElement.position(0, dim - 1)));
392 interactionVolume.setEdgePosition(edgeCoord, 0);
393 edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
394 interactionVolume.setEdgePosition(edgeCoord, 4);
395 edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
396 interactionVolume.setEdgePosition(edgeCoord, 3);
402 DimVector edgeCoord(geometry.global(referenceElement.position(3, dim - 1)));
403 interactionVolume.setEdgePosition(edgeCoord, 1);
404 edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
405 interactionVolume.setEdgePosition(edgeCoord, 2);
406 edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
407 interactionVolume.setEdgePosition(edgeCoord, 5);
413 DimVector edgeCoord(geometry.global(referenceElement.position(2, dim - 1)));
414 interactionVolume.setEdgePosition(edgeCoord, 1);
415 edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
416 interactionVolume.setEdgePosition(edgeCoord, 2);
417 edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
418 interactionVolume.setEdgePosition(edgeCoord, 3);
424 DimVector edgeCoord(geometry.global(referenceElement.position(1, dim - 1)));
425 interactionVolume.setEdgePosition(edgeCoord, 1);
426 edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
427 interactionVolume.setEdgePosition(edgeCoord, 4);
428 edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
429 interactionVolume.setEdgePosition(edgeCoord, 5);
435 DimVector edgeCoord(geometry.global(referenceElement.position(0, dim - 1)));
436 interactionVolume.setEdgePosition(edgeCoord, 1);
437 edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
438 interactionVolume.setEdgePosition(edgeCoord, 4);
439 edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
440 interactionVolume.setEdgePosition(edgeCoord, 3);
450 switch (interactionVolume.getElementNumber())
455 InteractionVolume hangingNodeVolume(problem_.gridView().grid());
457 std::vector<int> elemIdxOld;
458 for (
int i = 0; i < InteractionVolume::subVolumeTotalNum; i++)
460 if (interactionVolume.hasSubVolumeElement(i))
461 elemIdxOld.push_back(i);
464 int zeroFaceIdx = IndexTranslator::getFaceIndexFromElements(elemIdxOld[0], elemIdxOld[1]);
466 std::vector<int> elemIdxNew(2);
467 elemIdxNew[0] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elemIdxOld[0]);
468 elemIdxNew[1] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elemIdxOld[1]);
471 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[0]),
473 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[1]),
476 for (
int elem = 0; elem < InteractionVolume::subVolumeTotalNum; elem++)
478 for (
int i = 0; i < 3; i++)
480 int faceIdxOld = IndexTranslator::getFaceIndexFromSubVolume(elem, i);
481 int faceIdxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, faceIdxOld);
483 for (
int j = 0; j < 3; j++)
485 int elemNew = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elem);
486 int faceIdxNewTest = IndexTranslator::getFaceIndexFromSubVolume(elemNew, j);
488 if (faceIdxNew == faceIdxNewTest)
490 hangingNodeVolume.setNormal(interactionVolume.getNormal(elem, i),
492 hangingNodeVolume.setIndexOnElement(
493 interactionVolume.getIndexOnElement(elem, i), elemNew, j);
499 for (
int i = 0; i < InteractionVolume::fluxFacesTotalNum; i++)
501 int idxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, i);
502 hangingNodeVolume.setFacePosition(interactionVolume.getFacePosition(i), idxNew);
505 for (
int i = 0; i < InteractionVolume::fluxEdgesTotalNum; i++)
507 int idxNew = IndexTranslator::getNewEdgeIdxFromOldFaceIdxto0(zeroFaceIdx, i);
508 hangingNodeVolume.setEdgePosition(interactionVolume.getEdgePosition(i), idxNew);
511 interactionVolume = hangingNodeVolume;
513 interactionVolume.setHangingNodeType(InteractionVolume::twoSmallCells);
515 interactionVolume.setCenterPosition(centerPos);
517 auto element1 = interactionVolume.getSubVolumeElement(0);
519 for (
const auto& intersection :
intersections(problem_.gridView(), element1))
521 int idxInInside = intersection.indexInInside();
523 if (idxInInside == interactionVolume.getIndexOnElement(0, 2))
525 if (intersection.neighbor())
527 auto outside = intersection.outside();
528 if (element1.level() > outside.level())
530 interactionVolume.setSubVolumeElement(outside, 4);
531 interactionVolume.setSubVolumeElement(outside, 5);
535 else if (idxInInside == interactionVolume.getIndexOnElement(0, 1))
537 if (intersection.neighbor())
539 auto outside = intersection.outside();
540 if (element1.level() > outside.level())
542 interactionVolume.setSubVolumeElement(outside, 2);
543 interactionVolume.setSubVolumeElement(outside, 3);
548 auto element2 = interactionVolume.getSubVolumeElement(1);
549 auto element45 = interactionVolume.getSubVolumeElement(4);
550 auto element23 = interactionVolume.getSubVolumeElement(2);
552 for (
const auto& intersection1 :
intersections(problem_.gridView(), element45))
554 if (intersection1.neighbor())
556 auto element45Outside = intersection1.outside();
558 for (
const auto& intersection2 :
intersections(problem_.gridView(), element23))
560 if (intersection2.neighbor())
562 auto element23Outside = intersection2.outside();
564 if (element45Outside == element23Outside && element45Outside != element1
565 && element45Outside != element2)
567 interactionVolume.setSubVolumeElement(element45Outside, 6);
568 interactionVolume.setSubVolumeElement(element45Outside, 7);
569 DimVector normal = intersection2.centerUnitOuterNormal();
570 interactionVolume.setNormal(normal, 2, 2);
571 interactionVolume.setNormal(normal, 3, 2);
573 interactionVolume.setNormal(normal, 6, 0);
574 interactionVolume.setNormal(normal, 7, 0);
576 GlobalPosition globalPosFace = intersection2.geometry().center();
577 interactionVolume.setFacePosition(globalPosFace, 10);
578 interactionVolume.setFacePosition(globalPosFace, 11);
579 interactionVolume.setEdgePosition(centerPos, 4);
581 Scalar faceArea = intersection2.geometry().volume()/4.0;
583 interactionVolume.setFaceArea(faceArea, 10);
584 interactionVolume.setFaceArea(faceArea, 11);
586 normal = intersection1.centerUnitOuterNormal();
587 interactionVolume.setNormal(normal, 4, 2);
588 interactionVolume.setNormal(normal, 5, 1);
590 interactionVolume.setNormal(normal, 6, 1);
591 interactionVolume.setNormal(normal, 7, 2);
593 globalPosFace = intersection1.geometry().center();
594 interactionVolume.setFacePosition(globalPosFace, 5);
595 interactionVolume.setFacePosition(globalPosFace, 7);
596 interactionVolume.setEdgePosition(centerPos, 1);
598 faceArea = intersection1.geometry().volume()/4.0;
600 interactionVolume.setFaceArea(faceArea, 5);
601 interactionVolume.setFaceArea(faceArea, 7);
608 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
609 DimVector edgeCoord2(interactionVolume.getEdgePosition(1));
610 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
611 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
612 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
613 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
615 DimVector crossProductVector1(0);
616 DimVector crossProductVector2(0);
618 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
619 crossProductVector1 = centerPos - globalPosFace;
620 crossProductVector2 = edgeCoord1 - edgeCoord3;
621 Scalar faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
622 interactionVolume.setFaceArea(faceArea, 0);
624 globalPosFace = interactionVolume.getFacePosition(1);
625 crossProductVector1 = centerPos - globalPosFace;
626 crossProductVector2 = edgeCoord1 - edgeCoord4;
627 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
628 interactionVolume.setFaceArea(faceArea, 1);
630 globalPosFace = interactionVolume.getFacePosition(3);
631 crossProductVector1 = centerPos - globalPosFace;
632 crossProductVector2 = edgeCoord1 - edgeCoord6;
633 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
634 interactionVolume.setFaceArea(faceArea, 3);
636 globalPosFace = interactionVolume.getFacePosition(8);
637 crossProductVector1 = centerPos - globalPosFace;
638 crossProductVector2 = edgeCoord3 - edgeCoord6;
639 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
640 interactionVolume.setFaceArea(faceArea, 8);
642 globalPosFace = interactionVolume.getFacePosition(9);
643 crossProductVector1 = centerPos - globalPosFace;
644 crossProductVector2 = edgeCoord3 - edgeCoord4;
645 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
646 interactionVolume.setFaceArea(faceArea, 9);
653 InteractionVolume hangingNodeVolume(problem_.gridView().grid());
655 std::vector<int> elemIdxOld;
656 for (
int i = 0; i < InteractionVolume::subVolumeTotalNum; i++)
658 if (interactionVolume.hasSubVolumeElement(i))
659 elemIdxOld.push_back(i);
662 std::set<int> zeroFaceIdxVec;
663 for (
int i = 0; i < 4; i++)
665 for (
int j = 0; j < 4; j++)
667 int fIdx = IndexTranslator::getFaceIndexFromElements(elemIdxOld[i], elemIdxOld[j]);
669 zeroFaceIdxVec.insert(fIdx);
673 std::vector<int> elemIdxNew(4);
676 if (zeroFaceIdxVec.size() == 2)
678 hangingNodeVolume.setHangingNodeType(InteractionVolume::fourSmallCellsDiag);
680 if (zeroFaceIdxVec.find(0) == zeroFaceIdxVec.end())
681 zeroFaceIdx = *zeroFaceIdxVec.begin();
683 for (
int i = 0; i < 4; i++)
685 elemIdxNew[i] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elemIdxOld[i]);
686 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[i]),
690 else if (zeroFaceIdxVec.size() == 4)
692 std::set<int>::iterator it = zeroFaceIdxVec.begin();
693 for (; it != zeroFaceIdxVec.end(); ++it)
698 for (
int i = 0; i < 4; i++)
700 elemIdxNew[i] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx,
702 if (elemIdxNew[i] == 4 || elemIdxNew[i] == 5 || elemIdxNew[i] == 6 || elemIdxNew[i] == 7)
710 for (
int i = 0; i < 4; i++)
712 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[i]), elemIdxNew[i]);
721 std::set<int>::iterator it = zeroFaceIdxVec.begin();
722 std::cout<<
"zeroFaceIdxVec = ";
723 for(;it != zeroFaceIdxVec.end(); ++it)
730 for (
int elem = 0; elem < InteractionVolume::subVolumeTotalNum; elem++)
732 for (
int i = 0; i < 3; i++)
734 int faceIdxOld = IndexTranslator::getFaceIndexFromSubVolume(elem, i);
735 int faceIdxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, faceIdxOld);
737 for (
int j = 0; j < 3; j++)
739 int elemNew = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elem);
740 int faceIdxNewTest = IndexTranslator::getFaceIndexFromSubVolume(elemNew, j);
742 if (faceIdxNew == faceIdxNewTest)
744 hangingNodeVolume.setNormal(interactionVolume.getNormal(elem, i),
746 hangingNodeVolume.setIndexOnElement(interactionVolume.getIndexOnElement(elem, i), elemNew, j);
752 for (
int i = 0; i < InteractionVolume::fluxFacesTotalNum; i++)
754 int idxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, i);
755 hangingNodeVolume.setFacePosition(interactionVolume.getFacePosition(i), idxNew);
758 for (
int i = 0; i < InteractionVolume::fluxEdgesTotalNum; i++)
760 int idxNew = IndexTranslator::getNewEdgeIdxFromOldFaceIdxto0(zeroFaceIdx, i);
761 hangingNodeVolume.setEdgePosition(interactionVolume.getEdgePosition(i), idxNew);
764 interactionVolume = hangingNodeVolume;
766 interactionVolume.setCenterPosition(centerPos);
768 if (zeroFaceIdxVec.size() == 4)
770 auto element1 = interactionVolume.getSubVolumeElement(0);
771 auto element4 = interactionVolume.getSubVolumeElement(3);
773 auto outside1 = element1;
774 auto outside4 = element4;
776 for (
const auto& intersection1 :
intersections(problem_.gridView(), element1))
778 if (intersection1.neighbor())
780 if (intersection1.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
782 outside1 = intersection1.outside();
787 for (
const auto& intersection4 :
intersections(problem_.gridView(), element4))
789 if (intersection4.neighbor())
791 if (intersection4.indexInInside() == interactionVolume.getIndexOnElement(3, 2))
793 outside4 = intersection4.outside();
798 if (outside1 != outside4)
800 interactionVolume.setHangingNodeType(InteractionVolume::fourSmallCellsEdge);
802 else if (outside1 == outside4)
804 interactionVolume.setHangingNodeType(InteractionVolume::fourSmallCellsFace);
808 if (interactionVolume.getHangingNodeType() == InteractionVolume::fourSmallCellsFace)
810 auto element = interactionVolume.getSubVolumeElement(0);
812 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
813 DimVector edgeCoord2(interactionVolume.getEdgePosition(1));
814 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
815 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
816 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
817 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
819 DimVector crossProductVector1(0);
820 DimVector crossProductVector2(0);
822 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
823 crossProductVector1 = centerPos - globalPosFace;
824 crossProductVector2 = edgeCoord1 - edgeCoord3;
825 Scalar faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
826 interactionVolume.setFaceArea(faceArea, 0);
828 globalPosFace = interactionVolume.getFacePosition(1);
829 crossProductVector1 = centerPos - globalPosFace;
830 crossProductVector2 = edgeCoord1 - edgeCoord4;
831 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
832 interactionVolume.setFaceArea(faceArea, 1);
834 globalPosFace = interactionVolume.getFacePosition(2);
835 crossProductVector1 = centerPos - globalPosFace;
836 crossProductVector2 = edgeCoord1 - edgeCoord5;
837 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
838 interactionVolume.setFaceArea(faceArea, 2);
840 globalPosFace = interactionVolume.getFacePosition(3);
841 crossProductVector1 = centerPos - globalPosFace;
842 crossProductVector2 = edgeCoord1 - edgeCoord6;
843 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
844 interactionVolume.setFaceArea(faceArea, 3);
846 globalPosFace = interactionVolume.getFacePosition(8);
847 crossProductVector1 = centerPos - globalPosFace;
848 crossProductVector2 = edgeCoord3 - edgeCoord6;
849 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
850 interactionVolume.setFaceArea(faceArea, 8);
852 globalPosFace = interactionVolume.getFacePosition(9);
853 crossProductVector1 = centerPos - globalPosFace;
854 crossProductVector2 = edgeCoord3 - edgeCoord4;
855 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
856 interactionVolume.setFaceArea(faceArea, 9);
858 globalPosFace = interactionVolume.getFacePosition(10);
859 crossProductVector1 = centerPos - globalPosFace;
860 crossProductVector2 = edgeCoord4 - edgeCoord5;
861 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
862 interactionVolume.setFaceArea(faceArea, 10);
864 globalPosFace = interactionVolume.getFacePosition(11);
865 crossProductVector1 = centerPos - globalPosFace;
866 crossProductVector2 = edgeCoord5 - edgeCoord6;
867 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
868 interactionVolume.setFaceArea(faceArea, 11);
870 for (
const auto& intersection :
intersections(problem_.gridView(), element))
872 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
874 if (intersection.neighbor())
876 auto outside = intersection.outside();
877 if (element.level() > outside.level())
879 interactionVolume.setSubVolumeElement(outside, 4);
880 interactionVolume.setSubVolumeElement(outside, 5);
881 interactionVolume.setSubVolumeElement(outside, 6);
882 interactionVolume.setSubVolumeElement(outside, 7);
891 else if (interactionVolume.getHangingNodeType() == InteractionVolume::fourSmallCellsEdge)
893 auto element1 = interactionVolume.getSubVolumeElement(0);
894 auto element2 = interactionVolume.getSubVolumeElement(1);
895 auto element3 = interactionVolume.getSubVolumeElement(2);
896 auto element4 = interactionVolume.getSubVolumeElement(3);
898 for (
const auto& intersection1 :
intersections(problem_.gridView(), element1))
900 if (intersection1.neighbor())
902 if (intersection1.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
904 auto outside = intersection1.outside();
905 if (element1.level() > outside.level())
907 interactionVolume.setSubVolumeElement(outside, 4);
912 for (
const auto& intersection2 :
intersections(problem_.gridView(), element2))
914 if (intersection2.neighbor())
916 if (intersection2.indexInInside() == interactionVolume.getIndexOnElement(1, 2))
918 auto outside = intersection2.outside();
919 if (element2.level() > outside.level())
921 interactionVolume.setSubVolumeElement(outside, 5);
928 for (
const auto& intersection3 :
intersections(problem_.gridView(), element3))
930 if (intersection3.neighbor())
932 if (intersection3.indexInInside() == interactionVolume.getIndexOnElement(2, 2))
934 auto outside = intersection3.outside();
935 if (element3.level() > outside.level())
937 interactionVolume.setSubVolumeElement(outside, 6);
943 for (
const auto& intersection4 :
intersections(problem_.gridView(), element4))
945 if (intersection4.neighbor())
947 if (intersection4.indexInInside() == interactionVolume.getIndexOnElement(3, 2))
949 auto outside = intersection4.outside();
950 if (element4.level() > outside.level())
952 interactionVolume.setSubVolumeElement(outside, 7);
960 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
961 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
962 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
963 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
964 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
966 DimVector crossProductVector1(0);
967 DimVector crossProductVector2(0);
969 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
970 crossProductVector1 = centerPos - globalPosFace;
971 crossProductVector2 = edgeCoord1 - edgeCoord3;
972 Scalar faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
973 interactionVolume.setFaceArea(faceArea, 0);
975 globalPosFace = interactionVolume.getFacePosition(1);
976 crossProductVector1 = centerPos - globalPosFace;
977 crossProductVector2 = edgeCoord1 - edgeCoord4;
978 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
979 interactionVolume.setFaceArea(faceArea, 1);
981 globalPosFace = interactionVolume.getFacePosition(2);
982 crossProductVector1 = centerPos - globalPosFace;
983 crossProductVector2 = edgeCoord1 - edgeCoord5;
984 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
985 interactionVolume.setFaceArea(faceArea, 2);
987 globalPosFace = interactionVolume.getFacePosition(3);
988 crossProductVector1 = centerPos - globalPosFace;
989 crossProductVector2 = edgeCoord1 - edgeCoord6;
990 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
991 interactionVolume.setFaceArea(faceArea, 3);
993 globalPosFace = interactionVolume.getFacePosition(8);
994 crossProductVector1 = centerPos - globalPosFace;
995 crossProductVector2 = edgeCoord3 - edgeCoord6;
996 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
997 interactionVolume.setFaceArea(faceArea, 8);
999 globalPosFace = interactionVolume.getFacePosition(9);
1000 crossProductVector1 = centerPos - globalPosFace;
1001 crossProductVector2 = edgeCoord3 - edgeCoord4;
1002 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1003 interactionVolume.setFaceArea(faceArea, 9);
1005 globalPosFace = interactionVolume.getFacePosition(10);
1006 crossProductVector1 = centerPos - globalPosFace;
1007 crossProductVector2 = edgeCoord4 - edgeCoord5;
1008 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1009 interactionVolume.setFaceArea(faceArea, 10);
1011 globalPosFace = interactionVolume.getFacePosition(11);
1012 crossProductVector1 = centerPos - globalPosFace;
1013 crossProductVector2 = edgeCoord5 - edgeCoord6;
1014 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1015 interactionVolume.setFaceArea(faceArea, 11);
1017 auto element5 = interactionVolume.getSubVolumeElement(4);
1018 auto element6 = interactionVolume.getSubVolumeElement(5);
1019 auto element7 = interactionVolume.getSubVolumeElement(6);
1020 auto element8 = interactionVolume.getSubVolumeElement(7);
1022 if (element5 == element6)
1024 interactionVolume.setFacePosition(element5.geometry().center(), 4);
1025 interactionVolume.setFacePosition(element7.geometry().center(), 6);
1027 for (
const auto& intersection :
intersections(problem_.gridView(), element5))
1029 if (intersection.neighbor())
1031 auto outside = intersection.outside();
1033 if (outside == element7 || outside == element8)
1035 int indexInInside = intersection.indexInInside();
1036 interactionVolume.setIndexOnElement(indexInInside, 4, 2);
1037 interactionVolume.setIndexOnElement(indexInInside, 5, 1);
1038 DimVector normal = intersection.centerUnitOuterNormal();
1039 interactionVolume.setNormal(normal, 4, 2);
1040 interactionVolume.setNormal(normal, 5, 1);
1041 globalPosFace = intersection.geometry().center();
1042 interactionVolume.setFacePosition(globalPosFace, 5);
1043 interactionVolume.setFacePosition(globalPosFace, 7);
1044 int indexInOutside = intersection.indexInOutside();
1045 interactionVolume.setIndexOnElement(indexInOutside, 6, 1);
1046 interactionVolume.setIndexOnElement(indexInOutside, 7, 2);
1048 interactionVolume.setNormal(normal, 6, 1);
1049 interactionVolume.setNormal(normal, 7, 2);
1056 DimVector edgeCoord2(interactionVolume.getFacePosition(7));
1058 globalPosFace = interactionVolume.getFacePosition(5);
1059 crossProductVector1 = centerPos - globalPosFace;
1060 crossProductVector2 = edgeCoord2 - edgeCoord4;
1061 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1062 interactionVolume.setFaceArea(faceArea, 5);
1064 globalPosFace = interactionVolume.getFacePosition(7);
1065 crossProductVector1 = centerPos - globalPosFace;
1066 crossProductVector2 = edgeCoord2 - edgeCoord6;
1067 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1068 interactionVolume.setFaceArea(faceArea, 7);
1071 else if (element5 == element7)
1073 interactionVolume.setFacePosition(element6.geometry().center(), 5);
1074 interactionVolume.setFacePosition(element5.geometry().center(), 7);
1075 interactionVolume.setFacePosition(globalPosFace, 6);
1077 for (
const auto& intersection :
intersections(problem_.gridView(), element5))
1079 if (intersection.neighbor())
1081 auto outside = intersection.outside();
1083 if (outside == element6 || outside == element8)
1085 int indexInInside = intersection.indexInInside();
1086 interactionVolume.setIndexOnElement(indexInInside, 4, 1);
1087 interactionVolume.setIndexOnElement(indexInInside, 6, 2);
1088 DimVector normal = intersection.centerUnitOuterNormal();
1089 interactionVolume.setNormal(normal, 4, 1);
1090 interactionVolume.setNormal(normal, 6, 2);
1091 globalPosFace = intersection.geometry().center();
1092 interactionVolume.setFacePosition(globalPosFace, 4);
1093 interactionVolume.setFacePosition(globalPosFace, 6);
1094 int indexInOutside = intersection.indexInOutside();
1095 interactionVolume.setIndexOnElement(indexInOutside, 5, 2);
1096 interactionVolume.setIndexOnElement(indexInOutside, 7, 1);
1098 interactionVolume.setNormal(normal, 5, 2);
1099 interactionVolume.setNormal(normal, 7, 1);
1105 DimVector edgeCoord2(interactionVolume.getFacePosition(4));
1107 globalPosFace = interactionVolume.getFacePosition(4);
1108 crossProductVector1 = centerPos - globalPosFace;
1109 crossProductVector2 = edgeCoord2 - edgeCoord3;
1110 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1111 interactionVolume.setFaceArea(faceArea, 4);
1113 globalPosFace = interactionVolume.getFacePosition(6);
1114 crossProductVector1 = centerPos - globalPosFace;
1115 crossProductVector2 = edgeCoord2 - edgeCoord5;
1116 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1117 interactionVolume.setFaceArea(faceArea, 6);
1122 const ElementGeometry& geometry = element5.geometry();
1124 const auto referenceElement = ReferenceElements::general(geometry.type());
1126 int oldSubVolumElemIdx = IndexTranslator::getOldElemIdxFromNewFaceIdxto0(zeroFaceIdx, 4);
1127 int oldEdgeIdx = IndexTranslator::getOldEdgeIdxFromNewFaceIdxto0(zeroFaceIdx, 1);
1129 DimVector edgeCoord(0.0);
1130 switch (oldSubVolumElemIdx)
1137 edgeCoord = geometry.global(referenceElement.position(9, dim - 1));
1140 edgeCoord = geometry.global(referenceElement.position(3, dim - 1));
1143 edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
1154 edgeCoord = geometry.global(referenceElement.position(2, dim - 1));
1157 edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
1160 edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
1171 edgeCoord = geometry.global(referenceElement.position(1, dim - 1));
1174 edgeCoord = geometry.global(referenceElement.position(9, dim - 1));
1177 edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
1188 edgeCoord = geometry.global(referenceElement.position(0, dim - 1));
1191 edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
1194 edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
1205 edgeCoord = geometry.global(referenceElement.position(3, dim - 1));
1208 edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
1211 edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
1222 edgeCoord = geometry.global(referenceElement.position(2, dim - 1));
1225 edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
1228 edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
1239 edgeCoord = geometry.global(referenceElement.position(1, dim - 1));
1242 edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
1245 edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
1256 edgeCoord = geometry.global(referenceElement.position(0, dim - 1));
1259 edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
1262 edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
1270 interactionVolume.setEdgePosition(edgeCoord, 1);
1272 else if (interactionVolume.getHangingNodeType() == InteractionVolume::fourSmallCellsDiag)
1274 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
1275 DimVector edgeCoord2(interactionVolume.getEdgePosition(1));
1276 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
1277 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
1278 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
1279 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
1281 DimVector crossProductVector1(0);
1282 DimVector crossProductVector2(0);
1284 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
1285 crossProductVector1 = centerPos - globalPosFace;
1286 crossProductVector2 = edgeCoord1 - edgeCoord3;
1287 Scalar faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1288 interactionVolume.setFaceArea(faceArea, 0);
1290 globalPosFace = interactionVolume.getFacePosition(1);
1291 crossProductVector1 = centerPos - globalPosFace;
1292 crossProductVector2 = edgeCoord1 - edgeCoord4;
1293 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1294 interactionVolume.setFaceArea(faceArea, 1);
1296 globalPosFace = interactionVolume.getFacePosition(3);
1297 crossProductVector1 = centerPos - globalPosFace;
1298 crossProductVector2 = edgeCoord1 - edgeCoord6;
1299 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1300 interactionVolume.setFaceArea(faceArea, 3);
1302 globalPosFace = interactionVolume.getFacePosition(5);
1303 crossProductVector1 = centerPos - globalPosFace;
1304 crossProductVector2 = edgeCoord2 - edgeCoord4;
1305 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1306 interactionVolume.setFaceArea(faceArea, 5);
1308 globalPosFace = interactionVolume.getFacePosition(6);
1309 crossProductVector1 = centerPos - globalPosFace;
1310 crossProductVector2 = edgeCoord2 - edgeCoord5;
1311 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1312 interactionVolume.setFaceArea(faceArea, 6);
1314 globalPosFace = interactionVolume.getFacePosition(7);
1315 crossProductVector1 = centerPos - globalPosFace;
1316 crossProductVector2 = edgeCoord2 - edgeCoord6;
1317 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1318 interactionVolume.setFaceArea(faceArea, 7);
1320 globalPosFace = interactionVolume.getFacePosition(8);
1321 crossProductVector1 = centerPos - globalPosFace;
1322 crossProductVector2 = edgeCoord3 - edgeCoord6;
1323 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1324 interactionVolume.setFaceArea(faceArea, 8);
1326 globalPosFace = interactionVolume.getFacePosition(9);
1327 crossProductVector1 = centerPos - globalPosFace;
1328 crossProductVector2 = edgeCoord3 - edgeCoord4;
1329 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1330 interactionVolume.setFaceArea(faceArea, 9);
1332 globalPosFace = interactionVolume.getFacePosition(10);
1333 crossProductVector1 = centerPos - globalPosFace;
1334 crossProductVector2 = edgeCoord4 - edgeCoord5;
1335 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1336 interactionVolume.setFaceArea(faceArea, 10);
1338 globalPosFace = interactionVolume.getFacePosition(11);
1339 crossProductVector1 = centerPos - globalPosFace;
1340 crossProductVector2 = edgeCoord5 - edgeCoord6;
1341 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1342 interactionVolume.setFaceArea(faceArea, 11);
1344 auto element1 = interactionVolume.getSubVolumeElement(0);
1346 bool hasFaceOne =
false;
1347 bool hasFaceTwo =
false;
1348 for (
const auto& intersection :
intersections(problem_.gridView(), element1))
1350 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(0, 1))
1352 if (intersection.neighbor())
1354 auto outside = intersection.outside();
1355 if (element1.level() > outside.level())
1357 interactionVolume.setSubVolumeElement(outside, 2);
1358 interactionVolume.setSubVolumeElement(outside, 3);
1366 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
1368 if (intersection.neighbor())
1370 auto outside = intersection.outside();
1371 if (element1.level() > outside.level())
1373 interactionVolume.setSubVolumeElement(outside, 4);
1374 interactionVolume.setSubVolumeElement(outside, 5);
1390 InteractionVolume hangingNodeVolume(problem_.gridView().grid());
1392 std::vector<int> elemIdxOld;
1393 for (
int i = 0; i < InteractionVolume::subVolumeTotalNum; i++)
1395 if (interactionVolume.hasSubVolumeElement(i))
1396 elemIdxOld.push_back(i);
1399 std::set<int> zeroFaceIdxVec;
1400 for (
int i = 0; i < 6; i++)
1402 for (
int j = 0; j < 6; j++)
1404 int fIdx = IndexTranslator::getFaceIndexFromElements(elemIdxOld[i], elemIdxOld[j]);
1406 zeroFaceIdxVec.insert(fIdx);
1410 std::vector<int> elemIdxNew(6);
1411 int zeroFaceIdx = 0;
1414 std::set<int>::iterator it = zeroFaceIdxVec.begin();
1415 for (; it != zeroFaceIdxVec.end(); ++it)
1420 for (
int i = 0; i < 6; i++)
1422 elemIdxNew[i] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elemIdxOld[i]);
1423 if (elemIdxNew[i] == 6 || elemIdxNew[i] == 7)
1431 for (
int i = 0; i < 6; i++)
1433 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[i]),
1441 for (
int elem = 0; elem < InteractionVolume::subVolumeTotalNum; elem++)
1443 for (
int i = 0; i < 3; i++)
1445 int faceIdxOld = IndexTranslator::getFaceIndexFromSubVolume(elem, i);
1446 int faceIdxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, faceIdxOld);
1448 for (
int j = 0; j < 3; j++)
1450 int elemNew = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elem);
1451 int faceIdxNewTest = IndexTranslator::getFaceIndexFromSubVolume(elemNew, j);
1453 if (faceIdxNew == faceIdxNewTest)
1455 hangingNodeVolume.setNormal(interactionVolume.getNormal(elem, i),
1457 hangingNodeVolume.setIndexOnElement(
1458 interactionVolume.getIndexOnElement(elem, i), elemNew, j);
1464 for (
int i = 0; i < InteractionVolume::fluxFacesTotalNum; i++)
1466 int idxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, i);
1467 hangingNodeVolume.setFaceArea(interactionVolume.getFaceArea(i), idxNew);
1468 hangingNodeVolume.setFacePosition(interactionVolume.getFacePosition(i), idxNew);
1471 for (
int i = 0; i < InteractionVolume::fluxEdgesTotalNum; i++)
1473 int idxNew = IndexTranslator::getNewEdgeIdxFromOldFaceIdxto0(zeroFaceIdx, i);
1474 hangingNodeVolume.setEdgePosition(interactionVolume.getEdgePosition(i), idxNew);
1477 interactionVolume = hangingNodeVolume;
1479 interactionVolume.setCenterPosition(centerPos);
1481 interactionVolume.setHangingNodeType(InteractionVolume::sixSmallCells);
1483 auto element3 = interactionVolume.getSubVolumeElement(2);
1485 for (
const auto& intersection :
intersections(problem_.gridView(), element3))
1487 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(2, 2))
1489 if (intersection.neighbor())
1491 auto outside = intersection.outside();
1492 if (element3.level() > outside.level())
1494 interactionVolume.setSubVolumeElement(outside, 6);
1495 interactionVolume.setSubVolumeElement(outside, 7);
1503 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
1504 DimVector edgeCoord2(interactionVolume.getEdgePosition(1));
1505 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
1506 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
1507 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
1508 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
1510 DimVector crossProductVector1(0);
1511 DimVector crossProductVector2(0);
1513 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
1514 crossProductVector1 = centerPos - globalPosFace;
1515 crossProductVector2 = edgeCoord1 - edgeCoord3;
1516 Scalar faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1517 interactionVolume.setFaceArea(faceArea, 0);
1519 globalPosFace = interactionVolume.getFacePosition(1);
1520 crossProductVector1 = centerPos - globalPosFace;
1521 crossProductVector2 = edgeCoord1 - edgeCoord4;
1522 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1523 interactionVolume.setFaceArea(faceArea, 1);
1525 globalPosFace = interactionVolume.getFacePosition(2);
1526 crossProductVector1 = centerPos - globalPosFace;
1527 crossProductVector2 = edgeCoord1 - edgeCoord5;
1528 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1529 interactionVolume.setFaceArea(faceArea, 2);
1531 globalPosFace = interactionVolume.getFacePosition(3);
1532 crossProductVector1 = centerPos - globalPosFace;
1533 crossProductVector2 = edgeCoord1 - edgeCoord6;
1534 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1535 interactionVolume.setFaceArea(faceArea, 3);
1537 globalPosFace = interactionVolume.getFacePosition(4);
1538 crossProductVector1 = centerPos - globalPosFace;
1539 crossProductVector2 = edgeCoord2 - edgeCoord3;
1540 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1541 interactionVolume.setFaceArea(faceArea, 4);
1543 globalPosFace = interactionVolume.getFacePosition(5);
1544 crossProductVector1 = centerPos - globalPosFace;
1545 crossProductVector2 = edgeCoord2 - edgeCoord4;
1546 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1547 interactionVolume.setFaceArea(faceArea, 5);
1549 globalPosFace = interactionVolume.getFacePosition(7);
1550 crossProductVector1 = centerPos - globalPosFace;
1551 crossProductVector2 = edgeCoord2 - edgeCoord6;
1552 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1553 interactionVolume.setFaceArea(faceArea, 7);
1555 globalPosFace = interactionVolume.getFacePosition(8);
1556 crossProductVector1 = centerPos - globalPosFace;
1557 crossProductVector2 = edgeCoord3 - edgeCoord6;
1558 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1559 interactionVolume.setFaceArea(faceArea, 8);
1561 globalPosFace = interactionVolume.getFacePosition(9);
1562 crossProductVector1 = centerPos - globalPosFace;
1563 crossProductVector2 = edgeCoord3 - edgeCoord4;
1564 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1565 interactionVolume.setFaceArea(faceArea, 9);
1567 globalPosFace = interactionVolume.getFacePosition(10);
1568 crossProductVector1 = centerPos - globalPosFace;
1569 crossProductVector2 = edgeCoord4 - edgeCoord5;
1570 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1571 interactionVolume.setFaceArea(faceArea, 10);
1573 globalPosFace = interactionVolume.getFacePosition(11);
1574 crossProductVector1 = centerPos - globalPosFace;
1575 crossProductVector2 = edgeCoord5 - edgeCoord6;
1576 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1577 interactionVolume.setFaceArea(faceArea, 11);
1582 interactionVolume.printInteractionVolumeInfo();
1583 DUNE_THROW(Dune::NotImplemented,
"Hanging node shape not implemented");
1590template<
class TypeTag>
1593 std::vector < std::vector<int> > elemVertMap(problem_.gridView().size(dim), std::vector<int>(8, -1));
1596 for (
const auto& element : elements(problem_.gridView()))
1597 asImp_().storeSubVolumeElements(element, elemVertMap);
1599 for (
unsigned int i = 0; i < asImp_().interactionVolumes_.size(); i++)
1600 if (asImp_().interactionVolumes_[i].getElementNumber() == 0)
1601 asImp_().interactionVolumes_[i].printInteractionVolumeInfo();
1604 for (
const auto& element : elements(problem_.gridView()))
1605 asImp_().storeIntersectionInfo(element, elemVertMap);
1607 faceVertices_.clear();
1608 faceVertices_.resize(problem_.gridView().size(0));
1612 for (
const auto& vertex : vertices(problem_.gridView()))
1614 int vIdxGlobal = problem_.variables().index(vertex);
1616 InteractionVolume& interactionVolume = asImp_().interactionVolumes_[vIdxGlobal];
1618 if (interactionVolume.getElementNumber() == 8)
1620 asImp_().storeInnerInteractionVolume(interactionVolume, vertex);
1622 else if (interactionVolume.isBoundaryInteractionVolume())
1624 asImp_().storeBoundaryInteractionVolume(interactionVolume, vertex);
1629 storeHangingNodeInteractionVolume(interactionVolume, vertex);
1632 if (!interactionVolume.isBoundaryInteractionVolume())
1634 auto element1 = interactionVolume.getSubVolumeElement(0);
1635 auto element2 = interactionVolume.getSubVolumeElement(1);
1636 auto element3 = interactionVolume.getSubVolumeElement(2);
1637 auto element4 = interactionVolume.getSubVolumeElement(3);
1638 auto element5 = interactionVolume.getSubVolumeElement(4);
1639 auto element6 = interactionVolume.getSubVolumeElement(5);
1640 auto element7 = interactionVolume.getSubVolumeElement(6);
1641 auto element8 = interactionVolume.getSubVolumeElement(7);
1643 int globalIdx1 = problem_.variables().index(element1);
1644 int globalIdx2 = problem_.variables().index(element2);
1645 int globalIdx3 = problem_.variables().index(element3);
1646 int globalIdx4 = problem_.variables().index(element4);
1647 int globalIdx5 = problem_.variables().index(element5);
1648 int globalIdx6 = problem_.variables().index(element6);
1649 int globalIdx7 = problem_.variables().index(element7);
1650 int globalIdx8 = problem_.variables().index(element8);
1652 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(0, 0), globalIdx1, interactionVolume.getIndexOnElement(0, 0));
1653 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(0, 1), globalIdx1, interactionVolume.getIndexOnElement(0, 1));
1654 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(0, 2), globalIdx1, interactionVolume.getIndexOnElement(0, 2));
1655 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(1, 0), globalIdx2, interactionVolume.getIndexOnElement(1, 0));
1656 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(1, 1), globalIdx2, interactionVolume.getIndexOnElement(1, 1));
1657 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(1, 2), globalIdx2, interactionVolume.getIndexOnElement(1, 2));
1658 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(2, 0), globalIdx3, interactionVolume.getIndexOnElement(2, 0));
1659 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(3, 1), globalIdx4, interactionVolume.getIndexOnElement(3, 1));
1660 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(4, 0), globalIdx5, interactionVolume.getIndexOnElement(4, 0));
1661 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(5, 0), globalIdx6, interactionVolume.getIndexOnElement(5, 0));
1663 faceVertices_[globalIdx1][interactionVolume.getIndexOnElement(0, 0)].insert(vIdxGlobal);
1664 faceVertices_[globalIdx1][interactionVolume.getIndexOnElement(0, 1)].insert(vIdxGlobal);
1665 faceVertices_[globalIdx1][interactionVolume.getIndexOnElement(0, 2)].insert(vIdxGlobal);
1666 faceVertices_[globalIdx2][interactionVolume.getIndexOnElement(1, 0)].insert(vIdxGlobal);
1667 faceVertices_[globalIdx2][interactionVolume.getIndexOnElement(1, 1)].insert(vIdxGlobal);
1668 faceVertices_[globalIdx2][interactionVolume.getIndexOnElement(1, 2)].insert(vIdxGlobal);
1669 faceVertices_[globalIdx3][interactionVolume.getIndexOnElement(2, 0)].insert(vIdxGlobal);
1670 faceVertices_[globalIdx4][interactionVolume.getIndexOnElement(3, 1)].insert(vIdxGlobal);
1671 faceVertices_[globalIdx5][interactionVolume.getIndexOnElement(4, 0)].insert(vIdxGlobal);
1672 faceVertices_[globalIdx6][interactionVolume.getIndexOnElement(5, 0)].insert(vIdxGlobal);
1674 if (interactionVolume.getHangingNodeType() != InteractionVolume::twoSmallCells)
1676 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(2, 1), globalIdx3, interactionVolume.getIndexOnElement(2, 1));
1677 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(2, 2), globalIdx3, interactionVolume.getIndexOnElement(2, 2));
1678 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(3, 0), globalIdx4, interactionVolume.getIndexOnElement(3, 0));
1679 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(3, 2), globalIdx4, interactionVolume.getIndexOnElement(3, 2));
1680 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(4, 1), globalIdx5, interactionVolume.getIndexOnElement(4, 1));
1681 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(4, 2), globalIdx5, interactionVolume.getIndexOnElement(4, 2));
1682 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(5, 1), globalIdx6, interactionVolume.getIndexOnElement(5, 1));
1683 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(5, 2), globalIdx6, interactionVolume.getIndexOnElement(5, 2));
1684 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(6, 0), globalIdx7, interactionVolume.getIndexOnElement(6, 0));
1685 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(6, 1), globalIdx7, interactionVolume.getIndexOnElement(6, 1));
1686 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(6, 2), globalIdx7, interactionVolume.getIndexOnElement(6, 2));
1687 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(7, 0), globalIdx8, interactionVolume.getIndexOnElement(7, 0));
1688 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(7, 1), globalIdx8, interactionVolume.getIndexOnElement(7, 1));
1689 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(7, 2), globalIdx8, interactionVolume.getIndexOnElement(7, 2));
1691 faceVertices_[globalIdx3][interactionVolume.getIndexOnElement(2, 1)].insert(vIdxGlobal);
1692 faceVertices_[globalIdx3][interactionVolume.getIndexOnElement(2, 2)].insert(vIdxGlobal);
1693 faceVertices_[globalIdx4][interactionVolume.getIndexOnElement(3, 0)].insert(vIdxGlobal);
1694 faceVertices_[globalIdx4][interactionVolume.getIndexOnElement(3, 2)].insert(vIdxGlobal);
1695 faceVertices_[globalIdx5][interactionVolume.getIndexOnElement(4, 1)].insert(vIdxGlobal);
1696 faceVertices_[globalIdx5][interactionVolume.getIndexOnElement(4, 2)].insert(vIdxGlobal);
1697 faceVertices_[globalIdx6][interactionVolume.getIndexOnElement(5, 1)].insert(vIdxGlobal);
1698 faceVertices_[globalIdx6][interactionVolume.getIndexOnElement(5, 2)].insert(vIdxGlobal);
1699 faceVertices_[globalIdx7][interactionVolume.getIndexOnElement(6, 0)].insert(vIdxGlobal);
1700 faceVertices_[globalIdx7][interactionVolume.getIndexOnElement(6, 1)].insert(vIdxGlobal);
1701 faceVertices_[globalIdx7][interactionVolume.getIndexOnElement(6, 2)].insert(vIdxGlobal);
1702 faceVertices_[globalIdx8][interactionVolume.getIndexOnElement(7, 0)].insert(vIdxGlobal);
1703 faceVertices_[globalIdx8][interactionVolume.getIndexOnElement(7, 1)].insert(vIdxGlobal);
1704 faceVertices_[globalIdx8][interactionVolume.getIndexOnElement(7, 2)].insert(vIdxGlobal);
#define GET_PROP_TYPE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:283
Interactionvolume container for 3-d MPFA L-method.
Class including the information of an interaction volume of a MPFA 3D method that does not change wit...
Dune::FieldVector< Scalar, 3 > crossProduct(const Dune::FieldVector< Scalar, 3 > &vec1, const Dune::FieldVector< Scalar, 3 > &vec2)
Cross product of two vectors in three-dimensional Euclidean space.
Definition: math.hh:631
Dune::IteratorRange< typename MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper >::Intersections::const_iterator > intersections(const MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper > &glue)
Range generator to iterate with range-based for loops over all intersections as follows: for (const a...
Definition: glue.hh:62
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
Property tag Indices
Definition: porousmediumflow/sequential/properties.hh:59
Property tag MPFAInteractionVolumeContainer
Type of the data container for one interaction volume.
Definition: porousmediumflow/sequential/cellcentered/mpfa/properties.hh:97
Property tag MPFAInteractionVolume
Type of the data container for one interaction volume.
Definition: porousmediumflow/sequential/cellcentered/mpfa/properties.hh:96
Interactionvolume container for 3-d MPFA L-method.
Definition: 3dinteractionvolumecontainer.hh:46
typename GET_PROP_TYPE(TypeTag, MPFAInteractionVolume) InteractionVolume
Type for storing an MPFA-interaction-volume. (Usually of type FvMpfaL3dInteractionVolume or FvMpfaL3d...
Definition: 3dinteractionvolumecontainer.hh:93
InteractionVolume & interactionVolume(int vertexIdx)
Returns an interaction volume.
Definition: 3dinteractionvolumecontainer.hh:142
Interactionvolume container for 3-d MPFA L-method on an h-adaptive grid.
Definition: 3dinteractionvolumecontaineradaptive.hh:47
void storeInteractionVolumeInfo()
Stores interaction volumes for each grid vertex.
Definition: 3dinteractionvolumecontaineradaptive.hh:1591
std::set< int > & faceVerticeIndices(int eIdxGlobal, int fIdx)
Returns the set of vertices on an element face.
Definition: 3dinteractionvolumecontaineradaptive.hh:104
FvMpfaL3dInteractionVolumeContainerAdaptive(Problem &problem)
Constructs a FvMpfaL3dInteractionVolumeContainerAdaptive object.
Definition: 3dinteractionvolumecontaineradaptive.hh:113
Properties for a MPFA method.
Base file for properties related to sequential IMPET algorithms.