From c921bc8d78f1c653b2bf002b37b9803d8041c313 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Thu, 8 Dec 2016 13:18:33 +0100
Subject: [PATCH] [staggeredGrid][globalVolVars] Re-introduce boundary face
volVars
* Use dirichtlet values from problem or assume zero gradient and
use values from cell at boundary
---
.../staggered/globalvolumevariables.hh | 51 +++++++++++--------
1 file changed, 31 insertions(+), 20 deletions(-)
diff --git a/dumux/discretization/staggered/globalvolumevariables.hh b/dumux/discretization/staggered/globalvolumevariables.hh
index a41f5ce88b..8919c36dc0 100644
--- a/dumux/discretization/staggered/globalvolumevariables.hh
+++ b/dumux/discretization/staggered/globalvolumevariables.hh
@@ -63,6 +63,9 @@ class StaggeredGlobalVolumeVariables<TypeTag, /*enableGlobalVolVarsCache*/true>
static const int dim = GridView::dimension;
using Element = typename GridView::template Codim<0>::Entity;
+ using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+
+ enum { numEqCellCenter = GET_PROP_VALUE(TypeTag, NumEqCellCenter) };
public:
void update(Problem& problem, const SolutionVector& sol)
@@ -70,9 +73,9 @@ public:
problemPtr_ = &problem;
auto numScv = problem.model().globalFvGeometry().numScv();
-// auto numBoundaryScvf = problem.model().globalFvGeometry().numBoundaryScvf();
+ auto numBoundaryScvf = problem.model().globalFvGeometry().numBoundaryScvf();
- volumeVariables_.resize(numScv /*+ numBoundaryScvf*/);
+ volumeVariables_.resize(numScv + numBoundaryScvf);
for (const auto& element : elements(problem.gridView()))
{
auto fvGeometry = localView(problem.model().globalFvGeometry());
@@ -81,24 +84,32 @@ public:
for (auto&& scv : scvs(fvGeometry))
volumeVariables_[scv.index()].update(sol[cellCenterIdx][scv.dofIndex()], problem, element, scv);
-// // handle the boundary volume variables
-// for (auto&& scvf : scvfs(fvGeometry))
-// {
-// // if we are not on a boundary, skip the rest
-// if (!scvf.boundary())
-// continue;
-//
-// // check if boundary is a pure dirichlet boundary
-// const auto bcTypes = problem.boundaryTypes(element, scvf);
-// if (bcTypes.hasOnlyDirichlet())
-// {
-// const auto insideScvIdx = scvf.insideScvIdx();
-// const auto& insideScv = fvGeometry.scv(insideScvIdx);
-// const auto dirichletPriVars = problem.ccDirichlet(element, scvf);
-//
-// volumeVariables_[scvf.outsideScvIdx()].update(dirichletPriVars, problem, element, insideScv);
-// }
-// }
+ // handle the boundary volume variables
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ // if we are not on a boundary, skip the rest
+ if (!scvf.boundary())
+ continue;
+
+ const auto bcTypes = problem.boundaryTypes(element, scvf);
+ const auto insideScvIdx = scvf.insideScvIdx();
+ const auto& insideScv = fvGeometry.scv(insideScvIdx);
+
+ auto boundaryPriVars = CellCenterPrimaryVariables(0.0);
+
+ for(int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx)
+ {
+ if(bcTypes.isDirichlet(eqIdx))
+ boundaryPriVars[eqIdx] = problem.ccDirichlet(element, scvf)[eqIdx];
+ else if(bcTypes.isNeumann(eqIdx))
+ boundaryPriVars[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
+ DUNE_THROW(Dune::InvalidStateException, "Face at: " << scvf.center() << " has neither Dirichlet nor Neumann BC.");
+ }
+ volumeVariables_[scvf.outsideScvIdx()].update(boundaryPriVars, problem, element, insideScv);
+ }
}
}
--
GitLab