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