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 edgeCoord3(interactionVolume.getEdgePosition(2));
610 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
611 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
613 DimVector crossProductVector1(0);
614 DimVector crossProductVector2(0);
616 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
617 crossProductVector1 = centerPos - globalPosFace;
618 crossProductVector2 = edgeCoord1 - edgeCoord3;
619 Scalar faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
620 interactionVolume.setFaceArea(faceArea, 0);
622 globalPosFace = interactionVolume.getFacePosition(1);
623 crossProductVector1 = centerPos - globalPosFace;
624 crossProductVector2 = edgeCoord1 - edgeCoord4;
625 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
626 interactionVolume.setFaceArea(faceArea, 1);
628 globalPosFace = interactionVolume.getFacePosition(3);
629 crossProductVector1 = centerPos - globalPosFace;
630 crossProductVector2 = edgeCoord1 - edgeCoord6;
631 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
632 interactionVolume.setFaceArea(faceArea, 3);
634 globalPosFace = interactionVolume.getFacePosition(8);
635 crossProductVector1 = centerPos - globalPosFace;
636 crossProductVector2 = edgeCoord3 - edgeCoord6;
637 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
638 interactionVolume.setFaceArea(faceArea, 8);
640 globalPosFace = interactionVolume.getFacePosition(9);
641 crossProductVector1 = centerPos - globalPosFace;
642 crossProductVector2 = edgeCoord3 - edgeCoord4;
643 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
644 interactionVolume.setFaceArea(faceArea, 9);
651 InteractionVolume hangingNodeVolume(problem_.gridView().grid());
653 std::vector<int> elemIdxOld;
654 for (
int i = 0; i < InteractionVolume::subVolumeTotalNum; i++)
656 if (interactionVolume.hasSubVolumeElement(i))
657 elemIdxOld.push_back(i);
660 std::set<int> zeroFaceIdxVec;
661 for (
int i = 0; i < 4; i++)
663 for (
int j = 0; j < 4; j++)
665 int fIdx = IndexTranslator::getFaceIndexFromElements(elemIdxOld[i], elemIdxOld[j]);
667 zeroFaceIdxVec.insert(fIdx);
671 std::vector<int> elemIdxNew(4);
674 if (zeroFaceIdxVec.size() == 2)
676 hangingNodeVolume.setHangingNodeType(InteractionVolume::fourSmallCellsDiag);
678 if (zeroFaceIdxVec.find(0) == zeroFaceIdxVec.end())
679 zeroFaceIdx = *zeroFaceIdxVec.begin();
681 for (
int i = 0; i < 4; i++)
683 elemIdxNew[i] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elemIdxOld[i]);
684 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[i]),
688 else if (zeroFaceIdxVec.size() == 4)
690 std::set<int>::iterator it = zeroFaceIdxVec.begin();
691 for (; it != zeroFaceIdxVec.end(); ++it)
696 for (
int i = 0; i < 4; i++)
698 elemIdxNew[i] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx,
700 if (elemIdxNew[i] == 4 || elemIdxNew[i] == 5 || elemIdxNew[i] == 6 || elemIdxNew[i] == 7)
708 for (
int i = 0; i < 4; i++)
710 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[i]), elemIdxNew[i]);
719 std::set<int>::iterator it = zeroFaceIdxVec.begin();
720 std::cout<<
"zeroFaceIdxVec = ";
721 for(;it != zeroFaceIdxVec.end(); ++it)
728 for (
int elem = 0; elem < InteractionVolume::subVolumeTotalNum; elem++)
730 for (
int i = 0; i < 3; i++)
732 int faceIdxOld = IndexTranslator::getFaceIndexFromSubVolume(elem, i);
733 int faceIdxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, faceIdxOld);
735 for (
int j = 0; j < 3; j++)
737 int elemNew = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elem);
738 int faceIdxNewTest = IndexTranslator::getFaceIndexFromSubVolume(elemNew, j);
740 if (faceIdxNew == faceIdxNewTest)
742 hangingNodeVolume.setNormal(interactionVolume.getNormal(elem, i),
744 hangingNodeVolume.setIndexOnElement(interactionVolume.getIndexOnElement(elem, i), elemNew, j);
750 for (
int i = 0; i < InteractionVolume::fluxFacesTotalNum; i++)
752 int idxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, i);
753 hangingNodeVolume.setFacePosition(interactionVolume.getFacePosition(i), idxNew);
756 for (
int i = 0; i < InteractionVolume::fluxEdgesTotalNum; i++)
758 int idxNew = IndexTranslator::getNewEdgeIdxFromOldFaceIdxto0(zeroFaceIdx, i);
759 hangingNodeVolume.setEdgePosition(interactionVolume.getEdgePosition(i), idxNew);
762 interactionVolume = hangingNodeVolume;
764 interactionVolume.setCenterPosition(centerPos);
766 if (zeroFaceIdxVec.size() == 4)
768 auto element1 = interactionVolume.getSubVolumeElement(0);
769 auto element4 = interactionVolume.getSubVolumeElement(3);
771 auto outside1 = element1;
772 auto outside4 = element4;
774 for (
const auto& intersection1 : intersections(problem_.gridView(), element1))
776 if (intersection1.neighbor())
778 if (intersection1.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
780 outside1 = intersection1.outside();
785 for (
const auto& intersection4 : intersections(problem_.gridView(), element4))
787 if (intersection4.neighbor())
789 if (intersection4.indexInInside() == interactionVolume.getIndexOnElement(3, 2))
791 outside4 = intersection4.outside();
796 if (outside1 != outside4)
798 interactionVolume.setHangingNodeType(InteractionVolume::fourSmallCellsEdge);
800 else if (outside1 == outside4)
802 interactionVolume.setHangingNodeType(InteractionVolume::fourSmallCellsFace);
806 if (interactionVolume.getHangingNodeType() == InteractionVolume::fourSmallCellsFace)
808 auto element = interactionVolume.getSubVolumeElement(0);
810 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
811 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
812 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
813 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
814 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
816 DimVector crossProductVector1(0);
817 DimVector crossProductVector2(0);
819 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
820 crossProductVector1 = centerPos - globalPosFace;
821 crossProductVector2 = edgeCoord1 - edgeCoord3;
822 Scalar faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
823 interactionVolume.setFaceArea(faceArea, 0);
825 globalPosFace = interactionVolume.getFacePosition(1);
826 crossProductVector1 = centerPos - globalPosFace;
827 crossProductVector2 = edgeCoord1 - edgeCoord4;
828 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
829 interactionVolume.setFaceArea(faceArea, 1);
831 globalPosFace = interactionVolume.getFacePosition(2);
832 crossProductVector1 = centerPos - globalPosFace;
833 crossProductVector2 = edgeCoord1 - edgeCoord5;
834 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
835 interactionVolume.setFaceArea(faceArea, 2);
837 globalPosFace = interactionVolume.getFacePosition(3);
838 crossProductVector1 = centerPos - globalPosFace;
839 crossProductVector2 = edgeCoord1 - edgeCoord6;
840 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
841 interactionVolume.setFaceArea(faceArea, 3);
843 globalPosFace = interactionVolume.getFacePosition(8);
844 crossProductVector1 = centerPos - globalPosFace;
845 crossProductVector2 = edgeCoord3 - edgeCoord6;
846 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
847 interactionVolume.setFaceArea(faceArea, 8);
849 globalPosFace = interactionVolume.getFacePosition(9);
850 crossProductVector1 = centerPos - globalPosFace;
851 crossProductVector2 = edgeCoord3 - edgeCoord4;
852 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
853 interactionVolume.setFaceArea(faceArea, 9);
855 globalPosFace = interactionVolume.getFacePosition(10);
856 crossProductVector1 = centerPos - globalPosFace;
857 crossProductVector2 = edgeCoord4 - edgeCoord5;
858 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
859 interactionVolume.setFaceArea(faceArea, 10);
861 globalPosFace = interactionVolume.getFacePosition(11);
862 crossProductVector1 = centerPos - globalPosFace;
863 crossProductVector2 = edgeCoord5 - edgeCoord6;
864 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
865 interactionVolume.setFaceArea(faceArea, 11);
867 for (
const auto& intersection : intersections(problem_.gridView(), element))
869 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
871 if (intersection.neighbor())
873 auto outside = intersection.outside();
874 if (element.level() > outside.level())
876 interactionVolume.setSubVolumeElement(outside, 4);
877 interactionVolume.setSubVolumeElement(outside, 5);
878 interactionVolume.setSubVolumeElement(outside, 6);
879 interactionVolume.setSubVolumeElement(outside, 7);
888 else if (interactionVolume.getHangingNodeType() == InteractionVolume::fourSmallCellsEdge)
890 auto element1 = interactionVolume.getSubVolumeElement(0);
891 auto element2 = interactionVolume.getSubVolumeElement(1);
892 auto element3 = interactionVolume.getSubVolumeElement(2);
893 auto element4 = interactionVolume.getSubVolumeElement(3);
895 for (
const auto& intersection1 : intersections(problem_.gridView(), element1))
897 if (intersection1.neighbor())
899 if (intersection1.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
901 auto outside = intersection1.outside();
902 if (element1.level() > outside.level())
904 interactionVolume.setSubVolumeElement(outside, 4);
909 for (
const auto& intersection2 : intersections(problem_.gridView(), element2))
911 if (intersection2.neighbor())
913 if (intersection2.indexInInside() == interactionVolume.getIndexOnElement(1, 2))
915 auto outside = intersection2.outside();
916 if (element2.level() > outside.level())
918 interactionVolume.setSubVolumeElement(outside, 5);
925 for (
const auto& intersection3 : intersections(problem_.gridView(), element3))
927 if (intersection3.neighbor())
929 if (intersection3.indexInInside() == interactionVolume.getIndexOnElement(2, 2))
931 auto outside = intersection3.outside();
932 if (element3.level() > outside.level())
934 interactionVolume.setSubVolumeElement(outside, 6);
940 for (
const auto& intersection4 : intersections(problem_.gridView(), element4))
942 if (intersection4.neighbor())
944 if (intersection4.indexInInside() == interactionVolume.getIndexOnElement(3, 2))
946 auto outside = intersection4.outside();
947 if (element4.level() > outside.level())
949 interactionVolume.setSubVolumeElement(outside, 7);
957 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
958 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
959 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
960 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
961 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
963 DimVector crossProductVector1(0);
964 DimVector crossProductVector2(0);
966 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
967 crossProductVector1 = centerPos - globalPosFace;
968 crossProductVector2 = edgeCoord1 - edgeCoord3;
969 Scalar faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
970 interactionVolume.setFaceArea(faceArea, 0);
972 globalPosFace = interactionVolume.getFacePosition(1);
973 crossProductVector1 = centerPos - globalPosFace;
974 crossProductVector2 = edgeCoord1 - edgeCoord4;
975 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
976 interactionVolume.setFaceArea(faceArea, 1);
978 globalPosFace = interactionVolume.getFacePosition(2);
979 crossProductVector1 = centerPos - globalPosFace;
980 crossProductVector2 = edgeCoord1 - edgeCoord5;
981 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
982 interactionVolume.setFaceArea(faceArea, 2);
984 globalPosFace = interactionVolume.getFacePosition(3);
985 crossProductVector1 = centerPos - globalPosFace;
986 crossProductVector2 = edgeCoord1 - edgeCoord6;
987 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
988 interactionVolume.setFaceArea(faceArea, 3);
990 globalPosFace = interactionVolume.getFacePosition(8);
991 crossProductVector1 = centerPos - globalPosFace;
992 crossProductVector2 = edgeCoord3 - edgeCoord6;
993 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
994 interactionVolume.setFaceArea(faceArea, 8);
996 globalPosFace = interactionVolume.getFacePosition(9);
997 crossProductVector1 = centerPos - globalPosFace;
998 crossProductVector2 = edgeCoord3 - edgeCoord4;
999 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1000 interactionVolume.setFaceArea(faceArea, 9);
1002 globalPosFace = interactionVolume.getFacePosition(10);
1003 crossProductVector1 = centerPos - globalPosFace;
1004 crossProductVector2 = edgeCoord4 - edgeCoord5;
1005 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1006 interactionVolume.setFaceArea(faceArea, 10);
1008 globalPosFace = interactionVolume.getFacePosition(11);
1009 crossProductVector1 = centerPos - globalPosFace;
1010 crossProductVector2 = edgeCoord5 - edgeCoord6;
1011 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1012 interactionVolume.setFaceArea(faceArea, 11);
1014 auto element5 = interactionVolume.getSubVolumeElement(4);
1015 auto element6 = interactionVolume.getSubVolumeElement(5);
1016 auto element7 = interactionVolume.getSubVolumeElement(6);
1017 auto element8 = interactionVolume.getSubVolumeElement(7);
1019 if (element5 == element6)
1021 interactionVolume.setFacePosition(element5.geometry().center(), 4);
1022 interactionVolume.setFacePosition(element7.geometry().center(), 6);
1024 for (
const auto& intersection : intersections(problem_.gridView(), element5))
1026 if (intersection.neighbor())
1028 auto outside = intersection.outside();
1030 if (outside == element7 || outside == element8)
1032 int indexInInside = intersection.indexInInside();
1033 interactionVolume.setIndexOnElement(indexInInside, 4, 2);
1034 interactionVolume.setIndexOnElement(indexInInside, 5, 1);
1035 DimVector normal = intersection.centerUnitOuterNormal();
1036 interactionVolume.setNormal(normal, 4, 2);
1037 interactionVolume.setNormal(normal, 5, 1);
1038 globalPosFace = intersection.geometry().center();
1039 interactionVolume.setFacePosition(globalPosFace, 5);
1040 interactionVolume.setFacePosition(globalPosFace, 7);
1041 int indexInOutside = intersection.indexInOutside();
1042 interactionVolume.setIndexOnElement(indexInOutside, 6, 1);
1043 interactionVolume.setIndexOnElement(indexInOutside, 7, 2);
1045 interactionVolume.setNormal(normal, 6, 1);
1046 interactionVolume.setNormal(normal, 7, 2);
1053 DimVector edgeCoord2(interactionVolume.getFacePosition(7));
1055 globalPosFace = interactionVolume.getFacePosition(5);
1056 crossProductVector1 = centerPos - globalPosFace;
1057 crossProductVector2 = edgeCoord2 - edgeCoord4;
1058 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1059 interactionVolume.setFaceArea(faceArea, 5);
1061 globalPosFace = interactionVolume.getFacePosition(7);
1062 crossProductVector1 = centerPos - globalPosFace;
1063 crossProductVector2 = edgeCoord2 - edgeCoord6;
1064 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1065 interactionVolume.setFaceArea(faceArea, 7);
1068 else if (element5 == element7)
1070 interactionVolume.setFacePosition(element6.geometry().center(), 5);
1071 interactionVolume.setFacePosition(element5.geometry().center(), 7);
1072 interactionVolume.setFacePosition(globalPosFace, 6);
1074 for (
const auto& intersection : intersections(problem_.gridView(), element5))
1076 if (intersection.neighbor())
1078 auto outside = intersection.outside();
1080 if (outside == element6 || outside == element8)
1082 int indexInInside = intersection.indexInInside();
1083 interactionVolume.setIndexOnElement(indexInInside, 4, 1);
1084 interactionVolume.setIndexOnElement(indexInInside, 6, 2);
1085 DimVector normal = intersection.centerUnitOuterNormal();
1086 interactionVolume.setNormal(normal, 4, 1);
1087 interactionVolume.setNormal(normal, 6, 2);
1088 globalPosFace = intersection.geometry().center();
1089 interactionVolume.setFacePosition(globalPosFace, 4);
1090 interactionVolume.setFacePosition(globalPosFace, 6);
1091 int indexInOutside = intersection.indexInOutside();
1092 interactionVolume.setIndexOnElement(indexInOutside, 5, 2);
1093 interactionVolume.setIndexOnElement(indexInOutside, 7, 1);
1095 interactionVolume.setNormal(normal, 5, 2);
1096 interactionVolume.setNormal(normal, 7, 1);
1102 DimVector edgeCoord2(interactionVolume.getFacePosition(4));
1104 globalPosFace = interactionVolume.getFacePosition(4);
1105 crossProductVector1 = centerPos - globalPosFace;
1106 crossProductVector2 = edgeCoord2 - edgeCoord3;
1107 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1108 interactionVolume.setFaceArea(faceArea, 4);
1110 globalPosFace = interactionVolume.getFacePosition(6);
1111 crossProductVector1 = centerPos - globalPosFace;
1112 crossProductVector2 = edgeCoord2 - edgeCoord5;
1113 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1114 interactionVolume.setFaceArea(faceArea, 6);
1119 const ElementGeometry& geometry = element5.geometry();
1121 const auto referenceElement = ReferenceElements::general(geometry.type());
1123 int oldSubVolumElemIdx = IndexTranslator::getOldElemIdxFromNewFaceIdxto0(zeroFaceIdx, 4);
1124 int oldEdgeIdx = IndexTranslator::getOldEdgeIdxFromNewFaceIdxto0(zeroFaceIdx, 1);
1126 DimVector edgeCoord(0.0);
1127 switch (oldSubVolumElemIdx)
1134 edgeCoord = geometry.global(referenceElement.position(9, dim - 1));
1137 edgeCoord = geometry.global(referenceElement.position(3, dim - 1));
1140 edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
1151 edgeCoord = geometry.global(referenceElement.position(2, dim - 1));
1154 edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
1157 edgeCoord = geometry.global(referenceElement.position(11, dim - 1));
1168 edgeCoord = geometry.global(referenceElement.position(1, dim - 1));
1171 edgeCoord = geometry.global(referenceElement.position(9, dim - 1));
1174 edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
1185 edgeCoord = geometry.global(referenceElement.position(0, dim - 1));
1188 edgeCoord = geometry.global(referenceElement.position(8, dim - 1));
1191 edgeCoord = geometry.global(referenceElement.position(10, dim - 1));
1202 edgeCoord = geometry.global(referenceElement.position(3, dim - 1));
1205 edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
1208 edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
1219 edgeCoord = geometry.global(referenceElement.position(2, dim - 1));
1222 edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
1225 edgeCoord = geometry.global(referenceElement.position(7, dim - 1));
1236 edgeCoord = geometry.global(referenceElement.position(1, dim - 1));
1239 edgeCoord = geometry.global(referenceElement.position(5, dim - 1));
1242 edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
1253 edgeCoord = geometry.global(referenceElement.position(0, dim - 1));
1256 edgeCoord = geometry.global(referenceElement.position(4, dim - 1));
1259 edgeCoord = geometry.global(referenceElement.position(6, dim - 1));
1267 interactionVolume.setEdgePosition(edgeCoord, 1);
1269 else if (interactionVolume.getHangingNodeType() == InteractionVolume::fourSmallCellsDiag)
1271 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
1272 DimVector edgeCoord2(interactionVolume.getEdgePosition(1));
1273 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
1274 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
1275 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
1276 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
1278 DimVector crossProductVector1(0);
1279 DimVector crossProductVector2(0);
1281 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
1282 crossProductVector1 = centerPos - globalPosFace;
1283 crossProductVector2 = edgeCoord1 - edgeCoord3;
1284 Scalar faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1285 interactionVolume.setFaceArea(faceArea, 0);
1287 globalPosFace = interactionVolume.getFacePosition(1);
1288 crossProductVector1 = centerPos - globalPosFace;
1289 crossProductVector2 = edgeCoord1 - edgeCoord4;
1290 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1291 interactionVolume.setFaceArea(faceArea, 1);
1293 globalPosFace = interactionVolume.getFacePosition(3);
1294 crossProductVector1 = centerPos - globalPosFace;
1295 crossProductVector2 = edgeCoord1 - edgeCoord6;
1296 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1297 interactionVolume.setFaceArea(faceArea, 3);
1299 globalPosFace = interactionVolume.getFacePosition(5);
1300 crossProductVector1 = centerPos - globalPosFace;
1301 crossProductVector2 = edgeCoord2 - edgeCoord4;
1302 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1303 interactionVolume.setFaceArea(faceArea, 5);
1305 globalPosFace = interactionVolume.getFacePosition(6);
1306 crossProductVector1 = centerPos - globalPosFace;
1307 crossProductVector2 = edgeCoord2 - edgeCoord5;
1308 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1309 interactionVolume.setFaceArea(faceArea, 6);
1311 globalPosFace = interactionVolume.getFacePosition(7);
1312 crossProductVector1 = centerPos - globalPosFace;
1313 crossProductVector2 = edgeCoord2 - edgeCoord6;
1314 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1315 interactionVolume.setFaceArea(faceArea, 7);
1317 globalPosFace = interactionVolume.getFacePosition(8);
1318 crossProductVector1 = centerPos - globalPosFace;
1319 crossProductVector2 = edgeCoord3 - edgeCoord6;
1320 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1321 interactionVolume.setFaceArea(faceArea, 8);
1323 globalPosFace = interactionVolume.getFacePosition(9);
1324 crossProductVector1 = centerPos - globalPosFace;
1325 crossProductVector2 = edgeCoord3 - edgeCoord4;
1326 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1327 interactionVolume.setFaceArea(faceArea, 9);
1329 globalPosFace = interactionVolume.getFacePosition(10);
1330 crossProductVector1 = centerPos - globalPosFace;
1331 crossProductVector2 = edgeCoord4 - edgeCoord5;
1332 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1333 interactionVolume.setFaceArea(faceArea, 10);
1335 globalPosFace = interactionVolume.getFacePosition(11);
1336 crossProductVector1 = centerPos - globalPosFace;
1337 crossProductVector2 = edgeCoord5 - edgeCoord6;
1338 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1339 interactionVolume.setFaceArea(faceArea, 11);
1341 auto element1 = interactionVolume.getSubVolumeElement(0);
1343 bool hasFaceOne =
false;
1344 bool hasFaceTwo =
false;
1345 for (
const auto& intersection : intersections(problem_.gridView(), element1))
1347 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(0, 1))
1349 if (intersection.neighbor())
1351 auto outside = intersection.outside();
1352 if (element1.level() > outside.level())
1354 interactionVolume.setSubVolumeElement(outside, 2);
1355 interactionVolume.setSubVolumeElement(outside, 3);
1363 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
1365 if (intersection.neighbor())
1367 auto outside = intersection.outside();
1368 if (element1.level() > outside.level())
1370 interactionVolume.setSubVolumeElement(outside, 4);
1371 interactionVolume.setSubVolumeElement(outside, 5);
1387 InteractionVolume hangingNodeVolume(problem_.gridView().grid());
1389 std::vector<int> elemIdxOld;
1390 for (
int i = 0; i < InteractionVolume::subVolumeTotalNum; i++)
1392 if (interactionVolume.hasSubVolumeElement(i))
1393 elemIdxOld.push_back(i);
1396 std::set<int> zeroFaceIdxVec;
1397 for (
int i = 0; i < 6; i++)
1399 for (
int j = 0; j < 6; j++)
1401 int fIdx = IndexTranslator::getFaceIndexFromElements(elemIdxOld[i], elemIdxOld[j]);
1403 zeroFaceIdxVec.insert(fIdx);
1407 std::vector<int> elemIdxNew(6);
1408 int zeroFaceIdx = 0;
1411 std::set<int>::iterator it = zeroFaceIdxVec.begin();
1412 for (; it != zeroFaceIdxVec.end(); ++it)
1417 for (
int i = 0; i < 6; i++)
1419 elemIdxNew[i] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elemIdxOld[i]);
1420 if (elemIdxNew[i] == 6 || elemIdxNew[i] == 7)
1428 for (
int i = 0; i < 6; i++)
1430 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[i]),
1438 for (
int elem = 0; elem < InteractionVolume::subVolumeTotalNum; elem++)
1440 for (
int i = 0; i < 3; i++)
1442 int faceIdxOld = IndexTranslator::getFaceIndexFromSubVolume(elem, i);
1443 int faceIdxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, faceIdxOld);
1445 for (
int j = 0; j < 3; j++)
1447 int elemNew = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elem);
1448 int faceIdxNewTest = IndexTranslator::getFaceIndexFromSubVolume(elemNew, j);
1450 if (faceIdxNew == faceIdxNewTest)
1452 hangingNodeVolume.setNormal(interactionVolume.getNormal(elem, i),
1454 hangingNodeVolume.setIndexOnElement(
1455 interactionVolume.getIndexOnElement(elem, i), elemNew, j);
1461 for (
int i = 0; i < InteractionVolume::fluxFacesTotalNum; i++)
1463 int idxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, i);
1464 hangingNodeVolume.setFaceArea(interactionVolume.getFaceArea(i), idxNew);
1465 hangingNodeVolume.setFacePosition(interactionVolume.getFacePosition(i), idxNew);
1468 for (
int i = 0; i < InteractionVolume::fluxEdgesTotalNum; i++)
1470 int idxNew = IndexTranslator::getNewEdgeIdxFromOldFaceIdxto0(zeroFaceIdx, i);
1471 hangingNodeVolume.setEdgePosition(interactionVolume.getEdgePosition(i), idxNew);
1474 interactionVolume = hangingNodeVolume;
1476 interactionVolume.setCenterPosition(centerPos);
1478 interactionVolume.setHangingNodeType(InteractionVolume::sixSmallCells);
1480 auto element3 = interactionVolume.getSubVolumeElement(2);
1482 for (
const auto& intersection : intersections(problem_.gridView(), element3))
1484 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(2, 2))
1486 if (intersection.neighbor())
1488 auto outside = intersection.outside();
1489 if (element3.level() > outside.level())
1491 interactionVolume.setSubVolumeElement(outside, 6);
1492 interactionVolume.setSubVolumeElement(outside, 7);
1500 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
1501 DimVector edgeCoord2(interactionVolume.getEdgePosition(1));
1502 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
1503 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
1504 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
1505 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
1507 DimVector crossProductVector1(0);
1508 DimVector crossProductVector2(0);
1510 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
1511 crossProductVector1 = centerPos - globalPosFace;
1512 crossProductVector2 = edgeCoord1 - edgeCoord3;
1513 Scalar faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1514 interactionVolume.setFaceArea(faceArea, 0);
1516 globalPosFace = interactionVolume.getFacePosition(1);
1517 crossProductVector1 = centerPos - globalPosFace;
1518 crossProductVector2 = edgeCoord1 - edgeCoord4;
1519 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1520 interactionVolume.setFaceArea(faceArea, 1);
1522 globalPosFace = interactionVolume.getFacePosition(2);
1523 crossProductVector1 = centerPos - globalPosFace;
1524 crossProductVector2 = edgeCoord1 - edgeCoord5;
1525 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1526 interactionVolume.setFaceArea(faceArea, 2);
1528 globalPosFace = interactionVolume.getFacePosition(3);
1529 crossProductVector1 = centerPos - globalPosFace;
1530 crossProductVector2 = edgeCoord1 - edgeCoord6;
1531 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1532 interactionVolume.setFaceArea(faceArea, 3);
1534 globalPosFace = interactionVolume.getFacePosition(4);
1535 crossProductVector1 = centerPos - globalPosFace;
1536 crossProductVector2 = edgeCoord2 - edgeCoord3;
1537 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1538 interactionVolume.setFaceArea(faceArea, 4);
1540 globalPosFace = interactionVolume.getFacePosition(5);
1541 crossProductVector1 = centerPos - globalPosFace;
1542 crossProductVector2 = edgeCoord2 - edgeCoord4;
1543 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1544 interactionVolume.setFaceArea(faceArea, 5);
1546 globalPosFace = interactionVolume.getFacePosition(7);
1547 crossProductVector1 = centerPos - globalPosFace;
1548 crossProductVector2 = edgeCoord2 - edgeCoord6;
1549 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1550 interactionVolume.setFaceArea(faceArea, 7);
1552 globalPosFace = interactionVolume.getFacePosition(8);
1553 crossProductVector1 = centerPos - globalPosFace;
1554 crossProductVector2 = edgeCoord3 - edgeCoord6;
1555 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1556 interactionVolume.setFaceArea(faceArea, 8);
1558 globalPosFace = interactionVolume.getFacePosition(9);
1559 crossProductVector1 = centerPos - globalPosFace;
1560 crossProductVector2 = edgeCoord3 - edgeCoord4;
1561 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1562 interactionVolume.setFaceArea(faceArea, 9);
1564 globalPosFace = interactionVolume.getFacePosition(10);
1565 crossProductVector1 = centerPos - globalPosFace;
1566 crossProductVector2 = edgeCoord4 - edgeCoord5;
1567 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1568 interactionVolume.setFaceArea(faceArea, 10);
1570 globalPosFace = interactionVolume.getFacePosition(11);
1571 crossProductVector1 = centerPos - globalPosFace;
1572 crossProductVector2 = edgeCoord5 - edgeCoord6;
1573 faceArea =
crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1574 interactionVolume.setFaceArea(faceArea, 11);
1579 interactionVolume.printInteractionVolumeInfo();
1580 DUNE_THROW(Dune::NotImplemented,
"Hanging node shape not implemented");
1587template<
class TypeTag>
1590 std::vector < std::vector<int> > elemVertMap(problem_.gridView().size(dim), std::vector<int>(8, -1));
1593 for (
const auto& element : elements(problem_.gridView()))
1594 asImp_().storeSubVolumeElements(element, elemVertMap);
1596 for (
unsigned int i = 0; i < asImp_().interactionVolumes_.size(); i++)
1597 if (asImp_().interactionVolumes_[i].getElementNumber() == 0)
1598 asImp_().interactionVolumes_[i].printInteractionVolumeInfo();
1601 for (
const auto& element : elements(problem_.gridView()))
1602 asImp_().storeIntersectionInfo(element, elemVertMap);
1604 faceVertices_.clear();
1605 faceVertices_.resize(problem_.gridView().size(0));
1609 for (
const auto& vertex : vertices(problem_.gridView()))
1611 int vIdxGlobal = problem_.variables().index(vertex);
1613 InteractionVolume& interactionVolume = asImp_().interactionVolumes_[vIdxGlobal];
1615 if (interactionVolume.getElementNumber() == 8)
1617 asImp_().storeInnerInteractionVolume(interactionVolume, vertex);
1619 else if (interactionVolume.isBoundaryInteractionVolume())
1621 asImp_().storeBoundaryInteractionVolume(interactionVolume, vertex);
1626 storeHangingNodeInteractionVolume(interactionVolume, vertex);
1629 if (!interactionVolume.isBoundaryInteractionVolume())
1631 auto element1 = interactionVolume.getSubVolumeElement(0);
1632 auto element2 = interactionVolume.getSubVolumeElement(1);
1633 auto element3 = interactionVolume.getSubVolumeElement(2);
1634 auto element4 = interactionVolume.getSubVolumeElement(3);
1635 auto element5 = interactionVolume.getSubVolumeElement(4);
1636 auto element6 = interactionVolume.getSubVolumeElement(5);
1637 auto element7 = interactionVolume.getSubVolumeElement(6);
1638 auto element8 = interactionVolume.getSubVolumeElement(7);
1640 int globalIdx1 = problem_.variables().index(element1);
1641 int globalIdx2 = problem_.variables().index(element2);
1642 int globalIdx3 = problem_.variables().index(element3);
1643 int globalIdx4 = problem_.variables().index(element4);
1644 int globalIdx5 = problem_.variables().index(element5);
1645 int globalIdx6 = problem_.variables().index(element6);
1646 int globalIdx7 = problem_.variables().index(element7);
1647 int globalIdx8 = problem_.variables().index(element8);
1649 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(0, 0), globalIdx1, interactionVolume.getIndexOnElement(0, 0));
1650 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(0, 1), globalIdx1, interactionVolume.getIndexOnElement(0, 1));
1651 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(0, 2), globalIdx1, interactionVolume.getIndexOnElement(0, 2));
1652 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(1, 0), globalIdx2, interactionVolume.getIndexOnElement(1, 0));
1653 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(1, 1), globalIdx2, interactionVolume.getIndexOnElement(1, 1));
1654 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(1, 2), globalIdx2, interactionVolume.getIndexOnElement(1, 2));
1655 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(2, 0), globalIdx3, interactionVolume.getIndexOnElement(2, 0));
1656 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(3, 1), globalIdx4, interactionVolume.getIndexOnElement(3, 1));
1657 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(4, 0), globalIdx5, interactionVolume.getIndexOnElement(4, 0));
1658 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(5, 0), globalIdx6, interactionVolume.getIndexOnElement(5, 0));
1660 faceVertices_[globalIdx1][interactionVolume.getIndexOnElement(0, 0)].insert(vIdxGlobal);
1661 faceVertices_[globalIdx1][interactionVolume.getIndexOnElement(0, 1)].insert(vIdxGlobal);
1662 faceVertices_[globalIdx1][interactionVolume.getIndexOnElement(0, 2)].insert(vIdxGlobal);
1663 faceVertices_[globalIdx2][interactionVolume.getIndexOnElement(1, 0)].insert(vIdxGlobal);
1664 faceVertices_[globalIdx2][interactionVolume.getIndexOnElement(1, 1)].insert(vIdxGlobal);
1665 faceVertices_[globalIdx2][interactionVolume.getIndexOnElement(1, 2)].insert(vIdxGlobal);
1666 faceVertices_[globalIdx3][interactionVolume.getIndexOnElement(2, 0)].insert(vIdxGlobal);
1667 faceVertices_[globalIdx4][interactionVolume.getIndexOnElement(3, 1)].insert(vIdxGlobal);
1668 faceVertices_[globalIdx5][interactionVolume.getIndexOnElement(4, 0)].insert(vIdxGlobal);
1669 faceVertices_[globalIdx6][interactionVolume.getIndexOnElement(5, 0)].insert(vIdxGlobal);
1671 if (interactionVolume.getHangingNodeType() != InteractionVolume::twoSmallCells)
1673 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(2, 1), globalIdx3, interactionVolume.getIndexOnElement(2, 1));
1674 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(2, 2), globalIdx3, interactionVolume.getIndexOnElement(2, 2));
1675 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(3, 0), globalIdx4, interactionVolume.getIndexOnElement(3, 0));
1676 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(3, 2), globalIdx4, interactionVolume.getIndexOnElement(3, 2));
1677 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(4, 1), globalIdx5, interactionVolume.getIndexOnElement(4, 1));
1678 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(4, 2), globalIdx5, interactionVolume.getIndexOnElement(4, 2));
1679 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(5, 1), globalIdx6, interactionVolume.getIndexOnElement(5, 1));
1680 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(5, 2), globalIdx6, interactionVolume.getIndexOnElement(5, 2));
1681 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(6, 0), globalIdx7, interactionVolume.getIndexOnElement(6, 0));
1682 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(6, 1), globalIdx7, interactionVolume.getIndexOnElement(6, 1));
1683 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(6, 2), globalIdx7, interactionVolume.getIndexOnElement(6, 2));
1684 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(7, 0), globalIdx8, interactionVolume.getIndexOnElement(7, 0));
1685 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(7, 1), globalIdx8, interactionVolume.getIndexOnElement(7, 1));
1686 asImp_().addRealFluxFaceArea_(interactionVolume.getFaceArea(7, 2), globalIdx8, interactionVolume.getIndexOnElement(7, 2));
1688 faceVertices_[globalIdx3][interactionVolume.getIndexOnElement(2, 1)].insert(vIdxGlobal);
1689 faceVertices_[globalIdx3][interactionVolume.getIndexOnElement(2, 2)].insert(vIdxGlobal);
1690 faceVertices_[globalIdx4][interactionVolume.getIndexOnElement(3, 0)].insert(vIdxGlobal);
1691 faceVertices_[globalIdx4][interactionVolume.getIndexOnElement(3, 2)].insert(vIdxGlobal);
1692 faceVertices_[globalIdx5][interactionVolume.getIndexOnElement(4, 1)].insert(vIdxGlobal);
1693 faceVertices_[globalIdx5][interactionVolume.getIndexOnElement(4, 2)].insert(vIdxGlobal);
1694 faceVertices_[globalIdx6][interactionVolume.getIndexOnElement(5, 1)].insert(vIdxGlobal);
1695 faceVertices_[globalIdx6][interactionVolume.getIndexOnElement(5, 2)].insert(vIdxGlobal);
1696 faceVertices_[globalIdx7][interactionVolume.getIndexOnElement(6, 0)].insert(vIdxGlobal);
1697 faceVertices_[globalIdx7][interactionVolume.getIndexOnElement(6, 1)].insert(vIdxGlobal);
1698 faceVertices_[globalIdx7][interactionVolume.getIndexOnElement(6, 2)].insert(vIdxGlobal);
1699 faceVertices_[globalIdx8][interactionVolume.getIndexOnElement(7, 0)].insert(vIdxGlobal);
1700 faceVertices_[globalIdx8][interactionVolume.getIndexOnElement(7, 1)].insert(vIdxGlobal);
1701 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...
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
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:94
InteractionVolume & interactionVolume(int vertexIdx)
Returns an interaction volume.
Definition: 3dinteractionvolumecontainer.hh:143
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:1588
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.