24#ifndef DUMUX_FVMPFAL3D_INTERACTIONVOLUMECONTAINER_ADAPTIVE_HH
25#define DUMUX_FVMPFAL3D_INTERACTIONVOLUMECONTAINER_ADAPTIVE_HH
45template<
class TypeTag>
54 dim = GridView::dimension, dimWorld = GridView::dimensionworld
62 using Element =
typename GridView::Traits::template Codim<0>::Entity;
63 using Vertex =
typename GridView::Traits::template Codim<dim>::Entity;
64 using ElementGeometry =
typename Element::Geometry;
66 using GlobalPosition =
typename ElementGeometry::GlobalCoordinate;
67 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
69 using DimVector = Dune::FieldVector<Scalar, dim>;
70 using FaceVerticesVector = std::vector<Dune::FieldVector<std::set<int>, 2*dim> >;
74 pressureEqIdx = Indices::pressureEqIdx,
77 using IndexTranslator = IndexTranslatorAdaptive;
104 return faceVertices_[eIdxGlobal][fIdx];
116 DUNE_THROW(Dune::NotImplemented,
"Dimension not supported!");
124 Implementation &asImp_()
125 {
return *
static_cast<Implementation *
>(
this);}
128 const Implementation &asImp_()
const
129 {
return *
static_cast<const Implementation *
>(
this);}
132 FaceVerticesVector faceVertices_;
158template<
class TypeTag>
159void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeInnerInteractionVolume(InteractionVolume& interactionVolume,
160 const Vertex& vertex)
165 if (!interactionVolume.sameLevel())
168 std::vector < std::vector<int> > levelIdx(8, std::vector<int>(2));
169 for (
int i = 0; i < 8; i++)
172 levelIdx[i][1] = interactionVolume.getSubVolumeElement(i).level();
175 std::sort(levelIdx.begin(), levelIdx.end(), [](
const auto& a,
const auto& b) { return (a[1]<b[1]); });
182 for (
int i = 0; i < 8; i++)
184 int idx = levelIdx[i][0];
185 auto element = interactionVolume.getSubVolumeElement(idx);
187 const ElementGeometry& geometry = element.geometry();
189 const auto refElement = referenceElement(geometry);
195 DimVector edgeCoord(geometry.global(refElement.position(9, dim - 1)));
196 interactionVolume.setEdgePosition(edgeCoord, 2);
197 edgeCoord = geometry.global(refElement.position(3, dim - 1));
198 interactionVolume.setEdgePosition(edgeCoord, 0);
199 edgeCoord = geometry.global(refElement.position(11, dim - 1));
200 interactionVolume.setEdgePosition(edgeCoord, 5);
206 DimVector edgeCoord(geometry.global(refElement.position(2, dim - 1)));
207 interactionVolume.setEdgePosition(edgeCoord, 0);
208 edgeCoord = geometry.global(refElement.position(8, dim - 1));
209 interactionVolume.setEdgePosition(edgeCoord, 2);
210 edgeCoord = geometry.global(refElement.position(11, dim - 1));
211 interactionVolume.setEdgePosition(edgeCoord, 3);
217 DimVector edgeCoord(geometry.global(refElement.position(1, dim - 1)));
218 interactionVolume.setEdgePosition(edgeCoord, 0);
219 edgeCoord = geometry.global(refElement.position(9, dim - 1));
220 interactionVolume.setEdgePosition(edgeCoord, 4);
221 edgeCoord = geometry.global(refElement.position(10, dim - 1));
222 interactionVolume.setEdgePosition(edgeCoord, 5);
228 DimVector edgeCoord(geometry.global(refElement.position(0, dim - 1)));
229 interactionVolume.setEdgePosition(edgeCoord, 0);
230 edgeCoord = geometry.global(refElement.position(8, dim - 1));
231 interactionVolume.setEdgePosition(edgeCoord, 4);
232 edgeCoord = geometry.global(refElement.position(10, dim - 1));
233 interactionVolume.setEdgePosition(edgeCoord, 3);
239 DimVector edgeCoord(geometry.global(refElement.position(3, dim - 1)));
240 interactionVolume.setEdgePosition(edgeCoord, 1);
241 edgeCoord = geometry.global(refElement.position(5, dim - 1));
242 interactionVolume.setEdgePosition(edgeCoord, 2);
243 edgeCoord = geometry.global(refElement.position(7, dim - 1));
244 interactionVolume.setEdgePosition(edgeCoord, 5);
250 DimVector edgeCoord(geometry.global(refElement.position(2, dim - 1)));
251 interactionVolume.setEdgePosition(edgeCoord, 1);
252 edgeCoord = geometry.global(refElement.position(4, dim - 1));
253 interactionVolume.setEdgePosition(edgeCoord, 2);
254 edgeCoord = geometry.global(refElement.position(7, dim - 1));
255 interactionVolume.setEdgePosition(edgeCoord, 3);
261 DimVector edgeCoord(geometry.global(refElement.position(1, dim - 1)));
262 interactionVolume.setEdgePosition(edgeCoord, 1);
263 edgeCoord = geometry.global(refElement.position(5, dim - 1));
264 interactionVolume.setEdgePosition(edgeCoord, 4);
265 edgeCoord = geometry.global(refElement.position(6, dim - 1));
266 interactionVolume.setEdgePosition(edgeCoord, 5);
272 DimVector edgeCoord(geometry.global(refElement.position(0, dim - 1)));
273 interactionVolume.setEdgePosition(edgeCoord, 1);
274 edgeCoord = geometry.global(refElement.position(4, dim - 1));
275 interactionVolume.setEdgePosition(edgeCoord, 4);
276 edgeCoord = geometry.global(refElement.position(6, dim - 1));
277 interactionVolume.setEdgePosition(edgeCoord, 3);
283 ParentType::storeInnerInteractionVolume(interactionVolume, vertex,
false);
287 ParentType::storeInnerInteractionVolume(interactionVolume, vertex,
true);
313template<
class TypeTag>
314void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInteractionVolume(InteractionVolume& interactionVolume,
315 const Vertex& vertex)
317 const DimVector& centerPos = vertex.geometry().center();
319 interactionVolume.setCenterPosition(centerPos);
322 std::vector < std::vector<int> > levelIdx(8, std::vector<int>(2));
323 for (
int i = 0; i < 8; i++)
326 if (interactionVolume.hasSubVolumeElement(i))
327 levelIdx[i][1] = interactionVolume.getSubVolumeElement(i).level();
332 std::sort(levelIdx.begin(), levelIdx.end(), [](
const auto& a,
const auto& b) { return (a[1]<b[1]); });
339 for (
int i = 0; i < 8; i++)
341 if (levelIdx[i][1] < 0)
344 int idx = levelIdx[i][0];
346 auto element = interactionVolume.getSubVolumeElement(idx);
348 const ElementGeometry& geometry = element.geometry();
350 const auto refElement = referenceElement(geometry);
356 DimVector edgeCoord(geometry.global(refElement.position(9, dim - 1)));
357 interactionVolume.setEdgePosition(edgeCoord, 2);
358 edgeCoord = geometry.global(refElement.position(3, dim - 1));
359 interactionVolume.setEdgePosition(edgeCoord, 0);
360 edgeCoord = geometry.global(refElement.position(11, dim - 1));
361 interactionVolume.setEdgePosition(edgeCoord, 5);
367 DimVector edgeCoord(geometry.global(refElement.position(2, dim - 1)));
368 interactionVolume.setEdgePosition(edgeCoord, 0);
369 edgeCoord = geometry.global(refElement.position(8, dim - 1));
370 interactionVolume.setEdgePosition(edgeCoord, 2);
371 edgeCoord = geometry.global(refElement.position(11, dim - 1));
372 interactionVolume.setEdgePosition(edgeCoord, 3);
378 DimVector edgeCoord(geometry.global(refElement.position(1, dim - 1)));
379 interactionVolume.setEdgePosition(edgeCoord, 0);
380 edgeCoord = geometry.global(refElement.position(9, dim - 1));
381 interactionVolume.setEdgePosition(edgeCoord, 4);
382 edgeCoord = geometry.global(refElement.position(10, dim - 1));
383 interactionVolume.setEdgePosition(edgeCoord, 5);
389 DimVector edgeCoord(geometry.global(refElement.position(0, dim - 1)));
390 interactionVolume.setEdgePosition(edgeCoord, 0);
391 edgeCoord = geometry.global(refElement.position(8, dim - 1));
392 interactionVolume.setEdgePosition(edgeCoord, 4);
393 edgeCoord = geometry.global(refElement.position(10, dim - 1));
394 interactionVolume.setEdgePosition(edgeCoord, 3);
400 DimVector edgeCoord(geometry.global(refElement.position(3, dim - 1)));
401 interactionVolume.setEdgePosition(edgeCoord, 1);
402 edgeCoord = geometry.global(refElement.position(5, dim - 1));
403 interactionVolume.setEdgePosition(edgeCoord, 2);
404 edgeCoord = geometry.global(refElement.position(7, dim - 1));
405 interactionVolume.setEdgePosition(edgeCoord, 5);
411 DimVector edgeCoord(geometry.global(refElement.position(2, dim - 1)));
412 interactionVolume.setEdgePosition(edgeCoord, 1);
413 edgeCoord = geometry.global(refElement.position(4, dim - 1));
414 interactionVolume.setEdgePosition(edgeCoord, 2);
415 edgeCoord = geometry.global(refElement.position(7, dim - 1));
416 interactionVolume.setEdgePosition(edgeCoord, 3);
422 DimVector edgeCoord(geometry.global(refElement.position(1, dim - 1)));
423 interactionVolume.setEdgePosition(edgeCoord, 1);
424 edgeCoord = geometry.global(refElement.position(5, dim - 1));
425 interactionVolume.setEdgePosition(edgeCoord, 4);
426 edgeCoord = geometry.global(refElement.position(6, dim - 1));
427 interactionVolume.setEdgePosition(edgeCoord, 5);
433 DimVector edgeCoord(geometry.global(refElement.position(0, dim - 1)));
434 interactionVolume.setEdgePosition(edgeCoord, 1);
435 edgeCoord = geometry.global(refElement.position(4, dim - 1));
436 interactionVolume.setEdgePosition(edgeCoord, 4);
437 edgeCoord = geometry.global(refElement.position(6, dim - 1));
438 interactionVolume.setEdgePosition(edgeCoord, 3);
448 switch (interactionVolume.getElementNumber())
453 InteractionVolume hangingNodeVolume(problem_.gridView().grid());
455 std::vector<int> elemIdxOld;
456 for (
int i = 0; i < InteractionVolume::subVolumeTotalNum; i++)
458 if (interactionVolume.hasSubVolumeElement(i))
459 elemIdxOld.push_back(i);
462 int zeroFaceIdx = IndexTranslator::getFaceIndexFromElements(elemIdxOld[0], elemIdxOld[1]);
464 std::vector<int> elemIdxNew(2);
465 elemIdxNew[0] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elemIdxOld[0]);
466 elemIdxNew[1] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elemIdxOld[1]);
469 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[0]),
471 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[1]),
474 for (
int elem = 0; elem < InteractionVolume::subVolumeTotalNum; elem++)
476 for (
int i = 0; i < 3; i++)
478 int faceIdxOld = IndexTranslator::getFaceIndexFromSubVolume(elem, i);
479 int faceIdxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, faceIdxOld);
481 for (
int j = 0; j < 3; j++)
483 int elemNew = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elem);
484 int faceIdxNewTest = IndexTranslator::getFaceIndexFromSubVolume(elemNew, j);
486 if (faceIdxNew == faceIdxNewTest)
488 hangingNodeVolume.setNormal(interactionVolume.getNormal(elem, i),
490 hangingNodeVolume.setIndexOnElement(
491 interactionVolume.getIndexOnElement(elem, i), elemNew, j);
497 for (
int i = 0; i < InteractionVolume::fluxFacesTotalNum; i++)
499 int idxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, i);
500 hangingNodeVolume.setFacePosition(interactionVolume.getFacePosition(i), idxNew);
503 for (
int i = 0; i < InteractionVolume::fluxEdgesTotalNum; i++)
505 int idxNew = IndexTranslator::getNewEdgeIdxFromOldFaceIdxto0(zeroFaceIdx, i);
506 hangingNodeVolume.setEdgePosition(interactionVolume.getEdgePosition(i), idxNew);
509 interactionVolume = hangingNodeVolume;
511 interactionVolume.setHangingNodeType(InteractionVolume::twoSmallCells);
513 interactionVolume.setCenterPosition(centerPos);
515 auto element1 = interactionVolume.getSubVolumeElement(0);
517 for (
const auto& intersection : intersections(problem_.gridView(), element1))
519 int idxInInside = intersection.indexInInside();
521 if (idxInInside == interactionVolume.getIndexOnElement(0, 2))
523 if (intersection.neighbor())
525 auto outside = intersection.outside();
526 if (element1.level() > outside.level())
528 interactionVolume.setSubVolumeElement(outside, 4);
529 interactionVolume.setSubVolumeElement(outside, 5);
533 else if (idxInInside == interactionVolume.getIndexOnElement(0, 1))
535 if (intersection.neighbor())
537 auto outside = intersection.outside();
538 if (element1.level() > outside.level())
540 interactionVolume.setSubVolumeElement(outside, 2);
541 interactionVolume.setSubVolumeElement(outside, 3);
546 auto element2 = interactionVolume.getSubVolumeElement(1);
547 auto element45 = interactionVolume.getSubVolumeElement(4);
548 auto element23 = interactionVolume.getSubVolumeElement(2);
550 for (
const auto& intersection1 : intersections(problem_.gridView(), element45))
552 if (intersection1.neighbor())
554 auto element45Outside = intersection1.outside();
556 for (
const auto& intersection2 : intersections(problem_.gridView(), element23))
558 if (intersection2.neighbor())
560 auto element23Outside = intersection2.outside();
562 if (element45Outside == element23Outside && element45Outside != element1
563 && element45Outside != element2)
565 interactionVolume.setSubVolumeElement(element45Outside, 6);
566 interactionVolume.setSubVolumeElement(element45Outside, 7);
567 DimVector
normal = intersection2.centerUnitOuterNormal();
568 interactionVolume.setNormal(
normal, 2, 2);
569 interactionVolume.setNormal(
normal, 3, 2);
571 interactionVolume.setNormal(
normal, 6, 0);
572 interactionVolume.setNormal(
normal, 7, 0);
574 GlobalPosition globalPosFace = intersection2.geometry().center();
575 interactionVolume.setFacePosition(globalPosFace, 10);
576 interactionVolume.setFacePosition(globalPosFace, 11);
577 interactionVolume.setEdgePosition(centerPos, 4);
579 Scalar faceArea = intersection2.geometry().volume()/4.0;
581 interactionVolume.setFaceArea(faceArea, 10);
582 interactionVolume.setFaceArea(faceArea, 11);
584 normal = intersection1.centerUnitOuterNormal();
585 interactionVolume.setNormal(
normal, 4, 2);
586 interactionVolume.setNormal(
normal, 5, 1);
588 interactionVolume.setNormal(
normal, 6, 1);
589 interactionVolume.setNormal(
normal, 7, 2);
591 globalPosFace = intersection1.geometry().center();
592 interactionVolume.setFacePosition(globalPosFace, 5);
593 interactionVolume.setFacePosition(globalPosFace, 7);
594 interactionVolume.setEdgePosition(centerPos, 1);
596 faceArea = intersection1.geometry().volume()/4.0;
598 interactionVolume.setFaceArea(faceArea, 5);
599 interactionVolume.setFaceArea(faceArea, 7);
606 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
607 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
608 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
609 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
611 DimVector crossProductVector1(0);
612 DimVector crossProductVector2(0);
614 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
615 crossProductVector1 = centerPos - globalPosFace;
616 crossProductVector2 = edgeCoord1 - edgeCoord3;
617 Scalar faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
618 interactionVolume.setFaceArea(faceArea, 0);
620 globalPosFace = interactionVolume.getFacePosition(1);
621 crossProductVector1 = centerPos - globalPosFace;
622 crossProductVector2 = edgeCoord1 - edgeCoord4;
623 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
624 interactionVolume.setFaceArea(faceArea, 1);
626 globalPosFace = interactionVolume.getFacePosition(3);
627 crossProductVector1 = centerPos - globalPosFace;
628 crossProductVector2 = edgeCoord1 - edgeCoord6;
629 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
630 interactionVolume.setFaceArea(faceArea, 3);
632 globalPosFace = interactionVolume.getFacePosition(8);
633 crossProductVector1 = centerPos - globalPosFace;
634 crossProductVector2 = edgeCoord3 - edgeCoord6;
635 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
636 interactionVolume.setFaceArea(faceArea, 8);
638 globalPosFace = interactionVolume.getFacePosition(9);
639 crossProductVector1 = centerPos - globalPosFace;
640 crossProductVector2 = edgeCoord3 - edgeCoord4;
641 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
642 interactionVolume.setFaceArea(faceArea, 9);
649 InteractionVolume hangingNodeVolume(problem_.gridView().grid());
651 std::vector<int> elemIdxOld;
652 for (
int i = 0; i < InteractionVolume::subVolumeTotalNum; i++)
654 if (interactionVolume.hasSubVolumeElement(i))
655 elemIdxOld.push_back(i);
658 std::set<int> zeroFaceIdxVec;
659 for (
int i = 0; i < 4; i++)
661 for (
int j = 0; j < 4; j++)
663 int fIdx = IndexTranslator::getFaceIndexFromElements(elemIdxOld[i], elemIdxOld[j]);
665 zeroFaceIdxVec.insert(fIdx);
669 std::vector<int> elemIdxNew(4);
672 if (zeroFaceIdxVec.size() == 2)
674 hangingNodeVolume.setHangingNodeType(InteractionVolume::fourSmallCellsDiag);
676 if (zeroFaceIdxVec.find(0) == zeroFaceIdxVec.end())
677 zeroFaceIdx = *zeroFaceIdxVec.begin();
679 for (
int i = 0; i < 4; i++)
681 elemIdxNew[i] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elemIdxOld[i]);
682 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[i]),
686 else if (zeroFaceIdxVec.size() == 4)
688 std::set<int>::iterator it = zeroFaceIdxVec.begin();
689 for (; it != zeroFaceIdxVec.end(); ++it)
694 for (
int i = 0; i < 4; i++)
696 elemIdxNew[i] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx,
698 if (elemIdxNew[i] == 4 || elemIdxNew[i] == 5 || elemIdxNew[i] == 6 || elemIdxNew[i] == 7)
706 for (
int i = 0; i < 4; i++)
708 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[i]), elemIdxNew[i]);
717 std::set<int>::iterator it = zeroFaceIdxVec.begin();
718 std::cout<<
"zeroFaceIdxVec = ";
719 for(;it != zeroFaceIdxVec.end(); ++it)
726 for (
int elem = 0; elem < InteractionVolume::subVolumeTotalNum; elem++)
728 for (
int i = 0; i < 3; i++)
730 int faceIdxOld = IndexTranslator::getFaceIndexFromSubVolume(elem, i);
731 int faceIdxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, faceIdxOld);
733 for (
int j = 0; j < 3; j++)
735 int elemNew = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elem);
736 int faceIdxNewTest = IndexTranslator::getFaceIndexFromSubVolume(elemNew, j);
738 if (faceIdxNew == faceIdxNewTest)
740 hangingNodeVolume.setNormal(interactionVolume.getNormal(elem, i),
742 hangingNodeVolume.setIndexOnElement(interactionVolume.getIndexOnElement(elem, i), elemNew, j);
748 for (
int i = 0; i < InteractionVolume::fluxFacesTotalNum; i++)
750 int idxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, i);
751 hangingNodeVolume.setFacePosition(interactionVolume.getFacePosition(i), idxNew);
754 for (
int i = 0; i < InteractionVolume::fluxEdgesTotalNum; i++)
756 int idxNew = IndexTranslator::getNewEdgeIdxFromOldFaceIdxto0(zeroFaceIdx, i);
757 hangingNodeVolume.setEdgePosition(interactionVolume.getEdgePosition(i), idxNew);
760 interactionVolume = hangingNodeVolume;
762 interactionVolume.setCenterPosition(centerPos);
764 if (zeroFaceIdxVec.size() == 4)
766 auto element1 = interactionVolume.getSubVolumeElement(0);
767 auto element4 = interactionVolume.getSubVolumeElement(3);
769 auto outside1 = element1;
770 auto outside4 = element4;
772 for (
const auto& intersection1 : intersections(problem_.gridView(), element1))
774 if (intersection1.neighbor())
776 if (intersection1.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
778 outside1 = intersection1.outside();
783 for (
const auto& intersection4 : intersections(problem_.gridView(), element4))
785 if (intersection4.neighbor())
787 if (intersection4.indexInInside() == interactionVolume.getIndexOnElement(3, 2))
789 outside4 = intersection4.outside();
794 if (outside1 != outside4)
796 interactionVolume.setHangingNodeType(InteractionVolume::fourSmallCellsEdge);
798 else if (outside1 == outside4)
800 interactionVolume.setHangingNodeType(InteractionVolume::fourSmallCellsFace);
804 if (interactionVolume.getHangingNodeType() == InteractionVolume::fourSmallCellsFace)
806 auto element = interactionVolume.getSubVolumeElement(0);
808 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
809 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
810 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
811 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
812 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
814 DimVector crossProductVector1(0);
815 DimVector crossProductVector2(0);
817 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
818 crossProductVector1 = centerPos - globalPosFace;
819 crossProductVector2 = edgeCoord1 - edgeCoord3;
820 Scalar faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
821 interactionVolume.setFaceArea(faceArea, 0);
823 globalPosFace = interactionVolume.getFacePosition(1);
824 crossProductVector1 = centerPos - globalPosFace;
825 crossProductVector2 = edgeCoord1 - edgeCoord4;
826 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
827 interactionVolume.setFaceArea(faceArea, 1);
829 globalPosFace = interactionVolume.getFacePosition(2);
830 crossProductVector1 = centerPos - globalPosFace;
831 crossProductVector2 = edgeCoord1 - edgeCoord5;
832 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
833 interactionVolume.setFaceArea(faceArea, 2);
835 globalPosFace = interactionVolume.getFacePosition(3);
836 crossProductVector1 = centerPos - globalPosFace;
837 crossProductVector2 = edgeCoord1 - edgeCoord6;
838 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
839 interactionVolume.setFaceArea(faceArea, 3);
841 globalPosFace = interactionVolume.getFacePosition(8);
842 crossProductVector1 = centerPos - globalPosFace;
843 crossProductVector2 = edgeCoord3 - edgeCoord6;
844 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
845 interactionVolume.setFaceArea(faceArea, 8);
847 globalPosFace = interactionVolume.getFacePosition(9);
848 crossProductVector1 = centerPos - globalPosFace;
849 crossProductVector2 = edgeCoord3 - edgeCoord4;
850 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
851 interactionVolume.setFaceArea(faceArea, 9);
853 globalPosFace = interactionVolume.getFacePosition(10);
854 crossProductVector1 = centerPos - globalPosFace;
855 crossProductVector2 = edgeCoord4 - edgeCoord5;
856 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
857 interactionVolume.setFaceArea(faceArea, 10);
859 globalPosFace = interactionVolume.getFacePosition(11);
860 crossProductVector1 = centerPos - globalPosFace;
861 crossProductVector2 = edgeCoord5 - edgeCoord6;
862 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
863 interactionVolume.setFaceArea(faceArea, 11);
865 for (
const auto& intersection : intersections(problem_.gridView(), element))
867 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
869 if (intersection.neighbor())
871 auto outside = intersection.outside();
872 if (element.level() > outside.level())
874 interactionVolume.setSubVolumeElement(outside, 4);
875 interactionVolume.setSubVolumeElement(outside, 5);
876 interactionVolume.setSubVolumeElement(outside, 6);
877 interactionVolume.setSubVolumeElement(outside, 7);
886 else if (interactionVolume.getHangingNodeType() == InteractionVolume::fourSmallCellsEdge)
888 auto element1 = interactionVolume.getSubVolumeElement(0);
889 auto element2 = interactionVolume.getSubVolumeElement(1);
890 auto element3 = interactionVolume.getSubVolumeElement(2);
891 auto element4 = interactionVolume.getSubVolumeElement(3);
893 for (
const auto& intersection1 : intersections(problem_.gridView(), element1))
895 if (intersection1.neighbor())
897 if (intersection1.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
899 auto outside = intersection1.outside();
900 if (element1.level() > outside.level())
902 interactionVolume.setSubVolumeElement(outside, 4);
907 for (
const auto& intersection2 : intersections(problem_.gridView(), element2))
909 if (intersection2.neighbor())
911 if (intersection2.indexInInside() == interactionVolume.getIndexOnElement(1, 2))
913 auto outside = intersection2.outside();
914 if (element2.level() > outside.level())
916 interactionVolume.setSubVolumeElement(outside, 5);
923 for (
const auto& intersection3 : intersections(problem_.gridView(), element3))
925 if (intersection3.neighbor())
927 if (intersection3.indexInInside() == interactionVolume.getIndexOnElement(2, 2))
929 auto outside = intersection3.outside();
930 if (element3.level() > outside.level())
932 interactionVolume.setSubVolumeElement(outside, 6);
938 for (
const auto& intersection4 : intersections(problem_.gridView(), element4))
940 if (intersection4.neighbor())
942 if (intersection4.indexInInside() == interactionVolume.getIndexOnElement(3, 2))
944 auto outside = intersection4.outside();
945 if (element4.level() > outside.level())
947 interactionVolume.setSubVolumeElement(outside, 7);
955 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
956 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
957 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
958 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
959 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
961 DimVector crossProductVector1(0);
962 DimVector crossProductVector2(0);
964 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
965 crossProductVector1 = centerPos - globalPosFace;
966 crossProductVector2 = edgeCoord1 - edgeCoord3;
967 Scalar faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
968 interactionVolume.setFaceArea(faceArea, 0);
970 globalPosFace = interactionVolume.getFacePosition(1);
971 crossProductVector1 = centerPos - globalPosFace;
972 crossProductVector2 = edgeCoord1 - edgeCoord4;
973 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
974 interactionVolume.setFaceArea(faceArea, 1);
976 globalPosFace = interactionVolume.getFacePosition(2);
977 crossProductVector1 = centerPos - globalPosFace;
978 crossProductVector2 = edgeCoord1 - edgeCoord5;
979 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
980 interactionVolume.setFaceArea(faceArea, 2);
982 globalPosFace = interactionVolume.getFacePosition(3);
983 crossProductVector1 = centerPos - globalPosFace;
984 crossProductVector2 = edgeCoord1 - edgeCoord6;
985 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
986 interactionVolume.setFaceArea(faceArea, 3);
988 globalPosFace = interactionVolume.getFacePosition(8);
989 crossProductVector1 = centerPos - globalPosFace;
990 crossProductVector2 = edgeCoord3 - edgeCoord6;
991 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
992 interactionVolume.setFaceArea(faceArea, 8);
994 globalPosFace = interactionVolume.getFacePosition(9);
995 crossProductVector1 = centerPos - globalPosFace;
996 crossProductVector2 = edgeCoord3 - edgeCoord4;
997 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
998 interactionVolume.setFaceArea(faceArea, 9);
1000 globalPosFace = interactionVolume.getFacePosition(10);
1001 crossProductVector1 = centerPos - globalPosFace;
1002 crossProductVector2 = edgeCoord4 - edgeCoord5;
1003 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1004 interactionVolume.setFaceArea(faceArea, 10);
1006 globalPosFace = interactionVolume.getFacePosition(11);
1007 crossProductVector1 = centerPos - globalPosFace;
1008 crossProductVector2 = edgeCoord5 - edgeCoord6;
1009 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1010 interactionVolume.setFaceArea(faceArea, 11);
1012 auto element5 = interactionVolume.getSubVolumeElement(4);
1013 auto element6 = interactionVolume.getSubVolumeElement(5);
1014 auto element7 = interactionVolume.getSubVolumeElement(6);
1015 auto element8 = interactionVolume.getSubVolumeElement(7);
1017 if (element5 == element6)
1019 interactionVolume.setFacePosition(element5.geometry().center(), 4);
1020 interactionVolume.setFacePosition(element7.geometry().center(), 6);
1022 for (
const auto& intersection : intersections(problem_.gridView(), element5))
1024 if (intersection.neighbor())
1026 auto outside = intersection.outside();
1028 if (outside == element7 || outside == element8)
1030 int indexInInside = intersection.indexInInside();
1031 interactionVolume.setIndexOnElement(indexInInside, 4, 2);
1032 interactionVolume.setIndexOnElement(indexInInside, 5, 1);
1033 DimVector
normal = intersection.centerUnitOuterNormal();
1034 interactionVolume.setNormal(
normal, 4, 2);
1035 interactionVolume.setNormal(
normal, 5, 1);
1036 globalPosFace = intersection.geometry().center();
1037 interactionVolume.setFacePosition(globalPosFace, 5);
1038 interactionVolume.setFacePosition(globalPosFace, 7);
1039 int indexInOutside = intersection.indexInOutside();
1040 interactionVolume.setIndexOnElement(indexInOutside, 6, 1);
1041 interactionVolume.setIndexOnElement(indexInOutside, 7, 2);
1043 interactionVolume.setNormal(
normal, 6, 1);
1044 interactionVolume.setNormal(
normal, 7, 2);
1051 DimVector edgeCoord2(interactionVolume.getFacePosition(7));
1053 globalPosFace = interactionVolume.getFacePosition(5);
1054 crossProductVector1 = centerPos - globalPosFace;
1055 crossProductVector2 = edgeCoord2 - edgeCoord4;
1056 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1057 interactionVolume.setFaceArea(faceArea, 5);
1059 globalPosFace = interactionVolume.getFacePosition(7);
1060 crossProductVector1 = centerPos - globalPosFace;
1061 crossProductVector2 = edgeCoord2 - edgeCoord6;
1062 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1063 interactionVolume.setFaceArea(faceArea, 7);
1066 else if (element5 == element7)
1068 interactionVolume.setFacePosition(element6.geometry().center(), 5);
1069 interactionVolume.setFacePosition(element5.geometry().center(), 7);
1070 interactionVolume.setFacePosition(globalPosFace, 6);
1072 for (
const auto& intersection : intersections(problem_.gridView(), element5))
1074 if (intersection.neighbor())
1076 auto outside = intersection.outside();
1078 if (outside == element6 || outside == element8)
1080 int indexInInside = intersection.indexInInside();
1081 interactionVolume.setIndexOnElement(indexInInside, 4, 1);
1082 interactionVolume.setIndexOnElement(indexInInside, 6, 2);
1083 DimVector
normal = intersection.centerUnitOuterNormal();
1084 interactionVolume.setNormal(
normal, 4, 1);
1085 interactionVolume.setNormal(
normal, 6, 2);
1086 globalPosFace = intersection.geometry().center();
1087 interactionVolume.setFacePosition(globalPosFace, 4);
1088 interactionVolume.setFacePosition(globalPosFace, 6);
1089 int indexInOutside = intersection.indexInOutside();
1090 interactionVolume.setIndexOnElement(indexInOutside, 5, 2);
1091 interactionVolume.setIndexOnElement(indexInOutside, 7, 1);
1093 interactionVolume.setNormal(
normal, 5, 2);
1094 interactionVolume.setNormal(
normal, 7, 1);
1100 DimVector edgeCoord2(interactionVolume.getFacePosition(4));
1102 globalPosFace = interactionVolume.getFacePosition(4);
1103 crossProductVector1 = centerPos - globalPosFace;
1104 crossProductVector2 = edgeCoord2 - edgeCoord3;
1105 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1106 interactionVolume.setFaceArea(faceArea, 4);
1108 globalPosFace = interactionVolume.getFacePosition(6);
1109 crossProductVector1 = centerPos - globalPosFace;
1110 crossProductVector2 = edgeCoord2 - edgeCoord5;
1111 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1112 interactionVolume.setFaceArea(faceArea, 6);
1117 const ElementGeometry& geometry = element5.geometry();
1119 const auto refElement = referenceElement(geometry);
1121 int oldSubVolumElemIdx = IndexTranslator::getOldElemIdxFromNewFaceIdxto0(zeroFaceIdx, 4);
1122 int oldEdgeIdx = IndexTranslator::getOldEdgeIdxFromNewFaceIdxto0(zeroFaceIdx, 1);
1124 DimVector edgeCoord(0.0);
1125 switch (oldSubVolumElemIdx)
1132 edgeCoord = geometry.global(refElement.position(9, dim - 1));
1135 edgeCoord = geometry.global(refElement.position(3, dim - 1));
1138 edgeCoord = geometry.global(refElement.position(11, dim - 1));
1149 edgeCoord = geometry.global(refElement.position(2, dim - 1));
1152 edgeCoord = geometry.global(refElement.position(8, dim - 1));
1155 edgeCoord = geometry.global(refElement.position(11, dim - 1));
1166 edgeCoord = geometry.global(refElement.position(1, dim - 1));
1169 edgeCoord = geometry.global(refElement.position(9, dim - 1));
1172 edgeCoord = geometry.global(refElement.position(10, dim - 1));
1183 edgeCoord = geometry.global(refElement.position(0, dim - 1));
1186 edgeCoord = geometry.global(refElement.position(8, dim - 1));
1189 edgeCoord = geometry.global(refElement.position(10, dim - 1));
1200 edgeCoord = geometry.global(refElement.position(3, dim - 1));
1203 edgeCoord = geometry.global(refElement.position(5, dim - 1));
1206 edgeCoord = geometry.global(refElement.position(7, dim - 1));
1217 edgeCoord = geometry.global(refElement.position(2, dim - 1));
1220 edgeCoord = geometry.global(refElement.position(4, dim - 1));
1223 edgeCoord = geometry.global(refElement.position(7, dim - 1));
1234 edgeCoord = geometry.global(refElement.position(1, dim - 1));
1237 edgeCoord = geometry.global(refElement.position(5, dim - 1));
1240 edgeCoord = geometry.global(refElement.position(6, dim - 1));
1251 edgeCoord = geometry.global(refElement.position(0, dim - 1));
1254 edgeCoord = geometry.global(refElement.position(4, dim - 1));
1257 edgeCoord = geometry.global(refElement.position(6, dim - 1));
1265 interactionVolume.setEdgePosition(edgeCoord, 1);
1267 else if (interactionVolume.getHangingNodeType() == InteractionVolume::fourSmallCellsDiag)
1269 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
1270 DimVector edgeCoord2(interactionVolume.getEdgePosition(1));
1271 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
1272 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
1273 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
1274 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
1276 DimVector crossProductVector1(0);
1277 DimVector crossProductVector2(0);
1279 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
1280 crossProductVector1 = centerPos - globalPosFace;
1281 crossProductVector2 = edgeCoord1 - edgeCoord3;
1282 Scalar faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1283 interactionVolume.setFaceArea(faceArea, 0);
1285 globalPosFace = interactionVolume.getFacePosition(1);
1286 crossProductVector1 = centerPos - globalPosFace;
1287 crossProductVector2 = edgeCoord1 - edgeCoord4;
1288 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1289 interactionVolume.setFaceArea(faceArea, 1);
1291 globalPosFace = interactionVolume.getFacePosition(3);
1292 crossProductVector1 = centerPos - globalPosFace;
1293 crossProductVector2 = edgeCoord1 - edgeCoord6;
1294 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1295 interactionVolume.setFaceArea(faceArea, 3);
1297 globalPosFace = interactionVolume.getFacePosition(5);
1298 crossProductVector1 = centerPos - globalPosFace;
1299 crossProductVector2 = edgeCoord2 - edgeCoord4;
1300 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1301 interactionVolume.setFaceArea(faceArea, 5);
1303 globalPosFace = interactionVolume.getFacePosition(6);
1304 crossProductVector1 = centerPos - globalPosFace;
1305 crossProductVector2 = edgeCoord2 - edgeCoord5;
1306 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1307 interactionVolume.setFaceArea(faceArea, 6);
1309 globalPosFace = interactionVolume.getFacePosition(7);
1310 crossProductVector1 = centerPos - globalPosFace;
1311 crossProductVector2 = edgeCoord2 - edgeCoord6;
1312 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1313 interactionVolume.setFaceArea(faceArea, 7);
1315 globalPosFace = interactionVolume.getFacePosition(8);
1316 crossProductVector1 = centerPos - globalPosFace;
1317 crossProductVector2 = edgeCoord3 - edgeCoord6;
1318 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1319 interactionVolume.setFaceArea(faceArea, 8);
1321 globalPosFace = interactionVolume.getFacePosition(9);
1322 crossProductVector1 = centerPos - globalPosFace;
1323 crossProductVector2 = edgeCoord3 - edgeCoord4;
1324 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1325 interactionVolume.setFaceArea(faceArea, 9);
1327 globalPosFace = interactionVolume.getFacePosition(10);
1328 crossProductVector1 = centerPos - globalPosFace;
1329 crossProductVector2 = edgeCoord4 - edgeCoord5;
1330 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1331 interactionVolume.setFaceArea(faceArea, 10);
1333 globalPosFace = interactionVolume.getFacePosition(11);
1334 crossProductVector1 = centerPos - globalPosFace;
1335 crossProductVector2 = edgeCoord5 - edgeCoord6;
1336 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1337 interactionVolume.setFaceArea(faceArea, 11);
1339 auto element1 = interactionVolume.getSubVolumeElement(0);
1341 bool hasFaceOne =
false;
1342 bool hasFaceTwo =
false;
1343 for (
const auto& intersection : intersections(problem_.gridView(), element1))
1345 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(0, 1))
1347 if (intersection.neighbor())
1349 auto outside = intersection.outside();
1350 if (element1.level() > outside.level())
1352 interactionVolume.setSubVolumeElement(outside, 2);
1353 interactionVolume.setSubVolumeElement(outside, 3);
1361 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
1363 if (intersection.neighbor())
1365 auto outside = intersection.outside();
1366 if (element1.level() > outside.level())
1368 interactionVolume.setSubVolumeElement(outside, 4);
1369 interactionVolume.setSubVolumeElement(outside, 5);
1385 InteractionVolume hangingNodeVolume(problem_.gridView().grid());
1387 std::vector<int> elemIdxOld;
1388 for (
int i = 0; i < InteractionVolume::subVolumeTotalNum; i++)
1390 if (interactionVolume.hasSubVolumeElement(i))
1391 elemIdxOld.push_back(i);
1394 std::set<int> zeroFaceIdxVec;
1395 for (
int i = 0; i < 6; i++)
1397 for (
int j = 0; j < 6; j++)
1399 int fIdx = IndexTranslator::getFaceIndexFromElements(elemIdxOld[i], elemIdxOld[j]);
1401 zeroFaceIdxVec.insert(fIdx);
1405 std::vector<int> elemIdxNew(6);
1406 int zeroFaceIdx = 0;
1409 std::set<int>::iterator it = zeroFaceIdxVec.begin();
1410 for (; it != zeroFaceIdxVec.end(); ++it)
1415 for (
int i = 0; i < 6; i++)
1417 elemIdxNew[i] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elemIdxOld[i]);
1418 if (elemIdxNew[i] == 6 || elemIdxNew[i] == 7)
1426 for (
int i = 0; i < 6; i++)
1428 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[i]),
1436 for (
int elem = 0; elem < InteractionVolume::subVolumeTotalNum; elem++)
1438 for (
int i = 0; i < 3; i++)
1440 int faceIdxOld = IndexTranslator::getFaceIndexFromSubVolume(elem, i);
1441 int faceIdxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, faceIdxOld);
1443 for (
int j = 0; j < 3; j++)
1445 int elemNew = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elem);
1446 int faceIdxNewTest = IndexTranslator::getFaceIndexFromSubVolume(elemNew, j);
1448 if (faceIdxNew == faceIdxNewTest)
1450 hangingNodeVolume.setNormal(interactionVolume.getNormal(elem, i),
1452 hangingNodeVolume.setIndexOnElement(
1453 interactionVolume.getIndexOnElement(elem, i), elemNew, j);
1459 for (
int i = 0; i < InteractionVolume::fluxFacesTotalNum; i++)
1461 int idxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, i);
1462 hangingNodeVolume.setFaceArea(interactionVolume.getFaceArea(i), idxNew);
1463 hangingNodeVolume.setFacePosition(interactionVolume.getFacePosition(i), idxNew);
1466 for (
int i = 0; i < InteractionVolume::fluxEdgesTotalNum; i++)
1468 int idxNew = IndexTranslator::getNewEdgeIdxFromOldFaceIdxto0(zeroFaceIdx, i);
1469 hangingNodeVolume.setEdgePosition(interactionVolume.getEdgePosition(i), idxNew);
1472 interactionVolume = hangingNodeVolume;
1474 interactionVolume.setCenterPosition(centerPos);
1476 interactionVolume.setHangingNodeType(InteractionVolume::sixSmallCells);
1478 auto element3 = interactionVolume.getSubVolumeElement(2);
1480 for (
const auto& intersection : intersections(problem_.gridView(), element3))
1482 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(2, 2))
1484 if (intersection.neighbor())
1486 auto outside = intersection.outside();
1487 if (element3.level() > outside.level())
1489 interactionVolume.setSubVolumeElement(outside, 6);
1490 interactionVolume.setSubVolumeElement(outside, 7);
1498 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
1499 DimVector edgeCoord2(interactionVolume.getEdgePosition(1));
1500 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
1501 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
1502 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
1503 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
1505 DimVector crossProductVector1(0);
1506 DimVector crossProductVector2(0);
1508 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
1509 crossProductVector1 = centerPos - globalPosFace;
1510 crossProductVector2 = edgeCoord1 - edgeCoord3;
1511 Scalar faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1512 interactionVolume.setFaceArea(faceArea, 0);
1514 globalPosFace = interactionVolume.getFacePosition(1);
1515 crossProductVector1 = centerPos - globalPosFace;
1516 crossProductVector2 = edgeCoord1 - edgeCoord4;
1517 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1518 interactionVolume.setFaceArea(faceArea, 1);
1520 globalPosFace = interactionVolume.getFacePosition(2);
1521 crossProductVector1 = centerPos - globalPosFace;
1522 crossProductVector2 = edgeCoord1 - edgeCoord5;
1523 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1524 interactionVolume.setFaceArea(faceArea, 2);
1526 globalPosFace = interactionVolume.getFacePosition(3);
1527 crossProductVector1 = centerPos - globalPosFace;
1528 crossProductVector2 = edgeCoord1 - edgeCoord6;
1529 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1530 interactionVolume.setFaceArea(faceArea, 3);
1532 globalPosFace = interactionVolume.getFacePosition(4);
1533 crossProductVector1 = centerPos - globalPosFace;
1534 crossProductVector2 = edgeCoord2 - edgeCoord3;
1535 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1536 interactionVolume.setFaceArea(faceArea, 4);
1538 globalPosFace = interactionVolume.getFacePosition(5);
1539 crossProductVector1 = centerPos - globalPosFace;
1540 crossProductVector2 = edgeCoord2 - edgeCoord4;
1541 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1542 interactionVolume.setFaceArea(faceArea, 5);
1544 globalPosFace = interactionVolume.getFacePosition(7);
1545 crossProductVector1 = centerPos - globalPosFace;
1546 crossProductVector2 = edgeCoord2 - edgeCoord6;
1547 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1548 interactionVolume.setFaceArea(faceArea, 7);
1550 globalPosFace = interactionVolume.getFacePosition(8);
1551 crossProductVector1 = centerPos - globalPosFace;
1552 crossProductVector2 = edgeCoord3 - edgeCoord6;
1553 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1554 interactionVolume.setFaceArea(faceArea, 8);
1556 globalPosFace = interactionVolume.getFacePosition(9);
1557 crossProductVector1 = centerPos - globalPosFace;
1558 crossProductVector2 = edgeCoord3 - edgeCoord4;
1559 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1560 interactionVolume.setFaceArea(faceArea, 9);
1562 globalPosFace = interactionVolume.getFacePosition(10);
1563 crossProductVector1 = centerPos - globalPosFace;
1564 crossProductVector2 = edgeCoord4 - edgeCoord5;
1565 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1566 interactionVolume.setFaceArea(faceArea, 10);
1568 globalPosFace = interactionVolume.getFacePosition(11);
1569 crossProductVector1 = centerPos - globalPosFace;
1570 crossProductVector2 = edgeCoord5 - edgeCoord6;
1571 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1572 interactionVolume.setFaceArea(faceArea, 11);
1577 interactionVolume.printInteractionVolumeInfo();
1578 DUNE_THROW(Dune::NotImplemented,
"Hanging node shape not implemented");
1585template<
class TypeTag>
1588 std::vector < std::vector<int> > elemVertMap(problem_.gridView().size(dim), std::vector<int>(8, -1));
1591 for (
const auto& element : elements(problem_.gridView()))
1592 asImp_().storeSubVolumeElements(element, elemVertMap);
1594 for (
unsigned int i = 0; i < asImp_().interactionVolumes_.size(); i++)
1595 if (asImp_().interactionVolumes_[i].getElementNumber() == 0)
1596 asImp_().interactionVolumes_[i].printInteractionVolumeInfo();
1599 for (
const auto& element : elements(problem_.gridView()))
1600 asImp_().storeIntersectionInfo(element, elemVertMap);
1602 faceVertices_.clear();
1603 faceVertices_.resize(problem_.gridView().size(0));
1607 for (
const auto& vertex : vertices(problem_.gridView()))
1609 int vIdxGlobal = problem_.variables().index(vertex);
1611 InteractionVolume& interactionVolume = asImp_().interactionVolumes_[vIdxGlobal];
1613 if (interactionVolume.getElementNumber() == 8)
1615 asImp_().storeInnerInteractionVolume(interactionVolume, vertex);
1617 else if (interactionVolume.isBoundaryInteractionVolume())
1619 asImp_().storeBoundaryInteractionVolume(interactionVolume, vertex);
1624 storeHangingNodeInteractionVolume(interactionVolume, vertex);
1627 if (!interactionVolume.isBoundaryInteractionVolume())
1629 auto element1 = interactionVolume.getSubVolumeElement(0);
1630 auto element2 = interactionVolume.getSubVolumeElement(1);
1631 auto element3 = interactionVolume.getSubVolumeElement(2);
1632 auto element4 = interactionVolume.getSubVolumeElement(3);
1633 auto element5 = interactionVolume.getSubVolumeElement(4);
1634 auto element6 = interactionVolume.getSubVolumeElement(5);
1635 auto element7 = interactionVolume.getSubVolumeElement(6);
1636 auto element8 = interactionVolume.getSubVolumeElement(7);
1638 int globalIdx1 = problem_.variables().index(element1);
1639 int globalIdx2 = problem_.variables().index(element2);
1640 int globalIdx3 = problem_.variables().index(element3);
1641 int globalIdx4 = problem_.variables().index(element4);
1642 int globalIdx5 = problem_.variables().index(element5);
1643 int globalIdx6 = problem_.variables().index(element6);
1644 int globalIdx7 = problem_.variables().index(element7);
1645 int globalIdx8 = problem_.variables().index(element8);
1647 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(0, 0), globalIdx1, interactionVolume.getIndexOnElement(0, 0));
1648 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(0, 1), globalIdx1, interactionVolume.getIndexOnElement(0, 1));
1649 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(0, 2), globalIdx1, interactionVolume.getIndexOnElement(0, 2));
1650 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(1, 0), globalIdx2, interactionVolume.getIndexOnElement(1, 0));
1651 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(1, 1), globalIdx2, interactionVolume.getIndexOnElement(1, 1));
1652 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(1, 2), globalIdx2, interactionVolume.getIndexOnElement(1, 2));
1653 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(2, 0), globalIdx3, interactionVolume.getIndexOnElement(2, 0));
1654 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(3, 1), globalIdx4, interactionVolume.getIndexOnElement(3, 1));
1655 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(4, 0), globalIdx5, interactionVolume.getIndexOnElement(4, 0));
1656 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(5, 0), globalIdx6, interactionVolume.getIndexOnElement(5, 0));
1658 faceVertices_[globalIdx1][interactionVolume.getIndexOnElement(0, 0)].insert(vIdxGlobal);
1659 faceVertices_[globalIdx1][interactionVolume.getIndexOnElement(0, 1)].insert(vIdxGlobal);
1660 faceVertices_[globalIdx1][interactionVolume.getIndexOnElement(0, 2)].insert(vIdxGlobal);
1661 faceVertices_[globalIdx2][interactionVolume.getIndexOnElement(1, 0)].insert(vIdxGlobal);
1662 faceVertices_[globalIdx2][interactionVolume.getIndexOnElement(1, 1)].insert(vIdxGlobal);
1663 faceVertices_[globalIdx2][interactionVolume.getIndexOnElement(1, 2)].insert(vIdxGlobal);
1664 faceVertices_[globalIdx3][interactionVolume.getIndexOnElement(2, 0)].insert(vIdxGlobal);
1665 faceVertices_[globalIdx4][interactionVolume.getIndexOnElement(3, 1)].insert(vIdxGlobal);
1666 faceVertices_[globalIdx5][interactionVolume.getIndexOnElement(4, 0)].insert(vIdxGlobal);
1667 faceVertices_[globalIdx6][interactionVolume.getIndexOnElement(5, 0)].insert(vIdxGlobal);
1669 if (interactionVolume.getHangingNodeType() != InteractionVolume::twoSmallCells)
1671 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(2, 1), globalIdx3, interactionVolume.getIndexOnElement(2, 1));
1672 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(2, 2), globalIdx3, interactionVolume.getIndexOnElement(2, 2));
1673 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(3, 0), globalIdx4, interactionVolume.getIndexOnElement(3, 0));
1674 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(3, 2), globalIdx4, interactionVolume.getIndexOnElement(3, 2));
1675 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(4, 1), globalIdx5, interactionVolume.getIndexOnElement(4, 1));
1676 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(4, 2), globalIdx5, interactionVolume.getIndexOnElement(4, 2));
1677 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(5, 1), globalIdx6, interactionVolume.getIndexOnElement(5, 1));
1678 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(5, 2), globalIdx6, interactionVolume.getIndexOnElement(5, 2));
1679 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(6, 0), globalIdx7, interactionVolume.getIndexOnElement(6, 0));
1680 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(6, 1), globalIdx7, interactionVolume.getIndexOnElement(6, 1));
1681 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(6, 2), globalIdx7, interactionVolume.getIndexOnElement(6, 2));
1682 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(7, 0), globalIdx8, interactionVolume.getIndexOnElement(7, 0));
1683 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(7, 1), globalIdx8, interactionVolume.getIndexOnElement(7, 1));
1684 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(7, 2), globalIdx8, interactionVolume.getIndexOnElement(7, 2));
1686 faceVertices_[globalIdx3][interactionVolume.getIndexOnElement(2, 1)].insert(vIdxGlobal);
1687 faceVertices_[globalIdx3][interactionVolume.getIndexOnElement(2, 2)].insert(vIdxGlobal);
1688 faceVertices_[globalIdx4][interactionVolume.getIndexOnElement(3, 0)].insert(vIdxGlobal);
1689 faceVertices_[globalIdx4][interactionVolume.getIndexOnElement(3, 2)].insert(vIdxGlobal);
1690 faceVertices_[globalIdx5][interactionVolume.getIndexOnElement(4, 1)].insert(vIdxGlobal);
1691 faceVertices_[globalIdx5][interactionVolume.getIndexOnElement(4, 2)].insert(vIdxGlobal);
1692 faceVertices_[globalIdx6][interactionVolume.getIndexOnElement(5, 1)].insert(vIdxGlobal);
1693 faceVertices_[globalIdx6][interactionVolume.getIndexOnElement(5, 2)].insert(vIdxGlobal);
1694 faceVertices_[globalIdx7][interactionVolume.getIndexOnElement(6, 0)].insert(vIdxGlobal);
1695 faceVertices_[globalIdx7][interactionVolume.getIndexOnElement(6, 1)].insert(vIdxGlobal);
1696 faceVertices_[globalIdx7][interactionVolume.getIndexOnElement(6, 2)].insert(vIdxGlobal);
1697 faceVertices_[globalIdx8][interactionVolume.getIndexOnElement(7, 0)].insert(vIdxGlobal);
1698 faceVertices_[globalIdx8][interactionVolume.getIndexOnElement(7, 1)].insert(vIdxGlobal);
1699 faceVertices_[globalIdx8][interactionVolume.getIndexOnElement(7, 2)].insert(vIdxGlobal);
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...
Vector normal(const Vector &v)
Create a vector normal to the given one (v is expected to be non-zero)
Definition: geometry/normal.hh:36
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:640
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
Interactionvolume container for 3-d MPFA L-method.
Definition: 3dinteractionvolumecontainer.hh:46
GetPropType< TypeTag, Properties::MPFAInteractionVolume > InteractionVolume
Type for storing an MPFA-interaction-volume. (Usually of type FvMpfaL3dInteractionVolume or FvMpfaL3d...
Definition: 3dinteractionvolumecontainer.hh:92
InteractionVolume & interactionVolume(int vertexIdx)
Returns an interaction volume.
Definition: 3dinteractionvolumecontainer.hh:141
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:1586
std::set< int > & faceVerticeIndices(int eIdxGlobal, int fIdx)
Returns the set of vertices on an element face.
Definition: 3dinteractionvolumecontaineradaptive.hh:102
FvMpfaL3dInteractionVolumeContainerAdaptive(Problem &problem)
Constructs a FvMpfaL3dInteractionVolumeContainerAdaptive object.
Definition: 3dinteractionvolumecontaineradaptive.hh:111
Properties for a MPFA method.
Base file for properties related to sequential IMPET algorithms.