3.2-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
3dinteractionvolumecontainer.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_HH
25#define DUMUX_FVMPFAL3D_INTERACTIONVOLUMECONTAINER_HH
26
27// dumux environment
28#include <dumux/porousmediumflow/sequential/pressureproperties.hh>
29#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
30#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/linteractionvolume3d.hh>
31
32namespace Dumux {
33
44template<class TypeTag>
45class FvMpfaL3dInteractionVolumeContainer
46{
47 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
48 using Implementation = GetPropType<TypeTag, Properties::MPFAInteractionVolumeContainer>;
49
50 enum
51 {
52 dim = GridView::dimension, dimWorld = GridView::dimensionworld
53 };
54
55 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
56 using Problem = GetPropType<TypeTag, Properties::Problem>;
57
58 using ReferenceElements = Dune::ReferenceElements<Scalar, dim>;
59
60 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
61
62 using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
63 using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
64 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
65
66 using Element = typename GridView::Traits::template Codim<0>::Entity;
67 using Vertex = typename GridView::Traits::template Codim<dim>::Entity;
68 using ElementGeometry = typename Element::Geometry;
69
70 using Intersection = typename GridView::Intersection;
71 using IntersectionGeometry = typename Intersection::Geometry;
72
73 using GlobalPosition = typename ElementGeometry::GlobalCoordinate;
74 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
75
76 using DimVector = Dune::FieldVector<Scalar, dim>;
77
78 enum
79 {
80 pressureEqIdx = Indices::pressureEqIdx,
81 };
82
83 enum
84 {
85 innerEdgeFace = 2, innerSideFace = 1
86 };
87 enum
88 {
89 realFaceArea = 0, fluxFaceArea = 1
90 };
91
92public:
94 using InteractionVolume = GetPropType<TypeTag, Properties::MPFAInteractionVolume>;
95
96private:
97 using GlobalInteractionVolumeVector = std::vector<InteractionVolume>;
98 using FaceAreaVector = std::vector<Dune::FieldVector<Dune::FieldVector<Scalar, 2>, 2*dim> >;
99protected:
100 void storeSubVolumeElements(const Element& element, std::vector < std::vector<int> >& elemVertMap);
101 void storeIntersectionInfo(const Element& element, std::vector < std::vector<int> >& elemVertMap);
102 void storeInnerInteractionVolume(InteractionVolume& interactionVolume, const Vertex& vertex, bool sameLevel = true);
103 void storeBoundaryInteractionVolume(InteractionVolume& interactionVolume, const Vertex& vertex);
104private:
105 void storeInteractionVolumeInfo();
106public:
107
113 void update()
114 {
115 interactionVolumes_.clear();
116 realFluxFaceArea_.clear();
117
118 realFluxFaceArea_.resize(problem_.gridView().size(dim),
119 Dune::FieldVector<Dune::FieldVector<Scalar, 2>, 2 * dim>(Dune::FieldVector<Scalar, 2>(0.0)));
120 interactionVolumes_.resize(problem_.gridView().size(dim), InteractionVolume(problem_.gridView().grid()));
121
122 asImp_().storeInteractionVolumeInfo();
123 }
124
125
131 void initialize(bool solveTwice = true)
132 {
133 update();
134
135 return;
136 }
137
143 InteractionVolume& interactionVolume(int vertexIdx)
144 {
145 return interactionVolumes_[vertexIdx];
146 }
147
153 InteractionVolume& interactionVolume(int vertexIdx) const
154 {
155 return interactionVolumes_[vertexIdx];
156 }
157
159 GlobalInteractionVolumeVector& interactionVolumesGlobal()
160 {
161 return interactionVolumes_;
162 }
163
165 GlobalInteractionVolumeVector& interactionVolumesGlobal() const
166 {
167 return interactionVolumes_;
168 }
169
181 Scalar faceAreaFactor(InteractionVolume& interactionVolume, int elemGlobalIdx, int elemLocalIdx, int localFaceIdx)
182 {
183 Scalar factor = getRealFaceArea(interactionVolume, elemGlobalIdx, elemLocalIdx, localFaceIdx);
184 factor /= getRealFluxFaceArea(interactionVolume, elemGlobalIdx, elemLocalIdx, localFaceIdx);
185
186 return factor;
187 }
188
198 Scalar faceAreaFactor(int elemGlobalIdx, int indexInInside)
199 {
200 Scalar factor = getRealFaceArea(elemGlobalIdx, indexInInside);
201 factor /= getRealFluxFaceArea(elemGlobalIdx, indexInInside);
202
203 return factor;
204 }
205
216 Scalar getRealFluxFaceArea(InteractionVolume& interactionVolume, int elemGlobalIdx, int elemLocalIdx, int localFaceIdx)
217 {
218 Scalar factor = realFluxFaceArea_[elemGlobalIdx][interactionVolume.getIndexOnElement(elemLocalIdx, localFaceIdx)][fluxFaceArea];
219
220 return factor;
221 }
222
231 Scalar getRealFluxFaceArea(int elemGlobalIdx, int indexInInside)
232 {
233 Scalar factor = realFluxFaceArea_[elemGlobalIdx][indexInInside][fluxFaceArea];
234
235 return factor;
236 }
237
248 Scalar getRealFaceArea(InteractionVolume& interactionVolume, int elemGlobalIdx, int elemLocalIdx, int localFaceIdx)
249 {
250 Scalar factor = realFluxFaceArea_[elemGlobalIdx][interactionVolume.getIndexOnElement(elemLocalIdx, localFaceIdx)][realFaceArea];
251
252 return factor;
253 }
254
263 Scalar getRealFaceArea(int elemGlobalIdx, int indexInInside)
264 {
265 Scalar factor = realFluxFaceArea_[elemGlobalIdx][indexInInside][realFaceArea];
266
267 return factor;
268 }
269
275 FvMpfaL3dInteractionVolumeContainer(Problem& problem) :
276 problem_(problem)
277 {
278 if (dim != 3)
279 {
280 DUNE_THROW(Dune::NotImplemented, "Dimension not supported!");
281 }
282 }
283
284protected:
285 void addRealFluxFaceArea_(Scalar faceArea, int eIdxGlobal, int fIdx)
286 {
287 realFluxFaceArea_[eIdxGlobal][fIdx][fluxFaceArea] += faceArea;
288 }
289 void addRealFaceArea_(Scalar faceArea, int eIdxGlobal, int fIdx)
290 {
291 realFluxFaceArea_[eIdxGlobal][fIdx][realFaceArea] += faceArea;
292 }
293
294 Problem& problem_;
295
296 GlobalInteractionVolumeVector interactionVolumes_;
297 FaceAreaVector realFluxFaceArea_;
298private:
300 Implementation &asImp_()
301 { return *static_cast<Implementation *>(this); }
302
304 const Implementation &asImp_() const
305 { return *static_cast<const Implementation *>(this); }
306};
307
319template<class TypeTag>
320void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeSubVolumeElements(const Element& element,
321 std::vector < std::vector<int> >& elemVertMap)
322{
323 int eIdxGlobal = problem_.variables().index(element);
324 int vIdxGlobal = problem_.variables().vertexMapper().subIndex(element, 0, dim);
325 interactionVolumes_[vIdxGlobal].setSubVolumeElement(element, 7);
326 elemVertMap[vIdxGlobal][7] = eIdxGlobal;
327
328 vIdxGlobal = problem_.variables().vertexMapper().subIndex(element, 1, dim);
329 interactionVolumes_[vIdxGlobal].setSubVolumeElement(element, 6);
330 elemVertMap[vIdxGlobal][6] = eIdxGlobal;
331
332 vIdxGlobal = problem_.variables().vertexMapper().subIndex(element, 2, dim);
333 interactionVolumes_[vIdxGlobal].setSubVolumeElement(element, 5);
334 elemVertMap[vIdxGlobal][5] = eIdxGlobal;
335
336 vIdxGlobal = problem_.variables().vertexMapper().subIndex(element, 3, dim);
337 interactionVolumes_[vIdxGlobal].setSubVolumeElement(element, 4);
338 elemVertMap[vIdxGlobal][4] = eIdxGlobal;
339
340 vIdxGlobal = problem_.variables().vertexMapper().subIndex(element, 4, dim);
341 interactionVolumes_[vIdxGlobal].setSubVolumeElement(element, 3);
342 elemVertMap[vIdxGlobal][3] = eIdxGlobal;
343
344 vIdxGlobal = problem_.variables().vertexMapper().subIndex(element, 5, dim);
345 interactionVolumes_[vIdxGlobal].setSubVolumeElement(element, 2);
346 elemVertMap[vIdxGlobal][2] = eIdxGlobal;
347
348 vIdxGlobal = problem_.variables().vertexMapper().subIndex(element, 6, dim);
349 interactionVolumes_[vIdxGlobal].setSubVolumeElement(element, 1);
350 elemVertMap[vIdxGlobal][1] = eIdxGlobal;
351
352 vIdxGlobal = problem_.variables().vertexMapper().subIndex(element, 7, dim);
353 interactionVolumes_[vIdxGlobal].setSubVolumeElement(element, 0);
354 elemVertMap[vIdxGlobal][0] = eIdxGlobal;
355
356}
357
370template<class TypeTag>
371void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeIntersectionInfo(const Element& element,
372 std::vector < std::vector<int> >& elemVertMap)
373{
374 BoundaryTypes bcType;
375
376 int eIdxGlobal = problem_.variables().index(element);
377
378 const ElementGeometry& geometry = element.geometry();
379
380 const auto referenceElement = ReferenceElements::general(geometry.type());
381
382 int levelI = element.level();
383
384 // run through all intersections
385 for (const auto& intersection : intersections(problem_.gridView(), element))
386 {
387 int indexInInside = intersection.indexInInside();
388
389 DimVector normal = intersection.centerUnitOuterNormal();
390
391 const IntersectionGeometry& isGeometry = intersection.geometry();
392
393 Scalar faceVol = isGeometry.volume();
394
395 const DimVector& globalPosFace = isGeometry.center();
396
397 bool takeIntersection = true;
398 if (intersection.neighbor())
399 {
400 auto outside = intersection.outside();
401 int eIdxGlobalJ = problem_.variables().index(outside);
402
403 if (levelI == outside.level() && eIdxGlobal > eIdxGlobalJ)
404 takeIntersection = false;
405 if (levelI < outside.level())
406 takeIntersection = false;
407 }
408
409 if (takeIntersection)
410 {
411 addRealFaceArea_(faceVol, eIdxGlobal, indexInInside);
412 if (intersection.neighbor())
413 {
414 int eIdxGlobalJ = problem_.variables().index(intersection.outside());
415 addRealFaceArea_(faceVol, eIdxGlobalJ, intersection.indexInOutside());
416 }
417
418 for (int i = 0; i < isGeometry.corners(); i++)
419 {
420 int localVertIdx = referenceElement.subEntity(indexInInside, 1, i, dim);
421
422 int vIdxGlobal = problem_.variables().vertexMapper().subIndex(element, localVertIdx, dim);
423
424 InteractionVolume& interactionVolume = interactionVolumes_[vIdxGlobal];
425
426 if (elemVertMap[vIdxGlobal][0] == eIdxGlobal)
427 {
428 if (indexInInside == 1)
429 {
430 interactionVolume.setIndexOnElement(indexInInside, 0, 0);
431 interactionVolume.setNormal(normal, 0, 0);
432 interactionVolume.setFacePosition(globalPosFace, 0);
433
434 if (intersection.neighbor())
435 {
436 int indexInOutside = intersection.indexInOutside();
437
438 interactionVolume.setIndexOnElement(indexInOutside, 1, 1);
439 DimVector normalOutside = normal;
440 normalOutside *= -1;
441
442 interactionVolume.setNormal(normalOutside, 1, 1);
443 }
444 }
445 else if (indexInInside == 3)
446 {
447 interactionVolume.setIndexOnElement(indexInInside, 0, 1);
448 interactionVolume.setNormal(normal, 0, 1);
449 interactionVolume.setFacePosition(globalPosFace, 3);
450
451 if (intersection.neighbor())
452 {
453 int indexInOutside = intersection.indexInOutside();
454
455 interactionVolume.setIndexOnElement(indexInOutside, 2, 0);
456 DimVector normalOutside = normal;
457 normalOutside *= -1;
458
459 interactionVolume.setNormal(normalOutside, 2, 0);
460 }
461 }
462 else if (indexInInside == 5)
463 {
464 interactionVolume.setIndexOnElement(indexInInside, 0, 2);
465 interactionVolume.setNormal(normal, 0, 2);
466 interactionVolume.setFacePosition(globalPosFace, 8);
467
468 if (intersection.neighbor())
469 {
470 int indexInOutside = intersection.indexInOutside();
471
472 interactionVolume.setIndexOnElement(indexInOutside, 4, 0);
473 DimVector normalOutside = normal;
474 normalOutside *= -1;
475
476 interactionVolume.setNormal(normalOutside, 4, 0);
477 }
478 }
479 }
480 if (elemVertMap[vIdxGlobal][1] == eIdxGlobal)
481 {
482 if (indexInInside == 3)
483 {
484 interactionVolume.setIndexOnElement(indexInInside, 1, 0);
485 interactionVolume.setNormal(normal, 1, 0);
486 interactionVolume.setFacePosition(globalPosFace, 1);
487
488 if (intersection.neighbor())
489 {
490 int indexInOutside = intersection.indexInOutside();
491
492 interactionVolume.setIndexOnElement(indexInOutside, 3, 1);
493 DimVector normalOutside = normal;
494 normalOutside *= -1;
495
496 interactionVolume.setNormal(normalOutside, 3, 1);
497 }
498 }
499 else if (indexInInside == 0)
500 {
501 interactionVolume.setIndexOnElement(indexInInside, 1, 1);
502 interactionVolume.setNormal(normal, 1, 1);
503 interactionVolume.setFacePosition(globalPosFace, 0);
504
505 if (intersection.neighbor())
506 {
507 int indexInOutside = intersection.indexInOutside();
508
509 interactionVolume.setIndexOnElement(indexInOutside, 0, 0);
510 DimVector normalOutside = normal;
511 normalOutside *= -1;
512
513 interactionVolume.setNormal(normalOutside, 0, 0);
514 }
515 }
516 else if (indexInInside == 5)
517 {
518 interactionVolume.setIndexOnElement(indexInInside, 1, 2);
519 interactionVolume.setNormal(normal, 1, 2);
520 interactionVolume.setFacePosition(globalPosFace, 9);
521
522 if (intersection.neighbor())
523 {
524 int indexInOutside = intersection.indexInOutside();
525
526 interactionVolume.setIndexOnElement(indexInOutside, 5, 0);
527 DimVector normalOutside = normal;
528 normalOutside *= -1;
529
530 interactionVolume.setNormal(normalOutside, 5, 0);
531 }
532 }
533 }
534 if (elemVertMap[vIdxGlobal][2] == eIdxGlobal)
535 {
536 if (indexInInside == 2)
537 {
538 interactionVolume.setIndexOnElement(indexInInside, 2, 0);
539 interactionVolume.setNormal(normal, 2, 0);
540 interactionVolume.setFacePosition(globalPosFace, 3);
541
542 if (intersection.neighbor())
543 {
544 int indexInOutside = intersection.indexInOutside();
545
546 interactionVolume.setIndexOnElement(indexInOutside, 0, 1);
547 DimVector normalOutside = normal;
548 normalOutside *= -1;
549
550 interactionVolume.setNormal(normalOutside, 0, 1);
551 }
552 }
553 else if (indexInInside == 1)
554 {
555 interactionVolume.setIndexOnElement(indexInInside, 2, 1);
556 interactionVolume.setNormal(normal, 2, 1);
557 interactionVolume.setFacePosition(globalPosFace, 2);
558
559 if (intersection.neighbor())
560 {
561 int indexInOutside = intersection.indexInOutside();
562
563 interactionVolume.setIndexOnElement(indexInOutside, 3, 0);
564 DimVector normalOutside = normal;
565 normalOutside *= -1;
566
567 interactionVolume.setNormal(normalOutside, 3, 0);
568 }
569 }
570 else if (indexInInside == 5)
571 {
572 interactionVolume.setIndexOnElement(indexInInside, 2, 2);
573 interactionVolume.setNormal(normal, 2, 2);
574 interactionVolume.setFacePosition(globalPosFace, 11);
575
576 if (intersection.neighbor())
577 {
578 int indexInOutside = intersection.indexInOutside();
579
580 interactionVolume.setIndexOnElement(indexInOutside, 6, 0);
581 DimVector normalOutside = normal;
582 normalOutside *= -1;
583
584 interactionVolume.setNormal(normalOutside, 6, 0);
585 }
586 }
587 }
588 if (elemVertMap[vIdxGlobal][3] == eIdxGlobal)
589 {
590 if (indexInInside == 0)
591 {
592 interactionVolume.setIndexOnElement(indexInInside, 3, 0);
593 interactionVolume.setNormal(normal, 3, 0);
594 interactionVolume.setFacePosition(globalPosFace, 2);
595
596 if (intersection.neighbor())
597 {
598 int indexInOutside = intersection.indexInOutside();
599
600 interactionVolume.setIndexOnElement(indexInOutside, 2, 1);
601 DimVector normalOutside = normal;
602 normalOutside *= -1;
603
604 interactionVolume.setNormal(normalOutside, 2, 1);
605 }
606 }
607 else if (indexInInside == 2)
608 {
609 interactionVolume.setIndexOnElement(indexInInside, 3, 1);
610 interactionVolume.setNormal(normal, 3, 1);
611 interactionVolume.setFacePosition(globalPosFace, 1);
612
613 if (intersection.neighbor())
614 {
615 int indexInOutside = intersection.indexInOutside();
616
617 interactionVolume.setIndexOnElement(indexInOutside, 1, 0);
618 DimVector normalOutside = normal;
619 normalOutside *= -1;
620
621 interactionVolume.setNormal(normalOutside, 1, 0);
622 }
623 }
624 else if (indexInInside == 5)
625 {
626 interactionVolume.setIndexOnElement(indexInInside, 3, 2);
627 interactionVolume.setNormal(normal, 3, 2);
628 interactionVolume.setFacePosition(globalPosFace, 10);
629
630 if (intersection.neighbor())
631 {
632 int indexInOutside = intersection.indexInOutside();
633
634 interactionVolume.setIndexOnElement(indexInOutside, 7, 0);
635 DimVector normalOutside = normal;
636 normalOutside *= -1;
637
638 interactionVolume.setNormal(normalOutside, 7, 0);
639 }
640 }
641 }
642 if (elemVertMap[vIdxGlobal][4] == eIdxGlobal)
643 {
644 if (indexInInside == 4)
645 {
646 interactionVolume.setIndexOnElement(indexInInside, 4, 0);
647 interactionVolume.setNormal(normal, 4, 0);
648 interactionVolume.setFacePosition(globalPosFace, 8);
649
650 if (intersection.neighbor())
651 {
652 int indexInOutside = intersection.indexInOutside();
653
654 interactionVolume.setIndexOnElement(indexInOutside, 0, 2);
655 DimVector normalOutside = normal;
656 normalOutside *= -1;
657
658 interactionVolume.setNormal(normalOutside, 0, 2);
659 }
660 }
661 else if (indexInInside == 1)
662 {
663 interactionVolume.setIndexOnElement(indexInInside, 4, 1);
664 interactionVolume.setNormal(normal, 4, 1);
665 interactionVolume.setFacePosition(globalPosFace, 4);
666
667 if (intersection.neighbor())
668 {
669 int indexInOutside = intersection.indexInOutside();
670
671 interactionVolume.setIndexOnElement(indexInOutside, 5, 2);
672 DimVector normalOutside = normal;
673 normalOutside *= -1;
674
675 interactionVolume.setNormal(normalOutside, 5, 2);
676 }
677 }
678 else if (indexInInside == 3)
679 {
680 interactionVolume.setIndexOnElement(indexInInside, 4, 2);
681 interactionVolume.setNormal(normal, 4, 2);
682 interactionVolume.setFacePosition(globalPosFace, 7);
683
684 if (intersection.neighbor())
685 {
686 int indexInOutside = intersection.indexInOutside();
687
688 interactionVolume.setIndexOnElement(indexInOutside, 6, 1);
689 DimVector normalOutside = normal;
690 normalOutside *= -1;
691
692 interactionVolume.setNormal(normalOutside, 6, 1);
693 }
694 }
695 }
696 if (elemVertMap[vIdxGlobal][5] == eIdxGlobal)
697 {
698 if (indexInInside == 4)
699 {
700 interactionVolume.setIndexOnElement(indexInInside, 5, 0);
701 interactionVolume.setNormal(normal, 5, 0);
702 interactionVolume.setFacePosition(globalPosFace, 9);
703
704 if (intersection.neighbor())
705 {
706 int indexInOutside = intersection.indexInOutside();
707
708 interactionVolume.setIndexOnElement(indexInOutside, 1, 2);
709 DimVector normalOutside = normal;
710 normalOutside *= -1;
711
712 interactionVolume.setNormal(normalOutside, 1, 2);
713 }
714 }
715 else if (indexInInside == 3)
716 {
717 interactionVolume.setIndexOnElement(indexInInside, 5, 1);
718 interactionVolume.setNormal(normal, 5, 1);
719 interactionVolume.setFacePosition(globalPosFace, 5);
720
721 if (intersection.neighbor())
722 {
723 int indexInOutside = intersection.indexInOutside();
724
725 interactionVolume.setIndexOnElement(indexInOutside, 7, 2);
726 DimVector normalOutside = normal;
727 normalOutside *= -1;
728
729 interactionVolume.setNormal(normalOutside, 7, 2);
730 }
731
732 }
733 else if (indexInInside == 0)
734 {
735 interactionVolume.setIndexOnElement(indexInInside, 5, 2);
736 interactionVolume.setNormal(normal, 5, 2);
737 interactionVolume.setFacePosition(globalPosFace, 4);
738
739 if (intersection.neighbor())
740 {
741 int indexInOutside = intersection.indexInOutside();
742
743 interactionVolume.setIndexOnElement(indexInOutside, 4, 1);
744 DimVector normalOutside = normal;
745 normalOutside *= -1;
746
747 interactionVolume.setNormal(normalOutside, 4, 1);
748 }
749 }
750 }
751 if (elemVertMap[vIdxGlobal][6] == eIdxGlobal)
752 {
753 if (indexInInside == 4)
754 {
755 interactionVolume.setIndexOnElement(indexInInside, 6, 0);
756 interactionVolume.setNormal(normal, 6, 0);
757 interactionVolume.setFacePosition(globalPosFace, 11);
758
759 if (intersection.neighbor())
760 {
761 int indexInOutside = intersection.indexInOutside();
762
763 interactionVolume.setIndexOnElement(indexInOutside, 2, 2);
764 DimVector normalOutside = normal;
765 normalOutside *= -1;
766
767 interactionVolume.setNormal(normalOutside, 2, 2);
768 }
769 }
770 else if (indexInInside == 2)
771 {
772 interactionVolume.setIndexOnElement(indexInInside, 6, 1);
773 interactionVolume.setNormal(normal, 6, 1);
774 interactionVolume.setFacePosition(globalPosFace, 7);
775
776 if (intersection.neighbor())
777 {
778 int indexInOutside = intersection.indexInOutside();
779
780 interactionVolume.setIndexOnElement(indexInOutside, 4, 2);
781 DimVector normalOutside = normal;
782 normalOutside *= -1;
783
784 interactionVolume.setNormal(normalOutside, 4, 2);
785 }
786 }
787 else if (indexInInside == 1)
788 {
789 interactionVolume.setIndexOnElement(indexInInside, 6, 2);
790 interactionVolume.setNormal(normal, 6, 2);
791 interactionVolume.setFacePosition(globalPosFace, 6);
792
793 if (intersection.neighbor())
794 {
795 int indexInOutside = intersection.indexInOutside();
796
797 interactionVolume.setIndexOnElement(indexInOutside, 7, 1);
798 DimVector normalOutside = normal;
799 normalOutside *= -1;
800
801 interactionVolume.setNormal(normalOutside, 7, 1);
802 }
803 }
804 }
805 if (elemVertMap[vIdxGlobal][7] == eIdxGlobal)
806 {
807 if (indexInInside == 4)
808 {
809 interactionVolume.setIndexOnElement(indexInInside, 7, 0);
810 interactionVolume.setNormal(normal, 7, 0);
811 interactionVolume.setFacePosition(globalPosFace, 10);
812
813 if (intersection.neighbor())
814 {
815 int indexInOutside = intersection.indexInOutside();
816
817 interactionVolume.setIndexOnElement(indexInOutside, 3, 2);
818 DimVector normalOutside = normal;
819 normalOutside *= -1;
820
821 interactionVolume.setNormal(normalOutside, 3, 2);
822 }
823 }
824 else if (indexInInside == 0)
825 {
826 interactionVolume.setIndexOnElement(indexInInside, 7, 1);
827 interactionVolume.setNormal(normal, 7, 1);
828 interactionVolume.setFacePosition(globalPosFace, 6);
829
830 if (intersection.neighbor())
831 {
832 int indexInOutside = intersection.indexInOutside();
833
834 interactionVolume.setIndexOnElement(indexInOutside, 6, 2);
835 DimVector normalOutside = normal;
836 normalOutside *= -1;
837
838 interactionVolume.setNormal(normalOutside, 6, 2);
839 }
840 }
841 else if (indexInInside == 2)
842 {
843 interactionVolume.setIndexOnElement(indexInInside, 7, 2);
844 interactionVolume.setNormal(normal, 7, 2);
845 interactionVolume.setFacePosition(globalPosFace, 5);
846
847 if (intersection.neighbor())
848 {
849 int indexInOutside = intersection.indexInOutside();
850
851 interactionVolume.setIndexOnElement(indexInOutside, 5, 1);
852 DimVector normalOutside = normal;
853 normalOutside *= -1;
854
855 interactionVolume.setNormal(normalOutside, 5, 1);
856 }
857 }
858 }
859 if (intersection.boundary())
860 {
861 if (elemVertMap[vIdxGlobal][0] == eIdxGlobal)
862 {
863 if (indexInInside == 1)
864 {
865 problem_.boundaryTypes(bcType, intersection);
866 PrimaryVariables boundValues(0.0);
867
868 interactionVolume.setBoundary(bcType, 0);
869 if (bcType.isNeumann(pressureEqIdx))
870 {
871 problem_.neumann(boundValues, intersection);
872 boundValues *= faceVol/4.0;
873 interactionVolume.setNeumannCondition(boundValues, 0);
874 }
875 if (bcType.hasDirichlet())
876 {
877 problem_.dirichlet(boundValues, intersection);
878 interactionVolume.setDirichletCondition(boundValues, 0);
879 }
880 }
881 else if (indexInInside == 3)
882 {
883 problem_.boundaryTypes(bcType, intersection);
884 PrimaryVariables boundValues(0.0);
885
886 interactionVolume.setBoundary(bcType, 3);
887 if (bcType.isNeumann(pressureEqIdx))
888 {
889 problem_.neumann(boundValues, intersection);
890 boundValues *= faceVol/4.0;
891 interactionVolume.setNeumannCondition(boundValues, 3);
892 }
893 if (bcType.hasDirichlet())
894 {
895 problem_.dirichlet(boundValues, intersection);
896 interactionVolume.setDirichletCondition(boundValues, 3);
897 }
898 }
899 else if (indexInInside == 5)
900 {
901 problem_.boundaryTypes(bcType, intersection);
902 PrimaryVariables boundValues(0.0);
903
904 interactionVolume.setBoundary(bcType, 8);
905 if (bcType.isNeumann(pressureEqIdx))
906 {
907 problem_.neumann(boundValues, intersection);
908 boundValues *= faceVol/4.0;
909 interactionVolume.setNeumannCondition(boundValues, 8);
910 }
911 if (bcType.hasDirichlet())
912 {
913 problem_.dirichlet(boundValues, intersection);
914 interactionVolume.setDirichletCondition(boundValues, 8);
915 }
916 }
917 }
918 if (elemVertMap[vIdxGlobal][1] == eIdxGlobal)
919 {
920 if (indexInInside == 3)
921 {
922 problem_.boundaryTypes(bcType, intersection);
923 PrimaryVariables boundValues(0.0);
924
925 interactionVolume.setBoundary(bcType, 1);
926 if (bcType.isNeumann(pressureEqIdx))
927 {
928 problem_.neumann(boundValues, intersection);
929 boundValues *= faceVol/4.0;
930 interactionVolume.setNeumannCondition(boundValues, 1);
931 }
932 if (bcType.hasDirichlet())
933 {
934 problem_.dirichlet(boundValues, intersection);
935 interactionVolume.setDirichletCondition(boundValues, 1);
936 }
937 }
938 else if (indexInInside == 0)
939 {
940 problem_.boundaryTypes(bcType, intersection);
941 PrimaryVariables boundValues(0.0);
942
943 interactionVolume.setBoundary(bcType, 0);
944 if (bcType.isNeumann(pressureEqIdx))
945 {
946 problem_.neumann(boundValues, intersection);
947 boundValues *= faceVol/4.0;
948 interactionVolume.setNeumannCondition(boundValues, 0);
949 }
950 if (bcType.hasDirichlet())
951 {
952 problem_.dirichlet(boundValues, intersection);
953 interactionVolume.setDirichletCondition(boundValues, 0);
954 }
955 }
956 else if (indexInInside == 5)
957 {
958 problem_.boundaryTypes(bcType, intersection);
959 PrimaryVariables boundValues(0.0);
960
961 interactionVolume.setBoundary(bcType, 9);
962 if (bcType.isNeumann(pressureEqIdx))
963 {
964 problem_.neumann(boundValues, intersection);
965 boundValues *= faceVol/4.0;
966 interactionVolume.setNeumannCondition(boundValues, 9);
967 }
968 if (bcType.hasDirichlet())
969 {
970 problem_.dirichlet(boundValues, intersection);
971 interactionVolume.setDirichletCondition(boundValues, 9);
972 }
973 }
974 }
975 if (elemVertMap[vIdxGlobal][2] == eIdxGlobal)
976 {
977 if (indexInInside == 2)
978 {
979 problem_.boundaryTypes(bcType, intersection);
980 PrimaryVariables boundValues(0.0);
981
982 interactionVolume.setBoundary(bcType, 3);
983 if (bcType.isNeumann(pressureEqIdx))
984 {
985 problem_.neumann(boundValues, intersection);
986 boundValues *= faceVol/4.0;
987 interactionVolume.setNeumannCondition(boundValues, 3);
988 }
989 if (bcType.hasDirichlet())
990 {
991 problem_.dirichlet(boundValues, intersection);
992 interactionVolume.setDirichletCondition(boundValues, 3);
993 }
994 }
995 else if (indexInInside == 1)
996 {
997 problem_.boundaryTypes(bcType, intersection);
998 PrimaryVariables boundValues(0.0);
999
1000 interactionVolume.setBoundary(bcType, 2);
1001 if (bcType.isNeumann(pressureEqIdx))
1002 {
1003 problem_.neumann(boundValues, intersection);
1004 boundValues *= faceVol/4.0;
1005 interactionVolume.setNeumannCondition(boundValues, 2);
1006 }
1007 if (bcType.hasDirichlet())
1008 {
1009 problem_.dirichlet(boundValues, intersection);
1010 interactionVolume.setDirichletCondition(boundValues, 2);
1011 }
1012 }
1013 else if (indexInInside == 5)
1014 {
1015 problem_.boundaryTypes(bcType, intersection);
1016 PrimaryVariables boundValues(0.0);
1017
1018 interactionVolume.setBoundary(bcType, 11);
1019 if (bcType.isNeumann(pressureEqIdx))
1020 {
1021 problem_.neumann(boundValues, intersection);
1022 boundValues *= faceVol/4.0;
1023 interactionVolume.setNeumannCondition(boundValues, 11);
1024 }
1025 if (bcType.hasDirichlet())
1026 {
1027 problem_.dirichlet(boundValues, intersection);
1028 interactionVolume.setDirichletCondition(boundValues, 11);
1029 }
1030 }
1031 }
1032 if (elemVertMap[vIdxGlobal][3] == eIdxGlobal)
1033 {
1034 if (indexInInside == 0)
1035 {
1036 problem_.boundaryTypes(bcType, intersection);
1037 PrimaryVariables boundValues(0.0);
1038
1039 interactionVolume.setBoundary(bcType, 2);
1040 if (bcType.isNeumann(pressureEqIdx))
1041 {
1042 problem_.neumann(boundValues, intersection);
1043 boundValues *= faceVol/4.0;
1044 interactionVolume.setNeumannCondition(boundValues, 2);
1045 }
1046 if (bcType.hasDirichlet())
1047 {
1048 problem_.dirichlet(boundValues, intersection);
1049 interactionVolume.setDirichletCondition(boundValues, 2);
1050 }
1051 }
1052 else if (indexInInside == 2)
1053 {
1054 problem_.boundaryTypes(bcType, intersection);
1055 PrimaryVariables boundValues(0.0);
1056
1057 interactionVolume.setBoundary(bcType, 1);
1058 if (bcType.isNeumann(pressureEqIdx))
1059 {
1060 problem_.neumann(boundValues, intersection);
1061 boundValues *= faceVol/4.0;
1062 interactionVolume.setNeumannCondition(boundValues, 1);
1063 }
1064 if (bcType.hasDirichlet())
1065 {
1066 problem_.dirichlet(boundValues, intersection);
1067 interactionVolume.setDirichletCondition(boundValues, 1);
1068 }
1069 }
1070 else if (indexInInside == 5)
1071 {
1072 problem_.boundaryTypes(bcType, intersection);
1073 PrimaryVariables boundValues(0.0);
1074
1075 interactionVolume.setBoundary(bcType, 10);
1076 if (bcType.isNeumann(pressureEqIdx))
1077 {
1078 problem_.neumann(boundValues, intersection);
1079 boundValues *= faceVol/4.0;
1080 interactionVolume.setNeumannCondition(boundValues, 10);
1081 }
1082 if (bcType.hasDirichlet())
1083 {
1084 problem_.dirichlet(boundValues, intersection);
1085 interactionVolume.setDirichletCondition(boundValues, 10);
1086 }
1087 }
1088 }
1089 if (elemVertMap[vIdxGlobal][4] == eIdxGlobal)
1090 {
1091 if (indexInInside == 4)
1092 {
1093 problem_.boundaryTypes(bcType, intersection);
1094 PrimaryVariables boundValues(0.0);
1095
1096 interactionVolume.setBoundary(bcType, 8);
1097 if (bcType.isNeumann(pressureEqIdx))
1098 {
1099 problem_.neumann(boundValues, intersection);
1100 boundValues *= faceVol/4.0;
1101 interactionVolume.setNeumannCondition(boundValues, 8);
1102 }
1103 if (bcType.hasDirichlet())
1104 {
1105 problem_.dirichlet(boundValues, intersection);
1106 interactionVolume.setDirichletCondition(boundValues, 8);
1107 }
1108 }
1109 else if (indexInInside == 1)
1110 {
1111 problem_.boundaryTypes(bcType, intersection);
1112 PrimaryVariables boundValues(0.0);
1113
1114 interactionVolume.setBoundary(bcType, 4);
1115 if (bcType.isNeumann(pressureEqIdx))
1116 {
1117 problem_.neumann(boundValues, intersection);
1118 boundValues *= faceVol/4.0;
1119 interactionVolume.setNeumannCondition(boundValues, 4);
1120 }
1121 if (bcType.hasDirichlet())
1122 {
1123 problem_.dirichlet(boundValues, intersection);
1124 interactionVolume.setDirichletCondition(boundValues, 4);
1125 }
1126 }
1127 else if (indexInInside == 3)
1128 {
1129 problem_.boundaryTypes(bcType, intersection);
1130 PrimaryVariables boundValues(0.0);
1131
1132 interactionVolume.setBoundary(bcType, 7);
1133 if (bcType.isNeumann(pressureEqIdx))
1134 {
1135 problem_.neumann(boundValues, intersection);
1136 boundValues *= faceVol/4.0;
1137 interactionVolume.setNeumannCondition(boundValues, 7);
1138 }
1139 if (bcType.hasDirichlet())
1140 {
1141 problem_.dirichlet(boundValues, intersection);
1142 interactionVolume.setDirichletCondition(boundValues, 7);
1143 }
1144 }
1145 }
1146 if (elemVertMap[vIdxGlobal][5] == eIdxGlobal)
1147 {
1148 if (indexInInside == 4)
1149 {
1150 problem_.boundaryTypes(bcType, intersection);
1151 PrimaryVariables boundValues(0.0);
1152
1153 interactionVolume.setBoundary(bcType, 9);
1154 if (bcType.isNeumann(pressureEqIdx))
1155 {
1156 problem_.neumann(boundValues, intersection);
1157 boundValues *= faceVol/4.0;
1158 interactionVolume.setNeumannCondition(boundValues, 9);
1159 }
1160 if (bcType.hasDirichlet())
1161 {
1162 problem_.dirichlet(boundValues, intersection);
1163 interactionVolume.setDirichletCondition(boundValues, 9);
1164 }
1165 }
1166 else if (indexInInside == 3)
1167 {
1168 problem_.boundaryTypes(bcType, intersection);
1169 PrimaryVariables boundValues(0.0);
1170
1171 interactionVolume.setBoundary(bcType, 5);
1172 if (bcType.isNeumann(pressureEqIdx))
1173 {
1174 problem_.neumann(boundValues, intersection);
1175 boundValues *= faceVol/4.0;
1176 interactionVolume.setNeumannCondition(boundValues, 5);
1177 }
1178 if (bcType.hasDirichlet())
1179 {
1180 problem_.dirichlet(boundValues, intersection);
1181 interactionVolume.setDirichletCondition(boundValues, 5);
1182 }
1183 }
1184 else if (indexInInside == 0)
1185 {
1186 problem_.boundaryTypes(bcType, intersection);
1187 PrimaryVariables boundValues(0.0);
1188
1189 interactionVolume.setBoundary(bcType, 4);
1190 if (bcType.isNeumann(pressureEqIdx))
1191 {
1192 problem_.neumann(boundValues, intersection);
1193 boundValues *= faceVol/4.0;
1194 interactionVolume.setNeumannCondition(boundValues, 4);
1195 }
1196 if (bcType.hasDirichlet())
1197 {
1198 problem_.dirichlet(boundValues, intersection);
1199 interactionVolume.setDirichletCondition(boundValues, 4);
1200 }
1201 }
1202 }
1203 if (elemVertMap[vIdxGlobal][6] == eIdxGlobal)
1204 {
1205 if (indexInInside == 4)
1206 {
1207 problem_.boundaryTypes(bcType, intersection);
1208 PrimaryVariables boundValues(0.0);
1209
1210 interactionVolume.setBoundary(bcType, 11);
1211 if (bcType.isNeumann(pressureEqIdx))
1212 {
1213 problem_.neumann(boundValues, intersection);
1214 boundValues *= faceVol/4.0;
1215 interactionVolume.setNeumannCondition(boundValues, 11);
1216 }
1217 if (bcType.hasDirichlet())
1218 {
1219 problem_.dirichlet(boundValues, intersection);
1220 interactionVolume.setDirichletCondition(boundValues, 11);
1221 }
1222 }
1223 else if (indexInInside == 2)
1224 {
1225 problem_.boundaryTypes(bcType, intersection);
1226 PrimaryVariables boundValues(0.0);
1227
1228 interactionVolume.setBoundary(bcType, 7);
1229 if (bcType.isNeumann(pressureEqIdx))
1230 {
1231 problem_.neumann(boundValues, intersection);
1232 boundValues *= faceVol/4.0;
1233 interactionVolume.setNeumannCondition(boundValues, 7);
1234 }
1235 if (bcType.hasDirichlet())
1236 {
1237 problem_.dirichlet(boundValues, intersection);
1238 interactionVolume.setDirichletCondition(boundValues, 7);
1239 }
1240 }
1241 else if (indexInInside == 1)
1242 {
1243 problem_.boundaryTypes(bcType, intersection);
1244 PrimaryVariables boundValues(0.0);
1245
1246 interactionVolume.setBoundary(bcType, 6);
1247 if (bcType.isNeumann(pressureEqIdx))
1248 {
1249 problem_.neumann(boundValues, intersection);
1250 boundValues *= faceVol/4.0;
1251 interactionVolume.setNeumannCondition(boundValues, 6);
1252 }
1253 if (bcType.hasDirichlet())
1254 {
1255 problem_.dirichlet(boundValues, intersection);
1256 interactionVolume.setDirichletCondition(boundValues, 6);
1257 }
1258 }
1259 }
1260 if (elemVertMap[vIdxGlobal][7] == eIdxGlobal)
1261 {
1262 if (indexInInside == 4)
1263 {
1264 problem_.boundaryTypes(bcType, intersection);
1265 PrimaryVariables boundValues(0.0);
1266
1267 interactionVolume.setBoundary(bcType, 10);
1268 if (bcType.isNeumann(pressureEqIdx))
1269 {
1270 problem_.neumann(boundValues, intersection);
1271 boundValues *= faceVol/4.0;
1272 interactionVolume.setNeumannCondition(boundValues, 10);
1273 }
1274 if (bcType.hasDirichlet())
1275 {
1276 problem_.dirichlet(boundValues, intersection);
1277 interactionVolume.setDirichletCondition(boundValues, 10);
1278 }
1279 }
1280 else if (indexInInside == 0)
1281 {
1282 problem_.boundaryTypes(bcType, intersection);
1283 PrimaryVariables boundValues(0.0);
1284
1285 interactionVolume.setBoundary(bcType, 6);
1286 if (bcType.isNeumann(pressureEqIdx))
1287 {
1288 problem_.neumann(boundValues, intersection);
1289 boundValues *= faceVol/4.0;
1290 interactionVolume.setNeumannCondition(boundValues, 6);
1291 }
1292 if (bcType.hasDirichlet())
1293 {
1294 problem_.dirichlet(boundValues, intersection);
1295 interactionVolume.setDirichletCondition(boundValues, 6);
1296 }
1297 }
1298 else if (indexInInside == 2)
1299 {
1300 problem_.boundaryTypes(bcType, intersection);
1301 PrimaryVariables boundValues(0.0);
1302
1303 interactionVolume.setBoundary(bcType, 5);
1304 if (bcType.isNeumann(pressureEqIdx))
1305 {
1306 problem_.neumann(boundValues, intersection);
1307 boundValues *= faceVol/4.0;
1308 interactionVolume.setNeumannCondition(boundValues, 5);
1309 }
1310 if (bcType.hasDirichlet())
1311 {
1312 problem_.dirichlet(boundValues, intersection);
1313 interactionVolume.setDirichletCondition(boundValues, 5);
1314 }
1315 }
1316 }
1317 }
1318 }
1319 }
1320
1321 }
1322}
1323
1340template<class TypeTag>
1341void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeInnerInteractionVolume(InteractionVolume& interactionVolume,
1342 const Vertex& vertex, bool sameLevel)
1343{
1344 const DimVector& centerPos = vertex.geometry().center();
1345
1346 interactionVolume.setCenterPosition(centerPos);
1347
1348 if (sameLevel)
1349 {
1350 auto element1 = interactionVolume.getSubVolumeElement(0);
1351 auto element8 = interactionVolume.getSubVolumeElement(7);
1352
1353 const ElementGeometry& geometry1 = element1.geometry();
1354 const ElementGeometry& geometry8 = element8.geometry();
1355
1356 const auto referenceElement = ReferenceElements::general(geometry1.type());
1357
1358 DimVector edgeCoord(geometry1.global(referenceElement.position(9, dim - 1)));
1359 interactionVolume.setEdgePosition(edgeCoord, 2);
1360 edgeCoord = geometry1.global(referenceElement.position(3, dim - 1));
1361 interactionVolume.setEdgePosition(edgeCoord, 0);
1362 edgeCoord = geometry1.global(referenceElement.position(11, dim - 1));
1363 interactionVolume.setEdgePosition(edgeCoord, 5);
1364
1365 edgeCoord = geometry8.global(referenceElement.position(4, dim - 1));
1366 interactionVolume.setEdgePosition(edgeCoord, 4);
1367 edgeCoord = geometry8.global(referenceElement.position(6, dim - 1));
1368 interactionVolume.setEdgePosition(edgeCoord, 3);
1369 edgeCoord = geometry8.global(referenceElement.position(0, dim - 1));
1370 interactionVolume.setEdgePosition(edgeCoord, 1);
1371 }
1372
1373 DimVector edgeCoord1(interactionVolume.getEdgePosition(0));
1374 DimVector edgeCoord2(interactionVolume.getEdgePosition(1));
1375 DimVector edgeCoord3(interactionVolume.getEdgePosition(2));
1376 DimVector edgeCoord4(interactionVolume.getEdgePosition(3));
1377 DimVector edgeCoord5(interactionVolume.getEdgePosition(4));
1378 DimVector edgeCoord6(interactionVolume.getEdgePosition(5));
1379
1380 DimVector crossProductVector1(0);
1381 DimVector crossProductVector2(0);
1382
1383 GlobalPosition globalPosFace = interactionVolume.getFacePosition(0);
1384 crossProductVector1 = centerPos - globalPosFace;
1385 crossProductVector2 = edgeCoord1 - edgeCoord3;
1386 Scalar faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1387 interactionVolume.setFaceArea(faceArea, 0);
1388
1389 globalPosFace = interactionVolume.getFacePosition(1);
1390 crossProductVector1 = centerPos - globalPosFace;
1391 crossProductVector2 = edgeCoord1 - edgeCoord4;
1392
1393 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1394
1395 interactionVolume.setFaceArea(faceArea, 1);
1396
1397 globalPosFace = interactionVolume.getFacePosition(2);
1398 crossProductVector1 = centerPos - globalPosFace;
1399 crossProductVector2 = edgeCoord1 - edgeCoord5;
1400 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1401 interactionVolume.setFaceArea(faceArea, 2);
1402
1403 globalPosFace = interactionVolume.getFacePosition(3);
1404 crossProductVector1 = centerPos - globalPosFace;
1405 crossProductVector2 = edgeCoord1 - edgeCoord6;
1406 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1407 interactionVolume.setFaceArea(faceArea, 3);
1408
1409 globalPosFace = interactionVolume.getFacePosition(4);
1410 crossProductVector1 = centerPos - globalPosFace;
1411 crossProductVector2 = edgeCoord2 - edgeCoord3;
1412 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1413 interactionVolume.setFaceArea(faceArea, 4);
1414
1415 globalPosFace = interactionVolume.getFacePosition(5);
1416 crossProductVector1 = centerPos - globalPosFace;
1417 crossProductVector2 = edgeCoord2 - edgeCoord4;
1418 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1419 interactionVolume.setFaceArea(faceArea, 5);
1420
1421 globalPosFace = interactionVolume.getFacePosition(6);
1422 crossProductVector1 = centerPos - globalPosFace;
1423 crossProductVector2 = edgeCoord2 - edgeCoord5;
1424 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1425 interactionVolume.setFaceArea(faceArea, 6);
1426
1427 globalPosFace = interactionVolume.getFacePosition(7);
1428 crossProductVector1 = centerPos - globalPosFace;
1429 crossProductVector2 = edgeCoord2 - edgeCoord6;
1430 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1431 interactionVolume.setFaceArea(faceArea, 7);
1432
1433 globalPosFace = interactionVolume.getFacePosition(8);
1434 crossProductVector1 = centerPos - globalPosFace;
1435 crossProductVector2 = edgeCoord3 - edgeCoord6;
1436 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1437 interactionVolume.setFaceArea(faceArea, 8);
1438
1439 globalPosFace = interactionVolume.getFacePosition(9);
1440 crossProductVector1 = centerPos - globalPosFace;
1441 crossProductVector2 = edgeCoord3 - edgeCoord4;
1442 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1443 interactionVolume.setFaceArea(faceArea, 9);
1444
1445 globalPosFace = interactionVolume.getFacePosition(10);
1446 crossProductVector1 = centerPos - globalPosFace;
1447 crossProductVector2 = edgeCoord4 - edgeCoord5;
1448 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1449 interactionVolume.setFaceArea(faceArea, 10);
1450
1451 globalPosFace = interactionVolume.getFacePosition(11);
1452 crossProductVector1 = centerPos - globalPosFace;
1453 crossProductVector2 = edgeCoord5 - edgeCoord6;
1454 faceArea = crossProduct(crossProductVector1, crossProductVector2).two_norm()/2.0;
1455 interactionVolume.setFaceArea(faceArea, 11);
1456}
1457
1477template<class TypeTag>
1478void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeBoundaryInteractionVolume(InteractionVolume& interactionVolume,
1479 const Vertex& vertex)
1480{
1481 const DimVector& centerPos = vertex.geometry().center();
1482
1483 interactionVolume.setCenterPosition(centerPos);
1484
1485 //corner
1486 switch (interactionVolume.getElementNumber())
1487 {
1488 case 1:
1489 {
1490 if (interactionVolume.hasSubVolumeElement(0))
1491 {
1492 interactionVolume.setOutsideFace(1);
1493 interactionVolume.setOutsideFace(2);
1494 interactionVolume.setOutsideFace(4);
1495 interactionVolume.setOutsideFace(5);
1496 interactionVolume.setOutsideFace(6);
1497 interactionVolume.setOutsideFace(7);
1498 interactionVolume.setOutsideFace(9);
1499 interactionVolume.setOutsideFace(10);
1500 interactionVolume.setOutsideFace(11);
1501
1502 int eIdxGlobal = problem_.variables().index(interactionVolume.getSubVolumeElement(0));
1503 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 0, 0)/4.0,0);
1504 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 0, 1)/4.0,3);
1505 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 0, 2)/4.0,8);
1506 }
1507 if (interactionVolume.hasSubVolumeElement(1))
1508 {
1509 interactionVolume.setOutsideFace(2);
1510 interactionVolume.setOutsideFace(3);
1511 interactionVolume.setOutsideFace(4);
1512 interactionVolume.setOutsideFace(5);
1513 interactionVolume.setOutsideFace(6);
1514 interactionVolume.setOutsideFace(7);
1515 interactionVolume.setOutsideFace(8);
1516 interactionVolume.setOutsideFace(10);
1517 interactionVolume.setOutsideFace(11);
1518
1519 int eIdxGlobal = problem_.variables().index(interactionVolume.getSubVolumeElement(1));
1520 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 1, 0)/4.0,1);
1521 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 1, 1)/4.0,0);
1522 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 1, 2)/4.0,9);
1523 }
1524 if (interactionVolume.hasSubVolumeElement(2))
1525 {
1526 interactionVolume.setOutsideFace(0);
1527 interactionVolume.setOutsideFace(1);
1528 interactionVolume.setOutsideFace(4);
1529 interactionVolume.setOutsideFace(5);
1530 interactionVolume.setOutsideFace(6);
1531 interactionVolume.setOutsideFace(7);
1532 interactionVolume.setOutsideFace(8);
1533 interactionVolume.setOutsideFace(9);
1534 interactionVolume.setOutsideFace(10);
1535
1536 int eIdxGlobal = problem_.variables().index(interactionVolume.getSubVolumeElement(2));
1537 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 2, 0)/4.0,3);
1538 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 2, 1)/4.0,2);
1539 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 2, 2)/4.0,11);
1540 }
1541 if (interactionVolume.hasSubVolumeElement(3))
1542 {
1543 interactionVolume.setOutsideFace(0);
1544 interactionVolume.setOutsideFace(3);
1545 interactionVolume.setOutsideFace(4);
1546 interactionVolume.setOutsideFace(5);
1547 interactionVolume.setOutsideFace(6);
1548 interactionVolume.setOutsideFace(7);
1549 interactionVolume.setOutsideFace(8);
1550 interactionVolume.setOutsideFace(9);
1551 interactionVolume.setOutsideFace(11);
1552
1553 int eIdxGlobal = problem_.variables().index(interactionVolume.getSubVolumeElement(3));
1554 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 3, 0)/4.0,2);
1555 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 3, 1)/4.0,1);
1556 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 3, 2)/4.0,10);
1557 }
1558 if (interactionVolume.hasSubVolumeElement(4))
1559 {
1560 interactionVolume.setOutsideFace(0);
1561 interactionVolume.setOutsideFace(1);
1562 interactionVolume.setOutsideFace(2);
1563 interactionVolume.setOutsideFace(3);
1564 interactionVolume.setOutsideFace(5);
1565 interactionVolume.setOutsideFace(6);
1566 interactionVolume.setOutsideFace(9);
1567 interactionVolume.setOutsideFace(10);
1568 interactionVolume.setOutsideFace(11);
1569
1570 int eIdxGlobal = problem_.variables().index(interactionVolume.getSubVolumeElement(4));
1571 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 4, 0)/4.0,8);
1572 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 4, 1)/4.0,4);
1573 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 4, 2)/4.0,7);
1574 }
1575 if (interactionVolume.hasSubVolumeElement(5))
1576 {
1577 interactionVolume.setOutsideFace(0);
1578 interactionVolume.setOutsideFace(1);
1579 interactionVolume.setOutsideFace(2);
1580 interactionVolume.setOutsideFace(3);
1581 interactionVolume.setOutsideFace(6);
1582 interactionVolume.setOutsideFace(7);
1583 interactionVolume.setOutsideFace(8);
1584 interactionVolume.setOutsideFace(10);
1585 interactionVolume.setOutsideFace(11);
1586
1587 int eIdxGlobal = problem_.variables().index(interactionVolume.getSubVolumeElement(5));
1588 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 5, 0)/4.0,9);
1589 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 5, 1)/4.0,5);
1590 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 5, 2)/4.0,4);
1591 }
1592 if (interactionVolume.hasSubVolumeElement(6))
1593 {
1594 interactionVolume.setOutsideFace(0);
1595 interactionVolume.setOutsideFace(1);
1596 interactionVolume.setOutsideFace(2);
1597 interactionVolume.setOutsideFace(3);
1598 interactionVolume.setOutsideFace(4);
1599 interactionVolume.setOutsideFace(5);
1600 interactionVolume.setOutsideFace(8);
1601 interactionVolume.setOutsideFace(9);
1602 interactionVolume.setOutsideFace(10);
1603
1604 int eIdxGlobal = problem_.variables().index(interactionVolume.getSubVolumeElement(6));
1605 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 6, 0)/4.0,11);
1606 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 6, 1)/4.0,7);
1607 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 6, 2)/4.0,6);
1608 }
1609 if (interactionVolume.hasSubVolumeElement(7))
1610 {
1611 interactionVolume.setOutsideFace(0);
1612 interactionVolume.setOutsideFace(1);
1613 interactionVolume.setOutsideFace(2);
1614 interactionVolume.setOutsideFace(3);
1615 interactionVolume.setOutsideFace(4);
1616 interactionVolume.setOutsideFace(7);
1617 interactionVolume.setOutsideFace(8);
1618 interactionVolume.setOutsideFace(9);
1619 interactionVolume.setOutsideFace(11);
1620
1621 int eIdxGlobal = problem_.variables().index(interactionVolume.getSubVolumeElement(7));
1622 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 7, 0)/4.0,10);
1623 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 7, 1)/4.0,6);
1624 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal, 7, 2)/4.0,5);
1625 }
1626 break;
1627 }
1628 case 2:
1629 {
1630 // edge
1631 if (interactionVolume.hasSubVolumeElement(0))
1632 {
1633 int eIdxGlobal1 = problem_.variables().index(interactionVolume.getSubVolumeElement(0));
1634 if (interactionVolume.hasSubVolumeElement(1))
1635 {
1636 interactionVolume.setOutsideFace(4);
1637 interactionVolume.setOutsideFace(5);
1638 interactionVolume.setOutsideFace(6);
1639 interactionVolume.setOutsideFace(7);
1640 interactionVolume.setOutsideFace(2);
1641 interactionVolume.setOutsideFace(10);
1642 interactionVolume.setOutsideFace(11);
1643
1644 int eIdxGlobal2 = problem_.variables().index(interactionVolume.getSubVolumeElement(1));
1645 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 0, 1)/4.0,3);
1646 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 0, 2)/4.0,8);
1647 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 1, 0)/4.0,1);
1648 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 1, 2)/4.0,9);
1649
1650 return;
1651 }
1652 if (interactionVolume.hasSubVolumeElement(2))
1653 {
1654 interactionVolume.setOutsideFace(4);
1655 interactionVolume.setOutsideFace(5);
1656 interactionVolume.setOutsideFace(6);
1657 interactionVolume.setOutsideFace(7);
1658 interactionVolume.setOutsideFace(1);
1659 interactionVolume.setOutsideFace(9);
1660 interactionVolume.setOutsideFace(10);
1661
1662 int eIdxGlobal2 = problem_.variables().index(interactionVolume.getSubVolumeElement(2));
1663 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 0, 0)/4.0,0);
1664 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 0, 2)/4.0,8);
1665 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 2, 1)/4.0,2);
1666 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 2, 2)/4.0,11);
1667
1668 return;
1669 }
1670 if (interactionVolume.hasSubVolumeElement(4))
1671 {
1672 interactionVolume.setOutsideFace(1);
1673 interactionVolume.setOutsideFace(5);
1674 interactionVolume.setOutsideFace(9);
1675 interactionVolume.setOutsideFace(10);
1676 interactionVolume.setOutsideFace(2);
1677 interactionVolume.setOutsideFace(6);
1678 interactionVolume.setOutsideFace(11);
1679
1680 int eIdxGlobal2 = problem_.variables().index(interactionVolume.getSubVolumeElement(4));
1681 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 0, 0)/4.0,0);
1682 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 0, 1)/4.0,3);
1683 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 4, 1)/4.0,4);
1684 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 4, 2)/4.0,7);
1685
1686 return;
1687 }
1688 }
1689 if (interactionVolume.hasSubVolumeElement(7))
1690 {
1691 int eIdxGlobal1 = problem_.variables().index(interactionVolume.getSubVolumeElement(7));
1692 if (interactionVolume.hasSubVolumeElement(5))
1693 {
1694 interactionVolume.setOutsideFace(0);
1695 interactionVolume.setOutsideFace(1);
1696 interactionVolume.setOutsideFace(2);
1697 interactionVolume.setOutsideFace(3);
1698 interactionVolume.setOutsideFace(7);
1699 interactionVolume.setOutsideFace(8);
1700 interactionVolume.setOutsideFace(11);
1701
1702 int eIdxGlobal2 = problem_.variables().index(interactionVolume.getSubVolumeElement(5));
1703 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 7, 0)/4.0,10);
1704 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 7, 1)/4.0,6);
1705 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 5, 0)/4.0,9);
1706 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 5, 2)/4.0,4);
1707
1708 return;
1709 }
1710 if (interactionVolume.hasSubVolumeElement(6))
1711 {
1712 interactionVolume.setOutsideFace(0);
1713 interactionVolume.setOutsideFace(1);
1714 interactionVolume.setOutsideFace(2);
1715 interactionVolume.setOutsideFace(3);
1716 interactionVolume.setOutsideFace(4);
1717 interactionVolume.setOutsideFace(8);
1718 interactionVolume.setOutsideFace(9);
1719
1720 int eIdxGlobal2 = problem_.variables().index(interactionVolume.getSubVolumeElement(6));
1721 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 7, 0)/4.0,10);
1722 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 7, 2)/4.0,5);
1723 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 6, 0)/4.0,11);
1724 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 6, 1)/4.0,7);
1725
1726 return;
1727 }
1728 if (interactionVolume.hasSubVolumeElement(3))
1729 {
1730 interactionVolume.setOutsideFace(0);
1731 interactionVolume.setOutsideFace(4);
1732 interactionVolume.setOutsideFace(8);
1733 interactionVolume.setOutsideFace(9);
1734 interactionVolume.setOutsideFace(3);
1735 interactionVolume.setOutsideFace(7);
1736 interactionVolume.setOutsideFace(11);
1737
1738 int eIdxGlobal2 = problem_.variables().index(interactionVolume.getSubVolumeElement(3));
1739 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 7, 1)/4.0,6);
1740 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 7, 2)/4.0,5);
1741 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 3, 0)/4.0,2);
1742 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 3, 1)/4.0,1);
1743
1744 return;
1745 }
1746 }
1747 if (interactionVolume.hasSubVolumeElement(5))
1748 {
1749 int eIdxGlobal1 = problem_.variables().index(interactionVolume.getSubVolumeElement(5));
1750 if (interactionVolume.hasSubVolumeElement(1))
1751 {
1752 interactionVolume.setOutsideFace(3);
1753 interactionVolume.setOutsideFace(7);
1754 interactionVolume.setOutsideFace(8);
1755 interactionVolume.setOutsideFace(11);
1756 interactionVolume.setOutsideFace(2);
1757 interactionVolume.setOutsideFace(6);
1758 interactionVolume.setOutsideFace(10);
1759
1760 int eIdxGlobal2 = problem_.variables().index(interactionVolume.getSubVolumeElement(1));
1761 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 5, 1)/4.0,5);
1762 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 5, 2)/4.0,4);
1763 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 1, 0)/4.0,1);
1764 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 1, 1)/4.0,0);
1765
1766 return;
1767 }
1768 if (interactionVolume.hasSubVolumeElement(4))
1769 {
1770 interactionVolume.setOutsideFace(0);
1771 interactionVolume.setOutsideFace(1);
1772 interactionVolume.setOutsideFace(2);
1773 interactionVolume.setOutsideFace(3);
1774 interactionVolume.setOutsideFace(6);
1775 interactionVolume.setOutsideFace(10);
1776 interactionVolume.setOutsideFace(11);
1777
1778 int eIdxGlobal2 = problem_.variables().index(interactionVolume.getSubVolumeElement(4));
1779 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 5, 0)/4.0,9);
1780 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 5, 1)/4.0,5);
1781 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 4, 0)/4.0,8);
1782 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 4, 2)/4.0,7);
1783
1784 return;
1785 }
1786 }
1787 if (interactionVolume.hasSubVolumeElement(6))
1788 {
1789 int eIdxGlobal1 = problem_.variables().index(interactionVolume.getSubVolumeElement(6));
1790 if (interactionVolume.hasSubVolumeElement(4))
1791 {
1792 interactionVolume.setOutsideFace(1);
1793 interactionVolume.setOutsideFace(5);
1794 interactionVolume.setOutsideFace(9);
1795 interactionVolume.setOutsideFace(10);
1796 interactionVolume.setOutsideFace(0);
1797 interactionVolume.setOutsideFace(2);
1798 interactionVolume.setOutsideFace(3);
1799
1800 int eIdxGlobal2 = problem_.variables().index(interactionVolume.getSubVolumeElement(4));
1801 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 6, 0)/4.0,11);
1802 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 6, 2)/4.0,6);
1803 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 4, 0)/4.0,8);
1804 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 4, 1)/4.0,4);
1805
1806 return;
1807 }
1808 if (interactionVolume.hasSubVolumeElement(2))
1809 {
1810 interactionVolume.setOutsideFace(1);
1811 interactionVolume.setOutsideFace(5);
1812 interactionVolume.setOutsideFace(9);
1813 interactionVolume.setOutsideFace(10);
1814 interactionVolume.setOutsideFace(0);
1815 interactionVolume.setOutsideFace(4);
1816 interactionVolume.setOutsideFace(8);
1817
1818 int eIdxGlobal2 = problem_.variables().index(interactionVolume.getSubVolumeElement(2));
1819 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 6, 1)/4.0,7);
1820 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 6, 2)/4.0,6);
1821 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 2, 0)/4.0,3);
1822 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 2, 1)/4.0,2);
1823
1824 return;
1825 }
1826 }
1827 if (interactionVolume.hasSubVolumeElement(3))
1828 {
1829 int eIdxGlobal1 = problem_.variables().index(interactionVolume.getSubVolumeElement(3));
1830 if (interactionVolume.hasSubVolumeElement(1))
1831 {
1832 interactionVolume.setOutsideFace(4);
1833 interactionVolume.setOutsideFace(5);
1834 interactionVolume.setOutsideFace(6);
1835 interactionVolume.setOutsideFace(7);
1836 interactionVolume.setOutsideFace(3);
1837 interactionVolume.setOutsideFace(8);
1838 interactionVolume.setOutsideFace(11);
1839
1840 int eIdxGlobal2 = problem_.variables().index(interactionVolume.getSubVolumeElement(1));
1841 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 3, 0)/4.0,2);
1842 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 3, 2)/4.0,10);
1843 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 1, 1)/4.0,0);
1844 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 1, 2)/4.0,9);
1845
1846 return;
1847 }
1848 if (interactionVolume.hasSubVolumeElement(2))
1849 {
1850 interactionVolume.setOutsideFace(4);
1851 interactionVolume.setOutsideFace(5);
1852 interactionVolume.setOutsideFace(6);
1853 interactionVolume.setOutsideFace(7);
1854 interactionVolume.setOutsideFace(0);
1855 interactionVolume.setOutsideFace(8);
1856 interactionVolume.setOutsideFace(9);
1857
1858 int eIdxGlobal2 = problem_.variables().index(interactionVolume.getSubVolumeElement(2));
1859 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 3, 1)/4.0,1);
1860 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 3, 2)/4.0,10);
1861 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 2, 0)/4.0,3);
1862 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 2, 2)/4.0,11);
1863
1864 return;
1865 }
1866 }
1867 break;
1868 }
1869 case 4:
1870 {
1871 //side
1872 if (interactionVolume.hasSubVolumeElement(0))
1873 {
1874 int eIdxGlobal1 = problem_.variables().index(interactionVolume.getSubVolumeElement(0));
1875 if (interactionVolume.hasSubVolumeElement(1))
1876 {
1877 int eIdxGlobal2 = problem_.variables().index(interactionVolume.getSubVolumeElement(1));
1878 if (interactionVolume.hasSubVolumeElement(2) && interactionVolume.hasSubVolumeElement(3))
1879 {
1880 interactionVolume.setOutsideFace(4);
1881 interactionVolume.setOutsideFace(5);
1882 interactionVolume.setOutsideFace(6);
1883 interactionVolume.setOutsideFace(7);
1884
1885 int eIdxGlobal3 = problem_.variables().index(interactionVolume.getSubVolumeElement(2));
1886 int eIdxGlobal4 = problem_.variables().index(interactionVolume.getSubVolumeElement(3));
1887 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 0, 2)/4.0, 8);
1888 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 1, 2)/4.0, 9);
1889 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal3, 2, 2)/4.0, 10);
1890 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal4, 3, 2)/4.0, 11);
1891
1892 return;
1893 }
1894 if (interactionVolume.hasSubVolumeElement(4) && interactionVolume.hasSubVolumeElement(5))
1895 {
1896 interactionVolume.setOutsideFace(2);
1897 interactionVolume.setOutsideFace(6);
1898 interactionVolume.setOutsideFace(10);
1899 interactionVolume.setOutsideFace(11);
1900
1901 int eIdxGlobal3 = problem_.variables().index(interactionVolume.getSubVolumeElement(4));
1902 int eIdxGlobal4 = problem_.variables().index(interactionVolume.getSubVolumeElement(5));
1903 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 0, 1)/4.0, 3);
1904 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 1, 0)/4.0, 1);
1905 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal3, 4, 2)/4.0, 7);
1906 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal4, 5, 1)/4.0, 5);
1907
1908 return;
1909 }
1910 }
1911 if (interactionVolume.hasSubVolumeElement(2) && interactionVolume.hasSubVolumeElement(4)
1912 && interactionVolume.hasSubVolumeElement(6))
1913 {
1914 interactionVolume.setOutsideFace(1);
1915 interactionVolume.setOutsideFace(5);
1916 interactionVolume.setOutsideFace(9);
1917 interactionVolume.setOutsideFace(10);
1918
1919 int eIdxGlobal2 = problem_.variables().index(interactionVolume.getSubVolumeElement(2));
1920 int eIdxGlobal3 = problem_.variables().index(interactionVolume.getSubVolumeElement(4));
1921 int eIdxGlobal4 = problem_.variables().index(interactionVolume.getSubVolumeElement(6));
1922 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 0, 0)/4.0, 0);
1923 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 2, 1)/4.0, 2);
1924 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal3, 4, 1)/4.0, 4);
1925 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal4, 6, 2)/4.0, 6);
1926
1927 return;
1928 }
1929 }
1930 if (interactionVolume.hasSubVolumeElement(7))
1931 {
1932 int eIdxGlobal1 = problem_.variables().index(interactionVolume.getSubVolumeElement(7));
1933 if (interactionVolume.hasSubVolumeElement(5))
1934 {
1935 int eIdxGlobal2 = problem_.variables().index(interactionVolume.getSubVolumeElement(5));
1936 if (interactionVolume.hasSubVolumeElement(1) && interactionVolume.hasSubVolumeElement(3))
1937 {
1938 interactionVolume.setOutsideFace(3);
1939 interactionVolume.setOutsideFace(7);
1940 interactionVolume.setOutsideFace(8);
1941 interactionVolume.setOutsideFace(11);
1942
1943 int eIdxGlobal3 = problem_.variables().index(interactionVolume.getSubVolumeElement(1));
1944 int eIdxGlobal4 = problem_.variables().index(interactionVolume.getSubVolumeElement(3));
1945 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 7, 1)/4.0, 6);
1946 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 5, 2)/4.0, 4);
1947 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal3, 1, 1)/4.0, 0);
1948 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal4, 3, 0)/4.0, 2);
1949
1950 return;
1951 }
1952 if (interactionVolume.hasSubVolumeElement(4) && interactionVolume.hasSubVolumeElement(6))
1953 {
1954 interactionVolume.setOutsideFace(0);
1955 interactionVolume.setOutsideFace(1);
1956 interactionVolume.setOutsideFace(2);
1957 interactionVolume.setOutsideFace(3);
1958
1959 int eIdxGlobal3 = problem_.variables().index(interactionVolume.getSubVolumeElement(4));
1960 int eIdxGlobal4 = problem_.variables().index(interactionVolume.getSubVolumeElement(6));
1961 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 7, 0)/4.0, 10);
1962 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 5, 0)/4.0, 9);
1963 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal3, 4, 0)/4.0, 8);
1964 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal4, 6, 0)/4.0, 11);
1965
1966 return;
1967 }
1968 }
1969 if (interactionVolume.hasSubVolumeElement(6) && interactionVolume.hasSubVolumeElement(2)
1970 && interactionVolume.hasSubVolumeElement(3))
1971 {
1972 interactionVolume.setOutsideFace(0);
1973 interactionVolume.setOutsideFace(4);
1974 interactionVolume.setOutsideFace(8);
1975 interactionVolume.setOutsideFace(9);
1976
1977 int eIdxGlobal2 = problem_.variables().index(interactionVolume.getSubVolumeElement(6));
1978 int eIdxGlobal3 = problem_.variables().index(interactionVolume.getSubVolumeElement(2));
1979 int eIdxGlobal4 = problem_.variables().index(interactionVolume.getSubVolumeElement(3));
1980 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal1, 7, 2)/4.0, 5);
1981 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal2, 6, 1)/4.0, 7);
1982 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal3, 2, 0)/4.0, 3);
1983 interactionVolume.setFaceArea(getRealFaceArea(interactionVolume, eIdxGlobal4, 3, 1)/4.0, 1);
1984
1985 return;
1986 }
1987 }
1988 break;
1989 }
1990 default:
1991 DUNE_THROW(Dune::NotImplemented, "Boundary shape not implemented");
1992 break;
1993 }
1994}
1995
1996
2000template<class TypeTag>
2001void FvMpfaL3dInteractionVolumeContainer<TypeTag>::storeInteractionVolumeInfo()
2002{
2003 std::vector < std::vector<int> > elemVertMap(problem_.gridView().size(dim), std::vector<int>(8, -1));
2004
2005 //Add elements to the interaction volumes and store element-vertex map
2006 for (const auto& element : elements(problem_.gridView()))
2007 storeSubVolumeElements(element, elemVertMap);
2008
2009 for (unsigned int i = 0; i < interactionVolumes_.size(); i++)
2010 if (interactionVolumes_[i].getElementNumber() == 0)
2011 interactionVolumes_[i].printInteractionVolumeInfo();
2012
2013 // Store information related to DUNE intersections for all interaction volumes
2014 for (const auto& element : elements(problem_.gridView()))
2015 storeIntersectionInfo(element, elemVertMap);
2016
2017 // Complete storage of the interaction volumes using the previously stored information
2018 // about the orientation and relationship of the DUNE elements in the interaction volumes (see doc/docextra/3dmpfa)
2019 for (const auto& vertex : vertices(problem_.gridView()))
2020 {
2021 int vIdxGlobal = problem_.variables().index(vertex);
2022
2023 InteractionVolume& interactionVolume = interactionVolumes_[vIdxGlobal];
2024
2025 if (interactionVolume.getElementNumber() == 8)
2026 {
2027 storeInnerInteractionVolume(interactionVolume, vertex);
2028 }
2029 else if (interactionVolume.isBoundaryInteractionVolume())
2030 {
2031 storeBoundaryInteractionVolume(interactionVolume, vertex);
2032 }
2033 else
2034 {
2035 DUNE_THROW(Dune::NotImplemented,"Interaction volume is no boundary volume but consists of less than 8 elements");
2036 }
2037
2038 // Store information about the MPFA flux face areas for correct flux weighting of
2039 // fluxes though faces which intersect the domain boundary
2040 // (see M. Wolff, Y. Cao, B. Flemisch, R. Helmig, and B. Wohlmuth (2013a). Multi-point flux
2041 // approximation L-method in 3D: numerical convergence and application to two-phase
2042 // flow through porous media. In P. Bastian, J. Kraus, R. Scheichl, and M. Wheeler,
2043 // editors, Simulation of Flow in Porous Media - Applications in Energy and Environment. De Gruyter.)
2044 if (!interactionVolume.isBoundaryInteractionVolume())
2045 {
2046 auto element1 = interactionVolume.getSubVolumeElement(0);
2047 auto element2 = interactionVolume.getSubVolumeElement(1);
2048 auto element3 = interactionVolume.getSubVolumeElement(2);
2049 auto element4 = interactionVolume.getSubVolumeElement(3);
2050 auto element5 = interactionVolume.getSubVolumeElement(4);
2051 auto element6 = interactionVolume.getSubVolumeElement(5);
2052 auto element7 = interactionVolume.getSubVolumeElement(6);
2053 auto element8 = interactionVolume.getSubVolumeElement(7);
2054
2055 int eIdxGlobal1 = problem_.variables().index(element1);
2056 int eIdxGlobal2 = problem_.variables().index(element2);
2057 int eIdxGlobal3 = problem_.variables().index(element3);
2058 int eIdxGlobal4 = problem_.variables().index(element4);
2059 int eIdxGlobal5 = problem_.variables().index(element5);
2060 int eIdxGlobal6 = problem_.variables().index(element6);
2061 int eIdxGlobal7 = problem_.variables().index(element7);
2062 int eIdxGlobal8 = problem_.variables().index(element8);
2063
2064 addRealFluxFaceArea_(interactionVolume.getFaceArea(0, 0), eIdxGlobal1, interactionVolume.getIndexOnElement(0, 0));
2065 addRealFluxFaceArea_(interactionVolume.getFaceArea(0, 1), eIdxGlobal1, interactionVolume.getIndexOnElement(0, 1));
2066 addRealFluxFaceArea_(interactionVolume.getFaceArea(0, 2), eIdxGlobal1, interactionVolume.getIndexOnElement(0, 2));
2067 addRealFluxFaceArea_(interactionVolume.getFaceArea(1, 0), eIdxGlobal2, interactionVolume.getIndexOnElement(1, 0));
2068 addRealFluxFaceArea_(interactionVolume.getFaceArea(1, 1), eIdxGlobal2, interactionVolume.getIndexOnElement(1, 1));
2069 addRealFluxFaceArea_(interactionVolume.getFaceArea(1, 2), eIdxGlobal2, interactionVolume.getIndexOnElement(1, 2));
2070 addRealFluxFaceArea_(interactionVolume.getFaceArea(2, 0), eIdxGlobal3, interactionVolume.getIndexOnElement(2, 0));
2071 addRealFluxFaceArea_(interactionVolume.getFaceArea(2, 1), eIdxGlobal3, interactionVolume.getIndexOnElement(2, 1));
2072 addRealFluxFaceArea_(interactionVolume.getFaceArea(2, 2), eIdxGlobal3, interactionVolume.getIndexOnElement(2, 2));
2073 addRealFluxFaceArea_(interactionVolume.getFaceArea(3, 0), eIdxGlobal4, interactionVolume.getIndexOnElement(3, 0));
2074 addRealFluxFaceArea_(interactionVolume.getFaceArea(3, 1), eIdxGlobal4, interactionVolume.getIndexOnElement(3, 1));
2075 addRealFluxFaceArea_(interactionVolume.getFaceArea(3, 2), eIdxGlobal4, interactionVolume.getIndexOnElement(3, 2));
2076 addRealFluxFaceArea_(interactionVolume.getFaceArea(4, 0), eIdxGlobal5, interactionVolume.getIndexOnElement(4, 0));
2077 addRealFluxFaceArea_(interactionVolume.getFaceArea(4, 1), eIdxGlobal5, interactionVolume.getIndexOnElement(4, 1));
2078 addRealFluxFaceArea_(interactionVolume.getFaceArea(4, 2), eIdxGlobal5, interactionVolume.getIndexOnElement(4, 2));
2079 addRealFluxFaceArea_(interactionVolume.getFaceArea(5, 0), eIdxGlobal6, interactionVolume.getIndexOnElement(5, 0));
2080 addRealFluxFaceArea_(interactionVolume.getFaceArea(5, 1), eIdxGlobal6, interactionVolume.getIndexOnElement(5, 1));
2081 addRealFluxFaceArea_(interactionVolume.getFaceArea(5, 2), eIdxGlobal6, interactionVolume.getIndexOnElement(5, 2));
2082 addRealFluxFaceArea_(interactionVolume.getFaceArea(6, 0), eIdxGlobal7, interactionVolume.getIndexOnElement(6, 0));
2083 addRealFluxFaceArea_(interactionVolume.getFaceArea(6, 1), eIdxGlobal7, interactionVolume.getIndexOnElement(6, 1));
2084 addRealFluxFaceArea_(interactionVolume.getFaceArea(6, 2), eIdxGlobal7, interactionVolume.getIndexOnElement(6, 2));
2085 addRealFluxFaceArea_(interactionVolume.getFaceArea(7, 0), eIdxGlobal8, interactionVolume.getIndexOnElement(7, 0));
2086 addRealFluxFaceArea_(interactionVolume.getFaceArea(7, 1), eIdxGlobal8, interactionVolume.getIndexOnElement(7, 1));
2087 addRealFluxFaceArea_(interactionVolume.getFaceArea(7, 2), eIdxGlobal8, interactionVolume.getIndexOnElement(7, 2));
2088 }
2089
2090 }
2091}
2092} // end namespace Dumux
2093#endif
linteractionvolume3d.hh
Class including the information of an interaction volume of a MPFA 3D method that does not change wit...
Dumux::crossProduct
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
Dumux
Definition: adapt.hh:29
Dumux::GetProp
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type GetProp
get the type of a property (equivalent to old macro GET_PROP(...))
Definition: propertysystem.hh:140
Dumux::GetPropType
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
Dumux::FvMpfaL3dInteractionVolumeContainer
Interactionvolume container for 3-d MPFA L-method.
Definition: 3dinteractionvolumecontainer.hh:46
Dumux::FvMpfaL3dInteractionVolumeContainer::storeIntersectionInfo
void storeIntersectionInfo(const Element &element, std::vector< std::vector< int > > &elemVertMap)
Stores information with respect to DUNE intersections in the interaction volumes.
Definition: 3dinteractionvolumecontainer.hh:371
Dumux::FvMpfaL3dInteractionVolumeContainer::storeBoundaryInteractionVolume
void storeBoundaryInteractionVolume(InteractionVolume &interactionVolume, const Vertex &vertex)
Stores additional information for interaction volumes of boundary vertices.
Definition: 3dinteractionvolumecontainer.hh:1478
Dumux::FvMpfaL3dInteractionVolumeContainer::interactionVolumes_
GlobalInteractionVolumeVector interactionVolumes_
Definition: 3dinteractionvolumecontainer.hh:296
Dumux::FvMpfaL3dInteractionVolumeContainer::FvMpfaL3dInteractionVolumeContainer
FvMpfaL3dInteractionVolumeContainer(Problem &problem)
Constructs a FvMpfaL3dInteractionVolumeContainer object.
Definition: 3dinteractionvolumecontainer.hh:275
Dumux::FvMpfaL3dInteractionVolumeContainer::realFluxFaceArea_
FaceAreaVector realFluxFaceArea_
Definition: 3dinteractionvolumecontainer.hh:297
Dumux::FvMpfaL3dInteractionVolumeContainer::getRealFaceArea
Scalar getRealFaceArea(InteractionVolume &interactionVolume, int elemGlobalIdx, int elemLocalIdx, int localFaceIdx)
Returns the face area of the element.
Definition: 3dinteractionvolumecontainer.hh:248
Dumux::FvMpfaL3dInteractionVolumeContainer::InteractionVolume
GetPropType< TypeTag, Properties::MPFAInteractionVolume > InteractionVolume
Type for storing an MPFA-interaction-volume. (Usually of type FvMpfaL3dInteractionVolume or FvMpfaL3d...
Definition: 3dinteractionvolumecontainer.hh:94
Dumux::FvMpfaL3dInteractionVolumeContainer::interactionVolumesGlobal
GlobalInteractionVolumeVector & interactionVolumesGlobal() const
Returns the interaction volumes container.
Definition: 3dinteractionvolumecontainer.hh:165
Dumux::FvMpfaL3dInteractionVolumeContainer::faceAreaFactor
Scalar faceAreaFactor(InteractionVolume &interactionVolume, int elemGlobalIdx, int elemLocalIdx, int localFaceIdx)
Returns the area weighting factor for the fluxes.
Definition: 3dinteractionvolumecontainer.hh:181
Dumux::FvMpfaL3dInteractionVolumeContainer::problem_
Problem & problem_
Definition: 3dinteractionvolumecontainer.hh:294
Dumux::FvMpfaL3dInteractionVolumeContainer::faceAreaFactor
Scalar faceAreaFactor(int elemGlobalIdx, int indexInInside)
Returns the area weighting factor for the fluxes.
Definition: 3dinteractionvolumecontainer.hh:198
Dumux::FvMpfaL3dInteractionVolumeContainer::interactionVolumesGlobal
GlobalInteractionVolumeVector & interactionVolumesGlobal()
Returns the interaction volumes container.
Definition: 3dinteractionvolumecontainer.hh:159
Dumux::FvMpfaL3dInteractionVolumeContainer::initialize
void initialize(bool solveTwice=true)
Initializes the interaction volume container.
Definition: 3dinteractionvolumecontainer.hh:131
Dumux::FvMpfaL3dInteractionVolumeContainer::addRealFluxFaceArea_
void addRealFluxFaceArea_(Scalar faceArea, int eIdxGlobal, int fIdx)
Definition: 3dinteractionvolumecontainer.hh:285
Dumux::FvMpfaL3dInteractionVolumeContainer::interactionVolume
InteractionVolume & interactionVolume(int vertexIdx) const
Returns an interaction volume.
Definition: 3dinteractionvolumecontainer.hh:153
Dumux::FvMpfaL3dInteractionVolumeContainer::update
void update()
Updates the interaction volume container.
Definition: 3dinteractionvolumecontainer.hh:113
Dumux::FvMpfaL3dInteractionVolumeContainer::storeSubVolumeElements
void storeSubVolumeElements(const Element &element, std::vector< std::vector< int > > &elemVertMap)
Function for storing the elements of an interaction volume and constructing a map from a vertex to it...
Definition: 3dinteractionvolumecontainer.hh:320
Dumux::FvMpfaL3dInteractionVolumeContainer::addRealFaceArea_
void addRealFaceArea_(Scalar faceArea, int eIdxGlobal, int fIdx)
Definition: 3dinteractionvolumecontainer.hh:289
Dumux::FvMpfaL3dInteractionVolumeContainer::getRealFluxFaceArea
Scalar getRealFluxFaceArea(InteractionVolume &interactionVolume, int elemGlobalIdx, int elemLocalIdx, int localFaceIdx)
Returns the area trough which fluxes are calculated by the MPFA.
Definition: 3dinteractionvolumecontainer.hh:216
Dumux::FvMpfaL3dInteractionVolumeContainer::getRealFaceArea
Scalar getRealFaceArea(int elemGlobalIdx, int indexInInside)
Returns the face area of the element.
Definition: 3dinteractionvolumecontainer.hh:263
Dumux::FvMpfaL3dInteractionVolumeContainer::interactionVolume
InteractionVolume & interactionVolume(int vertexIdx)
Returns an interaction volume.
Definition: 3dinteractionvolumecontainer.hh:143
Dumux::FvMpfaL3dInteractionVolumeContainer::getRealFluxFaceArea
Scalar getRealFluxFaceArea(int elemGlobalIdx, int indexInInside)
Returns the area trough which fluxes are calculated by the MPFA.
Definition: 3dinteractionvolumecontainer.hh:231
Dumux::FvMpfaL3dInteractionVolumeContainer::storeInnerInteractionVolume
void storeInnerInteractionVolume(InteractionVolume &interactionVolume, const Vertex &vertex, bool sameLevel=true)
Stores additional information which can be constructed for interaction volumes of non-boundary vertic...
Definition: 3dinteractionvolumecontainer.hh:1341
properties.hh
Properties for a MPFA method.
pressureproperties.hh
Base file for properties related to sequential IMPET algorithms.
Generated by doxygen 1.9.3