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