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_;