3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
2dtransmissibilitycalculator.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2/*****************************************************************************
3 * See the file COPYING for full copying permissions. *
4 * *
5 * This program is free software: you can redistribute it and/or modify *
6 * it under the terms of the GNU General Public License as published by *
7 * the Free Software Foundation, either version 3 of the License, or *
8 * (at your option) any later version. *
9 * *
10 * This program is distributed in the hope that it will be useful, *
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13 * GNU General Public License for more details. *
14 * *
15 * You should have received a copy of the GNU General Public License *
16 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
17 *****************************************************************************/
23#ifndef DUMUX_FVMPFAL2D_TRANSMISSIBILITYCALCULATOR_HH
24#define DUMUX_FVMPFAL2D_TRANSMISSIBILITYCALCULATOR_HH
25
26// dumux environment
27#include <dumux/porousmediumflow/sequential/pressureproperties.hh>
28#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/properties.hh>
29#include <dumux/porousmediumflow/sequential/cellcentered/mpfa/linteractionvolume.hh>
30
31namespace Dumux {
32
42template<class TypeTag>
43class FvMpfaL2dTransmissibilityCalculator
44{
45 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
46
47 enum
48 {
49 dim = GridView::dimension, dimWorld = GridView::dimensionworld
50 };
51
52 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
53 using Problem = GetPropType<TypeTag, Properties::Problem>;
54
55 using Element = typename GridView::template Codim<0>::Entity;
56
57 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
58 using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
59
60 using DimVector = Dune::FieldVector<Scalar, dim>;
61
62 using InteractionVolume = FVMPFALInteractionVolume<TypeTag>;
63
64public:
65 using TransmissibilityType = Dune::FieldMatrix<Scalar, dim, 2*dim - dim + 1>;
66
68 enum
69 {
70 leftTriangle = -1,
71 noTransmissibility = 0,
72 rightTriangle = 1
73 };
74
76 int calculateTransmissibility(
77 TransmissibilityType& transmissibility,
78 InteractionVolume& interactionVolume,
79 std::vector<DimVector >& lambda,
80 int idx1, int idx2, int idx3, int idx4);
81
83 int calculateLeftHNTransmissibility(TransmissibilityType& transmissibility,
84 InteractionVolume& interactionVolume,
85 std::vector<DimVector >& lambda,
86 int idx1, int idx2, int idx3);
87
89 int calculateRightHNTransmissibility(TransmissibilityType& transmissibility,
90 InteractionVolume& interactionVolume,
91 std::vector<DimVector >& lambda,
92 int idx1, int idx2, int idx3);
93
99 FvMpfaL2dTransmissibilityCalculator(Problem& problem) :
100 problem_(problem), R_(0)
101 {
102 if (dim != 2)
103 {
104 DUNE_THROW(Dune::NotImplemented, "Dimension not supported!");
105 }
106
107 // evaluate matrix R
108 if (dim == 2)
109 {
110 R_[0][1] = 1;
111 R_[1][0] = -1;
112 }
113 }
114
115private:
116 Problem& problem_;
117 DimMatrix R_;
118};
119
135template<class TypeTag>
136int FvMpfaL2dTransmissibilityCalculator<TypeTag>::calculateTransmissibility(
137 TransmissibilityType& transmissibility,
138 InteractionVolume& interactionVolume,
139 std::vector<DimVector >& lambda,
140 int idx1, int idx2, int idx3, int idx4)
141{
142 auto element1 = interactionVolume.getSubVolumeElement(idx1);
143 auto element2 = interactionVolume.getSubVolumeElement(idx2);
144 auto element3 = interactionVolume.getSubVolumeElement(idx3);
145 auto element4 = interactionVolume.getSubVolumeElement(idx4);
146
147 if (element3 == element4 && element1.level() != element2.level())
148 {
149 return noTransmissibility;
150 }
151
152 // get global coordinate of cell centers
153 const GlobalPosition& globalPos1 = element1.geometry().center();
154 const GlobalPosition& globalPos2 = element2.geometry().center();
155 const GlobalPosition& globalPos3 = element3.geometry().center();
156 const GlobalPosition& globalPos4 = element4.geometry().center();
157
158 const GlobalPosition& globalPosCenter = interactionVolume.getCenterPosition();
159
160 const DimMatrix& K1 = problem_.spatialParams().intrinsicPermeability(element1);
161 const DimMatrix& K2 = problem_.spatialParams().intrinsicPermeability(element2);
162 const DimMatrix& K3 = problem_.spatialParams().intrinsicPermeability(element3);
163 const DimMatrix& K4 = problem_.spatialParams().intrinsicPermeability(element4);
164
165 const GlobalPosition& globalPosFace12 = interactionVolume.getFacePosition(idx1, 0);
166 const GlobalPosition& globalPosFace23 = interactionVolume.getFacePosition(idx2, 0);
167
168 // compute normal vectors nu1-nu7 in triangle R for first half edge
169 DimVector nu1R1(0);
170 R_.mv(globalPosFace12 - globalPos2, nu1R1);
171
172 DimVector nu2R1(0);
173 R_.mv(globalPos2 - globalPosFace23, nu2R1);
174
175 DimVector nu3R1(0);
176 R_.mv(globalPosFace23 - globalPos3, nu3R1);
177
178 DimVector nu4R1(0);
179 R_.mv(globalPos3 - globalPosCenter, nu4R1);
180
181 DimVector nu5R1(0);
182 R_.mv(globalPosCenter - globalPos1, nu5R1);
183
184 DimVector nu6R1(0);
185 R_.mv(globalPos1 - globalPosFace12, nu6R1);
186
187 DimVector nu7R1(0);
188 R_.mv(globalPosCenter - globalPos2, nu7R1);
189
190 // compute T, i.e., the area of quadrilateral made by normal vectors 'nu'
191 DimVector Rnu2R1(0);
192 R_.mv(nu2R1, Rnu2R1);
193 Scalar T1R1 = nu1R1 * Rnu2R1;
194
195 DimVector Rnu4R1(0);
196 R_.mv(nu4R1, Rnu4R1);
197 Scalar T2R1 = nu3R1 * Rnu4R1;
198
199 DimVector Rnu6R1(0);
200 R_.mv(nu6R1, Rnu6R1);
201 Scalar T3R1 = nu5R1 * Rnu6R1;
202
203 // compute components needed for flux calculation, denoted as 'omega' and 'chi'
204 DimVector K2nu1R1(0);
205 K2.mv(nu1R1, K2nu1R1);
206 DimVector K2nu2R1(0);
207 K2.mv(nu2R1, K2nu2R1);
208 DimVector K4nu3R1(0);
209 K3.mv(nu3R1, K4nu3R1);
210 DimVector K4nu4R1(0);
211 K3.mv(nu4R1, K4nu4R1);
212 DimVector K1nu5R1(0);
213 K1.mv(nu5R1, K1nu5R1);
214 DimVector K1nu6R1(0);
215 K1.mv(nu6R1, K1nu6R1);
216
217 DimVector Rnu1R1(0);
218 R_.mv(nu1R1, Rnu1R1);
219
220 DimVector &outerNormaln1R1 = interactionVolume.getNormal(idx2, 0);
221 DimVector &outerNormaln2 = interactionVolume.getNormal(idx1, 0);
222
223 Scalar omega111R1 = lambda[idx2][0] * (outerNormaln1R1 * K2nu1R1) * interactionVolume.getFaceArea(idx2, 0) / T1R1;
224 Scalar omega112R1 = lambda[idx2][0] * (outerNormaln1R1 * K2nu2R1) * interactionVolume.getFaceArea(idx2, 0) / T1R1;
225 Scalar omega211R1 = lambda[idx2][1] * (outerNormaln2 * K2nu1R1) * interactionVolume.getFaceArea(idx2, 1) / T1R1;
226 Scalar omega212R1 = lambda[idx2][1] * (outerNormaln2 * K2nu2R1) * interactionVolume.getFaceArea(idx2, 1) / T1R1;
227 Scalar omega123R1 = lambda[idx3][1] * (outerNormaln1R1 * K4nu3R1) * interactionVolume.getFaceArea(idx3, 1) / T2R1;
228 Scalar omega124R1 = lambda[idx3][1] * (outerNormaln1R1 * K4nu4R1) * interactionVolume.getFaceArea(idx3, 1) / T2R1;
229 Scalar omega235R1 = lambda[idx1][0] * (outerNormaln2 * K1nu5R1) * interactionVolume.getFaceArea(idx1, 0) / T3R1;
230 Scalar omega236R1 = lambda[idx1][0] * (outerNormaln2 * K1nu6R1) * interactionVolume.getFaceArea(idx1, 0) / T3R1;
231 Scalar chi711R1 = (nu7R1 * Rnu1R1) / T1R1;
232 Scalar chi712R1 = (nu7R1 * Rnu2R1) / T1R1;
233
234 // compute transmissibility matrix TR1 = CA^{-1}B+D
235 DimMatrix C(0), A(0);
236 Dune::FieldMatrix<Scalar, dim, 2 * dim - dim + 1> D(0), B(0);
237
238 // evaluate matrix C, D, A, B
239 C[0][0] = -omega111R1;
240 C[0][1] = -omega112R1;
241 C[1][0] = -omega211R1;
242 C[1][1] = -omega212R1;
243
244 D[0][0] = omega111R1 + omega112R1;
245 D[1][0] = omega211R1 + omega212R1;
246
247 A[0][0] = omega111R1 - omega124R1 - omega123R1 * chi711R1;
248 A[0][1] = omega112R1 - omega123R1 * chi712R1;
249 A[1][0] = omega211R1 - omega236R1 * chi711R1;
250 A[1][1] = omega212R1 - omega235R1 - omega236R1 * chi712R1;
251
252 B[0][0] = omega111R1 + omega112R1 + omega123R1 * (1.0 - chi711R1 - chi712R1);
253 B[0][1] = -omega123R1 - omega124R1;
254 B[1][0] = omega211R1 + omega212R1 + omega236R1 * (1.0 - chi711R1 - chi712R1);
255 B[1][2] = -omega235R1 - omega236R1;
256
257 // compute TR1
258 A.invert();
259 D += B.leftmultiply(C.rightmultiply(A));
260 Dune::FieldMatrix<Scalar, dim, 2 * dim - dim + 1> TR1(D);
261
262 // 2.use triangle L to compute the transmissibility of half edge
263 const GlobalPosition& globalPosFace14 = interactionVolume.getFacePosition(idx1, 1);
264
265 // compute normal vectors nu1-nu7 in triangle L for first half edge
266 DimVector nu1L1(0);
267 R_.mv(globalPosFace12 - globalPos1, nu1L1);
268
269 DimVector nu2L1(0);
270 R_.mv(globalPos1 - globalPosFace14, nu2L1);
271
272 DimVector nu3L1(0);
273 R_.mv(globalPosFace14 - globalPos4, nu3L1);
274
275 DimVector nu4L1(0);
276 R_.mv(globalPos4 - globalPosCenter, nu4L1);
277
278 DimVector nu5L1(0);
279 R_.mv(globalPosCenter - globalPos2, nu5L1);
280
281 DimVector nu6L1(0);
282 R_.mv(globalPos2 - globalPosFace12, nu6L1);
283
284 DimVector nu7L1(0);
285 R_.mv(globalPosCenter - globalPos1, nu7L1);
286
287 // compute T, i.e., the area of quadrilateral made by normal vectors 'nu'
288 DimVector Rnu2L1(0);
289 R_.mv(nu2L1, Rnu2L1);
290 Scalar T1L1 = nu1L1 * Rnu2L1;
291
292 DimVector Rnu4L1(0);
293 R_.mv(nu4L1, Rnu4L1);
294 Scalar T2L1 = nu3L1 * Rnu4L1;
295
296 DimVector Rnu6L1(0);
297 R_.mv(nu6L1, Rnu6L1);
298 Scalar T3L1 = nu5L1 * Rnu6L1;
299
300 // compute components needed for flux calculation, denoted as 'omega' and 'chi'
301 DimVector K1nu1L1(0);
302 K1.mv(nu1L1, K1nu1L1);
303 DimVector K1nu2L1(0);
304 K1.mv(nu2L1, K1nu2L1);
305 DimVector K3nu3L1(0);
306 K4.mv(nu3L1, K3nu3L1);
307 DimVector K3nu4L1(0);
308 K4.mv(nu4L1, K3nu4L1);
309 DimVector K2nu5L1(0);
310 K2.mv(nu5L1, K2nu5L1);
311 DimVector K2nu6L1(0);
312 K2.mv(nu6L1, K2nu6L1);
313
314 DimVector Rnu1L1(0);
315 R_.mv(nu1L1, Rnu1L1);
316
317 DimVector &outerNormaln1L1 = interactionVolume.getNormal(idx1, 1);
318
319 Scalar omega111L1 = lambda[idx1][1] * (outerNormaln1L1 * K1nu1L1) * interactionVolume.getFaceArea(idx1, 1) / T1L1;
320 Scalar omega112L1 = lambda[idx1][1] * (outerNormaln1L1 * K1nu2L1) * interactionVolume.getFaceArea(idx1, 1) / T1L1;
321 Scalar omega211L1 = lambda[idx1][0] * (outerNormaln2 * K1nu1L1) * interactionVolume.getFaceArea(idx1, 0) / T1L1;
322 Scalar omega212L1 = lambda[idx1][0] * (outerNormaln2 * K1nu2L1) * interactionVolume.getFaceArea(idx1, 0) / T1L1;
323 Scalar omega123L1 = lambda[idx4][0] * (outerNormaln1L1 * K3nu3L1) * interactionVolume.getFaceArea(idx4, 0) / T2L1;
324 Scalar omega124L1 = lambda[idx4][0] * (outerNormaln1L1 * K3nu4L1) * interactionVolume.getFaceArea(idx4, 0) / T2L1;
325 Scalar omega235L1 = lambda[idx2][1] * (outerNormaln2 * K2nu5L1) * interactionVolume.getFaceArea(idx2, 1) / T3L1;
326 Scalar omega236L1 = lambda[idx2][1] * (outerNormaln2 * K2nu6L1) * interactionVolume.getFaceArea(idx2, 1) / T3L1;
327 Scalar chi711L1 = (nu7L1 * Rnu1L1) / T1L1;
328 Scalar chi712L1 = (nu7L1 * Rnu2L1) / T1L1;
329
330 // compute transmissibility matrix TL1 = CA^{-1}B+D
331 C = 0;
332 A = 0;
333 D = 0;
334 B = 0;
335
336 // evaluate matrix C, D, A, B
337 C[0][0] = -omega111L1;
338 C[0][1] = -omega112L1;
339 C[1][0] = -omega211L1;
340 C[1][1] = -omega212L1;
341
342 D[0][0] = omega111L1 + omega112L1;
343 D[1][0] = omega211L1 + omega212L1;
344
345 A[0][0] = omega111L1 - omega124L1 - omega123L1 * chi711L1;
346 A[0][1] = omega112L1 - omega123L1 * chi712L1;
347 A[1][0] = omega211L1 - omega236L1 * chi711L1;
348 A[1][1] = omega212L1 - omega235L1 - omega236L1 * chi712L1;
349
350 B[0][0] = omega111L1 + omega112L1 + omega123L1 * (1.0 - chi711L1 - chi712L1);
351 B[0][1] = -omega123L1 - omega124L1;
352 B[1][0] = omega211L1 + omega212L1 + omega236L1 * (1.0 - chi711L1 - chi712L1);
353 B[1][2] = -omega235L1 - omega236L1;
354
355 // compute TL1
356 A.invert();
357 D += B.leftmultiply(C.rightmultiply(A));
358 Dune::FieldMatrix<Scalar, dim, 2 * dim - dim + 1> TL1(D);
359
360 //selection criterion
361 using std::abs;
362 Scalar sR = abs(TR1[1][2] - TR1[1][0]);
363 Scalar sL = abs(TL1[1][0] - TL1[1][2]);
364
365 // 3.decide which triangle (which transmissibility coefficients) to use
366 if (sR <= sL)
367 {
368 transmissibility = TR1;
369 return rightTriangle;
370 }
371 else
372 {
373 transmissibility = TL1;
374 return leftTriangle;
375 }
376}
377
392template<class TypeTag>
393int FvMpfaL2dTransmissibilityCalculator<TypeTag>::calculateLeftHNTransmissibility(
394 TransmissibilityType& transmissibilityLeft,
395 InteractionVolume& interactionVolume,
396 std::vector<DimVector >& lambda,
397 int idx1, int idx2, int idx3)
398{
399 auto element1 = interactionVolume.getSubVolumeElement(idx1);
400 auto element2 = interactionVolume.getSubVolumeElement(idx2);
401 auto element3 = interactionVolume.getSubVolumeElement(idx3);
402
403 if (element1.level() != element3.level())
404 {
405 return noTransmissibility;
406 }
407
408 // get global coordinate of cell centers
409 const GlobalPosition& globalPos1 = element1.geometry().center();
410 const GlobalPosition& globalPos2 = element2.geometry().center();
411 const GlobalPosition& globalPos3 = element3.geometry().center();
412
413 const GlobalPosition& globalPosCenter = interactionVolume.getCenterPosition();
414
415 const DimMatrix& K1 = problem_.spatialParams().intrinsicPermeability(element1);
416 const DimMatrix& K2 = problem_.spatialParams().intrinsicPermeability(element2);
417 const DimMatrix& K3 = problem_.spatialParams().intrinsicPermeability(element3);
418
419 const GlobalPosition& globalPosFace12 = interactionVolume.getFacePosition(idx1, 0);
420 DimVector &outerNormaln2 = interactionVolume.getNormal(idx1, 0);
421
422 // compute transmissibility matrix TR1 = CA^{-1}B+D
423 DimMatrix C(0), A(0);
424 TransmissibilityType D(0), B(0);
425
426 // 2.use triangle L to compute the transmissibility of half edge
427 const GlobalPosition& globalPosFace14 = interactionVolume.getFacePosition(idx1, 1);
428
429 // compute normal vectors nu1-nu7 in triangle L for first half edge
430 DimVector nu1L1(0);
431 R_.mv(globalPosFace12 - globalPos1, nu1L1);
432
433 DimVector nu2L1(0);
434 R_.mv(globalPos1 - globalPosFace14, nu2L1);
435
436 DimVector nu3L1(0);
437 R_.mv(globalPosFace14 - globalPos3, nu3L1);
438
439 DimVector nu4L1(0);
440 R_.mv(globalPos3 - globalPosCenter, nu4L1);
441
442 DimVector nu5L1(0);
443 R_.mv(globalPosCenter - globalPos2, nu5L1);
444
445 DimVector nu6L1(0);
446 R_.mv(globalPos2 - globalPosFace12, nu6L1);
447
448 DimVector nu7L1(0);
449 R_.mv(globalPosCenter - globalPos1, nu7L1);
450
451 // compute T, i.e., the area of quadrilateral made by normal vectors 'nu'
452 DimVector Rnu2L1(0);
453 R_.mv(nu2L1, Rnu2L1);
454 Scalar T1L1 = nu1L1 * Rnu2L1;
455
456 DimVector Rnu4L1(0);
457 R_.mv(nu4L1, Rnu4L1);
458 Scalar T2L1 = nu3L1 * Rnu4L1;
459
460 DimVector Rnu6L1(0);
461 R_.mv(nu6L1, Rnu6L1);
462 Scalar T3L1 = nu5L1 * Rnu6L1;
463
464 // compute components needed for flux calculation, denoted as 'omega' and 'chi'
465 DimVector K1nu1L1(0);
466 K1.mv(nu1L1, K1nu1L1);
467 DimVector K1nu2L1(0);
468 K1.mv(nu2L1, K1nu2L1);
469 DimVector K3nu3L1(0);
470 K3.mv(nu3L1, K3nu3L1);
471 DimVector K3nu4L1(0);
472 K3.mv(nu4L1, K3nu4L1);
473 DimVector K2nu5L1(0);
474 K2.mv(nu5L1, K2nu5L1);
475 DimVector K2nu6L1(0);
476 K2.mv(nu6L1, K2nu6L1);
477
478 DimVector Rnu1L1(0);
479 R_.mv(nu1L1, Rnu1L1);
480
481 DimVector &outerNormaln1L1 = interactionVolume.getNormal(idx1, 1);
482
483 Scalar omega111L1 = lambda[idx1][1] * (outerNormaln1L1 * K1nu1L1) * interactionVolume.getFaceArea(idx1, 1) / T1L1;
484 Scalar omega112L1 = lambda[idx1][1] * (outerNormaln1L1 * K1nu2L1) * interactionVolume.getFaceArea(idx1, 1) / T1L1;
485 Scalar omega211L1 = lambda[idx1][0] * (outerNormaln2 * K1nu1L1) * interactionVolume.getFaceArea(idx1, 0) / T1L1;
486 Scalar omega212L1 = lambda[idx1][0] * (outerNormaln2 * K1nu2L1) * interactionVolume.getFaceArea(idx1, 0) / T1L1;
487 Scalar omega123L1 = lambda[idx3][0] * (outerNormaln1L1 * K3nu3L1) * interactionVolume.getFaceArea(idx3, 0) / T2L1;
488 Scalar omega124L1 = lambda[idx3][0] * (outerNormaln1L1 * K3nu4L1) * interactionVolume.getFaceArea(idx3, 0) / T2L1;
489 Scalar omega235L1 = lambda[idx2][1] * (outerNormaln2 * K2nu5L1) * interactionVolume.getFaceArea(idx2, 1) / T3L1;
490 Scalar omega236L1 = lambda[idx2][1] * (outerNormaln2 * K2nu6L1) * interactionVolume.getFaceArea(idx2, 1) / T3L1;
491 Scalar chi711L1 = (nu7L1 * Rnu1L1) / T1L1;
492 Scalar chi712L1 = (nu7L1 * Rnu2L1) / T1L1;
493
494 // compute transmissibility matrix TL1 = CA^{-1}B+D
495 // evaluate matrix C, D, A, B
496 C[0][0] = -omega111L1;
497 C[0][1] = -omega112L1;
498 C[1][0] = -omega211L1;
499 C[1][1] = -omega212L1;
500
501 D[0][0] = omega111L1 + omega112L1;
502 D[1][0] = omega211L1 + omega212L1;
503
504 A[0][0] = omega111L1 - omega124L1 - omega123L1 * chi711L1;
505 A[0][1] = omega112L1 - omega123L1 * chi712L1;
506 A[1][0] = omega211L1 - omega236L1 * chi711L1;
507 A[1][1] = omega212L1 - omega235L1 - omega236L1 * chi712L1;
508
509 B[0][0] = omega111L1 + omega112L1 + omega123L1 * (1.0 - chi711L1 - chi712L1);
510 B[0][1] = -omega123L1 - omega124L1;
511 B[1][0] = omega211L1 + omega212L1 + omega236L1 * (1.0 - chi711L1 - chi712L1);
512 B[1][2] = -omega235L1 - omega236L1;
513
514 // compute TL1
515 A.invert();
516 D += B.leftmultiply(C.rightmultiply(A));
517
518 // 3.decide which triangle (which transmissibility coefficients) to use
519 transmissibilityLeft = D;
520 return leftTriangle;
521}
522
537template<class TypeTag>
538int FvMpfaL2dTransmissibilityCalculator<TypeTag>::calculateRightHNTransmissibility(
539 TransmissibilityType& transmissibilityRight,
540 InteractionVolume& interactionVolume,
541 std::vector<DimVector >& lambda,
542 int idx1, int idx2, int idx3)
543{
544 auto element1 = interactionVolume.getSubVolumeElement(idx1);
545 auto element2 = interactionVolume.getSubVolumeElement(idx2);
546 auto element3 = interactionVolume.getSubVolumeElement(idx3);
547
548 if (element2.level() != element3.level())
549 {
550 return noTransmissibility;
551 }
552
553 // get global coordinate of cell centers
554 const GlobalPosition& globalPos1 = element1.geometry().center();
555 const GlobalPosition& globalPos2 = element2.geometry().center();
556 const GlobalPosition& globalPos3 = element3.geometry().center();
557
558 const GlobalPosition& globalPosCenter = interactionVolume.getCenterPosition();
559
560 const DimMatrix& K1 = problem_.spatialParams().intrinsicPermeability(element1);
561 const DimMatrix& K2 = problem_.spatialParams().intrinsicPermeability(element2);
562 const DimMatrix& K3 = problem_.spatialParams().intrinsicPermeability(element3);
563
564 const GlobalPosition& globalPosFace12 = interactionVolume.getFacePosition(idx1, 0);
565 const GlobalPosition& globalPosFace23 = interactionVolume.getFacePosition(idx2, 0);
566
567 // compute normal vectors nu1-nu7 in triangle R for first half edge
568 DimVector nu1R1(0);
569 R_.mv(globalPosFace12 - globalPos2, nu1R1);
570
571 DimVector nu2R1(0);
572 R_.mv(globalPos2 - globalPosFace23, nu2R1);
573
574 DimVector nu3R1(0);
575 R_.mv(globalPosFace23 - globalPos3, nu3R1);
576
577 DimVector nu4R1(0);
578 R_.mv(globalPos3 - globalPosCenter, nu4R1);
579
580 DimVector nu5R1(0);
581 R_.mv(globalPosCenter - globalPos1, nu5R1);
582
583 DimVector nu6R1(0);
584 R_.mv(globalPos1 - globalPosFace12, nu6R1);
585
586 DimVector nu7R1(0);
587 R_.mv(globalPosCenter - globalPos2, nu7R1);
588
589 // compute T, i.e., the area of quadrilateral made by normal vectors 'nu'
590 DimVector Rnu2R1(0);
591 R_.mv(nu2R1, Rnu2R1);
592 Scalar T1R1 = nu1R1 * Rnu2R1;
593
594 DimVector Rnu4R1(0);
595 R_.mv(nu4R1, Rnu4R1);
596 Scalar T2R1 = nu3R1 * Rnu4R1;
597
598 DimVector Rnu6R1(0);
599 R_.mv(nu6R1, Rnu6R1);
600 Scalar T3R1 = nu5R1 * Rnu6R1;
601
602 // compute components needed for flux calculation, denoted as 'omega' and 'chi'
603 DimVector K2nu1R1(0);
604 K2.mv(nu1R1, K2nu1R1);
605 DimVector K2nu2R1(0);
606 K2.mv(nu2R1, K2nu2R1);
607 DimVector K3nu3R1(0);
608 K3.mv(nu3R1, K3nu3R1);
609 DimVector K3nu4R1(0);
610 K3.mv(nu4R1, K3nu4R1);
611 DimVector K1nu5R1(0);
612 K1.mv(nu5R1, K1nu5R1);
613 DimVector K1nu6R1(0);
614 K1.mv(nu6R1, K1nu6R1);
615
616 DimVector Rnu1R1(0);
617 R_.mv(nu1R1, Rnu1R1);
618
619 DimVector &outerNormaln1R1 = interactionVolume.getNormal(idx2, 0);
620 DimVector &outerNormaln2 = interactionVolume.getNormal(idx1, 0);
621
622 Scalar omega111R1 = lambda[idx2][0] * (outerNormaln1R1 * K2nu1R1) * interactionVolume.getFaceArea(idx2, 0) / T1R1;
623 Scalar omega112R1 = lambda[idx2][0] * (outerNormaln1R1 * K2nu2R1) * interactionVolume.getFaceArea(idx2, 0) / T1R1;
624 Scalar omega211R1 = lambda[idx2][1] * (outerNormaln2 * K2nu1R1) * interactionVolume.getFaceArea(idx2, 1) / T1R1;
625 Scalar omega212R1 = lambda[idx2][1] * (outerNormaln2 * K2nu2R1) * interactionVolume.getFaceArea(idx2, 1) / T1R1;
626 Scalar omega123R1 = lambda[idx3][1] * (outerNormaln1R1 * K3nu3R1) * interactionVolume.getFaceArea(idx3, 1) / T2R1;
627 Scalar omega124R1 = lambda[idx3][1] * (outerNormaln1R1 * K3nu4R1) * interactionVolume.getFaceArea(idx3, 1) / T2R1;
628 Scalar omega235R1 = lambda[idx1][0] * (outerNormaln2 * K1nu5R1) * interactionVolume.getFaceArea(idx1, 0) / T3R1;
629 Scalar omega236R1 = lambda[idx1][0] * (outerNormaln2 * K1nu6R1) * interactionVolume.getFaceArea(idx1, 0) / T3R1;
630 Scalar chi711R1 = (nu7R1 * Rnu1R1) / T1R1;
631 Scalar chi712R1 = (nu7R1 * Rnu2R1) / T1R1;
632
633 // compute transmissibility matrix TR1 = CA^{-1}B+D
634 DimMatrix C(0), A(0);
635 TransmissibilityType D(0), B(0);
636
637 // evaluate matrix C, D, A, B
638 C[0][0] = -omega111R1;
639 C[0][1] = -omega112R1;
640 C[1][0] = -omega211R1;
641 C[1][1] = -omega212R1;
642
643 D[0][0] = omega111R1 + omega112R1;
644 D[1][0] = omega211R1 + omega212R1;
645
646 A[0][0] = omega111R1 - omega124R1 - omega123R1 * chi711R1;
647 A[0][1] = omega112R1 - omega123R1 * chi712R1;
648 A[1][0] = omega211R1 - omega236R1 * chi711R1;
649 A[1][1] = omega212R1 - omega235R1 - omega236R1 * chi712R1;
650
651 B[0][0] = omega111R1 + omega112R1 + omega123R1 * (1.0 - chi711R1 - chi712R1);
652 B[0][1] = -omega123R1 - omega124R1;
653 B[1][0] = omega211R1 + omega212R1 + omega236R1 * (1.0 - chi711R1 - chi712R1);
654 B[1][2] = -omega235R1 - omega236R1;
655
656 // compute TR1
657 A.invert();
658 D += B.leftmultiply(C.rightmultiply(A));
659
660 transmissibilityRight = D;
661 return rightTriangle;
662}
663} // end namespace Dumux
664#endif
linteractionvolume.hh
Class including the information of an interaction volume of a MPFA L-method that does not change with...
Dumux::FvMpfaL2dTransmissibilityCalculator::calculateTransmissibility
int calculateTransmissibility(TransmissibilityType &transmissibility, InteractionVolume &interactionVolume, std::vector< DimVector > &lambda, int idx1, int idx2, int idx3, int idx4)
Calculates tranmissibility matrix.
Definition: 2dtransmissibilitycalculator.hh:136
Dumux::FvMpfaL2dTransmissibilityCalculator::calculateRightHNTransmissibility
int calculateRightHNTransmissibility(TransmissibilityType &transmissibility, InteractionVolume &interactionVolume, std::vector< DimVector > &lambda, int idx1, int idx2, int idx3)
Calculates tranmissibility matrix of right L-shape.
Definition: 2dtransmissibilitycalculator.hh:538
Dumux::FvMpfaL2dTransmissibilityCalculator::calculateLeftHNTransmissibility
int calculateLeftHNTransmissibility(TransmissibilityType &transmissibility, InteractionVolume &interactionVolume, std::vector< DimVector > &lambda, int idx1, int idx2, int idx3)
Calculates tranmissibility matrix of left L-shape.
Definition: 2dtransmissibilitycalculator.hh:393
Dumux
Definition: adapt.hh:29
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::FvMpfaL2dTransmissibilityCalculator
Provides methods for transmissibility calculation in 2-d.
Definition: 2dtransmissibilitycalculator.hh:44
Dumux::FvMpfaL2dTransmissibilityCalculator::FvMpfaL2dTransmissibilityCalculator
FvMpfaL2dTransmissibilityCalculator(Problem &problem)
Constructs a FvMpfaL2dTransmissibilityCalculator object.
Definition: 2dtransmissibilitycalculator.hh:99
Dumux::FvMpfaL2dTransmissibilityCalculator::TransmissibilityType
Dune::FieldMatrix< Scalar, dim, 2 *dim - dim+1 > TransmissibilityType
Type of the transmissibility matrix.
Definition: 2dtransmissibilitycalculator.hh:65
Dumux::FvMpfaL2dTransmissibilityCalculator::leftTriangle
@ leftTriangle
Left L-shape.
Definition: 2dtransmissibilitycalculator.hh:70
Dumux::FvMpfaL2dTransmissibilityCalculator::rightTriangle
@ rightTriangle
Right L-shape.
Definition: 2dtransmissibilitycalculator.hh:72
Dumux::FvMpfaL2dTransmissibilityCalculator::noTransmissibility
@ noTransmissibility
No transmissibility calculated.
Definition: 2dtransmissibilitycalculator.hh:71
Dumux::FVMPFALInteractionVolume
Class including the information of an interaction volume of a MPFA L-method that does not change with...
Definition: linteractionvolume.hh:41
Dumux::FVMPFALInteractionVolume::getNormal
DimVector & getNormal(int subVolumeIdx, int subVolumeFaceIdxInInside)
Get a flux face normal.
Definition: linteractionvolume.hh:364
Dumux::FVMPFALInteractionVolume::getSubVolumeElement
Element getSubVolumeElement(int subVolumeIdx)
Get an element of the interaction volume.
Definition: linteractionvolume.hh:281
Dumux::FVMPFALInteractionVolume::getFacePosition
DimVector & getFacePosition(int subVolumeIdx, int subVolumeFaceIdxInInside)
Get position of the center of a flux face.
Definition: linteractionvolume.hh:376
Dumux::FVMPFALInteractionVolume::getFaceArea
Scalar & getFaceArea(int subVolumeIdx, int subVolumeFaceIdxInInside)
Get a flux face area.
Definition: linteractionvolume.hh:388
Dumux::FVMPFALInteractionVolume::getCenterPosition
DimVector & getCenterPosition()
Get position vector of central vertex.
Definition: linteractionvolume.hh:240
properties.hh
Properties for a MPFA method.
pressureproperties.hh
Base file for properties related to sequential IMPET algorithms.
Generated by doxygen 1.9.3