From 52b82dc109a26592a88c7599cfa37558ef3e9e79 Mon Sep 17 00:00:00 2001
From: DennisGlaeser <dennis.glaeser@iws.uni-stuttgart.de>
Date: Wed, 8 Feb 2017 14:35:32 +0100
Subject: [PATCH] [mpfa] pass element to facetVolVars() method
---
.../cellcentered/mpfa/darcyslaw.hh | 21 +++++++--------
.../cellcentered/mpfa/fickslaw.hh | 26 +++++++++----------
.../cellcentered/mpfa/fourierslaw.hh | 15 +++++------
.../cellcentered/mpfa/interiorboundarydata.hh | 3 +--
dumux/discretization/upwindscheme.hh | 4 +--
5 files changed, 32 insertions(+), 37 deletions(-)
diff --git a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
index 485c63e96f..e9f35e251e 100644
--- a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
@@ -173,20 +173,17 @@ public:
const auto& tij = fluxVarsCache.advectionTij();
const bool isInteriorBoundary = enableInteriorBoundaries && fluxVarsCache.isInteriorBoundary();
- // For interior Neumann boundaries when using Tpfa for Neumann boundary conditions, we simply
- // return the user-specified flux. We assume phaseIdx = eqIdx here.
+ // For interior Neumann boundaries when using Tpfa on boundaries, return the user-specified flux
+ // We assume phaseIdx = eqIdx here.
if (isInteriorBoundary
&& useTpfaBoundary
&& fluxVarsCache.interiorBoundaryDataSelf().faceType() == MpfaFaceTypes::interiorNeumann)
return scvf.area()*
elemVolVars[scvf.insideScvIdx()].extrusionFactor()*
- problem.neumann(element,
- fvGeometry,
- elemVolVars,
- scvf)[phaseIdx];
+ problem.neumann(element, fvGeometry, elemVolVars, scvf)[phaseIdx];
// Calculate the interface density for gravity evaluation
- const auto rho = Implementation::interpolateDensity(fvGeometry, elemVolVars, scvf, fluxVarsCache, phaseIdx, isInteriorBoundary);
+ const auto rho = Implementation::interpolateDensity(element, fvGeometry, elemVolVars, scvf, fluxVarsCache, phaseIdx, isInteriorBoundary);
// calculate Tij*pj
Scalar flux(0.0);
@@ -214,13 +211,14 @@ public:
return useTpfaBoundary ? flux : flux + fluxVarsCache.advectionNeumannFlux(phaseIdx);
// Handle interior boundaries
- flux += Implementation::computeInteriorBoundaryContribution(problem, fvGeometry, fluxVarsCache, phaseIdx, rho);
+ flux += Implementation::computeInteriorBoundaryContribution(problem, element, fvGeometry, fluxVarsCache, phaseIdx, rho);
// return overall resulting flux
return useTpfaBoundary ? flux : flux + fluxVarsCache.advectionNeumannFlux(phaseIdx);
}
- static Scalar interpolateDensity(const FVElementGeometry& fvGeometry,
+ static Scalar interpolateDensity(const Element& element,
+ const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const FluxVariablesCache& fluxVarsCache,
@@ -238,7 +236,7 @@ public:
{
const auto& data = fluxVarsCache.interiorBoundaryDataSelf();
if (data.faceType() == MpfaFaceTypes::interiorDirichlet)
- return data.facetVolVars(fvGeometry).density(phaseIdx);
+ return data.facetVolVars(element, fvGeometry).density(phaseIdx);
}
// use arithmetic mean of the densities around the scvf
@@ -255,6 +253,7 @@ public:
}
static Scalar computeInteriorBoundaryContribution(const Problem& problem,
+ const Element& element,
const FVElementGeometry& fvGeometry,
const FluxVariablesCache& fluxVarsCache,
unsigned int phaseIdx, Scalar rho)
@@ -274,7 +273,7 @@ public:
{
if (data.faceType() == MpfaFaceTypes::interiorDirichlet)
{
- Scalar h = data.facetVolVars(fvGeometry).pressure(phaseIdx);
+ Scalar h = data.facetVolVars(element, fvGeometry).pressure(phaseIdx);
if (gravity)
{
diff --git a/dumux/discretization/cellcentered/mpfa/fickslaw.hh b/dumux/discretization/cellcentered/mpfa/fickslaw.hh
index 70bfc99d1b..710a58cdb6 100644
--- a/dumux/discretization/cellcentered/mpfa/fickslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/fickslaw.hh
@@ -188,14 +188,11 @@ public:
&& fluxVarsCache.interiorBoundaryDataSelf().faceType() == MpfaFaceTypes::interiorNeumann)
return scvf.area()*
elemVolVars[scvf.insideScvIdx()].extrusionFactor()*
- problem.neumann(element,
- fvGeometry,
- elemVolVars,
- scvf)[compIdx];
+ problem.neumann(element, fvGeometry, elemVolVars, scvf)[compIdx];
// get the scaling factor for the effective diffusive fluxes
- const auto effFactor = Implementation::computeEffectivityFactor(fvGeometry, elemVolVars, scvf, fluxVarsCache, phaseIdx, isInteriorBoundary);
+ const auto effFactor = Implementation::computeEffectivityFactor(element, fvGeometry, elemVolVars, scvf, fluxVarsCache, phaseIdx, isInteriorBoundary);
// if factor is zero, the flux will end up zero anyway
if (effFactor == 0.0)
@@ -209,7 +206,7 @@ public:
{ return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); };
// calculate the density at the interface
- const auto rho = Implementation::interpolateDensity(fvGeometry, elemVolVars, scvf, fluxVarsCache, getRho, isInteriorBoundary);
+ const auto rho = Implementation::interpolateDensity(element, fvGeometry, elemVolVars, scvf, fluxVarsCache, getRho, isInteriorBoundary);
// calculate Tij*xj
Scalar flux(0.0);
@@ -222,14 +219,15 @@ public:
return useTpfaBoundary ? flux*rho*effFactor : flux*rho*effFactor + fluxVarsCache.componentNeumannFlux(compIdx);
// Handle interior boundaries
- flux += Implementation::computeInteriorBoundaryContribution(fvGeometry, fluxVarsCache, getX, phaseIdx, compIdx);
+ flux += Implementation::computeInteriorBoundaryContribution(element, fvGeometry, fluxVarsCache, getX, phaseIdx, compIdx);
// return overall resulting flux
return useTpfaBoundary ? flux*rho*effFactor : flux*rho*effFactor + fluxVarsCache.componentNeumannFlux(compIdx);
}
template<typename GetRhoFunction>
- static Scalar interpolateDensity(const FVElementGeometry& fvGeometry,
+ static Scalar interpolateDensity(const Element& element,
+ const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const FluxVariablesCache& fluxVarsCache,
@@ -242,7 +240,7 @@ public:
{
const auto& data = fluxVarsCache.interiorBoundaryDataSelf();
if (data.faceType() == MpfaFaceTypes::interiorDirichlet)
- return getRho(data.facetVolVars(fvGeometry));
+ return getRho(data.facetVolVars(element, fvGeometry));
}
// use arithmetic mean of the densities around the scvf
@@ -261,7 +259,8 @@ public:
//! scaled to get the effective diffusivity. For this we use the effective diffusivity with
//! a diffusion coefficient of 1.0 as input. Then we scale the transmissibilites during flux
//! calculation (above) with the harmonic average of the two factors
- static Scalar computeEffectivityFactor(const FVElementGeometry& fvGeometry,
+ static Scalar computeEffectivityFactor(const Element& element,
+ const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const FluxVariablesCache& fluxVarsCache,
@@ -280,7 +279,7 @@ public:
insideVolVars.saturation(phaseIdx),
/*Diffusion coefficient*/ 1.0);
- const auto facetVolVars = data.facetVolVars(fvGeometry);
+ const auto facetVolVars = data.facetVolVars(element, fvGeometry);
const auto outsideFactor = EffDiffModel::effectiveDiffusivity(facetVolVars.porosity(),
facetVolVars.saturation(phaseIdx),
/*Diffusion coefficient*/ 1.0);
@@ -316,7 +315,8 @@ public:
}
template<typename GetXFunction>
- static Scalar computeInteriorBoundaryContribution(const FVElementGeometry& fvGeometry,
+ static Scalar computeInteriorBoundaryContribution(const Element& element,
+ const FVElementGeometry& fvGeometry,
const FluxVariablesCache& fluxVarsCache,
const GetXFunction& getX,
unsigned int phaseIdx, unsigned int compIdx)
@@ -332,7 +332,7 @@ public:
Scalar flux = 0.0;
for (auto&& data : fluxVarsCache.interiorBoundaryData())
if (data.faceType() == MpfaFaceTypes::interiorDirichlet)
- flux += tij[startIdx + data.localIndexInInteractionVolume()]*getX(data.facetVolVars(fvGeometry));
+ flux += tij[startIdx + data.localIndexInInteractionVolume()]*getX(data.facetVolVars(element, fvGeometry));
return flux;
}
diff --git a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
index bdabcd2a90..a85be8270a 100644
--- a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
@@ -147,17 +147,13 @@ public:
const auto& tij = fluxVarsCache.heatConductionTij();
const bool isInteriorBoundary = enableInteriorBoundaries && fluxVarsCache.isInteriorBoundary();
- // For interior Neumann boundaries when using Tpfa for Neumann boundary conditions, we simply
- // return the user-specified flux
+ // For interior Neumann boundaries when using Tpfa on boundaries, return the user-specified flux
if (isInteriorBoundary
&& useTpfaBoundary
&& fluxVarsCache.interiorBoundaryDataSelf().faceType() == MpfaFaceTypes::interiorNeumann)
return scvf.area()*
elemVolVars[scvf.insideScvIdx()].extrusionFactor()*
- problem.neumann(element,
- fvGeometry,
- elemVolVars,
- scvf)[energyEqIdx];
+ problem.neumann(element, fvGeometry, elemVolVars, scvf)[energyEqIdx];
// calculate Tij*tj
Scalar flux(0.0);
@@ -170,13 +166,14 @@ public:
return useTpfaBoundary ? flux : flux + fluxVarsCache.heatNeumannFlux();
// Handle interior boundaries
- flux += Implementation::computeInteriorBoundaryContribution(fvGeometry, fluxVarsCache);
+ flux += Implementation::computeInteriorBoundaryContribution(element, fvGeometry, fluxVarsCache);
// return overall resulting flux
return useTpfaBoundary ? flux : flux + fluxVarsCache.heatNeumannFlux();
}
- static Scalar computeInteriorBoundaryContribution(const FVElementGeometry& fvGeometry,
+ static Scalar computeInteriorBoundaryContribution(const Element& element,
+ const FVElementGeometry& fvGeometry,
const FluxVariablesCache& fluxVarsCache)
{
// obtain the transmissibilites associated with all pressures
@@ -190,7 +187,7 @@ public:
Scalar flux = 0.0;
for (auto&& data : fluxVarsCache.interiorBoundaryData())
if (data.faceType() == MpfaFaceTypes::interiorDirichlet)
- flux += tij[startIdx + data.localIndexInInteractionVolume()]*data.facetVolVars(fvGeometry).temperature();
+ flux += tij[startIdx + data.localIndexInInteractionVolume()]*data.facetVolVars(element, fvGeometry).temperature();
return flux;
}
diff --git a/dumux/discretization/cellcentered/mpfa/interiorboundarydata.hh b/dumux/discretization/cellcentered/mpfa/interiorboundarydata.hh
index 2e4d8067e6..17936b40a6 100644
--- a/dumux/discretization/cellcentered/mpfa/interiorboundarydata.hh
+++ b/dumux/discretization/cellcentered/mpfa/interiorboundarydata.hh
@@ -79,13 +79,12 @@ public:
{ return faceType_; }
//! returns the volume variables for interior dirichlet boundaries
- VolumeVariables facetVolVars(const FVElementGeometry& fvGeometry) const
+ VolumeVariables facetVolVars(const Element& element, const FVElementGeometry& fvGeometry) const
{
//! This can only be called for interior Dirichlet boundaries
assert(faceType_ == MpfaFaceTypes::interiorDirichlet && "requesting Dirichlet vol vars for a face which is"
"not marked as interior Dirichlet face.");
- auto element = problem_().model().globalFvGeometry().element(elementIndex());
auto priVars = problem_().dirichlet(element, fvGeometry.scvf(scvfIndex()));
VolumeVariables volVars;
diff --git a/dumux/discretization/upwindscheme.hh b/dumux/discretization/upwindscheme.hh
index 94d1cdc60d..a3b5b6c8bc 100644
--- a/dumux/discretization/upwindscheme.hh
+++ b/dumux/discretization/upwindscheme.hh
@@ -292,7 +292,7 @@ public:
if (facetCoupling || scvfFluxVarsCache.interiorBoundaryDataSelf().faceType() == MpfaFaceTypes::interiorDirichlet)
{
const auto& insideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().insideScvIdx()];
- const auto& outsideVolVars = scvfFluxVarsCache.interiorBoundaryDataSelf().facetVolVars(fluxVars.fvGeometry());
+ const auto& outsideVolVars = scvfFluxVarsCache.interiorBoundaryDataSelf().facetVolVars(fluxVars.element(), fluxVars.fvGeometry());
if (std::signbit(flux))
return flux*(upwindWeight*upwindTerm(outsideVolVars)
+ (1.0 - upwindWeight)*upwindTerm(insideVolVars));
@@ -361,7 +361,7 @@ public:
if (facetCoupling || scvfFluxVarsCache.interiorBoundaryDataSelf().faceType() == MpfaFaceTypes::interiorDirichlet)
{
const auto& insideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().insideScvIdx()];
- const auto& outsideVolVars = scvfFluxVarsCache.interiorBoundaryDataSelf().facetVolVars(fluxVars.fvGeometry());
+ const auto& outsideVolVars = scvfFluxVarsCache.interiorBoundaryDataSelf().facetVolVars(fluxVars.element(), fluxVars.fvGeometry());
if (std::signbit(flux))
return flux*(upwindWeight*upwindTerm(outsideVolVars)
+ (1.0 - upwindWeight)*upwindTerm(insideVolVars));
--
GitLab