Commit ac2829a7 authored by Timo Koch's avatar Timo Koch
Browse files

[disc] Use extrusion instead of raw volume/area throughout code

parent 78d4d3f1
......@@ -30,6 +30,7 @@
#include <dumux/common/properties.hh>
#include <dumux/assembly/fvlocalresidual.hh>
#include <dumux/discretization/extrusion.hh>
namespace Dumux {
......@@ -48,7 +49,9 @@ class BoxLocalResidual : public FVLocalResidual<TypeTag>
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using ElementBoundaryTypes = GetPropType<TypeTag, Properties::ElementBoundaryTypes>;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
using Extrusion = Extrusion_t<GridGeometry>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
......@@ -114,7 +117,7 @@ public:
auto neumannFluxes = problem.neumann(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
// multiply neumann fluxes with the area and the extrusion factor
neumannFluxes *= scvf.area()*elemVolVars[scv].extrusionFactor();
neumannFluxes *= Extrusion::area(scvf)*elemVolVars[scv].extrusionFactor();
// only add fluxes to equations for which Neumann is set
for (int eqIdx = 0; eqIdx < NumEqVector::dimension; ++eqIdx)
......
......@@ -28,6 +28,7 @@
#include <dumux/common/reservedblockvector.hh>
#include <dumux/common/properties.hh>
#include <dumux/assembly/fvlocalresidual.hh>
#include <dumux/discretization/extrusion.hh>
namespace Dumux {
......@@ -46,7 +47,9 @@ class CCLocalResidual : public FVLocalResidual<TypeTag>
using ElementBoundaryTypes = GetPropType<TypeTag, Properties::ElementBoundaryTypes>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
using Extrusion = Extrusion_t<GridGeometry>;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
public:
......@@ -100,7 +103,7 @@ public:
// multiply neumann fluxes with the area and the extrusion factor
const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
neumannFluxes *= scvf.area()*elemVolVars[scv].extrusionFactor();
neumannFluxes *= Extrusion::area(scvf)*elemVolVars[scv].extrusionFactor();
flux += neumannFluxes;
}
......
......@@ -31,6 +31,7 @@
#include <dumux/common/timeloop.hh>
#include <dumux/common/reservedblockvector.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/extrusion.hh>
namespace Dumux {
......@@ -51,8 +52,9 @@ class FVLocalResidual
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using SubControlVolume = typename GridGeometry::SubControlVolume;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using Extrusion = Extrusion_t<GridGeometry>;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using ElementBoundaryTypes = GetPropType<TypeTag, Properties::ElementBoundaryTypes>;
using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
......@@ -114,7 +116,7 @@ public:
auto localScvIdx = scv.localDofIndex();
const auto& volVars = elemVolVars[scv];
storage[localScvIdx] = asImp().computeStorage(problem, scv, volVars);
storage[localScvIdx] *= scv.volume() * volVars.extrusionFactor();
storage[localScvIdx] *= Extrusion::volume(scv) * volVars.extrusionFactor();
}
return storage;
......@@ -307,7 +309,7 @@ public:
storage *= curVolVars.extrusionFactor();
storage -= prevStorage;
storage *= scv.volume();
storage *= Extrusion::volume(scv);
storage /= timeLoop_->timeStepSize();
residual[scv.localDofIndex()] += storage;
......@@ -336,7 +338,7 @@ public:
//! Compute source with the model specific storage residual
const auto& curVolVars = curElemVolVars[scv];
NumEqVector source = asImp().computeSource(problem, element, fvGeometry, curElemVolVars, scv);
source *= scv.volume()*curVolVars.extrusionFactor();
source *= Extrusion::volume(scv)*curVolVars.extrusionFactor();
//! subtract source from local rate (sign convention in user interface)
residual[scv.localDofIndex()] -= source;
......
......@@ -27,6 +27,7 @@
#include <dumux/common/timeloop.hh>
#include <dumux/common/properties.hh>
#include <dumux/discretization/extrusion.hh>
namespace Dumux {
......@@ -52,6 +53,7 @@ class StaggeredLocalResidual
using SubControlVolume = typename GridGeometry::SubControlVolume;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using FaceSubControlVolume = typename GridGeometry::Traits::FaceSubControlVolume;
using Extrusion = Extrusion_t<GridGeometry>;
using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>;
using CellCenterResidual = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>;
......@@ -123,7 +125,7 @@ public:
// subtract the source term from the local rate
auto source = asImp_().computeSourceForCellCenter(problem, element, fvGeometry, curElemVolVars, curElemFaceVars, scv);
source *= scv.volume()*curExtrusionFactor;
source *= Extrusion::volume(scv)*curExtrusionFactor;
residual -= source;
}
......@@ -169,7 +171,7 @@ public:
storage = std::move(curCCStorage);
storage -= std::move(prevCCStorage);
storage *= scv.volume();
storage *= Extrusion::volume(scv);
storage /= timeLoop_->timeStepSize();
residual += storage;
......@@ -243,13 +245,12 @@ public:
const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
const auto extrusionFactor = elemVolVars[scv].extrusionFactor();
// contruct staggered scv (half of the element)
auto scvCenter = scvf.center() - scv.center();
scvCenter *= 0.5;
FaceSubControlVolume faceScv(scvCenter, 0.5*scv.volume());
// construct staggered scv (half of the element)
auto faceScvCenter = scvf.center() + scv.center();
faceScvCenter *= 0.5;
FaceSubControlVolume faceScv(faceScvCenter, 0.5*scv.volume());
// multiply by 0.5 because we only consider half of a staggered control volume here
source *= faceScv.volume()*extrusionFactor;
source *= Extrusion::volume(faceScv)*extrusionFactor;
residual -= source;
}
......@@ -286,12 +287,12 @@ public:
const auto extrusionFactor = curElemVolVars[scv].extrusionFactor();
// contruct staggered scv (half of the element)
auto scvCenter = scvf.center() - scv.center();
scvCenter *= 0.5;
FaceSubControlVolume faceScv(scvCenter, 0.5*scv.volume());
// construct staggered scv (half of the element)
auto faceScvCenter = scvf.center() + scv.center();
faceScvCenter *= 0.5;
FaceSubControlVolume faceScv(faceScvCenter, 0.5*scv.volume());
storage *= faceScv.volume()*extrusionFactor;
storage *= Extrusion::volume(faceScv)*extrusionFactor;
storage /= timeLoop_->timeStepSize();
residual += storage;
......
......@@ -34,6 +34,7 @@
#include <dumux/common/parameters.hh>
#include <dumux/common/boundarytypes.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/extrusion.hh>
#include <dumux/assembly/initialsolution.hh>
......@@ -58,6 +59,7 @@ class FVProblem
using GridView = typename GridGeometry::GridView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using Extrusion = Extrusion_t<GridGeometry>;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
......@@ -440,7 +442,7 @@ public:
// Add the contributions to the dof source values
// We divide by the volume. In the local residual this will be multiplied with the same
// factor again. That's because the user specifies absolute values in kg/s.
const auto volume = scv.volume()*elemVolVars[scv].extrusionFactor();
const auto volume = Extrusion::volume(scv)*elemVolVars[scv].extrusionFactor();
for (const auto& ps : pointSourceMap_.at(key))
{
......
......@@ -34,6 +34,7 @@
#include <dune/common/reservedvector.hh>
#include <dumux/common/typetraits/isvalid.hh>
#include <dumux/discretization/extrusion.hh>
#include "localassemblerhelper.hh"
namespace Dumux {
......@@ -57,6 +58,7 @@ class InteractionVolumeAssemblerBase
using Problem = P;
using FVElementGeometry = EG;
using ElementVolumeVariables = EV;
using Extrusion = Extrusion_t<typename EG::GridGeometry>;
using Helper = InteractionVolumeAssemblerHelper;
template< class IV >
......@@ -212,20 +214,20 @@ class InteractionVolumeAssemblerBase
rho /= numOutsideFaces + 1;
deltaG[localDofIdx] -= alpha_inside;
deltaG[localDofIdx] *= rho*curGlobalScvf.area();
deltaG[localDofIdx] *= rho*Extrusion::area(curGlobalScvf);
}
// use density resulting from Dirichlet BCs
else
rho = getRho(elemVolVars()[curGlobalScvf.outsideScvIdx()]);
// add "inside" & "outside" alphas to gravity containers
g[faceIdx] = alpha_inside*rho*curGlobalScvf.area();
g[faceIdx] = alpha_inside*rho*Extrusion::area(curGlobalScvf);
if (isSurfaceGrid)
{
unsigned int i = 0;
for (const auto& alpha : alpha_outside)
outsideG[faceIdx][i++] = alpha*rho*curGlobalScvf.area();
outsideG[faceIdx][i++] = alpha*rho*Extrusion::area(curGlobalScvf);
}
}
......
......@@ -32,6 +32,7 @@
#include <dumux/common/parameters.hh>
#include <dumux/common/properties.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/extrusion.hh>
namespace Dumux {
......@@ -63,6 +64,7 @@ class BoxDarcysLaw
{
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using Extrusion = Extrusion_t<GridGeometry>;
using GridView = typename GridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
......@@ -116,7 +118,7 @@ public:
gradP.axpy(-rho, problem.spatialParams().gravity(scvf.center()));
// apply the permeability and return the flux
return -1.0*vtmv(scvf.unitOuterNormal(), K, gradP)*scvf.area();
return -1.0*vtmv(scvf.unitOuterNormal(), K, gradP)*Extrusion::area(scvf);
}
// compute transmissibilities ti for analytical jacobians
......@@ -145,7 +147,7 @@ public:
std::vector<Scalar> ti(fvGeometry.numScv());
for (const auto& scv : scvs(fvGeometry))
ti[scv.indexInElement()] =
-1.0*scvf.area()*vtmv(scvf.unitOuterNormal(), K, fluxVarCache.gradN(scv.indexInElement()));
-1.0*Extrusion::area(scvf)*vtmv(scvf.unitOuterNormal(), K, fluxVarCache.gradN(scv.indexInElement()));
return ti;
}
......
......@@ -27,6 +27,7 @@
#include <dumux/flux/effectivestresslaw.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/extrusion.hh>
namespace Dumux {
......@@ -42,6 +43,7 @@ class EffectiveStressLaw<StressType, GridGeometry, DiscretizationMethod::box>
{
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using Extrusion = Extrusion_t<GridGeometry>;
using GridView = typename GridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
......@@ -74,7 +76,7 @@ public:
ForceVector scvfForce(0.0);
sigma.mv(scvf.unitOuterNormal(), scvfForce);
scvfForce *= scvf.area();
scvfForce *= Extrusion::area(scvf);
return scvfForce;
}
......
......@@ -31,6 +31,7 @@
#include <dumux/common/math.hh>
#include <dumux/common/properties.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/extrusion.hh>
#include <dumux/flux/fickiandiffusioncoefficients.hh>
#include <dumux/flux/referencesystemformulation.hh>
......@@ -52,9 +53,11 @@ class FicksLawImplementation<TypeTag, DiscretizationMethod::box, referenceSystem
using Problem = GetPropType<TypeTag, Properties::Problem>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename GridGeometry::SubControlVolume;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using Extrusion = Extrusion_t<GridGeometry>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
using FluxVarCache = GetPropType<TypeTag, Properties::FluxVariablesCache>;
......@@ -151,7 +154,7 @@ public:
ti[compIdx].resize(fvGeometry.numScv());
for (auto&& scv : scvs(fvGeometry))
ti[compIdx][scv.indexInElement()] = -rho*vtmv(scvf.unitOuterNormal(), D, fluxVarCache.gradN(scv.indexInElement()))*scvf.area();
ti[compIdx][scv.indexInElement()] = -rho*vtmv(scvf.unitOuterNormal(), D, fluxVarCache.gradN(scv.indexInElement()))*Extrusion::area(scvf);
}
return ti;
......@@ -203,7 +206,7 @@ private:
Dune::FieldVector<Scalar, dimWorld> gradX(0.0);
for (auto&& scv : scvs(fvGeometry))
gradX.axpy(massOrMoleFraction(scv), fluxVarsCache.gradN(scv.indexInElement()));
return -1.0*preFactor*vtmv(scvf.unitOuterNormal(), D, gradX)*scvf.area();
return -1.0*preFactor*vtmv(scvf.unitOuterNormal(), D, gradX)*Extrusion::area(scvf);
}
};
......
......@@ -28,6 +28,7 @@
#include <dumux/common/math.hh>
#include <dumux/common/properties.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/extrusion.hh>
namespace Dumux {
......@@ -44,8 +45,10 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethod::box>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using Extrusion = Extrusion_t<GridGeometry>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
......@@ -85,7 +88,7 @@ public:
gradTemp.axpy(elemVolVars[scv].temperature(), fluxVarsCache.gradN(scv.indexInElement()));
// comute the heat conduction flux
return -1.0*vtmv(scvf.unitOuterNormal(), lambda, gradTemp)*scvf.area();
return -1.0*vtmv(scvf.unitOuterNormal(), lambda, gradTemp)*Extrusion::area(scvf);
}
};
......
......@@ -31,7 +31,7 @@
#include <dumux/common/properties.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/extrusion.hh>
namespace Dumux {
......@@ -48,8 +48,10 @@ class FouriersLawNonEquilibriumImplementation<TypeTag, DiscretizationMethod::box
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using Extrusion = Extrusion_t<GridGeometry>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
......@@ -111,7 +113,7 @@ public:
}
// comute the heat conduction flux
return -1.0*vtmv(scvf.unitOuterNormal(), lambda, gradTemp)*scvf.area();
return -1.0*vtmv(scvf.unitOuterNormal(), lambda, gradTemp)*Extrusion::area(scvf);
}
};
......
......@@ -29,6 +29,7 @@
#include <dumux/flux/hookeslaw.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/extrusion.hh>
namespace Dumux {
......@@ -42,7 +43,8 @@ template<class ScalarType, class GridGeometry>
class HookesLaw<ScalarType, GridGeometry, DiscretizationMethod::box>
{
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using Extrusion = Extrusion_t<GridGeometry>;
using GridView = typename GridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
......@@ -74,7 +76,7 @@ public:
ForceVector scvfForce(0.0);
sigma.mv(scvf.unitOuterNormal(), scvfForce);
scvfForce *= scvf.area();
scvfForce *= Extrusion::area(scvf);
return scvfForce;
}
......
......@@ -31,6 +31,7 @@
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/extrusion.hh>
#include <dumux/flux/maxwellstefandiffusioncoefficients.hh>
#include <dumux/flux/fluxvariablescaching.hh>
......@@ -53,8 +54,10 @@ class MaxwellStefansLawImplementation<TypeTag, DiscretizationMethod::box, refere
using Problem = GetPropType<TypeTag, Properties::Problem>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using Extrusion = Extrusion_t<GridGeometry>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
......@@ -129,7 +132,7 @@ public:
normalX[compIdx] = gradX *scvf.unitOuterNormal();
}
reducedDiffusionMatrix.solve(reducedFlux,normalX);
reducedFlux *= -1.0*rho*scvf.area();
reducedFlux *= -1.0*rho*Extrusion::area(scvf);
for (int compIdx = 0; compIdx < numComponents-1; compIdx++)
{
......
......@@ -29,6 +29,7 @@
#include <dumux/common/properties.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/extrusion.hh>
#include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh>
namespace Dumux {
......@@ -131,6 +132,7 @@ class CCTpfaDarcysLaw<ScalarType, GridGeometry, /*isNetwork*/ false>
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename GridGeometry::SubControlVolume;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using Extrusion = Extrusion_t<GridGeometry>;
using GridView = typename GridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
......@@ -184,7 +186,7 @@ class CCTpfaDarcysLaw<ScalarType, GridGeometry, /*isNetwork*/ false>
//! compute alpha := n^T*K*g
const auto alpha_inside = vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g)*insideVolVars.extrusionFactor();
Scalar flux = tij*(pInside - pOutside) + rho*scvf.area()*alpha_inside;
Scalar flux = tij*(pInside - pOutside) + rho*Extrusion::area(scvf)*alpha_inside;
//! On interior faces we have to add K-weighted gravitational contributions
if (!scvf.boundary())
......@@ -233,7 +235,7 @@ class CCTpfaDarcysLaw<ScalarType, GridGeometry, /*isNetwork*/ false>
// on the boundary (dirichlet) we only need ti
if (scvf.boundary())
tij = scvf.area()*ti;
tij = Extrusion::area(scvf)*ti;
// otherwise we compute a tpfa harmonic mean
else
......@@ -252,7 +254,7 @@ class CCTpfaDarcysLaw<ScalarType, GridGeometry, /*isNetwork*/ false>
if (ti*tj <= 0.0)
tij = 0;
else
tij = scvf.area()*(ti * tj)/(ti + tj);
tij = Extrusion::area(scvf)*(ti * tj)/(ti + tj);
}
return tij;
......@@ -285,6 +287,7 @@ class CCTpfaDarcysLaw<ScalarType, GridGeometry, /*isNetwork*/ true>
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename GridGeometry::SubControlVolume;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using Extrusion = Extrusion_t<GridGeometry>;
using GridView = typename GridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
......@@ -367,7 +370,7 @@ public:
Scalar sumPTi(tij*pInside);
// add inside gravitational contribution
sumPTi += rho*scvf.area()
sumPTi += rho*Extrusion::area(scvf)
*insideVolVars.extrusionFactor()
*vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g);
......@@ -381,7 +384,7 @@ public:
sumPTi += outsideFluxVarsCache.advectionTij()*outsideVolVars.pressure(phaseIdx);
// add outside gravitational contribution
sumPTi += rho*scvf.area()
sumPTi += rho*Extrusion::area(scvf)
*outsideVolVars.extrusionFactor()
*vtmv(flippedScvf.unitOuterNormal(), outsideVolVars.permeability(), g);
}
......@@ -392,7 +395,7 @@ public:
//! precompute alpha := n^T*K*g
const auto alpha_inside = vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g)*insideVolVars.extrusionFactor();
Scalar flux = tij*(pInside - pOutside) + scvf.area()*rho*alpha_inside;
Scalar flux = tij*(pInside - pOutside) + Extrusion::area(scvf)*rho*alpha_inside;
//! On interior faces with one neighbor we have to add K-weighted gravitational contributions
if (!scvf.boundary() && scvf.numOutsideScvs() == 1)
......@@ -464,7 +467,7 @@ public:
// for the boundary (dirichlet) or at branching points we only need ti
if (scvf.boundary() || scvf.numOutsideScvs() > 1)
tij = scvf.area()*ti;
tij = Extrusion::area(scvf)*ti;
// otherwise we compute a tpfa harmonic mean
else
......@@ -483,7 +486,7 @@ public:
if (ti*tj <= 0.0)
tij = 0;
else
tij = scvf.area()*(ti * tj)/(ti + tj);
tij = Extrusion::area(scvf)*(ti * tj)/(ti + tj);
}
return tij;
......
......@@ -30,6 +30,7 @@
#include <dumux/common/properties.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/extrusion.hh>
#include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh>
#include <dumux/flux/fickiandiffusioncoefficients.hh>
......@@ -51,9 +52,11 @@ class FicksLawImplementation<TypeTag, DiscretizationMethod::cctpfa, referenceSys
using Implementation = FicksLawImplementation<TypeTag, DiscretizationMethod::cctpfa, referenceSystem>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using Extrusion = Extrusion_t<GridGeometry>;
using GridView = typename GridGeometry::GridView;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using Element = typename GridView::template Codim<0>::Entity;
using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
......@@ -198,7 +201,7 @@ public:
// for the boundary (dirichlet) or at branching points we only need ti
Scalar tij;
if (scvf.boundary() || scvf.numOutsideScvs() > 1)
tij = scvf.area()*ti;
tij = Extrusion::area(scvf)*ti;
// otherwise we compute a tpfa harmonic mean
else
......@@ -222,7 +225,7 @@ public:
if (ti*tj <= 0.0)
tij = 0;
else
tij = scvf.area()*(ti * tj)/(ti + tj);
tij = Extrusion::area(scvf)*(ti * tj)/(ti + tj);
}
return tij;
......
......@@ -33,6 +33,7 @@
#include <dumux/common/typetraits/typetraits.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/extrusion.hh>
#include <dumux/flux/cctpfa/darcyslaw.hh>
namespace Dumux {
......@@ -142,6 +143,7 @@ class CCTpfaForchheimersLaw<ScalarType, GridGeometry, /*isNetwork*/ false>
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename GridGeometry::SubControlVolume;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using Extrusion = Extrusion_t<GridGeometry>;