3.4
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
3dinteractionvolumecontaineradaptive.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_FVMPFAL3D_INTERACTIONVOLUMECONTAINER_ADAPTIVE_HH
25#define DUMUX_FVMPFAL3D_INTERACTIONVOLUMECONTAINER_ADAPTIVE_HH
26
27// dumux environment
32
33namespace Dumux {
34
45template<class TypeTag>
47{
48 friend class FvMpfaL3dInteractionVolumeContainer<TypeTag>;
51
52 enum
53 {
54 dim = GridView::dimension, dimWorld = GridView::dimensionworld
55 };
56
59
61
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;
65
66 using GlobalPosition = typename ElementGeometry::GlobalCoordinate;
67 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
68
69 using DimVector = Dune::FieldVector<Scalar, dim>;
70 using FaceVerticesVector = std::vector<Dune::FieldVector<std::set<int>, 2*dim> >;
71
72 enum
73 {
74 pressureEqIdx = Indices::pressureEqIdx,
75 };
76
77 using IndexTranslator = IndexTranslatorAdaptive;
78
79public:
83
84private:
85
86 void storeHangingNodeInteractionVolume(InteractionVolume& interactionVolume, const Vertex& vertex);
87 void storeInnerInteractionVolume(InteractionVolume& interactionVolume, const Vertex& vertex);
88
89protected:
91
92public:
93
102 std::set<int>& faceVerticeIndices(int eIdxGlobal, int fIdx)
103 {
104 return faceVertices_[eIdxGlobal][fIdx];
105 }
106
112 ParentType(problem), problem_(problem)
113 {
114 if (dim != 3)
115 {
116 DUNE_THROW(Dune::NotImplemented, "Dimension not supported!");
117 }
118 }
119
120private:
122
124 Implementation &asImp_()
125 { return *static_cast<Implementation *>(this);}
126
128 const Implementation &asImp_() const
129 { return *static_cast<const Implementation *>(this);}
130
131 Problem& problem_;
132 FaceVerticesVector faceVertices_;
133};
134
158template<class TypeTag>
159void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeInnerInteractionVolume(InteractionVolume& interactionVolume,
160 const Vertex& vertex)
161{
162 // if cells of different grid level appear, the geometric information is constructed going from the
163 // coarsest level to the finest level. This ensures that coarser information ins always overwritten by
164 // finer information.
165 if (!interactionVolume.sameLevel())
166 {
167 // sort the local element indices according to the grid level
168 std::vector < std::vector<int> > levelIdx(8, std::vector<int>(2));
169 for (int i = 0; i < 8; i++)
170 {
171 levelIdx[i][0] = i;
172 levelIdx[i][1] = interactionVolume.getSubVolumeElement(i).level();
173 }
174
175 std::sort(levelIdx.begin(), levelIdx.end(), [](const auto& a, const auto& b) { return (a[1]<b[1]); });
176
177 // Generate and store the geometric information going from the coarsest to the finest level.
178 // For the calculation we take advantage from the fact that the ordering inside the interaction volume
179 // with respect to the DUNE reference element is known due to the storage process of the elements in
180 // FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeSubVolumeElements
181 // (const Element& element, std::vector < std::vector<int> >& elemVertMap)
182 for (int i = 0; i < 8; i++)
183 {
184 int idx = levelIdx[i][0];
185 auto element = interactionVolume.getSubVolumeElement(idx);
186
187 const ElementGeometry& geometry = element.geometry();
188
189 const auto refElement = referenceElement(geometry);
190
191 switch (idx)
192 {
193 case 0:
194 {
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);
201
202 break;
203 }
204 case 1:
205 {
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);
212
213 break;
214 }
215 case 2:
216 {
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);
223
224 break;
225 }
226 case 3:
227 {
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);
234
235 break;
236 }
237 case 4:
238 {
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);
245
246 break;
247 }
248 case 5:
249 {
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);
256
257 break;
258 }
259 case 6:
260 {
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);
267
268 break;
269 }
270 case 7:
271 {
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);
278
279 break;
280 }
281 }
282 }
283 ParentType::storeInnerInteractionVolume(interactionVolume, vertex, false);
284 }
285 else
286 {
287 ParentType::storeInnerInteractionVolume(interactionVolume, vertex, true);
288 }
289}
290
313template<class TypeTag>
314void FvMpfaL3dInteractionVolumeContainerAdaptive<TypeTag>::storeHangingNodeInteractionVolume(InteractionVolume& interactionVolume,
315 const Vertex& vertex)
316{
317 const DimVector& centerPos = vertex.geometry().center();
318
319 interactionVolume.setCenterPosition(centerPos);
320
321 // sort the local element indices according to the grid level
322 std::vector < std::vector<int> > levelIdx(8, std::vector<int>(2));
323 for (int i = 0; i < 8; i++)
324 {
325 levelIdx[i][0] = i;
326 if (interactionVolume.hasSubVolumeElement(i))
327 levelIdx[i][1] = interactionVolume.getSubVolumeElement(i).level();
328 else
329 levelIdx[i][1] = -1;
330 }
331
332 std::sort(levelIdx.begin(), levelIdx.end(), [](const auto& a, const auto& b) { return (a[1]<b[1]); });
333
334 // Generate and store the geometric information going from the coarsest to the finest level.
335 // For the calculation we take advantage from the fact that the ordering inside the interaction volume
336 // with respect to the DUNE reference element is known due to the storage process of the elements in
337 // FvMpfaL3dInteractionVolumeContainer<TypeTag>::
338 // storeSubVolumeElements(const Element& element, std::vector < std::vector<int> >& elemVertMap)
339 for (int i = 0; i < 8; i++)
340 {
341 if (levelIdx[i][1] < 0)
342 continue;
343
344 int idx = levelIdx[i][0];
345
346 auto element = interactionVolume.getSubVolumeElement(idx);
347
348 const ElementGeometry& geometry = element.geometry();
349
350 const auto refElement = referenceElement(geometry);
351
352 switch (idx)
353 {
354 case 0:
355 {
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);
362
363 break;
364 }
365 case 1:
366 {
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);
373
374 break;
375 }
376 case 2:
377 {
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);
384
385 break;
386 }
387 case 3:
388 {
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);
395
396 break;
397 }
398 case 4:
399 {
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);
406
407 break;
408 }
409 case 5:
410 {
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);
417
418 break;
419 }
420 case 6:
421 {
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);
428
429 break;
430 }
431 case 7:
432 {
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);
439
440 break;
441 }
442 }
443 }
444
445 // Choose the type of the hanging-node-interaction-volume depending on the number of stored fine
446 // elements (see dissertation M. Wolff, http://elib.uni-stuttgart.de/opus/volltexte/2013/8661/)
447 // and add missing elements and geometric information
448 switch (interactionVolume.getElementNumber())
449 {
450 // hanging-node interaction volume of type 5 or 7
451 case 2:
452 {
453 InteractionVolume hangingNodeVolume(problem_.gridView().grid());
454
455 std::vector<int> elemIdxOld;
456 for (int i = 0; i < InteractionVolume::subVolumeTotalNum; i++)
457 {
458 if (interactionVolume.hasSubVolumeElement(i))
459 elemIdxOld.push_back(i);
460 }
461
462 int zeroFaceIdx = IndexTranslator::getFaceIndexFromElements(elemIdxOld[0], elemIdxOld[1]);
463
464 std::vector<int> elemIdxNew(2);
465 elemIdxNew[0] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elemIdxOld[0]);
466 elemIdxNew[1] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elemIdxOld[1]);
467
468
469 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[0]),
470 elemIdxNew[0]);
471 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[1]),
472 elemIdxNew[1]);
473
474 for (int elem = 0; elem < InteractionVolume::subVolumeTotalNum; elem++)
475 {
476 for (int i = 0; i < 3; i++)
477 {
478 int faceIdxOld = IndexTranslator::getFaceIndexFromSubVolume(elem, i);
479 int faceIdxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, faceIdxOld);
480
481 for (int j = 0; j < 3; j++)
482 {
483 int elemNew = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elem);
484 int faceIdxNewTest = IndexTranslator::getFaceIndexFromSubVolume(elemNew, j);
485
486 if (faceIdxNew == faceIdxNewTest)
487 {
488 hangingNodeVolume.setNormal(interactionVolume.getNormal(elem, i),
489 elemNew, j);
490 hangingNodeVolume.setIndexOnElement(
491 interactionVolume.getIndexOnElement(elem, i), elemNew, j);
492 }
493 }
494 }
495 }
496
497 for (int i = 0; i < InteractionVolume::fluxFacesTotalNum; i++)
498 {
499 int idxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, i);
500 hangingNodeVolume.setFacePosition(interactionVolume.getFacePosition(i), idxNew);
501 }
502
503 for (int i = 0; i < InteractionVolume::fluxEdgesTotalNum; i++)
504 {
505 int idxNew = IndexTranslator::getNewEdgeIdxFromOldFaceIdxto0(zeroFaceIdx, i);
506 hangingNodeVolume.setEdgePosition(interactionVolume.getEdgePosition(i), idxNew);
507 }
508
509 interactionVolume = hangingNodeVolume;
510
511 interactionVolume.setHangingNodeType(InteractionVolume::twoSmallCells);
512
513 interactionVolume.setCenterPosition(centerPos);
514
515 auto element1 = interactionVolume.getSubVolumeElement(0);
516
517 for (const auto& intersection : intersections(problem_.gridView(), element1))
518 {
519 int idxInInside = intersection.indexInInside();
520
521 if (idxInInside == interactionVolume.getIndexOnElement(0, 2))
522 {
523 if (intersection.neighbor())
524 {
525 auto outside = intersection.outside();
526 if (element1.level() > outside.level())
527 {
528 interactionVolume.setSubVolumeElement(outside, 4);
529 interactionVolume.setSubVolumeElement(outside, 5);
530 }
531 }
532 }
533 else if (idxInInside == interactionVolume.getIndexOnElement(0, 1))
534 {
535 if (intersection.neighbor())
536 {
537 auto outside = intersection.outside();
538 if (element1.level() > outside.level())
539 {
540 interactionVolume.setSubVolumeElement(outside, 2);
541 interactionVolume.setSubVolumeElement(outside, 3);
542 }
543 }
544 }
545 }
546 auto element2 = interactionVolume.getSubVolumeElement(1);
547 auto element45 = interactionVolume.getSubVolumeElement(4);
548 auto element23 = interactionVolume.getSubVolumeElement(2);
549
550 for (const auto& intersection1 : intersections(problem_.gridView(), element45))
551 {
552 if (intersection1.neighbor())
553 {
554 auto element45Outside = intersection1.outside();
555
556 for (const auto& intersection2 : intersections(problem_.gridView(), element23))
557 {
558 if (intersection2.neighbor())
559 {
560 auto element23Outside = intersection2.outside();
561
562 if (element45Outside == element23Outside && element45Outside != element1
563 && element45Outside != element2)
564 {
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);
570 normal *= -1;
571 interactionVolume.setNormal(normal, 6, 0);
572 interactionVolume.setNormal(normal, 7, 0);
573
574 GlobalPosition globalPosFace = intersection2.geometry().center();
575 interactionVolume.setFacePosition(globalPosFace, 10);
576 interactionVolume.setFacePosition(globalPosFace, 11);
577 interactionVolume.setEdgePosition(centerPos, 4);
578
579 Scalar faceArea = intersection2.geometry().volume()/4.0;
580
581 interactionVolume.setFaceArea(faceArea, 10);
582 interactionVolume.setFaceArea(faceArea, 11);
583
584 normal = intersection1.centerUnitOuterNormal();
585 interactionVolume.setNormal(normal, 4, 2);
586 interactionVolume.setNormal(normal, 5, 1);
587 normal *= -1;
588 interactionVolume.setNormal(normal, 6, 1);
589 interactionVolume.setNormal(normal, 7, 2);
590
591 globalPosFace = intersection1.geometry().center();
592 interactionVolume.setFacePosition(globalPosFace, 5);
593 interactionVolume.setFacePosition(globalPosFace, 7);
594 interactionVolume.setEdgePosition(centerPos, 1);
595
596 faceArea = intersection1.geometry().volume()/4.0;
597
598 interactionVolume.setFaceArea(faceArea, 5);
599 interactionVolume.setFaceArea(faceArea, 7);
600 }
601 }
602 }
603 }
604 }
605
606 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
607 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
608 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
609 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
610
611 DimVector crossProductVector1(0);
612 DimVector crossProductVector2(0);
613
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);
619
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);
625
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);
631
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);
637
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);
643
644 break;
645 }
646 //hanging-node interaction volume of type 1, 3 or 4
647 case 4:
648 {
649 InteractionVolume hangingNodeVolume(problem_.gridView().grid());
650
651 std::vector<int> elemIdxOld;
652 for (int i = 0; i < InteractionVolume::subVolumeTotalNum; i++)
653 {
654 if (interactionVolume.hasSubVolumeElement(i))
655 elemIdxOld.push_back(i);
656 }
657
658 std::set<int> zeroFaceIdxVec;
659 for (int i = 0; i < 4; i++)
660 {
661 for (int j = 0; j < 4; j++)
662 {
663 int fIdx = IndexTranslator::getFaceIndexFromElements(elemIdxOld[i], elemIdxOld[j]);
664 if (fIdx >= 0)
665 zeroFaceIdxVec.insert(fIdx);
666 }
667 }
668
669 std::vector<int> elemIdxNew(4);
670 int zeroFaceIdx = 0;
671
672 if (zeroFaceIdxVec.size() == 2)
673 {
674 hangingNodeVolume.setHangingNodeType(InteractionVolume::fourSmallCellsDiag);
675
676 if (zeroFaceIdxVec.find(0) == zeroFaceIdxVec.end())
677 zeroFaceIdx = *zeroFaceIdxVec.begin();
678
679 for (int i = 0; i < 4; i++)
680 {
681 elemIdxNew[i] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elemIdxOld[i]);
682 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[i]),
683 elemIdxNew[i]);
684 }
685 }
686 else if (zeroFaceIdxVec.size() == 4)
687 {
688 std::set<int>::iterator it = zeroFaceIdxVec.begin();
689 for (; it != zeroFaceIdxVec.end(); ++it)
690 {
691 zeroFaceIdx = *it;
692
693 bool isFace = true;
694 for (int i = 0; i < 4; i++)
695 {
696 elemIdxNew[i] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx,
697 elemIdxOld[i]);
698 if (elemIdxNew[i] == 4 || elemIdxNew[i] == 5 || elemIdxNew[i] == 6 || elemIdxNew[i] == 7)
699 {
700 isFace = false;
701 break;
702 }
703 }
704 if (isFace)
705 {
706 for (int i = 0; i < 4; i++)
707 {
708 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[i]), elemIdxNew[i]);
709 }
710
711 break;
712 }
713 }
714 }
715 else
716 {
717 std::set<int>::iterator it = zeroFaceIdxVec.begin();
718 std::cout<<"zeroFaceIdxVec = ";
719 for(;it != zeroFaceIdxVec.end(); ++it)
720 {
721 std::cout<<" "<<*it;
722 }
723 std::cout<<"\n";
724 }
725
726 for (int elem = 0; elem < InteractionVolume::subVolumeTotalNum; elem++)
727 {
728 for (int i = 0; i < 3; i++)
729 {
730 int faceIdxOld = IndexTranslator::getFaceIndexFromSubVolume(elem, i);
731 int faceIdxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, faceIdxOld);
732
733 for (int j = 0; j < 3; j++)
734 {
735 int elemNew = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elem);
736 int faceIdxNewTest = IndexTranslator::getFaceIndexFromSubVolume(elemNew, j);
737
738 if (faceIdxNew == faceIdxNewTest)
739 {
740 hangingNodeVolume.setNormal(interactionVolume.getNormal(elem, i),
741 elemNew, j);
742 hangingNodeVolume.setIndexOnElement(interactionVolume.getIndexOnElement(elem, i), elemNew, j);
743 }
744 }
745 }
746 }
747
748 for (int i = 0; i < InteractionVolume::fluxFacesTotalNum; i++)
749 {
750 int idxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, i);
751 hangingNodeVolume.setFacePosition(interactionVolume.getFacePosition(i), idxNew);
752 }
753
754 for (int i = 0; i < InteractionVolume::fluxEdgesTotalNum; i++)
755 {
756 int idxNew = IndexTranslator::getNewEdgeIdxFromOldFaceIdxto0(zeroFaceIdx, i);
757 hangingNodeVolume.setEdgePosition(interactionVolume.getEdgePosition(i), idxNew);
758 }
759
760 interactionVolume = hangingNodeVolume;
761
762 interactionVolume.setCenterPosition(centerPos);
763
764 if (zeroFaceIdxVec.size() == 4)
765 {
766 auto element1 = interactionVolume.getSubVolumeElement(0);
767 auto element4 = interactionVolume.getSubVolumeElement(3);
768
769 auto outside1 = element1;
770 auto outside4 = element4;
771
772 for (const auto& intersection1 : intersections(problem_.gridView(), element1))
773 {
774 if (intersection1.neighbor())
775 {
776 if (intersection1.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
777 {
778 outside1 = intersection1.outside();
779 break;
780 }
781 }
782 }
783 for (const auto& intersection4 : intersections(problem_.gridView(), element4))
784 {
785 if (intersection4.neighbor())
786 {
787 if (intersection4.indexInInside() == interactionVolume.getIndexOnElement(3, 2))
788 {
789 outside4 = intersection4.outside();
790 break;
791 }
792 }
793 }
794 if (outside1 != outside4)
795 {
796 interactionVolume.setHangingNodeType(InteractionVolume::fourSmallCellsEdge);
797 }
798 else if (outside1 == outside4)
799 {
800 interactionVolume.setHangingNodeType(InteractionVolume::fourSmallCellsFace);
801 }
802 }
803
804 if (interactionVolume.getHangingNodeType() == InteractionVolume::fourSmallCellsFace)
805 {
806 auto element = interactionVolume.getSubVolumeElement(0);
807
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));
813
814 DimVector crossProductVector1(0);
815 DimVector crossProductVector2(0);
816
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);
822
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);
828
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);
834
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);
840
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);
846
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);
852
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);
858
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);
864
865 for (const auto& intersection : intersections(problem_.gridView(), element))
866 {
867 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
868 {
869 if (intersection.neighbor())
870 {
871 auto outside = intersection.outside();
872 if (element.level() > outside.level())
873 {
874 interactionVolume.setSubVolumeElement(outside, 4);
875 interactionVolume.setSubVolumeElement(outside, 5);
876 interactionVolume.setSubVolumeElement(outside, 6);
877 interactionVolume.setSubVolumeElement(outside, 7);
878
879 break;
880 }
881 }
882 }
883
884 }
885 }
886 else if (interactionVolume.getHangingNodeType() == InteractionVolume::fourSmallCellsEdge)
887 {
888 auto element1 = interactionVolume.getSubVolumeElement(0);
889 auto element2 = interactionVolume.getSubVolumeElement(1);
890 auto element3 = interactionVolume.getSubVolumeElement(2);
891 auto element4 = interactionVolume.getSubVolumeElement(3);
892
893 for (const auto& intersection1 : intersections(problem_.gridView(), element1))
894 {
895 if (intersection1.neighbor())
896 {
897 if (intersection1.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
898 {
899 auto outside = intersection1.outside();
900 if (element1.level() > outside.level())
901 {
902 interactionVolume.setSubVolumeElement(outside, 4);
903 }
904 }
905 }
906 }
907 for (const auto& intersection2 : intersections(problem_.gridView(), element2))
908 {
909 if (intersection2.neighbor())
910 {
911 if (intersection2.indexInInside() == interactionVolume.getIndexOnElement(1, 2))
912 {
913 auto outside = intersection2.outside();
914 if (element2.level() > outside.level())
915 {
916 interactionVolume.setSubVolumeElement(outside, 5);
917
918 break;
919 }
920 }
921 }
922 }
923 for (const auto& intersection3 : intersections(problem_.gridView(), element3))
924 {
925 if (intersection3.neighbor())
926 {
927 if (intersection3.indexInInside() == interactionVolume.getIndexOnElement(2, 2))
928 {
929 auto outside = intersection3.outside();
930 if (element3.level() > outside.level())
931 {
932 interactionVolume.setSubVolumeElement(outside, 6);
933 break;
934 }
935 }
936 }
937 }
938 for (const auto& intersection4 : intersections(problem_.gridView(), element4))
939 {
940 if (intersection4.neighbor())
941 {
942 if (intersection4.indexInInside() == interactionVolume.getIndexOnElement(3, 2))
943 {
944 auto outside = intersection4.outside();
945 if (element4.level() > outside.level())
946 {
947 interactionVolume.setSubVolumeElement(outside, 7);
948
949 break;
950 }
951 }
952 }
953 }
954
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));
960
961 DimVector crossProductVector1(0);
962 DimVector crossProductVector2(0);
963
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);
969
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);
975
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);
981
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);
987
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);
993
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);
999
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);
1005
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);
1011
1012 auto element5 = interactionVolume.getSubVolumeElement(4);
1013 auto element6 = interactionVolume.getSubVolumeElement(5);
1014 auto element7 = interactionVolume.getSubVolumeElement(6);
1015 auto element8 = interactionVolume.getSubVolumeElement(7);
1016
1017 if (element5 == element6)
1018 {
1019 interactionVolume.setFacePosition(element5.geometry().center(), 4);
1020 interactionVolume.setFacePosition(element7.geometry().center(), 6);
1021
1022 for (const auto& intersection : intersections(problem_.gridView(), element5))
1023 {
1024 if (intersection.neighbor())
1025 {
1026 auto outside = intersection.outside();
1027
1028 if (outside == element7 || outside == element8)
1029 {
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);
1042 normal *= -1;
1043 interactionVolume.setNormal(normal, 6, 1);
1044 interactionVolume.setNormal(normal, 7, 2);
1045
1046 break;
1047 }
1048 }
1049 }
1050
1051 DimVector edgeCoord2(interactionVolume.getFacePosition(7));
1052
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);
1058 // interactionVolume.setFaceArea(0.0, 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);
1064 // interactionVolume.setFaceArea(0.0, 7);
1065 }
1066 else if (element5 == element7)
1067 {
1068 interactionVolume.setFacePosition(element6.geometry().center(), 5);
1069 interactionVolume.setFacePosition(element5.geometry().center(), 7);
1070 interactionVolume.setFacePosition(globalPosFace, 6);
1071
1072 for (const auto& intersection : intersections(problem_.gridView(), element5))
1073 {
1074 if (intersection.neighbor())
1075 {
1076 auto outside = intersection.outside();
1077
1078 if (outside == element6 || outside == element8)
1079 {
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);
1092 normal *= -1;
1093 interactionVolume.setNormal(normal, 5, 2);
1094 interactionVolume.setNormal(normal, 7, 1);
1095
1096 break;
1097 }
1098 }
1099 }
1100 DimVector edgeCoord2(interactionVolume.getFacePosition(4));
1101
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);
1107
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);
1113 // interactionVolume.setFaceArea(0.0, 4);
1114 // interactionVolume.setFaceArea(0.0, 6);
1115 }
1116
1117 const ElementGeometry& geometry = element5.geometry();
1118
1119 const auto refElement = referenceElement(geometry);
1120
1121 int oldSubVolumElemIdx = IndexTranslator::getOldElemIdxFromNewFaceIdxto0(zeroFaceIdx, 4);
1122 int oldEdgeIdx = IndexTranslator::getOldEdgeIdxFromNewFaceIdxto0(zeroFaceIdx, 1);
1123
1124 DimVector edgeCoord(0.0);
1125 switch (oldSubVolumElemIdx)
1126 {
1127 case 0:
1128 {
1129 switch (oldEdgeIdx)
1130 {
1131 case 2:
1132 edgeCoord = geometry.global(refElement.position(9, dim - 1));
1133 break;
1134 case 0:
1135 edgeCoord = geometry.global(refElement.position(3, dim - 1));
1136 break;
1137 case 5:
1138 edgeCoord = geometry.global(refElement.position(11, dim - 1));
1139 break;
1140 }
1141
1142 break;
1143 }
1144 case 1:
1145 {
1146 switch (oldEdgeIdx)
1147 {
1148 case 0:
1149 edgeCoord = geometry.global(refElement.position(2, dim - 1));
1150 break;
1151 case 2:
1152 edgeCoord = geometry.global(refElement.position(8, dim - 1));
1153 break;
1154 case 3:
1155 edgeCoord = geometry.global(refElement.position(11, dim - 1));
1156 break;
1157 }
1158
1159 break;
1160 }
1161 case 2:
1162 {
1163 switch (oldEdgeIdx)
1164 {
1165 case 0:
1166 edgeCoord = geometry.global(refElement.position(1, dim - 1));
1167 break;
1168 case 4:
1169 edgeCoord = geometry.global(refElement.position(9, dim - 1));
1170 break;
1171 case 5:
1172 edgeCoord = geometry.global(refElement.position(10, dim - 1));
1173 break;
1174 }
1175
1176 break;
1177 }
1178 case 3:
1179 {
1180 switch (oldEdgeIdx)
1181 {
1182 case 0:
1183 edgeCoord = geometry.global(refElement.position(0, dim - 1));
1184 break;
1185 case 4:
1186 edgeCoord = geometry.global(refElement.position(8, dim - 1));
1187 break;
1188 case 3:
1189 edgeCoord = geometry.global(refElement.position(10, dim - 1));
1190 break;
1191 }
1192
1193 break;
1194 }
1195 case 4:
1196 {
1197 switch (oldEdgeIdx)
1198 {
1199 case 1:
1200 edgeCoord = geometry.global(refElement.position(3, dim - 1));
1201 break;
1202 case 2:
1203 edgeCoord = geometry.global(refElement.position(5, dim - 1));
1204 break;
1205 case 5:
1206 edgeCoord = geometry.global(refElement.position(7, dim - 1));
1207 break;
1208 }
1209
1210 break;
1211 }
1212 case 5:
1213 {
1214 switch (oldEdgeIdx)
1215 {
1216 case 1:
1217 edgeCoord = geometry.global(refElement.position(2, dim - 1));
1218 break;
1219 case 2:
1220 edgeCoord = geometry.global(refElement.position(4, dim - 1));
1221 break;
1222 case 3:
1223 edgeCoord = geometry.global(refElement.position(7, dim - 1));
1224 break;
1225 }
1226
1227 break;
1228 }
1229 case 6:
1230 {
1231 switch (oldEdgeIdx)
1232 {
1233 case 1:
1234 edgeCoord = geometry.global(refElement.position(1, dim - 1));
1235 break;
1236 case 4:
1237 edgeCoord = geometry.global(refElement.position(5, dim - 1));
1238 break;
1239 case 5:
1240 edgeCoord = geometry.global(refElement.position(6, dim - 1));
1241 break;
1242 }
1243
1244 break;
1245 }
1246 case 7:
1247 {
1248 switch (oldEdgeIdx)
1249 {
1250 case 1:
1251 edgeCoord = geometry.global(refElement.position(0, dim - 1));
1252 break;
1253 case 4:
1254 edgeCoord = geometry.global(refElement.position(4, dim - 1));
1255 break;
1256 case 3:
1257 edgeCoord = geometry.global(refElement.position(6, dim - 1));
1258 break;
1259 }
1260
1261 break;
1262 }
1263 }
1264
1265 interactionVolume.setEdgePosition(edgeCoord, 1);
1266 }
1267 else if (interactionVolume.getHangingNodeType() == InteractionVolume::fourSmallCellsDiag)
1268 {
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));
1275
1276 DimVector crossProductVector1(0);
1277 DimVector crossProductVector2(0);
1278
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);
1284
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);
1290
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);
1296
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);
1302
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);
1308
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);
1314
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);
1320
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);
1326
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);
1332
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);
1338
1339 auto element1 = interactionVolume.getSubVolumeElement(0);
1340
1341 bool hasFaceOne = false;
1342 bool hasFaceTwo = false;
1343 for (const auto& intersection : intersections(problem_.gridView(), element1))
1344 {
1345 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(0, 1))
1346 {
1347 if (intersection.neighbor())
1348 {
1349 auto outside = intersection.outside();
1350 if (element1.level() > outside.level())
1351 {
1352 interactionVolume.setSubVolumeElement(outside, 2);
1353 interactionVolume.setSubVolumeElement(outside, 3);
1354
1355 hasFaceOne = true;
1356 if (hasFaceTwo)
1357 break;
1358 }
1359 }
1360 }
1361 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(0, 2))
1362 {
1363 if (intersection.neighbor())
1364 {
1365 auto outside = intersection.outside();
1366 if (element1.level() > outside.level())
1367 {
1368 interactionVolume.setSubVolumeElement(outside, 4);
1369 interactionVolume.setSubVolumeElement(outside, 5);
1370
1371 hasFaceTwo = true;
1372 if (hasFaceOne)
1373 break;
1374 }
1375 }
1376 }
1377 }
1378 }
1379
1380 break;
1381 }
1382 //hanging-node interaction volume of type 2 or 6
1383 case 6:
1384 {
1385 InteractionVolume hangingNodeVolume(problem_.gridView().grid());
1386
1387 std::vector<int> elemIdxOld;
1388 for (int i = 0; i < InteractionVolume::subVolumeTotalNum; i++)
1389 {
1390 if (interactionVolume.hasSubVolumeElement(i))
1391 elemIdxOld.push_back(i);
1392 }
1393
1394 std::set<int> zeroFaceIdxVec;
1395 for (int i = 0; i < 6; i++)
1396 {
1397 for (int j = 0; j < 6; j++)
1398 {
1399 int fIdx = IndexTranslator::getFaceIndexFromElements(elemIdxOld[i], elemIdxOld[j]);
1400 if (fIdx >= 0)
1401 zeroFaceIdxVec.insert(fIdx);
1402 }
1403 }
1404
1405 std::vector<int> elemIdxNew(6);
1406 int zeroFaceIdx = 0;
1407
1408 bool isFace = true;
1409 std::set<int>::iterator it = zeroFaceIdxVec.begin();
1410 for (; it != zeroFaceIdxVec.end(); ++it)
1411 {
1412 zeroFaceIdx = *it;
1413 isFace = true;
1414 // bool isFace = true;
1415 for (int i = 0; i < 6; i++)
1416 {
1417 elemIdxNew[i] = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elemIdxOld[i]);
1418 if (elemIdxNew[i] == 6 || elemIdxNew[i] == 7)
1419 {
1420 isFace = false;
1421 break;
1422 }
1423 }
1424 if (isFace)
1425 {
1426 for (int i = 0; i < 6; i++)
1427 {
1428 hangingNodeVolume.setSubVolumeElement(interactionVolume.getSubVolumeElement(elemIdxOld[i]),
1429 elemIdxNew[i]);
1430 }
1431
1432 break;
1433 }
1434 }
1435
1436 for (int elem = 0; elem < InteractionVolume::subVolumeTotalNum; elem++)
1437 {
1438 for (int i = 0; i < 3; i++)
1439 {
1440 int faceIdxOld = IndexTranslator::getFaceIndexFromSubVolume(elem, i);
1441 int faceIdxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, faceIdxOld);
1442
1443 for (int j = 0; j < 3; j++)
1444 {
1445 int elemNew = IndexTranslator::getNewElemIdxFromOldFaceIdxto0(zeroFaceIdx, elem);
1446 int faceIdxNewTest = IndexTranslator::getFaceIndexFromSubVolume(elemNew, j);
1447
1448 if (faceIdxNew == faceIdxNewTest)
1449 {
1450 hangingNodeVolume.setNormal(interactionVolume.getNormal(elem, i),
1451 elemNew, j);
1452 hangingNodeVolume.setIndexOnElement(
1453 interactionVolume.getIndexOnElement(elem, i), elemNew, j);
1454 }
1455 }
1456 }
1457 }
1458
1459 for (int i = 0; i < InteractionVolume::fluxFacesTotalNum; i++)
1460 {
1461 int idxNew = IndexTranslator::getNewFaceIdxFromOldIdxto0(zeroFaceIdx, i);
1462 hangingNodeVolume.setFaceArea(interactionVolume.getFaceArea(i), idxNew);
1463 hangingNodeVolume.setFacePosition(interactionVolume.getFacePosition(i), idxNew);
1464 }
1465
1466 for (int i = 0; i < InteractionVolume::fluxEdgesTotalNum; i++)
1467 {
1468 int idxNew = IndexTranslator::getNewEdgeIdxFromOldFaceIdxto0(zeroFaceIdx, i);
1469 hangingNodeVolume.setEdgePosition(interactionVolume.getEdgePosition(i), idxNew);
1470 }
1471
1472 interactionVolume = hangingNodeVolume;
1473
1474 interactionVolume.setCenterPosition(centerPos);
1475
1476 interactionVolume.setHangingNodeType(InteractionVolume::sixSmallCells);
1477
1478 auto element3 = interactionVolume.getSubVolumeElement(2);
1479
1480 for (const auto& intersection : intersections(problem_.gridView(), element3))
1481 {
1482 if (intersection.indexInInside() == interactionVolume.getIndexOnElement(2, 2))
1483 {
1484 if (intersection.neighbor())
1485 {
1486 auto outside = intersection.outside();
1487 if (element3.level() > outside.level())
1488 {
1489 interactionVolume.setSubVolumeElement(outside, 6);
1490 interactionVolume.setSubVolumeElement(outside, 7);
1491
1492 break;
1493 }
1494 }
1495 }
1496 }
1497
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));
1504
1505 DimVector crossProductVector1(0);
1506 DimVector crossProductVector2(0);
1507
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);
1513
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);
1519
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);
1525
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);
1531
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);
1537
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);
1543
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);
1549
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);
1555
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);
1561
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);
1567
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);
1573 break;
1574 }
1575 default:
1576 {
1577 interactionVolume.printInteractionVolumeInfo();
1578 DUNE_THROW(Dune::NotImplemented, "Hanging node shape not implemented");
1579 break;
1580 }
1581 }
1582}
1583
1585template<class TypeTag>
1587{
1588 std::vector < std::vector<int> > elemVertMap(problem_.gridView().size(dim), std::vector<int>(8, -1));
1589
1590 //Add elements to the interaction volumes and store element-vertex map
1591 for (const auto& element : elements(problem_.gridView()))
1592 asImp_().storeSubVolumeElements(element, elemVertMap);
1593
1594 for (unsigned int i = 0; i < asImp_().interactionVolumes_.size(); i++)
1595 if (asImp_().interactionVolumes_[i].getElementNumber() == 0)
1596 asImp_().interactionVolumes_[i].printInteractionVolumeInfo();
1597
1598 // Store information related to DUNE intersections for all interaction volumes
1599 for (const auto& element : elements(problem_.gridView()))
1600 asImp_().storeIntersectionInfo(element, elemVertMap);
1601
1602 faceVertices_.clear();
1603 faceVertices_.resize(problem_.gridView().size(0));
1604
1605 // Complete storage of the interaction volumes using the previously stored information
1606 // about the orientation and relationship of the DUNE elements in the interaction volumes (see doc/docextra/3dmpfa)
1607 for (const auto& vertex : vertices(problem_.gridView()))
1608 {
1609 int vIdxGlobal = problem_.variables().index(vertex);
1610
1611 InteractionVolume& interactionVolume = asImp_().interactionVolumes_[vIdxGlobal];
1612
1613 if (interactionVolume.getElementNumber() == 8)
1614 {
1615 asImp_().storeInnerInteractionVolume(interactionVolume, vertex);
1616 }
1617 else if (interactionVolume.isBoundaryInteractionVolume())
1618 {
1619 asImp_().storeBoundaryInteractionVolume(interactionVolume, vertex);
1620 }
1621 //hanging node!
1622 else
1623 {
1624 storeHangingNodeInteractionVolume(interactionVolume, vertex);
1625 }
1626
1627 if (!interactionVolume.isBoundaryInteractionVolume())
1628 {
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);
1637
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);
1646
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));
1657
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);
1668
1669 if (interactionVolume.getHangingNodeType() != InteractionVolume::twoSmallCells)
1670 {
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));
1685
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);
1700 }
1701 }
1702
1703 }
1704}
1705} // end namespace Dumux
1706#endif
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: 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:654
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property
Definition: propertysystem.hh:150
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.