3.5-git
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
multidomain/facet/cellcentered/mpfa/localassembler.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 2 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 *****************************************************************************/
24#ifndef DUMUX_MULTIDOMAIN_FACET_CC_MPFA_O_LOCAL_ASSEMBLER_HH
25#define DUMUX_MULTIDOMAIN_FACET_CC_MPFA_O_LOCAL_ASSEMBLER_HH
26
27#include <vector>
28
29#include <dune/common/exceptions.hh>
30#include <dune/common/float_cmp.hh>
31
32#include <dumux/common/math.hh>
39
40namespace Dumux {
41
52template< class P, class EG, class EV >
54: public InteractionVolumeAssemblerBase< P, EG, EV >
55{
59
60 template< class IV >
61 using Scalar = typename IV::Traits::MatVecTraits::FaceVector::value_type;
62
63public:
65 using ParentType::ParentType;
66
80 template< class DataHandle, class IV, class TensorFunc >
81 void assembleMatrices(DataHandle& handle, IV& iv, const TensorFunc& getT, Scalar<IV> wijZeroThresh = 0.0)
82 {
83 assembleLocalMatrices_(handle.A(), handle.AB(), handle.CA(), handle.T(), handle.omegas(), iv, getT, wijZeroThresh);
84
85 // maybe solve the local system
86 if (iv.numUnknowns() > 0)
87 Helper::solveLocalSystem(this->fvGeometry(), handle, iv);
88 }
89
98 template< class DataHandle, class IV, class GetU >
99 void assembleU(DataHandle& handle, const IV& iv, const GetU& getU)
100 {
101 auto& u = handle.uj();
102 Helper::resizeVector(u, iv.numKnowns());
103
104 // put the cell unknowns first ...
105 typename IV::Traits::IndexSet::LocalIndexType i = 0;
106 for (; i < iv.numScvs(); i++)
107 u[i] = getU( this->elemVolVars()[iv.localScv(i).gridScvIndex()] );
108
109 // ... then put facet unknowns ...
110 for (const auto& data : iv.interiorBoundaryData())
111 {
112 const auto& scvf = this->fvGeometry().scvf(data.scvfIndex());
113 const auto element = this->fvGeometry().gridGeometry().element(scvf.insideScvIdx());
114 u[i++] = getU( this->problem().couplingManager().getLowDimVolVars(element, scvf) );
115 }
116
117 // ... then put the Dirichlet unknowns
118 for (const auto& data : iv.dirichletData())
119 u[i++] = getU( this->elemVolVars()[data.volVarIndex()] );
120 }
121
130 template< class DataHandle, class IV, class GetRho >
131 void assembleGravity(DataHandle& handle, const IV& iv, const GetRho& getRho)
132 {
133 using GridView = typename IV::Traits::GridView;
134 static constexpr int dim = GridView::dimension;
135 static constexpr int dimWorld = GridView::dimensionworld;
136 static constexpr bool isSurfaceGrid = dim < dimWorld;
137
138 // resize the gravity vectors
139 auto& g = handle.g();
140 auto& deltaG = handle.deltaG();
141 auto& outsideG = handle.gOutside();
142 Helper::resizeVector(g, iv.numFaces());
143 Helper::resizeVector(deltaG, iv.numUnknowns());
144 if (isSurfaceGrid)
145 Helper::resizeVector(outsideG, iv.numFaces());
146
151 using Scalar = typename IV::Traits::MatVecTraits::TMatrix::value_type;
152 using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType;
153
154 // xi factor for coupling conditions
155 static const Scalar xi = getParamFromGroup<Scalar>(this->problem().paramGroup(), "FacetCoupling.Xi", 1.0);
156
157 for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx)
158 {
159 // gravitational acceleration on this face
160 const auto& curLocalScvf = iv.localScvf(faceIdx);
161 const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.gridScvfIndex());
162 const auto& gravity = this->problem().spatialParams().gravity(curGlobalScvf.ipGlobal());
163 const auto curIsInteriorBoundary = curLocalScvf.isOnInteriorBoundary();
164 const Scalar curXiFactor = curIsInteriorBoundary ? (curGlobalScvf.boundary() ? 1.0 : xi) : 1.0;
165
166 // get permeability tensor in "positive" sub volume
167 const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices();
168 const auto& posGlobalScv = this->fvGeometry().scv(iv.localScv(neighborScvIndices[0]).gridScvIndex());
169 const auto& posVolVars = this->elemVolVars()[posGlobalScv];
170 const auto alpha_inside = posVolVars.extrusionFactor()*vtmv(curGlobalScvf.unitOuterNormal(),
171 posVolVars.permeability(),
172 gravity);
173
174 const auto numOutsideFaces = !curGlobalScvf.boundary() ? curGlobalScvf.numOutsideScvs() : 0;
175 using OutsideAlphaStorage = std::conditional_t< isSurfaceGrid,
176 std::vector<Scalar>,
177 Dune::ReservedVector<Scalar, 1> >;
178 OutsideAlphaStorage alpha_outside; alpha_outside.resize(numOutsideFaces);
179 std::fill(alpha_outside.begin(), alpha_outside.end(), 0.0);
180 Scalar rho;
181
182 if (isSurfaceGrid && !curIsInteriorBoundary)
183 Helper::resizeVector(outsideG[faceIdx], numOutsideFaces);
184
185 if (!curLocalScvf.isDirichlet())
186 {
187 const auto localDofIdx = curLocalScvf.localDofIndex();
188 const Scalar curOneMinusXi = curIsInteriorBoundary ? -(1.0 - xi) : 1.0;
189
190 rho = getRho(posVolVars);
191 deltaG[localDofIdx] = 0.0;
192
193 if (!curGlobalScvf.boundary())
194 {
195 for (unsigned int idxInOutside = 0; idxInOutside < curGlobalScvf.numOutsideScvs(); ++idxInOutside)
196 {
197 // obtain outside tensor
198 const auto negLocalScvIdx = neighborScvIndices[idxInOutside+1];
199 const auto& negGlobalScv = this->fvGeometry().scv(iv.localScv(negLocalScvIdx).gridScvIndex());
200 const auto& negVolVars = this->elemVolVars()[negGlobalScv];
201 const auto& flipScvf = !isSurfaceGrid ? curGlobalScvf
202 : this->fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside);
203
204 alpha_outside[idxInOutside] = negVolVars.extrusionFactor()*vtmv(flipScvf.unitOuterNormal(),
205 negVolVars.permeability(),
206 gravity);
207 if (isSurfaceGrid)
208 alpha_outside[idxInOutside] *= -1.0;
209
210 if (!curIsInteriorBoundary)
211 rho += getRho(negVolVars);
212
213 deltaG[localDofIdx] += curOneMinusXi*alpha_outside[idxInOutside];
214 }
215 }
216
217 if (curIsInteriorBoundary)
218 {
219 const auto posElement = this->fvGeometry().gridGeometry().element(posGlobalScv.elementIndex());
220 const auto& facetVolVars = this->problem().couplingManager().getLowDimVolVars(posElement, curGlobalScvf);
221 rho += getRho(facetVolVars);
222 rho /= 2.0;
223 const auto alphaFacet = posVolVars.extrusionFactor()*vtmv(curGlobalScvf.unitOuterNormal(),
224 facetVolVars.permeability(),
225 gravity);
226 deltaG[localDofIdx] += alphaFacet;
227 }
228 else
229 rho /= numOutsideFaces + 1;
230
231 deltaG[localDofIdx] -= curXiFactor*alpha_inside;
232 deltaG[localDofIdx] *= rho*Extrusion::area(curGlobalScvf);
233 }
234 // use average density between facet and cell density on interior Dirichlet boundaries
235 else if (curIsInteriorBoundary)
236 {
237 const auto posElement = this->fvGeometry().gridGeometry().element(posGlobalScv.elementIndex());
238 const auto& facetVolVars = this->problem().couplingManager().getLowDimVolVars(posElement, curGlobalScvf);
239 rho = 0.5*(getRho(facetVolVars) + getRho(posVolVars));
240 }
241 // use density resulting from Dirichlet BCs
242 else
243 rho = getRho(this->elemVolVars()[curGlobalScvf.outsideScvIdx()]);
244
245 // add "inside" & "outside" alphas to gravity containers
246 g[faceIdx] = alpha_inside*rho*Extrusion::area(curGlobalScvf);
247
248 if (isSurfaceGrid && !curIsInteriorBoundary)
249 {
250 unsigned int i = 0;
251 for (const auto& alpha : alpha_outside)
252 outsideG[faceIdx][i++] = alpha*rho*Extrusion::area(curGlobalScvf);
253 }
254 }
255
256 // add iv-wide contributions to gravity vectors
257 handle.CA().umv(deltaG, g);
258 if (isSurfaceGrid)
259 {
260 using FaceVector = typename IV::Traits::MatVecTraits::FaceVector;
261 FaceVector AG;
262 Helper::resizeVector(AG, iv.numUnknowns());
263 handle.A().mv(deltaG, AG);
264
265 // compute gravitational accelerations
266 for (const auto& localFaceData : iv.localFaceData())
267 {
268 // continue only for "outside" faces
269 if (!localFaceData.isOutsideFace())
270 continue;
271
272 const auto localScvIdx = localFaceData.ivLocalInsideScvIndex();
273 const auto localScvfIdx = localFaceData.ivLocalScvfIndex();
274 const auto idxInOutside = localFaceData.scvfLocalOutsideScvfIndex();
275 const auto& posLocalScv = iv.localScv(localScvIdx);
276 const auto& wijk = handle.omegas()[localScvfIdx][idxInOutside+1];
277
278 // add contributions from all local directions
279 for (LocalIndexType localDir = 0; localDir < dim; localDir++)
280 {
281 // the scvf corresponding to this local direction in the scv
282 const auto& curLocalScvf = iv.localScvf(posLocalScv.localScvfIndex(localDir));
283 if (!curLocalScvf.isDirichlet())
284 outsideG[localScvfIdx][idxInOutside] -= wijk[localDir]*AG[curLocalScvf.localDofIndex()];
285 }
286 }
287 }
288 }
289
290private:
294 template< class IV, class TensorFunc >
295 void assembleLocalMatrices_(typename IV::Traits::MatVecTraits::AMatrix& A,
296 typename IV::Traits::MatVecTraits::BMatrix& B,
297 typename IV::Traits::MatVecTraits::CMatrix& C,
298 typename IV::Traits::MatVecTraits::DMatrix& D,
299 typename IV::Traits::MatVecTraits::OmegaStorage& wijk,
300 IV& iv, const TensorFunc& getT,
301 Scalar<IV> wijZeroThresh)
302 {
303 using Scalar = typename IV::Traits::MatVecTraits::TMatrix::value_type;
304 using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType;
305 static constexpr int dim = IV::Traits::GridView::dimension;
306 static constexpr int dimWorld = IV::Traits::GridView::dimensionworld;
307
308 // xi factor for coupling conditions
309 static const Scalar xi = getParamFromGroup<Scalar>(this->problem().paramGroup(), "FacetCoupling.Xi", 1.0);
310
311 // On surface grids only xi = 1.0 can be used, as the coupling condition
312 // for xi != 1.0 does not generalize for surface grids where the normal
313 // vectors of the inside/outside elements have different orientations.
314 if (dim < dimWorld)
315 if (Dune::FloatCmp::ne(xi, 1.0, 1e-6))
316 DUNE_THROW(Dune::InvalidStateException, "Xi != 1.0 cannot be used on surface grids");
317
318 std::vector< std::pair<LocalIndexType, LocalIndexType> > faceMarkers;
319 Helper::resizeVector(wijk, iv.numFaces());
320
321 // if only Dirichlet faces are present in the iv,
322 // the matrices A, B & C are undefined and D = T
323 if (iv.numUnknowns() == 0)
324 {
325 // resize & reset D matrix
326 Helper::resizeMatrix(D, iv.numFaces(), iv.numKnowns()); D = 0.0;
327
328 // Loop over all the faces, in this case these are all dirichlet boundaries
329 for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx)
330 {
331 const auto& curLocalScvf = iv.localScvf(faceIdx);
332 const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.gridScvfIndex());
333 const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices();
334
335 // get tensor in "positive" sub volume
336 const auto& posLocalScv = iv.localScv(neighborScvIndices[0]);
337 const auto& posGlobalScv = this->fvGeometry().scv(posLocalScv.gridScvIndex());
338 const auto& posVolVars = this->elemVolVars()[posGlobalScv];
339 const auto tensor = getT(posVolVars);
340
341 // the omega factors of the "positive" sub volume
342 Helper::resizeVector(wijk[faceIdx], /*no outside scvs present*/1);
343 wijk[faceIdx][0] = computeMpfaTransmissibility<EG>(posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor());
344
345 const auto posScvLocalDofIdx = posLocalScv.localDofIndex();
346 for (LocalIndexType localDir = 0; localDir < dim; localDir++)
347 {
348 const auto& otherLocalScvf = iv.localScvf( posLocalScv.localScvfIndex(localDir) );
349 const auto otherLocalDofIdx = otherLocalScvf.localDofIndex();
350 D[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir];
351 D[faceIdx][posScvLocalDofIdx] += wijk[faceIdx][0][localDir];
352 }
353 }
354 }
355 else
356 {
357 // resize & reset matrices
358 Helper::resizeMatrix(A, iv.numUnknowns(), iv.numUnknowns()); A = 0.0;
359 Helper::resizeMatrix(B, iv.numUnknowns(), iv.numKnowns()); B = 0.0;
360 Helper::resizeMatrix(C, iv.numFaces(), iv.numUnknowns()); C = 0.0;
361 Helper::resizeMatrix(D, iv.numFaces(), iv.numKnowns()); D = 0.0;
362
363 for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx)
364 {
365 const auto& curLocalScvf = iv.localScvf(faceIdx);
366 const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.gridScvfIndex());
367 const auto curIsDirichlet = curLocalScvf.isDirichlet();
368 const auto curLocalDofIdx = curLocalScvf.localDofIndex();
369 const auto curIsInteriorBoundary = curLocalScvf.isOnInteriorBoundary();
370 const Scalar curXiFactor = curIsInteriorBoundary ? (curGlobalScvf.boundary() ? 1.0 : xi) : 1.0;
371
372 // get tensor in "positive" sub volume
373 const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices();
374 const auto& posLocalScv = iv.localScv(neighborScvIndices[0]);
375 const auto& posGlobalScv = this->fvGeometry().scv(posLocalScv.gridScvIndex());
376 const auto& posVolVars = this->elemVolVars()[posGlobalScv];
377 const auto& posElement = iv.element(neighborScvIndices[0]);
378 const auto tensor = getT(posVolVars);
379
380 // the omega factors of the "positive" sub volume
381 Helper::resizeVector(wijk[faceIdx], neighborScvIndices.size());
382 wijk[faceIdx][0] = computeMpfaTransmissibility<EG>(posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor());
383
384 using std::abs;
385 bool isZeroWij = false;
386
387 // go over the coordinate directions in the positive sub volume
388 for (unsigned int localDir = 0; localDir < dim; localDir++)
389 {
390 const auto& otherLocalScvf = iv.localScvf( posLocalScv.localScvfIndex(localDir) );
391 const auto otherLocalDofIdx = otherLocalScvf.localDofIndex();
392
393 if (otherLocalDofIdx == curLocalDofIdx)
394 {
395 if (abs(wijk[faceIdx][0][localDir]) <= wijZeroThresh)
396 {
397 if (!curIsDirichlet)
398 {
399 isZeroWij = true;
400 faceMarkers.emplace_back( std::make_pair(curLocalDofIdx, faceIdx) );
401 }
402 }
403 }
404
405 // if we are not on a Dirichlet face, add entries associated with unknown face pressures
406 // i.e. in matrix C and maybe A (if current face is not a Dirichlet face)
407 if (!otherLocalScvf.isDirichlet())
408 {
409 C[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir];
410 if (!curIsDirichlet)
411 A[curLocalDofIdx][otherLocalDofIdx] -= curXiFactor*wijk[faceIdx][0][localDir];
412 }
413 // the current face is a Dirichlet face and creates entries in D & maybe B
414 else
415 {
416 D[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir];
417 if (!curIsDirichlet)
418 B[curLocalDofIdx][otherLocalDofIdx] += curXiFactor*wijk[faceIdx][0][localDir];
419 }
420
421 // add entries related to pressures at the scv centers (dofs)
422 const auto posScvLocalDofIdx = posLocalScv.localDofIndex();
423 D[faceIdx][posScvLocalDofIdx] += wijk[faceIdx][0][localDir];
424
425 if (!curIsDirichlet)
426 B[curLocalDofIdx][posScvLocalDofIdx] -= curXiFactor*wijk[faceIdx][0][localDir];
427 }
428
429 // handle the entries related to the coupled facet element
430 if (curIsInteriorBoundary && !curIsDirichlet)
431 {
432 const auto facetVolVars = this->problem().couplingManager().getLowDimVolVars(posElement, curGlobalScvf);
433 const auto facetTensor = getT(facetVolVars);
434
435 // On surface grids we use the square root of the extrusion factor as approximation of the aperture
436 using std::sqrt;
437 const auto wFacet = 2.0*Extrusion::area(curGlobalScvf)*posVolVars.extrusionFactor()
438 *vtmv(curGlobalScvf.unitOuterNormal(), facetTensor, curGlobalScvf.unitOuterNormal())
439 / (dim < dimWorld ? sqrt(facetVolVars.extrusionFactor()) : facetVolVars.extrusionFactor());
440
441 A[curLocalDofIdx][curLocalDofIdx] -= wFacet;
442 B[curLocalDofIdx][curLocalScvf.coupledFacetLocalDofIndex()] -= wFacet;
443
444 // check for zero transmissibilities (skip if inside has been zero already)
445 if (!isZeroWij && abs(wFacet) <= wijZeroThresh)
446 {
447 isZeroWij = true;
448 faceMarkers.emplace_back( std::make_pair(curLocalDofIdx, faceIdx) );
449 }
450 }
451
452 // If we are on an interior face (which isn't coupled via Dirichlet Bs), add values from negative sub volume
453 if (!curGlobalScvf.boundary() && !curIsDirichlet)
454 {
455 // the minus ensures the right sign of the coupling condition in the matrices
456 const Scalar curOneMinusXi = curIsInteriorBoundary ? -(1.0 - xi) : 1.0;
457
458 // loop over all the outside neighbors of this face and add entries
459 for (unsigned int idxInOutside = 0; idxInOutside < curGlobalScvf.numOutsideScvs(); ++idxInOutside)
460 {
461 const auto idxOnScvf = idxInOutside+1;
462 const auto& negLocalScv = iv.localScv( neighborScvIndices[idxOnScvf] );
463 const auto& negGlobalScv = this->fvGeometry().scv(negLocalScv.gridScvIndex());
464 const auto& negVolVars = this->elemVolVars()[negGlobalScv];
465 const auto negTensor = getT(negVolVars);
466
467 // On surface grids, use outside face for "negative" transmissibility calculation
468 const auto& scvf = dim < dimWorld ? this->fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside)
469 : curGlobalScvf;
470 wijk[faceIdx][idxOnScvf] = computeMpfaTransmissibility<EG>(negLocalScv, scvf, negTensor, negVolVars.extrusionFactor());
471
472 // flip sign on surface grids (since we used the "outside" normal)
473 if (dim < dimWorld)
474 wijk[faceIdx][idxOnScvf] *= -1.0;
475
476 // go over the coordinate directions in the negative sub volume
477 for (int localDir = 0; localDir < dim; localDir++)
478 {
479 const auto otherLocalScvfIdx = negLocalScv.localScvfIndex(localDir);
480 const auto& otherLocalScvf = iv.localScvf(otherLocalScvfIdx);
481 const auto otherLocalDofIdx = otherLocalScvf.localDofIndex();
482
483 // check for zero transmissibilities (skip if inside has been zero already)
484 if (otherLocalDofIdx == curLocalDofIdx && !isZeroWij)
485 if (abs(wijk[faceIdx][idxOnScvf][localDir]) <= wijZeroThresh)
486 faceMarkers.emplace_back( std::make_pair(curLocalDofIdx, faceIdx) );
487
488 if (!otherLocalScvf.isDirichlet())
489 A[curLocalDofIdx][otherLocalDofIdx] += curOneMinusXi*wijk[faceIdx][idxOnScvf][localDir];
490 else
491 B[curLocalDofIdx][otherLocalDofIdx] -= curOneMinusXi*wijk[faceIdx][idxOnScvf][localDir];
492
493 // add entries to matrix B
494 B[curLocalDofIdx][negLocalScv.localDofIndex()] += curOneMinusXi*wijk[faceIdx][idxOnScvf][localDir];
495 }
496 }
497 }
498 }
499 }
500 }
501};
502
503} // end namespace
504
505#endif
Define some often used mathematical functions.
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
A class that contains helper functions as well as functionality which is common to different mpfa sch...
The available mpfa schemes in Dumux.
Defines the general interface of classes used for the assembly of the local systems of equations invo...
Helper classes to compute the integration elements.
Dune::DenseMatrix< MAT >::value_type vtmv(const Dune::DenseVector< V1 > &v1, const Dune::DenseMatrix< MAT > &M, const Dune::DenseVector< V2 > &v2)
Evaluates the scalar product of a vector v2, projected by a matrix M, with a vector v1.
Definition: math.hh:863
Definition: adapt.hh:29
typename Extrusion< T >::type Extrusion_t
Convenience alias for obtaining the extrusion type.
Definition: extrusion.hh:177
Defines the general interface of the local assembler classes for the assembly of the interaction volu...
Definition: localassemblerbase.hh:57
const FVElementGeometry & fvGeometry() const
Definition: localassemblerbase.hh:87
const ElementVolumeVariables & elemVolVars() const
Definition: localassemblerbase.hh:88
const Problem & problem() const
Definition: localassemblerbase.hh:86
A class that contains helper functions as well as functionality which is common to different mpfa sch...
Definition: localassemblerhelper.hh:46
static void solveLocalSystem(const FVElementGeometry &fvGeometry, DataHandle &handle, IV &iv)
Solves a previously assembled iv-local system of equations and stores the resulting transmissibilitie...
Definition: localassemblerhelper.hh:82
static void resizeMatrix(Matrix &M, size_type rows, size_type cols)
resizes a matrix to the given sizes (specialization for dynamic matrix type)
Definition: localassemblerhelper.hh:222
static void resizeVector(Vector &v, size_type size)
resizes a vector to the given size (specialization for dynamic matrix type)
Definition: localassemblerhelper.hh:238
Specialization of the interaction volume-local assembler class for the schemes using an mpfa-o type a...
Definition: multidomain/facet/cellcentered/mpfa/localassembler.hh:55
void assembleU(DataHandle &handle, const IV &iv, const GetU &getU)
Assembles the vector of primary (cell) unknowns and (maybe) Dirichlet boundary conditions within an i...
Definition: multidomain/facet/cellcentered/mpfa/localassembler.hh:99
void assembleMatrices(DataHandle &handle, IV &iv, const TensorFunc &getT, Scalar< IV > wijZeroThresh=0.0)
Assembles the matrices involved in the flux expressions and the local system of equations within an i...
Definition: multidomain/facet/cellcentered/mpfa/localassembler.hh:81
void assembleGravity(DataHandle &handle, const IV &iv, const GetRho &getRho)
Assembles the gravitational flux contributions on the scvfs within an interaction volume.
Definition: multidomain/facet/cellcentered/mpfa/localassembler.hh:131
This file contains free functions to evaluate the transmissibilities associated with flux evaluations...