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