diff --git a/dumux/discretization/cellcentered/mpfa/helper.hh b/dumux/discretization/cellcentered/mpfa/helper.hh
index bc5a38b1ea486a452c70c7219ce86df30add9ef4..23fe038d1acbd21db62ddd97bf80e91a38c45050 100644
--- a/dumux/discretization/cellcentered/mpfa/helper.hh
+++ b/dumux/discretization/cellcentered/mpfa/helper.hh
@@ -75,7 +75,7 @@ public:
 
         unsigned int count = 0;
         ScvfVector scvfVector({nullptr});
-        for (const auto& scvf : scvfs(fvGeometry))
+        for (auto&& scvf : scvfs(fvGeometry))
         {
             if (scvf.vertexIndex() == vIdxGlobal)
             {
@@ -473,7 +473,7 @@ public:
         auto elementCenter = element.geometry().center();
 
         LocalIndexType count = 0;
-        for (const auto& scvf : scvfs(fvGeometry))
+        for (auto&& scvf : scvfs(fvGeometry))
         {
             if (scvf.vertexIndex() == vIdxGlobal)
             {
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh b/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh
index 34e5699952acf4ce3afb117137ef10935e84b679..f11cb30f2e2313d3415a4c7fdf39015c8d955bc3 100644
--- a/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/helper.hh
@@ -83,10 +83,6 @@ public:
 
         fillEntitySeeds_(scvSeeds, outerScvSeeds, globalScvfIndices, problem, element, fvGeometry, scvf);
 
-        // free unused memory
-        scvSeeds.shrink_to_fit();
-        outerScvSeeds.shrink_to_fit();
-
         // return interaction volume seed
         return InteractionVolumeSeed(std::move(scvSeeds), std::move(outerScvSeeds), std::move(globalScvfIndices));
     }
@@ -139,19 +135,21 @@ private:
         // make the first scv seed, we know this element will NOT be on the lowest local level
         auto scvfVector = Implementation::getScvFacesAtVertex(scvf.vertexIndex(), element, fvGeometry);
         auto localScvfIdx = Implementation::getLocalFaceIndex(scvf, scvfVector);
-        auto globalScvIdx = scvf.insideScvIdx();
         scvSeeds.emplace_back( GlobalIndexSet({scvfVector[0]->index(), scvfVector[1]->index()}),
-                               globalScvIdx,
+                               scvf.insideScvIdx(),
                                localScvfIdx );
 
         // get the surrounding elements and "outside" data
-        LocalIndexType otherScvfIdx = localScvfIdx == 1 ? 0 : 1;
+        LocalIndexType otherScvfIdx = 1-localScvfIdx;
         auto e2 = problem.model().globalFvGeometry().element(scvf.outsideScvIdx());
         auto e3 = problem.model().globalFvGeometry().element(scvfVector[otherScvfIdx]->outsideScvIdx());
+
         auto e2Geometry = localView(problem.model().globalFvGeometry());
         auto e3Geometry = localView(problem.model().globalFvGeometry());
+
         e2Geometry.bindElement(e2);
         e3Geometry.bindElement(e3);
+
         auto e2Scvfs = Implementation::getCommonAndNextScvFace(scvf, e2Geometry, /*clockwise?*/localScvfIdx == 1);
         auto e3Scvfs = Implementation::getCommonAndNextScvFace(*scvfVector[otherScvfIdx], e3Geometry, /*clockwise?*/localScvfIdx == 0);
 
@@ -165,16 +163,17 @@ private:
                                otherScvfIdx );
 
         // Outer seed for e3, we know the local common scvf index will be localScvfIdx in 2d
-        const auto& f3 = *e3Scvfs[localScvfIdx];
-        outerScvSeeds.emplace_back(f3.insideScvIdx(), f3.index());
+        outerScvSeeds.emplace_back(e3Scvfs[localScvfIdx]->insideScvIdx(),
+                                   e3Scvfs[localScvfIdx]->index());
 
         // Outer seed for outside of e2, we know the local scvf index here will be localScvfIdx in 2d
         auto e4 = problem.model().globalFvGeometry().element(e2Scvfs[localScvfIdx]->outsideScvIdx());
         auto e4Geometry = localView(problem.model().globalFvGeometry());
         e4Geometry.bindElement(e4);
         auto e4Scvfs = Implementation::getCommonAndNextScvFace(*e2Scvfs[localScvfIdx], e4Geometry, /*clockwise?*/localScvfIdx == 1);
-        const auto& f4 = *e4Scvfs[otherScvfIdx];
-        outerScvSeeds.emplace_back(f4.insideScvIdx(), f4.index());
+
+        outerScvSeeds.emplace_back(e4Scvfs[otherScvfIdx]->insideScvIdx(),
+                                   e4Scvfs[otherScvfIdx]->index());
     }
 };
 
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh
index 47c2967a6f6c1ab19262d0220ce665d0aee34bd8..5cf49681045d604dcfd5ef06777e22918b01304b 100644
--- a/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolume.hh
@@ -131,11 +131,10 @@ public:
             std::array<Matrix, 2> M;
             M[0] = assembleMatrix_(getTensor, interactionRegions_[0]);
             M[1] = assembleMatrix_(getTensor, interactionRegions_[1]);
-            const auto id = MpfaHelper::selectionCriterion(interactionRegions_[0], interactionRegions_[1], M[0], M[1]);
 
-            const auto& ir = interactionRegions_[id];
+            auto id = MpfaHelper::selectionCriterion(interactionRegions_[0], interactionRegions_[1], M[0], M[1]);
+            postSolve_(interactionRegions_[id]);
             T_ = M[id];
-            postSolve_(ir);
         }
     }
 
@@ -259,9 +258,9 @@ private:
         if (regionUnique_)
         {
             interactionRegions_.reserve(1);
-            const auto& scvSeed1 = seed.scvSeed(0);
-            const auto& scvSeed2 = seed.outerScvSeed(0);
-            const auto& scvSeed3 = seed.outerScvSeed(1);
+            auto&& scvSeed1 = seed.scvSeed(0);
+            auto&& scvSeed2 = seed.outerScvSeed(0);
+            auto&& scvSeed3 = seed.outerScvSeed(1);
             auto e1 = globalFvGeometry.element(scvSeed1.globalIndex());
             auto e2 = globalFvGeometry.element(scvSeed2.globalIndex());
             auto e3 = globalFvGeometry.element(scvSeed3.globalIndex());
@@ -271,10 +270,10 @@ private:
         else
         {
             interactionRegions_.reserve(2);
-            const auto& scvSeed1 = seed.scvSeed(0);
-            const auto& scvSeed2 = seed.scvSeed(1);
-            const auto& outerScvSeed1 = seed.outerScvSeed(0);
-            const auto& outerScvSeed2 = seed.outerScvSeed(1);
+            auto&& scvSeed1 = seed.scvSeed(0);
+            auto&& scvSeed2 = seed.scvSeed(1);
+            auto&& outerScvSeed1 = seed.outerScvSeed(0);
+            auto&& outerScvSeed2 = seed.outerScvSeed(1);
             auto e1 = globalFvGeometry.element(scvSeed1.globalIndex());
             auto e2 = globalFvGeometry.element(scvSeed2.globalIndex());
             auto e3 = globalFvGeometry.element(outerScvSeed1.globalIndex());
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolumeseed.hh b/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolumeseed.hh
index 808abbece40e3feaf9b6c02aca04a6fbb26eb59f..6579142df86878b1aa28f938bc19a2682067a9f4 100644
--- a/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolumeseed.hh
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/interactionvolumeseed.hh
@@ -75,9 +75,12 @@ public:
         for (auto&& localScvSeed : outerScvSeeds())
             globalIndices.push_back(localScvSeed.globalIndex());
 
-        // make the entries unique
-        std::sort(globalIndices.begin(), globalIndices.end());
-        globalIndices.erase(std::unique(globalIndices.begin(), globalIndices.end()), globalIndices.end());
+        // make the entries unique if interaction region is not uniquely defined
+        if (!isUnique())
+        {
+            std::sort(globalIndices.begin(), globalIndices.end());
+            globalIndices.erase(std::unique(globalIndices.begin(), globalIndices.end()), globalIndices.end());
+        }
 
         return globalIndices;
     }
diff --git a/dumux/discretization/cellcentered/mpfa/lmethod/localsubcontrolentityseeds.hh b/dumux/discretization/cellcentered/mpfa/lmethod/localsubcontrolentityseeds.hh
index 739ed82e5fc17e2514a17334203d230ea91ecd78..054323f6df8f5b00c75b2ed94b153fbe7eb55e6d 100644
--- a/dumux/discretization/cellcentered/mpfa/lmethod/localsubcontrolentityseeds.hh
+++ b/dumux/discretization/cellcentered/mpfa/lmethod/localsubcontrolentityseeds.hh
@@ -82,10 +82,10 @@ public:
       scvfIndexGlobal_(scvSeed.globalScvfIndices()[scvSeed.contiFaceLocalIdx()]) {}
 
 
-    const GlobalIndexType globalIndex() const
+    GlobalIndexType globalIndex() const
     { return scvIndexGlobal_; }
 
-    const GlobalIndexType globalScvfIndex() const
+    GlobalIndexType globalScvfIndex() const
     { return scvfIndexGlobal_; }
 
 private:
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/helper.hh b/dumux/discretization/cellcentered/mpfa/omethod/helper.hh
index 2e031e6ddbde4981dc11385c2d98724ec80e296a..12ae0c978118b4980f81a7b997219716f2a1378b 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/helper.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/helper.hh
@@ -159,8 +159,8 @@ private:
                                  const ScvfPointerVector& scvfVector,
                                  std::vector<ScvSeedType>& scvSeeds,
                                  std::vector<ScvfSeedType>& scvfSeeds,
-                                 const GlobalIndexType scvIdx0,
-                                 const bool clockWise = false)
+                                 GlobalIndexType scvIdx0,
+                                 bool clockWise = false)
     {
         // extract the actual local indices from the containers
         LocalIndexType localScvIdx = scvSeeds.size();
@@ -220,8 +220,8 @@ private:
             auto outsideScvfVector = Implementation::getCommonAndNextScvFace(curScvf, outsideFvGeometry, clockWise);
             GlobalIndexType commonFaceCoordIdx = clockWise ? 0 : 1;
             GlobalIndexType nextFaceCoordIdx = clockWise ? 1 : 0;
-            auto& commonScvf = *outsideScvfVector[commonFaceCoordIdx];
-            auto& nextScvf = *outsideScvfVector[nextFaceCoordIdx];
+            auto&& commonScvf = *outsideScvfVector[commonFaceCoordIdx];
+            auto&& nextScvf = *outsideScvfVector[nextFaceCoordIdx];
 
             // create local scv face entity of the current scvf
             LocalIndexType outsideLocalScvIdx = localScvIdx+1;
@@ -377,7 +377,7 @@ private:
 
         for (int coordDir = 0; coordDir < dim; ++coordDir)
         {
-            const auto& actualScvf = *scvfVector[coordDir];
+            auto&& actualScvf = *scvfVector[coordDir];
 
             // if scvf is on a boundary, we create the scvfSeed and make no neighbor
             if (actualScvf.boundary())
@@ -405,8 +405,8 @@ private:
                     // find scvf in outside corresponding to the actual scvf
                     auto outsideScvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, outsideElement, outsideFvGeometry);
                     auto commonFaceLocalIdx = Implementation::getCommonFaceLocalIndex(actualScvf, outsideScvfVector);
-                    const auto& outsideScvf = *outsideScvfVector[commonFaceLocalIdx];
-                    const auto outsideScvfIdx = outsideScvf.index();
+                    auto&& outsideScvf = *outsideScvfVector[commonFaceLocalIdx];
+                    auto outsideScvfIdx = outsideScvf.index();
 
                     // check if the outside scv already exists and get its local index
                     bool outsideScvExists = false;
@@ -415,8 +415,7 @@ private:
                     {
                         if (scvSeed.globalIndex() == outsideGlobalScvIdx)
                         {
-                            outsideScvExists = true;
-                            break;
+                            outsideScvExists = true; break;
                         }
                         // keep track of local index
                         outsideLocalScvIdx++;
@@ -657,8 +656,7 @@ private:
                 {
                     if (scvSeed.globalIndex() == outsideGlobalScvIdx)
                     {
-                        outsideExists = true;
-                        break;
+                        outsideExists = true; break;
                     }
                     // keep track of local index
                     outsideLocalScvIdx++;
@@ -678,7 +676,7 @@ private:
                     // find scvf in outside corresponding to the actual scvf
                     auto outsideScvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, outsideElement, outsideFvGeometry);
                     auto commonFaceLocalIdx = Implementation::getCommonFaceLocalIndex(actualScvf, outsideScvfVector);
-                    const auto& outsideScvf = *outsideScvfVector[commonFaceLocalIdx];
+                    auto&& outsideScvf = *outsideScvfVector[commonFaceLocalIdx];
 
                     // create scvf seed
                     scvfSeeds.emplace_back(actualScvf,
@@ -725,7 +723,7 @@ private:
                         // find scvf in outside corresponding to the actual scvf
                         auto outsideScvfVector = Implementation::getScvFacesAtVertex(vIdxGlobal, outsideElement, outsideFvGeometry);
                         auto commonFaceLocalIdx = Implementation::getCommonFaceLocalIndex(actualScvf, outsideScvfVector);
-                        const auto& outsideScvf = *outsideScvfVector[commonFaceLocalIdx];
+                        auto&& outsideScvf = *outsideScvfVector[commonFaceLocalIdx];
 
                         // make scv face seed
                         scvfSeeds.emplace_back(actualScvf,
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
index 36911bed391ed36c8b0c0e3db312ee3716662c33..5449a566333eaa0e7e829a1e9e52a39ed1847e0e 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
@@ -142,9 +142,9 @@ public:
         createLocalEntities_(seed);
 
         // reserve memory
-        const auto numLocalScvs = localScvs_.size();
-        const auto numLocalScvfs = localScvfs_.size();
-        const auto maxNumVolVars = numLocalScvs + numLocalScvfs;
+        auto numLocalScvs = localScvs_.size();
+        auto numLocalScvfs = localScvfs_.size();
+        auto maxNumVolVars = numLocalScvs + numLocalScvfs;
         volVarsStencil_ = seed.globalScvIndices(); // boundary vol vars are placed at the end
         volVarsStencil_.reserve(maxNumVolVars);
         volVarsPositions_.reserve(maxNumVolVars);
@@ -152,12 +152,12 @@ public:
         fluxFaceIndexSet_.reserve(numLocalScvfs);
 
         // the positions where the vol vars are defined at (required for the gravitational acceleration)
-        for (const auto& localScv : localScvs_)
+        for (auto&& localScv : localScvs_)
             volVarsPositions_.push_back(localScv.center());
 
         // eventually add dirichlet vol var indices and set up local index sets of flux and dirichlet faces
         LocalIndexType localScvfIdx = 0;
-        for (const auto& localScvf : localScvfs_)
+        for (auto&& localScvf : localScvfs_)
         {
             // eventually add vol var index and corresponding position
             if (localScvf.faceType() == MpfaFaceTypes::dirichlet)
@@ -203,20 +203,20 @@ public:
         T_ += D;
     }
 
-    void assembleNeumannFluxes(const unsigned int eqIdx)
+    void assembleNeumannFluxes(unsigned int eqIdx)
     {
         if (!onBoundary() || GET_PROP_VALUE(TypeTag, UseTpfaBoundary))
             return;
 
         LocalIndexType fluxFaceIdx = 0;
-        for (const auto localFluxFaceIdx : fluxFaceIndexSet_)
+        for (auto localFluxFaceIdx : fluxFaceIndexSet_)
         {
-            const auto& localScvf = localScvf_(localFluxFaceIdx);
+            auto&& localScvf = localScvf_(localFluxFaceIdx);
             const auto faceType = localScvf.faceType();
             if (faceType == MpfaFaceTypes::neumann)
             {
-                const auto& element = localElement_(localScvf.insideLocalScvIndex());
-                const auto& globalScvf = fvGeometry_().scvf(localScvf.insideGlobalScvfIndex());
+                auto&& element = localElement_(localScvf.insideLocalScvIndex());
+                auto&& globalScvf = fvGeometry_().scvf(localScvf.insideGlobalScvfIndex());
                 auto neumannFlux = problem_().neumann(element, this->fvGeometry_(), this->elemVolVars_(), globalScvf)[eqIdx];
                 neumannFlux *= globalScvf.area();
 
@@ -237,7 +237,7 @@ public:
         auto scvfGlobalIdx = scvf.index();
 
         LocalIndexType localIdx = 0;
-        for (const auto& localScvf : localScvfs_)
+        for (auto&& localScvf : localScvfs_)
         {
             if (localScvf.insideGlobalScvfIndex() == scvfGlobalIdx)
                 return LocalFaceData(localIdx, localScvf.insideLocalScvIndex(), false);
@@ -257,10 +257,10 @@ public:
     typename std::enable_if< (d < dw), LocalFaceData >::type
     getLocalFaceData(const SubControlVolumeFace& scvf) const
     {
-        const auto scvfGlobalIdx = scvf.index();
+        auto scvfGlobalIdx = scvf.index();
 
         LocalIndexType localFaceIdx = 0;
-        for (const auto& localScvf : localScvfs_)
+        for (auto&& localScvf : localScvfs_)
         {
             if (localScvf.insideGlobalScvfIndex() == scvfGlobalIdx)
                 return LocalFaceData(localFaceIdx, localScvf.insideLocalScvIndex(), false);
@@ -271,7 +271,7 @@ public:
                     return LocalFaceData(localFaceIdx, localScvf.outsideLocalScvIndex(), true);
 
                 // Now we have to find wich local scv is the one the global scvf is embedded in
-                const auto scvGlobalIdx = scvf.insideScvIdx();
+                auto scvGlobalIdx = scvf.insideScvIdx();
                 for (auto localScvIdx : localScvf.outsideLocalScvIndices())
                     if (localScv_(localScvIdx).globalIndex() == scvGlobalIdx)
                         return LocalFaceData(localFaceIdx, localScvIdx, true);
@@ -311,16 +311,16 @@ public:
         DynamicVector tij(volVarsStencil().size(), 0.0);
 
         // get the local scv and iterate over local coordinates
-        const std::size_t numLocalScvs = localScvs_.size();
-        const auto& localScv = localScv_(localFaceData.localScvIndex);
-        const auto& localScvf = localScvf_(localFaceData.localScvfIndex);
+        std::size_t numLocalScvs = localScvs_.size();
+        auto&& localScv = localScv_(localFaceData.localScvIndex);
+        auto&& localScvf = localScvf_(localFaceData.localScvfIndex);
 
         auto idxInOutside = this->findLocalIndex(localScvf.outsideLocalScvIndices(), localFaceData.localScvIndex);
-        const auto& wijk = wijk_[localFaceData.localScvfIndex][idxInOutside+1];
+        auto&& wijk = wijk_[localFaceData.localScvfIndex][idxInOutside+1];
         for (int localDir = 0; localDir < dim; localDir++)
         {
             auto localScvfIdx = localScv.localScvfIndex(localDir);
-            const auto& localScvf = localScvf_(localScvfIdx);
+            auto&& localScvf = localScvf_(localScvfIdx);
 
             // add entries from the face unknowns
             if (localScvf.faceType() != MpfaFaceTypes::dirichlet)
@@ -393,17 +393,17 @@ private:
         localScvs_.reserve(numScvs);
         localScvfs_.reserve(numScvfs);
 
-        for (const auto& scvSeed : seed.scvSeeds())
+        for (auto&& scvSeed : seed.scvSeeds())
         {
             auto element = problem_().model().globalFvGeometry().element(scvSeed.globalIndex());
             localScvs_.emplace_back(LocalScvType(problem_(), element, fvGeometry_(), scvSeed));
             localElements_.emplace_back(std::move(element));
         }
 
-        for (const auto& scvfSeed : seed.scvfSeeds())
+        for (auto&& scvfSeed : seed.scvfSeeds())
         {
             // we have to use the "inside" scv face here
-            const auto& scvf = fvGeometry_().scvf(scvfSeed.insideGlobalScvfIndex());
+            auto&& scvf = fvGeometry_().scvf(scvfSeed.insideGlobalScvfIndex());
             localScvfs_.emplace_back(LocalScvfType(scvfSeed, scvf));
         }
     }
@@ -422,16 +422,16 @@ private:
 
         // loop over the local faces
         LocalIndexType rowIdx = 0;
-        for (const auto& localScvf : localScvfs_)
+        for (auto&& localScvf : localScvfs_)
         {
             auto faceType = localScvf.faceType();
             bool hasUnknown = faceType != MpfaFaceTypes::dirichlet;
             LocalIndexType idxInFluxFaces = hasUnknown ? this->findLocalIndex(fluxScvfIndexSet_(), rowIdx) : -1;
 
             // get diffusion tensor in "positive" sub volume
-            const auto posLocalScvIdx = localScvf.insideLocalScvIndex();
-            const auto& posLocalScv = localScv_(posLocalScvIdx);
-            const auto& posGlobalScv = fvGeometry_().scv(posLocalScv.globalIndex());
+            auto posLocalScvIdx = localScvf.insideLocalScvIndex();
+            auto&& posLocalScv = localScv_(posLocalScvIdx);
+            auto&& posGlobalScv = fvGeometry_().scv(posLocalScv.globalIndex());
             auto element = localElement_(posLocalScvIdx);
             auto tensor = getTensor(element, elemVolVars_()[posGlobalScv], posGlobalScv);
 
@@ -443,7 +443,7 @@ private:
             for (int localDir = 0; localDir < dim; localDir++)
             {
                 auto curLocalScvfIdx = posLocalScv.localScvfIndex(localDir);
-                const auto& curLocalScvf = localScvf_(curLocalScvfIdx);
+                auto&& curLocalScvf = localScvf_(curLocalScvfIdx);
 
                 // First, add the entries associated with face pressures (unkown or dirichlet)
                 if (curLocalScvf.faceType() != MpfaFaceTypes::dirichlet)
@@ -482,8 +482,8 @@ private:
                 unsigned int indexInOutsideData = 0;
                 for (auto negLocalScvIdx : localScvf.outsideLocalScvIndices())
                 {
-                    const auto& negLocalScv = localScv_(negLocalScvIdx);
-                    const auto& negGlobalScv = fvGeometry_().scv(negLocalScv.globalIndex());
+                    auto&& negLocalScv = localScv_(negLocalScvIdx);
+                    auto&& negGlobalScv = fvGeometry_().scv(negLocalScv.globalIndex());
                     auto negElement = localElement_(negLocalScvIdx);
                     auto negTensor = getTensor(negElement, elemVolVars_()[negGlobalScv], negGlobalScv);
 
@@ -494,7 +494,7 @@ private:
                     if (dim < dimWorld)
                     {
                         // outside scvf
-                        auto&& outsideScvf = fvGeometry_().scvf(localScvf.outsideGlobalScvfIndices()[indexInOutsideData]);
+                        auto&& outsideScvf = fvGeometry_().scvf(localScvf.outsideGlobalScvfIndex(indexInOutsideData));
                         auto negNormal = outsideScvf.unitOuterNormal();
                         negNormal *= -1.0;
                         negWijk = calculateOmegas_(negLocalScv, negNormal, localScvf.area(), negTensor);
@@ -509,7 +509,7 @@ private:
                     for (int localDir = 0; localDir < dim; localDir++)
                     {
                         auto curLocalScvfIdx = negLocalScv.localScvfIndex(localDir);
-                        const auto& curLocalScvf = localScvf_(curLocalScvfIdx);
+                        auto&& curLocalScvf = localScvf_(curLocalScvfIdx);
 
                         if (curLocalScvf.faceType() != MpfaFaceTypes::dirichlet)
                         {
@@ -558,12 +558,12 @@ private:
 
         // Loop over all the faces, in this case these are all dirichlet boundaries
         LocalIndexType rowIdx = 0;
-        for (const auto& localScvf : localScvfs_)
+        for (auto&& localScvf : localScvfs_)
         {
             // get diffusion tensor in "positive" sub volume
-            const auto posLocalScvIdx = localScvf.insideLocalScvIndex();
-            const auto& posLocalScv = localScv_(posLocalScvIdx);
-            const auto& posGlobalScv = fvGeometry_().scv(posLocalScv.globalIndex());
+            auto posLocalScvIdx = localScvf.insideLocalScvIndex();
+            auto&& posLocalScv = localScv_(posLocalScvIdx);
+            auto&& posGlobalScv = fvGeometry_().scv(posLocalScv.globalIndex());
             auto element = localElement_(posLocalScvIdx);
             auto tensor = getTensor(element, elemVolVars_()[posGlobalScv], posGlobalScv);
 
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeseed.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeseed.hh
index bdcdc831fc42678eb57b8b7eca62ac20b8689857..3a4c5ccf35092ed127a68546801ad554267d6902 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeseed.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeseed.hh
@@ -82,7 +82,6 @@ public:
                 globalIndices.push_back(scvfIdxGlobal);
         }
 
-        globalIndices.shrink_to_fit();
         return globalIndices;
     }
 
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh
index ce4046ae6e6afebaff9ebc1fefc1d98f6078f497..c3d8b4c98e5785857375389fe67f20b9d1249c96 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh
@@ -141,41 +141,32 @@ public:
       area_(scvf.area())
     {}
 
-    const GlobalIndexType insideGlobalScvfIndex() const
+    GlobalIndexType insideGlobalScvfIndex() const
     { return scvfSeed_().insideGlobalScvfIndex(); }
 
+    GlobalIndexType insideGlobalScvIndex() const
+    { return scvfSeed_().insideGlobalScvIndex(); }
+
+    LocalIndexType insideLocalScvIndex() const
+    { return scvfSeed_().insideLocalScvIndex(); }
+
     const GlobalIndexSet& outsideGlobalScvfIndices() const
     { return scvfSeed_().outsideGlobalScvfIndices(); }
 
-    const GlobalIndexType outsideGlobalScvfIndex() const
-    {
-        assert(scvfSeed_().outsideGlobalScvfIndices().size() == 1 && "outside scvf index not uniquely defined");
-        return scvfSeed_().outsideGlobalScvfIndices()[0];
-    }
-
-    const LocalIndexType insideLocalScvIndex() const
-    { return scvfSeed_().insideLocalScvIndex(); }
+    const GlobalIndexSet& outsideGlobalScvIndices() const
+    { return scvfSeed_().outsideGlobalScvIndices(); }
 
     const LocalIndexSet& outsideLocalScvIndices() const
     { return scvfSeed_().outsideLocalScvIndices(); }
 
-    const LocalIndexType outsideLocalScvIndex() const
-    {
-        assert(scvfSeed_().outsideLocalScvIndices().size() == 1 && "outside local scv index not uniquely defined");
-        return scvfSeed_().outsideLocalScvIndices()[0];
-    }
+    GlobalIndexType outsideGlobalScvfIndex(unsigned int outsideIdx = 0) const
+    { return scvfSeed_().outsideGlobalScvfIndex(outsideIdx); }
 
-    const GlobalIndexType insideGlobalScvIndex() const
-    { return scvfSeed_().insideGlobalScvIndex(); }
+    GlobalIndexType outsideGlobalScvIndex(unsigned int outsideIdx = 0) const
+    { return scvfSeed_().outsideGlobalScvIndex(outsideIdx); }
 
-    const GlobalIndexSet& outsideGlobalScvIndices() const
-    { return scvfSeed_().outsideGlobalScvIndices(); }
-
-    const GlobalIndexType outsideGlobalScvIndex() const
-    {
-        assert(scvfSeed_().outsideGlobalScvIndices().size() == 1 && "outside scv index not uniquely defined");
-        return scvfSeed_().outsideGlobalScvIndices()[0];
-    }
+    LocalIndexType outsideLocalScvIndex(unsigned int outsideIdx = 0) const
+    { return scvfSeed_().outsideLocalScvIndex(outsideIdx); }
 
     MpfaFaceTypes faceType() const
     { return scvfSeed_().faceType(); }
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh
index 7f7cbba105739a833dd3cc323c489db72d930378..db5884a5ec3595fcd713d55d2cd725ab5af228c8 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentityseeds.hh
@@ -123,30 +123,33 @@ public:
         outsideGlobalScvfIndices_.reserve(size);
       }
 
-    const GlobalIndexType insideGlobalScvfIndex() const
+    GlobalIndexType insideGlobalScvfIndex() const
     { return insideScvfGlobalIdx_; }
 
-    const GlobalIndexSet& outsideGlobalScvfIndices() const
-    { return outsideGlobalScvfIndices_; }
-
-    const GlobalIndexType outsideGlobalScvfIndex() const
-    {
-        assert(outsideGlobalScvfIndices_.size() == 1 && "outside global scvf index not uniquely defined");
-        return outsideGlobalScvfIndices_[0];
-    }
+    GlobalIndexType insideGlobalScvIndex() const
+    { return insideScvGlobalIdx_; }
 
-    const LocalIndexType insideLocalScvIndex() const
+    LocalIndexType insideLocalScvIndex() const
     { return insideScvLocalIdx_; }
 
-    const LocalIndexSet& outsideLocalScvIndices() const
-    { return outsideLocalScvIndices_; }
+    GlobalIndexType outsideGlobalScvfIndex(unsigned int outsideIdx = 0) const
+    { return outsideGlobalScvfIndices_[outsideIdx]; }
 
-    const GlobalIndexType insideGlobalScvIndex() const
-    { return insideScvGlobalIdx_; }
+    GlobalIndexType outsideGlobalScvIndex(unsigned int outsideIdx = 0) const
+    { return outsideGlobalScvIndices_[outsideIdx]; }
+
+    LocalIndexType outsideLocalScvIndex(unsigned int outsideIdx = 0) const
+    { return outsideLocalScvIndices_[outsideIdx]; }
+
+    const GlobalIndexSet& outsideGlobalScvfIndices() const
+    { return outsideGlobalScvfIndices_; }
 
     const GlobalIndexSet& outsideGlobalScvIndices() const
     { return outsideGlobalScvIndices_; }
 
+    const LocalIndexSet& outsideLocalScvIndices() const
+    { return outsideLocalScvIndices_; }
+
     MpfaFaceTypes faceType() const
     { return faceType_; }
 
diff --git a/dumux/discretization/cellcentered/mpfa/omethodfps/helper.hh b/dumux/discretization/cellcentered/mpfa/omethodfps/helper.hh
index d0e6d3faa15dfac346e93162a9d7ee1b0aca231d..4ea6e4619d0e05582745d07360141c146a9c6d31 100644
--- a/dumux/discretization/cellcentered/mpfa/omethodfps/helper.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethodfps/helper.hh
@@ -30,10 +30,11 @@ namespace Dumux
 /*!
  * \ingroup Mpfa
  * \brief Helper class to get the required information on an interaction volume.
- *        Specialization for the Mpfa-O method in two dimensions.
+ *        We simply forward to the mpfa-o method here.
  */
 template<class TypeTag, int dim, int dimWorld>
-class MpfaMethodHelper<TypeTag, MpfaMethods::oMethodFps, dim, dimWorld> : public MpfaMethodHelper<TypeTag, MpfaMethods::oMethod, dim, dimWorld>
+class MpfaMethodHelper<TypeTag, MpfaMethods::oMethodFps, dim, dimWorld>
+: public MpfaMethodHelper<TypeTag, MpfaMethods::oMethod, dim, dimWorld>
 {};
 
 } // end namespace
diff --git a/dumux/discretization/cellcentered/mpfa/omethodfps/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethodfps/interactionvolume.hh
index bf0e537c826e0d426bec9591b76c194c9e8f326f..2262c42f330e1fcc66647d51a19fa8f226c0f977 100644
--- a/dumux/discretization/cellcentered/mpfa/omethodfps/interactionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethodfps/interactionvolume.hh
@@ -147,7 +147,7 @@ public:
         this->T_ += mc.BF;
     }
 
-    void assembleNeumannFluxes(const unsigned int eqIdx)
+    void assembleNeumannFluxes(unsigned int eqIdx)
     {
         ParentType::assembleNeumannFluxes(eqIdx);
 
@@ -156,7 +156,7 @@ public:
 
         auto& neumannFluxes = this->neumannFluxes_;
         neumannFluxes.resize(this->fluxScvfIndexSet_().size() + 1);
-        neumannFluxes[neumannFluxes.size()-1] = std::accumulate(neumannFluxes.begin(), neumannFluxes.end(), 0.0);
+        neumannFluxes[neumannFluxes.size()-1] = std::accumulate(neumannFluxes.begin(), neumannFluxes.end()-1, 0.0);
     }
 
 private:
@@ -184,15 +184,15 @@ private:
     template<typename GetTensorFunction>
     void assemblePositiveScv(const GetTensorFunction& getTensor,
                              const LocalScvfType& localScvf,
-                             const LocalIndexType localScvfIdx,
+                             LocalIndexType localScvfIdx,
                              LocalMatrixContainer& mc,
-                             const bool boundary = false)
+                             bool boundary = false)
     {
         // get diffusion tensor in "positive" sub volume
-        const auto localScvIdx = localScvf.insideLocalScvIndex();
-        const auto& localScv = this->localScv_(localScvIdx);
-        const auto& globalScv = this->fvGeometry_().scv(localScv.globalIndex());
-        const auto& element = this->localElement_(localScvIdx);
+        auto localScvIdx = localScvf.insideLocalScvIndex();
+        auto&& localScv = this->localScv_(localScvIdx);
+        auto&& globalScv = this->fvGeometry_().scv(localScv.globalIndex());
+        auto&& element = this->localElement_(localScvIdx);
         auto D = makeTensor_(getTensor(element, this->elemVolVars_()[globalScv], globalScv));
 
         // the local finite element basis
@@ -204,7 +204,7 @@ private:
         auto ipLocal = localScv.geometry().local(localScvf.ip());
 
         // find normal coordinate direction and integration point for divergence condition
-        LocalIndexType divEqNormalDir = normalDir == 1 ? 0 : 1;
+        LocalIndexType divEqNormalDir = 1 - normalDir;
         LocalPosition divEqIpLocal(0.0);
         divEqIpLocal[divEqNormalDir] = divEqNormalDir == 1 ? c_ : 1.0 - (1.0-c_)*p_;
         divEqIpLocal[normalDir] = divEqNormalDir == 1 ? c_ + (1.0-c_)*p_ : c_;
@@ -230,15 +230,15 @@ private:
     template<typename GetTensorFunction>
     void assembleNegativeScv(const GetTensorFunction& getTensor,
                              const LocalScvfType& localScvf,
-                             const LocalIndexType localScvfIdx,
+                             LocalIndexType localScvfIdx,
                              LocalMatrixContainer& mc)
     {
         // get diffusion tensor in "negative" sub volume
         for (auto localScvIdx : localScvf.outsideLocalScvIndices())
         {
-            const auto& localScv = this->localScv_(localScvIdx);
-            const auto& globalScv = this->fvGeometry_().scv(localScv.globalIndex());
-            const auto& element = this->localElement_(localScvIdx);;
+            auto&& localScv = this->localScv_(localScvIdx);
+            auto&& globalScv = this->fvGeometry_().scv(localScv.globalIndex());
+            auto&& element = this->localElement_(localScvIdx);;
             auto D = makeTensor_(getTensor(element, this->elemVolVars_()[globalScv], globalScv));
 
             // the local finite element bases of the scvs
@@ -250,7 +250,7 @@ private:
             auto ipLocal = localScv.geometry().local(localScvf.ip());
 
             // find normals and integration points in the two scvs for condition of zero divergence
-            LocalIndexType divEqNormalDir = normalDir == 1 ? 0 : 1;
+            LocalIndexType divEqNormalDir = 1 - normalDir;
             LocalPosition divEqIpLocal(0.0);
             divEqIpLocal[divEqNormalDir] = divEqNormalDir == 1 ? c_ : 1.0 - (1.0-c_)*p_;
             divEqIpLocal[normalDir] = divEqNormalDir == 1 ? c_ + (1.0-c_)*p_ : c_;
@@ -269,12 +269,12 @@ private:
     void addFaceFluxCoefficients_(const LocalScvType& localScv,
                                   const FeLocalBasis& localBasis,
                                   const Tensor& D,
-                                  const LocalIndexType rowIdx,
+                                  LocalIndexType rowIdx,
                                   const LocalPosition& ipLocal,
-                                  const LocalIndexType normalDir,
+                                  LocalIndexType normalDir,
                                   LocalMatrixContainer& mc,
-                                  const bool isFluxEq,
-                                  const bool isRHS = false)
+                                  bool isFluxEq,
+                                  bool isRHS = false)
     {
         // In case we're on a flux continuity face, get local index
         LocalIndexType eqSystemIdx = isFluxEq ? this->findLocalIndex(this->fluxScvfIndexSet_(), rowIdx) : -1;
@@ -307,7 +307,6 @@ private:
         for (int localDir = 0; localDir < dim; localDir++)
         {
             auto localScvfIdx = localScv.localScvfIndex(localDir);
-
             if (this->localScvf_(localScvfIdx).faceType() != MpfaFaceTypes::dirichlet)
             {
                 Scalar aij = factor*(localD[normalDir]*shapeJacobian[localDir+1][0]);
@@ -334,9 +333,9 @@ private:
                                      const FeLocalBasis& localBasis,
                                      const Tensor& D,
                                      const LocalPosition& ipLocal,
-                                     const LocalIndexType normalDir,
+                                     LocalIndexType normalDir,
                                      LocalMatrixContainer& mc,
-                                     const bool isBoundary = false)
+                                     bool isBoundary = false)
     {
         // fluxes on the auxiliary volume have to be scaled
         static const Scalar auxArea = 1.0 - c_;
diff --git a/dumux/discretization/cellcentered/mpfa/omethodfps/localsubcontrolentities.hh b/dumux/discretization/cellcentered/mpfa/omethodfps/localsubcontrolentities.hh
index 17057339ac5f5d57818e113d7cf83457fd4fb279..616d9466cbcc6f97fd4e85600ca1e6fe02c61361 100644
--- a/dumux/discretization/cellcentered/mpfa/omethodfps/localsubcontrolentities.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethodfps/localsubcontrolentities.hh
@@ -70,13 +70,13 @@ public:
     GlobalIndexType globalIndex() const
     { return scvSeed_().globalIndex(); }
 
-    GlobalIndexType localScvfIndex(const LocalIndexType coordDir) const
+    GlobalIndexType localScvfIndex(LocalIndexType coordDir) const
     {
         assert(coordDir < dim);
         return scvSeed_().localScvfIndices()[coordDir];
     }
 
-    LocalIndexType getScvfIdxInScv(const LocalIndexType localScvfIndex) const
+    LocalIndexType getScvfIdxInScv(LocalIndexType localScvfIndex) const
     {
         auto it = std::find(scvSeed_().localScvfIndices().begin(), scvSeed_().localScvfIndices().end(), localScvfIndex);
         assert(it != scvSeed_().localScvfIndices().end() && "Could not find the local coordinate of the local scvf");
diff --git a/dumux/porousmediumflow/implicit/fluxvariablescache.hh b/dumux/porousmediumflow/implicit/fluxvariablescache.hh
index ee9551e7d86b1394d33fd3be33ab76a69d07cd79..3d6efb86709ba3ac2298bcf816297d617f3ac0c0 100644
--- a/dumux/porousmediumflow/implicit/fluxvariablescache.hh
+++ b/dumux/porousmediumflow/implicit/fluxvariablescache.hh
@@ -32,8 +32,7 @@ namespace Dumux
 {
 // forward declaration
 template<class TypeTag, DiscretizationMethods Method>
-class PorousMediumFluxVariablesCacheImplementation
-{};
+class PorousMediumFluxVariablesCacheImplementation;
 
 namespace Properties
 {
@@ -164,7 +163,7 @@ private:
 
 // forward declaration of the base class of the mpfa flux variables cache
 template<class TypeTag, bool EnableAdvection, bool EnableMolecularDiffusion, bool EnableEnergyBalance>
-class MpfaPorousMediumFluxVariablesCache {};
+class MpfaPorousMediumFluxVariablesCache;
 
 //! Base class for the advective cache in mpfa methods
 template<class TypeTag>