26#ifndef DUMUX_DISCRETIZATION_CC_MPFA_LOCAL_ASSEMBLER_BASE_HH
27#define DUMUX_DISCRETIZATION_CC_MPFA_LOCAL_ASSEMBLER_BASE_HH
33#include <dune/common/exceptions.hh>
34#include <dune/common/reservedvector.hh>
53template<
class P,
class EG,
class EV >
58 using FVElementGeometry = EG;
59 using ElementVolumeVariables = EV;
63 using Scalar =
typename IV::Traits::MatVecTraits::FaceVector::value_type;
84 const Problem&
problem()
const {
return *problemPtr_; }
85 const FVElementGeometry&
fvGeometry()
const {
return *fvGeometryPtr_; }
86 const ElementVolumeVariables&
elemVolVars()
const {
return *elemVolVarsPtr_; }
103 template<
class DataHandle,
class IV,
class TensorFunc >
104 void assembleMatrices(DataHandle& handle, IV& iv,
const TensorFunc& getT, Scalar<IV> wijZeroThresh = 0.0)
106 DUNE_THROW(Dune::NotImplemented,
"Implementation does not provide an assembleMatrices() function");
120 template<
class DataHandle,
class IV,
class GetU >
121 void assembleU(DataHandle& handle,
const IV& iv,
const GetU& getU)
123 DUNE_THROW(Dune::NotImplemented,
"Implementation does not provide an assemble() function for the cell/Dirichlet unknowns");
134 template<
class DataHandle,
class IV,
class GetRho >
137 using GridView =
typename IV::Traits::GridView;
138 static constexpr int dim = GridView::dimension;
139 static constexpr int dimWorld = GridView::dimensionworld;
140 static constexpr bool isSurfaceGrid = dim < dimWorld;
143 auto& g = handle.g();
144 auto& deltaG = handle.deltaG();
145 auto& outsideG = handle.gOutside();
155 using Scalar =
typename IV::Traits::MatVecTraits::TMatrix::value_type;
156 using LocalIndexType =
typename IV::Traits::IndexSet::LocalIndexType;
158 for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx)
161 const auto& curLocalScvf = iv.localScvf(faceIdx);
162 const auto& curGlobalScvf =
fvGeometry().scvf(curLocalScvf.gridScvfIndex());
163 const auto& gravity =
problem().spatialParams().gravity(curGlobalScvf.ipGlobal());
166 const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices();
167 const auto& posGlobalScv =
fvGeometry().scv(iv.localScv(neighborScvIndices[0]).gridScvIndex());
168 const auto& posVolVars =
elemVolVars()[posGlobalScv];
169 const auto alpha_inside = posVolVars.extrusionFactor()*
vtmv(curGlobalScvf.unitOuterNormal(),
170 posVolVars.permeability(),
173 const auto numOutsideFaces = !curGlobalScvf.boundary() ? curGlobalScvf.numOutsideScvs() : 0;
174 using OutsideAlphaStorage = std::conditional_t< isSurfaceGrid,
176 Dune::ReservedVector<Scalar, 1> >;
177 OutsideAlphaStorage alpha_outside; alpha_outside.resize(numOutsideFaces);
178 std::fill(alpha_outside.begin(), alpha_outside.end(), 0.0);
184 if (!curLocalScvf.isDirichlet())
186 const auto localDofIdx = curLocalScvf.localDofIndex();
188 rho = getRho(posVolVars);
189 deltaG[localDofIdx] = 0.0;
191 if (!curGlobalScvf.boundary())
193 for (
unsigned int idxInOutside = 0; idxInOutside < curGlobalScvf.numOutsideScvs(); ++idxInOutside)
196 const auto negLocalScvIdx = neighborScvIndices[idxInOutside+1];
197 const auto& negGlobalScv =
fvGeometry().scv(iv.localScv(negLocalScvIdx).gridScvIndex());
198 const auto& negVolVars =
elemVolVars()[negGlobalScv];
199 const auto& flipScvf = !isSurfaceGrid ? curGlobalScvf
200 :
fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside);
202 alpha_outside[idxInOutside] = negVolVars.extrusionFactor()*
vtmv(flipScvf.unitOuterNormal(),
203 negVolVars.permeability(),
206 alpha_outside[idxInOutside] *= -1.0;
208 rho += getRho(negVolVars);
209 deltaG[localDofIdx] += alpha_outside[idxInOutside];
213 rho /= numOutsideFaces + 1;
214 deltaG[localDofIdx] -= alpha_inside;
215 deltaG[localDofIdx] *= rho*curGlobalScvf.area();
219 rho = getRho(
elemVolVars()[curGlobalScvf.outsideScvIdx()]);
222 g[faceIdx] = alpha_inside*rho*curGlobalScvf.area();
227 for (
const auto& alpha : alpha_outside)
228 outsideG[faceIdx][i++] = alpha*rho*curGlobalScvf.area();
233 handle.CA().umv(deltaG, g);
236 using FaceVector =
typename IV::Traits::MatVecTraits::FaceVector;
239 handle.A().mv(deltaG, AG);
242 for (
const auto& localFaceData : iv.localFaceData())
245 if (!localFaceData.isOutsideFace())
248 const auto localScvIdx = localFaceData.ivLocalInsideScvIndex();
249 const auto localScvfIdx = localFaceData.ivLocalScvfIndex();
250 const auto idxInOutside = localFaceData.scvfLocalOutsideScvfIndex();
251 const auto& posLocalScv = iv.localScv(localScvIdx);
252 const auto& wijk = handle.omegas()[localScvfIdx][idxInOutside+1];
255 for (LocalIndexType localDir = 0; localDir < dim; localDir++)
258 const auto& curLocalScvf = iv.localScvf(posLocalScv.localScvfIndex(localDir));
259 if (!curLocalScvf.isDirichlet())
260 outsideG[localScvfIdx][idxInOutside] -= wijk[localDir]*AG[curLocalScvf.localDofIndex()];
268 const Problem* problemPtr_;
269 const FVElementGeometry* fvGeometryPtr_;
270 const ElementVolumeVariables* elemVolVarsPtr_;
A helper function for class member function introspection.
A class that contains helper functions as well as functionality which is common to different mpfa sch...
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:840
Defines the general interface of the local assembler classes for the assembly of the interaction volu...
Definition: localassemblerbase.hh:56
const FVElementGeometry & fvGeometry() const
Definition: localassemblerbase.hh:85
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 m...
Definition: localassemblerbase.hh:104
const ElementVolumeVariables & elemVolVars() const
Definition: localassemblerbase.hh:86
InteractionVolumeAssemblerBase(const Problem &problem, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars)
The constructor. Sets pointers to the objects required for a subsequent call to assemble().
Definition: localassemblerbase.hh:74
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: localassemblerbase.hh:121
void assembleGravity(DataHandle &handle, const IV &iv, const GetRho &getRho)
Assembles the gravitational flux contributions on the scvfs within an interaction volume.
Definition: localassemblerbase.hh:135
const Problem & problem() const
Definition: localassemblerbase.hh:84
A class that contains helper functions as well as functionality which is common to different mpfa sch...
Definition: localassemblerhelper.hh:46
static void resizeVector(Vector &v, size_type size)
resizes a vector to the given size (specialization for dynamic matrix type)
Definition: localassemblerhelper.hh:238