diff --git a/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh b/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
index 2288f3b49d00aa2d83089d6d30224ae23928162d..84aca7ed72008982cba4e459729476ca3919caa1 100644
--- a/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
+++ b/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
@@ -80,49 +80,41 @@ public:
     template<typename SubControlVolumeFace>
     void insert(const SubControlVolumeFace& scvf)
     {
-        insert(scvf.boundary(),
-               scvf.index(),
+        insert(scvf.index(),
                scvf.insideScvIdx(),
-               scvf.outsideScvIndices());
+               scvf.boundary());
     }
 
     //! Inserts scvf data
-    template<typename OutsideGridIndexStorage>
-    void insert(const bool boundary,
-                const GridIndexType scvfIdx,
+    void insert(const GridIndexType scvfIdx,
                 const GridIndexType insideScvIdx,
-                const OutsideGridIndexStorage& outsideScvIndices)
+                const bool boundary)
     {
         // this should always be called only once per scvf
-        assert(std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx ) == scvfIndices_.end()
-               && "scvf has already been inserted!");
+        assert(std::count(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx ) == 0 && "scvf already inserted!");
 
         // the local index of the scvf data about to be inserted
         const LocalIndexType curScvfLocalIdx = scvfIndices_.size();
 
-        // add grid scvf data
-        ScvfNeighborIndexSet scvIndices;
-        scvIndices.push_back(insideScvIdx);
-        for (const auto& outsideIdx : outsideScvIndices)
-                scvIndices.push_back(outsideIdx);
-
         // if scvf is on boundary, increase counter
-        if (boundary)
-            numBoundaryScvfs_++;
+        if (boundary) numBoundaryScvfs_++;
 
         // insert data on the new scv
         scvfIndices_.push_back(scvfIdx);
         scvfIsOnBoundary_.push_back(boundary);
-        scvfNeighborScvIndices_.push_back(scvIndices);
 
-        // if entry for the inside scv exists, append scvf local index, create otherwise
+        // if entry for the inside scv exists append data, create otherwise
         auto it = std::find( scvIndices_.begin(), scvIndices_.end(), insideScvIdx );
         if (it != scvIndices_.end())
-            localScvfIndicesInScv_[ std::distance(scvIndices_.begin(), it) ].push_back(curScvfLocalIdx);
+        {
+            const auto scvIdxLocal = std::distance(scvIndices_.begin(), it);
+            scvfInsideScvIndices_.push_back(scvIdxLocal);
+            localScvfIndicesInScv_[scvIdxLocal].push_back(curScvfLocalIdx);
+        }
         else
         {
-            localScvfIndicesInScv_.push_back({});
-            localScvfIndicesInScv_.back().push_back(curScvfLocalIdx);
+            scvfInsideScvIndices_.push_back(scvIndices_.size());
+            localScvfIndicesInScv_.push_back({curScvfLocalIdx});
             scvIndices_.push_back(insideScvIdx);
         }
     }
@@ -136,49 +128,64 @@ public:
     //! returns the number of boundary scvfs around the node
     std::size_t numBoundaryScvfs() const { return numBoundaryScvfs_; }
 
-    //! returns the global scv indices connected to this dual grid node
+    //! returns the grid scv indices connected to this dual grid node
     const GridStencilType& globalScvIndices() const { return scvIndices_; }
 
-    //! returns the global scvf indices connected to this dual grid node
+    //! returns the grid scvf indices connected to this dual grid node
     const GridScvfStencilType& globalScvfIndices() const { return scvfIndices_; }
 
     //! returns whether or not the i-th scvf is on a domain boundary
-    bool scvfIsOnBoundary(unsigned int i) const { return scvfIsOnBoundary_[i]; }
+    bool scvfIsOnBoundary(unsigned int i) const
+    {
+        assert(i < numScvfs());
+        return scvfIsOnBoundary_[i];
+    }
 
-    //! returns the global scv idx of the i-th scv
-    GridIndexType scvIdxGlobal(unsigned int i) const { return scvIndices_[i]; }
+    //! returns the grid scv idx of the i-th scv
+    GridIndexType scvIdxGlobal(unsigned int i) const
+    {
+        assert(i < numScvs());
+        return scvIndices_[i];
+    }
 
     //! returns the index of the i-th scvf
-    GridIndexType scvfIdxGlobal(unsigned int i) const { return scvfIndices_[i]; }
+    GridIndexType scvfIdxGlobal(unsigned int i) const
+    {
+        assert(i < numScvfs());
+        return scvfIndices_[i];
+    }
 
-    //! returns the global index of the j-th scvf embedded in the i-th scv
+    //! returns the grid index of the j-th scvf embedded in the i-th scv
     GridIndexType scvfIdxGlobal(unsigned int i, unsigned int j) const
     {
-        assert(j < dim); // only dim faces can be embedded in an scv
+        assert(i < numScvs());
+        assert(j < localScvfIndicesInScv_[i].size());
         return scvfIndices_[ localScvfIndicesInScv_[i][j] ];
     }
 
     //! returns the node-local index of the j-th scvf embedded in the i-th scv
     LocalIndexType scvfIdxLocal(unsigned int i, unsigned int j) const
     {
-        assert(j < dim); // only dim faces can be embedded in an scv
+        assert(i < numScvs());
+        assert(j < localScvfIndicesInScv_[i].size());
         return localScvfIndicesInScv_[i][j];
     }
 
-    //! returns the indices of the neighboring scvs of the i-th scvf
-    const ScvfNeighborIndexSet& neighboringScvIndices(unsigned int i) const
-    { return scvfNeighborScvIndices_[i]; }
+    //! returns the node-local index of the inside scv of the i-th scvf
+    LocalIndexType insideScvIdxLocal(unsigned int i) const
+    {
+        assert(i < numScvfs());
+        return scvfInsideScvIndices_[i];
+    }
 
 private:
     GridStencilType scvIndices_;                                                       //!< The indices of the scvs around a dual grid node
     Dune::ReservedVector<DimIndexVector, maxNumElementsAtNode> localScvfIndicesInScv_; //!< Maps to each scv a list of scvf indices embedded in it
 
-    GridScvfStencilType scvfIndices_;                             //!< the indices of the scvfs around a dual grid node
-    std::size_t numBoundaryScvfs_;                                //!< stores how many boundary scvfs are embedded in this dual grid node
+    GridScvfStencilType scvfIndices_;                                //!< the indices of the scvfs around a dual grid node
+    std::size_t numBoundaryScvfs_;                                   //!< stores how many boundary scvfs are embedded in this dual grid node
     Dune::ReservedVector<bool, maxNumScvfsAtNode> scvfIsOnBoundary_; //!< Maps to each scvf a boolean to indicate if it is on the boundary
-
-    //! The neighboring scvs for the scvfs (order: 0 - inside, 1..n - outside)
-    Dune::ReservedVector<ScvfNeighborIndexSet, maxNumScvfsAtNode> scvfNeighborScvIndices_;
+    Dune::ReservedVector<LocalIndexType, maxNumScvfsAtNode> scvfInsideScvIndices_; //!< The inside local scv index for each scvf
 };
 
 /*!
diff --git a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
index e56f3b1dd7b7c797409f57fb93ecd1b58e289c73..a23d1674f7608d8b33e8531a429982d74023ac3a 100644
--- a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
@@ -74,6 +74,8 @@ class CCMpfaFVGridGeometry<GV, Traits, true>
     using ScvfOutsideGridIndexStorage = typename Traits::SubControlVolumeFace::Traits::OutsideGridIndexStorage;
 
 public:
+    //! export the flip scvf index set type
+    using FlipScvfIndexSet = std::vector<ScvfOutsideGridIndexStorage>;
     //! export the mpfa helper type
     using MpfaHelper = typename Traits::MpfaHelper;
     //! export the grid interaction volume index set type
@@ -302,32 +304,29 @@ public:
         }
 
         // Make the flip index set for network and surface grids
-        if (dim < dimWorld)
+        flipScvfIndices_.resize(scvfs_.size());
+        for (const auto& scvf : scvfs_)
         {
-            flipScvfIndices_.resize(scvfs_.size());
-            for (const auto& scvf : scvfs_)
-            {
-                if (scvf.boundary())
-                    continue;
+            if (scvf.boundary())
+                continue;
 
-                const auto numOutsideScvs = scvf.numOutsideScvs();
-                const auto vIdxGlobal = scvf.vertexIndex();
-                const auto insideScvIdx = scvf.insideScvIdx();
+            const auto numOutsideScvs = scvf.numOutsideScvs();
+            const auto vIdxGlobal = scvf.vertexIndex();
+            const auto insideScvIdx = scvf.insideScvIdx();
 
-                flipScvfIndices_[scvf.index()].resize(numOutsideScvs);
-                for (std::size_t i = 0; i < numOutsideScvs; ++i)
+            flipScvfIndices_[scvf.index()].resize(numOutsideScvs);
+            for (std::size_t i = 0; i < numOutsideScvs; ++i)
+            {
+                const auto outsideScvIdx = scvf.outsideScvIdx(i);
+                for (auto outsideScvfIndex : scvfIndicesOfScv_[outsideScvIdx])
                 {
-                    const auto outsideScvIdx = scvf.outsideScvIdx(i);
-                    for (auto outsideScvfIndex : scvfIndicesOfScv_[outsideScvIdx])
+                    const auto& outsideScvf = this->scvf(outsideScvfIndex);
+                    if (outsideScvf.vertexIndex() == vIdxGlobal &&
+                        MpfaHelper::vectorContainsValue(outsideScvf.outsideScvIndices(), insideScvIdx))
                     {
-                        const auto& outsideScvf = this->scvf(outsideScvfIndex);
-                        if (outsideScvf.vertexIndex() == vIdxGlobal &&
-                            MpfaHelper::vectorContainsValue(outsideScvf.outsideScvIndices(), insideScvIdx))
-                        {
-                            flipScvfIndices_[scvf.index()][i] = outsideScvfIndex;
-                            // there is always only one flip face in an outside element
-                            break;
-                        }
+                        flipScvfIndices_[scvf.index()][i] = outsideScvfIndex;
+                        // there is always only one flip face in an outside element
+                        break;
                     }
                 }
             }
@@ -360,15 +359,17 @@ public:
     //! Returns the grid interaction volume index set class.
     const GridIVIndexSets& gridInteractionVolumeIndexSets() const { return ivIndexSets_; }
 
+    //! Get the sub control volume face indices of an scv by global index
+    const std::vector<GridIndexType>& scvfIndicesOfScv(GridIndexType scvIdx) const { return scvfIndicesOfScv_[scvIdx]; }
+
+    //! Returns the flip scvf index set
+    const FlipScvfIndexSet& flipScvfIndexSet() const { return flipScvfIndices_; }
+
     //! Get the scvf on the same face but from the other side
     //! Note that e.g. the normals might be different in the case of surface grids
     const SubControlVolumeFace& flipScvf(GridIndexType scvfIdx, unsigned int outsideScvfIdx = 0) const
     { return scvfs_[flipScvfIndices_[scvfIdx][outsideScvfIdx]]; }
 
-    //! Get the sub control volume face indices of an scv by global index
-    const std::vector<GridIndexType>& scvfIndicesOfScv(GridIndexType scvIdx) const
-    { return scvfIndicesOfScv_[scvIdx]; }
-
 private:
     // connectivity map for efficient assembly
     ConnectivityMap connectivityMap_;
@@ -383,7 +384,7 @@ private:
     GridIndexType numBoundaryScvf_;
 
     // needed for embedded surface and network grids (dim < dimWorld)
-    std::vector<ScvfOutsideGridIndexStorage> flipScvfIndices_;
+    FlipScvfIndexSet flipScvfIndices_;
 
     // The grid interaction volume index set
     GridIVIndexSets ivIndexSets_;
@@ -419,6 +420,8 @@ class CCMpfaFVGridGeometry<GV, Traits, false>
     using ScvfOutsideGridIndexStorage = typename Traits::SubControlVolumeFace::Traits::OutsideGridIndexStorage;
 
 public:
+    //! export the flip scvf index set type
+    using FlipScvfIndexSet = std::vector<ScvfOutsideGridIndexStorage>;
     //! export the mpfa helper type
     using MpfaHelper = typename Traits::MpfaHelper;
     //! export the grid interaction volume index set type
@@ -616,7 +619,7 @@ public:
                                                     } ();
 
                     // insert the scvf data into the dual grid index set
-                    dualIdSet[vIdxGlobal].insert(boundary, numScvf_, eIdx, outsideScvIndices);
+                    dualIdSet[vIdxGlobal].insert(numScvf_, eIdx, boundary);
 
                     // store information on the scv face
                     scvfsIndexSet.push_back(numScvf_++);
@@ -635,40 +638,37 @@ public:
             neighborVolVarIndices_[eIdx] = neighborVolVarIndexSet;
         }
 
-        // Make the flip index set for network and surface grids
-        if (dim < dimWorld)
+        // Make the flip scvf index set
+        flipScvfIndices_.resize(numScvf_);
+        for (std::size_t scvIdx = 0; scvIdx < numScvs_; ++scvIdx)
         {
-            flipScvfIndices_.resize(numScvf_);
-            for (std::size_t scvIdx = 0; scvIdx < numScvs_; ++scvIdx)
+            const auto& scvfIndices = scvfIndicesOfScv_[scvIdx];
+            for (unsigned int i = 0; i < scvfIndices.size(); ++i)
             {
-                const auto& scvfIndices = scvfIndicesOfScv_[scvIdx];
-                for (unsigned int i = 0; i < scvfIndices.size(); ++i)
-                {
-                    // boundary scvf have no flip scvfs
-                    if (scvfIsOnBoundary[ scvfIndices[i] ])
-                        continue;
+                // boundary scvf have no flip scvfs
+                if (scvfIsOnBoundary[ scvfIndices[i] ])
+                    continue;
 
-                    const auto scvfIdx = scvfIndices[i];
-                    const auto vIdxGlobal = scvfVertexIndex[scvfIdx];
-                    const auto numOutsideScvs = neighborVolVarIndices_[scvIdx][i].size();
+                const auto scvfIdx = scvfIndices[i];
+                const auto vIdxGlobal = scvfVertexIndex[scvfIdx];
+                const auto numOutsideScvs = neighborVolVarIndices_[scvIdx][i].size();
 
-                    flipScvfIndices_[scvfIdx].resize(numOutsideScvs);
-                    for (unsigned int j = 0; j < numOutsideScvs; ++j)
+                flipScvfIndices_[scvfIdx].resize(numOutsideScvs);
+                for (unsigned int j = 0; j < numOutsideScvs; ++j)
+                {
+                    const auto outsideScvIdx = neighborVolVarIndices_[scvIdx][i][j];
+                    const auto& outsideScvfIndices = scvfIndicesOfScv_[outsideScvIdx];
+                    for (unsigned int k = 0; k < outsideScvfIndices.size(); ++k)
                     {
-                        const auto outsideScvIdx = neighborVolVarIndices_[scvIdx][i][j];
-                        const auto& outsideScvfIndices = scvfIndicesOfScv_[outsideScvIdx];
-                        for (unsigned int k = 0; k < outsideScvfIndices.size(); ++k)
+                        const auto outsideScvfIndex = outsideScvfIndices[k];
+                        const auto outsideScvfVertexIndex = scvfVertexIndex[outsideScvfIndex];
+                        const auto& outsideScvfNeighborIndices = neighborVolVarIndices_[outsideScvIdx][k];
+                        if (outsideScvfVertexIndex == vIdxGlobal &&
+                            MpfaHelper::vectorContainsValue(outsideScvfNeighborIndices, scvIdx))
                         {
-                            const auto outsideScvfIndex = outsideScvfIndices[k];
-                            const auto outsideScvfVertexIndex = scvfVertexIndex[outsideScvfIndex];
-                            const auto& outsideScvfNeighborIndices = neighborVolVarIndices_[outsideScvIdx][k];
-                            if (outsideScvfVertexIndex == vIdxGlobal &&
-                                MpfaHelper::vectorContainsValue(outsideScvfNeighborIndices, scvIdx))
-                            {
-                                flipScvfIndices_[scvfIdx][j] = outsideScvfIndex;
-                                // there is always only one flip face in an outside element
-                                break;
-                            }
+                            flipScvfIndices_[scvfIdx][j] = outsideScvfIndex;
+                            // there is always only one flip face in an outside element
+                            break;
                         }
                     }
                 }
@@ -702,6 +702,9 @@ public:
     const GridIndexType flipScvfIdx(GridIndexType scvfIdx, unsigned int outsideScvfIdx = 0) const
     { return flipScvfIndices_[scvfIdx][outsideScvfIdx]; }
 
+    //! Returns the flip scvf index set
+    const FlipScvfIndexSet& flipScvfIndexSet() const { return flipScvfIndices_; }
+
     //! Returns the connectivity map of which dofs
     //! have derivatives with respect to a given dof.
     const ConnectivityMap& connectivityMap() const { return connectivityMap_; }
@@ -723,7 +726,7 @@ private:
     GridIndexType numBoundaryScvf_;
 
     // needed for embedded surface and network grids (dim < dimWorld)
-    std::vector<ScvfOutsideGridIndexStorage> flipScvfIndices_;
+    FlipScvfIndexSet flipScvfIndices_;
 
     // The grid interaction volume index set
     GridIVIndexSets ivIndexSets_;
diff --git a/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh b/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh
index 7a6875ca5c1720aaaec215710a02812c877a4982..358e0b0fe1a2c7ec93329861d41835245293a866 100644
--- a/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh
+++ b/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh
@@ -93,9 +93,15 @@ public:
         {
             const auto vIdxGlobal = fvGridGeometry.vertexMapper().index(vertex);
             if (!fvGridGeometry.vertexUsesSecondaryInteractionVolume(vIdxGlobal))
-                PrimaryInteractionVolume::addInteractionVolumeIndexSets(primaryIVIndexSets_, scvfIndexMap_, (*dualGridIndexSet_)[vIdxGlobal]);
+                PrimaryInteractionVolume::addInteractionVolumeIndexSets(primaryIVIndexSets_,
+                                                                        scvfIndexMap_,
+                                                                        (*dualGridIndexSet_)[vIdxGlobal],
+                                                                        fvGridGeometry.flipScvfIndexSet());
             else
-                SecondaryInteractionVolume::addInteractionVolumeIndexSets(secondaryIVIndexSets_, scvfIndexMap_, (*dualGridIndexSet_)[vIdxGlobal]);
+                SecondaryInteractionVolume::addInteractionVolumeIndexSets(secondaryIVIndexSets_,
+                                                                          scvfIndexMap_,
+                                                                          (*dualGridIndexSet_)[vIdxGlobal],
+                                                                          fvGridGeometry.flipScvfIndexSet());
         }
     }
 
diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
index d5e46c9b0ff4a71744be6ffe45a8dd158bcac447..a4600b12b58ae2db93633b5c1936d066910cdbab 100644
--- a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
+++ b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
@@ -142,11 +142,15 @@ public:
 
     //! adds the iv index sets living around a vertex to a given container
     //! and stores the the corresponding index in a map for each scvf
-    template<class IvIndexSetContainer, class ScvfIndexMap, class NodalIndexSet>
+    template< class IvIndexSetContainer,
+              class ScvfIndexMap,
+              class NodalIndexSet,
+              class FlipScvfIndexSet >
     static void addInteractionVolumeIndexSets(IvIndexSetContainer& ivIndexSetContainer,
                                               ScvfIndexMap& scvfIndexMap,
-                                              const NodalIndexSet& nodalIndexSet)
-    { Impl::addInteractionVolumeIndexSets(ivIndexSetContainer, scvfIndexMap, nodalIndexSet); }
+                                              const NodalIndexSet& nodalIndexSet,
+                                              const FlipScvfIndexSet& flipScvfIndexSet)
+    { Impl::addInteractionVolumeIndexSets(ivIndexSetContainer, scvfIndexMap, nodalIndexSet, flipScvfIndexSet); }
 };
 
 } // end namespace Dumux
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
index 3091b94868fe550a7a639c9e2175ccbc543ed640..4629683fd96bf1c6e99e391cc92df4dd543bed6a 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
@@ -300,16 +300,20 @@ public:
 
     //! adds the iv index sets living around a vertex to a given container
     //! and stores the the corresponding index in a map for each scvf
-    template<class IvIndexSetContainer, class ScvfIndexMap, class NodalIndexSet>
+    template< class IvIndexSetContainer,
+              class ScvfIndexMap,
+              class NodalIndexSet,
+              class FlipScvfIndexSet >
     static void addInteractionVolumeIndexSets(IvIndexSetContainer& ivIndexSetContainer,
                                               ScvfIndexMap& scvfIndexMap,
-                                              const NodalIndexSet& nodalIndexSet)
+                                              const NodalIndexSet& nodalIndexSet,
+                                              const FlipScvfIndexSet& flipScvfIndexSet)
     {
         // the global index of the iv index set that is about to be created
         const auto curGlobalIndex = ivIndexSetContainer.size();
 
         // make the one index set for this node
-        ivIndexSetContainer.emplace_back(nodalIndexSet);
+        ivIndexSetContainer.emplace_back(nodalIndexSet, flipScvfIndexSet);
 
         // store the index mapping
         for (const auto scvfIdx : nodalIndexSet.globalScvfIndices())
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh
index 1d63d804f259675a3619111b2056aa24dcd9430f..86f11c045200f9529c045a7c49fd1635019b5f33 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh
@@ -56,10 +56,10 @@ public:
     using ScvfNeighborLocalIndexSet = typename DualGridNodalIndexSet::ScvfNeighborLocalIndexSet;
 
     //! The constructor
-    CCMpfaOInteractionVolumeIndexSet(const NodalIndexSet& nodalIndexSet) : nodalIndexSet_(nodalIndexSet)
+    template< class FlipScvfIndexSet >
+    CCMpfaOInteractionVolumeIndexSet(const NodalIndexSet& nodalIndexSet, const FlipScvfIndexSet& flipScvfIndexSet)
+    : nodalIndexSet_(nodalIndexSet)
     {
-        // determine the number of iv-local faces for memory reservation
-        // note that this might be a vast overestimation on surface grids!
         const auto numNodalScvfs = nodalIndexSet.numScvfs();
 
         // keeps track of which nodal scvfs have been handled already
@@ -77,6 +77,7 @@ public:
             // for scvfs touching the boundary there are no "outside" scvfs
             if (nodalIndexSet.scvfIsOnBoundary(i))
             {
+                scvfNeighborScvLocalIndices_.push_back({nodalIndexSet.insideScvIdxLocal(i)});
                 nodeToIvScvf_[i] = ivToNodeScvf_.size();
                 isHandled[i] = true;
                 ivToNodeScvf_.push_back(i);
@@ -84,79 +85,49 @@ public:
                 continue;
             }
 
-            // We insert a new iv-local face and find all "outside" scvfs that map
-            // to this face as well by comparing the set of neighboring scv indices.
-            const auto scvIndices = [&nodalIndexSet, i] ()
-                                    {
-                                        auto tmp = nodalIndexSet.neighboringScvIndices(i);
-                                        std::sort(tmp.begin(), tmp.end());
-                                        return tmp;
-                                    } ();
-            const auto numNeighborsI = scvIndices.size();
-
-            std::vector<LocalIndexType> outsideScvfs;
-            for (LocalIndexType j = i+1; j < numNodalScvfs; ++j)
-            {
-                // a face that has been handled already cannot be an "outside" face here
-                if (!isHandled[j])
-                {
-                    const auto scvIndices2 = [&nodalIndexSet, j] ()
-                                             {
-                                                 auto tmp = nodalIndexSet.neighboringScvIndices(j);
-                                                 std::sort(tmp.begin(), tmp.end());
-                                                 return tmp;
-                                             } ();
-
-                    // if the sizes aren't equal, this can't be an "outside" face
-                    if (scvIndices2.size() != numNeighborsI)
-                        continue;
-
-                    // if the neighboring scv indices are the same, this is an "outside" face
-                    if (std::equal(scvIndices.begin(), scvIndices.end(), scvIndices2.begin()))
-                        outsideScvfs.push_back(j);
-                }
-            }
-
-            // insert mappings
+            // insert a new iv-local face
             const auto curIvLocalScvfIdx = ivToNodeScvf_.size();
             nodeToIvScvf_[i] = curIvLocalScvfIdx;
             isHandled[i] = true;
-            for (const auto nodeLocalScvfIdx : outsideScvfs)
+
+            // construct local index sets
+            const auto& flipScvfIndices = flipScvfIndexSet[nodalIndexSet.scvfIdxGlobal(i)];
+            const auto numFlipIndices = flipScvfIndices.size();
+
+            ScvfNeighborLocalIndexSet neighborsLocal;
+            neighborsLocal.resize(numFlipIndices + 1);
+            neighborsLocal[0] = nodalIndexSet.insideScvIdxLocal(i);
+
+            // mappings for all flip scvf
+            for (unsigned int j = 0; j < numFlipIndices; ++j)
             {
-                nodeToIvScvf_[nodeLocalScvfIdx] = curIvLocalScvfIdx;
-                isHandled[nodeLocalScvfIdx] = true;
+                const auto outsideScvfIdx = flipScvfIndices[j];
+                for (unsigned int nodeLocalIdx = 0; nodeLocalIdx < nodalIndexSet.numScvfs(); ++nodeLocalIdx)
+                {
+                    if (nodalIndexSet.scvfIdxGlobal(nodeLocalIdx) == outsideScvfIdx)
+                    {
+                        neighborsLocal[j+1] = nodalIndexSet.insideScvIdxLocal(nodeLocalIdx);
+                        nodeToIvScvf_[nodeLocalIdx] = curIvLocalScvfIdx;
+                        isHandled[nodeLocalIdx] = true;
+                        break; // go to next outside scvf
+                    }
+                }
             }
+
+            scvfNeighborScvLocalIndices_.push_back(neighborsLocal);
             ivToNodeScvf_.push_back(i);
             numFaces_++;
         }
-
-        // compute local neighboring scv indices for the iv-local scvfs
-        scvfNeighborScvLocalIndices_.resize(numFaces_);
-        for (unsigned int i = 0; i < numFaces_; ++i)
-        {
-            const auto& neighborsGlobal = nodalIndexSet_.neighboringScvIndices(ivToNodeScvf_[i]);
-            const auto numNeighbors = nodalIndexSet_.scvfIsOnBoundary(ivToNodeScvf_[i]) ? 1 : neighborsGlobal.size();
-
-            scvfNeighborScvLocalIndices_[i].resize(numNeighbors);
-            for (unsigned int j = 0; j < numNeighbors; ++j)
-                scvfNeighborScvLocalIndices_[i][j] = findLocalScvIdx_(neighborsGlobal[j]);
-        }
     }
 
     //! returns the corresponding nodal index set
     const DualGridNodalIndexSet& nodalIndexSet() const { return nodalIndexSet_; }
 
-    //! returns a global scvf idx for a given iv_local scvf index
-    GridIndexType scvfIdxGlobal(LocalIndexType ivLocalScvfIdx) const
-    { return nodalIndexSet_.scvfIdxGlobal( ivToNodeScvf_[ivLocalScvfIdx] ); }
-
-    //! returns the iv-local scvf idx of the i-th scvfs embedded in a local scv
-    LocalIndexType scvfIdxLocal(LocalIndexType scvIdxLocal, unsigned int i) const
-    { return nodeToIvScvf_[ nodalIndexSet_.scvfIdxLocal(scvIdxLocal, i) ]; }
+    //! returns the global scv indices connected to this dual grid node
+    const GridStencilType& globalScvIndices() const { return nodalIndexSet_.globalScvIndices(); }
 
-    //! returns the local indices of the neighboring scvs of an scvf
-    const ScvfNeighborLocalIndexSet& neighboringLocalScvIndices(LocalIndexType ivLocalScvfIdx) const
-    { return scvfNeighborScvLocalIndices_[ivLocalScvfIdx]; }
+    //! returns the global scvf indices connected to this dual grid node
+    const GridScvfStencilType& globalScvfIndices() const { return nodalIndexSet_.globalScvfIndices(); }
 
     //! returns the number of faces in the interaction volume
     std::size_t numFaces() const { return numFaces_; }
@@ -164,11 +135,26 @@ public:
     //! returns the number of scvs in the interaction volume
     std::size_t numScvs() const { return nodalIndexSet_.numScvs(); }
 
-    //! returns the global scv indices connected to this dual grid node
-    const GridStencilType& globalScvIndices() const { return nodalIndexSet_.globalScvIndices(); }
+    //! returns a global scvf idx for a given iv-local scvf index
+    GridIndexType scvfIdxGlobal(LocalIndexType ivLocalScvfIdx) const
+    {
+        assert(ivLocalScvfIdx < numFaces());
+        return nodalIndexSet_.scvfIdxGlobal( ivToNodeScvf_[ivLocalScvfIdx] );
+    }
 
-    //! returns the global scvf indices connected to this dual grid node
-    const GridScvfStencilType& globalScvfIndices() const { return nodalIndexSet_.globalScvfIndices(); }
+    //! returns the iv-local scvf idx of the i-th scvf embedded in a local scv
+    LocalIndexType scvfIdxLocal(LocalIndexType scvIdxLocal, unsigned int i) const
+    {
+        assert(nodalIndexSet_.scvfIdxLocal(scvIdxLocal, i) < nodeToIvScvf_.size());
+        return nodeToIvScvf_[ nodalIndexSet_.scvfIdxLocal(scvIdxLocal, i) ];
+    }
+
+    //! returns the local indices of the neighboring scvs of an scvf
+    const ScvfNeighborLocalIndexSet& neighboringLocalScvIndices(LocalIndexType ivLocalScvfIdx) const
+    {
+        assert(ivLocalScvfIdx < numFaces());
+        return scvfNeighborScvLocalIndices_[ivLocalScvfIdx];
+    }
 
 private:
     //! returns the local scv index to a given global scv index