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