diff --git a/dumux/discretization/staggered/stencils.hh b/dumux/discretization/staggered/stencils.hh
index f8a0f42a3bb8c381e2d4efba28684b3294e31b6b..0cc0b8aaf7d20d0c12468c5c57cf4c9caa068fc8 100644
--- a/dumux/discretization/staggered/stencils.hh
+++ b/dumux/discretization/staggered/stencils.hh
@@ -50,43 +50,50 @@ public:
const Element& element,
const FVElementGeometry& fvGeometry)
{
- elementStencil_.clear();
+
+ const auto globalI = problem.elementMapper().index(element);
+
+ cellCenterToCellCenterStencil_.clear();
+ cellCenterToFaceStencil_.clear();
+
+ // the cell center dof indices
+ cellCenterToCellCenterStencil_.push_back(globalI);
// loop over sub control faces
for (auto&& scvf : scvfs(fvGeometry))
{
- FluxVariables fluxVars;
- const auto& centerStencil = fluxVars.computeCellCenterStencil(problem, element, fvGeometry, scvf);
- elementStencil_.insert(elementStencil_.end(), centerStencil.begin(), centerStencil.end());
+ if (!scvf.boundary())
+ cellCenterToCellCenterStencil_.push_back(scvf.outsideScvIdx());
+
+ cellCenterToFaceStencil_.push_back(scvf.dofIndexSelf());
}
- // make values in elementstencil unique
- std::sort(elementStencil_.begin(), elementStencil_.end());
- elementStencil_.erase(std::unique(elementStencil_.begin(), elementStencil_.end()), elementStencil_.end());
- auto globalI = problem.elementMapper().index(element);
- neighborStencil_ = elementStencil_;
+ // make values unique
+ std::sort(cellCenterToCellCenterStencil_.begin(), cellCenterToCellCenterStencil_.end());
+ cellCenterToCellCenterStencil_.erase(std::unique(cellCenterToCellCenterStencil_.begin(), cellCenterToCellCenterStencil_.end()), cellCenterToCellCenterStencil_.end());
- // remove the element itself and possible ghost neighbors from the neighbor stencil
- neighborStencil_.erase(std::remove_if(neighborStencil_.begin(), neighborStencil_.end(),
- [globalI](int i){ return (i == globalI); }),
- neighborStencil_.end());
+ std::sort(cellCenterToFaceStencil_.begin(), cellCenterToFaceStencil_.end());
+ cellCenterToFaceStencil_.erase(std::unique(cellCenterToFaceStencil_.begin(), cellCenterToFaceStencil_.end()), cellCenterToFaceStencil_.end());
}
+
//! The full element stencil (all element this element is interacting with)
- const Stencil& elementStencil() const
+ const Stencil& cellCenterToCellCenterStencil() const
{
- return elementStencil_;
+ return cellCenterToCellCenterStencil_;
}
- //! The full element stencil without this element
- const Stencil& neighborStencil() const
+ //! The full element stencil (all element this element is interacting with)
+ const Stencil& cellCenterToFaceStencil() const
{
- return neighborStencil_;
+ return cellCenterToFaceStencil_;
}
+
+
private:
- Stencil elementStencil_;
- Stencil neighborStencil_;
+ Stencil cellCenterToCellCenterStencil_;
+ Stencil cellCenterToFaceStencil_;
};
template<class TypeTag>
@@ -106,23 +113,45 @@ public:
void update(const Problem& problem,
const SubControlVolumeFace& scvf)
{
- faceStencil_.clear();
- FluxVariables fluxVars;
+ faceToCellCenterStencil_.clear();
+ faceToFaceStencil_.clear();
+
+ faceToCellCenterStencil_.push_back(scvf.insideScvIdx());
- const auto& faceStencil = fluxVars.computeFaceStencil(problem, scvf);
- faceStencil_.insert(faceStencil_.end(), faceStencil.begin(), faceStencil.end());
- std::sort(faceStencil_.begin(), faceStencil_.end());
- faceStencil_.erase(std::unique(faceStencil_.begin(), faceStencil_.end()), faceStencil_.end());
+ faceToFaceStencil_.push_back(scvf.dofIndexSelf());
+ faceToFaceStencil_.push_back(scvf.dofIndexOpposite());
+
+ for(const auto& data : scvf.pairData())
+ {
+ const auto& outerParallelElementDofIdx = data.outerParallelElementDofIdx;
+ const auto& outerParallelFaceDofIdx = data.outerParallelFaceDofIdx;
+ if(outerParallelElementDofIdx >= 0)
+ faceToCellCenterStencil_.push_back(outerParallelElementDofIdx);
+ if(outerParallelFaceDofIdx >= 0)
+ faceToFaceStencil_.push_back(outerParallelFaceDofIdx);
+
+ faceToFaceStencil_.push_back(data.normalPair.first);
+ if(!scvf.boundary())
+ faceToFaceStencil_.push_back(data.normalPair.second);
+ }
+
+ }
+
+ //! The full face stencil (all dofs this face is interacting with)
+ const Stencil& faceToCellCenterStencil() const
+ {
+ return faceToCellCenterStencil_;
}
//! The full face stencil (all dofs this face is interacting with)
- const Stencil& faceStencil() const
+ const Stencil& faceToFaceStencil() const
{
- return faceStencil_;
+ return faceToFaceStencil_;
}
private:
- Stencil faceStencil_;
+ Stencil faceToCellCenterStencil_;
+ Stencil faceToFaceStencil_;
};
@@ -150,11 +179,13 @@ public:
const IndexType numElements = problem.gridView().size(0);
const IndexType numFacets = problem.gridView().size(1);
const IndexType numBoundaryFacets = problem.gridView().grid().numBoundarySegments();
- elementStencils_.resize(numElements);
+ cellCenterStencils_.resize(numElements);
faceStencils_.resize(2*numFacets - numBoundaryFacets);
- std::vector<Stencil> fullFaceStencils;
- fullFaceStencils.resize(numFacets);
+ std::vector<Stencil> fullFaceToCellCenterStencils;
+ fullFaceToCellCenterStencils.resize(numFacets);
+ std::vector<Stencil> fullfaceToFaceStencils;
+ fullfaceToFaceStencils.resize(numFacets);
for (const auto& element : elements(problem.gridView()))
{
@@ -162,7 +193,7 @@ public:
// restrict the FvGeometry locally and bind to the element
auto fvGeometry = localView(problem.model().globalFvGeometry());
fvGeometry.bindElement(element); // for TPFA bind element is enough
- elementStencils_[eIdx].update(problem, element, fvGeometry);
+ cellCenterStencils_[eIdx].update(problem, element, fvGeometry);
// loop over sub control faces
for (auto&& scvf : scvfs(fvGeometry))
@@ -170,21 +201,28 @@ public:
faceStencils_[scvf.index()].update(problem, scvf);
const IndexType idx = scvf.dofIndexSelf() - numElements;
- const auto& faceStencil = faceStencils_[scvf.index()].faceStencil();
- fullFaceStencils[idx].insert(fullFaceStencils[idx].end(), faceStencil.begin(), faceStencil.end());
- std::sort(fullFaceStencils[idx].begin(), fullFaceStencils[idx].end());
- fullFaceStencils[idx].erase(std::unique(fullFaceStencils[idx].begin(), fullFaceStencils[idx].end()), fullFaceStencils[idx].end());
+
+ const auto& faceToCellCenterStencil = faceStencils_[scvf.index()].faceToCellCenterStencil();
+ fullFaceToCellCenterStencils[idx].insert(fullFaceToCellCenterStencils[idx].end(), faceToCellCenterStencil.begin(), faceToCellCenterStencil.end());
+ std::sort(fullFaceToCellCenterStencils[idx].begin(), fullFaceToCellCenterStencils[idx].end());
+ fullFaceToCellCenterStencils[idx].erase(std::unique(fullFaceToCellCenterStencils[idx].begin(), fullFaceToCellCenterStencils[idx].end()), fullFaceToCellCenterStencils[idx].end());
+
+ const auto& faceToFaceStencil = faceStencils_[scvf.index()].faceToFaceStencil();
+ fullfaceToFaceStencils[idx].insert(fullfaceToFaceStencils[idx].end(), faceToFaceStencil.begin(), faceToFaceStencil.end());
+ std::sort(fullfaceToFaceStencils[idx].begin(), fullfaceToFaceStencils[idx].end());
+ fullfaceToFaceStencils[idx].erase(std::unique(fullfaceToFaceStencils[idx].begin(), fullfaceToFaceStencils[idx].end()), fullfaceToFaceStencils[idx].end());
}
}
// TODO: is this a good idea?
- fullFaceDofStencils_ = std::make_unique<decltype(fullFaceStencils)>(fullFaceStencils);
+ fullFaceToCellCenterStencils_ = std::make_unique<decltype(fullFaceToCellCenterStencils)>(fullFaceToCellCenterStencils);
+ fullfaceToFaceStencils_ = std::make_unique<decltype(fullfaceToFaceStencils)>(fullfaceToFaceStencils);
}
//! overload for elements
auto& get(const Element& entity) const
{
- return elementStencils_[problemPtr_->elementMapper().index(entity)];
+ return cellCenterStencils_[problemPtr_->elementMapper().index(entity)];
}
//! overload for faces
@@ -197,28 +235,49 @@ public:
/*!
* \brief Returns the size of a complete face dof stencil
*/
- size_t completeFaceDofStencilSize(const int idx) const
+ size_t fullFaceToCellCenterStencilSize(const int idx) const
+ {
+// const IndexType numElements = problemPtr_->gridView().size(0);
+ assert(fullFaceToCellCenterStencils_ && "fullFaceToCellCenterStencils_ has already been called and deleted!");
+ return fullFaceToCellCenterStencils_.get()[0][idx/*-numElements*/].size();
+ // TODO: why does this not work?
+ }
+
+ /*!
+ * \brief Returns the size of a complete face dof stencil
+ */
+ size_t fullfaceToFaceStencilSize(const int idx) const
{
// const IndexType numElements = problemPtr_->gridView().size(0);
- assert(fullFaceDofStencils_ && "fullFaceDofStencils_ has already been called and deleted!");
- return fullFaceDofStencils_.get()[0][idx/*-numElements*/].size();
+ assert(fullfaceToFaceStencils_ && "fullfaceToFaceStencils_ has already been called and deleted!");
+ return fullfaceToFaceStencils_.get()[0][idx/*-numElements*/].size();
// TODO: why does this not work?
}
/*!
* \brief Returns a unique pointer to the complete face dof stencils which is used once for setting up the global matrix and deleted afterwards
*/
- auto getFullFaceDofStencilsPtr()
+ auto getFullFaceToCellCenterStencilsPtr()
+ {
+ assert(fullFaceToCellCenterStencils_ && "fullFaceToCellCenterStencils_ has already been used and deleted!");
+ return std::move(fullFaceToCellCenterStencils_);
+ }
+
+ /*!
+ * \brief Returns a unique pointer to the complete face dof stencils which is used once for setting up the global matrix and deleted afterwards
+ */
+ auto getFullfaceToFaceStencilsPtr()
{
- assert(fullFaceDofStencils_ && "fullFaceDofStencils_ has already been used and deleted!");
- return std::move(fullFaceDofStencils_);
+ assert(fullfaceToFaceStencils_ && "fullfaceToFaceStencils_ has already been used and deleted!");
+ return std::move(fullfaceToFaceStencils_);
}
private:
- std::vector<StaggeredCellCenterStencils<TypeTag>> elementStencils_;
+ std::vector<StaggeredCellCenterStencils<TypeTag>> cellCenterStencils_;
std::vector<StaggeredFaceStencils<TypeTag>> faceStencils_;
- std::unique_ptr<std::vector<Stencil>> fullFaceDofStencils_;
+ std::unique_ptr<std::vector<Stencil>> fullFaceToCellCenterStencils_;
+ std::unique_ptr<std::vector<Stencil>> fullfaceToFaceStencils_;
const Problem* problemPtr_;
diff --git a/dumux/freeflow/staggered/model.hh b/dumux/freeflow/staggered/model.hh
index 179c17bab9cf254bfe0ef4b10bfb08c6e0ffbcd3..08ef270ff2b222e8f21f18f2655fbdeac3503700 100644
--- a/dumux/freeflow/staggered/model.hh
+++ b/dumux/freeflow/staggered/model.hh
@@ -143,17 +143,33 @@ public:
/*!
* \brief Returns the size of a complete face dof stencil
*/
- size_t completeFaceDofStencilSize(const int idx) const
+ size_t fullFaceToCellCenterStencilSize(const int idx) const
{
- return this->stencilsVector_.completeFaceDofStencilSize(idx);
+ return this->stencilsVector_.fullFaceToCellCenterStencilSize(idx);
+ }
+
+ /*!
+ * \brief Returns the size of a complete face dof stencil
+ */
+ size_t fullfaceToFaceStencilSize(const int idx) const
+ {
+ return this->stencilsVector_.fullfaceToFaceStencilSize(idx);
+ }
+
+ /*!
+ * \brief Returns a unique pointer to the complete face dof stencils which is used once for setting up the global matrix and deleted afterwards
+ */
+ auto getFullFaceToCellCenterStencilsPtr()
+ {
+ return this->stencilsVector_.getFullFaceToCellCenterStencilsPtr();
}
/*!
* \brief Returns a unique pointer to the complete face dof stencils which is used once for setting up the global matrix and deleted afterwards
*/
- auto getFullFaceDofStencilsPtr()
+ auto getFullfaceToFaceStencilsPtr()
{
- return this->stencilsVector_.getFullFaceDofStencilsPtr();
+ return this->stencilsVector_.getFullfaceToFaceStencilsPtr();
}
};
diff --git a/dumux/implicit/staggered/assembler.hh b/dumux/implicit/staggered/assembler.hh
index af0a097e5f5ab0c7e7a344038dc7d290dba51154..3d9689cbcc76022a1038d15a9a0f5a8ad836086a 100644
--- a/dumux/implicit/staggered/assembler.hh
+++ b/dumux/implicit/staggered/assembler.hh
@@ -50,15 +50,18 @@ class StaggeredAssembler : public ImplicitAssembler<TypeTag>
{
// the global index of the element at hand
const auto globalI = this->elementMapper_().index(element);
- const auto& stencil = this->model_().stencils(element).elementStencil();
+ const auto& cellCenterToCellCenterStencil = this->model_().stencils(element).cellCenterToCellCenterStencil();
+ const auto& cellCenterToFaceStencil = this->model_().stencils(element).cellCenterToFaceStencil();
+ const auto size = cellCenterToCellCenterStencil.size() + cellCenterToFaceStencil.size();
- this->matrix().setrowsize(globalI, stencil.size());
+ this->matrix().setrowsize(globalI, size);
}
for(int facetIdx = 0; facetIdx < this->gridView_().size(1); ++facetIdx)
{
const int faceDofxIdx = facetIdx + this->gridView_().size(0);
- const int size = this->model_().completeFaceDofStencilSize(facetIdx);
+ auto size = this->model_().fullFaceToCellCenterStencilSize(facetIdx);
+ size += this->model_().fullfaceToFaceStencilSize(facetIdx);
this->matrix().setrowsize(faceDofxIdx, size);
}
this->matrix().endrowsizes();
@@ -70,22 +73,34 @@ class StaggeredAssembler : public ImplicitAssembler<TypeTag>
{
// the global index of the element at hand
const auto globalI = this->elementMapper_().index(element);
- const auto& stencil = this->model_().stencils(element).elementStencil();
+ const auto& cellCenterToCellCenterStencil = this->model_().stencils(element).cellCenterToCellCenterStencil();
+ const auto& cellCenterToFaceStencil = this->model_().stencils(element).cellCenterToFaceStencil();
- for (auto&& globalJ : stencil)
+ for (auto&& globalJ : cellCenterToCellCenterStencil)
+ this->matrix().addindex(globalI, globalJ);
+ for (auto&& globalJ : cellCenterToFaceStencil)
this->matrix().addindex(globalI, globalJ);
}
- auto ptr = this->model_().getFullFaceDofStencilsPtr();
+ auto ptr1 = this->model_().getFullFaceToCellCenterStencilsPtr();
+ auto ptr2 = this->model_().getFullfaceToFaceStencilsPtr();
- if(ptr)
+ if(ptr1 && ptr2)
std::cout << "success!!!" << std::endl;
- const auto& fullFaceDofStencils = ptr.get()[0];
+ const auto& fullFaceToCellCenterStencils = ptr1.get()[0];
+ const auto& fullfaceToFaceStencils = ptr2.get()[0];
int globalI = this->gridView_().size(0);
- for(const auto& stencil : fullFaceDofStencils)
+ for(const auto& stencil : fullFaceToCellCenterStencils)
+ {
+ for(auto&& globalJ : stencil)
+ this->matrix().addindex(globalI, globalJ);
+ ++globalI;
+ }
+ globalI = this->gridView_().size(0);
+ for(const auto& stencil : fullfaceToFaceStencils)
{
for(auto&& globalJ : stencil)
this->matrix().addindex(globalI, globalJ);
diff --git a/dumux/implicit/staggered/localjacobian.hh b/dumux/implicit/staggered/localjacobian.hh
index 719964d6996d46f36c3196872245a872267bb871..0059b955c2fa60c7694d319bf26c58cb806334fe 100644
--- a/dumux/implicit/staggered/localjacobian.hh
+++ b/dumux/implicit/staggered/localjacobian.hh
@@ -190,157 +190,157 @@ private:
const bool isGhost)
{
// get stencil informations
- const auto& neighborStencil = this->model_().stencils(element).neighborStencil();
- const auto numNeighbors = neighborStencil.size();
+// const auto& neighborStencil = this->model_().stencils(element).neighborStencil();
+// const auto numNeighbors = neighborStencil.size();
// the localresidual class used for the flux calculations
- LocalResidual localRes;
- localRes.init(this->problem_());
-
- // container to store the neighboring elements
- std::vector<Element> neighborElements;
- neighborElements.reserve(numNeighbors);
-
- // get the elements and calculate the flux into the element in the undeflected state
- Dune::BlockVector<PrimaryVariables> origFlux(numNeighbors);
- origFlux = 0.0;
-
- unsigned int j = 0;
- for (auto globalJ : neighborStencil)
- {
- auto elementJ = fvGeometry.globalFvGeometry().element(globalJ);
- neighborElements.push_back(elementJ);
-
- for (auto fluxVarIdx : assemblyMap_[globalI_][j])
- {
- auto&& scvf = fvGeometry.scvf(fluxVarIdx);
- origFlux[j] += localRes.evalFlux_(elementJ, fvGeometry, curElemVolVars, scvf, elemFluxVarsCache[scvf]);
- }
-
- ++j;
- }
-
- auto&& scv = fvGeometry.scv(globalI_);
- auto& curVolVars = getCurVolVars(curElemVolVars, scv);
- // save a copy of the original vol vars
- VolumeVariables origVolVars(curVolVars);
-
- // derivatives in the neighbors with repect to the current elements
- Dune::BlockVector<PrimaryVariables> neighborDeriv(numNeighbors);
- for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
- {
- // derivatives of element dof with respect to itself
- PrimaryVariables partialDeriv(0.0);
-
- if (isGhost)
- partialDeriv[pvIdx] = 1.0;
-
- neighborDeriv = 0.0;
- PrimaryVariables priVars(this->model_().curSol()[globalI_]);
-
- Scalar eps = this->numericEpsilon(scv, curVolVars, pvIdx);
- Scalar delta = 0;
-
- if (numericDifferenceMethod_ >= 0)
- {
- // we are not using backward differences, i.e. we need to
- // calculate f(x + \epsilon)
-
- // deflect primary variables
- priVars[pvIdx] += eps;
- delta += eps;
-
- // update the volume variables and bind the flux var cache again
- curVolVars.update(priVars, this->problem_(), element, scv);
- elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
-
- if (!isGhost)
- {
- // calculate the residual with the deflected primary variables
- this->localResidual().eval(element, fvGeometry, prevElemVolVars, curElemVolVars, elemBcTypes, elemFluxVarsCache);
-
- // store the residual and the storage term
- partialDeriv = this->localResidual().residual(0);
- }
-
- // calculate the fluxes into element with the deflected primary variables
- for (std::size_t k = 0; k < numNeighbors; ++k)
- {
- for (auto fluxVarIdx : assemblyMap_[globalI_][k])
- {
- auto&& scvf = fvGeometry.scvf(fluxVarIdx);
- neighborDeriv[k] += localRes.evalFlux_(neighborElements[k], fvGeometry, curElemVolVars, scvf, elemFluxVarsCache[scvf]);
- }
- }
- }
- else
- {
- // we are using backward differences, i.e. we don't need
- // to calculate f(x + \epsilon) and we can recycle the
- // (already calculated) residual f(x)
- if (!isGhost)
- partialDeriv = this->residual(0);
- neighborDeriv = origFlux;
- }
-
- if (numericDifferenceMethod_ <= 0)
- {
- // we are not using forward differences, i.e. we
- // need to calculate f(x - \epsilon)
-
- // deflect the primary variables
- priVars[pvIdx] -= delta + eps;
- delta += eps;
-
- // update the volume variables and bind the flux var cache again
- curVolVars.update(priVars, this->problem_(), element, scv);
- elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
-
- if (!isGhost)
- {
- // calculate the residual with the deflected primary variables
- this->localResidual().eval(element, fvGeometry, prevElemVolVars, curElemVolVars, elemBcTypes, elemFluxVarsCache);
-
- // subtract the residual from the derivative storage
- partialDeriv -= this->localResidual().residual(0);
- }
-
- // calculate the fluxes into element with the deflected primary variables
- for (std::size_t k = 0; k < numNeighbors; ++k)
- {
- for (auto fluxVarIdx : assemblyMap_[globalI_][k])
- {
- auto&& scvf = fvGeometry.scvf(fluxVarIdx);
- neighborDeriv[k] -= localRes.evalFlux_(neighborElements[k], fvGeometry, curElemVolVars, scvf, elemFluxVarsCache[scvf]);
- }
- }
- }
- else
- {
- // we are using forward differences, i.e. we don't need to
- // calculate f(x - \epsilon) and we can recycle the
- // (already calculated) residual f(x)
- if (!isGhost)
- partialDeriv -= this->residual(0);
- neighborDeriv -= origFlux;
- }
-
- // divide difference in residuals by the magnitude of the
- // deflections between the two function evaluation
- if (!isGhost)
- partialDeriv /= delta;
- neighborDeriv /= delta;
-
- // restore the original state of the scv's volume variables
- curVolVars = origVolVars;
-
- // update the global jacobian matrix with the current partial derivatives
- this->updateGlobalJacobian_(matrix, globalI_, globalI_, pvIdx, partialDeriv);
-
- j = 0;
- for (auto globalJ : neighborStencil)
- this->updateGlobalJacobian_(matrix, globalJ, globalI_, pvIdx, neighborDeriv[j++]);
- }
+// LocalResidual localRes;
+// localRes.init(this->problem_());
+//
+// // container to store the neighboring elements
+// std::vector<Element> neighborElements;
+// // neighborElements.reserve(numNeighbors);
+//
+// // get the elements and calculate the flux into the element in the undeflected state
+// Dune::BlockVector<PrimaryVariables> origFlux(numNeighbors);
+// origFlux = 0.0;
+//
+// unsigned int j = 0;
+// for (auto globalJ : neighborStencil)
+// {
+// auto elementJ = fvGeometry.globalFvGeometry().element(globalJ);
+// neighborElements.push_back(elementJ);
+//
+// for (auto fluxVarIdx : assemblyMap_[globalI_][j])
+// {
+// auto&& scvf = fvGeometry.scvf(fluxVarIdx);
+// origFlux[j] += localRes.evalFlux_(elementJ, fvGeometry, curElemVolVars, scvf, elemFluxVarsCache[scvf]);
+// }
+//
+// ++j;
+// }
+//
+// auto&& scv = fvGeometry.scv(globalI_);
+// auto& curVolVars = getCurVolVars(curElemVolVars, scv);
+// // save a copy of the original vol vars
+// VolumeVariables origVolVars(curVolVars);
+//
+// // derivatives in the neighbors with repect to the current elements
+// Dune::BlockVector<PrimaryVariables> neighborDeriv(numNeighbors);
+// for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
+// {
+// // derivatives of element dof with respect to itself
+// PrimaryVariables partialDeriv(0.0);
+//
+// if (isGhost)
+// partialDeriv[pvIdx] = 1.0;
+//
+// neighborDeriv = 0.0;
+// PrimaryVariables priVars(this->model_().curSol()[globalI_]);
+//
+// Scalar eps = this->numericEpsilon(scv, curVolVars, pvIdx);
+// Scalar delta = 0;
+//
+// if (numericDifferenceMethod_ >= 0)
+// {
+// // we are not using backward differences, i.e. we need to
+// // calculate f(x + \epsilon)
+//
+// // deflect primary variables
+// priVars[pvIdx] += eps;
+// delta += eps;
+//
+// // update the volume variables and bind the flux var cache again
+// curVolVars.update(priVars, this->problem_(), element, scv);
+// elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
+//
+// if (!isGhost)
+// {
+// // calculate the residual with the deflected primary variables
+// this->localResidual().eval(element, fvGeometry, prevElemVolVars, curElemVolVars, elemBcTypes, elemFluxVarsCache);
+//
+// // store the residual and the storage term
+// partialDeriv = this->localResidual().residual(0);
+// }
+//
+// // calculate the fluxes into element with the deflected primary variables
+// for (std::size_t k = 0; k < numNeighbors; ++k)
+// {
+// for (auto fluxVarIdx : assemblyMap_[globalI_][k])
+// {
+// auto&& scvf = fvGeometry.scvf(fluxVarIdx);
+// neighborDeriv[k] += localRes.evalFlux_(neighborElements[k], fvGeometry, curElemVolVars, scvf, elemFluxVarsCache[scvf]);
+// }
+// }
+// }
+// else
+// {
+// // we are using backward differences, i.e. we don't need
+// // to calculate f(x + \epsilon) and we can recycle the
+// // (already calculated) residual f(x)
+// if (!isGhost)
+// partialDeriv = this->residual(0);
+// neighborDeriv = origFlux;
+// }
+//
+// if (numericDifferenceMethod_ <= 0)
+// {
+// // we are not using forward differences, i.e. we
+// // need to calculate f(x - \epsilon)
+//
+// // deflect the primary variables
+// priVars[pvIdx] -= delta + eps;
+// delta += eps;
+//
+// // update the volume variables and bind the flux var cache again
+// curVolVars.update(priVars, this->problem_(), element, scv);
+// elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
+//
+// if (!isGhost)
+// {
+// // calculate the residual with the deflected primary variables
+// this->localResidual().eval(element, fvGeometry, prevElemVolVars, curElemVolVars, elemBcTypes, elemFluxVarsCache);
+//
+// // subtract the residual from the derivative storage
+// partialDeriv -= this->localResidual().residual(0);
+// }
+//
+// // calculate the fluxes into element with the deflected primary variables
+// for (std::size_t k = 0; k < numNeighbors; ++k)
+// {
+// for (auto fluxVarIdx : assemblyMap_[globalI_][k])
+// {
+// auto&& scvf = fvGeometry.scvf(fluxVarIdx);
+// neighborDeriv[k] -= localRes.evalFlux_(neighborElements[k], fvGeometry, curElemVolVars, scvf, elemFluxVarsCache[scvf]);
+// }
+// }
+// }
+// else
+// {
+// // we are using forward differences, i.e. we don't need to
+// // calculate f(x - \epsilon) and we can recycle the
+// // (already calculated) residual f(x)
+// if (!isGhost)
+// partialDeriv -= this->residual(0);
+// neighborDeriv -= origFlux;
+// }
+//
+// // divide difference in residuals by the magnitude of the
+// // deflections between the two function evaluation
+// if (!isGhost)
+// partialDeriv /= delta;
+// neighborDeriv /= delta;
+//
+// // restore the original state of the scv's volume variables
+// curVolVars = origVolVars;
+//
+// // update the global jacobian matrix with the current partial derivatives
+// this->updateGlobalJacobian_(matrix, globalI_, globalI_, pvIdx, partialDeriv);
+//
+// j = 0;
+// for (auto globalJ : neighborStencil)
+// this->updateGlobalJacobian_(matrix, globalJ, globalI_, pvIdx, neighborDeriv[j++]);
+// }
}
//! If the global vol vars caching is enabled we have to modify the global volvar object