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