diff --git a/dumux/discretization/staggered/elementvolumevariables.hh b/dumux/discretization/staggered/elementvolumevariables.hh
index cdf0ff46bbfcf42aca36c8e9c80a77375b567fbb..9e2cc021ceed2da95fcc04fa9cc5385fa7c8c1ea 100644
--- a/dumux/discretization/staggered/elementvolumevariables.hh
+++ b/dumux/discretization/staggered/elementvolumevariables.hh
@@ -93,8 +93,11 @@ class StaggeredElementVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/false
{
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using ElementSolution = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using GlobalVolumeVariables = typename GET_PROP_TYPE(TypeTag, GlobalVolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
@@ -103,6 +106,12 @@ class StaggeredElementVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/false
static const int dim = GridView::dimension;
using Element = typename GridView::template Codim<0>::Entity;
+ using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
+ typename DofTypeIndices::CellCenterIdx cellCenterIdx;
+ typename DofTypeIndices::FaceIdx faceIdx;
+
+ static constexpr auto numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter);
+
public:
//! Constructor
@@ -115,30 +124,35 @@ public:
const FVElementGeometry& fvGeometry,
const SolutionVector& sol)
{
- auto eIdx = globalVolVars().problem_().elementMapper().index(element);
+ clear();
- // stencil information
- const auto& neighborStencil = globalVolVars().problem_().model().stencils(element).neighborStencil();
- const auto numDofs = neighborStencil.size() + 1;
+ const auto& problem = globalVolVars().problem();
+ const auto& fvGridGeometry = fvGeometry.fvGridGeometry();
+ const auto globalI = fvGridGeometry.elementMapper().index(element);
+ const auto map = fvGridGeometry.connectivityMap();
+ const auto& connectivityMapI = map(cellCenterIdx, cellCenterIdx, globalI);
+ const auto numDofs = connectivityMapI.size();
+
+ auto&& scvI = fvGeometry.scv(globalI);
// resize local containers to the required size (for internal elements)
volumeVariables_.resize(numDofs);
volVarIndices_.resize(numDofs);
int localIdx = 0;
- // update the volume variables of the element at hand
- auto&& scvI = fvGeometry.scv(eIdx);
- volumeVariables_[localIdx].update(sol[eIdx], globalVolVars().problem_(), element, scvI);
- volVarIndices_[localIdx] = scvI.index();
- ++localIdx;
-
- // Update the volume variables of the neighboring elements
- for (auto globalJ : neighborStencil)
+ // Update the volume variables of the element at hand and the neighboring elements
+ for (auto globalJ : connectivityMapI)
{
- const auto& elementJ = fvGeometry.fvGridGeometry().element(globalJ);
+ const auto& elementJ = fvGridGeometry.element(globalJ);
auto&& scvJ = fvGeometry.scv(globalJ);
- volumeVariables_[localIdx].update(sol[globalJ], globalVolVars().problem_(), elementJ, scvJ);
- volVarIndices_[localIdx] = scvJ.index();
+ PrimaryVariables priVars(0.0);
+ priVars[cellCenterIdx] = sol[cellCenterIdx][globalJ];
+ ElementSolution elemSol{std::move(priVars)};
+ volumeVariables_[localIdx].update(elemSol,
+ problem,
+ elementJ,
+ scvJ);
+ volVarIndices_[localIdx] = scvJ.dofIndex();
++localIdx;
}
@@ -149,18 +163,33 @@ public:
if (!scvf.boundary())
continue;
- // check if boundary is a pure dirichlet boundary
- const auto bcTypes = globalVolVars().problem_().boundaryTypes(element, scvf);
- if (bcTypes.hasOnlyDirichlet())
- {
- const auto dirichletPriVars = globalVolVars().problem_().dirichlet(element, scvf);
+ const auto bcTypes = problem.boundaryTypes(element, scvf);
+
+ PrimaryVariables boundaryPriVars(0.0);
- volumeVariables_.resize(localIdx+1);
- volVarIndices_.resize(localIdx+1);
- volumeVariables_[localIdx].update(dirichletPriVars, globalVolVars().problem_(), element, scvI);
- volVarIndices_[localIdx] = scvf.outsideScvIdx();
- ++localIdx;
+ for(int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
+ {
+ if(bcTypes.isDirichlet(eqIdx) || bcTypes.isDirichletCell(eqIdx))
+ boundaryPriVars[cellCenterIdx][eqIdx] = problem.dirichlet(element, scvf)[cellCenterIdx][eqIdx];
+ else if(bcTypes.isNeumann(eqIdx) || bcTypes.isOutflow(eqIdx) || bcTypes.isSymmetry())
+ boundaryPriVars[cellCenterIdx][eqIdx] = sol[cellCenterIdx][scvf.insideScvIdx()][eqIdx];
+ //TODO: this assumes a zero-gradient for e.g. the pressure on the boundary
+ // could be made more general by allowing a non-zero-gradient, provided in problem file
+ else
+ if(eqIdx == Indices::pressureIdx)
+ DUNE_THROW(Dune::InvalidStateException, "Face at: " << scvf.center() << " has neither Dirichlet nor Neumann BC.");
}
+
+ volumeVariables_.resize(localIdx+1);
+ volVarIndices_.resize(localIdx+1);
+
+ ElementSolution elemSol{std::move(boundaryPriVars)};
+ volumeVariables_[localIdx].update(elemSol,
+ problem,
+ element,
+ scvI);
+ volVarIndices_[localIdx] = scvf.outsideScvIdx();
+ ++localIdx;
}
}
@@ -170,13 +199,21 @@ public:
const FVElementGeometry& fvGeometry,
const SolutionVector& sol)
{
- auto eIdx = globalVolVars().problem_().elementMapper().index(element);
+ clear();
+
+ const auto eIdx = fvGeometry.fvGridGeometry().elementMapper().index(element);
volumeVariables_.resize(1);
volVarIndices_.resize(1);
// update the volume variables of the element
auto&& scv = fvGeometry.scv(eIdx);
- volumeVariables_[0].update(sol[eIdx], globalVolVars().problem_(), element, scv);
+ PrimaryVariables priVars(0.0);
+ priVars[cellCenterIdx] = sol[cellCenterIdx][eIdx];
+ ElementSolution elemSol{std::move(priVars)};
+ volumeVariables_[0].update(elemSol,
+ globalVolVars().problem(),
+ element,
+ scv);
volVarIndices_[0] = scv.dofIndex();
}
@@ -196,6 +233,13 @@ public:
const GlobalVolumeVariables& globalVolVars() const
{ return *globalVolVarsPtr_; }
+ //! Clear all local storage
+ void clear()
+ {
+ volVarIndices_.clear();
+ volumeVariables_.clear();
+ }
+
private:
const GlobalVolumeVariables* globalVolVarsPtr_;
diff --git a/dumux/discretization/staggered/globalvolumevariables.hh b/dumux/discretization/staggered/globalvolumevariables.hh
index b9b0b3d89c226e014f653b06cf58a14643998c93..3514b8c35e38dbb0566edc8c1a6dc38fbf56358b 100644
--- a/dumux/discretization/staggered/globalvolumevariables.hh
+++ b/dumux/discretization/staggered/globalvolumevariables.hh
@@ -23,9 +23,6 @@
#ifndef DUMUX_DISCRETIZATION_STAGGERED_GLOBAL_VOLUMEVARIABLES_HH
#define DUMUX_DISCRETIZATION_STAGGERED_GLOBAL_VOLUMEVARIABLES_HH
-#include <dumux/common/basicproperties.hh>
-#include <dumux/discretization/staggered/elementvolumevariables.hh>
-
namespace Dumux
{
@@ -41,9 +38,6 @@ class StaggeredGlobalVolumeVariables
template<class TypeTag>
class StaggeredGlobalVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/true>
{
- // The local class needs to access and change volVars
- friend StaggeredElementVolumeVariables<TypeTag, true>;
-
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
@@ -154,15 +148,15 @@ private:
template<class TypeTag>
class StaggeredGlobalVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/false>
{
- // local class needs access to the problem
- friend StaggeredElementVolumeVariables<TypeTag, false>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
public:
- void update(Problem& problem, const SolutionVector& sol)
- { problemPtr_ = &problem; }
+ StaggeredGlobalVolumeVariables(const Problem& problem) : problemPtr_(&problem) {}
+
+ void update(const FVGridGeometry& fvGridGeometry, const SolutionVector& sol) {}
/*!
* \brief Return a local restriction of this global object
@@ -172,11 +166,12 @@ public:
friend inline ElementVolumeVariables localView(const StaggeredGlobalVolumeVariables& global)
{ return ElementVolumeVariables(global); }
-private:
- Problem& problem_() const
+ const Problem& problem() const
{ return *problemPtr_;}
- Problem* problemPtr_;
+private:
+
+ const Problem* problemPtr_;
};
} // end namespace