25#ifndef DUMUX_DISCRETIZATION_CC_MPFA_O_LOCAL_ASSEMBLER_HH
26#define DUMUX_DISCRETIZATION_CC_MPFA_O_LOCAL_ASSEMBLER_HH
45template<
class P,
class EG,
class EV >
53 using Scalar =
typename IV::Traits::MatVecTraits::FaceVector::value_type;
57 using ParentType::ParentType;
73 template<
class DataHandle,
class IV,
class TensorFunc >
74 void assembleMatrices(DataHandle& handle, IV& iv,
const TensorFunc& getT, Scalar<IV> wijZeroThresh = 0.0)
76 const auto zeroRows = assembleLocalMatrices_(handle.A(), handle.AB(), handle.CA(), handle.T(), handle.omegas(), iv, getT, wijZeroThresh);
79 if (iv.numUnknowns() > 0)
82 auto& B = handle.AB();
88 for (
const auto& zeroRowIndices : zeroRows)
90 const auto zeroRowDofIdx = zeroRowIndices.first;
92 row[zeroRowDofIdx] = 0.0;
93 A[zeroRowDofIdx] = 0.0;
94 A[zeroRowDofIdx][zeroRowDofIdx] = 1.0;
95 B[zeroRowDofIdx] = 0.0;
103 for (
const auto& zeroRowIndices : zeroRows)
105 const auto faceIdx = zeroRowIndices.second;
106 A[zeroRowIndices.first] = 0.0;
107 handle.CA()[faceIdx] = 0.0;
108 handle.T()[faceIdx] = 0.0;
111 static constexpr int dim = IV::Traits::GridView::dimension;
112 static constexpr int dimWorld = IV::Traits::GridView::dimensionworld;
113 if constexpr (dim < dimWorld)
114 std::for_each( handle.tijOutside()[faceIdx].begin(),
115 handle.tijOutside()[faceIdx].end(),
116 [] (
auto& outsideTij) { outsideTij = 0.0; } );
129 template<
class DataHandle,
class IV,
class GetU >
130 void assembleU(DataHandle& handle,
const IV& iv,
const GetU& getU)
132 auto& u = handle.uj();
136 typename IV::Traits::IndexSet::LocalIndexType i = 0;
137 for (; i < iv.numScvs(); i++)
138 u[i] = getU( this->
elemVolVars()[iv.localScv(i).gridScvIndex()] );
139 for (
const auto& data : iv.dirichletData())
140 u[i++] = getU( this->
elemVolVars()[data.volVarIndex()] );
174 template<
class IV,
class TensorFunc >
175 auto assembleLocalMatrices_(
typename IV::Traits::MatVecTraits::AMatrix& A,
176 typename IV::Traits::MatVecTraits::BMatrix& B,
177 typename IV::Traits::MatVecTraits::CMatrix& C,
178 typename IV::Traits::MatVecTraits::DMatrix& D,
179 typename IV::Traits::MatVecTraits::OmegaStorage& wijk,
180 IV& iv,
const TensorFunc& getT,
181 Scalar<IV> wijZeroThresh)
183 using LocalIndexType =
typename IV::Traits::IndexSet::LocalIndexType;
184 static constexpr int dim = IV::Traits::GridView::dimension;
185 static constexpr int dimWorld = IV::Traits::GridView::dimensionworld;
187 std::vector< std::pair<LocalIndexType, LocalIndexType> > faceMarkers;
192 if (iv.numUnknowns() == 0)
198 for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx)
200 const auto& curLocalScvf = iv.localScvf(faceIdx);
201 const auto& curGlobalScvf = this->
fvGeometry().scvf(curLocalScvf.gridScvfIndex());
202 const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices();
205 const auto& posLocalScv = iv.localScv(neighborScvIndices[0]);
206 const auto& posGlobalScv = this->
fvGeometry().scv(posLocalScv.gridScvIndex());
207 const auto& posVolVars = this->
elemVolVars()[posGlobalScv];
208 const auto tensor = getT(posVolVars);
212 wijk[faceIdx][0] = computeMpfaTransmissibility<EG>(this->
fvGeometry(), posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor());
214 const auto posScvLocalDofIdx = posLocalScv.localDofIndex();
215 for (LocalIndexType localDir = 0; localDir < dim; localDir++)
217 const auto& otherLocalScvf = iv.localScvf( posLocalScv.localScvfIndex(localDir) );
218 const auto otherLocalDofIdx = otherLocalScvf.localDofIndex();
219 D[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir];
220 D[faceIdx][posScvLocalDofIdx] += wijk[faceIdx][0][localDir];
232 for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx)
234 const auto& curLocalScvf = iv.localScvf(faceIdx);
235 const auto& curGlobalScvf = this->
fvGeometry().scvf(curLocalScvf.gridScvfIndex());
236 const auto curIsDirichlet = curLocalScvf.isDirichlet();
237 const auto curLocalDofIdx = curLocalScvf.localDofIndex();
240 const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices();
241 const auto& posLocalScv = iv.localScv(neighborScvIndices[0]);
242 const auto& posGlobalScv = this->
fvGeometry().scv(posLocalScv.gridScvIndex());
243 const auto& posVolVars = this->
elemVolVars()[posGlobalScv];
244 const auto tensor = getT(posVolVars);
248 wijk[faceIdx][0] = computeMpfaTransmissibility<EG>(this->
fvGeometry(), posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor());
251 bool insideZeroWij =
false;
254 for (
unsigned int localDir = 0; localDir < dim; localDir++)
256 const auto& otherLocalScvf = iv.localScvf( posLocalScv.localScvfIndex(localDir) );
257 const auto otherLocalDofIdx = otherLocalScvf.localDofIndex();
259 if (otherLocalDofIdx == curLocalDofIdx)
261 if (abs(wijk[faceIdx][0][localDir]) <= wijZeroThresh)
265 insideZeroWij =
true;
266 faceMarkers.emplace_back( std::make_pair(curLocalDofIdx, faceIdx) );
273 if (!otherLocalScvf.isDirichlet())
275 C[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir];
277 A[curLocalDofIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir];
282 D[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir];
284 B[curLocalDofIdx][otherLocalDofIdx] += wijk[faceIdx][0][localDir];
288 const auto posScvLocalDofIdx = posLocalScv.localDofIndex();
289 D[faceIdx][posScvLocalDofIdx] += wijk[faceIdx][0][localDir];
292 B[curLocalDofIdx][posScvLocalDofIdx] -= wijk[faceIdx][0][localDir];
296 if (!curGlobalScvf.boundary())
299 for (
unsigned int idxInOutside = 0; idxInOutside < curGlobalScvf.numOutsideScvs(); ++idxInOutside)
301 const auto idxOnScvf = idxInOutside+1;
302 const auto& negLocalScv = iv.localScv( neighborScvIndices[idxOnScvf] );
303 const auto& negGlobalScv = this->
fvGeometry().scv(negLocalScv.gridScvIndex());
304 const auto& negVolVars = this->
elemVolVars()[negGlobalScv];
305 const auto negTensor = getT(negVolVars);
308 const auto& scvf = dim < dimWorld ? this->
fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside)
310 wijk[faceIdx][idxOnScvf] = computeMpfaTransmissibility<EG>(this->
fvGeometry(), negLocalScv, scvf, negTensor, negVolVars.extrusionFactor());
314 wijk[faceIdx][idxOnScvf] *= -1.0;
317 for (
int localDir = 0; localDir < dim; localDir++)
319 const auto otherLocalScvfIdx = negLocalScv.localScvfIndex(localDir);
320 const auto& otherLocalScvf = iv.localScvf(otherLocalScvfIdx);
321 const auto otherLocalDofIdx = otherLocalScvf.localDofIndex();
324 if (otherLocalDofIdx == curLocalDofIdx && !insideZeroWij)
325 if (abs(wijk[faceIdx][idxOnScvf][localDir]) <= wijZeroThresh)
326 faceMarkers.emplace_back( std::make_pair(curLocalDofIdx, faceIdx) );
328 if (!otherLocalScvf.isDirichlet())
329 A[curLocalDofIdx][otherLocalDofIdx] += wijk[faceIdx][idxOnScvf][localDir];
331 B[curLocalDofIdx][otherLocalDofIdx] -= wijk[faceIdx][idxOnScvf][localDir];
334 B[curLocalDofIdx][negLocalScv.localDofIndex()] += wijk[faceIdx][idxOnScvf][localDir];
Defines the general interface of classes used for the assembly of the local systems of equations invo...
A class that contains helper functions as well as functionality which is common to different mpfa sch...
Adaption of the non-isothermal two-phase two-component flow model to problems with CO2.
Definition: adapt.hh:29
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
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: discretization/cellcentered/mpfa/omethod/localassembler.hh:48
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: discretization/cellcentered/mpfa/omethod/localassembler.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: discretization/cellcentered/mpfa/omethod/localassembler.hh:130
This file contains free functions to evaluate the transmissibilities associated with flux evaluations...