3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
3dvelocity.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_FVMPFAL2PVELOCITY2P_HH
25#define DUMUX_FVMPFAL2PVELOCITY2P_HH
26
27#include <dune/grid/common/gridenums.hh>
28#include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
29#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
30#include "3dtransmissibilitycalculator.hh"
31
32#include <dumux/common/deprecated.hh>
33
34namespace Dumux {
35
56template<class TypeTag> class FvMpfaL3dVelocity2p
57{
58 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
59 enum
60 {
61 dim = GridView::dimension, dimWorld = GridView::dimensionworld
62 };
63
64 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
65 using Problem = GetPropType<TypeTag, Properties::Problem>;
66
67 using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
68
69 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
70
71 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
72 using FluidState = GetPropType<TypeTag, Properties::FluidState>;
73
74 using BoundaryTypes = GetPropType<TypeTag, Properties::SequentialBoundaryTypes>;
75 using CellData = GetPropType<TypeTag, Properties::CellData>;
76
77 using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
78 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
79
80 using Element = typename GridView::Traits::template Codim<0>::Entity;
81 using Grid = typename GridView::Grid;
82 using IndexSet = typename GridView::IndexSet;
83
84 using Geometry = typename Element::Geometry;
85 using JacobianTransposed = typename Geometry::JacobianTransposed ;
86
87 using GridTypeIndices = GetPropType<TypeTag, Properties::GridTypeIndices>;
88
89 using InteractionVolume = GetPropType<TypeTag, Properties::MPFAInteractionVolume>;
90 using InteractionVolumeContainer = GetPropType<TypeTag, Properties::MPFAInteractionVolumeContainer>;
91 using TransmissibilityCalculator = FvMpfaL3dTransmissibilityCalculator<TypeTag>;
92 using TransmissibilityType = typename TransmissibilityCalculator::TransmissibilityType;
93
94 enum
95 {
96 pw = Indices::pressureW,
97 pn = Indices::pressureNw,
98 pGlobal = Indices::pressureGlobal,
99 sw = Indices::saturationW,
100 sn = Indices::saturationNw,
101 vw = Indices::velocityW,
102 vn = Indices::velocityNw,
103 vt = Indices::velocityTotal
104 };
105 enum
106 {
107 wPhaseIdx = Indices::wPhaseIdx,
108 nPhaseIdx = Indices::nPhaseIdx,
109 pressureIdx = Indices::pressureIdx,
110 saturationIdx = Indices::saturationIdx,
111 pressureEqIdx = Indices::pressureEqIdx,
112 satEqIdx = Indices::satEqIdx,
113 numPhases = getPropValue<TypeTag, Properties::NumPhases>()
114 };
115 enum
116 {
117 globalCorner = 2,
118 globalEdge = 3,
119 neumannNeumann = 0,
120 dirichletDirichlet = 1,
121 dirichletNeumann = 2,
122 neumannDirichlet = 3
123 };
124 enum
125 {
126 innerEdgeFace = 2, innerSideFace = 1
127 };
128
129 using GlobalPosition = typename Geometry::GlobalCoordinate;
130 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
131 using DimVector = Dune::FieldVector<Scalar, dim>;
132
133public:
139 FvMpfaL3dVelocity2p(Problem& problem) :
140 problem_(problem), gravity_(problem.gravity())
141 {
142 density_[wPhaseIdx] = 0.;
143 density_[nPhaseIdx] = 0.;
144 viscosity_[wPhaseIdx] = 0.;
145 viscosity_[nPhaseIdx] = 0.;
146
147 vtkOutputLevel_ = getParam<int>("Vtk.OutputLevel");
148 }
149
151 void calculateInnerInteractionVolumeVelocity(InteractionVolume& interactionVolume,
152 CellData & cellData1, CellData & cellData2, CellData & cellData3, CellData & cellData4,
153 CellData & cellData5, CellData & cellData6, CellData & cellData7, CellData & cellData8,
154 InteractionVolumeContainer& interactionVolumes,
155 TransmissibilityCalculator& transmissibilityCalculator, int fIdx = -1);
156 void calculateBoundaryInteractionVolumeVelocity(InteractionVolume& interactionVolume,
157 CellData& cellData, int elemIdx);
158
160 void initialize(bool solveTwice = true)
161 {
162 const auto element = *problem_.gridView().template begin<0>();
163 FluidState fluidState;
164 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
165 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
166 fluidState.setTemperature(problem_.temperature(element));
167 fluidState.setSaturation(wPhaseIdx, 1.);
168 fluidState.setSaturation(nPhaseIdx, 0.);
169 density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
170 density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
171 viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
172 viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
173
174 return;
175 }
176
187 template<class MultiWriter>
188 void addOutputVtkFields(MultiWriter &writer)
189 {
190 if (vtkOutputLevel_ > 0)
191 {
192 Dune::BlockVector < DimVector > &velocityWetting = *(writer.template allocateManagedBuffer<Scalar,
193 dim>(problem_.gridView().size(0)));
194 Dune::BlockVector < DimVector > &velocityNonwetting = *(writer.template allocateManagedBuffer<Scalar,
195 dim>(problem_.gridView().size(0)));
196
197 // compute update vector
198 for (const auto& element : elements(problem_.gridView()))
199 {
200 // cell index
201 int eIdxGlobal = problem_.variables().index(element);
202
203 CellData & cellData = problem_.variables().cellData(eIdxGlobal);
204
205 Dune::FieldVector < Scalar, 2 * dim > fluxW(0);
206 Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
207
208 // run through all intersections with neighbors and boundary
209 for (const auto& intersection : intersections(problem_.gridView(), element))
210 {
211 int isIndex = intersection.indexInInside();
212
213 fluxW[isIndex] += intersection.geometry().volume()
214 * (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
215 fluxNw[isIndex] += intersection.geometry().volume()
216 * (intersection.centerUnitOuterNormal() * cellData.fluxData().velocity(nPhaseIdx, isIndex));
217 }
218
219 DimVector refVelocity(0);
220 refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
221 refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
222 refVelocity[2] = 0.5 * (fluxW[5] - fluxW[4]);
223
224 const DimVector& localPos = referenceElement(element).position(0, 0);
225
226 // get the transposed Jacobian of the element mapping
227 const JacobianTransposed jacobianT = element.geometry().jacobianTransposed(localPos);
228
229 // calculate the element velocity by the Piola transformation
230 DimVector elementVelocity(0);
231 jacobianT.umtv(refVelocity, elementVelocity);
232 elementVelocity /= element.geometry().integrationElement(localPos);
233
234 velocityWetting[eIdxGlobal] = elementVelocity;
235
236 refVelocity = 0;
237 refVelocity[0] = 0.5 * (fluxNw[1] - fluxNw[0]);
238 refVelocity[1] = 0.5 * (fluxNw[3] - fluxNw[2]);
239 refVelocity[2] = 0.5 * (fluxNw[5] - fluxNw[4]);
240
241 // calculate the element velocity by the Piola transformation
242 elementVelocity = 0;
243 jacobianT.umtv(refVelocity, elementVelocity);
244 elementVelocity /= element.geometry().integrationElement(localPos);
245
246 velocityNonwetting[eIdxGlobal] = elementVelocity;
247 }
248
249 writer.attachCellData(velocityWetting, "wetting-velocity", dim);
250 writer.attachCellData(velocityNonwetting, "nonwetting-velocity", dim);
251 }
252
253 return;
254 }
255
256private:
257 Problem& problem_;
258 const Dune::FieldVector<Scalar, dimWorld>& gravity_;
259
260 Scalar density_[numPhases];
261 Scalar viscosity_[numPhases];
262
263 int vtkOutputLevel_;
264
265 static constexpr Scalar threshold_ = 1e-15;
267 static const int velocityType_ = getPropValue<TypeTag, Properties::VelocityFormulation>();
269 static const int pressureType_ = getPropValue<TypeTag, Properties::PressureFormulation>();
271 static const int saturationType_ = getPropValue<TypeTag, Properties::SaturationFormulation>();
272};
273// end of template
274
295template<class TypeTag>
296void FvMpfaL3dVelocity2p<TypeTag>::calculateInnerInteractionVolumeVelocity(InteractionVolume& interactionVolume,
297 CellData & cellData1, CellData & cellData2, CellData & cellData3, CellData & cellData4,
298 CellData & cellData5, CellData & cellData6, CellData & cellData7, CellData & cellData8,
299 InteractionVolumeContainer& interactionVolumes, TransmissibilityCalculator& transmissibilityCalculator, int fIdx)
300 {
301 auto element1 = interactionVolume.getSubVolumeElement(0);
302 auto element2 = interactionVolume.getSubVolumeElement(1);
303 auto element3 = interactionVolume.getSubVolumeElement(2);
304 auto element4 = interactionVolume.getSubVolumeElement(3);
305 auto element5 = interactionVolume.getSubVolumeElement(4);
306 auto element6 = interactionVolume.getSubVolumeElement(5);
307 auto element7 = interactionVolume.getSubVolumeElement(6);
308 auto element8 = interactionVolume.getSubVolumeElement(7);
309
310 // cell index
311 int eIdxGlobal1 = problem_.variables().index(element1);
312 int eIdxGlobal2 = problem_.variables().index(element2);
313 int eIdxGlobal3 = problem_.variables().index(element3);
314 int eIdxGlobal4 = problem_.variables().index(element4);
315 int eIdxGlobal5 = problem_.variables().index(element5);
316 int eIdxGlobal6 = problem_.variables().index(element6);
317 int eIdxGlobal7 = problem_.variables().index(element7);
318 int eIdxGlobal8 = problem_.variables().index(element8);
319
320 // pressures flux calculation
321 Dune::FieldVector<Scalar, 8> potW(0);
322 Dune::FieldVector<Scalar, 8> potNw(0);
323
324 potW[0] = cellData1.potential(wPhaseIdx);
325 potW[1] = cellData2.potential(wPhaseIdx);
326 potW[2] = cellData3.potential(wPhaseIdx);
327 potW[3] = cellData4.potential(wPhaseIdx);
328 potW[4] = cellData5.potential(wPhaseIdx);
329 potW[5] = cellData6.potential(wPhaseIdx);
330 potW[6] = cellData7.potential(wPhaseIdx);
331 potW[7] = cellData8.potential(wPhaseIdx);
332
333 potNw[0] = cellData1.potential(nPhaseIdx);
334 potNw[1] = cellData2.potential(nPhaseIdx);
335 potNw[2] = cellData3.potential(nPhaseIdx);
336 potNw[3] = cellData4.potential(nPhaseIdx);
337 potNw[4] = cellData5.potential(nPhaseIdx);
338 potNw[5] = cellData6.potential(nPhaseIdx);
339 potNw[6] = cellData7.potential(nPhaseIdx);
340 potNw[7] = cellData8.potential(nPhaseIdx);
341
342 // get mobilities of the phases
343 Dune::FieldVector<Scalar, numPhases> lambda1(cellData1.mobility(wPhaseIdx));
344 lambda1[nPhaseIdx] = cellData1.mobility(nPhaseIdx);
345
346 // compute total mobility of cell 1
347 Scalar lambdaTotal1 = lambda1[wPhaseIdx] + lambda1[nPhaseIdx];
348
349 // get mobilities of the phases
350 Dune::FieldVector<Scalar, numPhases> lambda2(cellData2.mobility(wPhaseIdx));
351 lambda2[nPhaseIdx] = cellData2.mobility(nPhaseIdx);
352
353 // compute total mobility of cell 2
354 Scalar lambdaTotal2 = lambda2[wPhaseIdx] + lambda2[nPhaseIdx];
355
356 // get mobilities of the phases
357 Dune::FieldVector<Scalar, numPhases> lambda3(cellData3.mobility(wPhaseIdx));
358 lambda3[nPhaseIdx] = cellData3.mobility(nPhaseIdx);
359
360 // compute total mobility of cell 3
361 Scalar lambdaTotal3 = lambda3[wPhaseIdx] + lambda3[nPhaseIdx];
362
363 // get mobilities of the phases
364 Dune::FieldVector<Scalar, numPhases> lambda4(cellData4.mobility(wPhaseIdx));
365 lambda4[nPhaseIdx] = cellData4.mobility(nPhaseIdx);
366
367 // compute total mobility of cell 4
368 Scalar lambdaTotal4 = lambda4[wPhaseIdx] + lambda4[nPhaseIdx];
369
370 // get mobilities of the phases
371 Dune::FieldVector<Scalar, numPhases> lambda5(cellData5.mobility(wPhaseIdx));
372 lambda5[nPhaseIdx] = cellData5.mobility(nPhaseIdx);
373
374 // compute total mobility of cell 5
375 Scalar lambdaTotal5 = lambda5[wPhaseIdx] + lambda5[nPhaseIdx];
376
377 // get mobilities of the phases
378 Dune::FieldVector<Scalar, numPhases> lambda6(cellData6.mobility(wPhaseIdx));
379 lambda6[nPhaseIdx] = cellData6.mobility(nPhaseIdx);
380
381 // compute total mobility of cell 6
382 Scalar lambdaTotal6 = lambda6[wPhaseIdx] + lambda6[nPhaseIdx];
383
384 // get mobilities of the phases
385 Dune::FieldVector<Scalar, numPhases> lambda7(cellData7.mobility(wPhaseIdx));
386 lambda7[nPhaseIdx] = cellData7.mobility(nPhaseIdx);
387
388 // compute total mobility of cell 7
389 Scalar lambdaTotal7 = lambda7[wPhaseIdx] + lambda7[nPhaseIdx];
390
391 // get mobilities of the phases
392 Dune::FieldVector<Scalar, numPhases> lambda8(cellData8.mobility(wPhaseIdx));
393 lambda8[nPhaseIdx] = cellData8.mobility(nPhaseIdx);
394
395 // compute total mobility of cell 8
396 Scalar lambdaTotal8 = lambda8[wPhaseIdx] + lambda8[nPhaseIdx];
397
398 std::vector<DimVector> lambda(8);
399 lambda[0][0] = lambdaTotal1;
400 lambda[0][1] = lambdaTotal1;
401 lambda[0][2] = lambdaTotal1;
402 lambda[1][0] = lambdaTotal2;
403 lambda[1][1] = lambdaTotal2;
404 lambda[1][2] = lambdaTotal2;
405 lambda[2][0] = lambdaTotal3;
406 lambda[2][1] = lambdaTotal3;
407 lambda[2][2] = lambdaTotal3;
408 lambda[3][0] = lambdaTotal4;
409 lambda[3][1] = lambdaTotal4;
410 lambda[3][2] = lambdaTotal4;
411 lambda[4][0] = lambdaTotal5;
412 lambda[4][1] = lambdaTotal5;
413 lambda[4][2] = lambdaTotal5;
414 lambda[5][0] = lambdaTotal6;
415 lambda[5][1] = lambdaTotal6;
416 lambda[5][2] = lambdaTotal6;
417 lambda[6][0] = lambdaTotal7;
418 lambda[6][1] = lambdaTotal7;
419 lambda[6][2] = lambdaTotal7;
420 lambda[7][0] = lambdaTotal8;
421 lambda[7][1] = lambdaTotal8;
422 lambda[7][2] = lambdaTotal8;
423
424 Scalar potentialDiffW0 = 0;
425 Scalar potentialDiffW1 = 0;
426 Scalar potentialDiffW2 = 0;
427 Scalar potentialDiffW3 = 0;
428 Scalar potentialDiffW4 = 0;
429 Scalar potentialDiffW5 = 0;
430 Scalar potentialDiffW6 = 0;
431 Scalar potentialDiffW7 = 0;
432 Scalar potentialDiffW8 = 0;
433 Scalar potentialDiffW9 = 0;
434 Scalar potentialDiffW10 = 0;
435 Scalar potentialDiffW11 = 0;
436
437 Scalar potentialDiffNw0 = 0;
438 Scalar potentialDiffNw1 = 0;
439 Scalar potentialDiffNw2 = 0;
440 Scalar potentialDiffNw3 = 0;
441 Scalar potentialDiffNw4 = 0;
442 Scalar potentialDiffNw5 = 0;
443 Scalar potentialDiffNw6 = 0;
444 Scalar potentialDiffNw7 = 0;
445 Scalar potentialDiffNw8 = 0;
446 Scalar potentialDiffNw9 = 0;
447 Scalar potentialDiffNw10 = 0;
448 Scalar potentialDiffNw11 = 0;
449
450 //flux vector
451 Dune::FieldVector<Scalar, 12> fluxW(0);
452 Dune::FieldVector<Scalar, 12> fluxNw(0);
453
454 DimVector Tu(0);
455 Dune::FieldVector<Scalar, 2 * dim - dim + 1> u(0);
456 TransmissibilityType T(0);
457
458 if (fIdx < 0 || fIdx == 0)
459 {
460 // calculate the flux through the subvolumeface 1 (subVolumeFaceIdx = 0)
461 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 0, 1, 2, 3,
462 4, 5);
463
464 if (caseL == 1)
465 {
466 u[0] = potW[0];
467 u[1] = potW[1];
468 u[2] = potW[2];
469 u[3] = potW[4];
470
471 T.mv(u, Tu);
472
473 fluxW[0] = Tu[0];
474 potentialDiffW0 = Tu[0];
475
476 u[0] = potNw[0];
477 u[1] = potNw[1];
478 u[2] = potNw[2];
479 u[3] = potNw[4];
480
481 T.mv(u, Tu);
482
483 fluxNw[0] = Tu[0];
484 potentialDiffNw0 = Tu[0];
485 }
486 else if (caseL == 2)
487 {
488 u[0] = potW[0];
489 u[1] = potW[1];
490 u[2] = potW[3];
491 u[3] = potW[5];
492
493 T.mv(u, Tu);
494
495 fluxW[0] = Tu[0];
496 potentialDiffW0 = Tu[0];
497
498 u[0] = potNw[0];
499 u[1] = potNw[1];
500 u[2] = potNw[3];
501 u[3] = potNw[5];
502
503 T.mv(u, Tu);
504
505 fluxNw[0] = Tu[0];
506 potentialDiffNw0 = Tu[0];
507 }
508 else if (caseL == 3)
509 {
510 u[0] = potW[0];
511 u[1] = potW[1];
512 u[2] = potW[3];
513 u[3] = potW[4];
514
515 T.mv(u, Tu);
516
517 fluxW[0] = Tu[0];
518 potentialDiffW0 = Tu[0];
519
520 u[0] = potNw[0];
521 u[1] = potNw[1];
522 u[2] = potNw[3];
523 u[3] = potNw[4];
524
525 T.mv(u, Tu);
526
527 fluxNw[0] = Tu[0];
528 potentialDiffNw0 = Tu[0];
529 }
530 else
531 {
532 u[0] = potW[0];
533 u[1] = potW[1];
534 u[2] = potW[2];
535 u[3] = potW[5];
536
537 T.mv(u, Tu);
538
539 fluxW[0] = Tu[0];
540 potentialDiffW0 = Tu[0];
541
542 u[0] = potNw[0];
543 u[1] = potNw[1];
544 u[2] = potNw[2];
545 u[3] = potNw[5];
546
547 T.mv(u, Tu);
548
549 fluxNw[0] = Tu[0];
550 potentialDiffNw0 = Tu[0];
551 }
552 }
553
554 if (fIdx < 0 || fIdx == 1)
555 {
556 // calculate the flux through the subvolumeface 2 (subVolumeFaceIdx = 1)
557 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 1, 3, 0, 2, 5,
558 7);
559
560 if (caseL == 1)
561 {
562 u[0] = potW[1];
563 u[1] = potW[3];
564 u[2] = potW[0];
565 u[3] = potW[5];
566
567 T.mv(u, Tu);
568
569 fluxW[1] = Tu[0];
570 potentialDiffW1 = Tu[0];
571
572 u[0] = potNw[1];
573 u[1] = potNw[3];
574 u[2] = potNw[0];
575 u[3] = potNw[5];
576
577 T.mv(u, Tu);
578
579 fluxNw[1] = Tu[0];
580 potentialDiffNw1 = Tu[0];
581 }
582 else if (caseL == 2)
583 {
584 u[0] = potW[1];
585 u[1] = potW[3];
586 u[2] = potW[2];
587 u[3] = potW[7];
588
589 T.mv(u, Tu);
590
591 fluxW[1] = Tu[0];
592 potentialDiffW1 = Tu[0];
593
594 u[0] = potNw[1];
595 u[1] = potNw[3];
596 u[2] = potNw[2];
597 u[3] = potNw[7];
598
599 T.mv(u, Tu);
600
601 fluxNw[1] = Tu[0];
602 potentialDiffNw1 = Tu[0];
603 }
604 else if (caseL == 3)
605 {
606 u[0] = potW[1];
607 u[1] = potW[3];
608 u[2] = potW[2];
609 u[3] = potW[5];
610
611 T.mv(u, Tu);
612
613 fluxW[1] = Tu[0];
614 potentialDiffW1 = Tu[0];
615
616 u[0] = potNw[1];
617 u[1] = potNw[3];
618 u[2] = potNw[2];
619 u[3] = potNw[5];
620
621 T.mv(u, Tu);
622
623 fluxNw[1] = Tu[0];
624 potentialDiffNw1 = Tu[0];
625 }
626 else
627 {
628 u[0] = potW[1];
629 u[1] = potW[3];
630 u[2] = potW[0];
631 u[3] = potW[7];
632
633 T.mv(u, Tu);
634
635 fluxW[1] = Tu[0];
636 potentialDiffW1 = Tu[0];
637
638 u[0] = potNw[1];
639 u[1] = potNw[3];
640 u[2] = potNw[0];
641 u[3] = potNw[7];
642
643 T.mv(u, Tu);
644
645 fluxNw[1] = Tu[0];
646 potentialDiffNw1 = Tu[0];
647 }
648 }
649
650 if (fIdx < 0 || fIdx == 2)
651 {
652 // calculate the flux through the subvolumeface 3 (subVolumeFaceIdx = 2)
653 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 3, 2, 1, 0, 7,
654 6);
655
656 if (caseL == 1)
657 {
658 u[0] = potW[3];
659 u[1] = potW[2];
660 u[2] = potW[1];
661 u[3] = potW[7];
662
663 T.mv(u, Tu);
664
665 fluxW[2] = Tu[0];
666 potentialDiffW2 = Tu[0];
667
668 u[0] = potNw[3];
669 u[1] = potNw[2];
670 u[2] = potNw[1];
671 u[3] = potNw[7];
672
673 T.mv(u, Tu);
674
675 fluxNw[2] = Tu[0];
676 potentialDiffNw2 = Tu[0];
677 }
678 else if (caseL == 2)
679 {
680 u[0] = potW[3];
681 u[1] = potW[2];
682 u[2] = potW[0];
683 u[3] = potW[6];
684
685 T.mv(u, Tu);
686
687 fluxW[2] = Tu[0];
688 potentialDiffW2 = Tu[0];
689
690 u[0] = potNw[3];
691 u[1] = potNw[2];
692 u[2] = potNw[0];
693 u[3] = potNw[6];
694
695 T.mv(u, Tu);
696
697 fluxNw[2] = Tu[0];
698 potentialDiffNw2 = Tu[0];
699 }
700 else if (caseL == 3)
701 {
702 u[0] = potW[3];
703 u[1] = potW[2];
704 u[2] = potW[0];
705 u[3] = potW[7];
706
707 T.mv(u, Tu);
708
709 fluxW[2] = Tu[0];
710 potentialDiffW2 = Tu[0];
711
712 u[0] = potNw[3];
713 u[1] = potNw[2];
714 u[2] = potNw[0];
715 u[3] = potNw[7];
716
717 T.mv(u, Tu);
718
719 fluxNw[2] = Tu[0];
720 potentialDiffNw2 = Tu[0];
721 }
722 else
723 {
724 u[0] = potW[3];
725 u[1] = potW[2];
726 u[2] = potW[1];
727 u[3] = potW[6];
728
729 T.mv(u, Tu);
730
731 fluxW[2] = Tu[0];
732 potentialDiffW2 = Tu[0];
733
734 u[0] = potNw[3];
735 u[1] = potNw[2];
736 u[2] = potNw[1];
737 u[3] = potNw[6];
738
739 T.mv(u, Tu);
740
741 fluxNw[2] = Tu[0];
742 potentialDiffNw2 = Tu[0];
743 }
744 }
745
746 if (fIdx < 0 || fIdx == 3)
747 {
748 // calculate the flux through the subvolumeface 4 (subVolumeFaceIdx = 3)
749 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 2, 0, 3, 1, 6,
750 4);
751
752 if (caseL == 1)
753 {
754 u[0] = potW[2];
755 u[1] = potW[0];
756 u[2] = potW[3];
757 u[3] = potW[6];
758
759 T.mv(u, Tu);
760
761 fluxW[3] = Tu[0];
762 potentialDiffW3 = Tu[0];
763
764 u[0] = potNw[2];
765 u[1] = potNw[0];
766 u[2] = potNw[3];
767 u[3] = potNw[6];
768
769 T.mv(u, Tu);
770
771 fluxNw[3] = Tu[0];
772 potentialDiffNw3 = Tu[0];
773 }
774 else if (caseL == 2)
775 {
776 u[0] = potW[2];
777 u[1] = potW[0];
778 u[2] = potW[1];
779 u[3] = potW[4];
780
781 T.mv(u, Tu);
782
783 fluxW[3] = Tu[0];
784 potentialDiffW3 = Tu[0];
785
786 u[0] = potNw[2];
787 u[1] = potNw[0];
788 u[2] = potNw[1];
789 u[3] = potNw[4];
790
791 T.mv(u, Tu);
792
793 fluxNw[3] = Tu[0];
794 potentialDiffNw3 = Tu[0];
795 }
796 else if (caseL == 3)
797 {
798 u[0] = potW[2];
799 u[1] = potW[0];
800 u[2] = potW[1];
801 u[3] = potW[6];
802
803 T.mv(u, Tu);
804
805 fluxW[3] = Tu[0];
806 potentialDiffW3 = Tu[0];
807
808 u[0] = potNw[2];
809 u[1] = potNw[0];
810 u[2] = potNw[1];
811 u[3] = potNw[6];
812
813 T.mv(u, Tu);
814
815 fluxNw[3] = Tu[0];
816 potentialDiffNw3 = Tu[0];
817 }
818 else
819 {
820 u[0] = potW[2];
821 u[1] = potW[0];
822 u[2] = potW[3];
823 u[3] = potW[4];
824
825 T.mv(u, Tu);
826
827 fluxW[3] = Tu[0];
828 potentialDiffW3 = Tu[0];
829
830 u[0] = potNw[2];
831 u[1] = potNw[0];
832 u[2] = potNw[3];
833 u[3] = potNw[4];
834
835 T.mv(u, Tu);
836
837 fluxNw[3] = Tu[0];
838 potentialDiffNw3 = Tu[0];
839 }
840 }
841
842 if (fIdx < 0 || fIdx == 4)
843 {
844 // calculate the flux through the subvolumeface 5 (subVolumeFaceIdx = 4)
845 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 5, 4, 7, 6, 1,
846 0);
847
848 if (caseL == 1)
849 {
850 u[0] = potW[5];
851 u[1] = potW[4];
852 u[2] = potW[7];
853 u[3] = potW[1];
854
855 T.mv(u, Tu);
856
857 fluxW[4] = Tu[0];
858 potentialDiffW4 = Tu[0];
859
860 u[0] = potNw[5];
861 u[1] = potNw[4];
862 u[2] = potNw[7];
863 u[3] = potNw[1];
864
865 T.mv(u, Tu);
866
867 fluxNw[4] = Tu[0];
868 potentialDiffNw4 = Tu[0];
869 }
870 else if (caseL == 2)
871 {
872 u[0] = potW[5];
873 u[1] = potW[4];
874 u[2] = potW[6];
875 u[3] = potW[0];
876
877 T.mv(u, Tu);
878
879 fluxW[4] = Tu[0];
880 potentialDiffW4 = Tu[0];
881
882 u[0] = potNw[5];
883 u[1] = potNw[4];
884 u[2] = potNw[6];
885 u[3] = potNw[0];
886
887 T.mv(u, Tu);
888
889 fluxNw[4] = Tu[0];
890 potentialDiffNw4 = Tu[0];
891 }
892 else if (caseL == 3)
893 {
894 u[0] = potW[5];
895 u[1] = potW[4];
896 u[2] = potW[6];
897 u[3] = potW[1];
898
899 T.mv(u, Tu);
900
901 fluxW[4] = Tu[0];
902 potentialDiffW4 = Tu[0];
903
904 u[0] = potNw[5];
905 u[1] = potNw[4];
906 u[2] = potNw[6];
907 u[3] = potNw[1];
908
909 T.mv(u, Tu);
910
911 fluxNw[4] = Tu[0];
912 potentialDiffNw4 = Tu[0];
913 }
914 else
915 {
916 u[0] = potW[5];
917 u[1] = potW[4];
918 u[2] = potW[7];
919 u[3] = potW[0];
920
921 T.mv(u, Tu);
922
923 fluxW[4] = Tu[0];
924 potentialDiffW4 = Tu[0];
925
926 u[0] = potNw[5];
927 u[1] = potNw[4];
928 u[2] = potNw[7];
929 u[3] = potNw[0];
930
931 T.mv(u, Tu);
932
933 fluxNw[4] = Tu[0];
934 potentialDiffNw4 = Tu[0];
935 }
936 }
937
938 if (fIdx < 0 || fIdx == 5)
939 {
940 // calculate the flux through the subvolumeface 6 (subVolumeFaceIdx = 5)
941 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 7, 5, 6, 4, 3,
942 1);
943
944 if (caseL == 1)
945 {
946 u[0] = potW[7];
947 u[1] = potW[5];
948 u[2] = potW[6];
949 u[3] = potW[3];
950
951 T.mv(u, Tu);
952
953 fluxW[5] = Tu[0];
954 potentialDiffW5 = Tu[0];
955
956 u[0] = potNw[7];
957 u[1] = potNw[5];
958 u[2] = potNw[6];
959 u[3] = potNw[3];
960
961 T.mv(u, Tu);
962
963 fluxNw[5] = Tu[0];
964 potentialDiffNw5 = Tu[0];
965 }
966 else if (caseL == 2)
967 {
968 u[0] = potW[7];
969 u[1] = potW[5];
970 u[2] = potW[4];
971 u[3] = potW[1];
972
973 T.mv(u, Tu);
974
975 fluxW[5] = Tu[0];
976 potentialDiffW5 = Tu[0];
977
978 u[0] = potNw[7];
979 u[1] = potNw[5];
980 u[2] = potNw[4];
981 u[3] = potNw[1];
982
983 T.mv(u, Tu);
984
985 fluxNw[5] = Tu[0];
986 potentialDiffNw5 = Tu[0];
987 }
988 else if (caseL == 3)
989 {
990 u[0] = potW[7];
991 u[1] = potW[5];
992 u[2] = potW[4];
993 u[3] = potW[3];
994
995 T.mv(u, Tu);
996
997 fluxW[5] = Tu[0];
998 potentialDiffW5 = Tu[0];
999
1000 u[0] = potNw[7];
1001 u[1] = potNw[5];
1002 u[2] = potNw[4];
1003 u[3] = potNw[3];
1004
1005 T.mv(u, Tu);
1006
1007 fluxNw[5] = Tu[0];
1008 potentialDiffNw5 = Tu[0];
1009 }
1010 else
1011 {
1012 u[0] = potW[7];
1013 u[1] = potW[5];
1014 u[2] = potW[6];
1015 u[3] = potW[1];
1016
1017 T.mv(u, Tu);
1018
1019 fluxW[5] = Tu[0];
1020 potentialDiffW5 = Tu[0];
1021
1022 u[0] = potNw[7];
1023 u[1] = potNw[5];
1024 u[2] = potNw[6];
1025 u[3] = potNw[1];
1026
1027 T.mv(u, Tu);
1028
1029 fluxNw[5] = Tu[0];
1030 potentialDiffNw5 = Tu[0];
1031 }
1032 }
1033
1034 if (fIdx < 0 || fIdx == 6)
1035 {
1036 // calculate the flux through the subvolumeface 7 (subVolumeFaceIdx = 6)
1037 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 6, 7, 4, 5, 2,
1038 3);
1039
1040 if (caseL == 1)
1041 {
1042 u[0] = potW[6];
1043 u[1] = potW[7];
1044 u[2] = potW[4];
1045 u[3] = potW[2];
1046
1047 T.mv(u, Tu);
1048
1049 fluxW[6] = Tu[0];
1050 potentialDiffW6 = Tu[0];
1051
1052 u[0] = potNw[6];
1053 u[1] = potNw[7];
1054 u[2] = potNw[4];
1055 u[3] = potNw[2];
1056
1057 T.mv(u, Tu);
1058
1059 fluxNw[6] = Tu[0];
1060 potentialDiffNw6 = Tu[0];
1061 }
1062 else if (caseL == 2)
1063 {
1064 u[0] = potW[6];
1065 u[1] = potW[7];
1066 u[2] = potW[5];
1067 u[3] = potW[3];
1068
1069 T.mv(u, Tu);
1070
1071 fluxW[6] = Tu[0];
1072 potentialDiffW6 = Tu[0];
1073
1074 u[0] = potNw[6];
1075 u[1] = potNw[7];
1076 u[2] = potNw[5];
1077 u[3] = potNw[3];
1078
1079 T.mv(u, Tu);
1080
1081 fluxNw[6] = Tu[0];
1082 potentialDiffNw6 = Tu[0];
1083 }
1084 else if (caseL == 3)
1085 {
1086 u[0] = potW[6];
1087 u[1] = potW[7];
1088 u[2] = potW[5];
1089 u[3] = potW[2];
1090
1091 T.mv(u, Tu);
1092
1093 fluxW[6] = Tu[0];
1094 potentialDiffW6 = Tu[0];
1095
1096 u[0] = potNw[6];
1097 u[1] = potNw[7];
1098 u[2] = potNw[5];
1099 u[3] = potNw[2];
1100
1101 T.mv(u, Tu);
1102
1103 fluxNw[6] = Tu[0];
1104 potentialDiffNw6 = Tu[0];
1105 }
1106 else
1107 {
1108 u[0] = potW[6];
1109 u[1] = potW[7];
1110 u[2] = potW[4];
1111 u[3] = potW[3];
1112
1113 T.mv(u, Tu);
1114
1115 fluxW[6] = Tu[0];
1116 potentialDiffW6 = Tu[0];
1117
1118 u[0] = potNw[6];
1119 u[1] = potNw[7];
1120 u[2] = potNw[4];
1121 u[3] = potNw[3];
1122
1123 T.mv(u, Tu);
1124
1125 fluxNw[6] = Tu[0];
1126 potentialDiffNw6 = Tu[0];
1127 }
1128 }
1129
1130 if (fIdx < 0 || fIdx == 7)
1131 {
1132 // calculate the flux through the subvolumeface 8 (subVolumeFaceIdx = 7)
1133 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 4, 6, 5, 7, 0,
1134 2);
1135
1136 if (caseL == 1)
1137 {
1138 u[0] = potW[4];
1139 u[1] = potW[6];
1140 u[2] = potW[5];
1141 u[3] = potW[0];
1142
1143 T.mv(u, Tu);
1144
1145 fluxW[7] = Tu[0];
1146 potentialDiffW7 = Tu[0];
1147
1148 u[0] = potNw[4];
1149 u[1] = potNw[6];
1150 u[2] = potNw[5];
1151 u[3] = potNw[0];
1152
1153 T.mv(u, Tu);
1154
1155 fluxNw[7] = Tu[0];
1156 potentialDiffNw7 = Tu[0];
1157 }
1158 else if (caseL == 2)
1159 {
1160 u[0] = potW[4];
1161 u[1] = potW[6];
1162 u[2] = potW[7];
1163 u[3] = potW[2];
1164
1165 T.mv(u, Tu);
1166
1167 fluxW[7] = Tu[0];
1168 potentialDiffW7 = Tu[0];
1169
1170 u[0] = potNw[4];
1171 u[1] = potNw[6];
1172 u[2] = potNw[7];
1173 u[3] = potNw[2];
1174
1175 T.mv(u, Tu);
1176
1177 fluxNw[7] = Tu[0];
1178 potentialDiffNw7 = Tu[0];
1179 }
1180 else if (caseL == 3)
1181 {
1182 u[0] = potW[4];
1183 u[1] = potW[6];
1184 u[2] = potW[7];
1185 u[3] = potW[0];
1186
1187 T.mv(u, Tu);
1188
1189 fluxW[7] = Tu[0];
1190 potentialDiffW7 = Tu[0];
1191
1192 u[0] = potNw[4];
1193 u[1] = potNw[6];
1194 u[2] = potNw[7];
1195 u[3] = potNw[0];
1196
1197 T.mv(u, Tu);
1198
1199 fluxNw[7] = Tu[0];
1200 potentialDiffNw7 = Tu[0];
1201 }
1202 else
1203 {
1204 u[0] = potW[4];
1205 u[1] = potW[6];
1206 u[2] = potW[5];
1207 u[3] = potW[2];
1208
1209 T.mv(u, Tu);
1210
1211 fluxW[7] = Tu[0];
1212 potentialDiffW7 = Tu[0];
1213
1214 u[0] = potNw[4];
1215 u[1] = potNw[6];
1216 u[2] = potNw[5];
1217 u[3] = potNw[2];
1218
1219 T.mv(u, Tu);
1220
1221 fluxNw[7] = Tu[0];
1222 potentialDiffNw7 = Tu[0];
1223 }
1224 }
1225
1226 if (fIdx < 0 || fIdx == 8)
1227 {
1228 // calculate the flux through the subvolumeface 9 (subVolumeFaceIdx = 8)
1229 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 4, 0, 6, 2, 5,
1230 1);
1231
1232 if (caseL == 1)
1233 {
1234 u[0] = potW[4];
1235 u[1] = potW[0];
1236 u[2] = potW[6];
1237 u[3] = potW[5];
1238
1239 T.mv(u, Tu);
1240
1241 fluxW[8] = Tu[0];
1242 potentialDiffW8 = Tu[0];
1243
1244 u[0] = potNw[4];
1245 u[1] = potNw[0];
1246 u[2] = potNw[6];
1247 u[3] = potNw[5];
1248
1249 T.mv(u, Tu);
1250
1251 fluxNw[8] = Tu[0];
1252 potentialDiffNw8 = Tu[0];
1253 }
1254 else if (caseL == 2)
1255 {
1256 u[0] = potW[4];
1257 u[1] = potW[0];
1258 u[2] = potW[2];
1259 u[3] = potW[1];
1260
1261 T.mv(u, Tu);
1262
1263 fluxW[8] = Tu[0];
1264 potentialDiffW8 = Tu[0];
1265
1266 u[0] = potNw[4];
1267 u[1] = potNw[0];
1268 u[2] = potNw[2];
1269 u[3] = potNw[1];
1270
1271 T.mv(u, Tu);
1272
1273 fluxNw[8] = Tu[0];
1274 potentialDiffNw8 = Tu[0];
1275 }
1276 else if (caseL == 3)
1277 {
1278 u[0] = potW[4];
1279 u[1] = potW[0];
1280 u[2] = potW[2];
1281 u[3] = potW[5];
1282
1283 T.mv(u, Tu);
1284
1285 fluxW[8] = Tu[0];
1286 potentialDiffW8 = Tu[0];
1287
1288 u[0] = potNw[4];
1289 u[1] = potNw[0];
1290 u[2] = potNw[2];
1291 u[3] = potNw[5];
1292
1293 T.mv(u, Tu);
1294
1295 fluxNw[8] = Tu[0];
1296 potentialDiffNw8 = Tu[0];
1297 }
1298 else
1299 {
1300 u[0] = potW[4];
1301 u[1] = potW[0];
1302 u[2] = potW[6];
1303 u[3] = potW[1];
1304
1305 T.mv(u, Tu);
1306
1307 fluxW[8] = Tu[0];
1308 potentialDiffW8 = Tu[0];
1309
1310 u[0] = potNw[4];
1311 u[1] = potNw[0];
1312 u[2] = potNw[6];
1313 u[3] = potNw[1];
1314
1315 T.mv(u, Tu);
1316
1317 fluxNw[8] = Tu[0];
1318 potentialDiffNw8 = Tu[0];
1319 }
1320 }
1321
1322 if (fIdx < 0 || fIdx == 9)
1323 {
1324 // calculate the flux through the subvolumeface 10 (subVolumeFaceIdx = 9)
1325 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 1, 5, 3, 7, 0,
1326 4);
1327
1328 if (caseL == 1)
1329 {
1330 u[0] = potW[1];
1331 u[1] = potW[5];
1332 u[2] = potW[3];
1333 u[3] = potW[0];
1334
1335 T.mv(u, Tu);
1336
1337 fluxW[9] = Tu[0];
1338 potentialDiffW9 = Tu[0];
1339
1340 u[0] = potNw[1];
1341 u[1] = potNw[5];
1342 u[2] = potNw[3];
1343 u[3] = potNw[0];
1344
1345 T.mv(u, Tu);
1346
1347 fluxNw[9] = Tu[0];
1348 potentialDiffNw9 = Tu[0];
1349 }
1350 else if (caseL == 2)
1351 {
1352 u[0] = potW[1];
1353 u[1] = potW[5];
1354 u[2] = potW[7];
1355 u[3] = potW[4];
1356
1357 T.mv(u, Tu);
1358
1359 fluxW[9] = Tu[0];
1360 potentialDiffW9 = Tu[0];
1361
1362 u[0] = potNw[1];
1363 u[1] = potNw[5];
1364 u[2] = potNw[7];
1365 u[3] = potNw[4];
1366
1367 T.mv(u, Tu);
1368
1369 fluxNw[9] = Tu[0];
1370 potentialDiffNw9 = Tu[0];
1371 }
1372 else if (caseL == 3)
1373 {
1374 u[0] = potW[1];
1375 u[1] = potW[5];
1376 u[2] = potW[7];
1377 u[3] = potW[0];
1378
1379 T.mv(u, Tu);
1380
1381 fluxW[9] = Tu[0];
1382 potentialDiffW9 = Tu[0];
1383
1384 u[0] = potNw[1];
1385 u[1] = potNw[5];
1386 u[2] = potNw[7];
1387 u[3] = potNw[0];
1388
1389 T.mv(u, Tu);
1390
1391 fluxNw[9] = Tu[0];
1392 potentialDiffNw9 = Tu[0];
1393 }
1394 else
1395 {
1396 u[0] = potW[1];
1397 u[1] = potW[5];
1398 u[2] = potW[3];
1399 u[3] = potW[4];
1400
1401 T.mv(u, Tu);
1402
1403 fluxW[9] = Tu[0];
1404 potentialDiffW9 = Tu[0];
1405
1406 u[0] = potNw[1];
1407 u[1] = potNw[5];
1408 u[2] = potNw[3];
1409 u[3] = potNw[4];
1410
1411 T.mv(u, Tu);
1412
1413 fluxNw[9] = Tu[0];
1414 potentialDiffNw9 = Tu[0];
1415 }
1416 }
1417
1418 if (fIdx < 0 || fIdx == 10)
1419 {
1420 // calculate the flux through the subvolumeface 11 (subVolumeFaceIdx = 10)
1421 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 7, 3, 5, 1, 6,
1422 2);
1423
1424 if (caseL == 1)
1425 {
1426 u[0] = potW[7];
1427 u[1] = potW[3];
1428 u[2] = potW[5];
1429 u[3] = potW[6];
1430
1431 T.mv(u, Tu);
1432
1433 fluxW[10] = Tu[0];
1434 potentialDiffW10 = Tu[0];
1435
1436 u[0] = potNw[7];
1437 u[1] = potNw[3];
1438 u[2] = potNw[5];
1439 u[3] = potNw[6];
1440
1441 T.mv(u, Tu);
1442
1443 fluxNw[10] = Tu[0];
1444 potentialDiffNw10 = Tu[0];
1445 }
1446 else if (caseL == 2)
1447 {
1448 u[0] = potW[7];
1449 u[1] = potW[3];
1450 u[2] = potW[1];
1451 u[3] = potW[2];
1452
1453 T.mv(u, Tu);
1454
1455 fluxW[10] = Tu[0];
1456 potentialDiffW10 = Tu[0];
1457
1458 u[0] = potNw[7];
1459 u[1] = potNw[3];
1460 u[2] = potNw[1];
1461 u[3] = potNw[2];
1462
1463 T.mv(u, Tu);
1464
1465 fluxNw[10] = Tu[0];
1466 potentialDiffNw10 = Tu[0];
1467 }
1468 else if (caseL == 3)
1469 {
1470 u[0] = potW[7];
1471 u[1] = potW[3];
1472 u[2] = potW[1];
1473 u[3] = potW[6];
1474
1475 T.mv(u, Tu);
1476
1477 fluxW[10] = Tu[0];
1478 potentialDiffW10 = Tu[0];
1479
1480 u[0] = potNw[7];
1481 u[1] = potNw[3];
1482 u[2] = potNw[1];
1483 u[3] = potNw[6];
1484
1485 T.mv(u, Tu);
1486
1487 fluxNw[10] = Tu[0];
1488 potentialDiffNw10 = Tu[0];
1489 }
1490 else
1491 {
1492 u[0] = potW[7];
1493 u[1] = potW[3];
1494 u[2] = potW[5];
1495 u[3] = potW[2];
1496
1497 T.mv(u, Tu);
1498
1499 fluxW[10] = Tu[0];
1500 potentialDiffW10 = Tu[0];
1501
1502 u[0] = potNw[7];
1503 u[1] = potNw[3];
1504 u[2] = potNw[5];
1505 u[3] = potNw[2];
1506
1507 T.mv(u, Tu);
1508
1509 fluxNw[10] = Tu[0];
1510 potentialDiffNw10 = Tu[0];
1511 }
1512 }
1513
1514 if (fIdx < 0 || fIdx == 11)
1515 {
1516 // calculate the flux through the subvolumeface 12 (subVolumeFaceIdx = 11)
1517 int caseL = transmissibilityCalculator.transmissibility(T, interactionVolume, lambda, 2, 6, 0, 4, 3,
1518 7);
1519
1520 if (caseL == 1)
1521 {
1522 u[0] = potW[2];
1523 u[1] = potW[6];
1524 u[2] = potW[0];
1525 u[3] = potW[3];
1526
1527 T.mv(u, Tu);
1528
1529 fluxW[11] = Tu[0];
1530 potentialDiffW11 = Tu[0];
1531
1532 u[0] = potNw[2];
1533 u[1] = potNw[6];
1534 u[2] = potNw[0];
1535 u[3] = potNw[3];
1536
1537 T.mv(u, Tu);
1538
1539 fluxNw[11] = Tu[0];
1540 potentialDiffNw11 = Tu[0];
1541 }
1542 else if (caseL == 2)
1543 {
1544 u[0] = potW[2];
1545 u[1] = potW[6];
1546 u[2] = potW[4];
1547 u[3] = potW[7];
1548
1549 T.mv(u, Tu);
1550
1551 fluxW[11] = Tu[0];
1552 potentialDiffW11 = Tu[0];
1553
1554 u[0] = potNw[2];
1555 u[1] = potNw[6];
1556 u[2] = potNw[4];
1557 u[3] = potNw[7];
1558
1559 T.mv(u, Tu);
1560
1561 fluxNw[11] = Tu[0];
1562 potentialDiffNw11 = Tu[0];
1563 }
1564 else if (caseL == 3)
1565 {
1566 u[0] = potW[2];
1567 u[1] = potW[6];
1568 u[2] = potW[4];
1569 u[3] = potW[3];
1570
1571 T.mv(u, Tu);
1572
1573 fluxW[11] = Tu[0];
1574 potentialDiffW11 = Tu[0];
1575
1576 u[0] = potNw[2];
1577 u[1] = potNw[6];
1578 u[2] = potNw[4];
1579 u[3] = potNw[3];
1580
1581 T.mv(u, Tu);
1582
1583 fluxNw[11] = Tu[0];
1584 potentialDiffNw11 = Tu[0];
1585 }
1586 else
1587 {
1588 u[0] = potW[2];
1589 u[1] = potW[6];
1590 u[2] = potW[0];
1591 u[3] = potW[7];
1592
1593 T.mv(u, Tu);
1594
1595 fluxW[11] = Tu[0];
1596 potentialDiffW11 = Tu[0];
1597
1598 u[0] = potNw[2];
1599 u[1] = potNw[6];
1600 u[2] = potNw[0];
1601 u[3] = potNw[7];
1602
1603 T.mv(u, Tu);
1604
1605 fluxNw[11] = Tu[0];
1606 potentialDiffNw11 = Tu[0];
1607 }
1608 }
1609
1610 //store potentials for further calculations (saturation, ...)
1611 cellData1.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(0, 0), potentialDiffW0);
1612 cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 0), potentialDiffNw0);
1613 cellData1.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(0, 1), -potentialDiffW3);
1614 cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 1), -potentialDiffNw3);
1615 cellData1.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(0, 2), -potentialDiffW8);
1616 cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 2), -potentialDiffNw8);
1617
1618 cellData2.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(1, 0), potentialDiffW1);
1619 cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 0), potentialDiffNw1);
1620 cellData2.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(1, 1), -potentialDiffW0);
1621 cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 1), -potentialDiffNw0);
1622 cellData2.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(1, 2), potentialDiffW9);
1623 cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 2), potentialDiffNw9);
1624
1625 cellData3.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(2, 0), potentialDiffW3);
1626 cellData3.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(2, 0), potentialDiffNw3);
1627 cellData3.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(2, 1), -potentialDiffW2);
1628 cellData3.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(2, 1), -potentialDiffNw2);
1629 cellData3.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(2, 2), potentialDiffW11);
1630 cellData3.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(2, 2), potentialDiffNw11);
1631
1632 cellData4.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(3, 0), potentialDiffW2);
1633 cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 0), potentialDiffNw2);
1634 cellData4.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(3, 1), -potentialDiffW1);
1635 cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 1), -potentialDiffNw1);
1636 cellData4.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(3, 2), -potentialDiffW10);
1637 cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 2), -potentialDiffNw10);
1638
1639 cellData5.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(4, 0), potentialDiffW8);
1640 cellData5.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(4, 0), potentialDiffNw8);
1641 cellData5.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(4, 1), -potentialDiffW4);
1642 cellData5.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(4, 1), -potentialDiffNw4);
1643 cellData5.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(4, 2), potentialDiffW7);
1644 cellData5.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(4, 2), potentialDiffNw7);
1645
1646 cellData6.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(5, 0), -potentialDiffW9);
1647 cellData6.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(5, 0), -potentialDiffNw9);
1648 cellData6.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(5, 1), -potentialDiffW5);
1649 cellData6.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(5, 1), -potentialDiffNw5);
1650 cellData6.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(5, 2), potentialDiffW4);
1651 cellData6.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(5, 2), potentialDiffNw4);
1652
1653 cellData7.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(6, 0), -potentialDiffW11);
1654 cellData7.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(6, 0), -potentialDiffNw11);
1655 cellData7.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(6, 1), -potentialDiffW7);
1656 cellData7.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(6, 1), -potentialDiffNw7);
1657 cellData7.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(6, 2), potentialDiffW6);
1658 cellData7.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(6, 2), potentialDiffNw6);
1659
1660 cellData8.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(7, 0), potentialDiffW10);
1661 cellData8.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(7, 0), potentialDiffNw10);
1662 cellData8.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(7, 1), -potentialDiffW6);
1663 cellData8.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(7, 1), -potentialDiffNw6);
1664 cellData8.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(7, 2), potentialDiffW5);
1665 cellData8.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(7, 2), potentialDiffNw5);
1666
1667 // compute mobilities of subvolumeface 1 (subVolumeFaceIdx = 0)
1668 Dune::FieldVector<Scalar, numPhases> lambda0Upw(0.0);
1669 lambda0Upw[wPhaseIdx] = (potentialDiffW0 >= 0) ? lambda1[wPhaseIdx] : lambda2[wPhaseIdx];
1670 lambda0Upw[nPhaseIdx] = (potentialDiffNw0 >= 0) ? lambda1[nPhaseIdx] : lambda2[nPhaseIdx];
1671
1672 // compute mobilities of subvolumeface 2 (subVolumeFaceIdx = 1)
1673 Dune::FieldVector<Scalar, numPhases> lambda1Upw(0.0);
1674 lambda1Upw[wPhaseIdx] = (potentialDiffW1 >= 0) ? lambda2[wPhaseIdx] : lambda4[wPhaseIdx];
1675 lambda1Upw[nPhaseIdx] = (potentialDiffNw1 >= 0) ? lambda2[nPhaseIdx] : lambda4[nPhaseIdx];
1676
1677 // compute mobilities of subvolumeface 3 (subVolumeFaceIdx = 2)
1678 Dune::FieldVector<Scalar, numPhases> lambda2Upw(0.0);
1679 lambda2Upw[wPhaseIdx] = (potentialDiffW2 >= 0) ? lambda4[wPhaseIdx] : lambda3[wPhaseIdx];
1680 lambda2Upw[nPhaseIdx] = (potentialDiffNw2 >= 0) ? lambda4[nPhaseIdx] : lambda3[nPhaseIdx];
1681
1682 // compute mobilities of subvolumeface 4 (subVolumeFaceIdx = 3)
1683 Dune::FieldVector<Scalar, numPhases> lambda3Upw(0.0);
1684 lambda3Upw[wPhaseIdx] = (potentialDiffW3 >= 0) ? lambda3[wPhaseIdx] : lambda1[wPhaseIdx];
1685 lambda3Upw[nPhaseIdx] = (potentialDiffNw3 >= 0) ? lambda3[nPhaseIdx] : lambda1[nPhaseIdx];
1686
1687 // compute mobilities of subvolumeface 5 (subVolumeFaceIdx = 4)
1688 Dune::FieldVector<Scalar, numPhases> lambda4Upw(0.0);
1689 lambda4Upw[wPhaseIdx] = (potentialDiffW4 >= 0) ? lambda6[wPhaseIdx] : lambda5[wPhaseIdx];
1690 lambda4Upw[nPhaseIdx] = (potentialDiffNw4 >= 0) ? lambda6[nPhaseIdx] : lambda5[nPhaseIdx];
1691
1692 // compute mobilities of subvolumeface 6 (subVolumeFaceIdx = 5)
1693 Dune::FieldVector<Scalar, numPhases> lambda5Upw(0.0);
1694 lambda5Upw[wPhaseIdx] = (potentialDiffW5 >= 0) ? lambda8[wPhaseIdx] : lambda6[wPhaseIdx];
1695 lambda5Upw[nPhaseIdx] = (potentialDiffNw5 >= 0) ? lambda8[nPhaseIdx] : lambda6[nPhaseIdx];
1696
1697 // compute mobilities of subvolumeface 7 (subVolumeFaceIdx = 6)
1698 Dune::FieldVector<Scalar, numPhases> lambda6Upw(0.0);
1699 lambda6Upw[wPhaseIdx] = (potentialDiffW6 >= 0) ? lambda7[wPhaseIdx] : lambda8[wPhaseIdx];
1700 lambda6Upw[nPhaseIdx] = (potentialDiffNw6 >= 0) ? lambda7[nPhaseIdx] : lambda8[nPhaseIdx];
1701
1702 // compute mobilities of subvolumeface 8 (subVolumeFaceIdx = 7)
1703 Dune::FieldVector<Scalar, numPhases> lambda7Upw(0.0);
1704 lambda7Upw[wPhaseIdx] = (potentialDiffW7 >= 0) ? lambda5[wPhaseIdx] : lambda7[wPhaseIdx];
1705 lambda7Upw[nPhaseIdx] = (potentialDiffNw7 >= 0) ? lambda5[nPhaseIdx] : lambda7[nPhaseIdx];
1706
1707 // compute mobilities of subvolumeface 9 (subVolumeFaceIdx = 8)
1708 Dune::FieldVector<Scalar, numPhases> lambda8Upw(0.0);
1709 lambda8Upw[wPhaseIdx] = (potentialDiffW8 >= 0) ? lambda5[wPhaseIdx] : lambda1[wPhaseIdx];
1710 lambda8Upw[nPhaseIdx] = (potentialDiffNw8 >= 0) ? lambda5[nPhaseIdx] : lambda1[nPhaseIdx];
1711
1712 // compute mobilities of subvolumeface 10 (subVolumeFaceIdx = 9)
1713 Dune::FieldVector<Scalar, numPhases> lambda9Upw(0.0);
1714 lambda9Upw[wPhaseIdx] = (potentialDiffW9 >= 0) ? lambda2[wPhaseIdx] : lambda6[wPhaseIdx];
1715 lambda9Upw[nPhaseIdx] = (potentialDiffNw9 >= 0) ? lambda2[nPhaseIdx] : lambda6[nPhaseIdx];
1716
1717 // compute mobilities of subvolumeface 11 (subVolumeFaceIdx = 10)
1718 Dune::FieldVector<Scalar, numPhases> lambda10Upw(0.0);
1719 lambda10Upw[wPhaseIdx] = (potentialDiffW10 >= 0) ? lambda8[wPhaseIdx] : lambda4[wPhaseIdx];
1720 lambda10Upw[nPhaseIdx] = (potentialDiffNw10 >= 0) ? lambda8[nPhaseIdx] : lambda4[nPhaseIdx];
1721
1722 // compute mobilities of subvolumeface 12 (subVolumeFaceIdx = 11)
1723 Dune::FieldVector<Scalar, numPhases> lambda11Upw(0.0);
1724 lambda11Upw[wPhaseIdx] = (potentialDiffW11 >= 0) ? lambda3[wPhaseIdx] : lambda7[wPhaseIdx];
1725 lambda11Upw[nPhaseIdx] = (potentialDiffNw11 >= 0) ? lambda3[nPhaseIdx] : lambda7[nPhaseIdx];
1726
1727 for (int i = 0; i < numPhases; i++)
1728 {
1729 // evaluate parts of velocity --> always take the normal for which the flux is calculated!
1730 DimVector vel12 = interactionVolume.getNormal(0, 0);
1731 DimVector vel21 = interactionVolume.getNormal(1, 1);
1732 DimVector vel24 = interactionVolume.getNormal(1, 0);
1733 DimVector vel42 = interactionVolume.getNormal(3, 1);
1734 DimVector vel43 = interactionVolume.getNormal(3, 0);
1735 DimVector vel34 = interactionVolume.getNormal(2, 1);
1736 DimVector vel31 = interactionVolume.getNormal(2, 0);
1737 DimVector vel13 = interactionVolume.getNormal(0, 1);
1738 DimVector vel65 = interactionVolume.getNormal(5, 2);
1739 DimVector vel56 = interactionVolume.getNormal(4, 1);
1740 DimVector vel86 = interactionVolume.getNormal(7, 2);
1741 DimVector vel68 = interactionVolume.getNormal(5, 1);
1742 DimVector vel78 = interactionVolume.getNormal(6, 2);
1743 DimVector vel87 = interactionVolume.getNormal(7, 1);
1744 DimVector vel57 = interactionVolume.getNormal(4, 2);
1745 DimVector vel75 = interactionVolume.getNormal(6, 1);
1746 DimVector vel51 = interactionVolume.getNormal(4, 0);
1747 DimVector vel15 = interactionVolume.getNormal(0, 2);
1748 DimVector vel26 = interactionVolume.getNormal(1, 2);
1749 DimVector vel62 = interactionVolume.getNormal(5, 0);
1750 DimVector vel84 = interactionVolume.getNormal(7, 0);
1751 DimVector vel48 = interactionVolume.getNormal(3, 2);
1752 DimVector vel37 = interactionVolume.getNormal(2, 2);
1753 DimVector vel73 = interactionVolume.getNormal(6, 0);
1754 vel21 *= -1;
1755 vel42 *= -1;
1756 vel34 *= -1;
1757 vel13 *= -1;
1758 vel56 *= -1;
1759 vel68 *= -1;
1760 vel87 *= -1;
1761 vel75 *= -1;
1762 vel15 *= -1;
1763 vel62 *= -1;
1764 vel48 *= -1;
1765 vel73 *= -1;
1766
1767 Dune::FieldVector<Scalar, 12> flux(0);
1768 switch (i)
1769 {
1770 case wPhaseIdx:
1771 {
1772 flux = fluxW;
1773 break;
1774 }
1775 case nPhaseIdx:
1776 {
1777 flux = fluxNw;
1778 break;
1779 }
1780 }
1781
1782 vel12 *= flux[0] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal1, 0, 0));
1783 vel21 *= flux[0] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal2, 1, 1));
1784 vel24 *= flux[1] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal2, 1, 0));
1785 vel42 *= flux[1] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal4, 3, 1));
1786 vel43 *= flux[2] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal4, 3, 0));
1787 vel34 *= flux[2] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal3, 2, 1));
1788 vel31 *= flux[3] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal3, 2, 0));
1789 vel13 *= flux[3] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal1, 0, 1));
1790 vel65 *= flux[4] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal6, 5, 2));
1791 vel56 *= flux[4] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal5, 4, 1));
1792 vel86 *= flux[5] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal8, 7, 2));
1793 vel68 *= flux[5] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal6, 5, 1));
1794 vel78 *= flux[6] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal7, 6, 2));
1795 vel87 *= flux[6] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal8, 7, 1));
1796 vel57 *= flux[7] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal5, 4, 2));
1797 vel75 *= flux[7] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal7, 6, 1));
1798 vel51 *= flux[8] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal5, 4, 0));
1799 vel15 *= flux[8] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal1, 0, 2));
1800 vel26 *= flux[9] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal2, 1, 2));
1801 vel62 *= flux[9] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal6, 5, 0));
1802 vel84 *= flux[10] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal8, 7, 0));
1803 vel48 *= flux[10] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal4, 3, 2));
1804 vel37 *= flux[11] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal3, 2, 2));
1805 vel73 *= flux[11] / (interactionVolumes.getRealFluxFaceArea(interactionVolume, eIdxGlobal7, 6, 0));
1806
1807
1808 Scalar lambdaT0 = lambda0Upw[wPhaseIdx] + lambda0Upw[nPhaseIdx];
1809 Scalar lambdaT1 = lambda1Upw[wPhaseIdx] + lambda1Upw[nPhaseIdx];
1810 Scalar lambdaT2 = lambda2Upw[wPhaseIdx] + lambda2Upw[nPhaseIdx];
1811 Scalar lambdaT3 = lambda3Upw[wPhaseIdx] + lambda3Upw[nPhaseIdx];
1812 Scalar lambdaT4 = lambda4Upw[wPhaseIdx] + lambda4Upw[nPhaseIdx];
1813 Scalar lambdaT5 = lambda5Upw[wPhaseIdx] + lambda5Upw[nPhaseIdx];
1814 Scalar lambdaT6 = lambda6Upw[wPhaseIdx] + lambda6Upw[nPhaseIdx];
1815 Scalar lambdaT7 = lambda7Upw[wPhaseIdx] + lambda7Upw[nPhaseIdx];
1816 Scalar lambdaT8 = lambda8Upw[wPhaseIdx] + lambda8Upw[nPhaseIdx];
1817 Scalar lambdaT9 = lambda9Upw[wPhaseIdx] + lambda9Upw[nPhaseIdx];
1818 Scalar lambdaT10 = lambda10Upw[wPhaseIdx] + lambda10Upw[nPhaseIdx];
1819 Scalar lambdaT11 = lambda11Upw[wPhaseIdx] + lambda11Upw[nPhaseIdx];
1820 Scalar fracFlow0 = (lambdaT0 > threshold_) ? lambda0Upw[i] / (lambdaT0) : 0.0;
1821 Scalar fracFlow1 = (lambdaT1 > threshold_) ? lambda1Upw[i] / (lambdaT1) : 0.0;
1822 Scalar fracFlow2 = (lambdaT2 > threshold_) ? lambda2Upw[i] / (lambdaT2) : 0.0;
1823 Scalar fracFlow3 = (lambdaT3 > threshold_) ? lambda3Upw[i] / (lambdaT3) : 0.0;
1824 Scalar fracFlow4 = (lambdaT4 > threshold_) ? lambda4Upw[i] / (lambdaT4) : 0.0;
1825 Scalar fracFlow5 = (lambdaT5 > threshold_) ? lambda5Upw[i] / (lambdaT5) : 0.0;
1826 Scalar fracFlow6 = (lambdaT6 > threshold_) ? lambda6Upw[i] / (lambdaT6) : 0.0;
1827 Scalar fracFlow7 = (lambdaT7 > threshold_) ? lambda7Upw[i] / (lambdaT7) : 0.0;
1828 Scalar fracFlow8 = (lambdaT8 > threshold_) ? lambda8Upw[i] / (lambdaT8) : 0.0;
1829 Scalar fracFlow9 = (lambdaT9 > threshold_) ? lambda9Upw[i] / (lambdaT9) : 0.0;
1830 Scalar fracFlow10 = (lambdaT10 > threshold_) ? lambda10Upw[i] / (lambdaT10) : 0.0;
1831 Scalar fracFlow11 = (lambdaT11 > threshold_) ? lambda11Upw[i] / (lambdaT11) : 0.0;
1832
1833 vel12 *= fracFlow0;
1834 vel21 *= fracFlow0;
1835 vel24 *= fracFlow1;
1836 vel42 *= fracFlow1;
1837 vel43 *= fracFlow2;
1838 vel34 *= fracFlow2;
1839 vel31 *= fracFlow3;
1840 vel13 *= fracFlow3;
1841 vel65 *= fracFlow4;
1842 vel56 *= fracFlow4;
1843 vel86 *= fracFlow5;
1844 vel68 *= fracFlow5;
1845 vel78 *= fracFlow6;
1846 vel87 *= fracFlow6;
1847 vel57 *= fracFlow7;
1848 vel75 *= fracFlow7;
1849 vel51 *= fracFlow8;
1850 vel15 *= fracFlow8;
1851 vel26 *= fracFlow9;
1852 vel62 *= fracFlow9;
1853 vel84 *= fracFlow10;
1854 vel48 *= fracFlow10;
1855 vel37 *= fracFlow11;
1856 vel73 *= fracFlow11;
1857
1858 //store velocities
1859 cellData1.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(0, 0), vel12);
1860 cellData1.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(0, 1), vel13);
1861 cellData1.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(0, 2), vel15);
1862 cellData2.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(1, 0), vel24);
1863 cellData2.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(1, 1), vel21);
1864 cellData2.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(1, 2), vel26);
1865 cellData3.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(2, 0), vel31);
1866 cellData3.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(2, 1), vel34);
1867 cellData3.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(2, 2), vel37);
1868 cellData4.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(3, 0), vel43);
1869 cellData4.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(3, 1), vel42);
1870 cellData4.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(3, 2), vel48);
1871 cellData5.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(4, 0), vel51);
1872 cellData5.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(4, 1), vel56);
1873 cellData5.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(4, 2), vel57);
1874 cellData6.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(5, 0), vel62);
1875 cellData6.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(5, 1), vel68);
1876 cellData6.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(5, 2), vel65);
1877 cellData7.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(6, 0), vel73);
1878 cellData7.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(6, 1), vel75);
1879 cellData7.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(6, 2), vel78);
1880 cellData8.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(7, 0), vel84);
1881 cellData8.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(7, 1), vel87);
1882 cellData8.fluxData().addVelocity(i, interactionVolume.getIndexOnElement(7, 2), vel86);
1883 }
1884 //set velocity marker
1885 cellData1.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(0, 0));
1886 cellData1.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(0, 1));
1887 cellData1.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(0, 2));
1888 cellData2.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(1, 0));
1889 cellData2.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(1, 1));
1890 cellData2.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(1, 2));
1891 cellData3.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(2, 0));
1892 cellData3.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(2, 1));
1893 cellData3.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(2, 2));
1894 cellData4.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(3, 0));
1895 cellData4.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(3, 1));
1896 cellData4.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(3, 2));
1897 cellData5.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(4, 0));
1898 cellData5.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(4, 1));
1899 cellData5.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(4, 2));
1900 cellData6.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(5, 0));
1901 cellData6.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(5, 1));
1902 cellData6.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(5, 2));
1903 cellData7.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(6, 0));
1904 cellData7.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(6, 1));
1905 cellData7.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(6, 2));
1906 cellData8.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(7, 0));
1907 cellData8.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(7, 1));
1908 cellData8.fluxData().setVelocityMarker(interactionVolume.getIndexOnElement(7, 2));
1909 }
1910
1920template<class TypeTag>
1921void FvMpfaL3dVelocity2p<TypeTag>::calculateBoundaryInteractionVolumeVelocity(InteractionVolume& interactionVolume, CellData& cellData, int elemIdx)
1922{
1923 auto element = interactionVolume.getSubVolumeElement(elemIdx);
1924
1925 // get global coordinate of cell centers
1926 const GlobalPosition& globalPos = element.geometry().center();
1927
1928 // permeability vector at boundary
1929 DimMatrix permeability(problem_.spatialParams().intrinsicPermeability(element));
1930
1931 //get mobilities of the phases
1932 Dune::FieldVector<Scalar, numPhases> lambda(cellData.mobility(wPhaseIdx));
1933 lambda[nPhaseIdx] = cellData.mobility(nPhaseIdx);
1934
1935 Scalar potW = cellData.potential(wPhaseIdx);
1936 Scalar potNw = cellData.potential(nPhaseIdx);
1937
1938 for (int fIdx = 0; fIdx < dim; fIdx++)
1939 {
1940 int intVolFaceIdx = interactionVolume.getFaceIndexFromSubVolume(elemIdx, fIdx);
1941
1942 if (interactionVolume.isBoundaryFace(intVolFaceIdx))
1943 {
1944 if (interactionVolume.getBoundaryType(intVolFaceIdx).isDirichlet(pressureEqIdx))
1945 {
1946 int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
1947
1948 const GlobalPosition& globalPosFace = interactionVolume.getFacePosition(elemIdx, fIdx);
1949
1950 DimVector distVec(globalPosFace - globalPos);
1951 Scalar dist = distVec.two_norm();
1952 DimVector& normal = interactionVolume.getNormal(elemIdx, fIdx);
1953
1954 // get pc and lambda at the boundary
1955 Scalar satWBound = cellData.saturation(wPhaseIdx);
1956 //check boundary sat at face 1
1957 if (interactionVolume.getBoundaryType(intVolFaceIdx).isDirichlet(satEqIdx))
1958 {
1959 Scalar satBound = interactionVolume.getDirichletValues(intVolFaceIdx)[saturationIdx];
1960 switch (saturationType_)
1961 {
1962 case sw:
1963 {
1964 satWBound = satBound;
1965 break;
1966 }
1967 case sn:
1968 {
1969 satWBound = 1 - satBound;
1970 break;
1971 }
1972 }
1973
1974 }
1975
1976 // old material law interface is deprecated: Replace this by
1977 // const auto& fluidMatrixInteraction = spatialParams.fluidMatrixInteractionAtPos(element.geometry().center());
1978 // after the release of 3.3, when the deprecated interface is no longer supported
1979 const auto fluidMatrixInteraction = Deprecated::makePcKrSw(Scalar{}, problem_.spatialParams(), element);
1980
1981 Scalar pcBound = fluidMatrixInteraction.pc(satWBound);
1982
1983 Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
1984 * (density_[nPhaseIdx] - density_[wPhaseIdx]);
1985
1986 pcBound += gravityDiffBound;
1987
1988 Dune::FieldVector<Scalar, numPhases> lambdaBound(fluidMatrixInteraction.krw(satWBound));
1989 lambdaBound[nPhaseIdx] = fluidMatrixInteraction.krn(satWBound);
1990 lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
1991 lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
1992
1993 Scalar gdeltaZ = (problem_.bBoxMax()-globalPosFace) * gravity_;
1994 Scalar potentialBoundW = interactionVolume.getDirichletValues(intVolFaceIdx)[pressureIdx] + density_[wPhaseIdx]*gdeltaZ;
1995 Scalar potentialBoundNw = potentialBoundW;
1996
1997 //calculate potential gradients
1998 switch (pressureType_)
1999 {
2000 case pw:
2001 {
2002 potentialBoundNw += pcBound;
2003 break;
2004 }
2005 case pn:
2006 {
2007 //calculate potential gradients
2008 potentialBoundW -= pcBound;
2009 break;
2010 }
2011 }
2012
2013 Scalar potentialDiffW = (potW - potentialBoundW) / dist;
2014 Scalar potentialDiffNw = (potNw - potentialBoundNw) / dist;
2015
2016 //store potentials for further calculations (saturation, ...)
2017 cellData.fluxData().addUpwindPotential(wPhaseIdx, boundaryFaceIdx, potentialDiffW);
2018 cellData.fluxData().addUpwindPotential(nPhaseIdx, boundaryFaceIdx, potentialDiffNw);
2019
2020 //calculated phase velocities from advective velocities -> capillary pressure velocity already added in pressure part!
2021 DimVector velocityW(0);
2022 DimVector velocityNw(0);
2023
2024 // calculate capillary pressure gradient
2025 DimVector potentialGradient = normal;
2026 potentialGradient *= (potW - potentialBoundW) / dist;
2027 permeability.mv(potentialGradient, velocityW);
2028
2029 potentialGradient = normal;
2030 potentialGradient *= (potNw - potentialBoundNw) / dist;
2031 permeability.mv(potentialGradient, velocityNw);
2032
2033 velocityW *= (potentialDiffW >= 0.) ? lambda[wPhaseIdx] : lambdaBound[wPhaseIdx];
2034 velocityNw *= (potentialDiffNw >= 0.) ? lambda[nPhaseIdx] : lambdaBound[nPhaseIdx];
2035
2036 //velocity is calculated from two vertices of one intersection!
2037 velocityW *= 0.25;
2038 velocityNw *= 0.25;
2039
2040 //store velocities
2041 cellData.fluxData().addVelocity(wPhaseIdx, boundaryFaceIdx, velocityW);
2042 cellData.fluxData().addVelocity(nPhaseIdx, boundaryFaceIdx, velocityNw);
2043 cellData.fluxData().setVelocityMarker(boundaryFaceIdx);
2044 }
2045 else if (interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressureEqIdx))
2046 {
2047 int boundaryFaceIdx = interactionVolume.getIndexOnElement(elemIdx, fIdx);
2048
2049 DimVector& normal = interactionVolume.getNormal(elemIdx, fIdx);
2050
2051 // get neumann boundary value
2052 PrimaryVariables boundValues(interactionVolume.getNeumannValues(intVolFaceIdx));
2053
2054 boundValues[wPhaseIdx] /= density_[wPhaseIdx];
2055 boundValues[nPhaseIdx] /= density_[nPhaseIdx];
2056
2057 DimVector velocityW(normal);
2058 DimVector velocityNw(normal);
2059
2060 velocityW *= boundValues[wPhaseIdx] / (4.0*interactionVolume.getFaceArea(elemIdx, fIdx));
2061 velocityNw *= boundValues[nPhaseIdx]
2062 / (4.0*interactionVolume.getFaceArea(elemIdx, fIdx));
2063
2064 //store potentials for further calculations (saturation, ...)
2065 cellData.fluxData().addUpwindPotential(wPhaseIdx, boundaryFaceIdx, boundValues[wPhaseIdx]);
2066 cellData.fluxData().addUpwindPotential(nPhaseIdx, boundaryFaceIdx, boundValues[nPhaseIdx]);
2067
2068 //store velocities
2069 cellData.fluxData().addVelocity(wPhaseIdx, boundaryFaceIdx, velocityW);
2070 cellData.fluxData().addVelocity(nPhaseIdx, boundaryFaceIdx, velocityNw);
2071 cellData.fluxData().setVelocityMarker(boundaryFaceIdx);
2072 }
2073 else
2074 {
2075 std::cout << "interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressureEqIdx)"
2076 << interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressureEqIdx) << "\n";
2077 DUNE_THROW(Dune::NotImplemented,
2078 "No valid boundary condition type defined for pressure equation!");
2079 }
2080 }
2081 }
2082}
2083
2084} // end namespace Dumux
2085#endif
deprecated.hh
Helpers for deprecation.
3dtransmissibilitycalculator.hh
Provides methods for transmissibility calculation in 3-d.
Dumux::normal
Vector normal(const Vector &v)
Create a vector normal to the given one (v is expected to be non-zero)
Definition: geometry/normal.hh:36
Dumux
Definition: adapt.hh:29
Dumux::GetProp
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type GetProp
get the type of a property (equivalent to old macro GET_PROP(...))
Definition: propertysystem.hh:140
Dumux::GetPropType
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
Dumux::IOName::viscosity
std::string viscosity(int phaseIdx) noexcept
I/O name of viscosity for multiphase systems.
Definition: name.hh:74
Dumux::IOName::permeability
std::string permeability() noexcept
I/O name of permeability.
Definition: name.hh:143
Dumux::IOName::density
std::string density(int phaseIdx) noexcept
I/O name of density for multiphase systems.
Definition: name.hh:65
Dumux::FvMpfaL3dTransmissibilityCalculator
Provides methods for transmissibility calculation in 3-d.
Definition: 3dtransmissibilitycalculator.hh:51
Dumux::FvMpfaL3dTransmissibilityCalculator::TransmissibilityType
Dune::FieldMatrix< Scalar, dim, 2 *dim - dim+1 > TransmissibilityType
Type of the transmissibility matrix.
Definition: 3dtransmissibilitycalculator.hh:78
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)
Dumux::FvMpfaL3dVelocity2p
Class for calculating 3-d velocities from cell-wise constant pressure values.
Definition: 3dvelocity.hh:57
Dumux::FvMpfaL3dVelocity2p::FvMpfaL3dVelocity2p
FvMpfaL3dVelocity2p(Problem &problem)
Constructs a FvMpfaL3dVelocity2p object.
Definition: 3dvelocity.hh:139
Dumux::FvMpfaL3dVelocity2p::initialize
void initialize(bool solveTwice=true)
Initializes the velocity model.
Definition: 3dvelocity.hh:160
Dumux::FvMpfaL3dVelocity2p::addOutputVtkFields
void addOutputVtkFields(MultiWriter &writer)
Adds velocity output to the output file.
Definition: 3dvelocity.hh:188
Dumux::FvMpfaL3dVelocity2p::calculateInnerInteractionVolumeVelocity
void calculateInnerInteractionVolumeVelocity(InteractionVolume &interactionVolume, CellData &cellData1, CellData &cellData2, CellData &cellData3, CellData &cellData4, CellData &cellData5, CellData &cellData6, CellData &cellData7, CellData &cellData8, InteractionVolumeContainer &interactionVolumes, TransmissibilityCalculator &transmissibilityCalculator, int fIdx=-1)
Calculate velocities for flux faces of an interaction volume.
Definition: 3dvelocity.hh:296
Dumux::FvMpfaL3dVelocity2p::calculateBoundaryInteractionVolumeVelocity
void calculateBoundaryInteractionVolumeVelocity(InteractionVolume &interactionVolume, CellData &cellData, int elemIdx)
Calculates the velocity at a boundary flux faces.
Definition: 3dvelocity.hh:1921
properties.hh
Specifies the properties for immiscible 2p diffusion/pressure models.
properties.hh
Properties for a MPFA method.
Generated by doxygen 1.9.3