Commit bcccf922 authored by Dennis Gläser's avatar Dennis Gläser Committed by Timo Koch
Browse files

[disc] reintroduce indexInElement() funciton with different meaning

Now, the indexInElement() function returns the index of the scv in the
element-local set of scvs. For schemes introducint additional scvs on
some elements this does not coincide with localDofIndex(). This commit
furthermore introduces the correct usage of the function throughout Dumux.
parent 4bb2a586
......@@ -407,7 +407,7 @@ public:
const SubControlVolume &scv) const
{
NumEqVector source(0);
auto scvIdx = scv.localDofIndex();
auto scvIdx = scv.indexInElement();
auto key = std::make_pair(fvGridGeometry_->elementMapper().index(element), scvIdx);
if (pointSourceMap_.count(key))
{
......
......@@ -303,7 +303,7 @@ public:
for (auto&& scv : scvs(fvGeometry))
{
if (intersectsPointGeometry(globalPos, scv.geometry()))
scvIndices.push_back(scv.localDofIndex());
scvIndices.push_back(scv.indexInElement());
}
// for all scvs that where tested positiv add the point sources
// to the element/scv to point source map
......
......@@ -98,10 +98,10 @@ public:
const auto& volVars = elemVolVars[scv];
if (enableGravity)
rho += volVars.density(phaseIdx)*shapeValues[scv.localDofIndex()][0];
rho += volVars.density(phaseIdx)*shapeValues[scv.indexInElement()][0];
// the global shape function gradient
gradP.axpy(volVars.pressure(phaseIdx), fluxVarCache.gradN(scv.localDofIndex()));
gradP.axpy(volVars.pressure(phaseIdx), fluxVarCache.gradN(scv.indexInElement()));
}
if (enableGravity)
......@@ -135,8 +135,8 @@ public:
std::vector<Scalar> ti(fvGeometry.numScv());
for (const auto& scv : scvs(fvGeometry))
ti[scv.localDofIndex()] =
-1.0*scvf.area()*vtmv(scvf.unitOuterNormal(), K, fluxVarCache.gradN(scv.localDofIndex()));
ti[scv.indexInElement()] =
-1.0*scvf.area()*vtmv(scvf.unitOuterNormal(), K, fluxVarCache.gradN(scv.indexInElement()));
return ti;
}
......
......@@ -63,7 +63,7 @@ public:
template<class SubControlVolume, typename std::enable_if_t<!std::is_integral<SubControlVolume>::value, int> = 0>
const VolumeVariables& operator [](const SubControlVolume& scv) const
{ return gridVolVars().volVars(eIdx_, scv.localDofIndex()); }
{ return gridVolVars().volVars(eIdx_, scv.indexInElement()); }
// For compatibility reasons with the case of not storing the vol vars.
// function to be called before assembling an element, preparing the vol vars within the stencil
......@@ -133,7 +133,7 @@ public:
// resize volume variables to the required size
volumeVariables_.resize(fvGeometry.numScv());
for (auto&& scv : scvs(fvGeometry))
volumeVariables_[scv.localDofIndex()].update(elemSol, gridVolVars().problem(), element, scv);
volumeVariables_[scv.indexInElement()].update(elemSol, gridVolVars().problem(), element, scv);
}
const VolumeVariables& operator [](std::size_t scvIdx) const
......@@ -144,11 +144,11 @@ public:
template<class SubControlVolume, typename std::enable_if_t<!std::is_integral<SubControlVolume>::value, int> = 0>
const VolumeVariables& operator [](const SubControlVolume& scv) const
{ return volumeVariables_[scv.localDofIndex()]; }
{ return volumeVariables_[scv.indexInElement()]; }
template<class SubControlVolume, typename std::enable_if_t<!std::is_integral<SubControlVolume>::value, int> = 0>
VolumeVariables& operator [](const SubControlVolume& scv)
{ return volumeVariables_[scv.localDofIndex()]; }
{ return volumeVariables_[scv.indexInElement()]; }
//! The global volume variables object we are a restriction of
const GridVolumeVariables& gridVolVars() const
......
......@@ -93,7 +93,7 @@ public:
// density interpolation
Scalar rho(0.0);
for (auto&& scv : scvs(fvGeometry))
rho += elemVolVars[scv].molarDensity(phaseIdx)*shapeValues[scv.localDofIndex()][0];
rho += elemVolVars[scv].molarDensity(phaseIdx)*shapeValues[scv.indexInElement()][0];
for (int compIdx = 0; compIdx < numComponents; compIdx++)
{
......@@ -119,7 +119,7 @@ public:
// the mole/mass fraction gradient
GlobalPosition gradX(0.0);
for (auto&& scv : scvs(fvGeometry))
gradX.axpy(elemVolVars[scv].moleFraction(phaseIdx, compIdx), fluxVarsCache.gradN(scv.localDofIndex()));
gradX.axpy(elemVolVars[scv].moleFraction(phaseIdx, compIdx), fluxVarsCache.gradN(scv.indexInElement()));
// compute the diffusive flux
componentFlux[compIdx] = -1.0*rho*vtmv(scvf.unitOuterNormal(), D, gradX)*scvf.area();
......@@ -142,7 +142,7 @@ public:
Scalar rho(0.0);
const auto& shapeValues = fluxVarCache.shapeValues();
for (auto&& scv : scvs(fvGeometry))
rho += elemVolVars[scv].molarDensity(phaseIdx)*shapeValues[scv.localDofIndex()][0];
rho += elemVolVars[scv].molarDensity(phaseIdx)*shapeValues[scv.indexInElement()][0];
const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
......@@ -171,7 +171,7 @@ public:
ti[compIdx].resize(fvGeometry.numScv());
for (auto&& scv : scvs(fvGeometry))
ti[compIdx][scv.localDofIndex()] = -rho*vtmv(scvf.unitOuterNormal(), D, fluxVarCache.gradN(scv.localDofIndex()))*scvf.area();
ti[compIdx][scv.indexInElement()] = -rho*vtmv(scvf.unitOuterNormal(), D, fluxVarCache.gradN(scv.indexInElement()))*scvf.area();
}
return ti;
......
......@@ -96,7 +96,7 @@ public:
// compute the temperature gradient with the shape functions
GlobalPosition gradTemp(0.0);
for (auto&& scv : scvs(fvGeometry))
gradTemp.axpy(elemVolVars[scv].temperature(), fluxVarsCache.gradN(scv.localDofIndex()));
gradTemp.axpy(elemVolVars[scv].temperature(), fluxVarsCache.gradN(scv.indexInElement()));
// comute the heat conduction flux
return -1.0*vtmv(scvf.unitOuterNormal(), lambda, gradTemp)*scvf.area();
......
......@@ -122,9 +122,9 @@ public:
{
// compute the temperature gradient with the shape functions
if (phaseIdx < numEnergyEqFluid)
gradTemp.axpy(elemVolVars[scv].temperature(phaseIdx), fluxVarsCache.gradN(scv.localDofIndex()));
gradTemp.axpy(elemVolVars[scv].temperature(phaseIdx), fluxVarsCache.gradN(scv.indexInElement()));
else
gradTemp.axpy(elemVolVars[scv].temperatureSolid(), fluxVarsCache.gradN(scv.localDofIndex()));
gradTemp.axpy(elemVolVars[scv].temperatureSolid(), fluxVarsCache.gradN(scv.indexInElement()));
}
// comute the heat conduction flux
......
......@@ -88,17 +88,17 @@ public:
// update the volvars of the element
volumeVariables_[eIdx].resize(fvGeometry.numScv());
for (auto&& scv : scvs(fvGeometry))
volumeVariables_[eIdx][scv.localDofIndex()].update(elemSol, problem(), element, scv);
volumeVariables_[eIdx][scv.indexInElement()].update(elemSol, problem(), element, scv);
}
}
template<class SubControlVolume, typename std::enable_if_t<!std::is_integral<SubControlVolume>::value, int> = 0>
const VolumeVariables& volVars(const SubControlVolume& scv) const
{ return volumeVariables_[scv.elementIndex()][scv.localDofIndex()]; }
{ return volumeVariables_[scv.elementIndex()][scv.indexInElement()]; }
template<class SubControlVolume, typename std::enable_if_t<!std::is_integral<SubControlVolume>::value, int> = 0>
VolumeVariables& volVars(const SubControlVolume& scv)
{ return volumeVariables_[scv.elementIndex()][scv.localDofIndex()]; }
{ return volumeVariables_[scv.elementIndex()][scv.indexInElement()]; }
const VolumeVariables& volVars(const std::size_t eIdx, const std::size_t scvIdx) const
{ return volumeVariables_[eIdx][scvIdx]; }
......
......@@ -103,14 +103,14 @@ public:
const auto& volVars = elemVolVars[scv];
// density interpolation
rho += volVars.molarDensity(phaseIdx)*shapeValues[scv.localDofIndex()][0];
rho += volVars.molarDensity(phaseIdx)*shapeValues[scv.indexInElement()][0];
// the ansatz function gradient
jacInvT.mv(shapeJacobian[scv.localDofIndex()][0], gradN[scv.localDofIndex()]);
jacInvT.mv(shapeJacobian[scv.localDofIndex()][0], gradN[scv.indexInElement()]);
//interpolate the mole fraction for the diffusion matrix
for (int compIdx = 0; compIdx < numComponents; compIdx++)
{
moleFrac[compIdx] += volVars.moleFraction(phaseIdx, compIdx)*shapeValues[scv.localDofIndex()][0];
moleFrac[compIdx] += volVars.moleFraction(phaseIdx, compIdx)*shapeValues[scv.indexInElement()][0];
}
}
......@@ -124,7 +124,7 @@ public:
const auto& volVars = elemVolVars[scv];
// the mole/mass fraction gradient
gradX.axpy(volVars.moleFraction(phaseIdx, compIdx), gradN[scv.localDofIndex()]);
gradX.axpy(volVars.moleFraction(phaseIdx, compIdx), gradN[scv.indexInElement()]);
}
normalX[compIdx] = gradX *scvf.unitOuterNormal();
......
......@@ -151,6 +151,13 @@ public:
return localDofIdx_;
}
//! The element-local index of this scv.
//! For the standard box scheme this is the local dof index.
LocalIndexType indexInElement() const
{
return localDofIdx_;
}
//! The index of the dof this scv is embedded in
GridIndexType dofIndex() const
{
......
......@@ -144,6 +144,13 @@ public:
return 0;
}
//! The element-local index of this scv.
//! In cell-centered schemes there is always only one scv per element.
LocalIndexType indexInElement() const
{
return 0;
}
// The position of the dof this scv is embedded in
const GlobalPosition& dofPosition() const
{
......
......@@ -195,7 +195,7 @@ public:
/ curElemVolVars[scvf.insideScvIdx()].viscosity();
for (const auto& scv : scvs(fvGeometry))
{
auto d = up*ti[scv.localDofIndex()];
auto d = up*ti[scv.indexInElement()];
A[insideScv.dofIndex()][scv.dofIndex()][conti0EqIdx][pressureIdx] += d;
A[outsideScv.dofIndex()][scv.dofIndex()][conti0EqIdx][pressureIdx] -= d;
}
......
......@@ -345,7 +345,7 @@ public:
for (const auto& scvJ : scvs(fvGeometry))
{
const auto globalJ = scvJ.dofIndex();
const auto localJ = scvJ.localDofIndex();
const auto localJ = scvJ.indexInElement();
// the transmissibily associated with the scvJ
const auto tj = ti[localJ];
......
......@@ -99,8 +99,7 @@ public:
// compute the gradN at for every scv/dof
gradN_.resize(fvGeometry.numScv());
for (const auto& scv: scvs(fvGeometry))
jacInvT_.mv(shapeJacobian_[scv.localDofIndex()][0], gradN_[scv.localDofIndex()]);
jacInvT_.mv(shapeJacobian_[scv.localDofIndex()][0], gradN_[scv.indexInElement()]);
}
const std::vector<ShapeJacobian>& shapeJacobian() const
......
......@@ -71,10 +71,9 @@ public:
for (auto&& scv : scvs(fvGeometry))
{
//here we need to set the constraints of the mpnc model into the residual
const auto localScvIdx = scv.localDofIndex();
// here we need to set the constraints of the mpnc model into the residual
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
residual[localScvIdx][phase0NcpIdx + phaseIdx] = elemVolVars[scv].phaseNcp(phaseIdx);
residual[scv.localDofIndex()][phase0NcpIdx + phaseIdx] = elemVolVars[scv].phaseNcp(phaseIdx);
}
return residual;
......
......@@ -297,8 +297,8 @@ public:
for (const auto& scv : scvs(fvGeometry))
{
// diffusive term
const auto diffDeriv = useMoles ? ti[compIdx][scv.localDofIndex()]
: ti[compIdx][scv.localDofIndex()]*FluidSystem::molarMass(compIdx);
const auto diffDeriv = useMoles ? ti[compIdx][scv.indexInElement()]
: ti[compIdx][scv.indexInElement()]*FluidSystem::molarMass(compIdx);
A[insideScv.dofIndex()][scv.dofIndex()][compIdx][compIdx] += diffDeriv;
A[outsideScv.dofIndex()][scv.dofIndex()][compIdx][compIdx] -= diffDeriv;
}
......
......@@ -239,7 +239,7 @@ public:
// calculate the subcontrolvolume velocity by the Piola transformation
Dune::FieldVector<CoordScalar, dimWorld> scvVelocity(0);
jacobianT2.mtv(scvVelocities[scv.localDofIndex()], scvVelocity);
jacobianT2.mtv(scvVelocities[scv.indexInElement()], scvVelocity);
scvVelocity /= geometry.integrationElement(localPos)*cellNum_[vIdxGlobal];
// add up the wetting phase subcontrolvolume velocities for each vertex
velocity[vIdxGlobal] += scvVelocity;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment