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