diff --git a/dumux/discretization/box/fvelementgeometry.hh b/dumux/discretization/box/fvelementgeometry.hh
index 926ef8c9ff2cc14c964afbdf95944d6399e866e9..3b06ab52300c9a3c9be04237a7ebe2067f52b025 100644
--- a/dumux/discretization/box/fvelementgeometry.hh
+++ b/dumux/discretization/box/fvelementgeometry.hh
@@ -148,6 +148,10 @@ public:
     const FVGridGeometry& fvGridGeometry() const
     { return *fvGridGeometryPtr_; }
 
+    //! Returns whether one of the geometry's scvfs lies on a boundary
+    bool hasBoundaryScvf() const
+    { return fvGridGeometry().hasBoundaryScvf(eIdx_); }
+
 private:
     const Element* elementPtr_;
     const FVGridGeometry* fvGridGeometryPtr_;
@@ -263,11 +267,16 @@ public:
     const FVGridGeometry& fvGridGeometry() const
     { return *fvGridGeometryPtr_; }
 
+    //! Returns whether one of the geometry's scvfs lies on a boundary
+    bool hasBoundaryScvf() const
+    { return hasBoundaryScvf_; }
+
 private:
 
     void makeElementGeometries(const Element& element)
     {
         auto eIdx = fvGridGeometry().elementMapper().index(element);
+        hasBoundaryScvf_ = false;
 
         // get the element geometry
         auto elementGeometry = element.geometry();
@@ -316,6 +325,7 @@ private:
             if (intersection.boundary() && !intersection.neighbor())
             {
                 const auto isGeometry = intersection.geometry();
+                hasBoundaryScvf_ = true;
 
                 for (unsigned int isScvfLocalIdx = 0; isScvfLocalIdx < isGeometry.corners(); ++isScvfLocalIdx)
                 {
@@ -348,6 +358,8 @@ private:
     //! vectors to store the geometries locally after binding an element
     std::vector<SubControlVolume> scvs_;
     std::vector<SubControlVolumeFace> scvfs_;
+
+    bool hasBoundaryScvf_ = false;
 };
 
 } // end namespace Dumux
diff --git a/dumux/discretization/box/fvgridgeometry.hh b/dumux/discretization/box/fvgridgeometry.hh
index 98ea2b78ecadcc33839ea7337be24681d6d25636..b6d0a3b8034b4caec8cfe74429c8110a898e78aa 100644
--- a/dumux/discretization/box/fvgridgeometry.hh
+++ b/dumux/discretization/box/fvgridgeometry.hh
@@ -156,6 +156,7 @@ public:
         auto numElements = this->gridView().size(0);
         scvs_.resize(numElements);
         scvfs_.resize(numElements);
+        hasBoundaryScvf_.resize(numElements, false);
 
         boundaryDofIndices_.assign(numDofs(), false);
 
@@ -215,6 +216,8 @@ public:
                 if (intersection.boundary() && !intersection.neighbor())
                 {
                     const auto isGeometry = intersection.geometry();
+                    hasBoundaryScvf_[eIdx] = true;
+
                     // count
                     numScvf_ += isGeometry.corners();
                     numBoundaryScvf_ += isGeometry.corners();
@@ -317,6 +320,10 @@ public:
     const std::unordered_map<std::size_t, std::size_t>& periodicVertexMap() const
     { return periodicVertexMap_; }
 
+    //! Returns whether one of the geometry's scvfs lies on a boundary
+    bool hasBoundaryScvf(std::size_t  eIdx) const
+    { return hasBoundaryScvf_[eIdx]; }
+
 private:
 
     const FeCache feCache_;
@@ -330,6 +337,7 @@ private:
 
     // vertices on the boudary
     std::vector<bool> boundaryDofIndices_;
+    std::vector<bool> hasBoundaryScvf_;
 
     // a map for periodic boundary vertices
     std::unordered_map<std::size_t, std::size_t> periodicVertexMap_;
diff --git a/dumux/discretization/cellcentered/gridvolumevariables.hh b/dumux/discretization/cellcentered/gridvolumevariables.hh
index b9b250ffd78686829fa1a33cfd07ec6a5724831d..21ad3d9edb5b61404b83ea8dbebdb64240af2861 100644
--- a/dumux/discretization/cellcentered/gridvolumevariables.hh
+++ b/dumux/discretization/cellcentered/gridvolumevariables.hh
@@ -80,25 +80,28 @@ public:
                 volumeVariables_[scv.dofIndex()].update(elemSol, problem(), element, scv);
             }
 
-            // handle the boundary volume variables
-            for (auto&& scvf : scvfs(fvGeometry))
+            if (fvGeometry.hasBoundaryScvf())
             {
-                // if we are not on a boundary, skip the rest
-                if (!scvf.boundary())
-                    continue;
-
-                // check if boundary is a pure dirichlet boundary
-                const auto bcTypes = problem().boundaryTypes(element, scvf);
-                if (bcTypes.hasOnlyDirichlet())
+                // handle the boundary volume variables
+                for (auto&& scvf : scvfs(fvGeometry))
                 {
-                    const auto insideScvIdx = scvf.insideScvIdx();
-                    const auto& insideScv = fvGeometry.scv(insideScvIdx);
-                    const auto dirichletPriVars = elementSolution<typename FVGridGeometry::LocalView>(problem().dirichlet(element, scvf));
-
-                    volumeVariables_[scvf.outsideScvIdx()].update(dirichletPriVars,
-                                                                  problem(),
-                                                                  element,
-                                                                  insideScv);
+                    // if we are not on a boundary, skip the rest
+                    if (!scvf.boundary())
+                        continue;
+
+                    // check if boundary is a pure dirichlet boundary
+                    const auto bcTypes = problem().boundaryTypes(element, scvf);
+                    if (bcTypes.hasOnlyDirichlet())
+                    {
+                        const auto insideScvIdx = scvf.insideScvIdx();
+                        const auto& insideScv = fvGeometry.scv(insideScvIdx);
+                        const auto dirichletPriVars = elementSolution<typename FVGridGeometry::LocalView>(problem().dirichlet(element, scvf));
+
+                        volumeVariables_[scvf.outsideScvIdx()].update(dirichletPriVars,
+                                                                      problem(),
+                                                                      element,
+                                                                      insideScv);
+                    }
                 }
             }
         }
diff --git a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh
index 926db5889c66147d7aa547fd2624361e9634d353..675980bf32850b70718c5f1f70b8803fb99f7284 100644
--- a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh
+++ b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh
@@ -163,28 +163,31 @@ public:
             volumeVariables_.reserve(numVolVars+maxNumBoundaryVolVars);
             volVarIndices_.reserve(numVolVars+maxNumBoundaryVolVars);
 
-            // treat the BCs inside the element
-            for (const auto& scvf : scvfs(fvGeometry))
+            if (fvGeometry.hasBoundaryScvf())
             {
-                // if we are not on a boundary, skip to the next scvf
-                if (!scvf.boundary())
-                    continue;
-
-                const auto bcTypes = problem.boundaryTypes(element, scvf);
-
-                // Only proceed on dirichlet boundaries. Fluxes across Neumann
-                // boundaries are never computed - the user-defined flux is taken.
-                if (bcTypes.hasOnlyDirichlet())
+                // treat the BCs inside the element
+                for (const auto& scvf : scvfs(fvGeometry))
                 {
-                    // boundary volume variables
-                    VolumeVariables dirichletVolVars;
-                    dirichletVolVars.update(elementSolution<FVElementGeometry>(problem.dirichlet(element, scvf)),
-                                            problem,
-                                            element,
-                                            scvI);
-
-                    volumeVariables_.emplace_back(std::move(dirichletVolVars));
-                    volVarIndices_.push_back(scvf.outsideScvIdx());
+                    // if we are not on a boundary, skip to the next scvf
+                    if (!scvf.boundary())
+                        continue;
+
+                    const auto bcTypes = problem.boundaryTypes(element, scvf);
+
+                    // Only proceed on dirichlet boundaries. Fluxes across Neumann
+                    // boundaries are never computed - the user-defined flux is taken.
+                    if (bcTypes.hasOnlyDirichlet())
+                    {
+                        // boundary volume variables
+                        VolumeVariables dirichletVolVars;
+                        dirichletVolVars.update(elementSolution<FVElementGeometry>(problem.dirichlet(element, scvf)),
+                                                problem,
+                                                element,
+                                                scvI);
+
+                        volumeVariables_.emplace_back(std::move(dirichletVolVars));
+                        volVarIndices_.push_back(scvf.outsideScvIdx());
+                    }
                 }
             }
 
diff --git a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
index 553f49de84bab780cc1f07aa7c9dae4df04c831e..56f8c37e31c04d439bd0ba9b8ae4691c17cd7a4a 100644
--- a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
@@ -154,6 +154,10 @@ public:
     const FVGridGeometry& fvGridGeometry() const
     { return *fvGridGeometryPtr_; }
 
+    //! Returns whether one of the geometry's scvfs lies on a boundary
+    bool hasBoundaryScvf() const
+    { return fvGridGeometry().hasBoundaryScvf(scvIndices_[0]); }
+
 private:
 
     std::array<GridIndexType, 1> scvIndices_;
@@ -309,6 +313,10 @@ public:
     const FVGridGeometry& fvGridGeometry() const
     { return *fvGridGeometryPtr_; }
 
+    //! Returns whether one of the geometry's scvfs lies on a boundary
+    bool hasBoundaryScvf() const
+    { return hasBoundaryScvf_; }
+
 private:
 
     //! Computes the number of neighboring scvfs that have to be prepared
@@ -386,6 +394,8 @@ private:
                 if (fvGridGeometry().isGhostVertex(vIdxGlobal))
                     continue;
 
+                hasBoundaryScvf_ = (hasBoundaryScvf_ || is.boundary());
+
                 scvfs_.emplace_back(MpfaHelper(),
                                     MpfaHelper::getScvfCorners(isPositions, numCorners, c),
                                     is.centerUnitOuterNormal(),
@@ -514,6 +524,8 @@ private:
         neighborScvfIndices_.clear();
         neighborScvs_.clear();
         neighborScvfs_.clear();
+
+        hasBoundaryScvf_ = false;
     }
 
     const FVGridGeometry* fvGridGeometryPtr_;
@@ -528,6 +540,8 @@ private:
     std::vector<GridIndexType> neighborScvfIndices_;
     std::vector<SubControlVolume> neighborScvs_;
     std::vector<SubControlVolumeFace> neighborScvfs_;
+
+    bool hasBoundaryScvf_ = false;
 };
 
 } // end namespace
diff --git a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
index 21a157d483fd56f8b6ce0a6e2c44fc7ae05d659e..a1f315aec636f5cf901abd20f2f102139e5bec44 100644
--- a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
@@ -188,6 +188,7 @@ public:
         scvs_.resize(numScvs);
         scvfs_.reserve(numScvf);
         scvfIndicesOfScv_.resize(numScvs);
+        hasBoundaryScvf_.resize(numScvs, false);
 
         // Some methods require to use a second type of interaction volume, e.g.
         // around vertices on the boundary or branching points (surface grids)
@@ -236,6 +237,9 @@ public:
                 const bool boundary = is.boundary();
                 const bool neighbor = is.neighbor();
 
+                if (boundary)
+                    hasBoundaryScvf_[eIdx] = true;
+
                 // for surface grids, skip the rest if handled already
                 if (dim < dimWorld && neighbor && outsideIndices[indexInInside].empty())
                     continue;
@@ -395,6 +399,10 @@ public:
     const SubControlVolumeFace& flipScvf(GridIndexType scvfIdx, unsigned int outsideScvfIdx = 0) const
     { return scvfs_[flipScvfIndices_[scvfIdx][outsideScvfIdx]]; }
 
+    //! Returns whether one of the geometry's scvfs lies on a boundary
+    bool hasBoundaryScvf(GridIndexType eIdx) const
+    { return hasBoundaryScvf_[eIdx]; }
+
 private:
     // connectivity map for efficient assembly
     ConnectivityMap connectivityMap_;
@@ -407,6 +415,7 @@ private:
     std::vector<std::vector<GridIndexType>> scvfIndicesOfScv_;
     std::vector<bool> secondaryInteractionVolumeVertices_;
     GridIndexType numBoundaryScvf_;
+    std::vector<bool> hasBoundaryScvf_;
 
     // needed for embedded surface and network grids (dim < dimWorld)
     FlipScvfIndexSet flipScvfIndices_;
diff --git a/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh
index 9612e6c1ea50e0020e9f6b669b233a7e9d0d72ea..15f70efdbe7c0446035c3cf71e6a6456e0925f6e 100644
--- a/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh
+++ b/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh
@@ -152,28 +152,31 @@ public:
             ++localIdx;
         }
 
-        // Update boundary volume variables
-        for (auto&& scvf : scvfs(fvGeometry))
+        if (fvGeometry.hasBoundaryScvf())
         {
-            // if we are not on a boundary, skip to the next scvf
-            if (!scvf.boundary())
+            // Update boundary volume variables
+            for (auto&& scvf : scvfs(fvGeometry))
+            {
+                // if we are not on a boundary, skip to the next scvf
+                if (!scvf.boundary())
                 continue;
 
-            // check if boundary is a pure dirichlet boundary
-            const auto bcTypes = problem.boundaryTypes(element, scvf);
-            if (bcTypes.hasOnlyDirichlet())
-            {
-                const auto dirichletPriVars = elementSolution<FVElementGeometry>(problem.dirichlet(element, scvf));
-
-                volumeVariables_.resize(localIdx+1);
-                volVarIndices_.resize(localIdx+1);
-                volumeVariables_[localIdx].update(dirichletPriVars,
-                                                  problem,
-                                                  element,
-                                                  scvI);
-                volVarIndices_[localIdx] = scvf.outsideScvIdx();
-                ++localIdx;
-            }
+                // check if boundary is a pure dirichlet boundary
+                const auto bcTypes = problem.boundaryTypes(element, scvf);
+                if (bcTypes.hasOnlyDirichlet())
+                {
+                    const auto dirichletPriVars = elementSolution<FVElementGeometry>(problem.dirichlet(element, scvf));
+
+                    volumeVariables_.resize(localIdx+1);
+                    volVarIndices_.resize(localIdx+1);
+                    volumeVariables_[localIdx].update(dirichletPriVars,
+                        problem,
+                        element,
+                        scvI);
+                        volVarIndices_[localIdx] = scvf.outsideScvIdx();
+                        ++localIdx;
+                    }
+                }
         }
 
         //! Check if user added additional DOF dependencies, i.e. the residual of DOF globalI depends
diff --git a/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
index c1e4b4388b1de6558445f1bf0c03ac832179b029..3ccaf0ef026651cf5dd796658dd1d49838ab4bad 100644
--- a/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
+++ b/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
@@ -156,6 +156,10 @@ public:
     const FVGridGeometry& fvGridGeometry() const
     { return *fvGridGeometryPtr_; }
 
+    //! Returns whether one of the geometry's scvfs lies on a boundary
+    bool hasBoundaryScvf() const
+    { return fvGridGeometry().hasBoundaryScvf(scvIndices_[0]); }
+
 private:
 
     const Element* elementPtr_;
@@ -351,6 +355,10 @@ public:
     const FVGridGeometry& fvGridGeometry() const
     { return *fvGridGeometryPtr_; }
 
+    //! Returns whether one of the geometry's scvfs lies on a boundary
+    bool hasBoundaryScvf() const
+    { return hasBoundaryScvf_; }
+
 private:
 
     IndexType findFlippedScvfIndex_(IndexType insideScvIdx, IndexType globalOutsideScvIdx)
@@ -412,11 +420,15 @@ private:
                 scvIndices.resize(scvfNeighborVolVarIndices.size() + 1);
                 scvIndices[0] = eIdx;
                 std::copy(scvfNeighborVolVarIndices.begin(), scvfNeighborVolVarIndices.end(), scvIndices.begin()+1);
+
+                const bool onBoundary = intersection.boundary() && !intersection.neighbor();
+                hasBoundaryScvf_ = (hasBoundaryScvf_ || onBoundary);
+
                 scvfs_.emplace_back(intersection,
                                     intersection.geometry(),
                                     scvFaceIndices[scvfCounter],
                                     scvIndices,
-                                    intersection.boundary() && !intersection.neighbor());
+                                    onBoundary);
                 scvfIndices_.emplace_back(scvFaceIndices[scvfCounter]);
                 scvfCounter++;
 
@@ -536,6 +548,8 @@ private:
         neighborScvs_.clear();
         neighborScvfs_.clear();
         flippedNeighborScvfIndices_.clear();
+
+        hasBoundaryScvf_ = false;
     }
 
     const Element* elementPtr_; //!< the element to which this fvgeometry is bound
@@ -555,6 +569,8 @@ private:
     std::vector<IndexType> neighborScvfIndices_;
     std::vector<SubControlVolumeFace> neighborScvfs_;
     std::vector<std::vector<IndexType>> flippedNeighborScvfIndices_;
+
+    bool hasBoundaryScvf_ = false;
 };
 
 } // end namespace Dumux
diff --git a/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh
index be19deafec40f9c5f412c2b2c0021a12dfec0e4f..9d59a04246533e3b67f09611aa6e304795290d69 100644
--- a/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh
+++ b/dumux/discretization/cellcentered/tpfa/fvgridgeometry.hh
@@ -172,6 +172,7 @@ public:
         scvs_.resize(numScvs);
         scvfs_.reserve(numScvf);
         scvfIndicesOfScv_.resize(numScvs);
+        hasBoundaryScvf_.resize(numScvs, false);
 
         // Build the scvs and scv faces
         IndexType scvfIdx = 0;
@@ -254,6 +255,8 @@ public:
                                         ScvfGridIndexStorage({eIdx, this->gridView().size(0) + numBoundaryScvf_++}),
                                         true);
                     scvfsIndexSet.push_back(scvfIdx++);
+
+                    hasBoundaryScvf_[eIdx] = true;
                 }
             }
 
@@ -314,6 +317,10 @@ public:
     const ConnectivityMap &connectivityMap() const
     { return connectivityMap_; }
 
+    //! Returns whether one of the geometry's scvfs lies on a boundary
+    bool hasBoundaryScvf(IndexType eIdx) const
+    { return hasBoundaryScvf_[eIdx]; }
+
 private:
     // find the scvf that has insideScvIdx in its outsideScvIdx list and outsideScvIdx as its insideScvIdx
     IndexType findFlippedScvfIndex_(IndexType insideScvIdx, IndexType outsideScvIdx)
@@ -338,6 +345,7 @@ private:
     std::vector<SubControlVolumeFace> scvfs_;
     std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
     IndexType numBoundaryScvf_;
+    std::vector<bool> hasBoundaryScvf_;
 
     //! needed for embedded surface and network grids (dim < dimWorld)
     std::vector<std::vector<IndexType>> flipScvfIndices_;
diff --git a/dumux/discretization/staggered/freeflow/elementvolumevariables.hh b/dumux/discretization/staggered/freeflow/elementvolumevariables.hh
index 420b66fed79e09848bd6367bf916b36d4e8a99af..70106b9fcaeae30a1611edff4a732fc69a4bf41b 100644
--- a/dumux/discretization/staggered/freeflow/elementvolumevariables.hh
+++ b/dumux/discretization/staggered/freeflow/elementvolumevariables.hh
@@ -91,7 +91,8 @@ public:
               const SolutionVector& sol)
     {
         // the last parameter {} is needed for the PassKey pattern which restricts access to the ElementVolVars class
-        gridVolVars().updateBoundary_(element, fvGeometry, sol, {});
+        if (fvGeometry.hasBoundaryScvf())
+            gridVolVars().updateBoundary_(element, fvGeometry, sol, {});
     }
 
     //! function to prepare the vol vars within the element
@@ -181,24 +182,27 @@ public:
             ++localIdx;
         }
 
-        // Update boundary volume variables
-        for (auto&& scvf : scvfs(fvGeometry))
+        if (fvGeometry.hasBoundaryScvf())
         {
-            // if we are not on a boundary, skip to the next scvf
-            if (!scvf.boundary())
-                continue;
-
-            volumeVariables_.resize(localIdx+1);
-            volVarIndices_.resize(localIdx+1);
-
-            auto boundaryPriVars = gridVolVars().getBoundaryPriVars(problem, sol, element, scvf);
-            auto elemSol = elementSolution<FVElementGeometry>(std::move(boundaryPriVars));
-            volumeVariables_[localIdx].update(elemSol,
-                                              problem,
-                                              element,
-                                              scvI);
-            volVarIndices_[localIdx] = scvf.outsideScvIdx();
-            ++localIdx;
+            // Update boundary volume variables
+            for (auto&& scvf : scvfs(fvGeometry))
+            {
+                // if we are not on a boundary, skip to the next scvf
+                if (!scvf.boundary())
+                    continue;
+
+                volumeVariables_.resize(localIdx+1);
+                volVarIndices_.resize(localIdx+1);
+
+                auto boundaryPriVars = gridVolVars().getBoundaryPriVars(problem, sol, element, scvf);
+                auto elemSol = elementSolution<FVElementGeometry>(std::move(boundaryPriVars));
+                volumeVariables_[localIdx].update(elemSol,
+                                                  problem,
+                                                  element,
+                                                  scvI);
+                volVarIndices_[localIdx] = scvf.outsideScvIdx();
+                ++localIdx;
+            }
         }
     }
 
diff --git a/dumux/discretization/staggered/fvelementgeometry.hh b/dumux/discretization/staggered/fvelementgeometry.hh
index 75d83d8ce3762b7308b1d2bf06c6437640072ca4..e34dea10ff41fcbc391180c990d08b08481fab77 100644
--- a/dumux/discretization/staggered/fvelementgeometry.hh
+++ b/dumux/discretization/staggered/fvelementgeometry.hh
@@ -208,6 +208,10 @@ public:
     const FVGridGeometry& fvGridGeometry() const
     { return *fvGridGeometryPtr_; }
 
+    //! Returns whether one of the geometry's scvfs lies on a boundary
+    bool hasBoundaryScvf() const
+    { return hasBoundaryScvf_; }
+
 private:
 
     //! create scvs and scvfs of the bound element
@@ -238,6 +242,9 @@ private:
                                     geometryHelper);
                 scvfIndices_.emplace_back(scvFaceIndices[scvfCounter]);
                 scvfCounter++;
+
+                if (intersection.boundary())
+                    hasBoundaryScvf_ = true;
             }
         }
     }
@@ -301,6 +308,8 @@ private:
         neighborScvfIndices_.clear();
         neighborScvs_.clear();
         neighborScvfs_.clear();
+
+        hasBoundaryScvf_ = false;
     }
 
     const Element* elementPtr_; //!< the element to which this fvgeometry is bound
@@ -318,6 +327,8 @@ private:
 
     std::vector<IndexType> neighborScvfIndices_;
     std::vector<SubControlVolumeFace> neighborScvfs_;
+
+    bool hasBoundaryScvf_ = false;
 };
 
 
diff --git a/dumux/discretization/staggered/fvgridgeometry.hh b/dumux/discretization/staggered/fvgridgeometry.hh
index c22674ed6f3aeb41bc038a5fd6bb1ba0a2941130..8535d30fd591b487cc7f623be22d500cf6c788ca 100644
--- a/dumux/discretization/staggered/fvgridgeometry.hh
+++ b/dumux/discretization/staggered/fvgridgeometry.hh
@@ -281,6 +281,7 @@ public:
         scvfs_.reserve(numScvf);
         scvfIndicesOfScv_.resize(numScvs);
         localToGlobalScvfIndices_.resize(numScvs);
+        hasBoundaryScvf_.resize(numScvs, false);
 
         // Build the scvs and scv faces
         IndexType scvfIdx = 0;
@@ -328,6 +329,8 @@ public:
                                         geometryHelper);
                     localToGlobalScvfIndices_[eIdx][localFaceIndex] = scvfIdx;
                     scvfsIndexSet.push_back(scvfIdx++);
+
+                    hasBoundaryScvf_[eIdx] = true;
                 }
             }
 
@@ -398,6 +401,10 @@ public:
         return FaceFVGridGeometry<ThisType>(this);
     }
 
+    //! Returns whether one of the geometry's scvfs lies on a boundary
+    bool hasBoundaryScvf(IndexType eIdx) const
+    { return hasBoundaryScvf_[eIdx]; }
+
 private:
 
     // mappers
@@ -409,6 +416,7 @@ private:
     std::vector<std::vector<IndexType>> scvfIndicesOfScv_;
     std::vector<std::vector<IndexType>> localToGlobalScvfIndices_;
     IndexType numBoundaryScvf_;
+    std::vector<bool> hasBoundaryScvf_;
 };
 
 /*!
diff --git a/test/discretization/box/CMakeLists.txt b/test/discretization/box/CMakeLists.txt
index 5bbe00764e5a6e9257a7cefbc4e632ca8dc9d390..3dc0efdbc24fbef5ac1db53190c3d27ec45c0606 100644
--- a/test/discretization/box/CMakeLists.txt
+++ b/test/discretization/box/CMakeLists.txt
@@ -1,12 +1,12 @@
 dune_add_test(NAME test_boxfvgeometry
               SOURCES test_boxfvgeometry.cc
               COMPILE_DEFINITIONS ENABLE_CACHING=false
-              LABELS unit)
+              LABELS unit discretization)
 
 dune_add_test(NAME test_boxfvgeometry_caching
               SOURCES test_boxfvgeometry.cc
               COMPILE_DEFINITIONS ENABLE_CACHING=true
-              LABELS unit)
+              LABELS unit discretization)
 
 #install sources
 install(FILES
diff --git a/test/discretization/box/test_boxfvgeometry.cc b/test/discretization/box/test_boxfvgeometry.cc
index 227191bc55f96971dc174bc99da585526df45b2f..b2130b25dacfd7e010c31bf630bb5a9e47ebba95 100644
--- a/test/discretization/box/test_boxfvgeometry.cc
+++ b/test/discretization/box/test_boxfvgeometry.cc
@@ -98,12 +98,21 @@ int main (int argc, char *argv[]) try
         if(0 != testForwardIterator(range2.begin(), range2.end(), op2))
             DUNE_THROW(Dune::Exception, "Iterator does not fulfill the forward iterator concept");
 
+        std::size_t boundaryCount = 0;
         for (auto&& scvf : scvfs(fvGeometry))
         {
             std::cout << "-- scvf " << scvf.index() << " ip at: " << scvf.ipGlobal() << " normal: " << scvf.unitOuterNormal();
-            if (scvf.boundary()) std::cout << " (on boundary).";
+            if (scvf.boundary())
+            {
+                ++boundaryCount;
+                std::cout << " (on boundary).";
+            }
             std::cout << std::endl;
         }
+
+        if ((boundaryCount>0) != fvGeometry.hasBoundaryScvf())
+            DUNE_THROW(Dune::InvalidStateException, "fvGeometry.hasBoundaryScvf() reports " << fvGeometry.hasBoundaryScvf()
+                            << " but the number of boundary scvfs is " << boundaryCount);
     }
 }
 // //////////////////////////////////
diff --git a/test/discretization/cellcentered/tpfa/CMakeLists.txt b/test/discretization/cellcentered/tpfa/CMakeLists.txt
index 813ea9e4e75ed418f6e2550d4f3251c849f44379..3acfe25990fa4461b896de4585c45a3b20096989 100644
--- a/test/discretization/cellcentered/tpfa/CMakeLists.txt
+++ b/test/discretization/cellcentered/tpfa/CMakeLists.txt
@@ -1,21 +1,21 @@
 dune_add_test(NAME test_tpfafvgeometry
               SOURCES test_tpfafvgeometry.cc
               COMPILE_DEFINITIONS ENABLE_CACHING=false
-              LABELS unit)
+              LABELS unit discretization)
 
 dune_add_test(NAME test_tpfafvgeometry_caching
               SOURCES test_tpfafvgeometry.cc
               COMPILE_DEFINITIONS ENABLE_CACHING=true
-              LABELS unit)
+              LABELS unit discretization)
 
 dune_add_test(NAME test_tpfafvgeometry_nonconforming
               SOURCES test_tpfafvgeometry_nonconforming.cc
               COMPILE_DEFINITIONS ENABLE_CACHING=false
               CMAKE_GUARD dune-alugrid_FOUND
-              LABELS unit)
+              LABELS unit discretization)
 
 dune_add_test(NAME test_cachedtpfafvgeometry_nonconforming
               SOURCES test_tpfafvgeometry_nonconforming.cc
               COMPILE_DEFINITIONS ENABLE_CACHING=true
               CMAKE_GUARD dune-alugrid_FOUND
-              LABELS unit)
+              LABELS unit discretization)
diff --git a/test/discretization/cellcentered/tpfa/test_tpfafvgeometry.cc b/test/discretization/cellcentered/tpfa/test_tpfafvgeometry.cc
index e7223e64bbef2895fc380130cb67e121f565ba7e..49d715f3ff45ce05bd0c132f6695d3a6fc85fa52 100644
--- a/test/discretization/cellcentered/tpfa/test_tpfafvgeometry.cc
+++ b/test/discretization/cellcentered/tpfa/test_tpfafvgeometry.cc
@@ -99,12 +99,21 @@ int main (int argc, char *argv[]) try
         if(0 != testForwardIterator(range2.begin(), range2.end(), op2))
             DUNE_THROW(Dune::Exception, "Iterator does not fulfill the forward iterator concept");
 
+        std::size_t boundaryCount = 0;
         for (auto&& scvf : scvfs(fvGeometry))
         {
             std::cout << "-- scvf " << scvf.index() << " ip at: " << scvf.ipGlobal();
-            if (scvf.boundary()) std::cout << " (on boundary).";
+            if (scvf.boundary())
+            {
+                ++boundaryCount;
+                std::cout << " (on boundary).";
+            }
             std::cout << std::endl;
         }
+
+        if ((boundaryCount>0) != fvGeometry.hasBoundaryScvf())
+            DUNE_THROW(Dune::InvalidStateException, "fvGeometry.hasBoundaryScvf() reports " << fvGeometry.hasBoundaryScvf()
+                            << " but the number of boundary scvfs is " << boundaryCount);
     }
 }
 // //////////////////////////////////
diff --git a/test/discretization/staggered/CMakeLists.txt b/test/discretization/staggered/CMakeLists.txt
index e5e3e4fcc918b20e68a7a3fd7c933299976b3bf1..cab4a791fdb6ab9db4f2a19032d9c16b2ec4f28f 100644
--- a/test/discretization/staggered/CMakeLists.txt
+++ b/test/discretization/staggered/CMakeLists.txt
@@ -1,10 +1,10 @@
 dune_add_test(NAME test_staggeredfvgeometry
               SOURCES test_staggeredfvgeometry.cc
-              LABELS unit)
+              LABELS unit discretization)
 
 dune_add_test(NAME test_staggered_free_flow_geometry
               SOURCES test_staggered_free_flow_geometry.cc
-              LABELS unit)
+              LABELS unit discretization)
 
 #install sources
 install(FILES
diff --git a/test/discretization/staggered/test_staggeredfvgeometry.cc b/test/discretization/staggered/test_staggeredfvgeometry.cc
index d3d45349e8bac2efb59c80a166d11b4253ac3793..130e0be2df6cc3d2028dfeaaaf67b4e1c9427413 100644
--- a/test/discretization/staggered/test_staggeredfvgeometry.cc
+++ b/test/discretization/staggered/test_staggeredfvgeometry.cc
@@ -133,12 +133,21 @@ int main (int argc, char *argv[]) try
         if(0 != testForwardIterator(range2.begin(), range2.end(), op2))
             DUNE_THROW(Dune::Exception, "Iterator does not fulfill the forward iterator concept");
 
+        std::size_t boundaryCount = 0;
         for (auto&& scvf : scvfs(fvGeometry))
         {
             std::cout << "-- scvf " << scvf.index() << " ip at: " << scvf.ipGlobal() << " normal: " << scvf.unitOuterNormal();
-            if (scvf.boundary()) std::cout << " (on boundary).";
+            if (scvf.boundary())
+            {
+                ++boundaryCount;
+                std::cout << " (on boundary).";
+            }
             std::cout << std::endl;
         }
+
+        if ((boundaryCount>0) != fvGeometry.hasBoundaryScvf())
+            DUNE_THROW(Dune::InvalidStateException, "fvGeometry.hasBoundaryScvf() reports " << fvGeometry.hasBoundaryScvf()
+                            << " but the number of boundary scvfs is " << boundaryCount);
     }
 }
 // //////////////////////////////////