3.1-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
3dtransmissibilitycalculator.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2/*****************************************************************************
3 * See the file COPYING for full copying permissions. *
4 * *
5 * This program is free software: you can redistribute it and/or modify *
6 * it under the terms of the GNU General Public License as published by *
7 * the Free Software Foundation, either version 3 of the License, or *
8 * (at your option) any later version. *
9 * *
10 * This program is distributed in the hope that it will be useful, *
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13 * GNU General Public License for more details. *
14 * *
15 * You should have received a copy of the GNU General Public License *
16 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
17 *****************************************************************************/
23#ifndef DUMUX_FVMPFAL3D_TRANSMISSIBILITYCALCULATOR_HH
24#define DUMUX_FVMPFAL3D_TRANSMISSIBILITYCALCULATOR_HH
25
26// dumux environment
27#include <dumux/porousmediumflow/sequential/pressureproperties.hh>
28#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
29
30namespace Dumux {
31
49template<class TypeTag>
50class FvMpfaL3dTransmissibilityCalculator
51{
52 using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
53
54 enum
55 {
56 dim = GridView::dimension, dimWorld = GridView::dimensionworld
57 };
58
59 enum TransCriterion
60 {
61 sDiff = 0, sSum = 1
62 };
63
64 using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
65 using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
66
67 using Element = typename GridView::template Codim<0>::Entity;
68
69 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
70 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
71
72 using DimVector = Dune::FieldVector<Scalar, dim>;
73
74 using InteractionVolume = typename GET_PROP_TYPE(TypeTag, MPFAInteractionVolume);
75
76
77public:
78 using TransmissibilityType = Dune::FieldMatrix<Scalar, dim, 2*dim - dim + 1>;
79
80 int chooseTransmissibility(TransmissibilityType& transmissibilityOne,
81 TransmissibilityType& transmissibilityTwo, int lTypeOne, int lTypeTwo);
82
83 int transmissibility(Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
84 InteractionVolume& interactionVolume,
85 std::vector<DimVector >& lambda,
86 int idx1, int idx2, int idx3, int idx4, int idx5, int idx6);
87
88 int transmissibility(Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
89 InteractionVolume& interactionVolume,
90 std::vector<DimVector >& lambda,
91 int idx1, int idx2, int idx3, int idx4, int idx5, int idx6,
92 Dune::FieldVector<bool, 4> &useCases);
93
94 int transmissibilityTPFA(Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
95 InteractionVolume& interactionVolume,
96 std::vector<DimVector >& lambda,
97 int idx1, int idx2);
98
99 int transmissibilityCaseOne(
100 Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
101 InteractionVolume& interactionVolume,
102 std::vector<DimVector >& lambda,
103 int idx1, int idx2, int idx3, int idx5);
104
105 int transmissibilityCaseTwo(
106 Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
107 InteractionVolume& interactionVolume,
108 std::vector<DimVector >& lambda,
109 int idx1, int idx2, int idx4, int idx6);
110 int transmissibilityCaseThree(
111 Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
112 InteractionVolume& interactionVolume,
113 std::vector<DimVector >& lambda,
114 int idx1, int idx2, int idx4, int idx5);
115
116 int transmissibilityCaseFour(
117 Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
118 InteractionVolume& interactionVolume,
119 std::vector<DimVector >& lambda,
120 int idx1, int idx2, int idx3, int idx6);
121
127 FvMpfaL3dTransmissibilityCalculator(Problem& problem) :
128 problem_(problem)
129 {
130 if (dim != 3)
131 {
132 DUNE_THROW(Dune::NotImplemented, "Dimension not supported!");
133 }
134
135 enableSimpleLStencil_ = getParam<bool>("MPFA.EnableSimpleLStencil", true);
136 enableComplexLStencil_ = getParam<bool>("MPFA.EnableComplexLStencil", true);
137 enableTPFA_= getParam<bool>("MPFA.EnableTPFA", false);
138
139 if (!enableSimpleLStencil_ && !enableComplexLStencil_)
140 {
141 std::cout<<"MPFA disabled: Use TPFA!\n";
142 enableTPFA_ = true;
143 }
144
145 transChoiceThreshold_ = getParam<Scalar>("MPFA.TransmissibilityCriterionThreshold", 1e-8);
146
147 transCriterion_ = getParam<int>("MPFA.TransmissibilityCriterion", 0);
148 if (transCriterion_ == sDiff)
149 {
150 std::cout << "MPFAL 3D: Using standard transmissibility criterion\n";
151 }
152 else if (transCriterion_ == sSum)
153 {
154 std::cout << "MPFAL 3D: Using accumulative transmissibility criterion\n";
155 }
156 else
157 {
158 transCriterion_ = sDiff;
159 std::cout << "MPFAL 3D: Defined transmissiblity criterion not implemented! Using standard criterion!\n";
160 }
161 }
162
163private:
164 Problem& problem_;
165 bool enableSimpleLStencil_;
166 bool enableComplexLStencil_;
167 bool enableTPFA_;
168 Scalar transChoiceThreshold_;
169 int transCriterion_;
170};
171
185template<class TypeTag>
186int FvMpfaL3dTransmissibilityCalculator<TypeTag>::chooseTransmissibility(TransmissibilityType& transmissibilityOne,
187 TransmissibilityType& transmissibilityTwo,
188 int lTypeOne, int lTypeTwo)
189{
190 if (transCriterion_ == sDiff)
191 {
192 using std::abs;
193 Scalar sOne = abs(transmissibilityOne[0][0] - transmissibilityOne[0][1]);
194 Scalar sTwo = abs(transmissibilityTwo[0][0] - transmissibilityTwo[0][1]);
195
196 //Decide whether to take case1 or case2
197 if (sOne < sTwo - transChoiceThreshold_)
198 {
199 return lTypeOne;
200 }
201 else
202 {
203 return lTypeTwo;
204 }
205 }
206 else if (transCriterion_ == sSum)
207 {
208
209 Scalar tSumOne = 0;
210 Scalar tSumTwo = 0;
211 using std::abs;
212
213 if (lTypeOne == 1)
214 tSumOne = abs(transmissibilityOne[0][0] + transmissibilityOne[0][2] + transmissibilityOne[0][3]);
215 else if (lTypeOne == 2)
216 tSumOne = abs(transmissibilityOne[0][1] + transmissibilityOne[0][2] + transmissibilityOne[0][3]);
217 else if (lTypeOne == 3)
218 tSumOne = abs(transmissibilityOne[0][0] + transmissibilityOne[0][3]);
219 else if (lTypeOne == 4)
220 tSumOne = abs(transmissibilityOne[0][0] + transmissibilityOne[0][2]);
221 else
222 DUNE_THROW(Dune::NotImplemented,"Transmissibility type not implemented");
223
224 if (lTypeTwo == 1)
225 tSumTwo = abs(transmissibilityTwo[0][0] + transmissibilityTwo[0][2] + transmissibilityTwo[0][3]);
226 else if (lTypeTwo == 2)
227 tSumTwo = abs(transmissibilityTwo[0][1] + transmissibilityTwo[0][2] + transmissibilityTwo[0][3]);
228 else if (lTypeTwo == 3)
229 tSumTwo = abs(transmissibilityTwo[0][0] + transmissibilityTwo[0][3]);
230 else if (lTypeTwo == 4)
231 tSumTwo = abs(transmissibilityTwo[0][0] + transmissibilityTwo[0][2]);
232 else
233 DUNE_THROW(Dune::NotImplemented,"Transmissibility type not implemented");
234
235
236 if (tSumOne > tSumTwo + transChoiceThreshold_)
237 {
238 return lTypeOne;
239 }
240 else
241 {
242 return lTypeTwo;
243 }
244 }
245 else
246 {
247 DUNE_THROW(Dune::NotImplemented,"Transmissibility criterion not implemented");
248 }
249
250}
251
267template<class TypeTag>
268int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibility(TransmissibilityType& transmissibility,
269 InteractionVolume& interactionVolume,
270 std::vector<DimVector>& lambda,
271 int idx1, int idx2, int idx3, int idx4,
272 int idx5, int idx6)
273{
274 int level1 = interactionVolume.getSubVolumeElement(idx1).level();
275 int level2 = interactionVolume.getSubVolumeElement(idx2).level();
276
277 if (enableTPFA_ && level1 == level2)
278 {
279 return transmissibilityTPFA(transmissibility, interactionVolume, lambda, idx1, idx2);
280 }
281 else if (enableTPFA_ && !enableSimpleLStencil_ && !enableComplexLStencil_)
282 {
283 return transmissibilityTPFA(transmissibility, interactionVolume, lambda, idx1, idx2);
284 }
285
286 if (enableSimpleLStencil_ && enableComplexLStencil_)
287 {
288 TransmissibilityType transTemp(0);
289
290 // decide which L-stencil (which transmissibility coefficients) to use
291 // calculate transmissibility differences
292 int transTypeTemp = transmissibilityCaseOne(transTemp, interactionVolume, lambda, idx1, idx2, idx3, idx5);
293 int transType = transmissibilityCaseTwo(transmissibility, interactionVolume, lambda, idx1, idx2, idx4, idx6);
294
295 transType = chooseTransmissibility(transTemp, transmissibility, 1, 2);
296
297 if (transType == 1)
298 transmissibility = transTemp;
299
300 TransmissibilityType transCaseThree(0);
301
302 int DUNE_UNUSED transTypeCaseThree = transmissibilityCaseThree(transCaseThree, interactionVolume,
303 lambda, idx1, idx2, idx4, idx5);
304 transTypeTemp = transmissibilityCaseFour(transTemp, interactionVolume, lambda, idx1, idx2, idx3, idx6);
305
306 transTypeTemp = chooseTransmissibility(transCaseThree, transTemp, 3, 4);
307
308 if (transTypeTemp == 3)
309 transTemp = transCaseThree;
310
311 transType = chooseTransmissibility(transmissibility, transTemp, transType, transTypeTemp);
312
313 if (transType == transTypeTemp)
314 transmissibility = transTemp;
315
316 return transType;
317 }
318 else if (enableSimpleLStencil_)
319 {
320 TransmissibilityType transTemp(0);
321
322 int DUNE_UNUSED transTypeTemp = transmissibilityCaseOne(transTemp, interactionVolume,
323 lambda, idx1, idx2, idx3, idx5);
324 int transType = transmissibilityCaseTwo(transmissibility, interactionVolume, lambda, idx1, idx2, idx4, idx6);
325
326 transType = chooseTransmissibility(transTemp, transmissibility, 1, 2);
327
328 if (transType == 1)
329 transmissibility = transTemp;
330
331 return transType;
332 }
333 else if (enableComplexLStencil_)
334 {
335 TransmissibilityType transTemp(0);
336
337 int DUNE_UNUSED transTypeTemp = transmissibilityCaseThree(transTemp, interactionVolume, lambda, idx1, idx2, idx4, idx5);
338 int transType = transmissibilityCaseFour(transmissibility, interactionVolume, lambda, idx1, idx2, idx3, idx6);
339
340 transType = chooseTransmissibility(transTemp, transmissibility, 3, 4);
341
342 if (transType == 3)
343 transmissibility = transTemp;
344
345 return transType;
346 }
347 else
348 {
349 DUNE_THROW(Dune::NotImplemented,"Stencil type not implemented");
350 }
351}
352
369template<class TypeTag>
370int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibility(TransmissibilityType& transmissibility,
371 InteractionVolume& interactionVolume,
372 std::vector<DimVector>& lambda, int idx1,
373 int idx2, int idx3, int idx4, int idx5,
374 int idx6, Dune::FieldVector<bool, 4> &useCases)
375{
376 int level1 = interactionVolume.getSubVolumeElement(idx1).level();
377 int level2 = interactionVolume.getSubVolumeElement(idx2).level();
378
379 if (enableTPFA_ && level1 == level2)
380 {
381 return transmissibilityTPFA(transmissibility, interactionVolume, lambda, idx1, idx2);
382 }
383 else if (enableTPFA_ && !enableSimpleLStencil_ && !enableComplexLStencil_)
384 {
385 return transmissibilityTPFA(transmissibility, interactionVolume, lambda, idx1, idx2);
386 }
387 else if (enableTPFA_ && !enableSimpleLStencil_ && !useCases[2] && !useCases[3])
388 {
389 return transmissibilityTPFA(transmissibility, interactionVolume, lambda, idx1, idx2);
390 }
391 else if (enableTPFA_ && !enableComplexLStencil_ && !useCases[0] && !useCases[1])
392 {
393 return transmissibilityTPFA(transmissibility, interactionVolume, lambda, idx1, idx2);
394 }
395
396 if (enableSimpleLStencil_ && enableComplexLStencil_)
397 {
398 TransmissibilityType transTemp(0);
399 int transTypeTemp = 0;
400 int transType = 0;
401
402 // decide which L-stencil (which transmissibility coefficients) to use
403 // calculate transmissibility differences
404
405 if (useCases[0])
406 transTypeTemp = transmissibilityCaseOne(transTemp, interactionVolume, lambda, idx1, idx2, idx3, idx5);
407
408 if (useCases[1])
409 transType = transmissibilityCaseTwo(transmissibility, interactionVolume, lambda, idx1, idx2, idx4, idx6);
410
411 if (useCases[0] && useCases[1])
412 {
413 transType = chooseTransmissibility(transTemp, transmissibility, 1, 2);
414
415 if (transType == 1)
416 transmissibility = transTemp;
417 }
418 else if (useCases[0])
419 {
420 transmissibility = transTemp;
421 transType = transTypeTemp;
422 }
423
424 TransmissibilityType transCaseThree(0);
425 int transTypeCaseThree = 0;
426
427 if (useCases[2])
428 transTypeCaseThree = transmissibilityCaseThree(transCaseThree, interactionVolume, lambda, idx1, idx2, idx4, idx5);
429
430 if (useCases[3])
431 transTypeTemp = transmissibilityCaseFour(transTemp, interactionVolume, lambda, idx1, idx2, idx3, idx6);
432
433 if (useCases[2] && useCases[3])
434 {
435 transTypeTemp = chooseTransmissibility(transCaseThree, transTemp, 3, 4);
436
437 if (transTypeTemp == 3)
438 transTemp = transCaseThree;
439 }
440 else if (useCases[2])
441 {
442 transTemp = transCaseThree;
443 transTypeTemp = transTypeCaseThree;
444 }
445
446 if ((useCases[0] || useCases[1]) && (useCases[2] || useCases[3]))
447 {
448 transType = chooseTransmissibility(transmissibility, transTemp, transType, transTypeTemp);
449
450 if (transType == transTypeTemp)
451 transmissibility = transTemp;
452 }
453 else if (useCases[2] || useCases[3])
454 {
455 transmissibility = transTemp;
456 transType = transTypeTemp;
457 }
458
459 return transType;
460 }
461 else if (enableSimpleLStencil_)
462 {
463 if (useCases[0] && useCases[1])
464 {
465 TransmissibilityType transTemp(0);
466
467 int DUNE_UNUSED transTypeTemp = transmissibilityCaseOne(transTemp, interactionVolume,
468 lambda, idx1, idx2, idx3, idx5);
469 int transType = transmissibilityCaseTwo(transmissibility, interactionVolume,
470 lambda, idx1, idx2, idx4, idx6);
471
472 transType = chooseTransmissibility(transTemp, transmissibility, 1, 2);
473
474 if (transType == 1)
475 transmissibility = transTemp;
476
477 return transType;
478 }
479 else if (useCases[0])
480 {
481 return transmissibilityCaseOne(transmissibility, interactionVolume, lambda, idx1, idx2, idx3, idx5);
482 }
483 else if (useCases[1])
484 {
485 return transmissibilityCaseTwo(transmissibility, interactionVolume, lambda, idx1, idx2, idx4, idx6);
486 }
487 else
488 {
489 DUNE_THROW(Dune::RangeError,"Only simple L-stencils allowed but no simple stencil chosen!");
490 }
491 }
492 else if (enableComplexLStencil_)
493 {
494 if (useCases[2] && useCases[3])
495 {
496 TransmissibilityType transTemp(0);
497
498 int DUNE_UNUSED transTypeTemp = transmissibilityCaseThree(transTemp, interactionVolume,
499 lambda, idx1, idx2, idx4, idx5);
500 int transType = transmissibilityCaseFour(transmissibility, interactionVolume,
501 lambda, idx1, idx2, idx3, idx6);
502
503 transType = chooseTransmissibility(transTemp, transmissibility, 3, 4);
504
505 if (transType == 3)
506 transmissibility = transTemp;
507
508 return transType;
509 }
510 else if (useCases[2])
511 {
512 return transmissibilityCaseThree(transmissibility, interactionVolume, lambda, idx1, idx2, idx4, idx5);
513 }
514 else if (useCases[3])
515 {
516 return transmissibilityCaseFour(transmissibility, interactionVolume, lambda, idx1, idx2, idx3, idx6);
517 }
518 else
519 {
520 DUNE_THROW(Dune::NotImplemented,"Transmissibility type not implemented");
521 }
522 }
523 else
524 {
525 DUNE_THROW(Dune::NotImplemented,"Stencil type not implemented");
526 }
527}
528
539template<class TypeTag>
540int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibilityTPFA(
541 Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
542 InteractionVolume& interactionVolume,
543 std::vector<DimVector >& lambda, int idx1, int idx2)
544{
545 auto element1 = interactionVolume.getSubVolumeElement(idx1);
546 auto element2 = interactionVolume.getSubVolumeElement(idx2);
547
548 // get global coordinate of cell centers
549 const GlobalPosition& globalPos1 = element1.geometry().center();
550 const GlobalPosition& globalPos2 = element2.geometry().center();
551
552 GlobalPosition distVec = globalPos1 - globalPos2;
553
554 Scalar dist = distVec.two_norm();
555
556 DimVector outerNormal(0);
557 Scalar lambda12 =0;
558 Scalar lambda21 =0;
559 Scalar faceArea = 0;
560
561 if ((idx1 == 0 && idx2 == 1) || (idx1 == 1 && idx2 == 3) || (idx1 == 3 && idx2 == 2) || (idx1 == 2 && idx2 == 0))
562 {
563
564 outerNormal = interactionVolume.getNormal(idx1, 0);
565
566 lambda12 = lambda[idx1][0];
567 lambda21 = lambda[idx2][1];
568
569 faceArea = interactionVolume.getFaceArea(idx1,0);
570 }
571 else if ((idx1 == 5 && idx2 == 4) || (idx1 == 4 && idx2 == 6) || (idx1 == 6 && idx2 == 7) || (idx1 == 7 && idx2 == 5))
572 {
573 outerNormal = interactionVolume.getNormal(idx1, 2);
574
575 lambda12 = lambda[idx1][2];
576 lambda21 = lambda[idx2][1];
577
578 faceArea = interactionVolume.getFaceArea(idx1,2);
579 }
580 else if ((idx1 == 4 && idx2 == 0) || (idx1 == 7 && idx2 == 3))
581 {
582 outerNormal = interactionVolume.getNormal(idx1, 0);
583
584 lambda12 = lambda[idx1][0];
585 lambda21 = lambda[idx2][2];
586
587 faceArea = interactionVolume.getFaceArea(idx1,0);
588 }
589 else
590 {
591 outerNormal = interactionVolume.getNormal(idx1, 2);
592
593 lambda12 = lambda[idx1][2];
594 lambda21 = lambda[idx2][0];
595
596 faceArea = interactionVolume.getFaceArea(idx1,2);
597 }
598
599 const DimMatrix& K1 = problem_.spatialParams().intrinsicPermeability(element1);
600 const DimMatrix& K2 = problem_.spatialParams().intrinsicPermeability(element2);
601
602 // compute vectorized permeabilities
603 DimVector permeability(0);
604 K1.mv(outerNormal, permeability);
605
606 Scalar perm1 = permeability * outerNormal;
607
608 permeability = 0;
609 K2.mv(outerNormal, permeability);
610
611 Scalar perm2 = permeability * outerNormal;
612
613 Scalar meanMob = 0;
614 if ((lambda12 + lambda21) > 0.0)
615 meanMob = 2*lambda12 * perm1 * lambda21 * perm2/(lambda12 * perm1 + lambda21 * perm2);
616
617 Scalar entry = meanMob / dist * faceArea;
618
619 transmissibility = 0;
620 transmissibility[0][0] = entry;
621 transmissibility[0][1] = -entry;
622
623 return 1;
624}
625
644template<class TypeTag>
645int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibilityCaseOne(
646 Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
647 InteractionVolume& interactionVolume,
648 std::vector<DimVector >& lambda, int idx1, int idx2, int idx3, int idx5)
649{
650 auto element1 = interactionVolume.getSubVolumeElement(idx1);
651 auto element2 = interactionVolume.getSubVolumeElement(idx2);
652 auto element3 = interactionVolume.getSubVolumeElement(idx3);
653 auto element5 = interactionVolume.getSubVolumeElement(idx5);
654
655 // get global coordinate of cell centers
656 const GlobalPosition& globalPos1 = element1.geometry().center();
657 const GlobalPosition& globalPos2 = element2.geometry().center();
658 const GlobalPosition& globalPos3 = element3.geometry().center();
659 const GlobalPosition& globalPos5 = element5.geometry().center();
660
661 const GlobalPosition& globalPosCenter = interactionVolume.getCenterPosition();
662
663 const DimMatrix& K1 = problem_.spatialParams().intrinsicPermeability(element1);
664 const DimMatrix& K2 = problem_.spatialParams().intrinsicPermeability(element2);
665 const DimMatrix& K3 = problem_.spatialParams().intrinsicPermeability(element3);
666 const DimMatrix& K5 = problem_.spatialParams().intrinsicPermeability(element5);
667
668 // 1.use centered L-stencil 1 to compute the transmissibility of 1/4 face
669 GlobalPosition edgeCoord1213(0);
670 GlobalPosition edgeCoord1215(0);
671 GlobalPosition edgeCoord1315(0);
672
673 GlobalPosition globalPosFace12(0);
674 GlobalPosition globalPosFace13(0);
675 GlobalPosition globalPosFace15(0);
676
677 DimVector outerNormaln1(0);
678 DimVector outerNormaln2(0);
679 DimVector outerNormaln3(0);
680
681 Scalar lambda12 = 0;
682 Scalar lambda21 = 0;
683 Scalar lambda13 = 0;
684 Scalar lambda31 = 0;
685 Scalar lambda15 = 0;
686 Scalar lambda51 = 0;
687 Scalar faceArea12 = 0;
688 Scalar faceArea21 = 0;
689 Scalar faceArea13 = 0;
690 Scalar faceArea31 = 0;
691 Scalar faceArea15 = 0;
692 Scalar faceArea51 = 0;
693
694 if ((idx1 == 0 && idx2 == 1) || (idx1 == 1 && idx2 == 3) || (idx1 == 3 && idx2 == 2) || (idx1 == 2 && idx2 == 0))
695 {
696 globalPosFace12 = interactionVolume.getFacePosition(idx1,0);
697 globalPosFace13 = interactionVolume.getFacePosition(idx1,1);
698 globalPosFace15 = interactionVolume.getFacePosition(idx1,2);
699
700 edgeCoord1213 = interactionVolume.getEdgePosition(idx1, 0, 1);
701 edgeCoord1215 = interactionVolume.getEdgePosition(idx1, 0, 0);
702 edgeCoord1315 = interactionVolume.getEdgePosition(idx1, 1, 1);
703
704 outerNormaln1 = interactionVolume.getNormal(idx1, 0);
705 outerNormaln2 = interactionVolume.getNormal(idx1, 1);
706 outerNormaln3 = interactionVolume.getNormal(idx1, 2);
707
708 lambda12 = lambda[idx1][0];
709 lambda21 = lambda[idx2][1];
710 lambda13 = lambda[idx1][1];
711 lambda31 = lambda[idx3][0];
712 lambda15 = lambda[idx1][2];
713 lambda51 = lambda[idx5][0];
714
715 faceArea12 = interactionVolume.getFaceArea(idx1,0);
716 faceArea21 = interactionVolume.getFaceArea(idx2,1);
717 faceArea13 = interactionVolume.getFaceArea(idx1,1);
718 faceArea31 = interactionVolume.getFaceArea(idx3,0);
719 faceArea15 = interactionVolume.getFaceArea(idx1,2);
720 faceArea51 = interactionVolume.getFaceArea(idx5,0);
721 }
722 else if ((idx1 == 5 && idx2 == 4) || (idx1 == 4 && idx2 == 6) || (idx1 == 6 && idx2 == 7) || (idx1 == 7 && idx2 == 5))
723 {
724 globalPosFace12 = interactionVolume.getFacePosition(idx1,2);
725 globalPosFace13 = interactionVolume.getFacePosition(idx1,1);
726 globalPosFace15 = interactionVolume.getFacePosition(idx1,0);
727
728 edgeCoord1213 = interactionVolume.getEdgePosition(idx1, 2, 1);
729 edgeCoord1215 = interactionVolume.getEdgePosition(idx1, 2, 0);
730 edgeCoord1315 = interactionVolume.getEdgePosition(idx1, 1, 1);
731
732 outerNormaln1 = interactionVolume.getNormal(idx1, 2);
733 outerNormaln2 = interactionVolume.getNormal(idx1, 1);
734 outerNormaln3 = interactionVolume.getNormal(idx1, 0);
735
736 lambda12 = lambda[idx1][2];
737 lambda21 = lambda[idx2][1];
738 lambda13 = lambda[idx1][1];
739 lambda31 = lambda[idx3][2];
740 lambda15 = lambda[idx1][0];
741 lambda51 = lambda[idx5][2];
742
743 faceArea12 = interactionVolume.getFaceArea(idx1,2);
744 faceArea21 = interactionVolume.getFaceArea(idx2,1);
745 faceArea13 = interactionVolume.getFaceArea(idx1,1);
746 faceArea31 = interactionVolume.getFaceArea(idx3,2);
747 faceArea15 = interactionVolume.getFaceArea(idx1,0);
748 faceArea51 = interactionVolume.getFaceArea(idx5,2);
749 }
750 else if ((idx1 == 4 && idx2 == 0) || (idx1 == 7 && idx2 == 3))
751 {
752 globalPosFace12 = interactionVolume.getFacePosition(idx1,0);
753 globalPosFace13 = interactionVolume.getFacePosition(idx1,2);
754 globalPosFace15 = interactionVolume.getFacePosition(idx1,1);
755
756 edgeCoord1213 = interactionVolume.getEdgePosition(idx1, 0, 1);
757 edgeCoord1215 = interactionVolume.getEdgePosition(idx1, 0, 0);
758 edgeCoord1315 = interactionVolume.getEdgePosition(idx1, 2, 1);
759
760 outerNormaln1 = interactionVolume.getNormal(idx1, 0);
761 outerNormaln2 = interactionVolume.getNormal(idx1, 2);
762 outerNormaln3 = interactionVolume.getNormal(idx1, 1);
763
764 lambda12 = lambda[idx1][0];
765 lambda21 = lambda[idx2][2];
766 lambda13 = lambda[idx1][2];
767 lambda31 = lambda[idx3][1];
768 lambda15 = lambda[idx1][1];
769 lambda51 = lambda[idx5][2];
770
771 faceArea12 = interactionVolume.getFaceArea(idx1,0);
772 faceArea21 = interactionVolume.getFaceArea(idx2,2);
773 faceArea13 = interactionVolume.getFaceArea(idx1,2);
774 faceArea31 = interactionVolume.getFaceArea(idx3,1);
775 faceArea15 = interactionVolume.getFaceArea(idx1,1);
776 faceArea51 = interactionVolume.getFaceArea(idx5,2);
777 }
778 else
779 {
780 globalPosFace12 = interactionVolume.getFacePosition(idx1,2);
781 globalPosFace13 = interactionVolume.getFacePosition(idx1,0);
782 globalPosFace15 = interactionVolume.getFacePosition(idx1,1);
783
784 edgeCoord1213 = interactionVolume.getEdgePosition(idx1, 2, 1);
785 edgeCoord1215 = interactionVolume.getEdgePosition(idx1, 2, 0);
786 edgeCoord1315 = interactionVolume.getEdgePosition(idx1, 0, 1);
787
788 outerNormaln1 = interactionVolume.getNormal(idx1, 2);
789 outerNormaln2 = interactionVolume.getNormal(idx1, 0);
790 outerNormaln3 = interactionVolume.getNormal(idx1, 1);
791
792 lambda12 = lambda[idx1][2];
793 lambda21 = lambda[idx2][0];
794 lambda13 = lambda[idx1][0];
795 lambda31 = lambda[idx3][1];
796 lambda15 = lambda[idx1][1];
797 lambda51 = lambda[idx5][0];
798
799 faceArea12 = interactionVolume.getFaceArea(idx1,2);
800 faceArea21 = interactionVolume.getFaceArea(idx2,0);
801 faceArea13 = interactionVolume.getFaceArea(idx1,0);
802 faceArea31 = interactionVolume.getFaceArea(idx3,1);
803 faceArea15 = interactionVolume.getFaceArea(idx1,1);
804 faceArea51 = interactionVolume.getFaceArea(idx5,0);
805 }
806
807 // compute normal vectors for case 1 (idx1, idx2, idx3, idx5)
808 DimVector crossProductVector1(globalPosFace13-globalPos1);
809 DimVector crossProductVector2(globalPosFace15-globalPos1);
810 DimVector nu11 = crossProduct(crossProductVector1, crossProductVector2);
811 crossProductVector1 = globalPosFace15-globalPos1;
812 crossProductVector2 = globalPosFace12-globalPos1;
813 DimVector nu12 = crossProduct(crossProductVector1, crossProductVector2);
814 crossProductVector1 = globalPosFace12-globalPos1;
815 crossProductVector2 = globalPosFace13-globalPos1;
816 DimVector nu13 = crossProduct(crossProductVector1, crossProductVector2);
817
818 crossProductVector1 = globalPosFace12-globalPos2;
819 crossProductVector2 = edgeCoord1215-globalPos2;
820 DimVector nu21 = crossProduct(crossProductVector1, crossProductVector2);
821 crossProductVector1 = edgeCoord1215-globalPos2;
822 crossProductVector2 = edgeCoord1213-globalPos2;
823 DimVector nu22 = crossProduct(crossProductVector1, crossProductVector2);
824 crossProductVector1 = edgeCoord1213-globalPos2;
825 crossProductVector2 = globalPosFace12-globalPos2;
826 DimVector nu23 = crossProduct(crossProductVector1, crossProductVector2);
827
828 crossProductVector1 = edgeCoord1213-globalPos3;
829 crossProductVector2 = edgeCoord1315-globalPos3;
830 DimVector nu31 = crossProduct(crossProductVector1, crossProductVector2);
831 crossProductVector1 = edgeCoord1315-globalPos3;
832 crossProductVector2 = globalPosFace13-globalPos3;
833 DimVector nu32 = crossProduct(crossProductVector1, crossProductVector2);
834 crossProductVector1 = globalPosFace13-globalPos3;
835 crossProductVector2 = edgeCoord1213-globalPos3;
836 DimVector nu33 = crossProduct(crossProductVector1, crossProductVector2);
837
838 crossProductVector1 = edgeCoord1215-globalPos5;
839 crossProductVector2 = globalPosFace15-globalPos5;
840 DimVector nu51 = crossProduct(crossProductVector1, crossProductVector2);
841 crossProductVector1 = globalPosFace15-globalPos5;
842 crossProductVector2 = edgeCoord1315-globalPos5;
843 DimVector nu52 = crossProduct(crossProductVector1, crossProductVector2);
844 crossProductVector1 = edgeCoord1315-globalPos5;
845 crossProductVector2 = edgeCoord1215-globalPos5;
846 DimVector nu53 = crossProduct(crossProductVector1, crossProductVector2);
847
848 // compute T, i.e., the volume of the subtetrahedron
849 crossProductVector1 = globalPosFace13-globalPos1;
850 crossProductVector2 = globalPosFace15-globalPos1;
851 Scalar T1 = (globalPosFace12-globalPos1) * crossProduct(crossProductVector1, crossProductVector2);
852 crossProductVector1 = edgeCoord1215-globalPos2;
853 crossProductVector2 = edgeCoord1213-globalPos2;
854 Scalar T2 = (globalPosFace12-globalPos2) * crossProduct(crossProductVector1, crossProductVector2);
855 crossProductVector1 = edgeCoord1213-globalPos3;
856 crossProductVector2 = edgeCoord1315-globalPos3;
857 Scalar T3 = (globalPosFace13-globalPos3) * crossProduct(crossProductVector1, crossProductVector2);
858 crossProductVector1 = edgeCoord1315-globalPos5;
859 crossProductVector2 = edgeCoord1215-globalPos5;
860 Scalar T5 = (globalPosFace15-globalPos5) * crossProduct(crossProductVector1, crossProductVector2);
861
862 // compute components needed for flux calculation, denoted as 'omega' and 'r'
863 DimVector K1nu11(0);
864 K1.mv(nu11, K1nu11);
865 DimVector K1nu12(0);
866 K1.mv(nu12, K1nu12);
867 DimVector K1nu13(0);
868 K1.mv(nu13, K1nu13);
869
870 DimVector K2nu21(0);
871 K2.mv(nu21, K2nu21);
872 DimVector K2nu22(0);
873 K2.mv(nu22, K2nu22);
874 DimVector K2nu23(0);
875 K2.mv(nu23, K2nu23);
876
877 DimVector K3nu31(0);
878 K3.mv(nu31, K3nu31);
879 DimVector K3nu32(0);
880 K3.mv(nu32, K3nu32);
881 DimVector K3nu33(0);
882 K3.mv(nu33, K3nu33);
883
884 DimVector K5nu51(0);
885 K5.mv(nu51, K5nu51);
886 DimVector K5nu52(0);
887 K5.mv(nu52, K5nu52);
888 DimVector K5nu53(0);
889 K5.mv(nu53, K5nu53);
890
891 Scalar omega111 = lambda12 * (outerNormaln1 * K1nu11) * faceArea12/T1;
892 Scalar omega112 = lambda12 * (outerNormaln1 * K1nu12) * faceArea12/T1;
893 Scalar omega113 = lambda12 * (outerNormaln1 * K1nu13) * faceArea12/T1;
894
895 Scalar omega121 = lambda21 * (outerNormaln1 * K2nu21) * faceArea21/T2;
896 Scalar omega122 = lambda21 * (outerNormaln1 * K2nu22) * faceArea21/T2;
897 Scalar omega123 = lambda21 * (outerNormaln1 * K2nu23) * faceArea21/T2;
898
899 Scalar omega211 = lambda13 * (outerNormaln2 * K1nu11) * faceArea13/T1;
900 Scalar omega212 = lambda13 * (outerNormaln2 * K1nu12) * faceArea13/T1;
901 Scalar omega213 = lambda13 * (outerNormaln2 * K1nu13) * faceArea13/T1;
902
903 Scalar omega231 = lambda31 * (outerNormaln2 * K3nu31) * faceArea31/T3;
904 Scalar omega232 = lambda31 * (outerNormaln2 * K3nu32) * faceArea31/T3;
905 Scalar omega233 = lambda31 * (outerNormaln2 * K3nu33) * faceArea31/T3;
906
907 Scalar omega311 = lambda15 * (outerNormaln3 * K1nu11) * faceArea15/T1;
908 Scalar omega312 = lambda15 * (outerNormaln3 * K1nu12) * faceArea15/T1;
909 Scalar omega313 = lambda15 * (outerNormaln3 * K1nu13) * faceArea15/T1;
910
911 Scalar omega351 = lambda51 * (outerNormaln3 * K5nu51) * faceArea51/T5;
912 Scalar omega352 = lambda51 * (outerNormaln3 * K5nu52) * faceArea51/T5;
913 Scalar omega353 = lambda51 * (outerNormaln3 * K5nu53) * faceArea51/T5;
914
915 Scalar r111 = (nu11 * (edgeCoord1213-globalPos1))/T1;
916 Scalar r112 = (nu11 * (edgeCoord1215-globalPos1))/T1;
917 Scalar r113 = (nu11 * (edgeCoord1315-globalPos1))/T1;
918
919 Scalar r121 = (nu12 * (edgeCoord1213-globalPos1))/T1;
920 Scalar r122 = (nu12 * (edgeCoord1215-globalPos1))/T1;
921 Scalar r123 = (nu12 * (edgeCoord1315-globalPos1))/T1;
922
923 Scalar r131 = (nu13 * (edgeCoord1213-globalPos1))/T1;
924 Scalar r132 = (nu13 * (edgeCoord1215-globalPos1))/T1;
925 Scalar r133 = (nu13 * (edgeCoord1315-globalPos1))/T1;
926
927 // compute transmissibility matrix T = CA^{-1}B+D
928 DimMatrix C(0), A(0);
929 Dune::FieldMatrix<Scalar,dim,2*dim-dim+1> D(0), B(0);
930
931 // evaluate matrix C, D, A, B
932 C[0][0] = -omega111;
933 C[0][1] = -omega112;
934 C[0][2] = -omega113;
935 C[1][0] = -omega211;
936 C[1][1] = -omega212;
937 C[1][2] = -omega213;
938 C[2][0] = -omega311;
939 C[2][1] = -omega312;
940 C[2][2] = -omega313;
941
942 D[0][0] = omega111 + omega112 + omega113;
943 D[1][0] = omega211 + omega212 + omega213;
944 D[2][0] = omega311 + omega312 + omega313;
945
946 A[0][0] = omega111 - omega122 - omega121*r111 - omega123*r112;
947 A[0][1] = omega112 - omega121*r121 - omega123*r122;
948 A[0][2] = omega113 - omega121*r131 - omega123*r132;
949 A[1][0] = omega211 - omega232*r111 - omega233*r113;
950 A[1][1] = omega212 - omega231 - omega232*r121 - omega233*r123;
951 A[1][2] = omega213 - omega232*r131 - omega233*r133;
952 A[2][0] = omega311 - omega351*r113 - omega352*r112;
953 A[2][1] = omega312 - omega351*r123 - omega352*r122;
954 A[2][2] = omega313 - omega353 - omega351*r133 - omega352*r132;
955
956 B[0][0] = omega111 + omega112 + omega113 + omega121*(1.0 - r111 - r121 -r131) + omega123*(1.0 - r112 - r122 - r132);
957 B[0][1] = -omega121 - omega122 - omega123;
958 B[1][0] = omega211 + omega212 + omega213 + omega232*(1.0 - r111 - r121 - r131) + omega233*(1.0 - r113 - r123 - r133);
959 B[1][2] = -omega231 - omega232 - omega233;
960 B[2][0] = omega311 + omega312 + omega313 + omega351*(1.0 - r113 - r123 - r133) + omega352*(1.0 - r112 -r122 -r132);
961 B[2][3] = -omega351 - omega352 -omega353;
962
963 // compute T
964 A.invert();
965 D += B.leftmultiply(C.rightmultiply(A));
966
967 transmissibility = D;
968
969 using std::isnan;
970 if (isnan(transmissibility.frobenius_norm()))
971 {
972 std::cout<<"idx: "<<idx1<<idx2<<idx3<<idx5<<"\n";
973
974 // interactionVolume.printInteractionVolumeInfo();
975
976 std::cout<<"case 1: transmissibility = "<<transmissibility<<"\n";
977 std::cout<<"globalPos1 = "<<globalPos1<<"\n";
978 std::cout<<"globalPos2 = "<<globalPos2<<"\n";
979 std::cout<<"globalPos3 = "<<globalPos3<<"\n";
980 std::cout<<"globalPos5 = "<<globalPos5<<"\n";
981 std::cout<<"globalPosCenter = "<<globalPosCenter<<"\n";
982 std::cout<<"outerNormaln1 = "<<outerNormaln1<<"\n";
983 std::cout<<"outerNormaln2 = "<<outerNormaln2<<"\n";
984 std::cout<<"outerNormaln3 = "<<outerNormaln3<<"\n";
985 std::cout<<"xbar_1 = "<<globalPosFace12<<"\n";
986 std::cout<<"xbar_2 = "<<globalPosFace13<<"\n";
987 std::cout<<"xbar_3 = "<<globalPosFace15<<"\n";
988 std::cout<<"xbar_4 = "<<edgeCoord1213<<"\n";
989 std::cout<<"xbar_5 = "<<edgeCoord1215<<"\n";
990 std::cout<<"xbar_6 = "<<edgeCoord1315<<"\n";
991 std::cout<<"perm1 = "<<K1<<"\n";
992 std::cout<<"perm2 = "<<K2<<"\n";
993 std::cout<<"perm3 = "<<K3<<"\n";
994 std::cout<<"perm5 = "<<K5<<"\n";
995 std::cout<<"lambda = ";
996 for (unsigned int i = 0; i < lambda.size(); i++)
997 {
998 std::cout<<lambda[i]<<" ";
999 }
1000 std::cout<<"\n";
1001 std::cout<<"\n";
1002 DUNE_THROW(Dune::MathError,"T is nan");
1003 }
1004
1005 return 1;
1006}
1007
1026template<class TypeTag>
1027int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibilityCaseTwo(
1028 Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
1029 InteractionVolume& interactionVolume,
1030 std::vector<DimVector >& lambda, int idx1, int idx2, int idx4, int idx6)
1031{
1032 auto element1 = interactionVolume.getSubVolumeElement(idx1);
1033 auto element2 = interactionVolume.getSubVolumeElement(idx2);
1034 auto element4 = interactionVolume.getSubVolumeElement(idx4);
1035 auto element6 = interactionVolume.getSubVolumeElement(idx6);
1036
1037 // get global coordinate of cell centers
1038 const GlobalPosition& globalPos1 = element1.geometry().center();
1039 const GlobalPosition& globalPos2 = element2.geometry().center();
1040 const GlobalPosition& globalPos4 = element4.geometry().center();
1041 const GlobalPosition& globalPos6 = element6.geometry().center();
1042
1043 const GlobalPosition& globalPosCenter = interactionVolume.getCenterPosition();
1044
1045 const DimMatrix& K1 = problem_.spatialParams().intrinsicPermeability(element1);
1046 const DimMatrix& K2 = problem_.spatialParams().intrinsicPermeability(element2);
1047 const DimMatrix& K4 = problem_.spatialParams().intrinsicPermeability(element4);
1048 const DimMatrix& K6 = problem_.spatialParams().intrinsicPermeability(element6);
1049
1050 // 1.use centered L-stencil 1 to compute the transmissibility of 1/4 face
1051 GlobalPosition globalPosFace12(0);
1052 GlobalPosition globalPosFace24(0);
1053 GlobalPosition globalPosFace26(0);
1054
1055 DimVector outerNormaln1(0);
1056 DimVector outerNormaln2(0);
1057 DimVector outerNormaln3(0);
1058
1059 GlobalPosition edgeCoord1224(0);
1060 GlobalPosition edgeCoord1226(0);
1061 GlobalPosition edgeCoord2426(0);
1062
1063 // get lambda and face area
1064 Scalar lambda12 = 0;
1065 Scalar lambda21 = 0;
1066 Scalar lambda24 = 0;
1067 Scalar lambda42 = 0;
1068 Scalar lambda26 = 0;
1069 Scalar lambda62 = 0;
1070 Scalar faceArea12 = 0;
1071 Scalar faceArea21 = 0;
1072 Scalar faceArea24 = 0;
1073 Scalar faceArea42 = 0;
1074 Scalar faceArea26 = 0;
1075 Scalar faceArea62 = 0;
1076
1077 if ((idx1 == 0 && idx2 == 1) || (idx1 == 1 && idx2 == 3) || (idx1 == 3 && idx2 == 2) || (idx1 == 2 && idx2 == 0))
1078 {
1079 globalPosFace12 = interactionVolume.getFacePosition(idx1,0);
1080
1081 outerNormaln1 = interactionVolume.getNormal(idx1, 0);
1082
1083 globalPosFace24 = interactionVolume.getFacePosition(idx2,0);
1084 globalPosFace26 = interactionVolume.getFacePosition(idx2,2);
1085
1086 edgeCoord1224 = interactionVolume.getEdgePosition(idx2, 1, 0);
1087 edgeCoord1226 = interactionVolume.getEdgePosition(idx2, 1, 1);
1088 edgeCoord2426 = interactionVolume.getEdgePosition(idx2, 0, 0);
1089
1090 outerNormaln2 = interactionVolume.getNormal(idx2, 0);
1091 outerNormaln3 = interactionVolume.getNormal(idx2, 2);
1092
1093 lambda12 = lambda[idx1][0];
1094 lambda21 = lambda[idx2][1];
1095 lambda24 = lambda[idx2][0];
1096 lambda42 = lambda[idx4][1];
1097 lambda26 = lambda[idx2][2];
1098 lambda62 = lambda[idx6][0];
1099
1100 faceArea12 = interactionVolume.getFaceArea(idx1,0);
1101 faceArea21 = interactionVolume.getFaceArea(idx2,1);
1102 faceArea24 = interactionVolume.getFaceArea(idx2,0);
1103 faceArea42 = interactionVolume.getFaceArea(idx4,1);
1104 faceArea26 = interactionVolume.getFaceArea(idx2,2);
1105 faceArea62 = interactionVolume.getFaceArea(idx6,0);
1106 }
1107 else if ((idx1 == 5 && idx2 == 4) || (idx1 == 4 && idx2 == 6) || (idx1 == 6 && idx2 == 7) || (idx1 == 7 && idx2 == 5))
1108 {
1109 globalPosFace12 = interactionVolume.getFacePosition(idx1,2);
1110
1111 outerNormaln1 = interactionVolume.getNormal(idx1, 2);
1112
1113 globalPosFace24 = interactionVolume.getFacePosition(idx2,2);
1114 globalPosFace26 = interactionVolume.getFacePosition(idx2,0);
1115
1116 edgeCoord1224 = interactionVolume.getEdgePosition(idx2, 1, 0);
1117 edgeCoord1226 = interactionVolume.getEdgePosition(idx2, 1, 1);
1118 edgeCoord2426 = interactionVolume.getEdgePosition(idx2, 2, 0);
1119
1120 outerNormaln2 = interactionVolume.getNormal(idx2, 2);
1121 outerNormaln3 = interactionVolume.getNormal(idx2, 0);
1122
1123 lambda12 = lambda[idx1][2];
1124 lambda21 = lambda[idx2][1];
1125 lambda24 = lambda[idx2][2];
1126 lambda42 = lambda[idx4][1];
1127 lambda26 = lambda[idx2][0];
1128 lambda62 = lambda[idx6][2];
1129
1130 faceArea12 = interactionVolume.getFaceArea(idx1,2);
1131 faceArea21 = interactionVolume.getFaceArea(idx2,1);
1132 faceArea24 = interactionVolume.getFaceArea(idx2,2);
1133 faceArea42 = interactionVolume.getFaceArea(idx4,1);
1134 faceArea26 = interactionVolume.getFaceArea(idx2,0);
1135 faceArea62 = interactionVolume.getFaceArea(idx6,2);
1136 }
1137 else if ((idx1 == 4 && idx2 == 0) || (idx1 == 7 && idx2 == 3))
1138 {
1139 globalPosFace12 = interactionVolume.getFacePosition(idx1,0);
1140
1141 outerNormaln1 = interactionVolume.getNormal(idx1, 0);
1142 outerNormaln2 = interactionVolume.getNormal(idx2, 1);
1143 outerNormaln3 = interactionVolume.getNormal(idx2, 0);
1144
1145 globalPosFace24 = interactionVolume.getFacePosition(idx2,1);
1146 globalPosFace26 = interactionVolume.getFacePosition(idx2,0);
1147
1148 edgeCoord1224 = interactionVolume.getEdgePosition(idx2, 2, 0);
1149 edgeCoord1226 = interactionVolume.getEdgePosition(idx2, 2, 1);
1150 edgeCoord2426 = interactionVolume.getEdgePosition(idx2, 1, 0);
1151
1152 lambda12 = lambda[idx1][0];
1153 lambda21 = lambda[idx2][2];
1154 lambda24 = lambda[idx2][1];
1155 lambda42 = lambda[idx4][0];
1156 lambda26 = lambda[idx2][0];
1157 lambda62 = lambda[idx6][1];
1158
1159 faceArea12 = interactionVolume.getFaceArea(idx1,0);
1160 faceArea21 = interactionVolume.getFaceArea(idx2,2);
1161 faceArea24 = interactionVolume.getFaceArea(idx2,1);
1162 faceArea42 = interactionVolume.getFaceArea(idx4,0);
1163 faceArea26 = interactionVolume.getFaceArea(idx2,0);
1164 faceArea62 = interactionVolume.getFaceArea(idx6,1);
1165 }
1166 else
1167 {
1168 globalPosFace12 = interactionVolume.getFacePosition(idx1,2);
1169
1170 outerNormaln1 = interactionVolume.getNormal(idx1, 2);
1171
1172 globalPosFace24 = interactionVolume.getFacePosition(idx2,1);
1173 globalPosFace26 = interactionVolume.getFacePosition(idx2,2);
1174
1175 edgeCoord1224 = interactionVolume.getEdgePosition(idx2, 0, 0);
1176 edgeCoord1226 = interactionVolume.getEdgePosition(idx2, 0, 1);
1177 edgeCoord2426 = interactionVolume.getEdgePosition(idx2, 1, 0);
1178
1179 outerNormaln2 = interactionVolume.getNormal(idx2, 1);
1180 outerNormaln3 = interactionVolume.getNormal(idx2, 2);
1181
1182 lambda12 = lambda[idx1][2];
1183 lambda21 = lambda[idx2][0];
1184 lambda24 = lambda[idx2][1];
1185 lambda42 = lambda[idx4][2];
1186 lambda26 = lambda[idx2][2];
1187 lambda62 = lambda[idx6][1];
1188
1189 faceArea12 = interactionVolume.getFaceArea(idx1,2);
1190 faceArea21 = interactionVolume.getFaceArea(idx2,0);
1191 faceArea24 = interactionVolume.getFaceArea(idx2,1);
1192 faceArea42 = interactionVolume.getFaceArea(idx4,2);
1193 faceArea26 = interactionVolume.getFaceArea(idx2,2);
1194 faceArea62 = interactionVolume.getFaceArea(idx6,1);
1195 }
1196
1197
1198 // compute transmissibility matrix TC1 = CA^{-1}B+D
1199 DimMatrix C(0), A(0);
1200 Dune::FieldMatrix<Scalar,dim,2*dim-dim+1> D(0), B(0);
1201
1202 // compute normal vectors for case 2 (idx1, idx2, idx4, idx6)
1203 DimVector crossProductVector1(edgeCoord1226-globalPos1);
1204 DimVector crossProductVector2(globalPosFace12-globalPos1);
1205 DimVector nu11 = crossProduct(crossProductVector1, crossProductVector2);
1206 crossProductVector1 = edgeCoord1224-globalPos1;
1207 crossProductVector2 = edgeCoord1226-globalPos1;
1208 DimVector nu12 = crossProduct(crossProductVector1, crossProductVector2);
1209 crossProductVector1 = globalPosFace12-globalPos1;
1210 crossProductVector2 = edgeCoord1224-globalPos1;
1211 DimVector nu13 = crossProduct(crossProductVector1, crossProductVector2);
1212
1213 crossProductVector1 = globalPosFace26-globalPos2;
1214 crossProductVector2 = globalPosFace24-globalPos2;
1215 DimVector nu21 = crossProduct(crossProductVector1, crossProductVector2);
1216 crossProductVector1 = globalPosFace12-globalPos2;
1217 crossProductVector2 = globalPosFace26-globalPos2;
1218 DimVector nu22 = crossProduct(crossProductVector1, crossProductVector2);
1219 crossProductVector1 = globalPosFace24-globalPos2;
1220 crossProductVector2 = globalPosFace12-globalPos2;
1221 DimVector nu23 = crossProduct(crossProductVector1, crossProductVector2);
1222
1223 crossProductVector1 = edgeCoord2426-globalPos4;
1224 crossProductVector2 = edgeCoord1224-globalPos4;
1225 DimVector nu41 = crossProduct(crossProductVector1, crossProductVector2);
1226 crossProductVector1 = globalPosFace24-globalPos4;
1227 crossProductVector2 = edgeCoord2426-globalPos4;
1228 DimVector nu42 = crossProduct(crossProductVector1, crossProductVector2);
1229 crossProductVector1 = edgeCoord1224-globalPos4;
1230 crossProductVector2 = globalPosFace24-globalPos4;
1231 DimVector nu43 = crossProduct(crossProductVector1, crossProductVector2);
1232
1233 crossProductVector1 = globalPosFace26-globalPos6;
1234 crossProductVector2 = edgeCoord1226-globalPos6;
1235 DimVector nu61 = crossProduct(crossProductVector1, crossProductVector2);
1236 crossProductVector1 = edgeCoord2426-globalPos6;
1237 crossProductVector2 = globalPosFace26-globalPos6;
1238 DimVector nu62 = crossProduct(crossProductVector1, crossProductVector2);
1239 crossProductVector1 = edgeCoord1226-globalPos6;
1240 crossProductVector2 = edgeCoord2426-globalPos6;
1241 DimVector nu63 = crossProduct(crossProductVector1, crossProductVector2);
1242
1243 // compute T, i.e., the volume of the subtetrahedron
1244 crossProductVector1 = edgeCoord1224-globalPos1;
1245 crossProductVector2 = edgeCoord1226-globalPos1;
1246 Scalar T1 = (globalPosFace12-globalPos1) * crossProduct(crossProductVector1, crossProductVector2);
1247 crossProductVector1 = globalPosFace26-globalPos2;
1248 crossProductVector2 = globalPosFace24-globalPos2;
1249 Scalar T2 = (globalPosFace12-globalPos2) * crossProduct(crossProductVector1, crossProductVector2);
1250 crossProductVector1 = edgeCoord2426-globalPos4;
1251 crossProductVector2 = edgeCoord1224-globalPos4;
1252 Scalar T4 = (globalPosFace24-globalPos4) * crossProduct(crossProductVector1, crossProductVector2);
1253 crossProductVector1 = edgeCoord1226-globalPos6;
1254 crossProductVector2 = edgeCoord2426-globalPos6;
1255 Scalar T6 = (globalPosFace26-globalPos6) * crossProduct(crossProductVector1, crossProductVector2);
1256
1257 // compute components needed for flux calculation, denoted as 'omega' and 'r'
1258 DimVector K1nu11(0);
1259 K1.mv(nu11, K1nu11);
1260 DimVector K1nu12(0);
1261 K1.mv(nu12, K1nu12);
1262 DimVector K1nu13(0);
1263 K1.mv(nu13, K1nu13);
1264
1265 DimVector K2nu21(0);
1266 K2.mv(nu21, K2nu21);
1267 DimVector K2nu22(0);
1268 K2.mv(nu22, K2nu22);
1269 DimVector K2nu23(0);
1270 K2.mv(nu23, K2nu23);
1271
1272 DimVector K4nu41(0);
1273 K4.mv(nu41, K4nu41);
1274 DimVector K4nu42(0);
1275 K4.mv(nu42, K4nu42);
1276 DimVector K4nu43(0);
1277 K4.mv(nu43, K4nu43);
1278
1279 DimVector K6nu61(0);
1280 K6.mv(nu61, K6nu61);
1281 DimVector K6nu62(0);
1282 K6.mv(nu62, K6nu62);
1283 DimVector K6nu63(0);
1284 K6.mv(nu63, K6nu63);
1285
1286 Scalar omega111 = lambda12 * (outerNormaln1 * K1nu11) * faceArea12/T1;
1287 Scalar omega112 = lambda12 * (outerNormaln1 * K1nu12) * faceArea12/T1;
1288 Scalar omega113 = lambda12 * (outerNormaln1 * K1nu13) * faceArea12/T1;
1289
1290 Scalar omega121 = lambda21 * (outerNormaln1 * K2nu21) * faceArea21/T2;
1291 Scalar omega122 = lambda21 * (outerNormaln1 * K2nu22) * faceArea21/T2;
1292 Scalar omega123 = lambda21 * (outerNormaln1 * K2nu23) * faceArea21/T2;
1293
1294 Scalar omega221 = lambda24 * (outerNormaln2 * K2nu21) * faceArea24/T2;
1295 Scalar omega222 = lambda24 * (outerNormaln2 * K2nu22) * faceArea24/T2;
1296 Scalar omega223 = lambda24 * (outerNormaln2 * K2nu23) * faceArea24/T2;
1297
1298 Scalar omega241 = lambda42 * (outerNormaln2 * K4nu41) * faceArea42/T4;
1299 Scalar omega242 = lambda42 * (outerNormaln2 * K4nu42) * faceArea42/T4;
1300 Scalar omega243 = lambda42 * (outerNormaln2 * K4nu43) * faceArea42/T4;
1301
1302 Scalar omega321 = lambda26 * (outerNormaln3 * K2nu21) * faceArea26/T2;
1303 Scalar omega322 = lambda26 * (outerNormaln3 * K2nu22) * faceArea26/T2;
1304 Scalar omega323 = lambda26 * (outerNormaln3 * K2nu23) * faceArea26/T2;
1305
1306 Scalar omega361 = lambda62 * (outerNormaln3 * K6nu61) * faceArea62/T6;
1307 Scalar omega362 = lambda62 * (outerNormaln3 * K6nu62) * faceArea62/T6;
1308 Scalar omega363 = lambda62 * (outerNormaln3 * K6nu63) * faceArea62/T6;
1309
1310 Scalar r211 = (nu21 * (edgeCoord1224-globalPos2))/T2;
1311 Scalar r212 = (nu21 * (edgeCoord1226-globalPos2))/T2;
1312 Scalar r213 = (nu21 * (edgeCoord2426-globalPos2))/T2;
1313
1314 Scalar r221 = (nu22 * (edgeCoord1224-globalPos2))/T2;
1315 Scalar r222 = (nu22 * (edgeCoord1226-globalPos2))/T2;
1316 Scalar r223 = (nu22 * (edgeCoord2426-globalPos2))/T2;
1317
1318 Scalar r231 = (nu23 * (edgeCoord1224-globalPos2))/T2;
1319 Scalar r232 = (nu23 * (edgeCoord1226-globalPos2))/T2;
1320 Scalar r233 = (nu23 * (edgeCoord2426-globalPos2))/T2;
1321
1322 // compute transmissibility matrix T = CA^{-1}B+D
1323 C = 0;
1324 A = 0;
1325 D = 0;
1326 B = 0;
1327
1328 // evaluate matrix C, D, A, B
1329 C[0][0] = -omega121;
1330 C[0][1] = -omega122;
1331 C[0][2] = -omega123;
1332 C[1][0] = -omega221;
1333 C[1][1] = -omega222;
1334 C[1][2] = -omega223;
1335 C[2][0] = -omega321;
1336 C[2][1] = -omega322;
1337 C[2][2] = -omega323;
1338
1339 D[0][1] = omega121 + omega122 + omega123;
1340 D[1][1] = omega221 + omega222 + omega223;
1341 D[2][1] = omega321 + omega322 + omega323;
1342
1343 A[0][0] = omega121 - omega112 - omega111*r211 - omega113*r212;
1344 A[0][1] = omega122 - omega111*r221 - omega113*r222;
1345 A[0][2] = omega123 - omega111*r231 - omega113*r232;
1346 A[1][0] = omega221 - omega242*r211 - omega243*r213;
1347 A[1][1] = omega222 - omega241 - omega242*r221 - omega243*r223;
1348 A[1][2] = omega223 - omega242*r231 - omega243*r233;
1349 A[2][0] = omega321 - omega361*r213 - omega362*r212;
1350 A[2][1] = omega322 - omega361*r223 - omega362*r222;
1351 A[2][2] = omega323 - omega363 - omega361*r233 - omega362*r232;
1352
1353 B[0][0] = -omega111 - omega112 - omega113;
1354 B[0][1] = omega121 + omega122 + omega123 + omega111*(1.0 - r211 - r221 -r231) + omega113*(1.0 - r212 - r222 - r232);
1355 B[1][1] = omega221 + omega222 + omega223 + omega242*(1.0 - r211 - r221 - r231) + omega243*(1.0 - r213 - r223 - r233);
1356 B[1][2] = -omega241 - omega242 - omega243;
1357 B[2][1] = omega321 + omega322 + omega323 + omega361*(1.0 - r213 - r223 - r233) + omega362*(1.0 - r212 -r222 -r232);
1358 B[2][3] = -omega361 - omega362 -omega363;
1359
1360 // compute T
1361 A.invert();
1362 D += B.leftmultiply(C.rightmultiply(A));
1363
1364 transmissibility = D;
1365
1366 using std::isnan;
1367 if (isnan(transmissibility.frobenius_norm()))
1368 {
1369 std::cout<<"idx: "<<idx1<<idx2<<idx4<<idx6<<"\n";
1370
1371 // interactionVolume.printInteractionVolumeInfo();
1372
1373 std::cout<<"case 2: transmissibility = "<<transmissibility<<"\n";
1374 std::cout<<"globalPos1 = "<<globalPos1<<"\n";
1375 std::cout<<"globalPos2 = "<<globalPos2<<"\n";
1376 std::cout<<"globalPos4 = "<<globalPos4<<"\n";
1377 std::cout<<"globalPos6 = "<<globalPos6<<"\n";
1378 std::cout<<"globalPosCenter = "<<globalPosCenter<<"\n";
1379 std::cout<<"outerNormaln1 = "<<outerNormaln1<<"\n";
1380 std::cout<<"outerNormaln2 = "<<outerNormaln2<<"\n";
1381 std::cout<<"outerNormaln3 = "<<outerNormaln3<<"\n";
1382 std::cout<<"xbar_1 = "<<globalPosFace12<<"\n";
1383 std::cout<<"xbar_2 = "<<globalPosFace24<<"\n";
1384 std::cout<<"xbar_3 = "<<globalPosFace26<<"\n";
1385 std::cout<<"xbar_6 = "<<edgeCoord2426<<"\n";
1386 std::cout<<"perm1 = "<<K1<<"\n";
1387 std::cout<<"perm2 = "<<K2<<"\n";
1388 std::cout<<"perm4 = "<<K4<<"\n";
1389 std::cout<<"perm6 = "<<K6<<"\n";
1390 std::cout<<"lambda = ";
1391 for (unsigned int i = 0; i < lambda.size(); i++)
1392 {
1393 std::cout<<lambda[i]<<" ";
1394 }
1395 std::cout<<"\n";
1396 std::cout<<"\n";
1397 DUNE_THROW(Dune::MathError,"T is nan");
1398 }
1399
1400 return 2;
1401}
1402
1421template<class TypeTag>
1422int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibilityCaseThree(
1423 Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
1424 InteractionVolume& interactionVolume,
1425 std::vector<DimVector >& lambda,int idx1, int idx2, int idx4, int idx5)
1426{
1427 auto element1 = interactionVolume.getSubVolumeElement(idx1);
1428 auto element2 = interactionVolume.getSubVolumeElement(idx2);
1429 auto element4 = interactionVolume.getSubVolumeElement(idx4);
1430 auto element5 = interactionVolume.getSubVolumeElement(idx5);
1431
1432 // get global coordinate of cell centers
1433 const GlobalPosition& globalPos1 = element1.geometry().center();
1434 const GlobalPosition& globalPos2 = element2.geometry().center();
1435 const GlobalPosition& globalPos4 = element4.geometry().center();
1436 const GlobalPosition& globalPos5 = element5.geometry().center();
1437
1438 const GlobalPosition& globalPosCenter = interactionVolume.getCenterPosition();
1439
1440 const DimMatrix& K1 = problem_.spatialParams().intrinsicPermeability(element1);
1441 const DimMatrix& K2 = problem_.spatialParams().intrinsicPermeability(element2);
1442 const DimMatrix& K4 = problem_.spatialParams().intrinsicPermeability(element4);
1443 const DimMatrix& K5 = problem_.spatialParams().intrinsicPermeability(element5);
1444
1445 // 1.use centered L-stencil 1 to compute the transmissibility of 1/4 face
1446 GlobalPosition globalPosFace12(0);
1447 GlobalPosition globalPosFace15(0);
1448 GlobalPosition globalPosFace24(0);
1449
1450 DimVector outerNormaln1(0);
1451 DimVector outerNormaln2(0);
1452 DimVector outerNormaln3(0);
1453
1454 GlobalPosition edgeCoord1215(0);
1455 GlobalPosition edgeCoord1315(0);
1456 GlobalPosition edgeCoord1224(0);
1457 GlobalPosition edgeCoord2426(0);
1458
1459 // get lambda and face area
1460 Scalar lambda12 = 0;
1461 Scalar lambda21 = 0;
1462 Scalar lambda24 = 0;
1463 Scalar lambda42 = 0;
1464 Scalar lambda15 = 0;
1465 Scalar lambda51 = 0;
1466 Scalar faceArea12 = 0;
1467 Scalar faceArea24 = 0;
1468 Scalar faceArea15 = 0;
1469 Scalar faceArea21 = 0;
1470 Scalar faceArea42 = 0;
1471 Scalar faceArea51 = 0;
1472
1473 if ((idx1 == 0 && idx2 == 1) || (idx1 == 1 && idx2 == 3) || (idx1 == 3 && idx2 == 2) || (idx1 == 2 && idx2 == 0))
1474 {
1475 globalPosFace12 = interactionVolume.getFacePosition(idx1,0);
1476 globalPosFace15 = interactionVolume.getFacePosition(idx1,2);
1477 globalPosFace24 = interactionVolume.getFacePosition(idx2,0);
1478
1479 outerNormaln1 = interactionVolume.getNormal(idx1, 0);
1480 outerNormaln2 = interactionVolume.getNormal(idx2, 0);
1481 outerNormaln3 = interactionVolume.getNormal(idx1, 2);
1482
1483 edgeCoord1215 = interactionVolume.getEdgePosition(idx1, 0, 0);
1484 edgeCoord1315 = interactionVolume.getEdgePosition(idx1, 1, 1);
1485 edgeCoord1224 = interactionVolume.getEdgePosition(idx2, 1, 0);
1486 edgeCoord2426 = interactionVolume.getEdgePosition(idx2, 0, 0);
1487
1488 lambda12 = lambda[idx1][0];
1489 lambda21 = lambda[idx2][1];
1490 lambda24 = lambda[idx2][0];
1491 lambda42 = lambda[idx4][1];
1492 lambda15 = lambda[idx1][2];
1493 lambda51 = lambda[idx5][0];
1494
1495 faceArea12 = interactionVolume.getFaceArea(idx1,0);
1496 faceArea21 = interactionVolume.getFaceArea(idx2,1);
1497 faceArea24 = interactionVolume.getFaceArea(idx2,0);
1498 faceArea42 = interactionVolume.getFaceArea(idx4,1);
1499 faceArea15 = interactionVolume.getFaceArea(idx1,2);
1500 faceArea51 = interactionVolume.getFaceArea(idx5,0);
1501 }
1502 else if ((idx1 == 5 && idx2 == 4) || (idx1 == 4 && idx2 == 6) || (idx1 == 6 && idx2 == 7) || (idx1 == 7 && idx2 == 5))
1503 {
1504 globalPosFace12 = interactionVolume.getFacePosition(idx1,2);
1505 globalPosFace15 = interactionVolume.getFacePosition(idx1,0);
1506 globalPosFace24 = interactionVolume.getFacePosition(idx2,2);
1507
1508 outerNormaln1 = interactionVolume.getNormal(idx1, 2);
1509 outerNormaln2 = interactionVolume.getNormal(idx2, 2);
1510 outerNormaln3 = interactionVolume.getNormal(idx1, 0);
1511
1512 edgeCoord1215 = interactionVolume.getEdgePosition(idx1, 2, 0);
1513 edgeCoord1315 = interactionVolume.getEdgePosition(idx1, 1, 1);
1514 edgeCoord1224 = interactionVolume.getEdgePosition(idx2, 1, 0);
1515 edgeCoord2426 = interactionVolume.getEdgePosition(idx2, 2, 0);
1516
1517 lambda12 = lambda[idx1][2];
1518 lambda21 = lambda[idx2][1];
1519 lambda24 = lambda[idx2][2];
1520 lambda42 = lambda[idx4][1];
1521 lambda15 = lambda[idx1][0];
1522 lambda51 = lambda[idx5][2];
1523
1524 faceArea12 = interactionVolume.getFaceArea(idx1,2);
1525 faceArea21 = interactionVolume.getFaceArea(idx2,1);
1526 faceArea24 = interactionVolume.getFaceArea(idx2,2);
1527 faceArea42 = interactionVolume.getFaceArea(idx4,1);
1528 faceArea15 = interactionVolume.getFaceArea(idx1,0);
1529 faceArea51 = interactionVolume.getFaceArea(idx5,2);
1530 }
1531 else if ((idx1 == 4 && idx2 == 0) || (idx1 == 7 && idx2 == 3))
1532 {
1533 globalPosFace12 = interactionVolume.getFacePosition(idx1,0);
1534 globalPosFace15 = interactionVolume.getFacePosition(idx1,1);
1535 globalPosFace24 = interactionVolume.getFacePosition(idx2,1);
1536
1537 outerNormaln1 = interactionVolume.getNormal(idx1, 0);
1538 outerNormaln2 = interactionVolume.getNormal(idx2, 1);
1539 outerNormaln3 = interactionVolume.getNormal(idx1, 1);
1540
1541 edgeCoord1215 = interactionVolume.getEdgePosition(idx1, 0, 0);
1542 edgeCoord1315 = interactionVolume.getEdgePosition(idx1, 2, 1);
1543 edgeCoord1224 = interactionVolume.getEdgePosition(idx2, 2, 0);
1544 edgeCoord2426 = interactionVolume.getEdgePosition(idx2, 1, 0);
1545
1546 lambda12 = lambda[idx1][0];
1547 lambda21 = lambda[idx2][2];
1548 lambda24 = lambda[idx2][1];
1549 lambda42 = lambda[idx4][0];
1550 lambda15 = lambda[idx1][1];
1551 lambda51 = lambda[idx5][2];
1552
1553 faceArea12 = interactionVolume.getFaceArea(idx1,0);
1554 faceArea21 = interactionVolume.getFaceArea(idx2,2);
1555 faceArea24 = interactionVolume.getFaceArea(idx2,1);
1556 faceArea42 = interactionVolume.getFaceArea(idx4,0);
1557 faceArea15 = interactionVolume.getFaceArea(idx1,1);
1558 faceArea51 = interactionVolume.getFaceArea(idx5,2);
1559 }
1560 else
1561 {
1562 globalPosFace12 = interactionVolume.getFacePosition(idx1,2);
1563 globalPosFace15 = interactionVolume.getFacePosition(idx1,1);
1564 globalPosFace24 = interactionVolume.getFacePosition(idx2,1);
1565
1566 outerNormaln1 = interactionVolume.getNormal(idx1, 2);
1567 outerNormaln2 = interactionVolume.getNormal(idx2, 1);
1568 outerNormaln3 = interactionVolume.getNormal(idx1, 1);
1569
1570 edgeCoord1215 = interactionVolume.getEdgePosition(idx1, 2, 0);
1571 edgeCoord1315 = interactionVolume.getEdgePosition(idx1, 0, 1);
1572 edgeCoord1224 = interactionVolume.getEdgePosition(idx2, 0, 0);
1573 edgeCoord2426 = interactionVolume.getEdgePosition(idx2, 1, 0);
1574
1575 lambda12 = lambda[idx1][2];
1576 lambda21 = lambda[idx2][0];
1577 lambda24 = lambda[idx2][1];
1578 lambda42 = lambda[idx4][2];
1579 lambda15 = lambda[idx1][1];
1580 lambda51 = lambda[idx5][0];
1581
1582 faceArea12 = interactionVolume.getFaceArea(idx1,2);
1583 faceArea21 = interactionVolume.getFaceArea(idx2,0);
1584 faceArea24 = interactionVolume.getFaceArea(idx2,1);
1585 faceArea42 = interactionVolume.getFaceArea(idx4,2);
1586 faceArea15 = interactionVolume.getFaceArea(idx1,1);
1587 faceArea51 = interactionVolume.getFaceArea(idx5,0);
1588 }
1589
1590 // 3.use non-centered L-stencil 1 to compute the transmissibility of 1/4 face
1591 // compute normal vectors for case 3 (idx1, idx2, idx4, idx5)
1592 DimVector crossProductVector1 = edgeCoord1224-globalPos1;
1593 DimVector crossProductVector2 = globalPosFace15-globalPos1;
1594 DimVector nu11 = crossProduct(crossProductVector1, crossProductVector2);
1595 crossProductVector1 = globalPosFace12-globalPos1;
1596 crossProductVector2 = edgeCoord1224-globalPos1;
1597 DimVector nu12 = crossProduct(crossProductVector1, crossProductVector2);
1598 crossProductVector1 = globalPosFace15-globalPos1;
1599 crossProductVector2 = globalPosFace12-globalPos1;
1600 DimVector nu13 = crossProduct(crossProductVector1, crossProductVector2);
1601
1602 crossProductVector1 = edgeCoord1215-globalPos2;
1603 crossProductVector2 = globalPosFace24-globalPos2;
1604 DimVector nu21 = crossProduct(crossProductVector1, crossProductVector2);
1605 crossProductVector1 = globalPosFace12-globalPos2;
1606 crossProductVector2 = edgeCoord1215-globalPos2;
1607 DimVector nu22 = crossProduct(crossProductVector1, crossProductVector2);
1608 crossProductVector1 = globalPosFace24-globalPos2;
1609 crossProductVector2 = globalPosFace12-globalPos2;
1610 DimVector nu23 = crossProduct(crossProductVector1, crossProductVector2);
1611
1612 crossProductVector1 = edgeCoord2426-globalPos4;
1613 crossProductVector2 = edgeCoord1224-globalPos4;
1614 DimVector nu41 = crossProduct(crossProductVector1, crossProductVector2);
1615 crossProductVector1 = globalPosFace24-globalPos4;
1616 crossProductVector2 = edgeCoord2426-globalPos4;
1617 DimVector nu42 = crossProduct(crossProductVector1, crossProductVector2);
1618 crossProductVector1 = edgeCoord1224-globalPos4;
1619 crossProductVector2 = globalPosFace24-globalPos4;
1620 DimVector nu43 = crossProduct(crossProductVector1, crossProductVector2);
1621
1622 crossProductVector1 = edgeCoord1315-globalPos5;
1623 crossProductVector2 = edgeCoord1215-globalPos5;
1624 DimVector nu51 = crossProduct(crossProductVector1, crossProductVector2);
1625 crossProductVector1 = globalPosFace15-globalPos5;
1626 crossProductVector2 = edgeCoord1315-globalPos5;
1627 DimVector nu52 = crossProduct(crossProductVector1, crossProductVector2);
1628 crossProductVector1 = edgeCoord1215-globalPos5;
1629 crossProductVector2 = globalPosFace15-globalPos5;
1630 DimVector nu53 = crossProduct(crossProductVector1, crossProductVector2);
1631
1632 // compute T, i.e., the volume of the subtetrahedron
1633 crossProductVector1 = edgeCoord1224-globalPos1;
1634 crossProductVector2 = globalPosFace15-globalPos1;
1635 Scalar T1 = (globalPosFace12-globalPos1) * crossProduct(crossProductVector1, crossProductVector2);
1636 crossProductVector1 = edgeCoord1215-globalPos2;
1637 crossProductVector2 = globalPosFace24-globalPos2;
1638 Scalar T2 = (globalPosFace12-globalPos2) * crossProduct(crossProductVector1, crossProductVector2);
1639 crossProductVector1 = edgeCoord2426-globalPos4;
1640 crossProductVector2 = edgeCoord1224-globalPos4;
1641 Scalar T4 = (globalPosFace24-globalPos4) * crossProduct(crossProductVector1, crossProductVector2);
1642 crossProductVector1 = edgeCoord1315-globalPos5;
1643 crossProductVector2 = edgeCoord1215-globalPos5;
1644 Scalar T5 = (globalPosFace15-globalPos5) * crossProduct(crossProductVector1, crossProductVector2);
1645
1646 // compute components needed for flux calculation, denoted as 'omega' and 'r'
1647 DimVector K1nu11(0);
1648 K1.mv(nu11, K1nu11);
1649 DimVector K1nu12(0);
1650 K1.mv(nu12, K1nu12);
1651 DimVector K1nu13(0);
1652 K1.mv(nu13, K1nu13);
1653
1654 DimVector K2nu21(0);
1655 K2.mv(nu21, K2nu21);
1656 DimVector K2nu22(0);
1657 K2.mv(nu22, K2nu22);
1658 DimVector K2nu23(0);
1659 K2.mv(nu23, K2nu23);
1660
1661 DimVector K4nu41(0);
1662 K4.mv(nu41, K4nu41);
1663 DimVector K4nu42(0);
1664 K4.mv(nu42, K4nu42);
1665 DimVector K4nu43(0);
1666 K4.mv(nu43, K4nu43);
1667
1668 DimVector K5nu51(0);
1669 K5.mv(nu51, K5nu51);
1670 DimVector K5nu52(0);
1671 K5.mv(nu52, K5nu52);
1672 DimVector K5nu53(0);
1673 K5.mv(nu53, K5nu53);
1674
1675
1676 Scalar omega111 = lambda12 * (outerNormaln1 * K1nu11) * faceArea12/T1;
1677 Scalar omega112 = lambda12 * (outerNormaln1 * K1nu12) * faceArea12/T1;
1678 Scalar omega113 = lambda12 * (outerNormaln1 * K1nu13) * faceArea12/T1;
1679
1680 Scalar omega121 = lambda21 * (outerNormaln1 * K2nu21) * faceArea21/T2;
1681 Scalar omega122 = lambda21 * (outerNormaln1 * K2nu22) * faceArea21/T2;
1682 Scalar omega123 = lambda21 * (outerNormaln1 * K2nu23) * faceArea21/T2;
1683
1684 Scalar omega221 = lambda24 * (outerNormaln2 * K2nu21) * faceArea24/T2;
1685 Scalar omega222 = lambda24 * (outerNormaln2 * K2nu22) * faceArea24/T2;
1686 Scalar omega223 = lambda24 * (outerNormaln2 * K2nu23) * faceArea24/T2;
1687
1688 Scalar omega241 = lambda42 * (outerNormaln2 * K4nu41) * faceArea42/T4;
1689 Scalar omega242 = lambda42 * (outerNormaln2 * K4nu42) * faceArea42/T4;
1690 Scalar omega243 = lambda42 * (outerNormaln2 * K4nu43) * faceArea42/T4;
1691
1692 Scalar omega311 = lambda15 * (outerNormaln3 * K1nu11) * faceArea15/T1;
1693 Scalar omega312 = lambda15 * (outerNormaln3 * K1nu12) * faceArea15/T1;
1694 Scalar omega313 = lambda15 * (outerNormaln3 * K1nu13) * faceArea15/T1;
1695
1696 Scalar omega351 = lambda51 * (outerNormaln3 * K5nu51) * faceArea51/T5;
1697 Scalar omega352 = lambda51 * (outerNormaln3 * K5nu52) * faceArea51/T5;
1698 Scalar omega353 = lambda51 * (outerNormaln3 * K5nu53) * faceArea51/T5;
1699
1700 Scalar r112 = (nu11 * (edgeCoord1215-globalPos1))/T1;
1701 Scalar r113 = (nu11 * (edgeCoord1315-globalPos1))/T1;
1702 Scalar r122 = (nu12 * (edgeCoord1215-globalPos1))/T1;
1703 Scalar r123 = (nu12 * (edgeCoord1315-globalPos1))/T1;
1704 Scalar r132 = (nu13 * (edgeCoord1215-globalPos1))/T1;
1705 Scalar r133 = (nu13 * (edgeCoord1315-globalPos1))/T1;
1706 Scalar r211 = (nu21 * (edgeCoord1224-globalPos2))/T2;
1707 Scalar r214 = (nu21 * (edgeCoord2426-globalPos2))/T2;
1708 Scalar r221 = (nu22 * (edgeCoord1224-globalPos2))/T2;
1709 Scalar r224 = (nu22 * (edgeCoord2426-globalPos2))/T2;
1710 Scalar r231 = (nu23 * (edgeCoord1224-globalPos2))/T2;
1711 Scalar r234 = (nu23 * (edgeCoord2426-globalPos2))/T2;
1712
1713 Scalar coef = 1.0/(1.0 - r231 * r132);
1714
1715 // compute transmissibility matrix T = CA^{-1}B+D
1716 DimMatrix C(0), A(0);
1717 Dune::FieldMatrix<Scalar,dim,2*dim-dim+1> D(0), B(0);
1718
1719 // evaluate matrix C, D, A, B
1720 C[0][0] = -omega111 - coef * omega113 * r231 * r112 - coef * omega113 * r211;
1721 C[0][1] = -coef * omega113 * r221;
1722 C[0][2] = -omega112 - coef * omega113 * r231 * r122;
1723 C[1][0] = -omega221 - coef * omega223 * r112 - coef * omega223 * r132 * r211;
1724 C[1][1] = -omega222 - coef * omega223 * r132 * r221;
1725 C[1][2] = -coef * omega223 * r122;
1726 C[2][0] = -omega311 - coef * omega313 * r231 * r112 - coef * omega313 * r211;
1727 C[2][1] = -coef * omega313 * r221;
1728 C[2][2] = -omega312 - coef * omega313 * r231 * r122;
1729
1730 D[0][0] = coef * omega113 * r231 * (r112 + r122 + r132 - 1.0) + omega111 + omega112 + omega113;
1731 D[0][1] = coef * omega113 * (r211 + r221 + r231 - 1.0);
1732 D[1][0] = coef * omega223 * (r112 + r122 + r132 - 1.0);
1733 D[1][1] = omega221 + omega222 + omega223 + coef * omega223 * r132 * (r211 + r221 + r231 - 1.0);
1734 D[2][0] = coef * omega313 * r231 * (r112 + r122 + r132 - 1.0) + omega311 + omega312 + omega313;
1735 D[2][1] = coef * omega313 * (r211 + r221 + r231 - 1.0);
1736
1737 A[0][0] = omega111 - omega121 + coef * omega113 * (r231 * r112 + r211) - coef * omega123 * (r112 + r132 * r211);
1738 A[0][1] = coef * omega113 * r221 - omega122 - coef * omega123 * r132 * r221;
1739 A[0][2] = omega112 + coef * omega113 * r231 * r122 - coef * omega123 * r122;
1740 A[1][0] = omega221 - omega243 * r214 + coef * (omega223 - omega243 * r234) * (r112 + r132 * r211) - coef * omega242 * (r231 * r112 + r211);
1741 A[1][1] = omega222 - omega243 * r224 - omega241 + coef * (omega223 - omega243 * r234) * r132 * r221 - coef * omega242 * r221;
1742 A[1][2] = coef * r122 * (omega223 - omega243 * r234 - omega242 * r231);
1743 A[2][0] = omega311 - omega353 * r113 + coef * (omega313 - omega353 * r133) * (r231 * r112 + r211) - coef * omega352 * (r112 + r132 * r211);
1744 A[2][1] = coef * r221 * (omega313 - omega353 * r133 - omega352 * r132);
1745 A[2][2] = omega312 - omega351 - omega353 * r123 + coef * (omega313 - omega353 * r133) * r231 * r122 - coef * omega352 * r122;
1746
1747 B[0][0] = coef * (omega113 * r231 - omega123) * (r112 + r122 + r132 - 1.0) + omega111 + omega112 + omega113;
1748 B[0][1] = coef * (omega113 - omega123 * r132) * (r211 + r221 + r231 - 1.0) - (omega121 + omega122 + omega123);
1749 B[1][0] = coef * (r112 + r122 + r132 - 1.0) * (omega223 - omega243 * r234 - omega242 * r231);
1750 B[1][1] = omega221 + omega222 + omega223 - omega243 * (r214 + r224 + r234 - 1.0) + coef * (r211 + r221 + r231 - 1.0)
1751 * (omega223 * r132 - omega243 * r234 * r132 - omega242);
1752 B[1][2] = -omega241 - omega242 - omega243;
1753 B[2][0] = coef * (r112 + r122 + r132 - 1.0) * (omega313 * r231 - omega353 * r133 * r231 - omega352) + omega311 + omega312 + omega313
1754 - omega353 * (r113 + r123 + r133 - 1.0);
1755 B[2][1] = coef * (r211 + r221 + r231 - 1.0) * (omega313 - omega353 * r133 - omega352 * r132);
1756 B[2][3] = -omega351 - omega352 - omega353;
1757
1758 // compute T
1759 A.invert();
1760 D += B.leftmultiply(C.rightmultiply(A));
1761
1762 transmissibility = D;
1763
1764 if (std::isnan(transmissibility.frobenius_norm()))
1765 {
1766 std::cout<<"case 3: transmissibility = "<<transmissibility<<"\n";
1767 std::cout<<"globalPos1 = "<<globalPos1<<"\n";
1768 std::cout<<"globalPos2 = "<<globalPos2<<"\n";
1769 std::cout<<"globalPos4 = "<<globalPos4<<"\n";
1770 std::cout<<"globalPos5 = "<<globalPos5<<"\n";
1771 std::cout<<"globalPosCenter = "<<globalPosCenter<<"\n";
1772 std::cout<<"outerNormaln1 = "<<outerNormaln1<<"\n";
1773 std::cout<<"outerNormaln2 = "<<outerNormaln2<<"\n";
1774 std::cout<<"outerNormaln3 = "<<outerNormaln3<<"\n";
1775 std::cout<<"xbar_1 = "<<globalPosFace12<<"\n";
1776 std::cout<<"xbar_2 = "<<globalPosFace24<<"\n";
1777 std::cout<<"xbar_3 = "<<globalPosFace15<<"\n";
1778 std::cout<<"xbar_5 = "<<edgeCoord1215<<"\n";
1779 std::cout<<"xbar_6 = "<<edgeCoord1315<<"\n";
1780 std::cout<<"xbar_7 = "<<edgeCoord2426<<"\n";
1781 std::cout<<"perm1 = "<<K1<<"\n";
1782 std::cout<<"perm2 = "<<K2<<"\n";
1783 std::cout<<"perm4 = "<<K4<<"\n";
1784 std::cout<<"perm5 = "<<K5<<"\n";
1785 std::cout<<"lambda = ";
1786 for (unsigned int i = 0; i < lambda.size(); i++)
1787 {
1788 std::cout<<lambda[i]<<" ";
1789 }
1790 std::cout<<"\n";
1791 std::cout<<"\n";
1792 DUNE_THROW(Dune::MathError,"T is nan");
1793 }
1794
1795 return 3;
1796}
1797
1816template<class TypeTag>
1817int FvMpfaL3dTransmissibilityCalculator<TypeTag>::transmissibilityCaseFour(
1818 Dune::FieldMatrix<Scalar,dim,2*dim-dim+1>& transmissibility,
1819 InteractionVolume& interactionVolume,
1820 std::vector<DimVector >& lambda,int idx1, int idx2, int idx3, int idx6)
1821{
1822 auto element1 = interactionVolume.getSubVolumeElement(idx1);
1823 auto element2 = interactionVolume.getSubVolumeElement(idx2);
1824 auto element3 = interactionVolume.getSubVolumeElement(idx3);
1825 auto element6 = interactionVolume.getSubVolumeElement(idx6);
1826
1827 // get global coordinate of cell centers
1828 const GlobalPosition& globalPos1 = element1.geometry().center();
1829 const GlobalPosition& globalPos2 = element2.geometry().center();
1830 const GlobalPosition& globalPos3 = element3.geometry().center();
1831 const GlobalPosition& globalPos6 = element6.geometry().center();
1832
1833 const GlobalPosition& globalPosCenter = interactionVolume.getCenterPosition();
1834
1835 const DimMatrix& K1 = problem_.spatialParams().intrinsicPermeability(element1);
1836 const DimMatrix& K2 = problem_.spatialParams().intrinsicPermeability(element2);
1837 const DimMatrix& K3 = problem_.spatialParams().intrinsicPermeability(element3);
1838 const DimMatrix& K6 = problem_.spatialParams().intrinsicPermeability(element6);
1839
1840 // 1.use centered L-stencil 1 to compute the transmissibility of 1/4 face
1841 GlobalPosition globalPosFace12(0);
1842 GlobalPosition globalPosFace13(0);
1843 GlobalPosition globalPosFace26(0);
1844
1845 DimVector outerNormaln1(0);
1846 DimVector outerNormaln2(0);
1847 DimVector outerNormaln3(0);
1848
1849 GlobalPosition edgeCoord1213(0);
1850 GlobalPosition edgeCoord1315(0);
1851 GlobalPosition edgeCoord1226(0);
1852 GlobalPosition edgeCoord2426(0);
1853
1854 // get lambda and face area
1855 Scalar lambda12 = 0;
1856 Scalar lambda21 = 0;
1857 Scalar lambda13 = 0;
1858 Scalar lambda31 = 0;
1859 Scalar lambda26 = 0;
1860 Scalar lambda62 = 0;
1861 Scalar faceArea12 = 0;
1862 Scalar faceArea21 = 0;
1863 Scalar faceArea13 = 0;
1864 Scalar faceArea31 = 0;
1865 Scalar faceArea26 = 0;
1866 Scalar faceArea62 = 0;
1867
1868 if ((idx1 == 0 && idx2 == 1) || (idx1 == 1 && idx2 == 3) || (idx1 == 3 && idx2 == 2) || (idx1 == 2 && idx2 == 0))
1869 {
1870 globalPosFace12 = interactionVolume.getFacePosition(idx1,0);
1871 globalPosFace13 = interactionVolume.getFacePosition(idx1,1);
1872 globalPosFace26 = interactionVolume.getFacePosition(idx2,2);
1873
1874 outerNormaln1 = interactionVolume.getNormal(idx1, 0);
1875 outerNormaln2 = interactionVolume.getNormal(idx1, 1);
1876 outerNormaln3 = interactionVolume.getNormal(idx2, 2);
1877
1878 edgeCoord1213 = interactionVolume.getEdgePosition(idx1, 0, 1);
1879 edgeCoord1315 = interactionVolume.getEdgePosition(idx1, 1, 1);
1880 edgeCoord1226 = interactionVolume.getEdgePosition(idx2, 1, 1);
1881 edgeCoord2426 = interactionVolume.getEdgePosition(idx2, 0, 0);
1882
1883 lambda12 = lambda[idx1][0];
1884 lambda21 = lambda[idx2][1];
1885 lambda13 = lambda[idx1][1];
1886 lambda31 = lambda[idx3][0];
1887 lambda26 = lambda[idx2][2];
1888 lambda62 = lambda[idx6][0];
1889
1890 faceArea12 = interactionVolume.getFaceArea(idx1,0);
1891 faceArea21 = interactionVolume.getFaceArea(idx2,1);
1892 faceArea13 = interactionVolume.getFaceArea(idx1,1);
1893 faceArea31 = interactionVolume.getFaceArea(idx3,0);
1894 faceArea26 = interactionVolume.getFaceArea(idx2,2);
1895 faceArea62 = interactionVolume.getFaceArea(idx6,0);
1896 }
1897 else if ((idx1 == 5 && idx2 == 4) || (idx1 == 4 && idx2 == 6) || (idx1 == 6 && idx2 == 7) || (idx1 == 7 && idx2 == 5))
1898 {
1899 globalPosFace12 = interactionVolume.getFacePosition(idx1,2);
1900 globalPosFace13 = interactionVolume.getFacePosition(idx1,1);
1901 globalPosFace26 = interactionVolume.getFacePosition(idx2,0);
1902
1903 outerNormaln1 = interactionVolume.getNormal(idx1, 2);
1904 outerNormaln2 = interactionVolume.getNormal(idx1, 1);
1905 outerNormaln3 = interactionVolume.getNormal(idx2, 0);
1906
1907 edgeCoord1213 = interactionVolume.getEdgePosition(idx1, 2, 1);
1908 edgeCoord1315 = interactionVolume.getEdgePosition(idx1, 1, 1);
1909 edgeCoord1226 = interactionVolume.getEdgePosition(idx2, 1, 1);
1910 edgeCoord2426 = interactionVolume.getEdgePosition(idx2, 2, 0);
1911
1912 lambda12 = lambda[idx1][2];
1913 lambda21 = lambda[idx2][1];
1914 lambda13 = lambda[idx1][1];
1915 lambda31 = lambda[idx3][2];
1916 lambda26 = lambda[idx2][0];
1917 lambda62 = lambda[idx6][2];
1918
1919 faceArea12 = interactionVolume.getFaceArea(idx1,2);
1920 faceArea21 = interactionVolume.getFaceArea(idx2,1);
1921 faceArea13 = interactionVolume.getFaceArea(idx1,1);
1922 faceArea31 = interactionVolume.getFaceArea(idx3,2);
1923 faceArea26 = interactionVolume.getFaceArea(idx2,0);
1924 faceArea62 = interactionVolume.getFaceArea(idx6,2);
1925 }
1926 else if ((idx1 == 4 && idx2 == 0) || (idx1 == 7 && idx2 == 3))
1927 {
1928 globalPosFace12 = interactionVolume.getFacePosition(idx1,0);
1929 globalPosFace13 = interactionVolume.getFacePosition(idx1,2);
1930 globalPosFace26 = interactionVolume.getFacePosition(idx2,0);
1931
1932 outerNormaln1 = interactionVolume.getNormal(idx1, 0);
1933 outerNormaln2 = interactionVolume.getNormal(idx1, 2);
1934 outerNormaln3 = interactionVolume.getNormal(idx2, 0);
1935
1936 edgeCoord1213 = interactionVolume.getEdgePosition(idx1, 0, 1);
1937 edgeCoord1315 = interactionVolume.getEdgePosition(idx1, 2, 1);
1938 edgeCoord1226 = interactionVolume.getEdgePosition(idx2, 2, 1);
1939 edgeCoord2426 = interactionVolume.getEdgePosition(idx2, 1, 0);
1940
1941 lambda12 = lambda[idx1][0];
1942 lambda21 = lambda[idx2][2];
1943 lambda13 = lambda[idx1][2];
1944 lambda31 = lambda[idx3][1];
1945 lambda26 = lambda[idx2][0];
1946 lambda62 = lambda[idx6][1];
1947
1948 faceArea12 = interactionVolume.getFaceArea(idx1,0);
1949 faceArea21 = interactionVolume.getFaceArea(idx2,2);
1950 faceArea13 = interactionVolume.getFaceArea(idx1,2);
1951 faceArea31 = interactionVolume.getFaceArea(idx3,1);
1952 faceArea26 = interactionVolume.getFaceArea(idx2,0);
1953 faceArea62 = interactionVolume.getFaceArea(idx6,1);
1954 }
1955 else
1956 {
1957 globalPosFace12 = interactionVolume.getFacePosition(idx1,2);
1958 globalPosFace13 = interactionVolume.getFacePosition(idx1,0);
1959 globalPosFace26 = interactionVolume.getFacePosition(idx2,2);
1960
1961 outerNormaln1 = interactionVolume.getNormal(idx1, 2);
1962 outerNormaln2 = interactionVolume.getNormal(idx1, 0);
1963 outerNormaln3 = interactionVolume.getNormal(idx2, 2);
1964
1965 edgeCoord1213 = interactionVolume.getEdgePosition(idx1, 2, 1);
1966 edgeCoord1315 = interactionVolume.getEdgePosition(idx1, 0, 1);
1967 edgeCoord1226 = interactionVolume.getEdgePosition(idx2, 0, 1);
1968 edgeCoord2426 = interactionVolume.getEdgePosition(idx2, 1, 0);
1969
1970 lambda12 = lambda[idx1][2];
1971 lambda21 = lambda[idx2][0];
1972 lambda13 = lambda[idx1][0];
1973 lambda31 = lambda[idx3][1];
1974 lambda26 = lambda[idx2][2];
1975 lambda62 = lambda[idx6][1];
1976
1977 faceArea12 = interactionVolume.getFaceArea(idx1,2);
1978 faceArea21 = interactionVolume.getFaceArea(idx2,0);
1979 faceArea13 = interactionVolume.getFaceArea(idx1,0);
1980 faceArea31 = interactionVolume.getFaceArea(idx3,1);
1981 faceArea26 = interactionVolume.getFaceArea(idx2,2);
1982 faceArea62 = interactionVolume.getFaceArea(idx6,1);
1983 }
1984
1985 // 4.use non-centered L-stencil 2 to compute the transmissibility of 1/4 face
1986 // compute normal vectors for case 4 (idx1, idx2, idx3, idx6)
1987 DimVector crossProductVector1 = globalPosFace13-globalPos1;
1988 DimVector crossProductVector2 = edgeCoord1226-globalPos1;
1989 DimVector nu11 = crossProduct(crossProductVector1, crossProductVector2);
1990 crossProductVector1 = edgeCoord1226-globalPos1;
1991 crossProductVector2 = globalPosFace12-globalPos1;
1992 DimVector nu12 = crossProduct(crossProductVector1, crossProductVector2);
1993 crossProductVector1 = globalPosFace12-globalPos1;
1994 crossProductVector2 = globalPosFace13-globalPos1;
1995 DimVector nu13 = crossProduct(crossProductVector1, crossProductVector2);
1996
1997 crossProductVector1 = globalPosFace26-globalPos2;
1998 crossProductVector2 = edgeCoord1213-globalPos2;
1999 DimVector nu21 = crossProduct(crossProductVector1, crossProductVector2);
2000 crossProductVector1 = edgeCoord1213-globalPos2;
2001 crossProductVector2 = globalPosFace12-globalPos2;
2002 DimVector nu22 = crossProduct(crossProductVector1, crossProductVector2);
2003 crossProductVector1 = globalPosFace12-globalPos2;
2004 crossProductVector2 = globalPosFace26-globalPos2;
2005 DimVector nu23 = crossProduct(crossProductVector1, crossProductVector2);
2006
2007 crossProductVector1 = edgeCoord1213-globalPos3;
2008 crossProductVector2 = edgeCoord1315-globalPos3;
2009 DimVector nu31 = crossProduct(crossProductVector1, crossProductVector2);
2010 crossProductVector1 = edgeCoord1315-globalPos3;
2011 crossProductVector2 = globalPosFace13-globalPos3;
2012 DimVector nu32 = crossProduct(crossProductVector1, crossProductVector2);
2013 crossProductVector1 = globalPosFace13-globalPos3;
2014 crossProductVector2 = edgeCoord1213-globalPos3;
2015 DimVector nu33 = crossProduct(crossProductVector1, crossProductVector2);
2016
2017 crossProductVector1 = edgeCoord1226-globalPos6;
2018 crossProductVector2 = edgeCoord2426-globalPos6;
2019 DimVector nu61 = crossProduct(crossProductVector1, crossProductVector2);
2020 crossProductVector1 = edgeCoord2426-globalPos6;
2021 crossProductVector2 = globalPosFace26-globalPos6;
2022 DimVector nu62 = crossProduct(crossProductVector1, crossProductVector2);
2023 crossProductVector1 = globalPosFace26-globalPos6;
2024 crossProductVector2 = edgeCoord1226-globalPos6;
2025 DimVector nu63 = crossProduct(crossProductVector1, crossProductVector2);
2026
2027 // compute T, i.e., the volume of the subtetrahedron
2028 crossProductVector1 = globalPosFace13-globalPos1;
2029 crossProductVector2 = edgeCoord1226-globalPos1;
2030 Scalar T1 = (globalPosFace12-globalPos1) * crossProduct(crossProductVector1, crossProductVector2);
2031 crossProductVector1 = globalPosFace26-globalPos2;
2032 crossProductVector2 = edgeCoord1213-globalPos2;
2033 Scalar T2 = (globalPosFace12-globalPos2) * crossProduct(crossProductVector1, crossProductVector2);
2034 crossProductVector1 = edgeCoord1213-globalPos3;
2035 crossProductVector2 = edgeCoord1315-globalPos3;
2036 Scalar T3 = (globalPosFace13-globalPos3) * crossProduct(crossProductVector1, crossProductVector2);
2037 crossProductVector1 = edgeCoord1226-globalPos6;
2038 crossProductVector2 = edgeCoord2426-globalPos6;
2039 Scalar T6 = (globalPosFace26-globalPos6) * crossProduct(crossProductVector1, crossProductVector2);
2040
2041 // compute components needed for flux calculation, denoted as 'omega' and 'r'
2042 DimVector K1nu11(0);
2043 K1.mv(nu11, K1nu11);
2044 DimVector K1nu12(0);
2045 K1.mv(nu12, K1nu12);
2046 DimVector K1nu13(0);
2047 K1.mv(nu13, K1nu13);
2048
2049 DimVector K2nu21(0);
2050 K2.mv(nu21, K2nu21);
2051 DimVector K2nu22(0);
2052 K2.mv(nu22, K2nu22);
2053 DimVector K2nu23(0);
2054 K2.mv(nu23, K2nu23);
2055
2056 DimVector K3nu31(0);
2057 K3.mv(nu31, K3nu31);
2058 DimVector K3nu32(0);
2059 K3.mv(nu32, K3nu32);
2060 DimVector K3nu33(0);
2061 K3.mv(nu33, K3nu33);
2062
2063 DimVector K6nu61(0);
2064 K6.mv(nu61, K6nu61);
2065 DimVector K6nu62(0);
2066 K6.mv(nu62, K6nu62);
2067 DimVector K6nu63(0);
2068 K6.mv(nu63, K6nu63);
2069
2070 // std::cout<<"outerNormaln2 = "<<outerNormaln2<<"\n";
2071 // std::cout<<"outerNormaln3 = "<<outerNormaln3<<"\n";
2072
2073 Scalar omega111 = lambda12 * (outerNormaln1 * K1nu11) * faceArea12/T1;
2074 Scalar omega112 = lambda12 * (outerNormaln1 * K1nu12) * faceArea12/T1;
2075 Scalar omega113 = lambda12 * (outerNormaln1 * K1nu13) * faceArea12/T1;
2076
2077 Scalar omega121 = lambda21 * (outerNormaln1 * K2nu21) * faceArea21/T2;
2078 Scalar omega122 = lambda21 * (outerNormaln1 * K2nu22) * faceArea21/T2;
2079 Scalar omega123 = lambda21 * (outerNormaln1 * K2nu23) * faceArea21/T2;
2080
2081 Scalar omega211 = lambda13 * (outerNormaln2 * K1nu11) * faceArea13/T1;
2082 Scalar omega212 = lambda13 * (outerNormaln2 * K1nu12) * faceArea13/T1;
2083 Scalar omega213 = lambda13 * (outerNormaln2 * K1nu13) * faceArea13/T1;
2084
2085 Scalar omega231 = lambda31 * (outerNormaln2 * K3nu31) * faceArea31/T3;
2086 Scalar omega232 = lambda31 * (outerNormaln2 * K3nu32) * faceArea31/T3;
2087 Scalar omega233 = lambda31 * (outerNormaln2 * K3nu33) * faceArea31/T3;
2088
2089 Scalar omega321 = lambda26 * (outerNormaln3 * K2nu21) * faceArea26/T2;
2090 Scalar omega322 = lambda26 * (outerNormaln3 * K2nu22) * faceArea26/T2;
2091 Scalar omega323 = lambda26 * (outerNormaln3 * K2nu23) * faceArea26/T2;
2092
2093 Scalar omega361 = lambda62 * (outerNormaln3 * K6nu61) * faceArea62/T6;
2094 Scalar omega362 = lambda62 * (outerNormaln3 * K6nu62) * faceArea62/T6;
2095 Scalar omega363 = lambda62 * (outerNormaln3 * K6nu63) * faceArea62/T6;
2096
2097 Scalar r111 = (nu11 * (edgeCoord1213-globalPos1))/T1;
2098 Scalar r114 = (nu11 * (edgeCoord1315-globalPos1))/T1;
2099 Scalar r121 = (nu12 * (edgeCoord1213-globalPos1))/T1;
2100 Scalar r124 = (nu12 * (edgeCoord1315-globalPos1))/T1;
2101 Scalar r131 = (nu13 * (edgeCoord1213-globalPos1))/T1;
2102 Scalar r134 = (nu13 * (edgeCoord1315-globalPos1))/T1;
2103
2104 Scalar r212 = (nu21 * (edgeCoord1226-globalPos2))/T2;
2105 Scalar r213 = (nu21 * (edgeCoord2426-globalPos2))/T2;
2106 Scalar r222 = (nu22 * (edgeCoord1226-globalPos2))/T2;
2107 Scalar r223 = (nu22 * (edgeCoord2426-globalPos2))/T2;
2108 Scalar r232 = (nu23 * (edgeCoord1226-globalPos2))/T2;
2109 Scalar r233 = (nu23 * (edgeCoord2426-globalPos2))/T2;
2110
2111 Scalar coef = 1.0/(1.0 - r232 * r131);
2112
2113 // compute transmissibility matrix TC1 = CA^{-1}B+D
2114 DimMatrix C(0), A(0);
2115 Dune::FieldMatrix<Scalar,dim,2*dim-dim+1> D(0), B(0);
2116
2117 // evaluate matrix C, D, A, B
2118 C[0][0] = -omega111 - coef * omega113 * r212 - coef * omega113 * r232 * r111;
2119 C[0][1] = -omega112 - coef * omega113 * r232 * r121;
2120 C[0][2] = -coef * omega113 * r222;
2121 C[1][0] = -omega211 - coef * omega213 * r212 - coef * omega213 * r232 * r111;
2122 C[1][1] = -omega212 - coef * omega213 * r232 * r121;
2123 C[1][2] = -coef * omega213 * r222;
2124 C[2][0] = -omega321 - coef * omega323 * r111 - coef * omega323 * r131 * r212;
2125 C[2][1] = -coef * omega323 * r121;
2126 C[2][2] = -omega322 - coef * omega323 * r131 * r222;
2127
2128 D[0][0] = coef * omega113 * r232 * (r111 + r121 + r131 - 1.0) + omega111 + omega112 + omega113;
2129 D[0][1] = coef * omega113 * (r212 + r222 + r232 - 1.0);
2130 D[1][0] = coef * omega213 * r232 * (r111 + r121 + r131 - 1.0) + omega211 + omega212 + omega213;
2131 D[1][1] = coef * omega213 * (r212 + r222 + r232 - 1.0);
2132 D[2][0] = coef * omega323 * (r111 + r121 + r131 - 1.0);
2133 D[2][1] = coef * omega323 * r131 * (r212 + r222 + r232 - 1.0) + omega321 + omega322 + omega323;
2134
2135 A[0][0] = omega111 - omega121 + coef * omega113 * (r212 + r232 * r111) - coef * omega123 * (r111 + r131 * r212);
2136 A[0][1] = omega112 + coef * omega113 * r232 * r121 - coef * omega123 * r121;
2137 A[0][2] = coef * omega113 * r222 - omega122 - coef * omega123 * r131 * r222;
2138 A[1][0] = omega211 - omega233 * r114 + coef * (omega213 - omega233 * r134) * (r212 + r232 * r111) - coef * omega232 * (r111 + r131 * r212);
2139 A[1][1] = omega212 - omega231 - omega233 * r124 + coef * (omega213 - omega233 * r134) * r232 * r121 - coef * omega232 * r121;
2140 A[1][2] = coef * r222 * (omega213 - omega233 * r134 - omega232 * r131);
2141 A[2][0] = omega321 - omega363 * r213 + coef * (omega323 - omega363 * r233) * (r111 + r131 * r212) - coef * omega362 * (r212 + r232 * r111);
2142 A[2][1] = coef * r121 * (omega323 - omega363 * r233 - omega362 * r232);
2143 A[2][2] = omega322 - omega361 - omega363 * r223 + coef * (omega323 - omega363 * r233) * r131 * r222 - coef * omega362 * r222;
2144
2145 B[0][0] = coef * (omega113 * r232 - omega123) * (r111 + r121 + r131 - 1.0) + omega111 + omega112 + omega113;
2146 B[0][1] = coef * (omega113 - omega123 * r131) * (r212 + r222 + r232 - 1.0) - (omega121 + omega122 + omega123);
2147 B[1][0] = omega211 + omega212 + omega213 - omega233 * (r114 + r124 + r134 - 1.0) + coef * (r111 + r121 + r131 - 1.0)
2148 * (omega213 * r232 - omega233 * r134 * r232 - omega232);
2149 B[1][1] = coef * (r212 + r222 + r232 - 1.0) * (omega213 - omega233 * r134 - omega232 * r131);
2150 B[1][2] = -omega231 - omega232 - omega233;
2151 B[2][0] = coef * (r111 + r121 + r131 - 1.0) * (omega323 - omega363 * r233 - omega362 * r232);
2152 B[2][1] = omega321 + omega322 + omega323 + coef * (r212 + r222 + r232 - 1.0) * (omega323 * r131 - omega363 * r233 * r131 - omega362)
2153 - omega363 * (r213 + r223 + r233 - 1.0);
2154 B[2][3] = -omega361 - omega362 - omega363;
2155
2156 // compute T
2157 A.invert();
2158 D += B.leftmultiply(C.rightmultiply(A));
2159
2160 transmissibility = D;
2161
2162 if (std::isnan(transmissibility.frobenius_norm()))
2163 {
2164 std::cout<<"case 4: transmissibility = "<<transmissibility<<"\n";
2165
2166 std::cout<<"globalPos1 = "<<globalPos1<<"\n";
2167 std::cout<<"globalPos2 = "<<globalPos2<<"\n";
2168 std::cout<<"globalPos3 = "<<globalPos3<<"\n";
2169 std::cout<<"globalPos6 = "<<globalPos6<<"\n";
2170
2171 std::cout<<"globalPosCenter = "<<globalPosCenter<<"\n";
2172 std::cout<<"outerNormaln1 = "<<outerNormaln1<<"\n";
2173 std::cout<<"outerNormaln2 = "<<outerNormaln2<<"\n";
2174 std::cout<<"outerNormaln3 = "<<outerNormaln3<<"\n";
2175 std::cout<<"xbar_1 = "<<globalPosFace12<<"\n";
2176 std::cout<<"xbar_2 = "<<globalPosFace13<<"\n";
2177 std::cout<<"xbar_3 = "<<globalPosFace26<<"\n";
2178 std::cout<<"xbar_4 = "<<edgeCoord1213<<"\n";
2179 std::cout<<"xbar_5 = "<<edgeCoord1226<<"\n";
2180 std::cout<<"xbar_6 = "<<edgeCoord2426<<"\n";
2181 std::cout<<"xbar_7 = "<<edgeCoord1315<<"\n";
2182 std::cout<<"perm1 = "<<K1<<"\n";
2183 std::cout<<"perm2 = "<<K2<<"\n";
2184 std::cout<<"perm3 = "<<K3<<"\n";
2185 std::cout<<"perm6 = "<<K6<<"\n";
2186 std::cout<<"lambda = ";
2187 for (unsigned int i = 0; i < lambda.size(); i++)
2188 {
2189 std::cout<<lambda[i]<<" ";
2190 }
2191 std::cout<<"\n";
2192 std::cout<<"\n";
2193 DUNE_THROW(Dune::MathError,"T is nan");
2194 }
2195
2196 return 4;
2197}
2198} // end namespace Dumux
2199#endif
GET_PROP_TYPE
#define GET_PROP_TYPE(TypeTag, PropTagName)
Definition: propertysystemmacros.hh:283
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
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
Dumux::Properties::MPFAInteractionVolume
Property tag MPFAInteractionVolume
Type of the data container for one interaction volume.
Definition: porousmediumflow/sequential/cellcentered/mpfa/properties.hh:96
Dumux::IOName::permeability
std::string permeability() noexcept
I/O name of permeability.
Definition: name.hh:143
Dumux::FvMpfaL3dTransmissibilityCalculator
Provides methods for transmissibility calculation in 3-d.
Definition: 3dtransmissibilitycalculator.hh:51
Dumux::FvMpfaL3dTransmissibilityCalculator::transmissibilityCaseFour
int transmissibilityCaseFour(Dune::FieldMatrix< Scalar, dim, 2 *dim-dim+1 > &transmissibility, InteractionVolume &interactionVolume, std::vector< DimVector > &lambda, int idx1, int idx2, int idx3, int idx6)
Calculates the transmissibility matrix for the flux face between the cells with the local index idx1 ...
Definition: 3dtransmissibilitycalculator.hh:1817
Dumux::FvMpfaL3dTransmissibilityCalculator::transmissibilityCaseThree
int transmissibilityCaseThree(Dune::FieldMatrix< Scalar, dim, 2 *dim-dim+1 > &transmissibility, InteractionVolume &interactionVolume, std::vector< DimVector > &lambda, int idx1, int idx2, int idx4, int idx5)
Calculates the transmissibility matrix for the flux face between the cells with the local index idx1 ...
Definition: 3dtransmissibilitycalculator.hh:1422
Dumux::FvMpfaL3dTransmissibilityCalculator::transmissibility
int transmissibility(Dune::FieldMatrix< Scalar, dim, 2 *dim-dim+1 > &transmissibility, InteractionVolume &interactionVolume, std::vector< DimVector > &lambda, int idx1, int idx2, int idx3, int idx4, int idx5, int idx6, Dune::FieldVector< bool, 4 > &useCases)
Dumux::FvMpfaL3dTransmissibilityCalculator::TransmissibilityType
Dune::FieldMatrix< Scalar, dim, 2 *dim - dim+1 > TransmissibilityType
Type of the transmissibility matrix.
Definition: 3dtransmissibilitycalculator.hh:78
Dumux::FvMpfaL3dTransmissibilityCalculator::transmissibilityTPFA
int transmissibilityTPFA(Dune::FieldMatrix< Scalar, dim, 2 *dim-dim+1 > &transmissibility, InteractionVolume &interactionVolume, std::vector< DimVector > &lambda, int idx1, int idx2)
Calculates a TPFA transmissibility matrix for the flux face between the cells with the local index id...
Definition: 3dtransmissibilitycalculator.hh:540
Dumux::FvMpfaL3dTransmissibilityCalculator::transmissibilityCaseOne
int transmissibilityCaseOne(Dune::FieldMatrix< Scalar, dim, 2 *dim-dim+1 > &transmissibility, InteractionVolume &interactionVolume, std::vector< DimVector > &lambda, int idx1, int idx2, int idx3, int idx5)
Calculates the transmissibility matrix for the flux face between the cells with the local index idx1 ...
Definition: 3dtransmissibilitycalculator.hh:645
Dumux::FvMpfaL3dTransmissibilityCalculator::chooseTransmissibility
int chooseTransmissibility(TransmissibilityType &transmissibilityOne, TransmissibilityType &transmissibilityTwo, int lTypeOne, int lTypeTwo)
Compares two transmissibility matrices according to a L-selection criterion.
Definition: 3dtransmissibilitycalculator.hh:186
Dumux::FvMpfaL3dTransmissibilityCalculator::FvMpfaL3dTransmissibilityCalculator
FvMpfaL3dTransmissibilityCalculator(Problem &problem)
Constructs a FvMpfaL3dTransmissibilityCalculator object.
Definition: 3dtransmissibilitycalculator.hh:127
Dumux::FvMpfaL3dTransmissibilityCalculator::transmissibilityCaseTwo
int transmissibilityCaseTwo(Dune::FieldMatrix< Scalar, dim, 2 *dim-dim+1 > &transmissibility, InteractionVolume &interactionVolume, std::vector< DimVector > &lambda, int idx1, int idx2, int idx4, int idx6)
Calculates the transmissibility matrix for the flux face between the cells with the local index idx1 ...
Definition: 3dtransmissibilitycalculator.hh:1027
Dumux::FvMpfaL3dTransmissibilityCalculator::transmissibility
int transmissibility(Dune::FieldMatrix< Scalar, dim, 2 *dim-dim+1 > &transmissibility, InteractionVolume &interactionVolume, std::vector< DimVector > &lambda, int idx1, int idx2, int idx3, int idx4, int idx5, int idx6)
properties.hh
Properties for a MPFA method.
pressureproperties.hh
Base file for properties related to sequential IMPET algorithms.
Generated by doxygen 1.9.3