3.2-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
3dpressure.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 *****************************************************************************/
26#ifndef DUMUX_FVMPFAL2PFABOUND3DPRESSURE2P_HH
27#define DUMUX_FVMPFAL2PFABOUND3DPRESSURE2P_HH
28
29// dumux environment
30#include <dumux/porousmediumflow/2p/sequential/diffusion/properties.hh>
31#include <dumux/porousmediumflow/sequential/cellcentered/pressure.hh>
32#include "3dinteractionvolumecontainer.hh"
33#include "3dtransmissibilitycalculator.hh"
34
35namespace Dumux {
36
72template<class TypeTag>
73class FvMpfaL3dPressure2p: public FVPressure<TypeTag>
74{
75 using ParentType = FVPressure<TypeTag>;
76 using Implementation = GetPropType<TypeTag, Properties::PressureModel>;
77 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
78
79 enum
80 {
81 dim = GridView::dimension, dimWorld = GridView::dimensionworld
82 };
83
84 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
85 using Problem = GetPropType<TypeTag, Properties::Problem>;
86
87 using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
88 using MaterialLaw = typename SpatialParams::MaterialLaw;
89
90 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
91
92 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
93 using FluidState = GetPropType<TypeTag, Properties::FluidState>;
94
95 using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
96 using CellData = GetPropType<TypeTag, Properties::CellData>;
97
98 using SolutionTypes = GetProp<TypeTag, Properties::SolutionTypes>;
99 using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
100 using ScalarSolutionType = typename SolutionTypes::ScalarSolution;
101
102 enum
103 {
104 pw = Indices::pressureW,
105 pn = Indices::pressureNw,
106 pGlobal = Indices::pressureGlobal,
107 sw = Indices::saturationW,
108 sn = Indices::saturationNw,
109 vw = Indices::velocityW,
110 vn = Indices::velocityNw,
111 vt = Indices::velocityTotal
112 };
113 enum
114 {
115 wPhaseIdx = Indices::wPhaseIdx,
116 nPhaseIdx = Indices::nPhaseIdx,
117 pressureIdx = Indices::pressureIdx,
118 saturationIdx = Indices::saturationIdx,
119 pressureEqIdx = Indices::pressureEqIdx,
120 satEqIdx = Indices::satEqIdx,
121 numPhases = getPropValue<TypeTag, Properties::NumPhases>()
122 };
123 enum
124 {
125 globalCorner = 2,
126 globalEdge = 3,
127 neumannNeumann = 0,
128 dirichletDirichlet = 1,
129 dirichletNeumann = 2,
130 neumannDirichlet = 3
131 };
132 enum
133 {
134 innerEdgeFace = 2,
135 innerSideFace = 1
136 };
137
138 using Element = typename GridView::Traits::template Codim<0>::Entity;
139 using Grid = typename GridView::Grid;
140 using Geometry = typename Element::Geometry;
141 using IntersectionIterator = typename GridView::IntersectionIterator;
142
143 using GlobalPosition = typename Geometry::GlobalCoordinate;
144 using GravityVector = Dune::FieldVector<Scalar, dimWorld>;
145 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
146
147 using DimVector = Dune::FieldVector<Scalar, dim>;
148
149 using InteractionVolumeContainer = GetPropType<TypeTag, Properties::MPFAInteractionVolumeContainer>;
150 using TransmissibilityCalculator = FvMpfaL3dTransmissibilityCalculator<TypeTag>;
151public:
153 using TransmissibilityType = typename TransmissibilityCalculator::TransmissibilityType;
155 using InteractionVolume = GetPropType<TypeTag, Properties::MPFAInteractionVolume>;
156protected:
158 friend class FVPressure<TypeTag>;
159 void initializeMatrix();
160 void initializeMatrixRowSize();
161 void initializeMatrixIndices();
162
164 void assemble();
165 void assembleInnerInteractionVolume(InteractionVolume& interactionVolume);
166 void assembleBoundaryInteractionVolume(InteractionVolume& interactionVolume);
167
168public:
169
171 void updateMaterialLaws();
172
178 void initialize(bool solveTwice = true)
179 {
180 const auto element = *problem_.gridView().template begin<0>();
181 FluidState fluidState;
182 fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
183 fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
184 fluidState.setTemperature(problem_.temperature(element));
185 fluidState.setSaturation(wPhaseIdx, 1.);
186 fluidState.setSaturation(nPhaseIdx, 0.);
187 density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
188 density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
189 viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
190 viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
191
192 interactionVolumes_.initialize();
193 ParentType::initialize();
194
195 updateMaterialLaws();
196
197 asImp_().assemble();
198 this->solve();
199
200 storePressureSolution();
201
202 return;
203 }
204
206 void storePressureSolution()
207 {
208 // iterate through leaf grid an evaluate c0 at cell center
209 for (const auto& element : elements(problem_.gridView()))
210 {
211 storePressureSolution(element);
212 }
213 }
214
220 void storePressureSolution(const Element& element)
221 {
222 int eIdxGlobal = problem_.variables().index(element);
223 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
224
225 switch (pressureType_)
226 {
227 case pw:
228 {
229 Scalar potW = this->pressure()[eIdxGlobal];
230
231 Scalar gravityDiff = (problem_.bBoxMax() - element.geometry().center()) * gravity_;
232 Scalar potPc = cellData.capillaryPressure()
233 + gravityDiff * (density_[nPhaseIdx] - density_[wPhaseIdx]);
234
235 cellData.setPotential(wPhaseIdx, potW);
236 cellData.setPotential(nPhaseIdx, potW + potPc);
237
238 Scalar pressW = potW - gravityDiff * density_[wPhaseIdx];
239
240 cellData.setPressure(wPhaseIdx, pressW);
241 cellData.setPressure(nPhaseIdx, pressW + cellData.capillaryPressure());
242
243 break;
244 }
245 case pn:
246 {
247 Scalar potNw = this->pressure()[eIdxGlobal];
248
249 Scalar gravityDiff = (problem_.bBoxMax() - element.geometry().center()) * gravity_;
250 Scalar potPc = cellData.capillaryPressure()
251 + gravityDiff * (density_[nPhaseIdx] - density_[wPhaseIdx]);
252
253 cellData.setPotential(nPhaseIdx, potNw);
254 cellData.setPotential(wPhaseIdx, potNw - potPc);
255
256 Scalar pressNw = potNw - gravityDiff * density_[nPhaseIdx];
257
258 cellData.setPressure(wPhaseIdx, pressNw - cellData.capillaryPressure());
259 cellData.setPressure(nPhaseIdx, pressNw);
260
261 break;
262 }
263 }
264 cellData.fluxData().resetVelocity();
265 }
266
272 void update()
273 {
274 int size = problem_.gridView().size(0);
275
276 //error bounds for error term for incompressible models to correct unphysical saturation
277 //over/undershoots due to saturation transport
278 timeStep_ = problem_.timeManager().timeStepSize();
279 maxError_ = 0.0;
280 for (int i = 0; i < size; i++)
281 {
282 CellData& cellData = problem_.variables().cellData(i);
283 Scalar sat = 0;
284 using std::max;
285 switch (saturationType_)
286 {
287 case sw:
288 sat = cellData.saturation(wPhaseIdx);
289 break;
290 case sn:
291 sat = cellData.saturation(nPhaseIdx);
292 break;
293 }
294 if (sat > 1.0)
295 {
296 maxError_ = max(maxError_, (sat - 1.0) / timeStep_);
297 }
298 if (sat < 0.0)
299 {
300 maxError_ = max(maxError_, (-sat) / timeStep_);
301 }
302
303 switch (pressureType_)
304 {
305 case pw:
306 this->pressure()[i] = cellData.pressure(wPhaseIdx);
307 break;
308 case pn:
309 this->pressure()[i] = cellData.pressure(nPhaseIdx);
310 break;
311 }
312 }
313
314 asImp_().assemble();
315
316 this->solve();
317
318 storePressureSolution();
319
320 return;
321 }
322
341 Scalar evaluateErrorTerm(CellData& cellData)
342 {
343 //error term for incompressible models to correct unphysical saturation over/undershoots due to saturation transport
344 // error reduction routine: volumetric error is damped and inserted to right hand side
345 Scalar sat = 0;
346 switch (saturationType_)
347 {
348 case sw:
349 sat = cellData.saturation(wPhaseIdx);
350 break;
351 case sn:
352 sat = cellData.saturation(nPhaseIdx);
353 break;
354 }
355
356 Scalar error = (sat > 1.0) ? sat - 1.0 : 0.0;
357 if (sat < 0.0)
358 {
359 error = sat;
360 }
361 error /= timeStep_;
362
363 using std::abs;
364 Scalar errorAbs = abs(error);
365
366 if ((errorAbs * timeStep_ > 1e-6) && (errorAbs > ErrorTermLowerBound_ * maxError_)
367 && (!problem_.timeManager().willBeFinished()))
368 {
369 return ErrorTermFactor_ * error;
370 }
371 return 0.0;
372 }
373
384 template<class MultiWriter>
385 void addOutputVtkFields(MultiWriter &writer)
386 {
387 int size = problem_.gridView().size(0);
388 ScalarSolutionType *potential = writer.allocateManagedBuffer(size);
389
390 (*potential) = this->pressure();
391
392 if (pressureType_ == pw)
393 {
394 writer.attachCellData(*potential, "wetting potential");
395 }
396
397 if (pressureType_ == pn)
398 {
399 writer.attachCellData(*potential, "nonwetting potential");
400 }
401
402 if (vtkOutputLevel_ > 0)
403 {
404 ScalarSolutionType *pressure = writer.allocateManagedBuffer(size);
405 ScalarSolutionType *pressureSecond = writer.allocateManagedBuffer(size);
406 ScalarSolutionType *potentialSecond = writer.allocateManagedBuffer(size);
407 ScalarSolutionType *pc = writer.allocateManagedBuffer(size);
408
409 for (const auto& element : elements(problem_.gridView()))
410 {
411 int idx = problem_.variables().index(element);
412 CellData& cellData = problem_.variables().cellData(idx);
413
414 (*pc)[idx] = cellData.capillaryPressure();
415
416 if (pressureType_ == pw)
417 {
418 (*pressure)[idx] = cellData.pressure(wPhaseIdx);
419 (*potentialSecond)[idx] = cellData.potential(nPhaseIdx);
420 (*pressureSecond)[idx] = cellData.pressure(nPhaseIdx);
421 }
422
423 if (pressureType_ == pn)
424 {
425 (*pressure)[idx] = cellData.pressure(nPhaseIdx);
426 (*potentialSecond)[idx] = cellData.potential(wPhaseIdx);
427 (*pressureSecond)[idx] = cellData.pressure(wPhaseIdx);
428 }
429 }
430
431 if (pressureType_ == pw)
432 {
433 writer.attachCellData(*pressure, "wetting pressure");
434 writer.attachCellData(*pressureSecond, "nonwetting pressure");
435 writer.attachCellData(*potentialSecond, "nonwetting potential");
436 }
437
438 if (pressureType_ == pn)
439 {
440 writer.attachCellData(*pressure, "nonwetting pressure");
441 writer.attachCellData(*pressureSecond, "wetting pressure");
442 writer.attachCellData(*potentialSecond, "wetting potential");
443 }
444
445 writer.attachCellData(*pc, "capillary pressure");
446 }
447 }
448
450 InteractionVolumeContainer& interactionVolumes()
451 {
452 return interactionVolumes_;
453 }
454
460 TransmissibilityCalculator& transmissibilityCalculator()
461 {
462 return transmissibilityCalculator_;
463 }
464
469 FvMpfaL3dPressure2p(Problem& problem) :
470 ParentType(problem), problem_(problem),
471 interactionVolumes_(problem), transmissibilityCalculator_(problem),
472 gravity_(problem.gravity()),
473 maxError_(0.), timeStep_(1.)
474 {
475 if (pressureType_ != pw && pressureType_ != pn)
476 {
477 DUNE_THROW(Dune::NotImplemented, "Pressure type not supported!");
478 }
479 if (saturationType_ != sw && saturationType_ != sn)
480 {
481 DUNE_THROW(Dune::NotImplemented, "Saturation type not supported!");
482 }
483 if (getPropValue<TypeTag, Properties::EnableCompressibility>())
484 {
485 DUNE_THROW(Dune::NotImplemented, "Compressibility not supported!");
486 }
487 if (dim != 3)
488 {
489 DUNE_THROW(Dune::NotImplemented, "Dimension not supported!");
490 }
491
492 ErrorTermFactor_ = getParam<Scalar>("Impet.ErrorTermFactor");
493 ErrorTermLowerBound_ = getParam<Scalar>("Impet.ErrorTermLowerBound");
494 ErrorTermUpperBound_ = getParam<Scalar>("Impet.ErrorTermUpperBound");
495
496 density_[wPhaseIdx] = 0.;
497 density_[nPhaseIdx] = 0.;
498 viscosity_[wPhaseIdx] = 0.;
499 viscosity_[nPhaseIdx] = 0.;
500
501 vtkOutputLevel_ = getParam<int>("Vtk.OutputLevel");
502 }
503
504private:
505 Problem& problem_;
506
507protected:
508 InteractionVolumeContainer interactionVolumes_;
510 TransmissibilityCalculator transmissibilityCalculator_;
511
512private:
514 Implementation &asImp_()
515 { return *static_cast<Implementation *>(this);}
516
518 const Implementation &asImp_() const
519 { return *static_cast<const Implementation *>(this);}
520
521 const GravityVector& gravity_;
522
523 Scalar maxError_;
524 Scalar timeStep_;
525 Scalar ErrorTermFactor_;
526 Scalar ErrorTermLowerBound_;
527 Scalar ErrorTermUpperBound_;
528
529 Scalar density_[numPhases];
530 Scalar viscosity_[numPhases];
531
532 int vtkOutputLevel_;
533
534 static constexpr Scalar threshold_ = 1e-15;
536 static const int pressureType_ = getPropValue<TypeTag, Properties::PressureFormulation>();
538 static const int saturationType_ = getPropValue<TypeTag, Properties::SaturationFormulation>();
540 static const int velocityType_ = getPropValue<TypeTag, Properties::VelocityFormulation>();
541};
542
544template<class TypeTag>
545void FvMpfaL3dPressure2p<TypeTag>::initializeMatrix()
546{
547 initializeMatrixRowSize();
548 this->A_.endrowsizes();
549 initializeMatrixIndices();
550 this->A_.endindices();
551}
552
554template<class TypeTag>
555void FvMpfaL3dPressure2p<TypeTag>::initializeMatrixRowSize()
556{
557 // determine matrix row sizes
558 for (const auto& element : elements(problem_.gridView()))
559 {
560 // cell index
561 int eIdxGlobalI = problem_.variables().index(element);
562
563 std::set<int> neighborIndices;
564
565 int numVertices = element.geometry().corners();
566
567 for (int vIdx = 0; vIdx < numVertices; vIdx++)
568 {
569 int vIdxGlobal = problem_.variables().vertexMapper().subIndex(element, vIdx, dim);
570
571 InteractionVolume& interactionVolume = interactionVolumes_.interactionVolume(vIdxGlobal);
572
573 for (int subVolumeIdx = 0; subVolumeIdx < InteractionVolume::subVolumeTotalNum; subVolumeIdx++)
574 {
575 if (interactionVolume.hasSubVolumeElement(subVolumeIdx))
576 {
577 int eIdxGlobalJ = problem_.variables().index(interactionVolume.getSubVolumeElement(subVolumeIdx));
578
579 neighborIndices.insert(eIdxGlobalJ);
580 }
581 }
582 }
583
584 this->A_.setrowsize(eIdxGlobalI, neighborIndices.size());
585 } // end of element loop
586
587 return;
588}
589
591template<class TypeTag>
592void FvMpfaL3dPressure2p<TypeTag>::initializeMatrixIndices()
593{
594 // determine position of matrix entries
595 for (const auto& element : elements(problem_.gridView()))
596 {
597 // cell index
598 int eIdxGlobalI = problem_.variables().index(element);
599
600 // add diagonal index
601 this->A_.addindex(eIdxGlobalI, eIdxGlobalI);
602
603 int numVertices = element.geometry().corners();
604
605 for (int vIdx = 0; vIdx < numVertices; vIdx++)
606 {
607 int vIdxGlobal = problem_.variables().vertexMapper().subIndex(element, vIdx, dim);
608
609 InteractionVolume& interactionVolume = interactionVolumes_.interactionVolume(vIdxGlobal);
610 for (int subVolumeIdx = 0; subVolumeIdx < InteractionVolume::subVolumeTotalNum; subVolumeIdx++)
611 {
612 if (interactionVolume.hasSubVolumeElement(subVolumeIdx))
613 {
614 int eIdxGlobalJ = problem_.variables().index(interactionVolume.getSubVolumeElement(subVolumeIdx));
615
616 this->A_.addindex(eIdxGlobalI, eIdxGlobalJ);
617 }
618 }
619 }
620 } // end of element loop
621 return;
622}
623
625template<class TypeTag>
626void FvMpfaL3dPressure2p<TypeTag>::assemble()
627{
628 // initialization: set global matrix this->A_ to zero
629 this->A_ = 0;
630 this->f_ = 0;
631
632 // run through all vertices
633 for (const auto& vertex : vertices(problem_.gridView()))
634 {
635 int vIdxGlobal = problem_.variables().index(vertex);
636
637 InteractionVolume& interactionVolume = interactionVolumes_.interactionVolume(vIdxGlobal);
638
639 // inner interactionvolume
640 if (interactionVolume.isInnerVolume())
641 {
642 assembleInnerInteractionVolume(interactionVolume);
643 }
644
645 // at least one face on boundary! (boundary interactionvolume)
646 else
647 {
648 assembleBoundaryInteractionVolume(interactionVolume);
649 } // end boundaries
650
651 } // end vertex iterator
652
653 return;
654}
655
657template<class TypeTag>
658void FvMpfaL3dPressure2p<TypeTag>::assembleInnerInteractionVolume(InteractionVolume& interactionVolume)
659{
660 auto element1 = interactionVolume.getSubVolumeElement(0);
661 auto element2 = interactionVolume.getSubVolumeElement(1);
662 auto element3 = interactionVolume.getSubVolumeElement(2);
663 auto element4 = interactionVolume.getSubVolumeElement(3);
664 auto element5 = interactionVolume.getSubVolumeElement(4);
665 auto element6 = interactionVolume.getSubVolumeElement(5);
666 auto element7 = interactionVolume.getSubVolumeElement(6);
667 auto element8 = interactionVolume.getSubVolumeElement(7);
668
669 // get global coordinate of cell centers
670 const GlobalPosition& globalPos1 = element1.geometry().center();
671 const GlobalPosition& globalPos2 = element2.geometry().center();
672 const GlobalPosition& globalPos3 = element3.geometry().center();
673 const GlobalPosition& globalPos4 = element4.geometry().center();
674 const GlobalPosition& globalPos5 = element5.geometry().center();
675 const GlobalPosition& globalPos6 = element6.geometry().center();
676 const GlobalPosition& globalPos7 = element7.geometry().center();
677 const GlobalPosition& globalPos8 = element8.geometry().center();
678
679 // cell volumes
680 Scalar volume1 = element1.geometry().volume();
681 Scalar volume2 = element2.geometry().volume();
682 Scalar volume3 = element3.geometry().volume();
683 Scalar volume4 = element4.geometry().volume();
684 Scalar volume5 = element5.geometry().volume();
685 Scalar volume6 = element6.geometry().volume();
686 Scalar volume7 = element7.geometry().volume();
687 Scalar volume8 = element8.geometry().volume();
688
689 // cell index
690 int eIdxGlobal1 = problem_.variables().index(element1);
691 int eIdxGlobal2 = problem_.variables().index(element2);
692 int eIdxGlobal3 = problem_.variables().index(element3);
693 int eIdxGlobal4 = problem_.variables().index(element4);
694 int eIdxGlobal5 = problem_.variables().index(element5);
695 int eIdxGlobal6 = problem_.variables().index(element6);
696 int eIdxGlobal7 = problem_.variables().index(element7);
697 int eIdxGlobal8 = problem_.variables().index(element8);
698
699 //get the cell Data
700 CellData& cellData1 = problem_.variables().cellData(eIdxGlobal1);
701 CellData& cellData2 = problem_.variables().cellData(eIdxGlobal2);
702 CellData& cellData3 = problem_.variables().cellData(eIdxGlobal3);
703 CellData& cellData4 = problem_.variables().cellData(eIdxGlobal4);
704 CellData& cellData5 = problem_.variables().cellData(eIdxGlobal5);
705 CellData& cellData6 = problem_.variables().cellData(eIdxGlobal6);
706 CellData& cellData7 = problem_.variables().cellData(eIdxGlobal7);
707 CellData& cellData8 = problem_.variables().cellData(eIdxGlobal8);
708
709 // get mobilities of the phases
710 Dune::FieldVector<Scalar, numPhases> lambda1(cellData1.mobility(wPhaseIdx));
711 lambda1[nPhaseIdx] = cellData1.mobility(nPhaseIdx);
712
713 // compute total mobility of cell 1
714 Scalar lambdaTotal1 = lambda1[wPhaseIdx] + lambda1[nPhaseIdx];
715
716 // get mobilities of the phases
717 Dune::FieldVector<Scalar, numPhases> lambda2(cellData2.mobility(wPhaseIdx));
718 lambda2[nPhaseIdx] = cellData2.mobility(nPhaseIdx);
719
720 // compute total mobility of cell 2
721 Scalar lambdaTotal2 = lambda2[wPhaseIdx] + lambda2[nPhaseIdx];
722
723 // get mobilities of the phases
724 Dune::FieldVector<Scalar, numPhases> lambda3(cellData3.mobility(wPhaseIdx));
725 lambda3[nPhaseIdx] = cellData3.mobility(nPhaseIdx);
726
727 // compute total mobility of cell 3
728 Scalar lambdaTotal3 = lambda3[wPhaseIdx] + lambda3[nPhaseIdx];
729
730 // get mobilities of the phases
731 Dune::FieldVector<Scalar, numPhases> lambda4(cellData4.mobility(wPhaseIdx));
732 lambda4[nPhaseIdx] = cellData4.mobility(nPhaseIdx);
733
734 // compute total mobility of cell 4
735 Scalar lambdaTotal4 = lambda4[wPhaseIdx] + lambda4[nPhaseIdx];
736
737 // get mobilities of the phases
738 Dune::FieldVector<Scalar, numPhases> lambda5(cellData5.mobility(wPhaseIdx));
739 lambda5[nPhaseIdx] = cellData5.mobility(nPhaseIdx);
740
741 // compute total mobility of cell 5
742 Scalar lambdaTotal5 = lambda5[wPhaseIdx] + lambda5[nPhaseIdx];
743
744 // get mobilities of the phases
745 Dune::FieldVector<Scalar, numPhases> lambda6(cellData6.mobility(wPhaseIdx));
746 lambda6[nPhaseIdx] = cellData6.mobility(nPhaseIdx);
747
748 // compute total mobility of cell 6
749 Scalar lambdaTotal6 = lambda6[wPhaseIdx] + lambda6[nPhaseIdx];
750
751 // get mobilities of the phases
752 Dune::FieldVector<Scalar, numPhases> lambda7(cellData7.mobility(wPhaseIdx));
753 lambda7[nPhaseIdx] = cellData7.mobility(nPhaseIdx);
754
755 // compute total mobility of cell 7
756 Scalar lambdaTotal7 = lambda7[wPhaseIdx] + lambda7[nPhaseIdx];
757
758 // get mobilities of the phases
759 Dune::FieldVector<Scalar, numPhases> lambda8(cellData8.mobility(wPhaseIdx));
760 lambda8[nPhaseIdx] = cellData8.mobility(nPhaseIdx);
761
762 // compute total mobility of cell 8
763 Scalar lambdaTotal8 = lambda8[wPhaseIdx] + lambda8[nPhaseIdx];
764
765 std::vector<DimVector> lambda(8);
766 lambda[0][0] = lambdaTotal1;
767 lambda[0][1] = lambdaTotal1;
768 lambda[0][2] = lambdaTotal1;
769 lambda[1][0] = lambdaTotal2;
770 lambda[1][1] = lambdaTotal2;
771 lambda[1][2] = lambdaTotal2;
772 lambda[2][0] = lambdaTotal3;
773 lambda[2][1] = lambdaTotal3;
774 lambda[2][2] = lambdaTotal3;
775 lambda[3][0] = lambdaTotal4;
776 lambda[3][1] = lambdaTotal4;
777 lambda[3][2] = lambdaTotal4;
778 lambda[4][0] = lambdaTotal5;
779 lambda[4][1] = lambdaTotal5;
780 lambda[4][2] = lambdaTotal5;
781 lambda[5][0] = lambdaTotal6;
782 lambda[5][1] = lambdaTotal6;
783 lambda[5][2] = lambdaTotal6;
784 lambda[6][0] = lambdaTotal7;
785 lambda[6][1] = lambdaTotal7;
786 lambda[6][2] = lambdaTotal7;
787 lambda[7][0] = lambdaTotal8;
788 lambda[7][1] = lambdaTotal8;
789 lambda[7][2] = lambdaTotal8;
790
791 // add capillary pressure and gravity terms to right-hand-side
792 // calculate capillary pressure velocity
793 Dune::FieldVector<Scalar, 8> pc(0);
794 pc[0] = cellData1.capillaryPressure();
795 pc[1] = cellData2.capillaryPressure();
796 pc[2] = cellData3.capillaryPressure();
797 pc[3] = cellData4.capillaryPressure();
798 pc[4] = cellData5.capillaryPressure();
799 pc[5] = cellData6.capillaryPressure();
800 pc[6] = cellData7.capillaryPressure();
801 pc[7] = cellData8.capillaryPressure();
802
803 Dune::FieldVector<Scalar, 8> gravityDiff(0);
804
805 gravityDiff[0] = (problem_.bBoxMax() - globalPos1) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
806 gravityDiff[1] = (problem_.bBoxMax() - globalPos2) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
807 gravityDiff[2] = (problem_.bBoxMax() - globalPos3) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
808 gravityDiff[3] = (problem_.bBoxMax() - globalPos4) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
809 gravityDiff[4] = (problem_.bBoxMax() - globalPos5) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
810 gravityDiff[5] = (problem_.bBoxMax() - globalPos6) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
811 gravityDiff[6] = (problem_.bBoxMax() - globalPos7) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
812 gravityDiff[7] = (problem_.bBoxMax() - globalPos8) * gravity_ * (density_[nPhaseIdx] - density_[wPhaseIdx]);
813
814 pc += gravityDiff;
815
816 Dune::FieldVector<Dune::FieldVector<Scalar, 3>, 8> pcFlux(Dune::FieldVector<Scalar, 3>(0));
817
818 Scalar pcPotential0 = 0;
819 Scalar pcPotential1 = 0;
820 Scalar pcPotential2 = 0;
821 Scalar pcPotential3 = 0;
822 Scalar pcPotential4 = 0;
823 Scalar pcPotential5 = 0;
824 Scalar pcPotential6 = 0;
825 Scalar pcPotential7 = 0;
826 Scalar pcPotential8 = 0;
827 Scalar pcPotential9 = 0;
828 Scalar pcPotential10 = 0;
829 Scalar pcPotential11 = 0;
830
831 // interactionVolume.printInteractionVolumeInfo();
832 // evaluate right hand side
833 PrimaryVariables source(0.0);
834 problem_.source(source, element1);
835 this->f_[eIdxGlobal1] += volume1 / (8.0)
836 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
837 problem_.source(source, element2);
838 this->f_[eIdxGlobal2] += volume2 / (8.0)
839 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
840 problem_.source(source, element3);
841 this->f_[eIdxGlobal3] += volume3 / (8.0)
842 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
843 problem_.source(source, element4);
844 this->f_[eIdxGlobal4] += volume4 / (8.0)
845 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
846 problem_.source(source, element5);
847 this->f_[eIdxGlobal5] += volume5 / (8.0)
848 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
849 problem_.source(source, element6);
850 this->f_[eIdxGlobal6] += volume6 / (8.0)
851 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
852 problem_.source(source, element7);
853 this->f_[eIdxGlobal7] += volume7 / (8.0)
854 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
855 problem_.source(source, element8);
856 this->f_[eIdxGlobal8] += volume8 / (8.0)
857 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
858
859 this->f_[eIdxGlobal1] += evaluateErrorTerm(cellData1) * volume1 / (8.0);
860 this->f_[eIdxGlobal2] += evaluateErrorTerm(cellData2) * volume2 / (8.0);
861 this->f_[eIdxGlobal3] += evaluateErrorTerm(cellData3) * volume3 / (8.0);
862 this->f_[eIdxGlobal4] += evaluateErrorTerm(cellData4) * volume4 / (8.0);
863 this->f_[eIdxGlobal5] += evaluateErrorTerm(cellData5) * volume5 / (8.0);
864 this->f_[eIdxGlobal6] += evaluateErrorTerm(cellData6) * volume6 / (8.0);
865 this->f_[eIdxGlobal7] += evaluateErrorTerm(cellData7) * volume7 / (8.0);
866 this->f_[eIdxGlobal8] += evaluateErrorTerm(cellData8) * volume8 / (8.0);
867
868 DimVector Tu(0);
869 Dune::FieldVector<Scalar, 2 * dim - dim + 1> u(0);
870 TransmissibilityType T(0);
871 TransmissibilityType TSecond(0);
872
873 // calculate the flux through the subvolumeface 1 (subVolumeFaceIdx = 0)
874 int caseL = transmissibilityCalculator_.transmissibility(T, interactionVolume, lambda, 0, 1, 2, 3, 4,
875 5);
876
877 TSecond = T;
878 T *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal1, 0, 0);
879 TSecond *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal2, 1, 1);
880
881 if (caseL == 1)
882 {
883 this->A_[eIdxGlobal1][eIdxGlobal1] += T[0][0];
884 this->A_[eIdxGlobal1][eIdxGlobal2] += T[0][1];
885 this->A_[eIdxGlobal1][eIdxGlobal3] += T[0][2];
886 this->A_[eIdxGlobal1][eIdxGlobal5] += T[0][3];
887
888 this->A_[eIdxGlobal2][eIdxGlobal1] -= TSecond[0][0];
889 this->A_[eIdxGlobal2][eIdxGlobal2] -= TSecond[0][1];
890 this->A_[eIdxGlobal2][eIdxGlobal3] -= TSecond[0][2];
891 this->A_[eIdxGlobal2][eIdxGlobal5] -= TSecond[0][3];
892
893 u[0] = pc[0];
894 u[1] = pc[1];
895 u[2] = pc[2];
896 u[3] = pc[4];
897
898 T.mv(u, Tu);
899
900 pcFlux[0][0] = Tu[0];
901 pcPotential0 = Tu[0];
902
903 TSecond.mv(u, Tu);
904 pcFlux[1][1]= Tu[0];
905 }
906 else if (caseL == 2)
907 {
908 this->A_[eIdxGlobal1][eIdxGlobal1] += T[0][0];
909 this->A_[eIdxGlobal1][eIdxGlobal2] += T[0][1];
910 this->A_[eIdxGlobal1][eIdxGlobal4] += T[0][2];
911 this->A_[eIdxGlobal1][eIdxGlobal6] += T[0][3];
912
913 this->A_[eIdxGlobal2][eIdxGlobal1] -= TSecond[0][0];
914 this->A_[eIdxGlobal2][eIdxGlobal2] -= TSecond[0][1];
915 this->A_[eIdxGlobal2][eIdxGlobal4] -= TSecond[0][2];
916 this->A_[eIdxGlobal2][eIdxGlobal6] -= TSecond[0][3];
917
918 u[0] = pc[0];
919 u[1] = pc[1];
920 u[2] = pc[3];
921 u[3] = pc[5];
922
923 T.mv(u, Tu);
924 pcFlux[0][0] = Tu[0];
925 pcPotential0 = Tu[0];
926
927 TSecond.mv(u, Tu);
928 pcFlux[1][1]= Tu[0];
929 }
930 else if (caseL == 3)
931 {
932 this->A_[eIdxGlobal1][eIdxGlobal1] += T[0][0];
933 this->A_[eIdxGlobal1][eIdxGlobal2] += T[0][1];
934 this->A_[eIdxGlobal1][eIdxGlobal4] += T[0][2];
935 this->A_[eIdxGlobal1][eIdxGlobal5] += T[0][3];
936
937 this->A_[eIdxGlobal2][eIdxGlobal1] -= TSecond[0][0];
938 this->A_[eIdxGlobal2][eIdxGlobal2] -= TSecond[0][1];
939 this->A_[eIdxGlobal2][eIdxGlobal4] -= TSecond[0][2];
940 this->A_[eIdxGlobal2][eIdxGlobal5] -= TSecond[0][3];
941
942 u[0] = pc[0];
943 u[1] = pc[1];
944 u[2] = pc[3];
945 u[3] = pc[4];
946
947 T.mv(u, Tu);
948 pcFlux[0][0] = Tu[0];
949 pcPotential0 = Tu[0];
950
951 TSecond.mv(u, Tu);
952 pcFlux[1][1]= Tu[0];
953 }
954 else
955 {
956 this->A_[eIdxGlobal1][eIdxGlobal1] += T[0][0];
957 this->A_[eIdxGlobal1][eIdxGlobal2] += T[0][1];
958 this->A_[eIdxGlobal1][eIdxGlobal3] += T[0][2];
959 this->A_[eIdxGlobal1][eIdxGlobal6] += T[0][3];
960
961 this->A_[eIdxGlobal2][eIdxGlobal1] -= TSecond[0][0];
962 this->A_[eIdxGlobal2][eIdxGlobal2] -= TSecond[0][1];
963 this->A_[eIdxGlobal2][eIdxGlobal3] -= TSecond[0][2];
964 this->A_[eIdxGlobal2][eIdxGlobal6] -= TSecond[0][3];
965
966 u[0] = pc[0];
967 u[1] = pc[1];
968 u[2] = pc[2];
969 u[3] = pc[5];
970
971 T.mv(u, Tu);
972 pcFlux[0][0] = Tu[0];
973 pcPotential0 = Tu[0];
974
975 TSecond.mv(u, Tu);
976 pcFlux[1][1]= Tu[0];
977 }
978
979 // calculate the flux through the subvolumeface 2 (subVolumeFaceIdx = 1)
980 caseL = transmissibilityCalculator_.transmissibility(T, interactionVolume, lambda, 1, 3, 0, 2, 5, 7);
981
982 TSecond = T;
983 T *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal2, 1, 0);
984 TSecond *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal4, 3, 1);
985
986 if (caseL == 1)
987 {
988 this->A_[eIdxGlobal2][eIdxGlobal2] += T[0][0];
989 this->A_[eIdxGlobal2][eIdxGlobal4] += T[0][1];
990 this->A_[eIdxGlobal2][eIdxGlobal1] += T[0][2];
991 this->A_[eIdxGlobal2][eIdxGlobal6] += T[0][3];
992
993 this->A_[eIdxGlobal4][eIdxGlobal2] -= TSecond[0][0];
994 this->A_[eIdxGlobal4][eIdxGlobal4] -= TSecond[0][1];
995 this->A_[eIdxGlobal4][eIdxGlobal1] -= TSecond[0][2];
996 this->A_[eIdxGlobal4][eIdxGlobal6] -= TSecond[0][3];
997
998 u[0] = pc[1];
999 u[1] = pc[3];
1000 u[2] = pc[0];
1001 u[3] = pc[5];
1002
1003 T.mv(u, Tu);
1004 pcFlux[1][0] = Tu[0];
1005 pcPotential1 = Tu[0];
1006
1007 TSecond.mv(u, Tu);
1008 pcFlux[3][1] = Tu[0];
1009 }
1010 else if (caseL == 2)
1011 {
1012 this->A_[eIdxGlobal2][eIdxGlobal2] += T[0][0];
1013 this->A_[eIdxGlobal2][eIdxGlobal4] += T[0][1];
1014 this->A_[eIdxGlobal2][eIdxGlobal3] += T[0][2];
1015 this->A_[eIdxGlobal2][eIdxGlobal8] += T[0][3];
1016
1017 this->A_[eIdxGlobal4][eIdxGlobal2] -= TSecond[0][0];
1018 this->A_[eIdxGlobal4][eIdxGlobal4] -= TSecond[0][1];
1019 this->A_[eIdxGlobal4][eIdxGlobal3] -= TSecond[0][2];
1020 this->A_[eIdxGlobal4][eIdxGlobal8] -= TSecond[0][3];
1021
1022 u[0] = pc[1];
1023 u[1] = pc[3];
1024 u[2] = pc[2];
1025 u[3] = pc[7];
1026
1027 T.mv(u, Tu);
1028 pcFlux[1][0] = Tu[0];
1029 pcPotential1 = Tu[0];
1030
1031 TSecond.mv(u, Tu);
1032 pcFlux[3][1] = Tu[0];
1033 }
1034 else if (caseL == 3)
1035 {
1036 this->A_[eIdxGlobal2][eIdxGlobal2] += T[0][0];
1037 this->A_[eIdxGlobal2][eIdxGlobal4] += T[0][1];
1038 this->A_[eIdxGlobal2][eIdxGlobal3] += T[0][2];
1039 this->A_[eIdxGlobal2][eIdxGlobal6] += T[0][3];
1040
1041 this->A_[eIdxGlobal4][eIdxGlobal2] -= TSecond[0][0];
1042 this->A_[eIdxGlobal4][eIdxGlobal4] -= TSecond[0][1];
1043 this->A_[eIdxGlobal4][eIdxGlobal3] -= TSecond[0][2];
1044 this->A_[eIdxGlobal4][eIdxGlobal6] -= TSecond[0][3];
1045
1046 u[0] = pc[1];
1047 u[1] = pc[3];
1048 u[2] = pc[2];
1049 u[3] = pc[5];
1050
1051 T.mv(u, Tu);
1052 pcFlux[1][0] = Tu[0];
1053 pcPotential1 = Tu[0];
1054
1055 TSecond.mv(u, Tu);
1056 pcFlux[3][1] = Tu[0];
1057 }
1058 else
1059 {
1060 this->A_[eIdxGlobal2][eIdxGlobal2] += T[0][0];
1061 this->A_[eIdxGlobal2][eIdxGlobal4] += T[0][1];
1062 this->A_[eIdxGlobal2][eIdxGlobal1] += T[0][2];
1063 this->A_[eIdxGlobal2][eIdxGlobal8] += T[0][3];
1064
1065 this->A_[eIdxGlobal4][eIdxGlobal2] -= TSecond[0][0];
1066 this->A_[eIdxGlobal4][eIdxGlobal4] -= TSecond[0][1];
1067 this->A_[eIdxGlobal4][eIdxGlobal1] -= TSecond[0][2];
1068 this->A_[eIdxGlobal4][eIdxGlobal8] -= TSecond[0][3];
1069
1070 u[0] = pc[1];
1071 u[1] = pc[3];
1072 u[2] = pc[0];
1073 u[3] = pc[7];
1074
1075 T.mv(u, Tu);
1076 pcFlux[1][0] = Tu[0];
1077 pcPotential1 = Tu[0];
1078
1079 TSecond.mv(u, Tu);
1080 pcFlux[3][1] = Tu[0];
1081 }
1082
1083 // calculate the flux through the subvolumeface 3 (subVolumeFaceIdx = 2)
1084 caseL = transmissibilityCalculator_.transmissibility(T, interactionVolume, lambda, 3, 2, 1, 0, 7, 6);
1085
1086 TSecond = T;
1087 T *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal4, 3, 0);
1088 TSecond *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal3, 2, 1);
1089
1090 if (caseL == 1)
1091 {
1092 this->A_[eIdxGlobal4][eIdxGlobal4] += T[0][0];
1093 this->A_[eIdxGlobal4][eIdxGlobal3] += T[0][1];
1094 this->A_[eIdxGlobal4][eIdxGlobal2] += T[0][2];
1095 this->A_[eIdxGlobal4][eIdxGlobal8] += T[0][3];
1096
1097 this->A_[eIdxGlobal3][eIdxGlobal4] -= TSecond[0][0];
1098 this->A_[eIdxGlobal3][eIdxGlobal3] -= TSecond[0][1];
1099 this->A_[eIdxGlobal3][eIdxGlobal2] -= TSecond[0][2];
1100 this->A_[eIdxGlobal3][eIdxGlobal8] -= TSecond[0][3];
1101
1102 u[0] = pc[3];
1103 u[1] = pc[2];
1104 u[2] = pc[1];
1105 u[3] = pc[7];
1106
1107 T.mv(u, Tu);
1108 pcPotential2 = Tu[0];
1109 pcFlux[3][0] = Tu[0];
1110
1111 TSecond.mv(u, Tu);
1112 pcFlux[2][1] = Tu[0];
1113 }
1114 else if (caseL == 2)
1115 {
1116 this->A_[eIdxGlobal4][eIdxGlobal4] += T[0][0];
1117 this->A_[eIdxGlobal4][eIdxGlobal3] += T[0][1];
1118 this->A_[eIdxGlobal4][eIdxGlobal1] += T[0][2];
1119 this->A_[eIdxGlobal4][eIdxGlobal7] += T[0][3];
1120
1121 this->A_[eIdxGlobal3][eIdxGlobal4] -= TSecond[0][0];
1122 this->A_[eIdxGlobal3][eIdxGlobal3] -= TSecond[0][1];
1123 this->A_[eIdxGlobal3][eIdxGlobal1] -= TSecond[0][2];
1124 this->A_[eIdxGlobal3][eIdxGlobal7] -= TSecond[0][3];
1125
1126 u[0] = pc[3];
1127 u[1] = pc[2];
1128 u[2] = pc[0];
1129 u[3] = pc[6];
1130
1131 T.mv(u, Tu);
1132 pcPotential2 = Tu[0];
1133 pcFlux[3][0] = Tu[0];
1134
1135 TSecond.mv(u, Tu);
1136 pcFlux[2][1] = Tu[0];
1137 }
1138 else if (caseL == 3)
1139 {
1140 this->A_[eIdxGlobal4][eIdxGlobal4] += T[0][0];
1141 this->A_[eIdxGlobal4][eIdxGlobal3] += T[0][1];
1142 this->A_[eIdxGlobal4][eIdxGlobal1] += T[0][2];
1143 this->A_[eIdxGlobal4][eIdxGlobal8] += T[0][3];
1144
1145 this->A_[eIdxGlobal3][eIdxGlobal4] -= TSecond[0][0];
1146 this->A_[eIdxGlobal3][eIdxGlobal3] -= TSecond[0][1];
1147 this->A_[eIdxGlobal3][eIdxGlobal1] -= TSecond[0][2];
1148 this->A_[eIdxGlobal3][eIdxGlobal8] -= TSecond[0][3];
1149
1150 u[0] = pc[3];
1151 u[1] = pc[2];
1152 u[2] = pc[0];
1153 u[3] = pc[7];
1154
1155 T.mv(u, Tu);
1156 pcPotential2 = Tu[0];
1157 pcFlux[3][0] = Tu[0];
1158
1159 TSecond.mv(u, Tu);
1160 pcFlux[2][1] = Tu[0];
1161 }
1162 else
1163 {
1164 this->A_[eIdxGlobal4][eIdxGlobal4] += T[0][0];
1165 this->A_[eIdxGlobal4][eIdxGlobal3] += T[0][1];
1166 this->A_[eIdxGlobal4][eIdxGlobal2] += T[0][2];
1167 this->A_[eIdxGlobal4][eIdxGlobal7] += T[0][3];
1168
1169 this->A_[eIdxGlobal3][eIdxGlobal4] -= TSecond[0][0];
1170 this->A_[eIdxGlobal3][eIdxGlobal3] -= TSecond[0][1];
1171 this->A_[eIdxGlobal3][eIdxGlobal2] -= TSecond[0][2];
1172 this->A_[eIdxGlobal3][eIdxGlobal7] -= TSecond[0][3];
1173
1174 u[0] = pc[3];
1175 u[1] = pc[2];
1176 u[2] = pc[1];
1177 u[3] = pc[6];
1178
1179 T.mv(u, Tu);
1180 pcPotential2 = Tu[0];
1181 pcFlux[3][0] = Tu[0];
1182
1183 TSecond.mv(u, Tu);
1184 pcFlux[2][1] = Tu[0];
1185 }
1186
1187 // calculate the flux through the subvolumeface 4 (subVolumeFaceIdx = 3)
1188 caseL = transmissibilityCalculator_.transmissibility(T, interactionVolume, lambda, 2, 0, 3, 1, 6, 4);
1189
1190 TSecond = T;
1191 T *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal3, 2, 0);
1192 TSecond *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal1, 0, 1);
1193
1194 if (caseL == 1)
1195 {
1196 this->A_[eIdxGlobal3][eIdxGlobal3] += T[0][0];
1197 this->A_[eIdxGlobal3][eIdxGlobal1] += T[0][1];
1198 this->A_[eIdxGlobal3][eIdxGlobal4] += T[0][2];
1199 this->A_[eIdxGlobal3][eIdxGlobal7] += T[0][3];
1200
1201 this->A_[eIdxGlobal1][eIdxGlobal3] -= TSecond[0][0];
1202 this->A_[eIdxGlobal1][eIdxGlobal1] -= TSecond[0][1];
1203 this->A_[eIdxGlobal1][eIdxGlobal4] -= TSecond[0][2];
1204 this->A_[eIdxGlobal1][eIdxGlobal7] -= TSecond[0][3];
1205
1206 u[0] = pc[2];
1207 u[1] = pc[0];
1208 u[2] = pc[3];
1209 u[3] = pc[6];
1210
1211 T.mv(u, Tu);
1212 pcFlux[2][0] = Tu[0];
1213 pcPotential3 = Tu[0];
1214
1215 TSecond.mv(u, Tu);
1216 pcFlux[0][1] = Tu[0];
1217 }
1218 else if (caseL == 2)
1219 {
1220 this->A_[eIdxGlobal3][eIdxGlobal3] += T[0][0];
1221 this->A_[eIdxGlobal3][eIdxGlobal1] += T[0][1];
1222 this->A_[eIdxGlobal3][eIdxGlobal2] += T[0][2];
1223 this->A_[eIdxGlobal3][eIdxGlobal5] += T[0][3];
1224
1225 this->A_[eIdxGlobal1][eIdxGlobal3] -= TSecond[0][0];
1226 this->A_[eIdxGlobal1][eIdxGlobal1] -= TSecond[0][1];
1227 this->A_[eIdxGlobal1][eIdxGlobal2] -= TSecond[0][2];
1228 this->A_[eIdxGlobal1][eIdxGlobal5] -= TSecond[0][3];
1229
1230 u[0] = pc[2];
1231 u[1] = pc[0];
1232 u[2] = pc[1];
1233 u[3] = pc[4];
1234
1235 T.mv(u, Tu);
1236 pcFlux[2][0] = Tu[0];
1237 pcPotential3 = Tu[0];
1238
1239 TSecond.mv(u, Tu);
1240 pcFlux[0][1] = Tu[0];
1241 }
1242 else if (caseL == 3)
1243 {
1244 this->A_[eIdxGlobal3][eIdxGlobal3] += T[0][0];
1245 this->A_[eIdxGlobal3][eIdxGlobal1] += T[0][1];
1246 this->A_[eIdxGlobal3][eIdxGlobal2] += T[0][2];
1247 this->A_[eIdxGlobal3][eIdxGlobal7] += T[0][3];
1248
1249 this->A_[eIdxGlobal1][eIdxGlobal3] -= TSecond[0][0];
1250 this->A_[eIdxGlobal1][eIdxGlobal1] -= TSecond[0][1];
1251 this->A_[eIdxGlobal1][eIdxGlobal2] -= TSecond[0][2];
1252 this->A_[eIdxGlobal1][eIdxGlobal7] -= TSecond[0][3];
1253
1254 u[0] = pc[2];
1255 u[1] = pc[0];
1256 u[2] = pc[1];
1257 u[3] = pc[6];
1258
1259 T.mv(u, Tu);
1260 pcFlux[2][0] = Tu[0];
1261 pcPotential3 = Tu[0];
1262
1263 TSecond.mv(u, Tu);
1264 pcFlux[0][1] = Tu[0];
1265 }
1266 else
1267 {
1268 this->A_[eIdxGlobal3][eIdxGlobal3] += T[0][0];
1269 this->A_[eIdxGlobal3][eIdxGlobal1] += T[0][1];
1270 this->A_[eIdxGlobal3][eIdxGlobal4] += T[0][2];
1271 this->A_[eIdxGlobal3][eIdxGlobal5] += T[0][3];
1272
1273 this->A_[eIdxGlobal1][eIdxGlobal3] -= TSecond[0][0];
1274 this->A_[eIdxGlobal1][eIdxGlobal1] -= TSecond[0][1];
1275 this->A_[eIdxGlobal1][eIdxGlobal4] -= TSecond[0][2];
1276 this->A_[eIdxGlobal1][eIdxGlobal5] -= TSecond[0][3];
1277
1278 u[0] = pc[2];
1279 u[1] = pc[0];
1280 u[2] = pc[3];
1281 u[3] = pc[4];
1282
1283 T.mv(u, Tu);
1284 pcFlux[2][0] = Tu[0];
1285 pcPotential3 = Tu[0];
1286
1287 TSecond.mv(u, Tu);
1288 pcFlux[0][1] = Tu[0];
1289 }
1290
1291 // calculate the flux through the subvolumeface 5 (subVolumeFaceIdx = 4)
1292 caseL = transmissibilityCalculator_.transmissibility(T, interactionVolume, lambda, 5, 4, 7, 6, 1, 0);
1293
1294 TSecond = T;
1295 T *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal6, 5, 2);
1296 TSecond *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal5, 4, 1);
1297
1298 if (caseL == 1)
1299 {
1300 this->A_[eIdxGlobal6][eIdxGlobal6] += T[0][0];
1301 this->A_[eIdxGlobal6][eIdxGlobal5] += T[0][1];
1302 this->A_[eIdxGlobal6][eIdxGlobal8] += T[0][2];
1303 this->A_[eIdxGlobal6][eIdxGlobal2] += T[0][3];
1304
1305 this->A_[eIdxGlobal5][eIdxGlobal6] -= TSecond[0][0];
1306 this->A_[eIdxGlobal5][eIdxGlobal5] -= TSecond[0][1];
1307 this->A_[eIdxGlobal5][eIdxGlobal8] -= TSecond[0][2];
1308 this->A_[eIdxGlobal5][eIdxGlobal2] -= TSecond[0][3];
1309
1310 u[0] = pc[5];
1311 u[1] = pc[4];
1312 u[2] = pc[7];
1313 u[3] = pc[1];
1314
1315 T.mv(u, Tu);
1316 pcFlux[5][2] = Tu[0];
1317 pcPotential4 = Tu[0];
1318
1319 TSecond.mv(u, Tu);
1320 pcFlux[4][1] = Tu[0];
1321 }
1322 else if (caseL == 2)
1323 {
1324 this->A_[eIdxGlobal6][eIdxGlobal6] += T[0][0];
1325 this->A_[eIdxGlobal6][eIdxGlobal5] += T[0][1];
1326 this->A_[eIdxGlobal6][eIdxGlobal7] += T[0][2];
1327 this->A_[eIdxGlobal6][eIdxGlobal1] += T[0][3];
1328
1329 this->A_[eIdxGlobal5][eIdxGlobal6] -= TSecond[0][0];
1330 this->A_[eIdxGlobal5][eIdxGlobal5] -= TSecond[0][1];
1331 this->A_[eIdxGlobal5][eIdxGlobal7] -= TSecond[0][2];
1332 this->A_[eIdxGlobal5][eIdxGlobal1] -= TSecond[0][3];
1333
1334 u[0] = pc[5];
1335 u[1] = pc[4];
1336 u[2] = pc[6];
1337 u[3] = pc[0];
1338
1339 T.mv(u, Tu);
1340 pcFlux[5][2] = Tu[0];
1341 pcPotential4 = Tu[0];
1342
1343 TSecond.mv(u, Tu);
1344 pcFlux[4][1] = Tu[0];
1345 }
1346 else if (caseL == 3)
1347 {
1348 this->A_[eIdxGlobal6][eIdxGlobal6] += T[0][0];
1349 this->A_[eIdxGlobal6][eIdxGlobal5] += T[0][1];
1350 this->A_[eIdxGlobal6][eIdxGlobal7] += T[0][2];
1351 this->A_[eIdxGlobal6][eIdxGlobal2] += T[0][3];
1352
1353 this->A_[eIdxGlobal5][eIdxGlobal6] -= TSecond[0][0];
1354 this->A_[eIdxGlobal5][eIdxGlobal5] -= TSecond[0][1];
1355 this->A_[eIdxGlobal5][eIdxGlobal7] -= TSecond[0][2];
1356 this->A_[eIdxGlobal5][eIdxGlobal2] -= TSecond[0][3];
1357
1358 u[0] = pc[5];
1359 u[1] = pc[4];
1360 u[2] = pc[6];
1361 u[3] = pc[1];
1362
1363 T.mv(u, Tu);
1364 pcFlux[5][2] = Tu[0];
1365 pcPotential4 = Tu[0];
1366
1367 TSecond.mv(u, Tu);
1368 pcFlux[4][1] = Tu[0];
1369 }
1370 else
1371 {
1372 this->A_[eIdxGlobal6][eIdxGlobal6] += T[0][0];
1373 this->A_[eIdxGlobal6][eIdxGlobal5] += T[0][1];
1374 this->A_[eIdxGlobal6][eIdxGlobal8] += T[0][2];
1375 this->A_[eIdxGlobal6][eIdxGlobal1] += T[0][3];
1376
1377 this->A_[eIdxGlobal5][eIdxGlobal6] -= TSecond[0][0];
1378 this->A_[eIdxGlobal5][eIdxGlobal5] -= TSecond[0][1];
1379 this->A_[eIdxGlobal5][eIdxGlobal8] -= TSecond[0][2];
1380 this->A_[eIdxGlobal5][eIdxGlobal1] -= TSecond[0][3];
1381
1382 u[0] = pc[5];
1383 u[1] = pc[4];
1384 u[2] = pc[7];
1385 u[3] = pc[0];
1386
1387 T.mv(u, Tu);
1388 pcFlux[5][2] = Tu[0];
1389 pcPotential4 = Tu[0];
1390
1391 TSecond.mv(u, Tu);
1392 pcFlux[4][1] = Tu[0];
1393 }
1394
1395 // calculate the flux through the subvolumeface 6 (subVolumeFaceIdx = 5)
1396 caseL = transmissibilityCalculator_.transmissibility(T, interactionVolume, lambda, 7, 5, 6, 4, 3, 1);
1397
1398 TSecond = T;
1399 T *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal8, 7, 2);
1400 TSecond *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal6, 5, 1);
1401
1402
1403 if (caseL == 1)
1404 {
1405 this->A_[eIdxGlobal8][eIdxGlobal8] += T[0][0];
1406 this->A_[eIdxGlobal8][eIdxGlobal6] += T[0][1];
1407 this->A_[eIdxGlobal8][eIdxGlobal7] += T[0][2];
1408 this->A_[eIdxGlobal8][eIdxGlobal4] += T[0][3];
1409
1410 this->A_[eIdxGlobal6][eIdxGlobal8] -= TSecond[0][0];
1411 this->A_[eIdxGlobal6][eIdxGlobal6] -= TSecond[0][1];
1412 this->A_[eIdxGlobal6][eIdxGlobal7] -= TSecond[0][2];
1413 this->A_[eIdxGlobal6][eIdxGlobal4] -= TSecond[0][3];
1414
1415 u[0] = pc[7];
1416 u[1] = pc[5];
1417 u[2] = pc[6];
1418 u[3] = pc[3];
1419
1420 T.mv(u, Tu);
1421 pcFlux[7][2] = Tu[0];
1422 pcPotential5 = Tu[0];
1423
1424 TSecond.mv(u, Tu);
1425 pcFlux[5][1] = Tu[0];
1426 }
1427 else if (caseL == 2)
1428 {
1429 this->A_[eIdxGlobal8][eIdxGlobal8] += T[0][0];
1430 this->A_[eIdxGlobal8][eIdxGlobal6] += T[0][1];
1431 this->A_[eIdxGlobal8][eIdxGlobal5] += T[0][2];
1432 this->A_[eIdxGlobal8][eIdxGlobal2] += T[0][3];
1433
1434 this->A_[eIdxGlobal6][eIdxGlobal8] -= TSecond[0][0];
1435 this->A_[eIdxGlobal6][eIdxGlobal6] -= TSecond[0][1];
1436 this->A_[eIdxGlobal6][eIdxGlobal5] -= TSecond[0][2];
1437 this->A_[eIdxGlobal6][eIdxGlobal2] -= TSecond[0][3];
1438
1439 u[0] = pc[7];
1440 u[1] = pc[5];
1441 u[2] = pc[4];
1442 u[3] = pc[1];
1443
1444 T.mv(u, Tu);
1445 pcFlux[7][2] = Tu[0];
1446 pcPotential5 = Tu[0];
1447
1448 TSecond.mv(u, Tu);
1449 pcFlux[5][1] = Tu[0];
1450 }
1451 else if (caseL == 3)
1452 {
1453 this->A_[eIdxGlobal8][eIdxGlobal8] += T[0][0];
1454 this->A_[eIdxGlobal8][eIdxGlobal6] += T[0][1];
1455 this->A_[eIdxGlobal8][eIdxGlobal5] += T[0][2];
1456 this->A_[eIdxGlobal8][eIdxGlobal4] += T[0][3];
1457
1458 this->A_[eIdxGlobal6][eIdxGlobal8] -= TSecond[0][0];
1459 this->A_[eIdxGlobal6][eIdxGlobal6] -= TSecond[0][1];
1460 this->A_[eIdxGlobal6][eIdxGlobal5] -= TSecond[0][2];
1461 this->A_[eIdxGlobal6][eIdxGlobal4] -= TSecond[0][3];
1462
1463 u[0] = pc[7];
1464 u[1] = pc[5];
1465 u[2] = pc[4];
1466 u[3] = pc[3];
1467
1468 T.mv(u, Tu);
1469 pcFlux[7][2] = Tu[0];
1470 pcPotential5 = Tu[0];
1471
1472 TSecond.mv(u, Tu);
1473 pcFlux[5][1] = Tu[0];
1474 }
1475 else
1476 {
1477 this->A_[eIdxGlobal8][eIdxGlobal8] += T[0][0];
1478 this->A_[eIdxGlobal8][eIdxGlobal6] += T[0][1];
1479 this->A_[eIdxGlobal8][eIdxGlobal7] += T[0][2];
1480 this->A_[eIdxGlobal8][eIdxGlobal2] += T[0][3];
1481
1482 this->A_[eIdxGlobal6][eIdxGlobal8] -= TSecond[0][0];
1483 this->A_[eIdxGlobal6][eIdxGlobal6] -= TSecond[0][1];
1484 this->A_[eIdxGlobal6][eIdxGlobal7] -= TSecond[0][2];
1485 this->A_[eIdxGlobal6][eIdxGlobal2] -= TSecond[0][3];
1486
1487 u[0] = pc[7];
1488 u[1] = pc[5];
1489 u[2] = pc[6];
1490 u[3] = pc[1];
1491
1492 T.mv(u, Tu);
1493 pcFlux[7][2] = Tu[0];
1494 pcPotential5 = Tu[0];
1495
1496 TSecond.mv(u, Tu);
1497 pcFlux[5][1] = Tu[0];
1498 }
1499
1500 // calculate the flux through the subvolumeface 7 (subVolumeFaceIdx = 6)
1501 caseL = transmissibilityCalculator_.transmissibility(T, interactionVolume, lambda, 6, 7, 4, 5, 2, 3);
1502
1503 TSecond = T;
1504 T *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal7, 6, 2);
1505 TSecond *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal8, 7, 1);
1506
1507 if (caseL == 1)
1508 {
1509 this->A_[eIdxGlobal7][eIdxGlobal7] += T[0][0];
1510 this->A_[eIdxGlobal7][eIdxGlobal8] += T[0][1];
1511 this->A_[eIdxGlobal7][eIdxGlobal5] += T[0][2];
1512 this->A_[eIdxGlobal7][eIdxGlobal3] += T[0][3];
1513
1514 this->A_[eIdxGlobal8][eIdxGlobal7] -= TSecond[0][0];
1515 this->A_[eIdxGlobal8][eIdxGlobal8] -= TSecond[0][1];
1516 this->A_[eIdxGlobal8][eIdxGlobal5] -= TSecond[0][2];
1517 this->A_[eIdxGlobal8][eIdxGlobal3] -= TSecond[0][3];
1518
1519 u[0] = pc[6];
1520 u[1] = pc[7];
1521 u[2] = pc[4];
1522 u[3] = pc[2];
1523
1524 T.mv(u, Tu);
1525 pcFlux[6][2] = Tu[0];
1526 pcPotential6 = Tu[0];
1527
1528 TSecond.mv(u, Tu);
1529 pcFlux[7][1] = Tu[0];
1530 }
1531 else if (caseL == 2)
1532 {
1533 this->A_[eIdxGlobal7][eIdxGlobal7] += T[0][0];
1534 this->A_[eIdxGlobal7][eIdxGlobal8] += T[0][1];
1535 this->A_[eIdxGlobal7][eIdxGlobal6] += T[0][2];
1536 this->A_[eIdxGlobal7][eIdxGlobal4] += T[0][3];
1537
1538 this->A_[eIdxGlobal8][eIdxGlobal7] -= TSecond[0][0];
1539 this->A_[eIdxGlobal8][eIdxGlobal8] -= TSecond[0][1];
1540 this->A_[eIdxGlobal8][eIdxGlobal6] -= TSecond[0][2];
1541 this->A_[eIdxGlobal8][eIdxGlobal4] -= TSecond[0][3];
1542
1543 u[0] = pc[6];
1544 u[1] = pc[7];
1545 u[2] = pc[5];
1546 u[3] = pc[3];
1547
1548 T.mv(u, Tu);
1549 pcFlux[6][2] = Tu[0];
1550 pcPotential6 = Tu[0];
1551
1552 TSecond.mv(u, Tu);
1553 pcFlux[7][1] = Tu[0];
1554 }
1555 else if (caseL == 3)
1556 {
1557 this->A_[eIdxGlobal7][eIdxGlobal7] += T[0][0];
1558 this->A_[eIdxGlobal7][eIdxGlobal8] += T[0][1];
1559 this->A_[eIdxGlobal7][eIdxGlobal6] += T[0][2];
1560 this->A_[eIdxGlobal7][eIdxGlobal3] += T[0][3];
1561
1562 this->A_[eIdxGlobal8][eIdxGlobal7] -= TSecond[0][0];
1563 this->A_[eIdxGlobal8][eIdxGlobal8] -= TSecond[0][1];
1564 this->A_[eIdxGlobal8][eIdxGlobal6] -= TSecond[0][2];
1565 this->A_[eIdxGlobal8][eIdxGlobal3] -= TSecond[0][3];
1566
1567 u[0] = pc[6];
1568 u[1] = pc[7];
1569 u[2] = pc[5];
1570 u[3] = pc[2];
1571
1572 T.mv(u, Tu);
1573 pcFlux[6][2] = Tu[0];
1574 pcPotential6 = Tu[0];
1575
1576 TSecond.mv(u, Tu);
1577 pcFlux[7][1] = Tu[0];
1578 }
1579 else
1580 {
1581 this->A_[eIdxGlobal7][eIdxGlobal7] += T[0][0];
1582 this->A_[eIdxGlobal7][eIdxGlobal8] += T[0][1];
1583 this->A_[eIdxGlobal7][eIdxGlobal5] += T[0][2];
1584 this->A_[eIdxGlobal7][eIdxGlobal4] += T[0][3];
1585
1586 this->A_[eIdxGlobal8][eIdxGlobal7] -= TSecond[0][0];
1587 this->A_[eIdxGlobal8][eIdxGlobal8] -= TSecond[0][1];
1588 this->A_[eIdxGlobal8][eIdxGlobal5] -= TSecond[0][2];
1589 this->A_[eIdxGlobal8][eIdxGlobal4] -= TSecond[0][3];
1590
1591 u[0] = pc[6];
1592 u[1] = pc[7];
1593 u[2] = pc[4];
1594 u[3] = pc[3];
1595
1596 T.mv(u, Tu);
1597 pcFlux[6][2] = Tu[0];
1598 pcPotential6 = Tu[0];
1599
1600 TSecond.mv(u, Tu);
1601 pcFlux[7][1] = Tu[0];
1602 }
1603
1604 // calculate the flux through the subvolumeface 8 (subVolumeFaceIdx = 7)
1605 caseL = transmissibilityCalculator_.transmissibility(T, interactionVolume, lambda, 4, 6, 5, 7, 0, 2);
1606
1607 TSecond = T;
1608 T *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal5, 4, 2);
1609 TSecond *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal7, 6, 1);
1610
1611 if (caseL == 1)
1612 {
1613 this->A_[eIdxGlobal5][eIdxGlobal5] += T[0][0];
1614 this->A_[eIdxGlobal5][eIdxGlobal7] += T[0][1];
1615 this->A_[eIdxGlobal5][eIdxGlobal6] += T[0][2];
1616 this->A_[eIdxGlobal5][eIdxGlobal1] += T[0][3];
1617
1618 this->A_[eIdxGlobal7][eIdxGlobal5] -= TSecond[0][0];
1619 this->A_[eIdxGlobal7][eIdxGlobal7] -= TSecond[0][1];
1620 this->A_[eIdxGlobal7][eIdxGlobal6] -= TSecond[0][2];
1621 this->A_[eIdxGlobal7][eIdxGlobal1] -= TSecond[0][3];
1622
1623 u[0] = pc[4];
1624 u[1] = pc[6];
1625 u[2] = pc[5];
1626 u[3] = pc[0];
1627
1628 T.mv(u, Tu);
1629 pcFlux[4][2] = Tu[0];
1630 pcPotential7 = Tu[0];
1631
1632 TSecond.mv(u, Tu);
1633 pcFlux[6][1] = Tu[0];
1634 }
1635 else if (caseL == 2)
1636 {
1637 this->A_[eIdxGlobal5][eIdxGlobal5] += T[0][0];
1638 this->A_[eIdxGlobal5][eIdxGlobal7] += T[0][1];
1639 this->A_[eIdxGlobal5][eIdxGlobal8] += T[0][2];
1640 this->A_[eIdxGlobal5][eIdxGlobal3] += T[0][3];
1641
1642 this->A_[eIdxGlobal7][eIdxGlobal5] -= TSecond[0][0];
1643 this->A_[eIdxGlobal7][eIdxGlobal7] -= TSecond[0][1];
1644 this->A_[eIdxGlobal7][eIdxGlobal8] -= TSecond[0][2];
1645 this->A_[eIdxGlobal7][eIdxGlobal3] -= TSecond[0][3];
1646
1647 u[0] = pc[4];
1648 u[1] = pc[6];
1649 u[2] = pc[7];
1650 u[3] = pc[2];
1651
1652 T.mv(u, Tu);
1653 pcFlux[4][2] = Tu[0];
1654 pcPotential7 = Tu[0];
1655
1656 TSecond.mv(u, Tu);
1657 pcFlux[6][1] = Tu[0];
1658 }
1659 else if (caseL == 3)
1660 {
1661 this->A_[eIdxGlobal5][eIdxGlobal5] += T[0][0];
1662 this->A_[eIdxGlobal5][eIdxGlobal7] += T[0][1];
1663 this->A_[eIdxGlobal5][eIdxGlobal8] += T[0][2];
1664 this->A_[eIdxGlobal5][eIdxGlobal1] += T[0][3];
1665
1666 this->A_[eIdxGlobal7][eIdxGlobal5] -= TSecond[0][0];
1667 this->A_[eIdxGlobal7][eIdxGlobal7] -= TSecond[0][1];
1668 this->A_[eIdxGlobal7][eIdxGlobal8] -= TSecond[0][2];
1669 this->A_[eIdxGlobal7][eIdxGlobal1] -= TSecond[0][3];
1670
1671 u[0] = pc[4];
1672 u[1] = pc[6];
1673 u[2] = pc[7];
1674 u[3] = pc[0];
1675
1676 T.mv(u, Tu);
1677 pcFlux[4][2] = Tu[0];
1678 pcPotential7 = Tu[0];
1679
1680 TSecond.mv(u, Tu);
1681 pcFlux[6][1] = Tu[0];
1682 }
1683 else
1684 {
1685 this->A_[eIdxGlobal5][eIdxGlobal5] += T[0][0];
1686 this->A_[eIdxGlobal5][eIdxGlobal7] += T[0][1];
1687 this->A_[eIdxGlobal5][eIdxGlobal6] += T[0][2];
1688 this->A_[eIdxGlobal5][eIdxGlobal3] += T[0][3];
1689
1690 this->A_[eIdxGlobal7][eIdxGlobal5] -= TSecond[0][0];
1691 this->A_[eIdxGlobal7][eIdxGlobal7] -= TSecond[0][1];
1692 this->A_[eIdxGlobal7][eIdxGlobal6] -= TSecond[0][2];
1693 this->A_[eIdxGlobal7][eIdxGlobal3] -= TSecond[0][3];
1694
1695 u[0] = pc[4];
1696 u[1] = pc[6];
1697 u[2] = pc[5];
1698 u[3] = pc[2];
1699
1700 T.mv(u, Tu);
1701 pcFlux[4][2] = Tu[0];
1702 pcPotential7 = Tu[0];
1703
1704 TSecond.mv(u, Tu);
1705 pcFlux[6][1] = Tu[0];
1706 }
1707
1708 // calculate the flux through the subvolumeface 9 (subVolumeFaceIdx = 8)
1709 caseL = transmissibilityCalculator_.transmissibility(T, interactionVolume, lambda, 4, 0, 6, 2, 5, 1);
1710
1711 TSecond = T;
1712 T *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal5, 4, 0);
1713 TSecond *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal1, 0, 2);
1714
1715 if (caseL == 1)
1716 {
1717 this->A_[eIdxGlobal5][eIdxGlobal5] += T[0][0];
1718 this->A_[eIdxGlobal5][eIdxGlobal1] += T[0][1];
1719 this->A_[eIdxGlobal5][eIdxGlobal7] += T[0][2];
1720 this->A_[eIdxGlobal5][eIdxGlobal6] += T[0][3];
1721
1722 this->A_[eIdxGlobal1][eIdxGlobal5] -= TSecond[0][0];
1723 this->A_[eIdxGlobal1][eIdxGlobal1] -= TSecond[0][1];
1724 this->A_[eIdxGlobal1][eIdxGlobal7] -= TSecond[0][2];
1725 this->A_[eIdxGlobal1][eIdxGlobal6] -= TSecond[0][3];
1726
1727 u[0] = pc[4];
1728 u[1] = pc[0];
1729 u[2] = pc[6];
1730 u[3] = pc[5];
1731
1732 T.mv(u, Tu);
1733 pcFlux[4][0] = Tu[0];
1734 pcPotential8 = Tu[0];
1735
1736 TSecond.mv(u, Tu);
1737 pcFlux[0][2] = Tu[0];
1738 }
1739 else if (caseL == 2)
1740 {
1741 this->A_[eIdxGlobal5][eIdxGlobal5] += T[0][0];
1742 this->A_[eIdxGlobal5][eIdxGlobal1] += T[0][1];
1743 this->A_[eIdxGlobal5][eIdxGlobal3] += T[0][2];
1744 this->A_[eIdxGlobal5][eIdxGlobal2] += T[0][3];
1745
1746 this->A_[eIdxGlobal1][eIdxGlobal5] -= TSecond[0][0];
1747 this->A_[eIdxGlobal1][eIdxGlobal1] -= TSecond[0][1];
1748 this->A_[eIdxGlobal1][eIdxGlobal3] -= TSecond[0][2];
1749 this->A_[eIdxGlobal1][eIdxGlobal2] -= TSecond[0][3];
1750
1751 u[0] = pc[4];
1752 u[1] = pc[0];
1753 u[2] = pc[2];
1754 u[3] = pc[1];
1755
1756 T.mv(u, Tu);
1757 pcFlux[4][0] = Tu[0];
1758 pcPotential8 = Tu[0];
1759
1760 TSecond.mv(u, Tu);
1761 pcFlux[0][2] = Tu[0];
1762 }
1763 else if (caseL == 3)
1764 {
1765 this->A_[eIdxGlobal5][eIdxGlobal5] += T[0][0];
1766 this->A_[eIdxGlobal5][eIdxGlobal1] += T[0][1];
1767 this->A_[eIdxGlobal5][eIdxGlobal3] += T[0][2];
1768 this->A_[eIdxGlobal5][eIdxGlobal6] += T[0][3];
1769
1770 this->A_[eIdxGlobal1][eIdxGlobal5] -= TSecond[0][0];
1771 this->A_[eIdxGlobal1][eIdxGlobal1] -= TSecond[0][1];
1772 this->A_[eIdxGlobal1][eIdxGlobal3] -= TSecond[0][2];
1773 this->A_[eIdxGlobal1][eIdxGlobal6] -= TSecond[0][3];
1774
1775 u[0] = pc[4];
1776 u[1] = pc[0];
1777 u[2] = pc[2];
1778 u[3] = pc[5];
1779
1780 T.mv(u, Tu);
1781 pcFlux[4][0] = Tu[0];
1782 pcPotential8 = Tu[0];
1783
1784 TSecond.mv(u, Tu);
1785 pcFlux[0][2] = Tu[0];
1786 }
1787 else
1788 {
1789 this->A_[eIdxGlobal5][eIdxGlobal5] += T[0][0];
1790 this->A_[eIdxGlobal5][eIdxGlobal1] += T[0][1];
1791 this->A_[eIdxGlobal5][eIdxGlobal7] += T[0][2];
1792 this->A_[eIdxGlobal5][eIdxGlobal2] += T[0][3];
1793
1794 this->A_[eIdxGlobal1][eIdxGlobal5] -= TSecond[0][0];
1795 this->A_[eIdxGlobal1][eIdxGlobal1] -= TSecond[0][1];
1796 this->A_[eIdxGlobal1][eIdxGlobal7] -= TSecond[0][2];
1797 this->A_[eIdxGlobal1][eIdxGlobal2] -= TSecond[0][3];
1798
1799 u[0] = pc[4];
1800 u[1] = pc[0];
1801 u[2] = pc[6];
1802 u[3] = pc[1];
1803
1804 T.mv(u, Tu);
1805 pcFlux[4][0] = Tu[0];
1806 pcPotential8 = Tu[0];
1807
1808 TSecond.mv(u, Tu);
1809 pcFlux[0][2] = Tu[0];
1810 }
1811
1812 // calculate the flux through the subvolumeface 10 (subVolumeFaceIdx = 9)
1813 caseL = transmissibilityCalculator_.transmissibility(T, interactionVolume, lambda, 1, 5, 3, 7, 0, 4);
1814
1815 TSecond = T;
1816 T *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal2, 1, 2);
1817 TSecond *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal6, 5, 0);
1818
1819 if (caseL == 1)
1820 {
1821 this->A_[eIdxGlobal2][eIdxGlobal2] += T[0][0];
1822 this->A_[eIdxGlobal2][eIdxGlobal6] += T[0][1];
1823 this->A_[eIdxGlobal2][eIdxGlobal4] += T[0][2];
1824 this->A_[eIdxGlobal2][eIdxGlobal1] += T[0][3];
1825
1826 this->A_[eIdxGlobal6][eIdxGlobal2] -= TSecond[0][0];
1827 this->A_[eIdxGlobal6][eIdxGlobal6] -= TSecond[0][1];
1828 this->A_[eIdxGlobal6][eIdxGlobal4] -= TSecond[0][2];
1829 this->A_[eIdxGlobal6][eIdxGlobal1] -= TSecond[0][3];
1830
1831 u[0] = pc[1];
1832 u[1] = pc[5];
1833 u[2] = pc[3];
1834 u[3] = pc[0];
1835
1836 T.mv(u, Tu);
1837 pcFlux[1][2] = Tu[0];
1838 pcPotential9 = Tu[0];
1839
1840 TSecond.mv(u, Tu);
1841 pcFlux[5][0] = Tu[0];
1842 }
1843 else if (caseL == 2)
1844 {
1845 this->A_[eIdxGlobal2][eIdxGlobal2] += T[0][0];
1846 this->A_[eIdxGlobal2][eIdxGlobal6] += T[0][1];
1847 this->A_[eIdxGlobal2][eIdxGlobal8] += T[0][2];
1848 this->A_[eIdxGlobal2][eIdxGlobal5] += T[0][3];
1849
1850 this->A_[eIdxGlobal6][eIdxGlobal2] -= TSecond[0][0];
1851 this->A_[eIdxGlobal6][eIdxGlobal6] -= TSecond[0][1];
1852 this->A_[eIdxGlobal6][eIdxGlobal8] -= TSecond[0][2];
1853 this->A_[eIdxGlobal6][eIdxGlobal5] -= TSecond[0][3];
1854
1855 u[0] = pc[1];
1856 u[1] = pc[5];
1857 u[2] = pc[7];
1858 u[3] = pc[4];
1859
1860 T.mv(u, Tu);
1861 pcFlux[1][2] = Tu[0];
1862 pcPotential9 = Tu[0];
1863
1864 TSecond.mv(u, Tu);
1865 pcFlux[5][0] = Tu[0];
1866 }
1867 else if (caseL == 3)
1868 {
1869 this->A_[eIdxGlobal2][eIdxGlobal2] += T[0][0];
1870 this->A_[eIdxGlobal2][eIdxGlobal6] += T[0][1];
1871 this->A_[eIdxGlobal2][eIdxGlobal8] += T[0][2];
1872 this->A_[eIdxGlobal2][eIdxGlobal1] += T[0][3];
1873
1874 this->A_[eIdxGlobal6][eIdxGlobal2] -= TSecond[0][0];
1875 this->A_[eIdxGlobal6][eIdxGlobal6] -= TSecond[0][1];
1876 this->A_[eIdxGlobal6][eIdxGlobal8] -= TSecond[0][2];
1877 this->A_[eIdxGlobal6][eIdxGlobal1] -= TSecond[0][3];
1878
1879 u[0] = pc[1];
1880 u[1] = pc[5];
1881 u[2] = pc[7];
1882 u[3] = pc[0];
1883
1884 T.mv(u, Tu);
1885 pcFlux[1][2] = Tu[0];
1886 pcPotential9 = Tu[0];
1887
1888 TSecond.mv(u, Tu);
1889 pcFlux[5][0] = Tu[0];
1890 }
1891 else
1892 {
1893 this->A_[eIdxGlobal2][eIdxGlobal2] += T[0][0];
1894 this->A_[eIdxGlobal2][eIdxGlobal6] += T[0][1];
1895 this->A_[eIdxGlobal2][eIdxGlobal4] += T[0][2];
1896 this->A_[eIdxGlobal2][eIdxGlobal5] += T[0][3];
1897
1898 this->A_[eIdxGlobal6][eIdxGlobal2] -= TSecond[0][0];
1899 this->A_[eIdxGlobal6][eIdxGlobal6] -= TSecond[0][1];
1900 this->A_[eIdxGlobal6][eIdxGlobal4] -= TSecond[0][2];
1901 this->A_[eIdxGlobal6][eIdxGlobal5] -= TSecond[0][3];
1902
1903 u[0] = pc[1];
1904 u[1] = pc[5];
1905 u[2] = pc[3];
1906 u[3] = pc[4];
1907
1908 T.mv(u, Tu);
1909 pcFlux[1][2] = Tu[0];
1910 pcPotential9 = Tu[0];
1911
1912 TSecond.mv(u, Tu);
1913 pcFlux[5][0] = Tu[0];
1914 }
1915
1916 // calculate the flux through the subvolumeface 11 (subVolumeFaceIdx = 10)
1917 caseL = transmissibilityCalculator_.transmissibility(T, interactionVolume, lambda, 7, 3, 5, 1, 6, 2);
1918
1919 TSecond = T;
1920 T *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal8, 7, 0);
1921 TSecond *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal4, 3, 2);
1922
1923 if (caseL == 1)
1924 {
1925 this->A_[eIdxGlobal8][eIdxGlobal8] += T[0][0];
1926 this->A_[eIdxGlobal8][eIdxGlobal4] += T[0][1];
1927 this->A_[eIdxGlobal8][eIdxGlobal6] += T[0][2];
1928 this->A_[eIdxGlobal8][eIdxGlobal7] += T[0][3];
1929
1930 this->A_[eIdxGlobal4][eIdxGlobal8] -= TSecond[0][0];
1931 this->A_[eIdxGlobal4][eIdxGlobal4] -= TSecond[0][1];
1932 this->A_[eIdxGlobal4][eIdxGlobal6] -= TSecond[0][2];
1933 this->A_[eIdxGlobal4][eIdxGlobal7] -= TSecond[0][3];
1934
1935 u[0] = pc[7];
1936 u[1] = pc[3];
1937 u[2] = pc[5];
1938 u[3] = pc[6];
1939
1940 T.mv(u, Tu);
1941 pcFlux[7][0] = Tu[0];
1942 pcPotential10 = Tu[0];
1943
1944 TSecond.mv(u, Tu);
1945 pcFlux[3][2] = Tu[0];
1946 }
1947 else if (caseL == 2)
1948 {
1949 this->A_[eIdxGlobal8][eIdxGlobal8] += T[0][0];
1950 this->A_[eIdxGlobal8][eIdxGlobal4] += T[0][1];
1951 this->A_[eIdxGlobal8][eIdxGlobal2] += T[0][2];
1952 this->A_[eIdxGlobal8][eIdxGlobal3] += T[0][3];
1953
1954 this->A_[eIdxGlobal4][eIdxGlobal8] -= TSecond[0][0];
1955 this->A_[eIdxGlobal4][eIdxGlobal4] -= TSecond[0][1];
1956 this->A_[eIdxGlobal4][eIdxGlobal2] -= TSecond[0][2];
1957 this->A_[eIdxGlobal4][eIdxGlobal3] -= TSecond[0][3];
1958
1959 u[0] = pc[7];
1960 u[1] = pc[3];
1961 u[2] = pc[1];
1962 u[3] = pc[2];
1963
1964 T.mv(u, Tu);
1965 pcFlux[7][0] = Tu[0];
1966 pcPotential10 = Tu[0];
1967
1968 TSecond.mv(u, Tu);
1969 pcFlux[3][2] = Tu[0];
1970 }
1971 else if (caseL == 3)
1972 {
1973 this->A_[eIdxGlobal8][eIdxGlobal8] += T[0][0];
1974 this->A_[eIdxGlobal8][eIdxGlobal4] += T[0][1];
1975 this->A_[eIdxGlobal8][eIdxGlobal2] += T[0][2];
1976 this->A_[eIdxGlobal8][eIdxGlobal7] += T[0][3];
1977
1978 this->A_[eIdxGlobal4][eIdxGlobal8] -= TSecond[0][0];
1979 this->A_[eIdxGlobal4][eIdxGlobal4] -= TSecond[0][1];
1980 this->A_[eIdxGlobal4][eIdxGlobal2] -= TSecond[0][2];
1981 this->A_[eIdxGlobal4][eIdxGlobal7] -= TSecond[0][3];
1982
1983 u[0] = pc[7];
1984 u[1] = pc[3];
1985 u[2] = pc[1];
1986 u[3] = pc[6];
1987
1988 T.mv(u, Tu);
1989 pcFlux[7][0] = Tu[0];
1990 pcPotential10 = Tu[0];
1991
1992 TSecond.mv(u, Tu);
1993 pcFlux[3][2] = Tu[0];
1994 }
1995 else
1996 {
1997 this->A_[eIdxGlobal8][eIdxGlobal8] += T[0][0];
1998 this->A_[eIdxGlobal8][eIdxGlobal4] += T[0][1];
1999 this->A_[eIdxGlobal8][eIdxGlobal6] += T[0][2];
2000 this->A_[eIdxGlobal8][eIdxGlobal3] += T[0][3];
2001
2002 this->A_[eIdxGlobal4][eIdxGlobal8] -= TSecond[0][0];
2003 this->A_[eIdxGlobal4][eIdxGlobal4] -= TSecond[0][1];
2004 this->A_[eIdxGlobal4][eIdxGlobal6] -= TSecond[0][2];
2005 this->A_[eIdxGlobal4][eIdxGlobal3] -= TSecond[0][3];
2006
2007 u[0] = pc[7];
2008 u[1] = pc[3];
2009 u[2] = pc[5];
2010 u[3] = pc[2];
2011
2012 T.mv(u, Tu);
2013 pcFlux[7][0] = Tu[0];
2014 pcPotential10 = Tu[0];
2015
2016 TSecond.mv(u, Tu);
2017 pcFlux[3][2] = Tu[0];
2018 }
2019
2020 // calculate the flux through the subvolumeface 12 (subVolumeFaceIdx = 11)
2021 caseL = transmissibilityCalculator_.transmissibility(T, interactionVolume, lambda, 2, 6, 0, 4, 3, 7);
2022
2023 TSecond = T;
2024 T *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal3, 2, 2);
2025 TSecond *= interactionVolumes_.faceAreaFactor(interactionVolume, eIdxGlobal7, 6, 0);
2026
2027 if (caseL == 1)
2028 {
2029 this->A_[eIdxGlobal3][eIdxGlobal3] += T[0][0];
2030 this->A_[eIdxGlobal3][eIdxGlobal7] += T[0][1];
2031 this->A_[eIdxGlobal3][eIdxGlobal1] += T[0][2];
2032 this->A_[eIdxGlobal3][eIdxGlobal4] += T[0][3];
2033
2034 this->A_[eIdxGlobal7][eIdxGlobal3] -= TSecond[0][0];
2035 this->A_[eIdxGlobal7][eIdxGlobal7] -= TSecond[0][1];
2036 this->A_[eIdxGlobal7][eIdxGlobal1] -= TSecond[0][2];
2037 this->A_[eIdxGlobal7][eIdxGlobal4] -= TSecond[0][3];
2038
2039 u[0] = pc[2];
2040 u[1] = pc[6];
2041 u[2] = pc[0];
2042 u[3] = pc[3];
2043
2044 T.mv(u, Tu);
2045 pcFlux[2][2] = Tu[0];
2046 pcPotential11 = Tu[0];
2047
2048 TSecond.mv(u, Tu);
2049 pcFlux[6][0] = Tu[0];
2050 }
2051 else if (caseL == 2)
2052 {
2053 this->A_[eIdxGlobal3][eIdxGlobal3] += T[0][0];
2054 this->A_[eIdxGlobal3][eIdxGlobal7] += T[0][1];
2055 this->A_[eIdxGlobal3][eIdxGlobal5] += T[0][2];
2056 this->A_[eIdxGlobal3][eIdxGlobal8] += T[0][3];
2057
2058 this->A_[eIdxGlobal7][eIdxGlobal3] -= TSecond[0][0];
2059 this->A_[eIdxGlobal7][eIdxGlobal7] -= TSecond[0][1];
2060 this->A_[eIdxGlobal7][eIdxGlobal5] -= TSecond[0][2];
2061 this->A_[eIdxGlobal7][eIdxGlobal8] -= TSecond[0][3];
2062
2063 u[0] = pc[2];
2064 u[1] = pc[6];
2065 u[2] = pc[4];
2066 u[3] = pc[7];
2067
2068 T.mv(u, Tu);
2069 pcFlux[2][2] = Tu[0];
2070 pcPotential11 = Tu[0];
2071
2072 TSecond.mv(u, Tu);
2073 pcFlux[6][0] = Tu[0];
2074 }
2075 else if (caseL == 3)
2076 {
2077 this->A_[eIdxGlobal3][eIdxGlobal3] += T[0][0];
2078 this->A_[eIdxGlobal3][eIdxGlobal7] += T[0][1];
2079 this->A_[eIdxGlobal3][eIdxGlobal5] += T[0][2];
2080 this->A_[eIdxGlobal3][eIdxGlobal4] += T[0][3];
2081
2082 this->A_[eIdxGlobal7][eIdxGlobal3] -= TSecond[0][0];
2083 this->A_[eIdxGlobal7][eIdxGlobal7] -= TSecond[0][1];
2084 this->A_[eIdxGlobal7][eIdxGlobal5] -= TSecond[0][2];
2085 this->A_[eIdxGlobal7][eIdxGlobal4] -= TSecond[0][3];
2086
2087 u[0] = pc[2];
2088 u[1] = pc[6];
2089 u[2] = pc[4];
2090 u[3] = pc[3];
2091
2092 T.mv(u, Tu);
2093 pcFlux[2][2] = Tu[0];
2094 pcPotential11 = Tu[0];
2095
2096 TSecond.mv(u, Tu);
2097 pcFlux[6][0] = Tu[0];
2098 }
2099 else
2100 {
2101 this->A_[eIdxGlobal3][eIdxGlobal3] += T[0][0];
2102 this->A_[eIdxGlobal3][eIdxGlobal7] += T[0][1];
2103 this->A_[eIdxGlobal3][eIdxGlobal1] += T[0][2];
2104 this->A_[eIdxGlobal3][eIdxGlobal8] += T[0][3];
2105
2106 this->A_[eIdxGlobal7][eIdxGlobal3] -= TSecond[0][0];
2107 this->A_[eIdxGlobal7][eIdxGlobal7] -= TSecond[0][1];
2108 this->A_[eIdxGlobal7][eIdxGlobal1] -= TSecond[0][2];
2109 this->A_[eIdxGlobal7][eIdxGlobal8] -= TSecond[0][3];
2110
2111 u[0] = pc[2];
2112 u[1] = pc[6];
2113 u[2] = pc[0];
2114 u[3] = pc[7];
2115
2116 T.mv(u, Tu);
2117 pcFlux[2][2] = Tu[0];
2118 pcPotential11 = Tu[0];
2119
2120 TSecond.mv(u, Tu);
2121 pcFlux[6][0] = Tu[0];
2122 }
2123
2124 if (pc[0] == 0 && pc[1] == 0 && pc[2] == 0 && pc[3] == 0 && pc[4] == 0 && pc[5] == 0 && pc[6] == 0
2125 && pc[7] == 0)
2126 {
2127 return;
2128 }
2129
2130 // compute mobilities of subvolumeface 1 (subVolumeFaceIdx = 0)
2131 Dune::FieldVector<Scalar, numPhases> lambda0Upw(0.0);
2132 lambda0Upw[wPhaseIdx] = (pcPotential0 >= 0) ? lambda1[wPhaseIdx] : lambda2[wPhaseIdx];
2133 lambda0Upw[nPhaseIdx] = (pcPotential0 >= 0) ? lambda1[nPhaseIdx] : lambda2[nPhaseIdx];
2134
2135 // compute mobilities of subvolumeface 2 (subVolumeFaceIdx = 1)
2136 Dune::FieldVector<Scalar, numPhases> lambda1Upw(0.0);
2137 lambda1Upw[wPhaseIdx] = (pcPotential1 >= 0) ? lambda2[wPhaseIdx] : lambda4[wPhaseIdx];
2138 lambda1Upw[nPhaseIdx] = (pcPotential1 >= 0) ? lambda2[nPhaseIdx] : lambda4[nPhaseIdx];
2139
2140 // compute mobilities of subvolumeface 3 (subVolumeFaceIdx = 2)
2141 Dune::FieldVector<Scalar, numPhases> lambda2Upw(0.0);
2142 lambda2Upw[wPhaseIdx] = (pcPotential2 >= 0) ? lambda4[wPhaseIdx] : lambda3[wPhaseIdx];
2143 lambda2Upw[nPhaseIdx] = (pcPotential2 >= 0) ? lambda4[nPhaseIdx] : lambda3[nPhaseIdx];
2144
2145 // compute mobilities of subvolumeface 4 (subVolumeFaceIdx = 3)
2146 Dune::FieldVector<Scalar, numPhases> lambda3Upw(0.0);
2147 lambda3Upw[wPhaseIdx] = (pcPotential3 >= 0) ? lambda3[wPhaseIdx] : lambda1[wPhaseIdx];
2148 lambda3Upw[nPhaseIdx] = (pcPotential3 >= 0) ? lambda3[nPhaseIdx] : lambda1[nPhaseIdx];
2149
2150 // compute mobilities of subvolumeface 5 (subVolumeFaceIdx = 4)
2151 Dune::FieldVector<Scalar, numPhases> lambda4Upw(0.0);
2152 lambda4Upw[wPhaseIdx] = (pcPotential4 >= 0) ? lambda6[wPhaseIdx] : lambda5[wPhaseIdx];
2153 lambda4Upw[nPhaseIdx] = (pcPotential4 >= 0) ? lambda6[nPhaseIdx] : lambda5[nPhaseIdx];
2154
2155 // compute mobilities of subvolumeface 6 (subVolumeFaceIdx = 5)
2156 Dune::FieldVector<Scalar, numPhases> lambda5Upw(0.0);
2157 lambda5Upw[wPhaseIdx] = (pcPotential5 >= 0) ? lambda8[wPhaseIdx] : lambda6[wPhaseIdx];
2158 lambda5Upw[nPhaseIdx] = (pcPotential5 >= 0) ? lambda8[nPhaseIdx] : lambda6[nPhaseIdx];
2159
2160 // compute mobilities of subvolumeface 7 (subVolumeFaceIdx = 6)
2161 Dune::FieldVector<Scalar, numPhases> lambda6Upw(0.0);
2162 lambda6Upw[wPhaseIdx] = (pcPotential6 >= 0) ? lambda7[wPhaseIdx] : lambda8[wPhaseIdx];
2163 lambda6Upw[nPhaseIdx] = (pcPotential6 >= 0) ? lambda7[nPhaseIdx] : lambda8[nPhaseIdx];
2164
2165 // compute mobilities of subvolumeface 8 (subVolumeFaceIdx = 7)
2166 Dune::FieldVector<Scalar, numPhases> lambda7Upw(0.0);
2167 lambda7Upw[wPhaseIdx] = (pcPotential7 >= 0) ? lambda5[wPhaseIdx] : lambda7[wPhaseIdx];
2168 lambda7Upw[nPhaseIdx] = (pcPotential7 >= 0) ? lambda5[nPhaseIdx] : lambda7[nPhaseIdx];
2169
2170 // compute mobilities of subvolumeface 9 (subVolumeFaceIdx = 8)
2171 Dune::FieldVector<Scalar, numPhases> lambda8Upw(0.0);
2172 lambda8Upw[wPhaseIdx] = (pcPotential8 >= 0) ? lambda5[wPhaseIdx] : lambda1[wPhaseIdx];
2173 lambda8Upw[nPhaseIdx] = (pcPotential8 >= 0) ? lambda5[nPhaseIdx] : lambda1[nPhaseIdx];
2174
2175 // compute mobilities of subvolumeface 10 (subVolumeFaceIdx = 9)
2176 Dune::FieldVector<Scalar, numPhases> lambda9Upw(0.0);
2177 lambda9Upw[wPhaseIdx] = (pcPotential9 >= 0) ? lambda2[wPhaseIdx] : lambda6[wPhaseIdx];
2178 lambda9Upw[nPhaseIdx] = (pcPotential9 >= 0) ? lambda2[nPhaseIdx] : lambda6[nPhaseIdx];
2179
2180 // compute mobilities of subvolumeface 11 (subVolumeFaceIdx = 10)
2181 Dune::FieldVector<Scalar, numPhases> lambda10Upw(0.0);
2182 lambda10Upw[wPhaseIdx] = (pcPotential10 >= 0) ? lambda8[wPhaseIdx] : lambda4[wPhaseIdx];
2183 lambda10Upw[nPhaseIdx] = (pcPotential10 >= 0) ? lambda8[nPhaseIdx] : lambda4[nPhaseIdx];
2184
2185 // compute mobilities of subvolumeface 12 (subVolumeFaceIdx = 11)
2186 Dune::FieldVector<Scalar, numPhases> lambda11Upw(0.0);
2187 lambda11Upw[wPhaseIdx] = (pcPotential11 >= 0) ? lambda3[wPhaseIdx] : lambda7[wPhaseIdx];
2188 lambda11Upw[nPhaseIdx] = (pcPotential11 >= 0) ? lambda3[nPhaseIdx] : lambda7[nPhaseIdx];
2189
2190 for (int i = 0; i < numPhases; i++)
2191 {
2192 Scalar lambdaT0 = lambda0Upw[wPhaseIdx] + lambda0Upw[nPhaseIdx];
2193 Scalar lambdaT1 = lambda1Upw[wPhaseIdx] + lambda1Upw[nPhaseIdx];
2194 Scalar lambdaT2 = lambda2Upw[wPhaseIdx] + lambda2Upw[nPhaseIdx];
2195 Scalar lambdaT3 = lambda3Upw[wPhaseIdx] + lambda3Upw[nPhaseIdx];
2196 Scalar lambdaT4 = lambda4Upw[wPhaseIdx] + lambda4Upw[nPhaseIdx];
2197 Scalar lambdaT5 = lambda5Upw[wPhaseIdx] + lambda5Upw[nPhaseIdx];
2198 Scalar lambdaT6 = lambda6Upw[wPhaseIdx] + lambda6Upw[nPhaseIdx];
2199 Scalar lambdaT7 = lambda7Upw[wPhaseIdx] + lambda7Upw[nPhaseIdx];
2200 Scalar lambdaT8 = lambda8Upw[wPhaseIdx] + lambda8Upw[nPhaseIdx];
2201 Scalar lambdaT9 = lambda9Upw[wPhaseIdx] + lambda9Upw[nPhaseIdx];
2202 Scalar lambdaT10 = lambda10Upw[wPhaseIdx] + lambda10Upw[nPhaseIdx];
2203 Scalar lambdaT11 = lambda11Upw[wPhaseIdx] + lambda11Upw[nPhaseIdx];
2204 Scalar fracFlow0 = (lambdaT0 > threshold_) ? lambda0Upw[i] / (lambdaT0) : 0.0;
2205 Scalar fracFlow1 = (lambdaT1 > threshold_) ? lambda1Upw[i] / (lambdaT1) : 0.0;
2206 Scalar fracFlow2 = (lambdaT2 > threshold_) ? lambda2Upw[i] / (lambdaT2) : 0.0;
2207 Scalar fracFlow3 = (lambdaT3 > threshold_) ? lambda3Upw[i] / (lambdaT3) : 0.0;
2208 Scalar fracFlow4 = (lambdaT4 > threshold_) ? lambda4Upw[i] / (lambdaT4) : 0.0;
2209 Scalar fracFlow5 = (lambdaT5 > threshold_) ? lambda5Upw[i] / (lambdaT5) : 0.0;
2210 Scalar fracFlow6 = (lambdaT6 > threshold_) ? lambda6Upw[i] / (lambdaT6) : 0.0;
2211 Scalar fracFlow7 = (lambdaT7 > threshold_) ? lambda7Upw[i] / (lambdaT7) : 0.0;
2212 Scalar fracFlow8 = (lambdaT8 > threshold_) ? lambda8Upw[i] / (lambdaT8) : 0.0;
2213 Scalar fracFlow9 = (lambdaT9 > threshold_) ? lambda9Upw[i] / (lambdaT9) : 0.0;
2214 Scalar fracFlow10 = (lambdaT10 > threshold_) ? lambda10Upw[i] / (lambdaT10) : 0.0;
2215 Scalar fracFlow11 = (lambdaT11 > threshold_) ? lambda11Upw[i] / (lambdaT11) : 0.0;
2216
2217 switch (pressureType_)
2218 {
2219 case pw:
2220 {
2221 if (i == nPhaseIdx)
2222 {
2223 // add capillary pressure term to right hand side
2224 this->f_[eIdxGlobal1] -= (fracFlow0 * pcFlux[0][0] - fracFlow3 * pcFlux[0][1] - fracFlow8 * pcFlux[0][2]);
2225 this->f_[eIdxGlobal2] -= (fracFlow1 * pcFlux[1][0] - fracFlow0 * pcFlux[1][1] + fracFlow9 * pcFlux[1][2]);
2226 this->f_[eIdxGlobal3] -= (fracFlow3 * pcFlux[2][0] - fracFlow2 * pcFlux[2][1] + fracFlow11 * pcFlux[2][2]);
2227 this->f_[eIdxGlobal4] -= (fracFlow2 * pcFlux[3][0] - fracFlow1 * pcFlux[3][1] - fracFlow10 * pcFlux[3][2]);
2228 this->f_[eIdxGlobal5] -= (fracFlow8 * pcFlux[4][0] - fracFlow4 * pcFlux[4][1] + fracFlow7 * pcFlux[4][2]);
2229 this->f_[eIdxGlobal6] -= (-fracFlow9 * pcFlux[5][0] - fracFlow5 * pcFlux[5][1] + fracFlow4 * pcFlux[5][2]);
2230 this->f_[eIdxGlobal7] -= (-fracFlow11 * pcFlux[6][0] - fracFlow7 * pcFlux[6][1] + fracFlow6 * pcFlux[6][2]);
2231 this->f_[eIdxGlobal8] -= (fracFlow10 * pcFlux[7][0] - fracFlow6 * pcFlux[7][1] + fracFlow5 * pcFlux[7][2]);
2232 }
2233 break;
2234 }
2235 case pn:
2236 {
2237 if (i == wPhaseIdx)
2238 {
2239 // add capillary pressure term to right hand side
2240 this->f_[eIdxGlobal1] += (fracFlow0 * pcFlux[0][0] - fracFlow3 * pcFlux[0][1] - fracFlow8 * pcFlux[0][2]);
2241 this->f_[eIdxGlobal2] += (fracFlow1 * pcFlux[1][0] - fracFlow0 * pcFlux[1][1] + fracFlow9 * pcFlux[1][2]);
2242 this->f_[eIdxGlobal3] += (fracFlow3 * pcFlux[2][0] - fracFlow2 * pcFlux[2][1] + fracFlow11 * pcFlux[2][2]);
2243 this->f_[eIdxGlobal4] += (fracFlow2 * pcFlux[3][0] - fracFlow1 * pcFlux[3][1] - fracFlow10 * pcFlux[3][2]);
2244 this->f_[eIdxGlobal5] += (fracFlow8 * pcFlux[4][0] - fracFlow4 * pcFlux[4][1] + fracFlow7 * pcFlux[4][2]);
2245 this->f_[eIdxGlobal6] += (-fracFlow9 * pcFlux[5][0] - fracFlow5 * pcFlux[5][1] + fracFlow4 * pcFlux[5][2]);
2246 this->f_[eIdxGlobal7] += (-fracFlow11 * pcFlux[6][0] - fracFlow7 * pcFlux[6][1] + fracFlow6 * pcFlux[6][2]);
2247 this->f_[eIdxGlobal8] += (fracFlow10 * pcFlux[7][0] - fracFlow6 * pcFlux[7][1] + fracFlow5 * pcFlux[7][2]);
2248 }
2249 break;
2250 }
2251 }
2252 }
2253}
2254
2255// TODO doc me!
2256template<class TypeTag>
2257void FvMpfaL3dPressure2p<TypeTag>::assembleBoundaryInteractionVolume(InteractionVolume& interactionVolume)
2258{
2259 for (int elemIdx = 0; elemIdx < 8; elemIdx++)
2260 {
2261 if (!interactionVolume.hasSubVolumeElement(elemIdx))
2262 {
2263 continue;
2264 }
2265 bool isOutside = false;
2266 for (int fIdx = 0; fIdx < dim; fIdx++)
2267 {
2268 int intVolFaceIdx = interactionVolume.getFaceIndexFromSubVolume(elemIdx, fIdx);
2269 if (interactionVolume.isOutsideFace(intVolFaceIdx))
2270 {
2271 isOutside = true;
2272 break;
2273 }
2274 }
2275 if (isOutside)
2276 {
2277 continue;
2278 }
2279
2280 auto element = interactionVolume.getSubVolumeElement(elemIdx);
2281
2282 // get global coordinate of cell centers
2283 const GlobalPosition& globalPos = element.geometry().center();
2284
2285 // cell volumes
2286 Scalar volume = element.geometry().volume();
2287
2288 // cell index
2289 int eIdxGlobal = problem_.variables().index(element);
2290
2291 //get the cell Data
2292 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
2293
2294 // permeability vector at boundary
2295 DimMatrix permeability(problem_.spatialParams().intrinsicPermeability(element));
2296
2297 // evaluate right hand side
2298 PrimaryVariables source(0);
2299 problem_.source(source, element);
2300 this->f_[eIdxGlobal] += volume / (8.0)
2301 * (source[wPhaseIdx] / density_[wPhaseIdx] + source[nPhaseIdx] / density_[nPhaseIdx]);
2302
2303 this->f_[eIdxGlobal] += evaluateErrorTerm(cellData) * volume / (8.0);
2304
2305 //get mobilities of the phases
2306 Dune::FieldVector<Scalar, numPhases> lambda(cellData.mobility(wPhaseIdx));
2307 lambda[nPhaseIdx] = cellData.mobility(nPhaseIdx);
2308
2309 Scalar pc = cellData.capillaryPressure();
2310
2311 Scalar gravityDiff = (problem_.bBoxMax() - globalPos) * gravity_
2312 * (density_[nPhaseIdx] - density_[wPhaseIdx]);
2313
2314 pc += gravityDiff; //minus because of gravity definition!
2315
2316 for (int fIdx = 0; fIdx < dim; fIdx++)
2317 {
2318 int intVolFaceIdx = interactionVolume.getFaceIndexFromSubVolume(elemIdx, fIdx);
2319
2320 if (interactionVolume.isBoundaryFace(intVolFaceIdx))
2321 {
2322 if (interactionVolume.getBoundaryType(intVolFaceIdx).isDirichlet(pressureEqIdx))
2323 {
2324 const GlobalPosition& globalPosFace = interactionVolume.getFacePosition(elemIdx, fIdx);
2325
2326 DimVector distVec(globalPosFace - globalPos);
2327 Scalar dist = distVec.two_norm();
2328 DimVector& normal = interactionVolume.getNormal(elemIdx, fIdx);
2329
2330 Scalar faceArea = interactionVolume.getFaceArea(elemIdx, fIdx);
2331
2332 // get pc and lambda at the boundary
2333 Scalar satWBound = cellData.saturation(wPhaseIdx);
2334 //check boundary sat at face 1
2335 if (interactionVolume.getBoundaryType(intVolFaceIdx).isDirichlet(satEqIdx))
2336 {
2337 Scalar satBound = interactionVolume.getDirichletValues(intVolFaceIdx)[saturationIdx];
2338 switch (saturationType_)
2339 {
2340 case sw:
2341 {
2342 satWBound = satBound;
2343 break;
2344 }
2345 case sn:
2346 {
2347 satWBound = 1 - satBound;
2348 break;
2349 }
2350 }
2351
2352 }
2353
2354 Scalar pcBound = MaterialLaw::pc(
2355 problem_.spatialParams().materialLawParams(element), satWBound);
2356
2357 Scalar gravityDiffBound = (problem_.bBoxMax() - globalPosFace) * gravity_
2358 * (density_[nPhaseIdx] - density_[wPhaseIdx]);
2359
2360 pcBound += gravityDiffBound;
2361
2362 Dune::FieldVector<Scalar, numPhases>
2363 lambdaBound(MaterialLaw::krw(problem_.spatialParams().materialLawParams(element),satWBound));
2364 lambdaBound[nPhaseIdx] = MaterialLaw::krn(
2365 problem_.spatialParams().materialLawParams(element), satWBound);
2366 lambdaBound[wPhaseIdx] /= viscosity_[wPhaseIdx];
2367 lambdaBound[nPhaseIdx] /= viscosity_[nPhaseIdx];
2368
2369 Scalar potentialBound = interactionVolume.getDirichletValues(intVolFaceIdx)[pressureIdx];
2370 Scalar gdeltaZ = (problem_.bBoxMax()-globalPosFace) * gravity_;
2371
2372 //calculate potential gradients
2373 Scalar potentialDiffW = 0;
2374 Scalar potentialDiffNw = 0;
2375 switch (pressureType_)
2376 {
2377 case pw:
2378 {
2379 potentialBound += density_[wPhaseIdx]*gdeltaZ;
2380 potentialDiffW = (cellData.potential(wPhaseIdx) - potentialBound) / dist;
2381 potentialDiffNw = (cellData.potential(nPhaseIdx) - potentialBound - pcBound) / dist;
2382 break;
2383 }
2384 case pn:
2385 {
2386 potentialBound += density_[nPhaseIdx]*gdeltaZ;
2387 potentialDiffW = (cellData.potential(wPhaseIdx) - potentialBound + pcBound) / dist;
2388 potentialDiffNw = (cellData.potential(nPhaseIdx) - potentialBound) / dist;
2389 break;
2390 }
2391 }
2392
2393 Scalar lambdaTotal = (potentialDiffW >= 0.) ? lambda[wPhaseIdx] : lambdaBound[wPhaseIdx];
2394 lambdaTotal += (potentialDiffNw >= 0.) ? lambda[nPhaseIdx] : lambdaBound[nPhaseIdx];
2395
2396 DimVector permTimesNormal(0);
2397 permeability.mv(normal, permTimesNormal);
2398 Scalar scalarPerm = permTimesNormal.two_norm();
2399
2400 //calculate current matrix entry
2401 Scalar entry = lambdaTotal * scalarPerm / dist * faceArea;
2402
2403 // set diagonal entry and right hand side entry
2404 this->A_[eIdxGlobal][eIdxGlobal] += entry;
2405 this->f_[eIdxGlobal] += entry * potentialBound;
2406
2407 if (pc == 0 && pcBound == 0)
2408 {
2409 continue;
2410 }
2411
2412 //calculate right hand side
2413 Scalar pcFlux = 0;
2414
2415 switch (pressureType_)
2416 {
2417 case pw:
2418 {
2419 //add capillary pressure term to right hand side
2420 pcFlux = 0.5 * (lambda[nPhaseIdx] + lambdaBound[nPhaseIdx])
2421 * scalarPerm * (pc - pcBound) / dist * faceArea;
2422
2423 break;
2424 }
2425 case pn:
2426 {
2427 //add capillary pressure term to right hand side
2428 pcFlux = 0.5 * (lambda[wPhaseIdx] + lambdaBound[wPhaseIdx])
2429 * scalarPerm * (pc - pcBound) / dist * faceArea;
2430
2431 break;
2432 }
2433 }
2434
2435 for (int i = 0; i < numPhases; i++)
2436 {
2437 switch (pressureType_)
2438 {
2439 case pw:
2440 {
2441 if (i == nPhaseIdx)
2442 {
2443 //add capillary pressure term to right hand side
2444 this->f_[eIdxGlobal] -= pcFlux;
2445 }
2446 break;
2447 }
2448 case pn:
2449 {
2450 if (i == wPhaseIdx)
2451 {
2452 //add capillary pressure term to right hand side
2453 this->f_[eIdxGlobal] += pcFlux;
2454 }
2455 break;
2456 }
2457 }
2458 }
2459 }
2460 else if (interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressureEqIdx))
2461 {
2462 Scalar J = interactionVolume.getNeumannValues(intVolFaceIdx)[wPhaseIdx]
2463 / density_[wPhaseIdx];
2464 J += interactionVolume.getNeumannValues(intVolFaceIdx)[nPhaseIdx] / density_[nPhaseIdx];
2465 this->f_[eIdxGlobal] -= J;
2466 }
2467 else
2468 {
2469 std::cout << "interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressureEqIdx)"
2470 << interactionVolume.getBoundaryType(intVolFaceIdx).isNeumann(pressureEqIdx) << "\n";
2471 DUNE_THROW(Dune::NotImplemented,
2472 "No valid boundary condition type defined for pressure equation!");
2473 }
2474 }
2475 }
2476 }
2477}
2478
2480template<class TypeTag>
2481void FvMpfaL3dPressure2p<TypeTag>::updateMaterialLaws()
2482{
2483 // iterate through leaf grid an evaluate c0 at cell center
2484 for (const auto& element : elements(problem_.gridView()))
2485 {
2486 int eIdxGlobal = problem_.variables().index(element);
2487
2488 CellData& cellData = problem_.variables().cellData(eIdxGlobal);
2489
2490 Scalar satW = cellData.saturation(wPhaseIdx);
2491
2492 Scalar pc = MaterialLaw::pc(problem_.spatialParams().materialLawParams(element), satW);
2493
2494 cellData.setCapillaryPressure(pc);
2495
2496 // initialize mobilities
2497 Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(element), satW)
2498 / viscosity_[wPhaseIdx];
2499 Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(element), satW)
2500 / viscosity_[nPhaseIdx];
2501
2502 // initialize mobilities
2503 cellData.setMobility(wPhaseIdx, mobilityW);
2504 cellData.setMobility(nPhaseIdx, mobilityNw);
2505
2506 //initialize fractional flow functions
2507 cellData.setFracFlowFunc(wPhaseIdx, mobilityW / (mobilityW + mobilityNw));
2508 cellData.setFracFlowFunc(nPhaseIdx, mobilityNw / (mobilityW + mobilityNw));
2509 }
2510 return;
2511}
2512
2513} // end namespace Dumux
2514#endif
3dinteractionvolumecontainer.hh
Interactionvolume container for 3-d MPFA L-method.
3dtransmissibilitycalculator.hh
Provides methods for transmissibility calculation in 3-d.
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::FvMpfaL3dPressure2p
3-d finite volume MPFA L-method discretization of a two-phase flow pressure equation of the sequentia...
Definition: 3dpressure.hh:74
Dumux::FvMpfaL3dPressure2p::evaluateErrorTerm
Scalar evaluateErrorTerm(CellData &cellData)
Volume correction term to correct for unphysical saturation overshoots/undershoots.
Definition: 3dpressure.hh:341
Dumux::FvMpfaL3dPressure2p::initializeMatrixRowSize
void initializeMatrixRowSize()
Initializes the row size of the sparse matrix for the pressure solution.
Definition: 3dpressure.hh:555
Dumux::FvMpfaL3dPressure2p::assembleInnerInteractionVolume
void assembleInnerInteractionVolume(InteractionVolume &interactionVolume)
assembles the matrix entries of one interaction volume into the global matrix
Definition: 3dpressure.hh:658
Dumux::FvMpfaL3dPressure2p::TransmissibilityType
typename TransmissibilityCalculator::TransmissibilityType TransmissibilityType
Type including methods for calculation of MPFA transmissibilities.
Definition: 3dpressure.hh:153
Dumux::FvMpfaL3dPressure2p::assemble
void assemble()
function which assembles the system of equations to be solved
Definition: 3dpressure.hh:626
Dumux::FvMpfaL3dPressure2p::update
void update()
Pressure update.
Definition: 3dpressure.hh:272
Dumux::FvMpfaL3dPressure2p::interactionVolumes_
InteractionVolumeContainer interactionVolumes_
Definition: 3dpressure.hh:508
Dumux::FvMpfaL3dPressure2p::InteractionVolume
GetPropType< TypeTag, Properties::MPFAInteractionVolume > InteractionVolume
Type for storing interaction volume information.
Definition: 3dpressure.hh:155
Dumux::FvMpfaL3dPressure2p::assembleBoundaryInteractionVolume
void assembleBoundaryInteractionVolume(InteractionVolume &interactionVolume)
Definition: 3dpressure.hh:2257
Dumux::FvMpfaL3dPressure2p::transmissibilityCalculator_
TransmissibilityCalculator transmissibilityCalculator_
The transmissibility calculator including methods for the MPFA transmissibility calculation.
Definition: 3dpressure.hh:510
Dumux::FvMpfaL3dPressure2p::initialize
void initialize(bool solveTwice=true)
Initializes the pressure model.
Definition: 3dpressure.hh:178
Dumux::FvMpfaL3dPressure2p::interactionVolumes
InteractionVolumeContainer & interactionVolumes()
Returns the global container of the stored interaction volumes.
Definition: 3dpressure.hh:450
Dumux::FvMpfaL3dPressure2p::transmissibilityCalculator
TransmissibilityCalculator & transmissibilityCalculator()
Returns the transmissibility calculator.
Definition: 3dpressure.hh:460
Dumux::FvMpfaL3dPressure2p::initializeMatrix
void initializeMatrix()
Initializes the sparse matrix for the pressure solution.
Definition: 3dpressure.hh:545
Dumux::FvMpfaL3dPressure2p::FvMpfaL3dPressure2p
FvMpfaL3dPressure2p(Problem &problem)
Constructs a FvMpfaL3dPressure2p object.
Definition: 3dpressure.hh:469
Dumux::FvMpfaL3dPressure2p::storePressureSolution
void storePressureSolution()
Globally stores the pressure solution.
Definition: 3dpressure.hh:206
Dumux::FvMpfaL3dPressure2p::addOutputVtkFields
void addOutputVtkFields(MultiWriter &writer)
Adds pressure output to the output file.
Definition: 3dpressure.hh:385
Dumux::FvMpfaL3dPressure2p::storePressureSolution
void storePressureSolution(const Element &element)
Stores the pressure solution of a cell.
Definition: 3dpressure.hh:220
Dumux::FvMpfaL3dPressure2p::updateMaterialLaws
void updateMaterialLaws()
constitutive functions are initialized and stored in the variables object
Definition: 3dpressure.hh:2481
Dumux::FvMpfaL3dPressure2p::initializeMatrixIndices
void initializeMatrixIndices()
Initializes the indices of the sparse matrix for the pressure solution.
Definition: 3dpressure.hh:592
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::FVPressure
The finite volume base class for the solution of a pressure equation.
Definition: sequential/cellcentered/pressure.hh:49
Dumux::FVPressure::initialize
void initialize()
Initialize pressure model.
Definition: sequential/cellcentered/pressure.hh:213
Dumux::FVPressure::pressure
PressureSolution & pressure()
Returns the vector containing the pressure solution.
Definition: sequential/cellcentered/pressure.hh:120
Dumux::FVPressure::solve
void solve()
Solves the global system of equations to get the spatial distribution of the pressure.
Definition: sequential/cellcentered/pressure.hh:527
properties.hh
Specifies the properties for immiscible 2p diffusion/pressure models.
pressure.hh
Finite Volume Diffusion Model.
Generated by doxygen 1.9.3