diff --git a/dumux/implicit/assembler.hh b/dumux/implicit/assembler.hh
index 2193aaaf9634a69ab04af3c58169763a74886878..dd5ffdd40a789e49969d7bdaf7cf3c8d21e95846 100644
--- a/dumux/implicit/assembler.hh
+++ b/dumux/implicit/assembler.hh
@@ -152,31 +152,9 @@ protected:
     {
         asImp_().resetSystem_();
 
-        // reassemble the elements...
-        for (const auto& element : elements(gridView_())) {
-            if (element.partitionType() == Dune::GhostEntity)
-            {
-                asImp_().assembleGhostElement_(element);
-            }
-            else
-            {
-                asImp_().assembleElement_(element);
-            }
-        }
-    }
-    // "assemble" a ghost element
-    void assembleGhostElement_(const Element &element)
-    {
-        int globalI = elementMapper_().index(element);
-
-        // update the right hand side
-        residual_[globalI] = 0.0;
-
-        // update the diagonal entry
-        typedef typename JacobianMatrix::block_type BlockType;
-        BlockType &J = (*matrix_)[globalI][globalI];
-        for (int j = 0; j < BlockType::rows; ++j)
-            J[j][j] = 1.0;
+        // assemble the elements...
+        for (const auto& element : elements(gridView_()))
+            asImp_().assembleElement_(element);
     }
 
 protected:
diff --git a/dumux/implicit/cellcentered/elementboundarytypes.hh b/dumux/implicit/cellcentered/elementboundarytypes.hh
index de787f440a86e9afd2755cb54c807ea453c35ca6..db59ac6170b8119b7756bc273bae5a51beacf7cd 100644
--- a/dumux/implicit/cellcentered/elementboundarytypes.hh
+++ b/dumux/implicit/cellcentered/elementboundarytypes.hh
@@ -87,12 +87,12 @@ public:
 
         (*this)[0].reset();
 
-        for (auto&& scv : fvGeometry.scvs())
+        for (const auto& scv : fvGeometry.scvs())
         {
             if (!problem.model().onBoundary(scv))
                 return;
 
-            for (auto&& scvFace : fvGeometry.scvfs())
+            for (const auto& scvFace : fvGeometry.scvfs())
             {
                 if (!scvFace.boundary())
                     continue;
diff --git a/dumux/implicit/cellcentered/localjacobian.hh b/dumux/implicit/cellcentered/localjacobian.hh
index efb921092294e6a43381904aba68b1f73d0c49e3..0cf4858490f5cdb28f612a61af0071f7f63630f3 100644
--- a/dumux/implicit/cellcentered/localjacobian.hh
+++ b/dumux/implicit/cellcentered/localjacobian.hh
@@ -80,6 +80,7 @@ class CCLocalJacobian : public ImplicitLocalJacobian<TypeTag>
     typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
     typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
     typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
     typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
     typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
 
@@ -114,22 +115,31 @@ public:
         assemblyMap_.resize(problem.gridView().size(0));
         for (const auto& element : elements(problem.gridView()))
         {
+            // bind FVGeometry to the element to prepare all the geometries of the stencil
+            problem.model().fvGeometries_().bind(element);
+
             auto globalI = problem.elementMapper().index(element);
             const auto& neighborStencil = this->model_().stencils(element).neighborStencil();
 
             assemblyMap_[globalI].reserve(neighborStencil.size());
-
-            for (auto&& globalJ : neighborStencil)
+            for (auto globalJ : neighborStencil)
             {
+                const auto& fvGeometry = this->model_().fvGeometries(globalJ);
+                const auto& elementJ = this->model_().fvGeometries().element(globalJ);
+
                 // find the flux vars needed for the calculation of the flux into element
                 std::vector<IndexType> fluxVarIndices;
-                for (auto&& scvFaceJ : this->model_().fvGeometries(globalJ).scvfs())
+                for (const auto& scvFaceJ : fvGeometry.scvfs())
                 {
                     auto fluxVarsIdx = scvFaceJ.index();
 
                     // if globalI is in flux var stencil, add to list
-                    if (this->model_().fluxVars(fluxVarsIdx).stencil().count(globalI))
-                        fluxVarIndices.push_back(fluxVarsIdx);
+                    FluxVariables fluxVars;
+                    fluxVars.init(problem, elementJ, scvFaceJ);
+                    const auto& fluxStencil = fluxVars.stencil();
+                    for (auto globalIdx : fluxStencil)
+                        if (globalIdx == globalI)
+                            fluxVarIndices.push_back(fluxVarsIdx);
                 }
 
                 assemblyMap_[globalI].emplace_back(std::move(fluxVarIndices));
@@ -145,20 +155,35 @@ public:
      */
     void assemble(const Element& element, JacobianMatrix& matrix)
     {
+        const bool isGhost = (element.partitionType() == Dune::GhostEntity);
+
+        // prepare the volvars/fvGeometries in case caching is disabled
+        this->model_().fvGeometries_().bind(element);
+        this->model_().curVolVars_().bind(element);
+        this->model_().prevVolVars_().bindElement(element);
+
         // set the current grid element and update the element's
         // finite volume geometry
         globalI_ = this->problem_().elementMapper().index(element);
 
-        bcTypes_.update(this->problem_(), element, fvElemGeom_());
+        const auto& fvGeometry = fvElemGeom_();
+        bcTypes_.update(this->problem_(), element, fvGeometry);
 
         // calculate the local residual
-        this->localResidual().eval(element, bcTypes_);
-        this->residual_ = this->localResidual().residual();
+        if (isGhost)
+        {
+            this->residual_ = 0.0;
+        }
+        else
+        {
+            this->localResidual().eval(element, bcTypes_);
+            this->residual_ = this->localResidual().residual();
+        }
 
-        this->model_().updatePVWeights(fvElemGeom_());
+        this->model_().updatePVWeights(fvGeometry);
 
         // calculate derivatives of all dofs in stencil with respect to the dofs in the element
-        evalPartialDerivatives_(element, matrix);
+        evalPartialDerivatives_(element, matrix, isGhost);
 
         // TODO: calculate derivatives in the case of an extended source stencil
         // const auto& extendedSourceStencil = model_().stencils(element).extendedSourceStencil();
@@ -175,14 +200,14 @@ public:
         // }
     }
 
-    void evalPartialDerivatives_(const Element& element, JacobianMatrix& matrix)
+    void evalPartialDerivatives_(const Element& element, JacobianMatrix& matrix, const bool isGhost)
     {
         // get stencil informations
         const auto& neighborStencil = this->model_().stencils(element).neighborStencil();
         const auto numNeighbors = neighborStencil.size();
 
         // derivatives in the element and the neighbors
-        PrimaryVariables partialDeriv(0.0);
+        PrimaryVariables partialDeriv(isGhost ? 1.0 : 0.0);
         Dune::BlockVector<PrimaryVariables> neighborDeriv(numNeighbors);
 
         // the localresidual class used for the flux calculations
@@ -232,11 +257,14 @@ public:
                 // update the volume variables
                 this->model_().curVolVars_(scv).update(priVars, this->problem_(), element, scv);
 
-                // calculate the residual with the deflected primary variables
-                this->localResidual().eval(element, bcTypes_);
+                if (!isGhost)
+                {
+                    // calculate the residual with the deflected primary variables
+                    this->localResidual().eval(element, bcTypes_);
 
-                // store the residual and the storage term
-                partialDeriv = this->localResidual().residual(0);
+                    // 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)
@@ -250,7 +278,8 @@ public:
                 // 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)
-                partialDeriv = this->residual_[0];
+                if (!isGhost)
+                    partialDeriv = this->residual_[0];
                 neighborDeriv = origFlux;
             }
 
@@ -266,11 +295,14 @@ public:
                 // update the volume variables
                 this->model_().curVolVars_(scv).update(priVars, this->problem_(), element, scv);
 
-                // calculate the residual with the deflected primary variables
-                this->localResidual().eval(element, bcTypes_);
+                if (!isGhost)
+                {
+                    // calculate the residual with the deflected primary variables
+                    this->localResidual().eval(element, bcTypes_);
 
-                // subtract the residual from the derivative storage
-                partialDeriv -= this->localResidual().residual(0);
+                    // 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)
@@ -284,25 +316,19 @@ public:
                 // 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)
-                partialDeriv -= this->residual_[0];
+                if (!isGhost)
+                    partialDeriv -= this->residual_[0];
                 neighborDeriv -= origFlux;
             }
 
             // divide difference in residuals by the magnitude of the
             // deflections between the two function evaluation
-            partialDeriv /= delta;
+            if (!isGhost)
+                partialDeriv /= delta;
             neighborDeriv /= delta;
 
             // restore the original state of the scv's volume variables
-            this->model_().curVolVars_(scv) = origVolVars;
-
-#if HAVE_VALGRIND
-            for (unsigned i = 0; i < partialDeriv.size(); ++i)
-                Valgrind::CheckDefined(partialDeriv[i]);
-            for (int i = 0; i < neighborDeriv.size(); ++idx)
-                for(unsigned k = 0; k < neighborDeriv[i].size(); ++k)
-                    Valgrind::CheckDefined(neighborDeriv[i][k]);
-#endif
+            this->model_().curVolVars_(scv) = std::move(origVolVars);
 
             // update the global jacobian matrix with the current partial derivatives
             this->updateGlobalJacobian_(matrix, globalI_, globalI_, pvIdx, partialDeriv);
diff --git a/dumux/implicit/cellcentered/localresidual.hh b/dumux/implicit/cellcentered/localresidual.hh
index a543c8932cda9c1b9eb25348e0ebdffc7a9c6c85..c9e0b970b9ac76e3f7b0568949a70dfd7af3a5e8 100644
--- a/dumux/implicit/cellcentered/localresidual.hh
+++ b/dumux/implicit/cellcentered/localresidual.hh
@@ -51,6 +51,10 @@ class CCLocalResidual : public ImplicitLocalResidual<TypeTag>
         numEq = GET_PROP_VALUE(TypeTag, NumEq)
     };
 
+    enum { constantBC = GET_PROP_VALUE(TypeTag, ConstantBoundaryConditions) };
+
+    enum { enableVolVarsCache = GET_PROP_VALUE(TypeTag, EnableVolumeVariablesCache) };
+
     typedef typename GridView::template Codim<0>::Entity Element;
     typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
     typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
@@ -69,9 +73,7 @@ protected:
     PrimaryVariables computeFlux_(const SubControlVolumeFace &scvf)
     {
         if (!scvf.boundary() /*TODO: || GET_PROP_VALUE(TypeTag, BoundaryReconstruction)*/)
-        {
             return this->asImp_().computeFlux(scvf);
-        }
         else
             return this->asImp_().evalBoundary_(scvf);
     }
@@ -215,6 +217,15 @@ protected:
         // temporary vector to store the Dirichlet boundary fluxes
         PrimaryVariables flux(0);
 
+        if (!constantBC || !enableVolVarsCache)
+        {
+            // update corresponding boundary volume variables before flux calculation
+            const auto insideScvIdx = scvf.insideScvIdx();
+            const auto& insideScv = this->model_().fvGeometries().subControlVolume(insideScvIdx);
+            const auto dirichletPriVars = this->problem_().dirichlet(this->element_(), scvf);
+            this->model_().curVolVars_(scvf.outsideScvIdx()).update(dirichletPriVars, this->problem_(), this->element_(), insideScv);
+        }
+
         auto dirichletFlux = this->asImp_().computeFlux(scvf);
 
         // add fluxes to the residual
@@ -237,7 +248,7 @@ protected:
         PrimaryVariables dirichletValues = this->problem_().dirichlet(this->element_(), scvf);
 
         // get the primary variables
-        const auto& priVars = this->problem_().model().curVolVars(scvf.insideScvIdx()).priVars();
+        const auto& priVars = this->model_().curVolVars(scvf.insideScvIdx()).priVars();
 
         // set Dirichlet conditions in a strong sense
         for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
diff --git a/dumux/implicit/cellcentered/stencils.hh b/dumux/implicit/cellcentered/stencils.hh
index ad6f15b3e37b75c54782b34c0b6ac93096159a39..c0dcfea856e0e80ad985dbe6443adccbc6efa822 100644
--- a/dumux/implicit/cellcentered/stencils.hh
+++ b/dumux/implicit/cellcentered/stencils.hh
@@ -37,20 +37,42 @@ class CCElementStencils
 {
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
     using IndexType = typename GridView::IndexSet::IndexType;
     using Element = typename GridView::template Codim<0>::Entity;
     // TODO a separate stencil class all stencils can derive from?
-    using Stencil = std::set<IndexType>;
+    using Stencil = std::vector<IndexType>;
 public:
     void update(const Problem& problem, const Element& element)
     {
-        for (auto&& scvf : problem.model().fvGeometries(element).scvfs())
+        elementStencil_.clear();
+
+        const auto& fvGeometry = problem.model().fvGeometries(element);
+        // loop over sub control faces
+        for (const auto& scvf : fvGeometry.scvfs())
         {
-            auto&& fluxStencil = problem.model().fluxVars(scvf).stencil();
-            elementStencil_.insert(fluxStencil.begin(), fluxStencil.end());
+            FluxVariables fluxVars;
+            const auto& stencil = fluxVars.computeStencil(problem, scvf);
+
+            elementStencil_.insert(elementStencil_.end(), stencil.begin(), stencil.end());
         }
+        // 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_;
-        neighborStencil_.erase(problem.elementMapper().index(element));
+
+        //remove the element itself and possible ghost neighbors from the stencil
+        auto pred = [&problem, globalI](const int i) -> bool
+        {
+            if (i == globalI)
+                return true;
+            if (problem.model().fvGeometries().element(i).partitionType() == Dune::GhostEntity)
+                return true;
+            return false;
+        };
+        neighborStencil_.erase(std::remove_if(neighborStencil_.begin(), neighborStencil_.end(), pred), neighborStencil_.end());
     }
 
     //! The full element stencil (all element this element is interacting with)
@@ -59,7 +81,7 @@ public:
         return elementStencil_;
     }
 
-    //! //! The full element stencil without this element
+    //! The full element stencil without this element
     const Stencil& neighborStencil() const
     {
         return neighborStencil_;
@@ -79,14 +101,18 @@ class CCStencilsVector
 {
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using StencilType = CCElementStencils<TypeTag>;
 
 public:
-    void update(const Problem& problem)
+    void update(Problem& problem)
     {
         problemPtr_ = &problem;
         elementStencils_.resize(problem.gridView().size(0));
         for (const auto& element : elements(problem.gridView()))
         {
+            // bind the FvGeometry to the element before using it
+            problem.model().fvGeometries_().bind(element);
+
             auto eIdx = problem.elementMapper().index(element);
             elementStencils_[eIdx].update(problem, element);
         }
@@ -94,7 +120,7 @@ public:
 
     //! overload for elements
     template <class Entity>
-    typename std::enable_if<Entity::codimension == 0, const CCElementStencils<TypeTag>&>::type
+    typename std::enable_if<Entity::codimension == 0, const StencilType&>::type
     get(const Entity& entity) const
     {
         return elementStencils_[problemPtr_->elementMapper().index(entity)];
diff --git a/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh b/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh
index 83a77b0bc517c68712152eaaa29a03e8eefa19b0..96081fd91c65e23e87a0e70946938753160690ea 100644
--- a/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh
+++ b/dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh
@@ -39,8 +39,13 @@ namespace Dumux
  *        This builds up the sub control volumes and sub control volume faces
  *        for each element.
  */
-template<class TypeTag>
+template<class TypeTag, bool enableFVElementGeometryCache>
 class CCTpfaFVElementGeometryVector
+{};
+
+// specialization in case the FVElementGeometries are stored
+template<class TypeTag>
+class CCTpfaFVElementGeometryVector<TypeTag, true>
 {
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
@@ -73,7 +78,7 @@ public:
     //! Get a sub control volume face with a global scvf index
     const SubControlVolumeFace& subControlVolumeFace(IndexType scvfIdx) const
     {
-        return *scvfs_[scvfIdx];
+        return scvfs_[scvfIdx];
     }
 
     //! The total number of sub control volumes
@@ -82,12 +87,18 @@ public:
         return scvs_.size();
     }
 
-    //! The total number of sun control volume faces
+    //! The total number of sub control volume faces
     std::size_t numScvf() const
     {
         return scvfs_.size();
     }
 
+    //! The total number of boundary sub control volume faces
+    std::size_t numBoundaryScvf() const
+    {
+        return numBoundaryScvf_;
+    }
+
     // Get an element from a sub control volume contained in it
     Element element(const SubControlVolume& scv) const
     { return elementMap_.element(scv.elementIndex()); }
@@ -106,6 +117,7 @@ public:
 
         // Build the SCV and SCV faces
         IndexType scvfIdx = 0;
+        numBoundaryScvf_ = 0;
         elementMap_.resize(gridView_.size(0));
         scvs_.resize(gridView_.size(0));
         for (const auto& element : elements(gridView_))
@@ -120,25 +132,27 @@ public:
             std::vector<IndexType> scvfsIndexSet;
             for (const auto& intersection : intersections(gridView_, element))
             {
-                if (!intersection.boundary())
+                // inner sub control volume faces
+                if (intersection.neighbor())
                 {
                     auto nIdx = problem.elementMapper().index(intersection.outside());
-                    scvfs_.push_back(std::make_shared<SubControlVolumeFace>(intersection.geometry(),
-                                                                            intersection.geometry().center(),
-                                                                            intersection.centerUnitOuterNormal(),
-                                                                            scvfIdx,
-                                                                            std::vector<IndexType>({eIdx, nIdx}),
-                                                                            false));
+                    scvfs_.push_back(SubControlVolumeFace(intersection.geometry(),
+                                                          intersection.geometry().center(),
+                                                          intersection.centerUnitOuterNormal(),
+                                                          scvfIdx,
+                                                          std::vector<IndexType>({eIdx, nIdx}),
+                                                          false));
                     scvfsIndexSet.push_back(scvfIdx++);
                 }
-                else
+                // boundary sub control volume faces
+                else if (intersection.boundary())
                 {
-                    scvfs_.push_back(std::make_shared<SubControlVolumeFace>(intersection.geometry(),
-                                                                            intersection.geometry().center(),
-                                                                            intersection.centerUnitOuterNormal(),
-                                                                            scvfIdx,
-                                                                            std::vector<IndexType>({eIdx}),
-                                                                            true));
+                    scvfs_.push_back(SubControlVolumeFace(intersection.geometry(),
+                                                          intersection.geometry().center(),
+                                                          intersection.centerUnitOuterNormal(),
+                                                          scvfIdx,
+                                                          std::vector<IndexType>({eIdx, gridView_.size(0) + numBoundaryScvf_++}),
+                                                          true));
                     scvfsIndexSet.push_back(scvfIdx++);
                 }
             }
@@ -148,11 +162,262 @@ public:
         }
     }
 
+    // Binding of an element, called by the local jacobian to prepare element assembly
+    // For compatibility reasons with the FVGeometry cache disabled
+    void bind(const Element& element) {}
+    void bindElement(const Element& element) {}
+
 private:
     GridView gridView_;
     Dumux::ElementMap<GridView> elementMap_;
     std::vector<std::shared_ptr<SubControlVolume>> scvs_;
-    std::vector<std::shared_ptr<SubControlVolumeFace>> scvfs_;
+    std::vector<SubControlVolumeFace> scvfs_;
+    std::vector<FVElementGeometry> fvGeometries_;
+    IndexType numBoundaryScvf_;
+};
+
+// specialization in case the FVElementGeometries are not stored
+template<class TypeTag>
+class CCTpfaFVElementGeometryVector<TypeTag, false>
+{
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using IndexType = typename GridView::IndexSet::IndexType;
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+    using Element = typename GridView::template Codim<0>::Entity;
+
+public:
+    //! Constructor
+    CCTpfaFVElementGeometryVector(const GridView gridView)
+    : gridView_(gridView), elementMap_(gridView), eIdxBound_(-1) {}
+
+    /* \brief Get the finite volume geometry of an element
+     * \note The finite volume geometry offers iterators over the sub control volumes
+     *       and the sub control volume faces of an element.
+     */
+    const FVElementGeometry& fvGeometry(IndexType eIdx) const
+    {
+        return fvGeometries_[getStencilScvIdx_(eIdx)];
+    }
+
+    const FVElementGeometry& fvGeometry(const Element& element) const
+    {
+        return fvGeometries_[getStencilScvIdx_(problem_().elementMapper().index(element))];
+    }
+
+    //! Get a sub control volume with a global scv index
+    const SubControlVolume& subControlVolume(IndexType scvIdx) const
+    {
+        return stencilScvs_[getStencilScvIdx_(scvIdx)];
+    }
+
+    //! Get a sub control volume face with a global scvf index
+    const SubControlVolumeFace& subControlVolumeFace(IndexType scvfIdx) const
+    {
+        return stencilScvfs_[getStencilScvFaceIdx_(scvfIdx)];
+    }
+
+    //! The total number of sub control volumes
+    std::size_t numScv() const
+    {
+        return numScvs_;
+    }
+
+    //! The total number of sub control volume faces
+    std::size_t numScvf() const
+    {
+        return numScvf_;
+    }
+
+    //! The total number of boundary sub control volume faces
+    std::size_t numBoundaryScvf() const
+    {
+        return numBoundaryScvf_;
+    }
+
+    // Get an element from a sub control volume contained in it
+    Element element(const SubControlVolume& scv) const
+    { return elementMap_.element(scv.elementIndex()); }
+
+    // Get an element from a global element index
+    Element element(IndexType eIdx) const
+    { return elementMap_.element(eIdx); }
+
+    //! update all fvElementGeometries (do this again after grid adaption)
+    void update(const Problem& problem)
+    {
+        problemPtr_ = &problem;
+        elementMap_.clear();
+        scvFaceToElementMap_.clear();
+        eIdxBound_ = -1;
+
+        // Build the SCV and SCV faces
+        numScvf_ = 0;
+        numBoundaryScvf_ = 0;
+        numScvs_ = gridView_.size(0);
+        elementMap_.resize(numScvs_);
+        scvFaceIndices_.resize(numScvs_);
+        neighborVolVarIndices_.resize(numScvs_);
+        for (const auto& element : elements(gridView_))
+        {
+            auto eIdx = problem.elementMapper().index(element);
+
+            // fill the element map with seeds
+            elementMap_[eIdx] = element.seed();
+
+            // the element-wise index sets for finite volume geometry
+            std::vector<IndexType> scvfsIndexSet;
+            std::vector<IndexType> neighborVolVarIndexSet;
+            for (const auto& intersection : intersections(gridView_, element))
+            {
+                // inner sub control volume faces
+                if (intersection.neighbor())
+                {
+                    scvFaceToElementMap_.push_back(eIdx);
+                    scvfsIndexSet.push_back(numScvf_++);
+                    auto nIdx = problem.elementMapper().index(intersection.outside());
+                    neighborVolVarIndexSet.push_back(nIdx);
+                }
+                // boundary sub control volume faces
+                else if (intersection.boundary())
+                {
+                    scvFaceToElementMap_.push_back(eIdx);
+                    scvfsIndexSet.push_back(numScvf_++);
+                    neighborVolVarIndexSet.push_back(numScvs_ + numBoundaryScvf_++);
+                }
+            }
+
+            // store the sets of indices in the data container
+            scvFaceIndices_[eIdx] = scvfsIndexSet;
+            neighborVolVarIndices_[eIdx] = neighborVolVarIndexSet;
+        }
+    }
+
+    // Binding of an element preparing the geometries of the whole stencil
+    // called by the local jacobian to prepare element assembly
+    void bind(const Element& element)
+    {
+        const auto eIdx = problem_().elementMapper().index(element);
+        if (eIdxBound_ == eIdx && fvGeometries_.size() > 1)
+            return;
+
+        release_();
+        eIdxBound_ = eIdx;
+
+        makeElementGeometries_(element);
+        for (const auto& intersection : intersections(gridView_, element))
+        {
+            if (intersection.neighbor())
+                makeElementGeometries_(intersection.outside());
+        }
+    }
+
+    // Binding of an element preparing the geometries only inside the element
+    void bindElement(const Element& element)
+    {
+        const auto eIdx = problem_().elementMapper().index(element);
+        if (eIdx == eIdxBound_ || std::find(stencilScvIndices_.begin(), stencilScvIndices_.end(), eIdx) != stencilScvIndices_.end() )
+            return;
+
+        release_();
+        eIdxBound_ = eIdx;
+
+        makeElementGeometries_(element);
+    }
+
+
+private:
+
+    void release_()
+    {
+        stencilScvs_.clear();
+        stencilScvIndices_.clear();
+        stencilScvfs_.clear();
+        stencilScvFaceIndices_.clear();
+        fvGeometries_.clear();
+    }
+
+    void makeElementGeometries_(const Element& element)
+    {
+        const auto eIdx = problem_().elementMapper().index(element);
+
+        std::vector<SubControlVolume> scv({SubControlVolume(element.geometry(), eIdx)});
+        stencilScvs_.push_back(scv[0]);
+        stencilScvIndices_.push_back(eIdx);
+
+        const auto& scvFaceIndices = scvFaceIndices_[eIdx];
+        const auto& neighborVolVarIndices = neighborVolVarIndices_[eIdx];
+
+        int scvfCounter = 0;
+        for (const auto& intersection : intersections(gridView_, element))
+        {
+            if (intersection.neighbor() || intersection.boundary())
+            {
+                stencilScvfs_.push_back(SubControlVolumeFace(intersection.geometry(),
+                                                             intersection.geometry().center(),
+                                                             intersection.centerUnitOuterNormal(),
+                                                             scvFaceIndices[scvfCounter],
+                                                             std::vector<IndexType>({eIdx, neighborVolVarIndices[scvfCounter]}),
+                                                             intersection.boundary()));
+
+                stencilScvFaceIndices_.push_back(scvFaceIndices[scvfCounter]);
+                scvfCounter++;
+            }
+        }
+
+        // Compute the finite volume element geometries
+        fvGeometries_.push_back(FVElementGeometry(*this, {eIdx}, scvFaceIndices));
+    }
+
+    const int getStencilScvIdx_(const int scvIdx) const
+    {
+        auto it = std::find(stencilScvIndices_.begin(), stencilScvIndices_.end(), scvIdx);
+
+        if (it != stencilScvIndices_.end())
+            return std::distance(stencilScvIndices_.begin(), it);
+        else
+            DUNE_THROW(Dune::InvalidStateException,
+                       "Could not find the sub control volume for scvIdx = " << scvIdx <<
+                       ", make sure to properly bind the FVGeometries to the element before using them.");
+    }
+
+    const int getStencilScvFaceIdx_(const int scvfIdx) const
+    {
+        auto it = std::find(stencilScvFaceIndices_.begin(), stencilScvFaceIndices_.end(), scvfIdx);
+
+        if (it != stencilScvFaceIndices_.end())
+            return std::distance(stencilScvFaceIndices_.begin(), it);
+        else
+            DUNE_THROW(Dune::InvalidStateException,
+                       "Could not find the sub control volume face for scvfIdx = " << scvfIdx <<
+                       ", make sure to properly bind the FVGeometries to the to the element before using them.");
+    }
+
+    const Problem& problem_() const
+    { return *problemPtr_; }
+
+    const Problem* problemPtr_;
+    GridView gridView_;
+
+    // Information on the global number of geometries
+    IndexType numScvs_;
+    IndexType numScvf_;
+    IndexType numBoundaryScvf_;
+
+    // vectors that stores global data
+    Dumux::ElementMap<GridView> elementMap_;
+    std::vector<std::vector<IndexType>> scvFaceIndices_;
+    std::vector<std::vector<IndexType>> neighborVolVarIndices_;
+    std::vector<IndexType> scvFaceToElementMap_;
+
+    // vectors to store the geometries temporarily after binding an element
+    IndexType eIdxBound_;
+    std::vector<SubControlVolume> stencilScvs_;
+    std::vector<IndexType> stencilScvIndices_;
+    std::vector<SubControlVolumeFace> stencilScvfs_;
+    std::vector<IndexType> stencilScvFaceIndices_;
     std::vector<FVElementGeometry> fvGeometries_;
 };
 
diff --git a/dumux/implicit/cellcentered/tpfa/propertydefaults.hh b/dumux/implicit/cellcentered/tpfa/propertydefaults.hh
index ed73e4960b8dd3e8f30bfcc558abe32494ab4e7c..348ad210a3d03da3f9c6cd0e0c0042c58f1bd2a6 100644
--- a/dumux/implicit/cellcentered/tpfa/propertydefaults.hh
+++ b/dumux/implicit/cellcentered/tpfa/propertydefaults.hh
@@ -29,6 +29,7 @@
 
 #include <dumux/implicit/propertydefaults.hh>
 #include <dumux/implicit/fvelementgeometry.hh>
+ #include <dumux/implicit/fluxvariablescache.hh>
 #include <dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh>
 #include <dumux/implicit/cellcentered/properties.hh>
 #include <dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh>
@@ -41,7 +42,7 @@ template<class TypeTag> class CCElementBoundaryTypes;
 
 namespace Properties {
 //! Set the default for the FVElementGeometry vector
-SET_TYPE_PROP(CCTpfaModel, FVElementGeometryVector, CCTpfaFVElementGeometryVector<TypeTag>);
+SET_TYPE_PROP(CCTpfaModel, FVElementGeometryVector, CCTpfaFVElementGeometryVector<TypeTag, GET_PROP_VALUE(TypeTag, EnableFVElementGeometryCache)>);
 
 //! The sub control volume
 SET_PROP(CCTpfaModel, SubControlVolume)
@@ -64,6 +65,9 @@ public:
     typedef Dumux::SubControlVolumeFace<ScvfGeometry, IndexType> type;
 };
 
+//! The flux variables cache class
+SET_TYPE_PROP(CCTpfaModel, FluxVariablesCache, Dumux::CCTpfaFluxVariablesCache<TypeTag>);
+
 //! The darcy flux variables
 SET_TYPE_PROP(CCTpfaModel, AdvectionType, Dumux::CCTpfaDarcysLaw<TypeTag>);
 
diff --git a/dumux/implicit/fluxvariables.hh b/dumux/implicit/fluxvariables.hh
index 3e1d6814715c4488b1fcda97fd5d06ac0bfe12ac..e7a451316110d2099af6d8d566d2d4bc0e461aeb 100644
--- a/dumux/implicit/fluxvariables.hh
+++ b/dumux/implicit/fluxvariables.hh
@@ -37,6 +37,82 @@ NEW_PROP_TAG(NumComponents);
 /*!
  * \ingroup ImplicitModel
  * \brief Base class for the flux variables
+ *        Actual flux variables inherit from this class
+ */
+template<class TypeTag, class Implementation>
+class FluxVariablesBase
+{
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Element = typename GridView::template Codim<0>::Entity;
+    using IndexType = typename GridView::IndexSet::IndexType;
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Stencil = std::vector<IndexType>;
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+
+    enum{ enableFluxVarsCache = GET_PROP_VALUE(TypeTag, EnableFluxVariablesCache) };
+
+public:
+    FluxVariablesBase() : problemPtr_(nullptr), scvFacePtr_(nullptr)
+    {}
+
+    void init(const Problem& problem,
+              const Element& element,
+              const SubControlVolumeFace &scvFace)
+    {
+        problemPtr_ = &problem;
+        scvFacePtr_ = &scvFace;
+
+        // update the stencil if needed
+        if (!enableFluxVarsCache)
+            stencil_ = asImp_().computeStencil(problem, scvFace);
+    }
+
+    // when caching is enabled, get the stencil from the cache class
+    template <typename T = TypeTag>
+    const typename std::enable_if<GET_PROP_VALUE(T, EnableFluxVariablesCache), Stencil>::type& stencil() const
+    { return problem().model().fluxVarsCache(scvFace()).stencil(); }
+
+    // when caching is disabled, return the private stencil variable. The update(...) routine has to be called beforehand.
+    template <typename T = TypeTag>
+    const typename std::enable_if<!GET_PROP_VALUE(T, EnableFluxVariablesCache), Stencil>::type& stencil()
+    { return stencil_; }
+
+    const Problem& problem() const
+    { return *problemPtr_; }
+
+    const SubControlVolumeFace& scvFace() const
+    {
+        return *scvFacePtr_;
+    }
+
+    Stencil computeStencil(const Problem& problem, const SubControlVolumeFace& scvFace)
+    { DUNE_THROW(Dune::InvalidStateException, "computeStencil() routine is not provided by the implementation."); }
+
+private:
+
+    Implementation &asImp_()
+    {
+        assert(static_cast<Implementation*>(this) != 0);
+        return *static_cast<Implementation*>(this);
+    }
+
+    const Implementation &asImp_() const
+    {
+        assert(static_cast<const Implementation*>(this) != 0);
+        return *static_cast<const Implementation*>(this);
+    }
+
+    const Problem *problemPtr_;              //! Pointer to the problem
+    const SubControlVolumeFace *scvFacePtr_; //! Pointer to the sub control volume face for which the flux variables are created
+    Stencil stencil_;                        //! The flux stencil
+};
+
+
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief the flux variables class
  *        specializations are provided for combinations of physical processes
  */
 template<class TypeTag, bool enableAdvection, bool enableMolecularDiffusion, bool enableEnergyBalance>
@@ -45,65 +121,64 @@ class FluxVariables {};
 
 // specialization for pure advective flow (e.g. 1p/2p/3p immiscible darcy flow)
 template<class TypeTag>
-class FluxVariables<TypeTag, true, false, false>
+class FluxVariables<TypeTag, true, false, false> : public FluxVariablesBase<TypeTag, FluxVariables<TypeTag, true, false, false>>
 {
+    using ParentType = FluxVariablesBase<TypeTag, FluxVariables<TypeTag, true, false, false>>;
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using Element = typename GridView::template Codim<0>::Entity;
     using IndexType = typename GridView::IndexSet::IndexType;
-    using Stencil = std::set<IndexType>;
-    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Stencil = std::vector<IndexType>;
     using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
     using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
 
+    enum
+    {
+        isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox),
+        enableFluxVarsCache = GET_PROP_VALUE(TypeTag, EnableFluxVariablesCache),
+        constantBC = GET_PROP_VALUE(TypeTag, ConstantBoundaryConditions)
+    };
+
 public:
-    void update(const Problem& problem,
-                const Element& element,
-                const SubControlVolumeFace &scvf)
-    {
-        if (scvf.boundary())
-        {
-            if(!boundaryVolVars_)
-                boundaryVolVars_ = Dune::Std::make_unique<VolumeVariables>();
-            advection_.update(problem, element, scvf, boundaryVolVars_.get());
-        }
-        else
-        {
-            advection_.update(problem, element, scvf);
-        }
-    }
 
-    void beginFluxComputation()
+    void initAndComputeFluxes(const Problem& problem,
+                              const Element& element,
+                              const SubControlVolumeFace &scvFace)
     {
-        advection_.beginFluxComputation();
+        ParentType::init(problem, element, scvFace);
     }
 
-    const AdvectionType& advection() const
+    template<typename FunctionType>
+    Scalar advectiveFlux(const int phaseIdx, const FunctionType upwindFunction)
     {
-        return advection_;
-    }
+        Scalar flux = AdvectionType::flux(this->problem(), this->scvFace(), phaseIdx);
 
-    const Stencil& stencil() const
-    {
-        return advection_.stencil();
+        const auto& insideVolVars = this->problem().model().curVolVars(this->scvFace().insideScvIdx());
+        const auto& outsideVolVars = this->problem().model().curVolVars(this->scvFace().outsideScvIdx());
+
+        if (std::signbit(flux))
+            return flux*upwindFunction(outsideVolVars, insideVolVars);
+        else
+            return flux*upwindFunction(insideVolVars, outsideVolVars);
     }
 
-private:
-    AdvectionType advection_;
-    // boundary volume variables in case of Dirichlet boundaries
-    std::unique_ptr<VolumeVariables> boundaryVolVars_;
+    Stencil computeStencil(const Problem& problem, const SubControlVolumeFace& scvFace)
+    { return AdvectionType::stencil(problem, scvFace); }
 };
 
 
 // specialization for isothermal advection molecularDiffusion equations
 template<class TypeTag>
-class FluxVariables<TypeTag, true, true, false>
+class FluxVariables<TypeTag, true, true, false> : public FluxVariablesBase<TypeTag, FluxVariables<TypeTag, true, true, false>>
 {
+    using ParentType = FluxVariablesBase<TypeTag, FluxVariables<TypeTag, true, true, false>>;
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using Element = typename GridView::template Codim<0>::Entity;
     using IndexType = typename GridView::IndexSet::IndexType;
-    using Stencil = std::set<IndexType>;
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Stencil = std::vector<IndexType>;
     using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
     using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
     using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
@@ -116,93 +191,51 @@ class FluxVariables<TypeTag, true, true, false>
     };
 
 public:
-    void update(const Problem& problem,
-                const Element& element,
-                const SubControlVolumeFace &scvf)
-    {
-        if (scvf.boundary())
-        {
-            if(!boundaryVolVars_)
-                boundaryVolVars_ = Dune::Std::make_unique<VolumeVariables>();
-
-            advection_.update(problem, element, scvf, boundaryVolVars_.get());
-            for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-                for (int compIdx = 0; compIdx < numComponents; ++compIdx)
-                {
-                    if (phaseIdx != compIdx)
-                        molecularDiffusion(phaseIdx, compIdx).update(problem, element, scvf, phaseIdx, compIdx, boundaryVolVars_.get());
-                }
-        }
-        else
-        {
-            advection_.update(problem, element, scvf);
-            for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-                for (int compIdx = 0; compIdx < numComponents; ++compIdx)
-                {
-                    if (phaseIdx != compIdx)
-                        molecularDiffusion(phaseIdx, compIdx).update(problem, element, scvf, phaseIdx, compIdx);
-                }
-        }
-    }
 
-    void beginFluxComputation()
+    void initAndComputeFluxes(const Problem& problem,
+                              const Element& element,
+                              const SubControlVolumeFace &scvFace)
     {
-        advection_.beginFluxComputation();
+        ParentType::init(problem, element, scvFace);
+
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-        {
-            for (int compIdx = 0; compIdx < numComponents; ++compIdx)
-            {
-                if (phaseIdx != compIdx)
-                    molecularDiffusion(phaseIdx, compIdx).beginFluxComputation(true);
-            }
-        }
+            advectiveFluxes_[phaseIdx] = AdvectionType::flux(problem, scvFace, phaseIdx);
     }
 
-    // TODO update transmissibilities?
-
-    const AdvectionType& advection() const
+    Stencil computeStencil(const Problem& problem, const SubControlVolumeFace& scvFace)
     {
-        return advection_;
-    }
+        // unifiy advective and diffusive stencil
+        Stencil stencil = AdvectionType::stencil(problem, scvFace);
+        Stencil diffusionStencil = MolecularDiffusionType::stencil(problem, scvFace);
 
-    AdvectionType& advection()
-    {
-        return advection_;
-    }
+        stencil.insert(stencil.end(), diffusionStencil.begin(), diffusionStencil.end());
+        std::sort(stencil.begin(), stencil.end());
+        stencil.erase(std::unique(stencil.begin(), stencil.end()), stencil.end());
 
-    const MolecularDiffusionType& molecularDiffusion(const int phaseIdx, const int compIdx) const
-    {
-        if (compIdx < phaseIdx)
-            return molecularDiffusion_[phaseIdx][compIdx];
-        else if (compIdx > phaseIdx)
-            return molecularDiffusion_[phaseIdx][compIdx-1];
-        else
-            DUNE_THROW(Dune::InvalidStateException, "Diffusion flux called for phaseIdx = compIdx");
+        return stencil;
     }
 
-    MolecularDiffusionType& molecularDiffusion(const int phaseIdx, const int compIdx)
+    template<typename FunctionType>
+    Scalar advectiveFlux(const int phaseIdx, const FunctionType upwindFunction)
     {
-        if (compIdx < phaseIdx)
-            return molecularDiffusion_[phaseIdx][compIdx];
-        else if (compIdx > phaseIdx)
-            return molecularDiffusion_[phaseIdx][compIdx-1];
+        const auto& insideVolVars = this->problem().model().curVolVars(this->scvFace().insideScvIdx());
+        const auto& outsideVolVars = this->problem().model().curVolVars(this->scvFace().outsideScvIdx());
+
+        if (std::signbit(advectiveFluxes_[phaseIdx]))
+            return advectiveFluxes_[phaseIdx]*upwindFunction(outsideVolVars, insideVolVars);
         else
-            DUNE_THROW(Dune::InvalidStateException, "Diffusion flux called for phaseIdx = compIdx");
+            return advectiveFluxes_[phaseIdx]*upwindFunction(insideVolVars, outsideVolVars);
     }
 
-    Stencil stencil() const
+    Scalar molecularDiffusionFlux(const int phaseIdx, const int compIdx)
     {
-        // std::set<IndexType> stencil(advection().stencil().begin(), advection().stencil().end());
-        // TODO: insert all molecularDiffusion stencils of all components in all phases
-        // stencil.insert();
-        return advection().stencil();
+        Scalar flux = MolecularDiffusionType::flux(this->problem(), this->scvFace(), phaseIdx, compIdx);
+        return flux;
     }
 
 private:
-    AdvectionType advection_;
-    std::array< std::array<MolecularDiffusionType, numComponents-1>, numPhases> molecularDiffusion_;
-    // boundary volume variables in case of Dirichlet boundaries
-    std::unique_ptr<VolumeVariables> boundaryVolVars_;
+    // storage for calculated advective fluxes to not having to calculate them again
+    std::array<Scalar, numPhases> advectiveFluxes_;
 };
 
 
diff --git a/dumux/implicit/fluxvariablescache.hh b/dumux/implicit/fluxvariablescache.hh
new file mode 100644
index 0000000000000000000000000000000000000000..153b2499126846fe2efe224040d394667826d2fd
--- /dev/null
+++ b/dumux/implicit/fluxvariablescache.hh
@@ -0,0 +1,118 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for the flux variables
+ */
+#ifndef DUMUX_IMPLICIT_FLUXVARIABLESCACHE_HH
+#define DUMUX_IMPLICIT_FLUXVARIABLESCACHE_HH
+
+#include <dumux/implicit/properties.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+NEW_PROP_TAG(NumPhases);
+NEW_PROP_TAG(NumComponents);
+}
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief The flux variables cache classes
+ *        stores the transmissibilities and stencils
+ */
+// specialization for the Box Method
+template<class TypeTag>
+class BoxFluxVariablesCache
+{
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+    using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
+    using Element = typename GridView::template Codim<0>::Entity;
+    using IndexType = typename GridView::IndexSet::IndexType;
+    using Stencil = std::vector<IndexType>;
+    using TransmissibilityVector = std::vector<IndexType>;
+
+public:
+    void update(const Problem& problem,
+                const Element& element,
+                const SubControlVolumeFace &scvFace)
+    {
+        stencil_ = AdvectionType::stencil(problem, scvFace);
+        volVarsStencil_ = AdvectionType::volVarsStencil(problem, element, scvFace);
+        tij_ = AdvectionType::calculateTransmissibilities(problem, scvFace);
+    }
+
+    const Stencil& stencil() const
+    { return stencil_; }
+
+    const Stencil& darcyStencil() const
+    { return volVarsStencil_; }
+
+    const TransmissibilityVector& tij() const
+    { return tij_; }
+
+private:
+    Stencil volVarsStencil_;
+    Stencil stencil_;
+    TransmissibilityVector tij_;
+};
+
+// specialization for the cell centered tpfa method
+template<class TypeTag>
+class CCTpfaFluxVariablesCache
+{
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+    using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
+    using Element = typename GridView::template Codim<0>::Entity;
+    using IndexType = typename GridView::IndexSet::IndexType;
+    using Stencil = std::vector<IndexType>;
+
+public:
+    void update(const Problem& problem,
+                const Element& element,
+                const SubControlVolumeFace &scvFace)
+    {
+        FluxVariables fluxVars;
+        stencil_ = fluxVars.computeStencil(problem, scvFace);
+        tij_ = AdvectionType::calculateTransmissibilities(problem, scvFace);
+    }
+
+    const Stencil& stencil() const
+    { return stencil_; }
+
+    const Scalar& tij() const
+    { return tij_; }
+
+private:
+    Stencil stencil_;
+    Scalar tij_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/implicit/fluxvariablescachevector.hh b/dumux/implicit/fluxvariablescachevector.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e426ac6dd24da5c4e82b797f6e29b2d4bb14458c
--- /dev/null
+++ b/dumux/implicit/fluxvariablescachevector.hh
@@ -0,0 +1,81 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Base class for the volume variables vector
+ */
+#ifndef DUMUX_IMPLICIT_FLUXVARSCACHEVECTOR_HH
+#define DUMUX_IMPLICIT_FLUXVARSCACHEVECTOR_HH
+
+#include <dumux/implicit/properties.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup ImplicitModel
+ * \brief Base class for the flux variables cache vector, we store one cache per face
+ */
+template<class TypeTag>
+class FluxVariablesCacheVector
+{
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using IndexType = typename GridView::IndexSet::IndexType;
+    using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+
+public:
+    template <typename T = TypeTag>
+    typename std::enable_if<GET_PROP_VALUE(T, EnableFluxVariablesCache)>::type update(Problem& problem)
+    {
+        fluxVarsCache_.resize(problem.model().fvGeometries().numScvf());
+        for (const auto& element : elements(problem.gridView()))
+        {
+            // bind the geometries and volume variables to the element (all the elements in stencil)
+            problem.model().fvGeometries_().bind(element);
+            problem.model().curVolVars_().bind(element);
+
+            const auto& fvGeometry = problem.model().fvGeometries(element);
+            for (const auto& scvf : fvGeometry.scvfs())
+            {
+                (*this)[scvf.index()].update(problem, element, scvf);
+            }
+        }
+    }
+
+    template <typename T = TypeTag>
+    typename std::enable_if<!GET_PROP_VALUE(T, EnableFluxVariablesCache)>::type update(Problem& problem)
+    {}
+
+    template <typename T = TypeTag>
+    const typename std::enable_if<GET_PROP_VALUE(T, EnableFluxVariablesCache), FluxVariablesCache>::type& operator [](IndexType scvfIdx) const
+    { return fluxVarsCache_[scvfIdx]; }
+
+    template <typename T = TypeTag>
+    typename std::enable_if<GET_PROP_VALUE(T, EnableFluxVariablesCache), FluxVariablesCache>::type& operator [](IndexType scvfIdx)
+    { return fluxVarsCache_[scvfIdx]; }
+
+private:
+    std::vector<FluxVariablesCache> fluxVarsCache_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/implicit/localresidual.hh b/dumux/implicit/localresidual.hh
index 1ab9f085ae0e88c3d7ddb231277a27d5c835cd8a..4f9502b9bea1a46d915800fd694d2860482c0f56 100644
--- a/dumux/implicit/localresidual.hh
+++ b/dumux/implicit/localresidual.hh
@@ -91,6 +91,11 @@ public:
      */
     void eval(const Element &element)
     {
+        // make sure FVElementGeometry and volume variables are bound to the element
+        model_().fvGeometries_().bind(element);
+        model_().curVolVars_().bind(element);
+        model_().prevVolVars_().bindElement(element);
+
         ElementBoundaryTypes bcTypes;
         bcTypes.update(problem_(), element, fvGeometry_());
 
@@ -110,6 +115,11 @@ public:
     {
         elemPtr_ = &element;
 
+        // make sure FVElementGeometry and volume variables are bound to the element
+        model_().fvGeometries_().bindElement(element);
+        model_().curVolVars_().bindElement(element);
+        model_().prevVolVars_().bindElement(element);
+
         ElementBoundaryTypes bcTypes;
         bcTypes.update(problem_(), element, fvGeometry_());
         bcTypesPtr_ = &bcTypes;
@@ -129,6 +139,11 @@ public:
     {
         elemPtr_ = &element;
 
+        // make sure FVElementGeometry and volume variables are bound to the element
+        model_().fvGeometries_().bind(element);
+        model_().curVolVars_().bind(element);
+        model_().prevVolVars_().bindElement(element);
+
         ElementBoundaryTypes bcTypes;
         bcTypes.update(problem_(), element, fvGeometry_());
 
@@ -158,6 +173,7 @@ public:
     void eval(const Element &element,
               const ElementBoundaryTypes &bcTypes)
     {
+        // At this point the fv geometry has to be bound already
         elemPtr_ = &element;
         bcTypesPtr_ = &bcTypes;
 
@@ -170,19 +186,7 @@ public:
         storageTerm_ = 0.0;
 
         asImp_().evalFluxes_();
-
-#if !defined NDEBUG && HAVE_VALGRIND
-        for (int i=0; i < fvGeometry_().numScv; i++)
-            Valgrind::CheckDefined(residual_[i]);
-#endif // HAVE_VALGRIND
-
-       asImp_().evalVolumeTerms_();
-
-#if !defined NDEBUG && HAVE_VALGRIND
-        for (int i=0; i < fvGeometry_().numScv; i++) {
-            Valgrind::CheckDefined(residual_[i]);
-        }
-#endif // HAVE_VALGRIND
+        asImp_().evalVolumeTerms_();
     }
 
     /*!
@@ -250,7 +254,8 @@ protected:
     void evalFluxes_()
     {
         // calculate the mass flux over the scv faces and subtract
-        for (auto&& scvf : fvGeometry_().scvfs())
+        const auto& fvGeometry = fvGeometry_();
+        for (const auto& scvf : fvGeometry.scvfs())
         {
             auto flux = asImp_().computeFlux_(scvf);
 
@@ -306,7 +311,8 @@ protected:
 
         // calculate the amount of conservation each quantity inside
         // all sub control volumes
-        for (auto&& scv : fvGeometry_().scvs())
+        const auto& fvGeometry = fvGeometry_();
+        for (const auto& scv : fvGeometry.scvs())
         {
             auto scvIdx = scv.indexInElement();
 
@@ -318,7 +324,8 @@ protected:
     PrimaryVariables evalSource_()
     {
         PrimaryVariables source(0);
-        for (auto&& scv : fvGeometry_().scvs())
+        const auto& fvGeometry = fvGeometry_();
+        for (const auto& scv : fvGeometry)
         {
             source += this->problem_().source(element_(), scv);
 
@@ -339,7 +346,8 @@ protected:
     evalVolumeTerms_()
     {
         // evaluate the volume terms (storage + source terms)
-        for (auto&& scv : fvGeometry_().scvs())
+        const auto& fvGeometry = fvGeometry_();
+        for (const auto& scv : fvGeometry.scvs())
         {
             auto scvIdx = scv.indexInElement();
             auto curExtrusionFactor = model_().curVolVars(scv).extrusionFactor();
@@ -362,7 +370,8 @@ protected:
     evalVolumeTerms_()
     {
         // evaluate the volume terms (storage + source terms)
-        for (auto&& scv : fvGeometry_().scvs())
+        const auto& fvGeometry = fvGeometry_();
+        for (const auto& scv : fvGeometry.scvs())
         {
             auto scvIdx = scv.indexInElement();
             auto prevExtrusionFactor = model_().prevVolVars(scv).extrusionFactor();
@@ -439,7 +448,7 @@ protected:
      * \brief Returns a reference to the current element's finite
      *        volume geometry.
      */
-    const FVElementGeometry &fvGeometry_() const
+    const FVElementGeometry& fvGeometry_() const
     {
         return model_().fvGeometries(element_());
     }
diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh
index 1979e3eaa8eb32db3b790192a7ac44565f0f1df0..d666958efe7910cfe7de58e2e6b8740157cb18d3 100644
--- a/dumux/implicit/model.hh
+++ b/dumux/implicit/model.hh
@@ -27,13 +27,19 @@
 #include <dune/istl/bvector.hh>
 
 #include "properties.hh"
+#include "localresidual.hh"
 #include <dumux/implicit/adaptive/gridadaptproperties.hh>
 #include <dumux/common/valgrind.hh>
 #include <dumux/parallel/vertexhandles.hh>
+ #include <dumux/porousmediumflow/compositional/primaryvariableswitch.hh>
 
 namespace Dumux
 {
 
+// forward declaration of the friend localresidual classes
+template <class TypeTag> class CCLocalResidual;
+template <class TypeTag> class BoxLocalResidual;
+
 /*!
  * \ingroup ImplicitModel
  * \brief The base class for the vertex centered finite volume
@@ -42,8 +48,17 @@ namespace Dumux
 template<class TypeTag>
 class ImplicitModel
 {
+    friend CCLocalResidual<TypeTag>;
+    friend BoxLocalResidual<TypeTag>;
+    friend ImplicitLocalResidual<TypeTag>;
+    friend PrimaryVariableSwitch<TypeTag>;
+    friend typename GET_PROP_TYPE(TypeTag, Problem);
     friend typename GET_PROP_TYPE(TypeTag, LocalJacobian);
-    friend typename GET_PROP_TYPE(TypeTag, LocalResidual);
+    friend typename GET_PROP_TYPE(TypeTag, StencilsVector);
+    friend typename GET_PROP_TYPE(TypeTag, CurrentVolumeVariablesVector);
+    friend typename GET_PROP_TYPE(TypeTag, PreviousVolumeVariablesVector);
+    friend typename GET_PROP_TYPE(TypeTag, FluxVariablesCacheVector);
+
     typedef typename GET_PROP_TYPE(TypeTag, Model) Implementation;
     typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
     typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
@@ -53,12 +68,13 @@ class ImplicitModel
     typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
     typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
     typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler;
-    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariablesVector) VolumeVariablesVector;
+    typedef typename GET_PROP_TYPE(TypeTag, CurrentVolumeVariablesVector) CurrentVolumeVariablesVector;
+    typedef typename GET_PROP_TYPE(TypeTag, PreviousVolumeVariablesVector) PreviousVolumeVariablesVector;
     typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
     typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
     typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
-    typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, FluxVariablesVector) FluxVariablesVector;
+    typedef typename GET_PROP_TYPE(TypeTag, FluxVariablesCache) FluxVariablesCache;
+    typedef typename GET_PROP_TYPE(TypeTag, FluxVariablesCacheVector) FluxVariablesCacheVector;
     typedef typename GET_PROP_TYPE(TypeTag, StencilsVector) StencilsVector;
 
     enum {
@@ -103,8 +119,8 @@ public:
 
         updateBoundaryIndices_();
 
-        fvGeometries_ = std::make_shared<FVElementGeometryVector>(gridView_());
-        fvGeometries_->update(problem_());
+        fvGeometryVector_ = std::make_shared<FVElementGeometryVector>(gridView_());
+        fvGeometryVector_->update(problem_());
 
         int numDofs = asImp_().numDofs();
         uCur_.resize(numDofs);
@@ -116,15 +132,14 @@ public:
         asImp_().applyInitialSolution_();
 
         // resize and update the volVars with the initial solution
-        curVolVarsVector_.resize(fvGeometries().numScv());
         curVolVarsVector_.update(problem_(), curSol());
 
-        // update the flux vars (precompute transmissibilities)
-        fluxVarsVector_.update(problem_());
-
         // update stencils
         stencilsVector_.update(problem_());
 
+        // update the flux variables caches
+        fluxVarsCacheVector_.update(problem_());
+
         // initialize assembler and create matrix
         localJacobian_.init(problem_());
         jacAsm_ = std::make_shared<JacobianAssembler>();
@@ -740,7 +755,7 @@ public:
      */
     bool onBoundary(const SubControlVolume &scv) const
     {
-        return asImp_().onBoundary(onBoundary(scv.dofIndex()));
+        return asImp_().onBoundary(scv.dofIndex());
     }
 
     /*!
@@ -769,11 +784,11 @@ public:
     //                                         fvGeometry, scvIdx, fluidState);
     // }
 
-    const FluxVariables& fluxVars(const SubControlVolumeFace& scvf) const
-    { return fluxVarsVector_[scvf.index()]; }
+    const FluxVariablesCache& fluxVarsCache(const SubControlVolumeFace& scvf) const
+    { return fluxVarsCacheVector_[scvf.index()]; }
 
-    const FluxVariables& fluxVars(unsigned int scvfIdx) const
-    { return fluxVarsVector_[scvfIdx]; }
+    const FluxVariablesCache& fluxVarsCache(unsigned int scvfIdx) const
+    { return fluxVarsCacheVector_[scvfIdx]; }
 
     const VolumeVariables& curVolVars(const SubControlVolume& scv) const
     { return curVolVarsVector_[scv.index()]; }
@@ -788,37 +803,42 @@ public:
     { return prevVolVarsVector_[scvIdx]; }
 
     const FVElementGeometryVector& fvGeometries() const
-    { return *fvGeometries_; }
+    { return *fvGeometryVector_; }
 
     const FVElementGeometry& fvGeometries(const Element& element) const
-    { return fvGeometries_->fvGeometry(elementMapper().index(element)); }
+    { return fvGeometryVector_->fvGeometry(elementMapper().index(element)); }
 
     const FVElementGeometry& fvGeometries(unsigned int eIdx) const
-    { return fvGeometries_->fvGeometry(eIdx); }
+    { return fvGeometryVector_->fvGeometry(eIdx); }
 
 protected:
-    // only to be called by friends and deriving classes
-    VolumeVariablesVector& curVolVars_()
+
+    CurrentVolumeVariablesVector& curVolVars_()
     { return curVolVarsVector_; }
 
-    // only to be called by friends and deriving classes
     VolumeVariables& curVolVars_(const SubControlVolume& scv)
     { return curVolVarsVector_[scv.index()]; }
 
     VolumeVariables& curVolVars_(unsigned int scvIdx)
     { return curVolVarsVector_[scvIdx]; }
 
+    PreviousVolumeVariablesVector& prevVolVars_()
+    { return prevVolVarsVector_; }
+
     VolumeVariables& prevVolVars_(const SubControlVolume& scv)
     { return prevVolVarsVector_[scv.index()]; }
 
     VolumeVariables& prevVolVars_(unsigned int scvIdx)
     { return prevVolVarsVector_[scvIdx]; }
 
-    FluxVariables& fluxVars_(const SubControlVolumeFace& scvf)
-    { return fluxVarsVector_[scvf.index()]; }
+    FluxVariablesCache& fluxVarsCache_(const SubControlVolumeFace& scvf)
+    { return fluxVarsCacheVector_[scvf.index()]; }
 
-    FluxVariables& fluxVars_(unsigned int scvfIdx)
-    { return fluxVarsVector_[scvfIdx]; }
+    FluxVariablesCache& fluxVarsCache_(unsigned int scvfIdx)
+    { return fluxVarsCacheVector_[scvfIdx]; }
+
+    FVElementGeometryVector& fvGeometries_()
+    { return *fvGeometryVector_; }
 
     /*!
      * \brief A reference to the problem on which the model is applied.
@@ -859,8 +879,12 @@ protected:
         // condition at the center of each sub control volume
         for (const auto& element : elements(gridView_()))
         {
-            // deal with the current element
-            for(auto&& scv : fvGeometries(element).scvs())
+            // deal with the current element, bind FVGeometry to it
+            fvGeometries_().bindElement(element);
+            const auto& fvGeometry = fvGeometries(element);
+
+            // loop over sub control volumes
+            for(const auto& scv : fvGeometry.scvs())
             {
                 // let the problem do the dirty work of nailing down
                 // the initial solution.
@@ -968,17 +992,17 @@ protected:
     SolutionVector uPrev_;
 
     // the current and previous variables (primary and secondary variables)
-    VolumeVariablesVector curVolVarsVector_;
-    VolumeVariablesVector prevVolVarsVector_;
+    CurrentVolumeVariablesVector curVolVarsVector_;
+    PreviousVolumeVariablesVector prevVolVarsVector_;
 
-    // the flux variables vector
-    FluxVariablesVector fluxVarsVector_;
+    // the flux variables cache vector vector
+    FluxVariablesCacheVector fluxVarsCacheVector_;
 
     // the stencils vector
     StencilsVector stencilsVector_;
 
     // the finite volume element geometries
-    std::shared_ptr<FVElementGeometryVector> fvGeometries_;
+    std::shared_ptr<FVElementGeometryVector> fvGeometryVector_;
 
     Dune::BlockVector<Dune::FieldVector<Scalar, 1> > boxVolume_;
 
diff --git a/dumux/implicit/properties.hh b/dumux/implicit/properties.hh
index 41e22f039f8d15d6a7f5d14afcfe8c26109c37e7..cb9e36510f2484f1befcc8204260b6268303c621 100644
--- a/dumux/implicit/properties.hh
+++ b/dumux/implicit/properties.hh
@@ -61,6 +61,7 @@ NEW_PROP_TAG(SubControlVolume);//!< The type of the sub control volume
 NEW_PROP_TAG(SubControlVolumeFace); //!< The type of the sub control volume face
 NEW_PROP_TAG(FVElementGeometry); //!< The type of the finite volume geometry (iterators over scvs, scvfs)
 NEW_PROP_TAG(FVElementGeometryVector); //!< The type of the finite volume geometry vector
+NEW_PROP_TAG(EnableFVElementGeometryCache); //! specifies if geometric data is be saved (faster, but more memory consuming)
 
 NEW_PROP_TAG(JacobianAssembler); //!< Assembles the global jacobian matrix
 NEW_PROP_TAG(JacobianMatrix); //!< Type of the global jacobian matrix
@@ -71,11 +72,16 @@ NEW_PROP_TAG(PrimaryVariables); //!< A vector of primary variables within a sub-
 NEW_PROP_TAG(SolutionVector); //!< Vector containing all primary variable vector of the grid
 NEW_PROP_TAG(ElementSolutionVector); //!< A vector of primary variables within a sub-control volume
 
+NEW_PROP_TAG(EnableVolumeVariablesCache); //!< If disabled, the volume variables are not stored (reduces memory, but is slower)
 NEW_PROP_TAG(VolumeVariables);  //!< The secondary variables within a sub-control volume
-NEW_PROP_TAG(VolumeVariablesVector);  //!< The type for a container of volume variables
+NEW_PROP_TAG(CurrentVolumeVariablesVector);  //!< The type for a container of current volume variables
+NEW_PROP_TAG(PreviousVolumeVariablesVector);  //!< The type for a container of previous volume variables
 NEW_PROP_TAG(FluxVariables); //!< Container storing the different types of flux variables
-NEW_PROP_TAG(FluxVariablesVector); //!< The global vector of flux variable containers
-NEW_PROP_TAG(BoundaryVariables); //!< Data required to calculate fluxes over boundary faces (outflow)
+NEW_PROP_TAG(FluxVariablesCache); //!< Stores data associated with flux vars (if enabled)
+NEW_PROP_TAG(FluxVariablesCacheVector); //!< The global vector of flux variable containers
+NEW_PROP_TAG(EnableFluxVariablesCache); //! specifies if data on flux vars should be saved (faster, but more memory consuming)
+NEW_PROP_TAG(BoundaryVariables); //!< Data required to calculate fluxes over boundary faces in cc models(outflow)
+NEW_PROP_TAG(ConstantBoundaryConditions); //!< boundary data is stored in case the BC are constant
 
 // Specify the forms of fluxes that should be considered in the model
 // also, specify their corresponding flux variables
diff --git a/dumux/implicit/propertydefaults.hh b/dumux/implicit/propertydefaults.hh
index 4a413661a11aff0f20026384d3ae3d7e8009bc9a..375b0c46ef061087aa5f49d5f6a88535648d6973 100644
--- a/dumux/implicit/propertydefaults.hh
+++ b/dumux/implicit/propertydefaults.hh
@@ -42,7 +42,8 @@
 #include "volumevariables.hh"
 #include "volumevariablesvector.hh"
 #include "fluxvariables.hh"
-#include "fluxvariablesvector.hh"
+#include "fluxvariablescache.hh"
+#include "fluxvariablescachevector.hh"
 #include "fvelementgeometry.hh"
 
 namespace Dumux {
@@ -95,11 +96,14 @@ SET_TYPE_PROP(ImplicitBase, FVElementGeometry, Dumux::FVElementGeometry<TypeTag>
 //! The volume variable class, to be overloaded by the model
 SET_TYPE_PROP(ImplicitBase, VolumeVariables, Dumux::ImplicitVolumeVariables<TypeTag>);
 
-//! The global volume variables vector class
-SET_TYPE_PROP(ImplicitBase, VolumeVariablesVector, Dumux::VolumeVariablesVector<TypeTag>);
+//! The global current volume variables vector class
+SET_TYPE_PROP(ImplicitBase, CurrentVolumeVariablesVector, Dumux::VolumeVariablesVector<TypeTag, false, GET_PROP_VALUE(TypeTag, EnableVolumeVariablesCache)>);
 
-//! The global volume variables vector class
-SET_TYPE_PROP(ImplicitBase, FluxVariablesVector, Dumux::FluxVariablesVector<TypeTag>);
+//! The global previous volume variables vector class
+SET_TYPE_PROP(ImplicitBase, PreviousVolumeVariablesVector, Dumux::VolumeVariablesVector<TypeTag, true, GET_PROP_VALUE(TypeTag, EnableVolumeVariablesCache)>);
+
+//! The global flux variables cache vector class
+SET_TYPE_PROP(ImplicitBase, FluxVariablesCacheVector, Dumux::FluxVariablesCacheVector<TypeTag>);
 
 // //! The class that contains the different flux variables (i.e. darcy, diffusion, energy)
 SET_PROP(ImplicitBase, FluxVariables)
@@ -142,6 +146,18 @@ SET_BOOL_PROP(ImplicitBase, ImplicitEnableJacobianRecycling, false);
 //! disable partial reassembling by default
 SET_BOOL_PROP(ImplicitBase, ImplicitEnablePartialReassemble, false);
 
+//! We do not store the FVGeometry by default
+SET_BOOL_PROP(ImplicitBase, EnableFVElementGeometryCache, false);
+
+//! We do not store the volume variables by default
+SET_BOOL_PROP(ImplicitBase, EnableVolumeVariablesCache, false);
+
+//! disable flux variables data caching by default
+SET_BOOL_PROP(ImplicitBase, EnableFluxVariablesCache, false);
+
+//! boundary conditions are not stationary by default
+SET_BOOL_PROP(ImplicitBase, ConstantBoundaryConditions, false);
+
 //! Set the type of a global jacobian matrix from the solution types
 SET_PROP(ImplicitBase, JacobianMatrix)
 {
diff --git a/dumux/implicit/subcontrolvolumeface.hh b/dumux/implicit/subcontrolvolumeface.hh
index ed2f9292eb67a4fbda62658fb7688de4222442d4..d9261e819a3da2488c4439afbbcf78f87d4f11a1 100644
--- a/dumux/implicit/subcontrolvolumeface.hh
+++ b/dumux/implicit/subcontrolvolumeface.hh
@@ -109,7 +109,7 @@ public:
     }
 
     //! index of the outside sub control volume for spatial param evaluation
-    // This results in undefined behaviour if boundary is false
+    // This results in undefined behaviour if boundary is true
     IndexType outsideScvIdx() const
     {
         return scvIndices_[1];
diff --git a/dumux/implicit/volumevariablesvector.hh b/dumux/implicit/volumevariablesvector.hh
index 6a3db2f03e4a01c1b7afff814451afb0441f157f..45d99cd44a8cae41bcda967b7a4dd335c5f75b78 100644
--- a/dumux/implicit/volumevariablesvector.hh
+++ b/dumux/implicit/volumevariablesvector.hh
@@ -32,20 +32,56 @@ namespace Dumux
  * \ingroup ImplicitModel
  * \brief Base class for the volume variables vector
  */
-template<class TypeTag>
-class VolumeVariablesVector : public std::vector<typename GET_PROP_TYPE(TypeTag, VolumeVariables)>
+template<class TypeTag, bool useOldSol, bool enableVolVarsCache>
+class VolumeVariablesVector
+{};
+
+// specialization in case of storing the volume variables
+template<class TypeTag, bool useOldSol>
+class VolumeVariablesVector<TypeTag, useOldSol,/*enableVolVarCaching*/true> : public std::vector<typename GET_PROP_TYPE(TypeTag, VolumeVariables)>
 {
+    friend VolumeVariablesVector<TypeTag, !useOldSol, true>;
+    friend ImplicitModel<TypeTag>;
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
     using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
     using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using IndexType = typename GridView::IndexSet::IndexType;
+    using Element = typename GridView::template Codim<0>::Entity;
+
+    enum{ isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+
+    VolumeVariablesVector<TypeTag, useOldSol, true>& operator= (const VolumeVariablesVector<TypeTag, useOldSol, true>& other) = default;
+
+    VolumeVariablesVector<TypeTag, useOldSol, true>& operator= (const VolumeVariablesVector<TypeTag, !useOldSol, true>& other)
+    {
+        // do the copy
+        numScvs_ = other.numScvs_;
+        numBoundaryScvs_ = other.numBoundaryScvs_;
+        volumeVariables_ = other.volumeVariables_;
+        boundaryVolumeVariables_ = other.boundaryVolumeVariables_;
+
+        // return the existing object
+        return *this;
+    };
 
 public:
-    void update(const Problem& problem, const SolutionVector& sol)
+    void update(Problem& problem, const SolutionVector& sol)
     {
+        numScvs_ = problem.model().fvGeometries().numScv();
+        if (!isBox)
+            numBoundaryScvs_ = problem.model().fvGeometries().numBoundaryScvf();
+        else
+            numBoundaryScvs_ = 0;
+
+        volumeVariables_.resize(numScvs_);
+        boundaryVolumeVariables_.resize(numBoundaryScvs_);
         for (const auto& element : elements(problem.gridView()))
         {
-            for (auto&& scv : problem.model().fvGeometries(element).scvs())
+            problem.model().fvGeometries_().bindElement(element);
+            const auto& fvGeometry = problem.model().fvGeometries(element);
+            for (const auto& scv : fvGeometry.scvs())
             {
                 (*this)[scv.index()].update(sol[scv.dofIndex()],
                                             problem,
@@ -53,8 +89,317 @@ public:
                                             scv);
             }
 
+            if (!isBox)
+            {
+                for (const auto& scvFace : fvGeometry.scvfs())
+                {
+                    // if we are not on a boundary, skip the rest
+                    if (!scvFace.boundary())
+                        continue;
+
+                    // When complex boundary handling is inactive, we only use BC vol vars on pure Dirichlet boundaries
+                    const auto bcTypes = problem.boundaryTypes(element, scvFace);
+                    if (/*TODO !GET_PROP_VALUE(TypeTag, BoundaryReconstruction) && */!(bcTypes.hasDirichlet() && !bcTypes.hasNeumann()))
+                        continue;
+
+                    const auto insideScvIdx = scvFace.insideScvIdx();
+                    const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx);
+                    const auto dirichletPriVars = problem.dirichlet(element, scvFace);
+
+                    (*this)[scvFace.outsideScvIdx()].update(dirichletPriVars, problem, element, insideScv);
+                }
+            }
+        }
+    }
+
+    const VolumeVariables& operator [](IndexType scvIdx) const
+    {
+        assert(scvIdx < numScvs_ + numBoundaryScvs_);
+
+        if (scvIdx < numScvs_)
+            return volumeVariables_[scvIdx];
+        else
+            return boundaryVolumeVariables_[scvIdx-numScvs_];
+    }
+
+    VolumeVariables& operator [](IndexType scvIdx)
+    {
+        assert(scvIdx < numScvs_ + numBoundaryScvs_);
+
+        if (scvIdx < numScvs_)
+            return volumeVariables_[scvIdx];
+        else
+            return boundaryVolumeVariables_[scvIdx-numScvs_];
+    }
+
+    // For compatibility reasons with the case of not storing the vol vars.
+    // function to be called before assembling an element, preparing the vol vars within the stencil
+    void bind(const Element& element) const {}
+    // function to prepare the vol vars within the element
+    void bindElement(const Element& element) const {}
+
+private:
+    IndexType numScvs_;
+    IndexType numBoundaryScvs_;
+    std::vector<VolumeVariables> volumeVariables_;
+    std::vector<VolumeVariables> boundaryVolumeVariables_;
+};
+
+
+// Specialization when the current volume variables are not stored
+template<class TypeTag>
+class VolumeVariablesVector<TypeTag, /*isOldSol*/false, /*enableVolVarCaching*/false>
+{
+    // prev vol vars have to be a friend class in order for the assignment operator to work
+    friend VolumeVariablesVector<TypeTag, true, false>;
+    friend ImplicitModel<TypeTag>;
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using IndexType = typename GridView::IndexSet::IndexType;
+    using Element = typename GridView::template Codim<0>::Entity;
+
+    VolumeVariablesVector& operator= (const VolumeVariablesVector<TypeTag, /*isOldSol*/false, /*enableVolVarCaching*/false>& other) = default;
+
+    void operator= (const VolumeVariablesVector<TypeTag, /*isOldSol*/true, /*enableVolVarCaching*/false>& other)
+    {
+        eIdxBound_ = -1;
+    }
+
+    VolumeVariablesVector() : problemPtr_(nullptr), eIdxBound_(-1) {}
+
+
+public:
+
+    void update(Problem& problem, const SolutionVector& sol)
+    {
+        problemPtr_ = &problem;
+    }
+
+
+    // Binding of an element, prepares the volume variables within the element stencil
+    // called by the local jacobian to prepare element assembly
+    // specialization for cc models
+    template <typename T = TypeTag>
+    typename std::enable_if<!GET_PROP_VALUE(TypeTag, ImplicitIsBox)>::type
+    bind(const Element& element)
+    {
+        const auto eIdx = problem_().elementMapper().index(element);
+        if (eIdx == eIdxBound_ && volVarIndices_.size() > 1)
+            return;
+
+        eIdxBound_ = eIdx;
+
+        // make sure the FVElementGeometry is bound to the element
+        problem_().model().fvGeometries_().bind(element);
+        const auto& fvGeometry = problem_().model().fvGeometries(eIdx);
+
+        // get stencil information
+        const auto& elementStencil = problem_().model().stencils(element).elementStencil();
+        const auto& neighborStencil = problem_().model().stencils(element).neighborStencil();
+        const auto numDofs = elementStencil.size();
+
+        // resize volume variables to the required size
+        volumeVariables_.resize(numDofs);
+        volVarIndices_.resize(numDofs);
+        int localIdx = 0;
+
+        // update the volume variables of the element at hand
+        const auto& scvI = problem_().model().fvGeometries().subControlVolume(eIdx);
+        const auto& solI = problem_().model().curSol()[eIdx];
+        volumeVariables_[localIdx].update(solI, problem_(), element, scvI);
+        volVarIndices_[localIdx] = scvI.index();
+        localIdx++;
+
+        // Update the volume variables of the neighboring elements and the boundary
+        for (auto globalJ : neighborStencil)
+        {
+            const auto& elementJ = problem_().model().fvGeometries().element(globalJ);
+            const auto& scvJ = problem_().model().fvGeometries().subControlVolume(globalJ);
+            const auto& solJ = problem_().model().curSol()[globalJ];
+            volumeVariables_[localIdx].update(solJ, problem_(), elementJ, scvJ);
+            volVarIndices_[localIdx] = scvJ.index();
+            localIdx++;
         }
+
+        // Check for boundary volume variables
+        for (const auto& scvFace : fvGeometry.scvfs())
+        {
+            // if we are not on a boundary, skip the rest
+            if (!scvFace.boundary())
+                continue;
+
+            // When complex boundary handling is inactive, we only use BC vol vars on pure Dirichlet boundaries
+            const auto bcTypes = problem_().boundaryTypes(element, scvFace);
+            if (/*TODO !GET_PROP_VALUE(TypeTag, BoundaryReconstruction) && */!(bcTypes.hasDirichlet() && !bcTypes.hasNeumann()))
+                continue;
+
+            volumeVariables_.resize(localIdx+1);
+            volVarIndices_.resize(localIdx+1);
+            const auto dirichletPriVars = problem_().dirichlet(element, scvFace);
+            volumeVariables_[localIdx].update(dirichletPriVars, problem_(), element, scvI);
+            volVarIndices_[localIdx] = scvFace.outsideScvIdx();
+            localIdx++;
+        }
+    }
+
+    // Binding of an element, prepares only the volume variables of the element
+    // specialization for cc models
+    template <typename T = TypeTag>
+    typename std::enable_if<!GET_PROP_VALUE(TypeTag, ImplicitIsBox)>::type
+    bindElement(const Element& element)
+    {
+        const auto eIdx = problem_().elementMapper().index(element);
+        if (eIdx == eIdxBound_ || std::find(volVarIndices_.begin(), volVarIndices_.end(), eIdx) != volVarIndices_.end())
+            return;
+
+        volumeVariables_.resize(1);
+        volVarIndices_.resize(1);
+        eIdxBound_ = eIdx;
+
+        // make sure the FVElementGeometry is bound to the element
+        problem_().model().fvGeometries_().bindElement(element);
+
+        // update the volume variables of the element at hand
+        const auto& scv = problem_().model().fvGeometries().subControlVolume(eIdx);
+        const auto& sol = problem_().model().curSol()[eIdx];
+        volumeVariables_[0].update(sol, problem_(), element, scv);
+        volVarIndices_[0] = scv.index();
+    }
+
+    const VolumeVariables& operator [](IndexType scvIdx) const
+    {
+        return volumeVariables_[getLocalIdx_(scvIdx)];
+    }
+
+    VolumeVariables& operator [](IndexType scvIdx)
+    {
+        return volumeVariables_[getLocalIdx_(scvIdx)];
+    }
+
+private:
+
+    void release_()
+    {
+        volumeVariables_.clear();
+        volVarIndices_.clear();
+    }
+
+    const int getLocalIdx_(const int volVarIdx) const
+    {
+        auto it = std::find(volVarIndices_.begin(), volVarIndices_.end(), volVarIdx);
+
+        if (it != volVarIndices_.end())
+            return std::distance(volVarIndices_.begin(), it);
+        else
+            DUNE_THROW(Dune::InvalidStateException, "Could not find the current volume variables for volVarIdx = " << volVarIdx <<
+                                                    ", make sure to properly bind the volume variables to the element before using them");
+    }
+
+    Problem& problem_() const
+    { return *problemPtr_;}
+
+    Problem* problemPtr_;
+    IndexType eIdxBound_;
+    std::vector<IndexType> volVarIndices_;
+    std::vector<VolumeVariables> volumeVariables_;
+};
+
+// Specialization when the previous volume variables are not stored
+template<class TypeTag>
+class VolumeVariablesVector<TypeTag, /*isOldSol*/true, /*enableVolVarCaching*/false>
+{
+    // current vol vars have to be a friend class in order for the assignment operator to work
+    friend VolumeVariablesVector<TypeTag, false, false>;
+    friend ImplicitModel<TypeTag>;
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using IndexType = typename GridView::IndexSet::IndexType;
+    using Element = typename GridView::template Codim<0>::Entity;
+
+    VolumeVariablesVector& operator= (const VolumeVariablesVector<TypeTag, /*isOldSol*/false, /*enableVolVarCaching*/false>& other)
+    {
+        release_();
+        problemPtr_ = other.problemPtr_;
+        eIdxBound_ = -1;
+        return *this;
+    }
+
+    VolumeVariablesVector() : problemPtr_(nullptr), eIdxBound_(-1) {}
+
+
+public:
+
+    void update(const Problem& problem, const SolutionVector& sol)
+    {
+        problemPtr_ = &problem;
+    }
+
+    // Binding of an element, prepares the volume variables of only the element
+    // specialization for cc models
+    template <typename T = TypeTag>
+    typename std::enable_if<!GET_PROP_VALUE(TypeTag, ImplicitIsBox)>::type
+    bindElement(const Element& element)
+    {
+        const auto eIdx = problem_().elementMapper().index(element);
+        if (eIdx == eIdxBound_)
+            return;
+
+        volumeVariables_.resize(1);
+        volVarIndices_.resize(1);
+        eIdxBound_ = eIdx;
+
+        // make sure FVElementGeometry is bound to the element
+        problem_().model().fvGeometries_().bindElement(element);
+
+        // update the volume variables of the element at hand
+        const auto& scv = problem_().model().fvGeometries().subControlVolume(eIdx);
+        const auto& sol = problem_().model().prevSol()[eIdx];
+        volumeVariables_[0].update(sol, problem_(), element, scv);
+        volVarIndices_[0] = scv.index();
     }
+
+    const VolumeVariables& operator [](IndexType scvIdx) const
+    {
+        return volumeVariables_[getLocalIdx_(scvIdx)];
+    }
+
+    VolumeVariables& operator [](IndexType scvIdx)
+    {
+        return volumeVariables_[getLocalIdx_(scvIdx)];
+    }
+
+private:
+
+    void release_()
+    {
+        volumeVariables_.clear();
+        volVarIndices_.clear();
+    }
+
+    const int getLocalIdx_(const int volVarIdx) const
+    {
+        auto it = std::find(volVarIndices_.begin(), volVarIndices_.end(), volVarIdx);
+
+        if (it != volVarIndices_.end())
+            return std::distance(volVarIndices_.begin(), it);
+        else
+            DUNE_THROW(Dune::InvalidStateException, "Could not find the previous volume variables for volVarIdx = " << volVarIdx <<
+                                                    ", make sure to properly bind the volume variables to the element before using them");
+    }
+
+    Problem& problem_() const
+    { return *problemPtr_;}
+
+    Problem* problemPtr_;
+    IndexType eIdxBound_;
+    std::vector<IndexType> volVarIndices_;
+    std::vector<VolumeVariables> volumeVariables_;
 };
 
 } // end namespace
diff --git a/dumux/porousmediumflow/1p/implicit/model.hh b/dumux/porousmediumflow/1p/implicit/model.hh
index 50c585d679ba999ecaa9febda235ff64b95860e1..6b6684a34fd512565fb42e58a384bb82ac41f929 100644
--- a/dumux/porousmediumflow/1p/implicit/model.hh
+++ b/dumux/porousmediumflow/1p/implicit/model.hh
@@ -107,7 +107,12 @@ public:
                 int eIdx = this->elementMapper().index(element);
                 (*rank)[eIdx] = this->gridView_().comm().rank();
 
-                for (auto&& scv : this->fvGeometries(element).scvs())
+                // make sure FVElementGeometry and the volume variables are bound to the element
+                this->fvGeometries_().bindElement(element);
+                this->curVolVars_().bindElement(element);
+
+                const auto& fvGeometry = this->fvGeometries(element);
+                for (const auto& scv : fvGeometry.scvs())
                 {
                     const auto& spatialParams = this->problem_().spatialParams();
                     auto dofIdxGlobal = scv.dofIndex();
diff --git a/dumux/porousmediumflow/2p/implicit/model.hh b/dumux/porousmediumflow/2p/implicit/model.hh
index d0ad7eebdde6bba2cfdf2b42f9b17f4e3f42c700..4014bece2f3c33e73d39261fb4b399c6ed7add3a 100644
--- a/dumux/porousmediumflow/2p/implicit/model.hh
+++ b/dumux/porousmediumflow/2p/implicit/model.hh
@@ -144,7 +144,12 @@ public:
                 int eIdx = this->elementMapper().index(element);
                 (*rank)[eIdx] = this->gridView_().comm().rank();
 
-                for (auto&& scv : this->fvGeometries(element).scvs())
+                // make sure FVElementGeometry & vol vars are bound to the element
+                this->fvGeometries_().bindElement(element);
+                this->curVolVars_().bindElement(element);
+
+                const auto& fvGeometry = this->fvGeometries(element);
+                for (const auto& scv : fvGeometry.scvs())
                 {
                     auto dofIdxGlobal = scv.dofIndex();
                     const auto& volVars = this->curVolVars(scv);
diff --git a/dumux/porousmediumflow/3p/implicit/model.hh b/dumux/porousmediumflow/3p/implicit/model.hh
index 0048e07f59366c2b71de4ad0f488f2e221f7c27f..697e796f03371988ae21e1671ef62103c26e4ac7 100644
--- a/dumux/porousmediumflow/3p/implicit/model.hh
+++ b/dumux/porousmediumflow/3p/implicit/model.hh
@@ -139,9 +139,12 @@ public:
         {
             if(element.partitionType() == Dune::InteriorEntity)
             {
-                const auto& fvGeometry = this->fvGeometries(element);
+                // make sure FVElementGeometry & vol vars are bound to the element
+                this->fvGeometries_().bindElement(element);
+                this->curVolVars_().bindElement(element);
 
-                for (auto&& scv : fvGeometry.scvs())
+                const auto& fvGeometry = this->fvGeometries(element);
+                for (const auto& scv : fvGeometry.scvs())
                 {
                     auto eIdx = scv.elementIndex();
                     auto dofIdxGlobal = scv.dofIndex();
diff --git a/dumux/porousmediumflow/3p3c/implicit/localresidual.hh b/dumux/porousmediumflow/3p3c/implicit/localresidual.hh
index ffe92cbabc18dfd923c54bb094aaf5beb714acd4..9f6c94b9c208d590cf3497fd0379a8386d965e47 100644
--- a/dumux/porousmediumflow/3p3c/implicit/localresidual.hh
+++ b/dumux/porousmediumflow/3p3c/implicit/localresidual.hh
@@ -47,6 +47,7 @@ class ThreePThreeCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual
     typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
     typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
     typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
     typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
 
     enum {
@@ -117,8 +118,8 @@ public:
      */
     PrimaryVariables computeFlux(const SubControlVolumeFace& scvFace)
     {
-        auto& fluxVars = this->model_().fluxVars_(scvFace);
-        fluxVars.beginFluxComputation();
+        FluxVariables fluxVars;
+        fluxVars.initAndComputeFluxes(asImp_()->problem_(), this->element_(), scvFace);
 
         // get upwind weights into local scope
         auto massWeight = massWeight_;
@@ -141,17 +142,17 @@ public:
 
                 // get equation index
                 auto eqIdx = conti0EqIdx + compIdx;
-                flux[eqIdx] += fluxVars.advection().flux(phaseIdx, upwindRule);
+                flux[eqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindRule);
             }
         }
 
         // diffusive fluxes
-        Scalar jGW = fluxVars.molecularDiffusion(wPhaseIdx, gCompIdx).flux();
-        Scalar jNW = fluxVars.molecularDiffusion(wPhaseIdx, nCompIdx).flux();
+        Scalar jGW = fluxVars.molecularDiffusionFlux(wPhaseIdx, gCompIdx);
+        Scalar jNW = fluxVars.molecularDiffusionFlux(wPhaseIdx, nCompIdx);
         Scalar jWW = -(jGW+jNW);
 
-        Scalar jWG = fluxVars.molecularDiffusion(gPhaseIdx, wCompIdx).flux();
-        Scalar jNG = fluxVars.molecularDiffusion(gPhaseIdx, nCompIdx).flux();
+        Scalar jWG = fluxVars.molecularDiffusionFlux(gPhaseIdx, wCompIdx);
+        Scalar jNG = fluxVars.molecularDiffusionFlux(gPhaseIdx, nCompIdx);
         Scalar jGG = -(jWG+jNG);
 
         // Scalar jWN = fluxVars.molecularDiffusion().flux(nPhaseIdx, wCompIdx);
diff --git a/dumux/porousmediumflow/3p3c/implicit/model.hh b/dumux/porousmediumflow/3p3c/implicit/model.hh
index 16fbf7bfec32a89918e7cbf1a93f19f401bceabf..5fdaea88fcaf7164984f4eb0c5668e9dba1f0087 100644
--- a/dumux/porousmediumflow/3p3c/implicit/model.hh
+++ b/dumux/porousmediumflow/3p3c/implicit/model.hh
@@ -159,7 +159,12 @@ public:
         {
             for (const auto& element : elements(this->problem_().gridView()))
             {
-                for (auto&& scv : this->fvGeometries(element).scvs())
+                // make sure FVElementGeometry & vol vars are bound to the element
+                this->fvGeometries_().bindElement(element);
+                this->curVolVars_().bindElement(element);
+
+                const auto& fvGeometry = this->fvGeometries(element);
+                for (const auto& scv : fvGeometry.scvs())
                 {
                     auto dofIdxGlobal = scv.dofIndex();
                     if (priVarSwitch_().wasSwitched(dofIdxGlobal))
@@ -274,9 +279,12 @@ public:
                 int eIdx = this->problem_().elementMapper().index(element);
                 (*rank)[eIdx] = this->gridView_().comm().rank();
 
-                const auto& fvGeometry = this->fvGeometries(element);
+                // make sure FVElementGeometry & vol vars are bound to the element
+                this->fvGeometries_().bindElement(element);
+                this->curVolVars_().bindElement(element);
 
-                for (auto&& scv : fvGeometry.scvs())
+                const auto& fvGeometry = this->fvGeometries(element);
+                for (const auto& scv : fvGeometry.scvs())
                 {
                     const auto& volVars = this->curVolVars(scv);
                     int dofIdxGlobal = scv.dofIndex();
diff --git a/dumux/porousmediumflow/compositional/primaryvariableswitch.hh b/dumux/porousmediumflow/compositional/primaryvariableswitch.hh
index 799a2a5ecfff8f68725919852745e17b8df49214..1eae984eeb1bf9df68ba44626660a10025462c51 100644
--- a/dumux/porousmediumflow/compositional/primaryvariableswitch.hh
+++ b/dumux/porousmediumflow/compositional/primaryvariableswitch.hh
@@ -45,11 +45,11 @@ class PrimaryVariableSwitch
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
     using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
-    using VolumeVariablesVector = typename GET_PROP_TYPE(TypeTag, VolumeVariablesVector);
+    using CurrentVolumeVariablesVector = typename GET_PROP_TYPE(TypeTag, CurrentVolumeVariablesVector);
 
 public:
 
-    void init(const Problem& problem)
+    void init(Problem& problem)
     {
         auto numDofs = problem.model().numDofs();
         phasePresence_.resize(numDofs);
@@ -57,6 +57,9 @@ public:
 
         for (const auto& element : elements(problem.gridView()))
         {
+            // make sure FVElementGeometry is bound to the element
+            problem.model().fvGeometries_().bindElement(element);
+
             const auto& fvGeometry = problem.model().fvGeometries(element);
             for (auto&& scv : fvGeometry.scvs())
             {
@@ -113,15 +116,19 @@ public:
      * \param problem The problem
      * \param curSol The current solution to be updated / modified
      */
-    bool update(const Problem& problem, SolutionVector& curSol,
-                VolumeVariablesVector& volVarsVector)
+    bool update(Problem& problem, SolutionVector& curSol,
+                CurrentVolumeVariablesVector& volVarsVector)
     {
         bool switched = false;
         visited_.assign(phasePresence_.size(), false);
         for (const auto& element : elements(problem.gridView()))
         {
+            // make sure FVElementGeometry is bound to the element
+            problem.model().fvGeometries_().bindElement(element);
+            volVarsVector.bindElement(element);
+
             const auto& fvGeometry = problem.model().fvGeometries(element);
-            for (auto&& scv : fvGeometry.scvs())
+            for (const auto& scv : fvGeometry.scvs())
             {
                 auto dofIdxGlobal = scv.dofIndex();
                 if (!visited_[dofIdxGlobal])
diff --git a/dumux/porousmediumflow/immiscible/localresidual.hh b/dumux/porousmediumflow/immiscible/localresidual.hh
index 2fe3347df61b6722f2ddccf06f2cc027d3379c01..ee999d5aaedbaf090add9ff343cde9738be2a992 100644
--- a/dumux/porousmediumflow/immiscible/localresidual.hh
+++ b/dumux/porousmediumflow/immiscible/localresidual.hh
@@ -97,8 +97,8 @@ public:
      */
     PrimaryVariables computeFlux(const SubControlVolumeFace& scvf)
     {
-        auto& fluxVars = this->model_().fluxVars_(scvf);
-        fluxVars.beginFluxComputation();
+        FluxVariables fluxVars;
+        fluxVars.initAndComputeFluxes(asImp_()->problem_(), this->element_(), scvf);
 
         // copy weight to local scope for use in lambda expression
         auto w = upwindWeight_;
@@ -112,7 +112,7 @@ public:
                                      + (dn.density(phaseIdx)*dn.mobility(phaseIdx))*(1-w); };
 
             auto eqIdx = conti0EqIdx + phaseIdx;
-            flux[eqIdx] = fluxVars.advection().flux(phaseIdx, upwindRule);
+            flux[eqIdx] = fluxVars.advectiveFlux(phaseIdx, upwindRule);
         }
 
         return flux;
diff --git a/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh b/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh
index 41f8016f2c92c7a0631b2ec34c59a5a3e284a5ff..49a35b14d0539d1b0bb67b9f312f5d0c81e1ce36 100644
--- a/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh
+++ b/dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyslaw.hh
@@ -53,16 +53,16 @@ NEW_PROP_TAG(ProblemEnableGravity);
 template <class TypeTag>
 class CCTpfaDarcysLaw
 {
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
     typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
     typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
     typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
-    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
     typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
     typedef typename GridView::IndexSet::IndexType IndexType;
-    using Element = typename GridView::template Codim<0>::Entity;
+    typedef std::vector<IndexType> Stencil;
 
+    using Element = typename GridView::template Codim<0>::Entity;
     enum { dim = GridView::dimension} ;
     enum { dimWorld = GridView::dimensionworld} ;
     enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
@@ -72,234 +72,137 @@ class CCTpfaDarcysLaw
 
 public:
 
-    void update(const Problem &problem,
-                const Element& element,
-                const SubControlVolumeFace &scvFace)
+    static Scalar flux(const Problem& problem,
+                       const SubControlVolumeFace& scvFace,
+                       const IndexType phaseIdx)
     {
-        problemPtr_ = &problem;
-        scvFacePtr_ = &scvFace;
-        enableGravity_ = GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity);
-
-        updateTransmissibilities_();
-        updateStencil_();
-    }
-
-    void update(const Problem& problem,
-                const Element& element,
-                const SubControlVolumeFace &scvFace,
-                VolumeVariables* boundaryVolVars)
-    {
-        boundaryVolVars_ = boundaryVolVars;
-        update(problem, element, scvFace);
-    }
-
-    void updateTransmissibilities(const Problem &problem, const SubControlVolumeFace &scvFace)
-    {
-        updateTransmissibilities_();
-    }
+        const auto& tij = getTransmissibilities(problem, scvFace);
 
-    void beginFluxComputation(bool boundaryVolVarsUpdated = false)
-    {
         // Get the inside volume variables
-        const auto insideScvIdx = scvFace_().insideScvIdx();
-        const auto& insideScv = problem_().model().fvGeometries().subControlVolume(insideScvIdx);
-        const auto* insideVolVars = &problem_().model().curVolVars(insideScv);
+        const auto insideScvIdx = scvFace.insideScvIdx();
+        const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx);
+        const auto& insideVolVars = problem.model().curVolVars(insideScv);
 
         // and the outside volume variables
-        const VolumeVariables* outsideVolVars;
-        if (!scvFace_().boundary())
-            outsideVolVars = &problem_().model().curVolVars(scvFace_().outsideScvIdx());
-        else
-        {
-            outsideVolVars = boundaryVolVars_;
-            if (!boundaryVolVarsUpdated)
-            {
-                // update the boudary volvars for Dirichlet boundaries
-                const auto element = problem_().model().fvGeometries().element(insideScv);
-                const auto dirichletPriVars = problem_().dirichlet(element, scvFace_());
-                boundaryVolVars_->update(dirichletPriVars, problem_(), element, insideScv);
-            }
-        }
+        const auto& outsideVolVars = problem.model().curVolVars(scvFace.outsideScvIdx());
 
-        // loop over all phases to compute the volume flux
-        for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
+        auto hInside = insideVolVars.pressure(phaseIdx);
+        auto hOutside = outsideVolVars.pressure(phaseIdx);
+
+        if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
         {
-            auto hInside = insideVolVars->pressure(phaseIdx);
-            auto hOutside = outsideVolVars->pressure(phaseIdx);
+            // do averaging for the density
+            const auto rhoInside = insideVolVars.density(phaseIdx);
+            const auto rhoOutide = outsideVolVars.density(phaseIdx);
+            const auto rho = (rhoInside + rhoOutide)*0.5;
 
-            if (enableGravity_)
-            {
-                // do averaging for the density
-                const auto rhoInside = insideVolVars->density(phaseIdx);
-                const auto rhoOutide = outsideVolVars->density(phaseIdx);
-                const auto rho = (rhoInside + rhoOutide)*0.5;
-
-
-                // ask for the gravitational acceleration in the inside neighbor
-                const auto xInside = insideScv.center();
-                const auto gInside = problem_().gravityAtPos(xInside);
-
-                hInside -= rho*(gInside*xInside);
-
-                // and the outside neighbor
-                if (scvFace_().boundary())
-                {
-                    const auto xOutside = scvFace_().center();
-                    const auto gOutside = problem_().gravityAtPos(xOutside);
-                    hOutside -= rho*(gOutside*xOutside);
-                }
-                else
-                {
-                    const auto outsideScvIdx = scvFace_().outsideScvIdx();
-                    const auto& outsideScv = problem_().model().fvGeometries().subControlVolume(outsideScvIdx);
-                    const auto xOutside = outsideScv.center();
-                    const auto gOutside = problem_().gravityAtPos(xOutside);
-                    hOutside -= rho*(gOutside*xOutside);
-                }
-            }
 
-            kGradPNormal_[phaseIdx] = tij_*(hInside - hOutside);
+            // ask for the gravitational acceleration in the inside neighbor
+            const auto xInside = insideScv.center();
+            const auto gInside = problem.gravityAtPos(xInside);
 
-            if (std::signbit(kGradPNormal_[phaseIdx]))
+            hInside -= rho*(gInside*xInside);
+
+            // and the outside neighbor
+            if (scvFace.boundary())
             {
-                if (scvFace_().boundary())
-                    upWindIndices_[phaseIdx] = std::make_pair(-1, scvFace_().insideScvIdx());
-                else
-                    upWindIndices_[phaseIdx] = std::make_pair(scvFace_().outsideScvIdx(), scvFace_().insideScvIdx());
+                const auto xOutside = scvFace.center();
+                const auto gOutside = problem.gravityAtPos(xOutside);
+                hOutside -= rho*(gOutside*xOutside);
             }
             else
             {
-                if (scvFace_().boundary())
-                    upWindIndices_[phaseIdx] = std::make_pair(scvFace_().insideScvIdx(), -1);
-                else
-                    upWindIndices_[phaseIdx] = std::make_pair(scvFace_().insideScvIdx(), scvFace_().outsideScvIdx());
+                const auto outsideScvIdx = scvFace.outsideScvIdx();
+                const auto& outsideScv = problem.model().fvGeometries().subControlVolume(outsideScvIdx);
+                const auto xOutside = outsideScv.center();
+                const auto gOutside = problem.gravityAtPos(xOutside);
+                hOutside -= rho*(gOutside*xOutside);
             }
         }
+
+        return tij*(hInside - hOutside);
     }
 
-    /*!
-     * \brief A function to calculate the mass flux over a sub control volume face
-     *
-     * \param phaseIdx The index of the phase of which the flux is to be calculated
-     * \param upwindFunction A function which does the upwinding
-     */
-    template<typename FunctionType>
-    Scalar flux(IndexType phaseIdx, FunctionType upwindFunction) const
+    static Stencil stencil(const Problem& problem, const SubControlVolumeFace& scvFace)
     {
-        return kGradPNormal_[phaseIdx]*upwindFunction(upVolVars(phaseIdx), dnVolVars(phaseIdx));
+        std::vector<IndexType> stencil;
+        if (!scvFace.boundary())
+        {
+            stencil.push_back(scvFace.insideScvIdx());
+            stencil.push_back(scvFace.outsideScvIdx());
+        }
+        else
+            stencil.push_back(scvFace.insideScvIdx());
+
+        return stencil;
     }
 
-    const std::set<IndexType>& stencil() const
+    template <typename T = TypeTag>
+    static const typename std::enable_if<GET_PROP_VALUE(T, EnableFluxVariablesCache), Scalar>::type& getTransmissibilities(const Problem& problem, const SubControlVolumeFace& scvFace)
     {
-        return stencil_;
+        return problem.model().fluxVarsCache(scvFace).tij();
     }
 
-    const VolumeVariables& upVolVars(IndexType phaseIdx) const
+    template <typename T = TypeTag>
+    static const typename std::enable_if<!GET_PROP_VALUE(T, EnableFluxVariablesCache), Scalar>::type getTransmissibilities(const Problem& problem, const SubControlVolumeFace& scvFace)
     {
-        if(upWindIndices_[phaseIdx].first != -1)
-            return problem_().model().curVolVars(upWindIndices_[phaseIdx].first);
-        else
-            return *boundaryVolVars_;
+        return calculateTransmissibilities(problem, scvFace);
     }
 
-    const VolumeVariables& dnVolVars(IndexType phaseIdx) const
+    static Scalar calculateTransmissibilities(const Problem& problem, const SubControlVolumeFace& scvFace)
     {
-        if(upWindIndices_[phaseIdx].second != -1)
-            return problem_().model().curVolVars(upWindIndices_[phaseIdx].second);
-        else
-            return *boundaryVolVars_;
-    }
+        Scalar tij;
 
-private:
+        const auto insideScvIdx = scvFace.insideScvIdx();
+        const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx);
+        const auto insideK = problem.spatialParams().intrinsicPermeability(insideScv);
+        Scalar ti = calculateOmega_(problem, scvFace, insideK, insideScv);
 
-    void updateTransmissibilities_()
-    {
-        if (!scvFace_().boundary())
+        if (!scvFace.boundary())
         {
-            const auto insideScvIdx = scvFace_().insideScvIdx();
-            const auto& insideScv = problem_().model().fvGeometries().subControlVolume(insideScvIdx);
-            const auto insideK = problem_().spatialParams().intrinsicPermeability(insideScv);
-            Scalar ti = calculateOmega_(insideK, insideScv);
-
-            const auto outsideScvIdx = scvFace_().outsideScvIdx();
-            const auto& outsideScv = problem_().model().fvGeometries().subControlVolume(outsideScvIdx);
-            const auto outsideK = problem_().spatialParams().intrinsicPermeability(outsideScv);
-            Scalar tj = -1.0*calculateOmega_(outsideK, outsideScv);
+            const auto outsideScvIdx = scvFace.outsideScvIdx();
+            const auto& outsideScv = problem.model().fvGeometries().subControlVolume(outsideScvIdx);
+            const auto outsideK = problem.spatialParams().intrinsicPermeability(outsideScv);
+            Scalar tj = -1.0*calculateOmega_(problem, scvFace, outsideK, outsideScv);
 
-            tij_ = scvFace_().area()*(ti * tj)/(ti + tj);
+            tij = scvFace.area()*(ti * tj)/(ti + tj);
         }
         else
         {
-            const auto insideScvIdx = scvFace_().insideScvIdx();
-            const auto& insideScv = problem_().model().fvGeometries().subControlVolume(insideScvIdx);
-            const auto insideK = problem_().spatialParams().intrinsicPermeability(insideScv);
-            Scalar ti = calculateOmega_(insideK, insideScv);
-
-            tij_ = scvFace_().area()*ti;
+            tij = scvFace.area()*ti;
         }
-    }
 
-    void updateStencil_()
-    {
-        // fill the stencil
-        if (!scvFace_().boundary())
-            stencil_.insert({scvFace_().insideScvIdx(), scvFace_().outsideScvIdx()});
-        else
-            // fill the stencil
-            stencil_.insert(scvFace_().insideScvIdx());
+        return tij;
     }
 
-    Scalar calculateOmega_(const DimWorldMatrix &K, const SubControlVolume &scv) const
+private:
+
+    static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, const DimWorldMatrix &K, const SubControlVolume &scv)
     {
         GlobalPosition Knormal;
-        K.mv(scvFace_().unitOuterNormal(), Knormal);
+        K.mv(scvFace.unitOuterNormal(), Knormal);
 
-        auto distanceVector = scvFace_().center();
+        auto distanceVector = scvFace.center();
         distanceVector -= scv.center();
         distanceVector /= distanceVector.two_norm2();
 
         Scalar omega = Knormal * distanceVector;
-        omega *= problem_().model().curVolVars(scv).extrusionFactor();
+        omega *= problem.model().curVolVars(scv).extrusionFactor();
 
         return omega;
     }
 
-    Scalar calculateOmega_(Scalar K, const SubControlVolume &scv) const
+    static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, Scalar K, const SubControlVolume &scv)
     {
-        auto distanceVector = scvFace_().center();
+        auto distanceVector = scvFace.center();
         distanceVector -= scv.center();
         distanceVector /= distanceVector.two_norm2();
 
-        Scalar omega = K * (distanceVector * scvFace_().unitOuterNormal());
-        omega *= problem_().model().curVolVars(scv).extrusionFactor();
+        Scalar omega = K * (distanceVector * scvFace.unitOuterNormal());
+        omega *= problem.model().curVolVars(scv).extrusionFactor();
 
         return omega;
     }
-
-    const Problem &problem_() const
-    {
-        return *problemPtr_;
-    }
-
-    const SubControlVolumeFace& scvFace_() const
-    {
-        return *scvFacePtr_;
-    }
-
-    const Problem *problemPtr_;              //! Pointer to the problem
-    const SubControlVolumeFace *scvFacePtr_; //! Pointer to the sub control volume face for which the flux variables are created
-    bool enableGravity_;                     //! If we have a problem considering gravitational effects
-    std::set<IndexType> stencil_;         //! Indices of the cells of which the pressure is needed for the flux calculation
-
-    //! Boundary volume variables (they only get updated on Dirichlet boundaries)
-    VolumeVariables* boundaryVolVars_;
-
-    //! The upstream (first) and downstream (second) volume variable indices
-    std::array<std::pair<IndexType, IndexType>, numPhases> upWindIndices_;
-
-    //! Precomputed values
-    Scalar tij_; //! transmissibility for the flux calculation tij(ui - uj)
-    std::array<Scalar, numPhases> kGradPNormal_; //! K(grad(p) - rho*g)*n
 };
 
 } // end namespace
diff --git a/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh b/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh
index 0636512688bccdf1f8a3a286fb64a29bd9232f2c..309e438e81f99f624b0eae5ae033e3e3c7d5253b 100644
--- a/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh
+++ b/dumux/porousmediumflow/implicit/cellcentered/tpfa/fickslaw.hh
@@ -56,10 +56,10 @@ class CCTpfaFicksLaw
     typedef typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel) EffDiffModel;
     typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
     typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
-    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
     typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
     typedef typename GridView::IndexSet::IndexType IndexType;
+    typedef typename std::vector<IndexType> Stencil;
+
     using Element = typename GridView::template Codim<0>::Entity;
 
     enum { dim = GridView::dimension} ;
@@ -71,174 +71,106 @@ class CCTpfaFicksLaw
 
 public:
 
-    void update(const Problem& problem,
-                const Element& element,
-                const SubControlVolumeFace &scvFace,
-                int phaseIdx, int compIdx)
-    {
-        problemPtr_ = &problem;
-        scvFacePtr_ = &scvFace;
-
-        phaseIdx_ = phaseIdx;
-        compIdx_ = compIdx;
-
-        updateStencil_();
-
-        // TODO for non solution dependent diffusion tensors...
-    }
-
-    void update(const Problem &problem,
-                const Element &element,
-                const SubControlVolumeFace &scvFace,
-                int phaseIdx, int compIdx,
-                VolumeVariables* boundaryVolVars)
-    {
-        boundaryVolVars_ = boundaryVolVars;
-        update(problem, element, scvFace, phaseIdx, compIdx);
-    }
-
-    void beginFluxComputation(bool boundaryVolVarsUpdated = false)
+    static Scalar flux(const Problem& problem,
+                       const SubControlVolumeFace& scvFace,
+                       const int phaseIdx,
+                       const int compIdx)
     {
         // diffusion tensors are always solution dependent
-        updateTransmissibilities_();
+        Scalar tij = calculateTransmissibility_(problem, scvFace, phaseIdx, compIdx);
 
         // Get the inside volume variables
-        const auto insideScvIdx = scvFace_().insideScvIdx();
-        const auto& insideScv = problem_().model().fvGeometries().subControlVolume(insideScvIdx);
-        const auto* insideVolVars = &problem_().model().curVolVars(insideScv);
+        const auto insideScvIdx = scvFace.insideScvIdx();
+        const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx);
+        const auto& insideVolVars = problem.model().curVolVars(insideScv);
 
         // and the outside volume variables
-        const VolumeVariables* outsideVolVars;
-        if (!scvFace_().boundary())
-            outsideVolVars = &problem_().model().curVolVars(scvFace_().outsideScvIdx());
-        else
-        {
-            outsideVolVars = boundaryVolVars_;
-            if (!boundaryVolVarsUpdated)
-            {
-                // update the boudary volvars for Dirichlet boundaries
-                const auto element = problem_().model().fvGeometries().element(insideScv);
-                const auto dirichletPriVars = problem_().dirichlet(element, scvFace_());
-                boundaryVolVars_->update(dirichletPriVars, problem_(), element, insideScv);
-            }
-        }
+        const auto& outsideVolVars = problem.model().curVolVars(scvFace.outsideScvIdx());
 
         // compute the diffusive flux
-        const auto xInside = insideVolVars->moleFraction(phaseIdx_, compIdx_);
-        const auto xOutside = outsideVolVars->moleFraction(phaseIdx_, compIdx_);
-        const auto rho = 0.5*(insideVolVars->molarDensity(phaseIdx_) + outsideVolVars->molarDensity(phaseIdx_));
+        const auto xInside = insideVolVars.moleFraction(phaseIdx, compIdx);
+        const auto xOutside = outsideVolVars.moleFraction(phaseIdx, compIdx);
+        const auto rho = 0.5*(insideVolVars.molarDensity(phaseIdx) + outsideVolVars.molarDensity(phaseIdx));
 
-        rhoDGradXNormal_ = rho*tij_*(xInside - xOutside);
+        return rho*tij*(xInside - xOutside);
     }
 
-
-
-    /*!
-     * \brief A function to calculate the mass flux over a sub control volume face
-     *
-     * \param phaseIdx The index of the phase of which the flux is to be calculated
-     * \param compIdx The index of the transported component
-     */
-    Scalar flux() const
+    static Stencil stencil(const Problem& problem, const SubControlVolumeFace& scvFace)
     {
-        return rhoDGradXNormal_;
-    }
+        std::vector<IndexType> stencil;
+        stencil.clear();
+        if (!scvFace.boundary())
+        {
+            stencil.push_back(scvFace.insideScvIdx());
+            stencil.push_back(scvFace.outsideScvIdx());
+        }
+        else
+            stencil.push_back(scvFace.insideScvIdx());
 
-    std::set<IndexType> stencil() const
-    {
-        return stencil_;
+        return stencil;
     }
 
-protected:
+private:
 
 
-    void updateTransmissibilities_()
+    static Scalar calculateTransmissibility_(const Problem& problem, const SubControlVolumeFace& scvFace, const int phaseIdx, const int compIdx)
     {
-        const auto insideScvIdx = scvFace_().insideScvIdx();
-        const auto& insideScv = problem_().model().fvGeometries().subControlVolume(insideScvIdx);
-        const auto& insideVolVars = problem_().model().curVolVars(insideScvIdx);
+        Scalar tij;
+
+        const auto insideScvIdx = scvFace.insideScvIdx();
+        const auto& insideScv = problem.model().fvGeometries().subControlVolume(insideScvIdx);
+        const auto& insideVolVars = problem.model().curVolVars(insideScvIdx);
 
-        auto insideD = insideVolVars.diffusionCoefficient(phaseIdx_, compIdx_);
-        insideD = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx_), insideD);
-        Scalar ti = calculateOmega_(insideD, insideScv);
+        auto insideD = insideVolVars.diffusionCoefficient(phaseIdx, compIdx);
+        insideD = EffDiffModel::effectiveDiffusivity(insideVolVars.porosity(), insideVolVars.saturation(phaseIdx), insideD);
+        Scalar ti = calculateOmega_(problem, scvFace, insideD, insideScv);
 
-        if (!scvFace_().boundary())
+        if (!scvFace.boundary())
         {
-            const auto outsideScvIdx = scvFace_().outsideScvIdx();
-            const auto& outsideScv = problem_().model().fvGeometries().subControlVolume(outsideScvIdx);
-            const auto& outsideVolVars = problem_().model().curVolVars(outsideScvIdx);
+            const auto outsideScvIdx = scvFace.outsideScvIdx();
+            const auto& outsideScv = problem.model().fvGeometries().subControlVolume(outsideScvIdx);
+            const auto& outsideVolVars = problem.model().curVolVars(outsideScvIdx);
 
-            auto outsideD = outsideVolVars.diffusionCoefficient(phaseIdx_, compIdx_);
-            outsideD = EffDiffModel::effectiveDiffusivity(outsideVolVars.porosity(), outsideVolVars.saturation(phaseIdx_), outsideD);
-            Scalar tj = -1.0*calculateOmega_(outsideD, outsideScv);
+            auto outsideD = outsideVolVars.diffusionCoefficient(phaseIdx, compIdx);
+            outsideD = EffDiffModel::effectiveDiffusivity(outsideVolVars.porosity(), outsideVolVars.saturation(phaseIdx), outsideD);
+            Scalar tj = -1.0*calculateOmega_(problem, scvFace, outsideD, outsideScv);
 
-            tij_ = scvFace_().area()*(ti * tj)/(ti + tj);
+            tij = scvFace.area()*(ti * tj)/(ti + tj);
         }
         else
         {
-            tij_ = scvFace_().area()*ti;
+            tij = scvFace.area()*ti;
         }
-    }
 
-    void updateStencil_()
-    {
-        // fill the stencil
-        if (!scvFace_().boundary())
-            stencil_= {scvFace_().insideScvIdx(), scvFace_().outsideScvIdx()};
-        else
-            // fill the stencil
-            stencil_ = {scvFace_().insideScvIdx()};
+        return tij;
     }
 
-    Scalar calculateOmega_(const DimWorldMatrix &D, const SubControlVolume &scv) const
+    static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, const DimWorldMatrix &D, const SubControlVolume &scv)
     {
         GlobalPosition Dnormal;
-        D.mv(scvFace_().unitOuterNormal(), Dnormal);
+        D.mv(scvFace.unitOuterNormal(), Dnormal);
 
-        auto distanceVector = scvFace_().center();
+        auto distanceVector = scvFace.center();
         distanceVector -= scv.center();
         distanceVector /= distanceVector.two_norm2();
 
         Scalar omega = Dnormal * distanceVector;
-        omega *= problem_().model().curVolVars(scv).extrusionFactor();
+        omega *= problem.model().curVolVars(scv).extrusionFactor();
 
         return omega;
     }
 
-    Scalar calculateOmega_(Scalar D, const SubControlVolume &scv) const
+    static Scalar calculateOmega_(const Problem& problem, const SubControlVolumeFace& scvFace, Scalar D, const SubControlVolume &scv)
     {
-        auto distanceVector = scvFace_().center();
+        auto distanceVector = scvFace.center();
         distanceVector -= scv.center();
         distanceVector /= distanceVector.two_norm2();
 
-        Scalar omega = D * (distanceVector * scvFace_().unitOuterNormal());
-        omega *= problem_().model().curVolVars(scv).extrusionFactor();
+        Scalar omega = D * (distanceVector * scvFace.unitOuterNormal());
+        omega *= problem.model().curVolVars(scv).extrusionFactor();
 
         return omega;
     }
-
-    const Problem &problem_() const
-    {
-        return *problemPtr_;
-    }
-
-    const SubControlVolumeFace& scvFace_() const
-    {
-        return *scvFacePtr_;
-    }
-
-    const Problem *problemPtr_;
-    const SubControlVolumeFace *scvFacePtr_; //!< Pointer to the sub control volume face for which the flux variables are created
-    std::set<IndexType> stencil_;         //!< Indices of the cells of which the pressure is needed for the flux calculation
-
-    //! Boundary volume variables (they only get updated on Dirichlet boundaries)
-    VolumeVariables* boundaryVolVars_;
-
-    IndexType phaseIdx_;
-    IndexType compIdx_;
-    //! Precomputed values
-    Scalar tij_; //!< transmissibility for the flux calculation
-    Scalar rhoDGradXNormal_; //! rho*D(grad(x))*n
 };
 
 } // end namespace
diff --git a/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh b/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh
index d639edadbe4afd983ddc57f7b925ebe7096bff2e..786220b7d3272b2fe1f61ff40bf395ce5e0e2a21 100644
--- a/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh
+++ b/test/porousmediumflow/3p3c/implicit/infiltrationproblem.hh
@@ -345,8 +345,11 @@ public:
 
         for (const auto& element : elements(this->gridView()))
         {
-            auto fvGeometry = this->model().fvGeometries(element);
-            for (auto&& scv : fvGeometry.scvs())
+            // make sure the FVElementGeometry is bound to the element
+            this->model().fvGeometries_().bindElement(element);
+            const auto& fvGeometry = this->model().fvGeometries(element);
+
+            for (const auto& scv : fvGeometry.scvs())
             {
                 auto dofIdxGlobal = scv.dofIndex();
                 (*Kxx)[dofIdxGlobal] = this->spatialParams().intrinsicPermeability(scv);