diff --git a/dumux/freeflow/stokes/model.hh b/dumux/freeflow/stokes/model.hh
index c430a35301217c682278f1074284371f6337fb1f..215ac1c2ae897c7ff012990539ba92f0e199efa2 100644
--- a/dumux/freeflow/stokes/model.hh
+++ b/dumux/freeflow/stokes/model.hh
@@ -104,11 +104,8 @@ public:
         ElementVolumeVariables elemVolVars;
 
         // Loop over elements
-        for (const auto& element : elements(this->problem_.gridView()))
+        for (const auto& element : elements(this->problem_.gridView(), Dune::Partitions::interior))
         {
-            if (element.partitionType() != Dune::InteriorEntity)
-                continue;
-
             fvGeometry.update(this->gridView_(), element);
             elemVolVars.update(this->problem_(), element, fvGeometry);
             this->localResidual().evalFluxes(element, elemVolVars);
diff --git a/dumux/geomechanics/el1p2c/model.hh b/dumux/geomechanics/el1p2c/model.hh
index 9998cc7fcf6b44105f1c8b9df69f3a2fd9b23078..a04dbe03a79f5dcb57671af90d24a64239c39d2d 100644
--- a/dumux/geomechanics/el1p2c/model.hh
+++ b/dumux/geomechanics/el1p2c/model.hh
@@ -246,132 +246,129 @@ public:
 
         // initialize start and end of element iterator
         // loop over all elements (cells)
-        for (const auto& element : elements(this->gridView_()))
+        for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
         {
-            if(element.partitionType() == Dune::InteriorEntity)
-            {
-                unsigned int eIdx = this->problem_().model().elementMapper().index(element);
-                rank[eIdx] = this->gridView_().comm().rank();
+            unsigned int eIdx = this->problem_().model().elementMapper().index(element);
+            rank[eIdx] = this->gridView_().comm().rank();
 
-                fvGeometry.update(this->gridView_(), element);
-                elemBcTypes.update(this->problem_(), element, fvGeometry);
-                elemVolVars.update(this->problem_(), element, fvGeometry, false);
+            fvGeometry.update(this->gridView_(), element);
+            elemBcTypes.update(this->problem_(), element, fvGeometry);
+            elemVolVars.update(this->problem_(), element, fvGeometry, false);
 
-                // loop over all local vertices of the cell
-                int numScv = element.subEntities(dim);
+            // loop over all local vertices of the cell
+            int numScv = element.subEntities(dim);
 
-                for (int scvIdx = 0; scvIdx < numScv; ++scvIdx)
-                {
-                    unsigned int vIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dim);
-
-                    pressure[vIdxGlobal] = elemVolVars[scvIdx].pressure();
-                    moleFraction0[vIdxGlobal] = elemVolVars[scvIdx].moleFraction(0);
-                    moleFraction1[vIdxGlobal] = elemVolVars[scvIdx].moleFraction(1);
-                    massFraction0[vIdxGlobal] = elemVolVars[scvIdx].massFraction(0);
-                    massFraction1[vIdxGlobal] = elemVolVars[scvIdx].massFraction(1);
-                    // in case of rock mechanics sign convention solid displacement is
-                    // defined to be negative if it points in positive coordinate direction
-                    if(rockMechanicsSignConvention_){
-                        DimVector tmpDispl;
-                        tmpDispl = Scalar(0);
-                        tmpDispl -= elemVolVars[scvIdx].displacement();
-                        displacement[vIdxGlobal] = tmpDispl;
-                        }
-
-                    else
-                        displacement[vIdxGlobal] = elemVolVars[scvIdx].displacement();
-
-                    density[vIdxGlobal] = elemVolVars[scvIdx].density();
-                    viscosity[vIdxGlobal] = elemVolVars[scvIdx].viscosity();
-                    porosity[vIdxGlobal] = elemVolVars[scvIdx].porosity();
-                    Kx[vIdxGlobal] =    this->problem_().spatialParams().intrinsicPermeability(
-                                    element, fvGeometry, scvIdx)[0][0];
-                    // calculate cell quantities by adding up scv quantities and dividing through numScv
-                    cellPorosity[eIdx] += elemVolVars[scvIdx].porosity()    / numScv;
-                    cellKx[eIdx] += this->problem_().spatialParams().intrinsicPermeability(
-                                    element, fvGeometry, scvIdx)[0][0] / numScv;
-                    cellPressure[eIdx] += elemVolVars[scvIdx].pressure()    / numScv;
-                };
-
-                // calculate cell quantities for variables which are defined at the integration point
-                Scalar tmpEffPoro;
-                DimMatrix tmpEffStress;
-                tmpEffStress = Scalar(0);
-                tmpEffPoro = Scalar(0);
-
-                // loop over all scv-faces of the cell
-                for (int fIdx = 0; fIdx < fvGeometry.numScvf; fIdx++) {
-
-                    //prepare the flux calculations (set up and prepare geometry, FE gradients)
-                    FluxVariables fluxVars;
-                    fluxVars.update(this->problem_(),
-                                    element, fvGeometry,
-                                    fIdx,
-                                    elemVolVars);
-
-                    // divide by number of scv-faces and sum up edge values
-                    tmpEffPoro = fluxVars.effPorosity() / fvGeometry.numScvf;
-                    tmpEffStress = fluxVars.sigma();
-                    tmpEffStress /= fvGeometry.numScvf;
-
-                    effPorosity[eIdx] += tmpEffPoro;
-
-                    // in case of rock mechanics sign convention compressive stresses
-                    // are defined to be positive
-                    if(rockMechanicsSignConvention_){
-                        effStressX[eIdx] -= tmpEffStress[0];
-                        if (dim >= 2) {
-                            effStressY[eIdx] -= tmpEffStress[1];
-                        }
-                        if (dim >= 3) {
-                            effStressZ[eIdx] -= tmpEffStress[2];
-                        }
-                    }
-                    else{
-                        effStressX[eIdx] += tmpEffStress[0];
-                        if (dim >= 2) {
-                            effStressY[eIdx] += tmpEffStress[1];
-                        }
-                        if (dim >= 3) {
-                            effStressZ[eIdx] += tmpEffStress[2];
-                        }
+            for (int scvIdx = 0; scvIdx < numScv; ++scvIdx)
+            {
+                unsigned int vIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dim);
+
+                pressure[vIdxGlobal] = elemVolVars[scvIdx].pressure();
+                moleFraction0[vIdxGlobal] = elemVolVars[scvIdx].moleFraction(0);
+                moleFraction1[vIdxGlobal] = elemVolVars[scvIdx].moleFraction(1);
+                massFraction0[vIdxGlobal] = elemVolVars[scvIdx].massFraction(0);
+                massFraction1[vIdxGlobal] = elemVolVars[scvIdx].massFraction(1);
+                // in case of rock mechanics sign convention solid displacement is
+                // defined to be negative if it points in positive coordinate direction
+                if(rockMechanicsSignConvention_){
+                    DimVector tmpDispl;
+                    tmpDispl = Scalar(0);
+                    tmpDispl -= elemVolVars[scvIdx].displacement();
+                    displacement[vIdxGlobal] = tmpDispl;
                     }
-                }
 
-                // calculate total stresses
+                else
+                    displacement[vIdxGlobal] = elemVolVars[scvIdx].displacement();
+
+                density[vIdxGlobal] = elemVolVars[scvIdx].density();
+                viscosity[vIdxGlobal] = elemVolVars[scvIdx].viscosity();
+                porosity[vIdxGlobal] = elemVolVars[scvIdx].porosity();
+                Kx[vIdxGlobal] =    this->problem_().spatialParams().intrinsicPermeability(
+                                element, fvGeometry, scvIdx)[0][0];
+                // calculate cell quantities by adding up scv quantities and dividing through numScv
+                cellPorosity[eIdx] += elemVolVars[scvIdx].porosity()    / numScv;
+                cellKx[eIdx] += this->problem_().spatialParams().intrinsicPermeability(
+                                element, fvGeometry, scvIdx)[0][0] / numScv;
+                cellPressure[eIdx] += elemVolVars[scvIdx].pressure()    / numScv;
+            };
+
+            // calculate cell quantities for variables which are defined at the integration point
+            Scalar tmpEffPoro;
+            DimMatrix tmpEffStress;
+            tmpEffStress = Scalar(0);
+            tmpEffPoro = Scalar(0);
+
+            // loop over all scv-faces of the cell
+            for (int fIdx = 0; fIdx < fvGeometry.numScvf; fIdx++) {
+
+                //prepare the flux calculations (set up and prepare geometry, FE gradients)
+                FluxVariables fluxVars;
+                fluxVars.update(this->problem_(),
+                                element, fvGeometry,
+                                fIdx,
+                                elemVolVars);
+
+                // divide by number of scv-faces and sum up edge values
+                tmpEffPoro = fluxVars.effPorosity() / fvGeometry.numScvf;
+                tmpEffStress = fluxVars.sigma();
+                tmpEffStress /= fvGeometry.numScvf;
+
+                effPorosity[eIdx] += tmpEffPoro;
+
                 // in case of rock mechanics sign convention compressive stresses
-                // are defined to be positive and total stress is calculated by adding the pore pressure
+                // are defined to be positive
                 if(rockMechanicsSignConvention_){
-                    totalStressX[eIdx][0] = effStressX[eIdx][0]    + cellPressure[eIdx];
-                    totalStressX[eIdx][1] = effStressX[eIdx][1];
-                    totalStressX[eIdx][2] = effStressX[eIdx][2];
+                    effStressX[eIdx] -= tmpEffStress[0];
                     if (dim >= 2) {
-                        totalStressY[eIdx][0] = effStressY[eIdx][0];
-                        totalStressY[eIdx][1] = effStressY[eIdx][1]    + cellPressure[eIdx];
-                        totalStressY[eIdx][2] = effStressY[eIdx][2];
+                        effStressY[eIdx] -= tmpEffStress[1];
                     }
                     if (dim >= 3) {
-                        totalStressZ[eIdx][0] = effStressZ[eIdx][0];
-                        totalStressZ[eIdx][1] = effStressZ[eIdx][1];
-                        totalStressZ[eIdx][2] = effStressZ[eIdx][2]    + cellPressure[eIdx];
+                        effStressZ[eIdx] -= tmpEffStress[2];
                     }
                 }
                 else{
-                    totalStressX[eIdx][0] = effStressX[eIdx][0]    - cellPressure[eIdx];
-                    totalStressX[eIdx][1] = effStressX[eIdx][1];
-                    totalStressX[eIdx][2] = effStressX[eIdx][2];
+                    effStressX[eIdx] += tmpEffStress[0];
                     if (dim >= 2) {
-                        totalStressY[eIdx][0] = effStressY[eIdx][0];
-                        totalStressY[eIdx][1] = effStressY[eIdx][1]    - cellPressure[eIdx];
-                        totalStressY[eIdx][2] = effStressY[eIdx][2];
+                        effStressY[eIdx] += tmpEffStress[1];
                     }
                     if (dim >= 3) {
-                        totalStressZ[eIdx][0] = effStressZ[eIdx][0];
-                        totalStressZ[eIdx][1] = effStressZ[eIdx][1];
-                        totalStressZ[eIdx][2] = effStressZ[eIdx][2]    - cellPressure[eIdx];
+                        effStressZ[eIdx] += tmpEffStress[2];
                     }
                 }
             }
+
+            // calculate total stresses
+            // in case of rock mechanics sign convention compressive stresses
+            // are defined to be positive and total stress is calculated by adding the pore pressure
+            if(rockMechanicsSignConvention_){
+                totalStressX[eIdx][0] = effStressX[eIdx][0]    + cellPressure[eIdx];
+                totalStressX[eIdx][1] = effStressX[eIdx][1];
+                totalStressX[eIdx][2] = effStressX[eIdx][2];
+                if (dim >= 2) {
+                    totalStressY[eIdx][0] = effStressY[eIdx][0];
+                    totalStressY[eIdx][1] = effStressY[eIdx][1]    + cellPressure[eIdx];
+                    totalStressY[eIdx][2] = effStressY[eIdx][2];
+                }
+                if (dim >= 3) {
+                    totalStressZ[eIdx][0] = effStressZ[eIdx][0];
+                    totalStressZ[eIdx][1] = effStressZ[eIdx][1];
+                    totalStressZ[eIdx][2] = effStressZ[eIdx][2]    + cellPressure[eIdx];
+                }
+            }
+            else{
+                totalStressX[eIdx][0] = effStressX[eIdx][0]    - cellPressure[eIdx];
+                totalStressX[eIdx][1] = effStressX[eIdx][1];
+                totalStressX[eIdx][2] = effStressX[eIdx][2];
+                if (dim >= 2) {
+                    totalStressY[eIdx][0] = effStressY[eIdx][0];
+                    totalStressY[eIdx][1] = effStressY[eIdx][1]    - cellPressure[eIdx];
+                    totalStressY[eIdx][2] = effStressY[eIdx][2];
+                }
+                if (dim >= 3) {
+                    totalStressZ[eIdx][0] = effStressZ[eIdx][0];
+                    totalStressZ[eIdx][1] = effStressZ[eIdx][1];
+                    totalStressZ[eIdx][2] = effStressZ[eIdx][2]    - cellPressure[eIdx];
+                }
+            }
         }
 
         // calculate principal stresses i.e. the eigenvalues of the total stress tensor
diff --git a/dumux/geomechanics/el2p/basemodel.hh b/dumux/geomechanics/el2p/basemodel.hh
index 76b1719632deb813a9df2db08e4d9d2cb0326505..bb2850a4987fcbdd4ab58a9e4b0521fd1fa4fae7 100644
--- a/dumux/geomechanics/el2p/basemodel.hh
+++ b/dumux/geomechanics/el2p/basemodel.hh
@@ -242,21 +242,18 @@ public:
     {
         storage = 0;
 
-        for (const auto& element : elements(gridView_()))
+        for (const auto& element : elements(gridView_(), Dune::Partitions::interior))
         {
-            if(element.partitionType() == Dune::InteriorEntity)
-            {
-                localResidual().evalStorage(element);
+            localResidual().evalStorage(element);
 
-                if (isBox)
-                {
-                    for (int i = 0; i < element.subEntities(dim); ++i)
-                        storage += localResidual().storageTerm()[i];
-                }
-                else
-                {
-                    storage += localResidual().storageTerm()[0];
-                }
+            if (isBox)
+            {
+                for (int i = 0; i < element.subEntities(dim); ++i)
+                    storage += localResidual().storageTerm()[i];
+            }
+            else
+            {
+                storage += localResidual().storageTerm()[0];
             }
         }
 
diff --git a/dumux/geomechanics/el2p/model.hh b/dumux/geomechanics/el2p/model.hh
index 7cca7f6ca6e49ce8968d8812ebd2e82a1b1d47cf..9837b85c788c0739eebf01c1297403f5469cfa3b 100644
--- a/dumux/geomechanics/el2p/model.hh
+++ b/dumux/geomechanics/el2p/model.hh
@@ -326,10 +326,8 @@ public:
         const typename GridFunctionSpace::Ordering& ordering = gridFunctionSpace.ordering();
         // initialize start and end of element iterator
         // loop over all elements (cells)
-        for (const auto& element : elements(this->gridView_())) {
-            if(element.partitionType() == Dune::InteriorEntity)
-            {
-
+        for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
+        {
             // get FE function spaces to calculate gradients (gradient data of momentum balance
             // equation is not stored in fluxvars since it is not evaluated at box integration point)
             // copy the values of the sol vector to the localFunctionSpace values of the current element
@@ -527,7 +525,6 @@ public:
                 }
             }
         }
-        }
 
         // calculate principal stresses i.e. the eigenvalues of the total stress tensor
         Scalar a1, a2, a3;
diff --git a/dumux/geomechanics/elastic/model.hh b/dumux/geomechanics/elastic/model.hh
index 973b01e34926c5c2e9bdae7655bd3aa4b103370a..f24a2b00936145eaeb3e67b536f645aa9e5e9039 100644
--- a/dumux/geomechanics/elastic/model.hh
+++ b/dumux/geomechanics/elastic/model.hh
@@ -118,10 +118,8 @@ public:
         VolumeVariables volVars;
         ElementBoundaryTypes elemBcTypes;
 
-        for (const auto& element : elements(this->gridView_()))
+        for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
         {
-            if(element.partitionType() == Dune::InteriorEntity)
-            {
             int eIdx = this->problem_().model().elementMapper().index(element);
             rank[eIdx] = this->gridView_().comm().rank();
 
@@ -183,7 +181,6 @@ public:
                     sigmaz[eIdx] += stress[2];
                 }
             }
-            }
         }
 
 
diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh
index c8bec031c46a065c6e859b5b8d208df22e046cb8..40b2d0cec6d2da4667bc65d77c0171a0e4789c3e 100644
--- a/dumux/implicit/model.hh
+++ b/dumux/implicit/model.hh
@@ -268,24 +268,21 @@ public:
     {
         storage = 0;
 
-        for (const auto& element : elements(gridView_()))
+        for (const auto& element : elements(gridView_(), Dune::Partitions::interior))
         {
-            if(element.partitionType() == Dune::InteriorEntity)
-            {
-                localResidual().evalStorage(element);
+            localResidual().evalStorage(element);
 
-                if (isBox)
-                {
-                    for (int i = 0; i < element.subEntities(dim); ++i)
-                    {
-                        storage += localResidual().storageTerm()[i];
-                    }
-                }
-                else
+            if (isBox)
+            {
+                for (int i = 0; i < element.subEntities(dim); ++i)
                 {
-                    storage += localResidual().storageTerm()[0];
+                    storage += localResidual().storageTerm()[i];
                 }
             }
+            else
+            {
+                storage += localResidual().storageTerm()[0];
+            }
         }
 
         if (gridView_().comm().size() > 1)
diff --git a/dumux/porousmediumflow/1p/implicit/model.hh b/dumux/porousmediumflow/1p/implicit/model.hh
index 55561343dcc8327f4eac20907730425632abe8e4..d07a500d62d0d65ade7e7aeffcde89ec5aed385b 100644
--- a/dumux/porousmediumflow/1p/implicit/model.hh
+++ b/dumux/porousmediumflow/1p/implicit/model.hh
@@ -100,33 +100,30 @@ public:
         unsigned numElements = this->gridView_().size(0);
         ScalarField *rank = writer.allocateManagedBuffer(numElements);
 
-        for (const auto& element : elements(this->gridView_()))
+        for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
         {
-            if(element.partitionType() == Dune::InteriorEntity)
+            int eIdx = this->problem_().model().elementMapper().index(element);
+            (*rank)[eIdx] = this->gridView_().comm().rank();
+
+            FVElementGeometry fvGeometry;
+            fvGeometry.update(this->gridView_(), element);
+
+            ElementVolumeVariables elemVolVars;
+            elemVolVars.update(this->problem_(),
+                               element,
+                               fvGeometry,
+                               false /* oldSol? */);
+
+            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+            {
+                int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
+                (*p)[dofIdxGlobal] = elemVolVars[scvIdx].pressure();
+            }
+
+            // velocity output
+            if (velocityOutput.enableOutput())
             {
-                int eIdx = this->problem_().model().elementMapper().index(element);
-                (*rank)[eIdx] = this->gridView_().comm().rank();
-
-                FVElementGeometry fvGeometry;
-                fvGeometry.update(this->gridView_(), element);
-
-                ElementVolumeVariables elemVolVars;
-                elemVolVars.update(this->problem_(),
-                                   element,
-                                   fvGeometry,
-                                   false /* oldSol? */);
-
-                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
-                {
-                    int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
-                    (*p)[dofIdxGlobal] = elemVolVars[scvIdx].pressure();
-                }
-
-                // velocity output
-                if (velocityOutput.enableOutput())
-                {
-                    velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, element, /*phaseIdx=*/0);
-                }
+                velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, element, /*phaseIdx=*/0);
             }
         }
 
diff --git a/dumux/porousmediumflow/1p2c/implicit/model.hh b/dumux/porousmediumflow/1p2c/implicit/model.hh
index edfa45fe0691ea83048af081b3d85066480c4db1..c772fbae84781bfaec8406010187016b0c719d4c 100644
--- a/dumux/porousmediumflow/1p2c/implicit/model.hh
+++ b/dumux/porousmediumflow/1p2c/implicit/model.hh
@@ -129,40 +129,37 @@ public:
         unsigned numElements = this->gridView_().size(0);
         ScalarField &rank = *writer.allocateManagedBuffer(numElements);
 
-        for (const auto& element : elements(this->gridView_()))
+        for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
         {
-            if(element.partitionType() == Dune::InteriorEntity)
+           int eIdx = this->problem_().model().elementMapper().index(element);
+
+            rank[eIdx] = this->gridView_().comm().rank();
+
+            FVElementGeometry fvGeometry;
+            fvGeometry.update(this->gridView_(), element);
+
+            ElementVolumeVariables elemVolVars;
+            elemVolVars.update(this->problem_(),
+                               element,
+                               fvGeometry,
+                               false /* oldSol? */);
+
+            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
             {
-                int eIdx = this->problem_().model().elementMapper().index(element);
-
-                rank[eIdx] = this->gridView_().comm().rank();
-
-                FVElementGeometry fvGeometry;
-                fvGeometry.update(this->gridView_(), element);
-
-                ElementVolumeVariables elemVolVars;
-                elemVolVars.update(this->problem_(),
-                                   element,
-                                   fvGeometry,
-                                   false /* oldSol? */);
-
-                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
-                {
-                    int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
-
-                    pressure[dofIdxGlobal] = elemVolVars[scvIdx].pressure();
-                    delp[dofIdxGlobal] = elemVolVars[scvIdx].pressure() - 1e5;
-                    moleFraction0[dofIdxGlobal] = elemVolVars[scvIdx].moleFraction(0);
-                    moleFraction1[dofIdxGlobal] = elemVolVars[scvIdx].moleFraction(1);
-                    massFraction0[dofIdxGlobal] = elemVolVars[scvIdx].massFraction(0);
-                    massFraction1[dofIdxGlobal] = elemVolVars[scvIdx].massFraction(1);
-                    rho[dofIdxGlobal] = elemVolVars[scvIdx].density();
-                    mu[dofIdxGlobal] = elemVolVars[scvIdx].viscosity();
-                }
-
-                // velocity output
-                velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, element, phaseIdx);
+                int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
+
+                pressure[dofIdxGlobal] = elemVolVars[scvIdx].pressure();
+                delp[dofIdxGlobal] = elemVolVars[scvIdx].pressure() - 1e5;
+                moleFraction0[dofIdxGlobal] = elemVolVars[scvIdx].moleFraction(0);
+                moleFraction1[dofIdxGlobal] = elemVolVars[scvIdx].moleFraction(1);
+                massFraction0[dofIdxGlobal] = elemVolVars[scvIdx].massFraction(0);
+                massFraction1[dofIdxGlobal] = elemVolVars[scvIdx].massFraction(1);
+                rho[dofIdxGlobal] = elemVolVars[scvIdx].density();
+                mu[dofIdxGlobal] = elemVolVars[scvIdx].viscosity();
             }
+
+            // velocity output
+            velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, element, phaseIdx);
         }
 
         writer.attachDofData(pressure, "P", isBox);
diff --git a/dumux/porousmediumflow/2p/implicit/model.hh b/dumux/porousmediumflow/2p/implicit/model.hh
index 024e032b46d5828d5c1a7fcf750df862d75386dc..5ea13e75476ec0f22cba24a769cf7ff2e334f09d 100644
--- a/dumux/porousmediumflow/2p/implicit/model.hh
+++ b/dumux/porousmediumflow/2p/implicit/model.hh
@@ -138,44 +138,41 @@ public:
         unsigned numElements = this->gridView_().size(0);
         ScalarField *rank = writer.allocateManagedBuffer(numElements);
 
-        for (const auto& element : elements(this->gridView_()))
+        for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
         {
-            if(element.partitionType() == Dune::InteriorEntity)
+            int eIdx = this->elementMapper().index(element);
+
+            (*rank)[eIdx] = this->gridView_().comm().rank();
+
+            FVElementGeometry fvGeometry;
+            fvGeometry.update(this->gridView_(), element);
+
+            ElementVolumeVariables elemVolVars;
+            elemVolVars.update(this->problem_(),
+                               element,
+                               fvGeometry,
+                               false /* oldSol? */);
+
+            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
             {
-                int eIdx = this->elementMapper().index(element);
-
-                (*rank)[eIdx] = this->gridView_().comm().rank();
-
-                FVElementGeometry fvGeometry;
-                fvGeometry.update(this->gridView_(), element);
-
-                ElementVolumeVariables elemVolVars;
-                elemVolVars.update(this->problem_(),
-                                   element,
-                                   fvGeometry,
-                                   false /* oldSol? */);
-
-                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
-                {
-                    int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
-
-                    (*pw)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(wPhaseIdx);
-                    (*pn)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(nPhaseIdx);
-                    (*pc)[dofIdxGlobal] = elemVolVars[scvIdx].capillaryPressure();
-                    (*sw)[dofIdxGlobal] = elemVolVars[scvIdx].saturation(wPhaseIdx);
-                    (*sn)[dofIdxGlobal] = elemVolVars[scvIdx].saturation(nPhaseIdx);
-                    (*rhoW)[dofIdxGlobal] = elemVolVars[scvIdx].density(wPhaseIdx);
-                    (*rhoN)[dofIdxGlobal] = elemVolVars[scvIdx].density(nPhaseIdx);
-                    (*mobW)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(wPhaseIdx);
-                    (*mobN)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(nPhaseIdx);
-                    (*poro)[dofIdxGlobal] = elemVolVars[scvIdx].porosity();
-                    (*Te)[dofIdxGlobal] = elemVolVars[scvIdx].temperature();
-                }
-
-                // velocity output
-                velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, element, wPhaseIdx);
-                velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, nPhaseIdx);
+                int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
+
+                (*pw)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(wPhaseIdx);
+                (*pn)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(nPhaseIdx);
+                (*pc)[dofIdxGlobal] = elemVolVars[scvIdx].capillaryPressure();
+                (*sw)[dofIdxGlobal] = elemVolVars[scvIdx].saturation(wPhaseIdx);
+                (*sn)[dofIdxGlobal] = elemVolVars[scvIdx].saturation(nPhaseIdx);
+                (*rhoW)[dofIdxGlobal] = elemVolVars[scvIdx].density(wPhaseIdx);
+                (*rhoN)[dofIdxGlobal] = elemVolVars[scvIdx].density(nPhaseIdx);
+                (*mobW)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(wPhaseIdx);
+                (*mobN)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(nPhaseIdx);
+                (*poro)[dofIdxGlobal] = elemVolVars[scvIdx].porosity();
+                (*Te)[dofIdxGlobal] = elemVolVars[scvIdx].temperature();
             }
+
+            // velocity output
+            velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, element, wPhaseIdx);
+            velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, nPhaseIdx);
         }
 
         writer.attachDofData(*sn, "Sn", isBox);
diff --git a/dumux/porousmediumflow/2p2c/implicit/model.hh b/dumux/porousmediumflow/2p2c/implicit/model.hh
index efbd00e4d23138d1e4001fa605c4dc7d99c343cd..3a1701f3e0376ec89bef52aeb384694cda981de2 100644
--- a/dumux/porousmediumflow/2p2c/implicit/model.hh
+++ b/dumux/porousmediumflow/2p2c/implicit/model.hh
@@ -201,16 +201,12 @@ public:
     {
         storage = 0;
 
-        for (const auto& element : elements(this->gridView_())) {
-            if(element.partitionType() == Dune::InteriorEntity)
-            {
-
-
-                this->localResidual().evalPhaseStorage(element, phaseIdx);
+        for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
+        {
+            this->localResidual().evalPhaseStorage(element, phaseIdx);
 
-                for (unsigned int i = 0; i < this->localResidual().storageTerm().size(); ++i)
-                    storage += this->localResidual().storageTerm()[i];
-            }
+            for (unsigned int i = 0; i < this->localResidual().storageTerm().size(); ++i)
+                storage += this->localResidual().storageTerm()[i];
         }
         if (this->gridView_().comm().size() > 1)
             storage = this->gridView_().comm().sum(storage);
@@ -321,62 +317,59 @@ public:
         unsigned numElements = this->gridView_().size(0);
         ScalarField *rank = writer.allocateManagedBuffer(numElements);
 
-        for (const auto& element : elements(this->gridView_()))
+        for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
         {
-            if(element.partitionType() == Dune::InteriorEntity)
-            {
-                int eIdx = this->elementMapper().index(element);
-                (*rank)[eIdx] = this->gridView_().comm().rank();
+            int eIdx = this->elementMapper().index(element);
+            (*rank)[eIdx] = this->gridView_().comm().rank();
 
-                FVElementGeometry fvGeometry;
-                fvGeometry.update(this->gridView_(), element);
+            FVElementGeometry fvGeometry;
+            fvGeometry.update(this->gridView_(), element);
 
-                ElementVolumeVariables elemVolVars;
-                elemVolVars.update(this->problem_(),
-                                   element,
-                                   fvGeometry,
-                                   false /* oldSol? */);
+            ElementVolumeVariables elemVolVars;
+            elemVolVars.update(this->problem_(),
+                               element,
+                               fvGeometry,
+                               false /* oldSol? */);
 
-                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
-                {
-                    int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
-
-                    (*sN)[dofIdxGlobal]    = elemVolVars[scvIdx].saturation(nPhaseIdx);
-                    (*sW)[dofIdxGlobal]    = elemVolVars[scvIdx].saturation(wPhaseIdx);
-                    (*pn)[dofIdxGlobal]    = elemVolVars[scvIdx].pressure(nPhaseIdx);
-                    (*pw)[dofIdxGlobal]    = elemVolVars[scvIdx].pressure(wPhaseIdx);
-                    (*pc)[dofIdxGlobal]    = elemVolVars[scvIdx].capillaryPressure();
-                    (*rhoW)[dofIdxGlobal]  = elemVolVars[scvIdx].density(wPhaseIdx);
-                    (*rhoN)[dofIdxGlobal]  = elemVolVars[scvIdx].density(nPhaseIdx);
-                    (*mobW)[dofIdxGlobal]  = elemVolVars[scvIdx].mobility(wPhaseIdx);
-                    (*mobN)[dofIdxGlobal]  = elemVolVars[scvIdx].mobility(nPhaseIdx);
-                    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-                        for (int compIdx = 0; compIdx < numComponents; ++compIdx)
-                        {
-                            (*massFrac[phaseIdx][compIdx])[dofIdxGlobal]
-                                = elemVolVars[scvIdx].massFraction(phaseIdx, compIdx);
-
-                            Valgrind::CheckDefined((*massFrac[phaseIdx][compIdx])[dofIdxGlobal][0]);
-                        }
-                    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-                        for (int compIdx = 0; compIdx < numComponents; ++compIdx)
-                        {
-                            (*moleFrac[phaseIdx][compIdx])[dofIdxGlobal]
-                                = elemVolVars[scvIdx].moleFraction(phaseIdx, compIdx);
+            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+            {
+                int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
+
+                (*sN)[dofIdxGlobal]    = elemVolVars[scvIdx].saturation(nPhaseIdx);
+                (*sW)[dofIdxGlobal]    = elemVolVars[scvIdx].saturation(wPhaseIdx);
+                (*pn)[dofIdxGlobal]    = elemVolVars[scvIdx].pressure(nPhaseIdx);
+                (*pw)[dofIdxGlobal]    = elemVolVars[scvIdx].pressure(wPhaseIdx);
+                (*pc)[dofIdxGlobal]    = elemVolVars[scvIdx].capillaryPressure();
+                (*rhoW)[dofIdxGlobal]  = elemVolVars[scvIdx].density(wPhaseIdx);
+                (*rhoN)[dofIdxGlobal]  = elemVolVars[scvIdx].density(nPhaseIdx);
+                (*mobW)[dofIdxGlobal]  = elemVolVars[scvIdx].mobility(wPhaseIdx);
+                (*mobN)[dofIdxGlobal]  = elemVolVars[scvIdx].mobility(nPhaseIdx);
+                for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                    for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+                    {
+                        (*massFrac[phaseIdx][compIdx])[dofIdxGlobal]
+                            = elemVolVars[scvIdx].massFraction(phaseIdx, compIdx);
 
-                            Valgrind::CheckDefined((*moleFrac[phaseIdx][compIdx])[dofIdxGlobal][0]);
-                        }
-                    (*poro)[dofIdxGlobal]  = elemVolVars[scvIdx].porosity();
-                    (*temperature)[dofIdxGlobal] = elemVolVars[scvIdx].temperature();
-                    (*phasePresence)[dofIdxGlobal]
-                        = staticDat_[dofIdxGlobal].phasePresence;
-                }
+                        Valgrind::CheckDefined((*massFrac[phaseIdx][compIdx])[dofIdxGlobal][0]);
+                    }
+                for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                    for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+                    {
+                        (*moleFrac[phaseIdx][compIdx])[dofIdxGlobal]
+                            = elemVolVars[scvIdx].moleFraction(phaseIdx, compIdx);
 
-                // velocity output
-                velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, element, wPhaseIdx);
-                velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, nPhaseIdx);
+                        Valgrind::CheckDefined((*moleFrac[phaseIdx][compIdx])[dofIdxGlobal][0]);
+                    }
+                (*poro)[dofIdxGlobal]  = elemVolVars[scvIdx].porosity();
+                (*temperature)[dofIdxGlobal] = elemVolVars[scvIdx].temperature();
+                (*phasePresence)[dofIdxGlobal]
+                    = staticDat_[dofIdxGlobal].phasePresence;
             }
 
+            // velocity output
+            velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, element, wPhaseIdx);
+            velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, nPhaseIdx);
+
         } // loop over elements
 
         writer.attachDofData(*sN,     "Sn", isBox);
diff --git a/dumux/porousmediumflow/2pnc/implicit/model.hh b/dumux/porousmediumflow/2pnc/implicit/model.hh
index 93838e82e81a648a55e6c2d560b34cf58178c5e7..2d05d5dfa1af84cc7ee9ee17cc6283d3bbcc8eb4 100644
--- a/dumux/porousmediumflow/2pnc/implicit/model.hh
+++ b/dumux/porousmediumflow/2pnc/implicit/model.hh
@@ -196,15 +196,12 @@ public:
     {
         storage = 0;
 
-        for (const auto& element : elements(this->gridView_()))
+        for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
         {
-            if(element.partitionType() == Dune::InteriorEntity)
-            {
-                this->localResidual().evalPhaseStorage(element, phaseIdx);
+            this->localResidual().evalPhaseStorage(element, phaseIdx);
 
-                for (unsigned int i = 0; i < this->localResidual().storageTerm().size(); ++i)
-                    storage += this->localResidual().storageTerm()[i];
-            }
+            for (unsigned int i = 0; i < this->localResidual().storageTerm().size(); ++i)
+                storage += this->localResidual().storageTerm()[i];
         }
 
         this->gridView_().comm().sum(storage);
diff --git a/dumux/porousmediumflow/3p/implicit/model.hh b/dumux/porousmediumflow/3p/implicit/model.hh
index 29975622d2715e98cfa62d6d63c926a7783a7117..1b25cd90fecfeebf594a14f953711e297cf8b8f2 100644
--- a/dumux/porousmediumflow/3p/implicit/model.hh
+++ b/dumux/porousmediumflow/3p/implicit/model.hh
@@ -135,43 +135,40 @@ public:
         unsigned numElements = this->gridView_().size(0);
         ScalarField *rank = writer.allocateManagedBuffer (numElements);
 
-        for (const auto& element : elements(this->gridView_()))
+        for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
         {
-            if(element.partitionType() == Dune::InteriorEntity)
-            {
-                int eIdx = this->problem_().elementMapper().index(element);
-                (*rank)[eIdx] = this->gridView_().comm().rank();
-
-                FVElementGeometry fvGeometry;
-                fvGeometry.update(this->gridView_(), element);
+            int eIdx = this->problem_().elementMapper().index(element);
+            (*rank)[eIdx] = this->gridView_().comm().rank();
 
+            FVElementGeometry fvGeometry;
+            fvGeometry.update(this->gridView_(), element);
 
-                ElementVolumeVariables elemVolVars;
-                elemVolVars.update(this->problem_(),
-                                   element,
-                                   fvGeometry,
-                                   false /* oldSol? */);
 
-                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
-                {
-                    int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
+            ElementVolumeVariables elemVolVars;
+            elemVolVars.update(this->problem_(),
+                               element,
+                               fvGeometry,
+                               false /* oldSol? */);
 
-                    for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
-                        (*saturation[phaseIdx])[dofIdxGlobal] = elemVolVars[scvIdx].saturation(phaseIdx);
-                        (*pressure[phaseIdx])[dofIdxGlobal] = elemVolVars[scvIdx].pressure(phaseIdx);
-                        (*density[phaseIdx])[dofIdxGlobal] = elemVolVars[scvIdx].density(phaseIdx);
-                    }
+            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+            {
+                int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
 
-                    (*poro)[dofIdxGlobal] = elemVolVars[scvIdx].porosity();
-                    (*perm)[dofIdxGlobal] = elemVolVars[scvIdx].permeability();
-                    (*temperature)[dofIdxGlobal] = elemVolVars[scvIdx].temperature();
+                for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+                    (*saturation[phaseIdx])[dofIdxGlobal] = elemVolVars[scvIdx].saturation(phaseIdx);
+                    (*pressure[phaseIdx])[dofIdxGlobal] = elemVolVars[scvIdx].pressure(phaseIdx);
+                    (*density[phaseIdx])[dofIdxGlobal] = elemVolVars[scvIdx].density(phaseIdx);
                 }
 
-                // velocity output
-                velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, element, wPhaseIdx);
-                velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, nPhaseIdx);
-                velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, gPhaseIdx);
+                (*poro)[dofIdxGlobal] = elemVolVars[scvIdx].porosity();
+                (*perm)[dofIdxGlobal] = elemVolVars[scvIdx].permeability();
+                (*temperature)[dofIdxGlobal] = elemVolVars[scvIdx].temperature();
             }
+
+            // velocity output
+            velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, element, wPhaseIdx);
+            velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, nPhaseIdx);
+            velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, gPhaseIdx);
         }
 
         writer.attachDofData(*saturation[wPhaseIdx], "sw", isBox);
diff --git a/dumux/porousmediumflow/3p3c/implicit/model.hh b/dumux/porousmediumflow/3p3c/implicit/model.hh
index ff78a61f8fce83f6273d7f0157f7a463799bdd42..daa976dcc510fc7b91385fafe18b67c6c0cbe286 100644
--- a/dumux/porousmediumflow/3p3c/implicit/model.hh
+++ b/dumux/porousmediumflow/3p3c/implicit/model.hh
@@ -201,14 +201,12 @@ public:
     {
         storage = 0;
 
-        for (const auto& element : elements(this->gridView_())) {
-            if(element.partitionType() == Dune::InteriorEntity)
-            {
-                this->localResidual().evalPhaseStorage(element, phaseIdx);
+        for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
+        {
+            this->localResidual().evalPhaseStorage(element, phaseIdx);
 
-                for (unsigned int i = 0; i < this->localResidual().storageTerm().size(); ++i)
-                    storage += this->localResidual().storageTerm()[i];
-            }
+            for (unsigned int i = 0; i < this->localResidual().storageTerm().size(); ++i)
+                storage += this->localResidual().storageTerm()[i];
         }
         if (this->gridView_().comm().size() > 1)
             storage = this->gridView_().comm().sum(storage);
@@ -321,53 +319,50 @@ public:
         unsigned numElements = this->gridView_().size(0);
         ScalarField *rank = writer.allocateManagedBuffer (numElements);
 
-        for (const auto& element : elements(this->gridView_()))
+        for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
         {
-            if(element.partitionType() == Dune::InteriorEntity)
-            {
-                int eIdx = this->problem_().elementMapper().index(element);
-                (*rank)[eIdx] = this->gridView_().comm().rank();
+            int eIdx = this->problem_().elementMapper().index(element);
+            (*rank)[eIdx] = this->gridView_().comm().rank();
 
-                FVElementGeometry fvGeometry;
-                fvGeometry.update(this->gridView_(), element);
+            FVElementGeometry fvGeometry;
+            fvGeometry.update(this->gridView_(), element);
 
 
-                ElementVolumeVariables elemVolVars;
-                elemVolVars.update(this->problem_(),
-                                   element,
-                                   fvGeometry,
-                                   false /* oldSol? */);
+            ElementVolumeVariables elemVolVars;
+            elemVolVars.update(this->problem_(),
+                               element,
+                               fvGeometry,
+                               false /* oldSol? */);
 
-                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
-                {
-                    int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
+            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+            {
+                int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
 
-                    for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
-                        (*saturation[phaseIdx])[dofIdxGlobal] = elemVolVars[scvIdx].saturation(phaseIdx);
-                        (*pressure[phaseIdx])[dofIdxGlobal] = elemVolVars[scvIdx].pressure(phaseIdx);
-                        (*density[phaseIdx])[dofIdxGlobal] = elemVolVars[scvIdx].density(phaseIdx);
-                    }
+                for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+                    (*saturation[phaseIdx])[dofIdxGlobal] = elemVolVars[scvIdx].saturation(phaseIdx);
+                    (*pressure[phaseIdx])[dofIdxGlobal] = elemVolVars[scvIdx].pressure(phaseIdx);
+                    (*density[phaseIdx])[dofIdxGlobal] = elemVolVars[scvIdx].density(phaseIdx);
+                }
 
-                    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-                        for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
-                            (*moleFraction[phaseIdx][compIdx])[dofIdxGlobal] =
-                                elemVolVars[scvIdx].moleFraction(phaseIdx,
-                                                                  compIdx);
+                for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+                    for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
+                        (*moleFraction[phaseIdx][compIdx])[dofIdxGlobal] =
+                            elemVolVars[scvIdx].moleFraction(phaseIdx,
+                                                              compIdx);
 
-                            Valgrind::CheckDefined((*moleFraction[phaseIdx][compIdx])[dofIdxGlobal]);
-                        }
+                        Valgrind::CheckDefined((*moleFraction[phaseIdx][compIdx])[dofIdxGlobal]);
                     }
-
-                    (*poro)[dofIdxGlobal] = elemVolVars[scvIdx].porosity();
-                    (*temperature)[dofIdxGlobal] = elemVolVars[scvIdx].temperature();
-                    (*phasePresence)[dofIdxGlobal] = staticDat_[dofIdxGlobal].phasePresence;
                 }
 
-                // velocity output
-                velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, element, wPhaseIdx);
-                velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, nPhaseIdx);
-                velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, gPhaseIdx);
+                (*poro)[dofIdxGlobal] = elemVolVars[scvIdx].porosity();
+                (*temperature)[dofIdxGlobal] = elemVolVars[scvIdx].temperature();
+                (*phasePresence)[dofIdxGlobal] = staticDat_[dofIdxGlobal].phasePresence;
             }
+
+            // velocity output
+            velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, element, wPhaseIdx);
+            velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, nPhaseIdx);
+            velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, gPhaseIdx);
         }
 
         writer.attachDofData(*saturation[wPhaseIdx], "Sw", isBox);
diff --git a/dumux/porousmediumflow/3pwateroil/model.hh b/dumux/porousmediumflow/3pwateroil/model.hh
index 11e374b31b6f04e5ad917d67848cc1bc1fcff119..a17c7ef2aa16a37cd98d37af1bcd9697376e942f 100644
--- a/dumux/porousmediumflow/3pwateroil/model.hh
+++ b/dumux/porousmediumflow/3pwateroil/model.hh
@@ -156,16 +156,14 @@ public:
 
         if (isBox)
         {
-            VertexIterator vIt = this->gridView_().template begin<dim> ();
-            const VertexIterator &vEndIt = this->gridView_().template end<dim> ();
-            for (; vIt != vEndIt; ++vIt)
+            for (const auto& vertex : vertices(this->gridView_()))
             {
-                int globalIdx = this->dofMapper().index(*vIt);
-                const GlobalPosition &globalPos = vIt->geometry().corner(0);
+                int globalIdx = this->dofMapper().index(vertex);
+                const GlobalPosition &globalPos = vertex.geometry().corner(0);
 
                 // initialize phase presence
                 staticDat_[globalIdx].phasePresence
-                    = this->problem_().initialPhasePresence(*vIt, globalIdx,
+                    = this->problem_().initialPhasePresence(vertex, globalIdx,
                                                         globalPos);
                 staticDat_[globalIdx].wasSwitched = false;
 
@@ -175,12 +173,10 @@ public:
         }
         else
         {
-            ElementIterator eIt = this->gridView_().template begin<0> ();
-            const ElementIterator &eEndIt = this->gridView_().template end<0> ();
-            for (; eIt != eEndIt; ++eIt)
+            for (const auto& element : elements(this->gridView_()))
             {
-                int globalIdx = this->dofMapper().index(*eIt);
-                const GlobalPosition &globalPos = eIt->geometry().center();
+                int globalIdx = this->dofMapper().index(element);
+                const GlobalPosition &globalPos = element.geometry().center();
 
                 // initialize phase presence
                 staticDat_[globalIdx].phasePresence
@@ -205,16 +201,12 @@ public:
     {
         storage = 0;
 
-        ElementIterator eIt = this->gridView_().template begin<0>();
-        const ElementIterator eEndIt = this->gridView_().template end<0>();
-        for (; eIt != eEndIt; ++eIt) {
-            if(eIt->partitionType() == Dune::InteriorEntity)
-            {
-                this->localResidual().evalPhaseStorage(*eIt, phaseIdx);
+        for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
+        {
+            this->localResidual().evalPhaseStorage(element, phaseIdx);
 
-                for (unsigned int i = 0; i < this->localResidual().storageTerm().size(); ++i)
-                    storage += this->localResidual().storageTerm()[i];
-            }
+            for (unsigned int i = 0; i < this->localResidual().storageTerm().size(); ++i)
+                storage += this->localResidual().storageTerm()[i];
         }
         if (this->gridView_().comm().size() > 1)
             storage = this->gridView_().comm().sum(storage);
@@ -335,26 +327,24 @@ public:
         unsigned numElements = this->gridView_().size(0);
         ScalarField *rank = writer.allocateManagedBuffer (numElements);
 
-        ElementIterator eIt = this->gridView_().template begin<0>();
-        ElementIterator eEndIt = this->gridView_().template end<0>();
-        for (; eIt != eEndIt; ++eIt)
+        for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
         {
-            int idx = this->problem_().elementMapper().index(*eIt);
+            int idx = this->problem_().elementMapper().index(element);
             (*rank)[idx] = this->gridView_().comm().rank();
 
             FVElementGeometry fvGeometry;
-            fvGeometry.update(this->gridView_(), *eIt);
+            fvGeometry.update(this->gridView_(), element);
 
 
             ElementVolumeVariables elemVolVars;
             elemVolVars.update(this->problem_(),
-                               *eIt,
+                               element,
                                fvGeometry,
                                false /* oldSol? */);
 
             for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
             {
-                int globalIdx = this->dofMapper().subIndex(*eIt, scvIdx, dofCodim);
+                int globalIdx = this->dofMapper().subIndex(element, scvIdx, dofCodim);
 
                 for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
                     (*saturation[phaseIdx])[globalIdx] = elemVolVars[scvIdx].fluidState().saturation(phaseIdx);
@@ -382,9 +372,9 @@ public:
             }
 
             // velocity output
-            velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, *eIt, wPhaseIdx);
-            velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, nPhaseIdx);
-            velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, gPhaseIdx);
+            velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, element, wPhaseIdx);
+            velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, nPhaseIdx);
+            velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, gPhaseIdx);
 
         }
 
@@ -493,14 +483,12 @@ public:
 
         FVElementGeometry fvGeometry;
         static VolumeVariables volVars;
-        ElementIterator eIt = this->gridView_().template begin<0> ();
-        const ElementIterator &eEndIt = this->gridView_().template end<0> ();
-        for (; eIt != eEndIt; ++eIt)
+        for (const auto& element : elements(this->gridView_()))
         {
-            fvGeometry.update(this->gridView_(), *eIt);
+            fvGeometry.update(this->gridView_(), element);
             for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
             {
-                int globalIdx = this->dofMapper().subIndex(*eIt, scvIdx, dofCodim);
+                int globalIdx = this->dofMapper().subIndex(element, scvIdx, dofCodim);
 
                 if (staticDat_[globalIdx].visited)
                     continue;
@@ -508,7 +496,7 @@ public:
                 staticDat_[globalIdx].visited = true;
                 volVars.update(curGlobalSol[globalIdx],
                                this->problem_(),
-                               *eIt,
+                               element,
                                fvGeometry,
                                scvIdx,
                                false);
diff --git a/dumux/porousmediumflow/mpnc/implicit/vtkwriter.hh b/dumux/porousmediumflow/mpnc/implicit/vtkwriter.hh
index 9fd84683c1b96710bb7539ace121fb153f5760ad..539abc38c3996e26085dec146d8a6fe9fe39ac01 100644
--- a/dumux/porousmediumflow/mpnc/implicit/vtkwriter.hh
+++ b/dumux/porousmediumflow/mpnc/implicit/vtkwriter.hh
@@ -75,38 +75,35 @@ public:
         ElementVolumeVariables elemVolVars;
         ElementBoundaryTypes elemBcTypes;
 
-        for (const auto& element : elements(problem_.gridView()))
+        for (const auto& element : elements(problem_.gridView(), Dune::Partitions::interior))
         {
-            if(element.partitionType() == Dune::InteriorEntity)
-            {
-                fvGeometry.update(problem_.gridView(), element);
-                elemBcTypes.update(problem_, element);
-                this->problem_.model().setHints(element, elemVolVars);
-                elemVolVars.update(problem_,
-                                   element,
-                                   fvGeometry,
-                                   false);
-                this->problem_.model().updateCurHints(element, elemVolVars);
+            fvGeometry.update(problem_.gridView(), element);
+            elemBcTypes.update(problem_, element);
+            this->problem_.model().setHints(element, elemVolVars);
+            elemVolVars.update(problem_,
+                               element,
+                               fvGeometry,
+                               false);
+            this->problem_.model().updateCurHints(element, elemVolVars);
 
-                // tell the sub-writers to do what ever they need to with
-                // their internal buffers when a given element is seen.
-                commonWriter_.processElement(element,
-                                             fvGeometry,
-                                             elemVolVars,
-                                             elemBcTypes);
-                massWriter_.processElement(element,
-                                           fvGeometry,
-                                           elemVolVars,
-                                           elemBcTypes);
-                energyWriter_.processElement(element,
-                                             fvGeometry,
-                                             elemVolVars,
-                                             elemBcTypes);
-                customWriter_.processElement(element,
-                                             fvGeometry,
-                                             elemVolVars,
-                                             elemBcTypes);
-            }
+            // tell the sub-writers to do what ever they need to with
+            // their internal buffers when a given element is seen.
+            commonWriter_.processElement(element,
+                                         fvGeometry,
+                                         elemVolVars,
+                                         elemBcTypes);
+            massWriter_.processElement(element,
+                                       fvGeometry,
+                                       elemVolVars,
+                                       elemBcTypes);
+            energyWriter_.processElement(element,
+                                         fvGeometry,
+                                         elemVolVars,
+                                         elemBcTypes);
+            customWriter_.processElement(element,
+                                         fvGeometry,
+                                         elemVolVars,
+                                         elemBcTypes);
         }
 
         // write everything to the output file
diff --git a/dumux/porousmediumflow/nonisothermal/implicit/model.hh b/dumux/porousmediumflow/nonisothermal/implicit/model.hh
index 1372b4eba0c8b67a71cbd9a08adf4d80c1846b33..ffb4af3cfe5d5e1286cde6d8e917f0ab389edfa3 100644
--- a/dumux/porousmediumflow/nonisothermal/implicit/model.hh
+++ b/dumux/porousmediumflow/nonisothermal/implicit/model.hh
@@ -95,24 +95,21 @@ public:
         typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
         ScalarField &temperature = *writer.allocateManagedBuffer(numDofs);
 
-        for (const auto& element : elements(this->gridView_()))
+        for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
         {
-            if(element.partitionType() == Dune::InteriorEntity)
-            {
-                FVElementGeometry fvGeometry;
-                fvGeometry.update(this->gridView_(), element);
+            FVElementGeometry fvGeometry;
+            fvGeometry.update(this->gridView_(), element);
 
-                ElementVolumeVariables elemVolVars;
-                elemVolVars.update(this->problem_(),
-                                   element,
-                                   fvGeometry,
-                                   false /* oldSol? */);
+            ElementVolumeVariables elemVolVars;
+            elemVolVars.update(this->problem_(),
+                               element,
+                               fvGeometry,
+                               false /* oldSol? */);
 
-                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
-                {
-                    int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
-                    temperature[dofIdxGlobal] = elemVolVars[scvIdx].temperature();
-                }
+            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+            {
+                int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
+                temperature[dofIdxGlobal] = elemVolVars[scvIdx].temperature();
             }
         }
 
diff --git a/dumux/porousmediumflow/richards/implicit/model.hh b/dumux/porousmediumflow/richards/implicit/model.hh
index 89a26369fbdc1a4afbc6c402e27d36efb8bf7085..1aebffb040b0cbaaf7b76a440980c9f0546ed4a5 100644
--- a/dumux/porousmediumflow/richards/implicit/model.hh
+++ b/dumux/porousmediumflow/richards/implicit/model.hh
@@ -174,48 +174,45 @@ public:
         unsigned numElements = this->gridView_().size(0);
         ScalarField *rank = writer.allocateManagedBuffer (numElements);
 
-        for (const auto& element : elements(this->gridView_()))
+        for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
         {
-            if(element.partitionType() == Dune::InteriorEntity)
+            int eIdx = this->problem_().model().elementMapper().index(element);
+
+            (*rank)[eIdx] = this->gridView_().comm().rank();
+
+            FVElementGeometry fvGeometry;
+            fvGeometry.update(this->gridView_(), element);
+
+            ElementVolumeVariables elemVolVars;
+            elemVolVars.update(this->problem_(),
+                               element,
+                               fvGeometry,
+                               false /* oldSol? */);
+
+            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
             {
-                int eIdx = this->problem_().model().elementMapper().index(element);
-
-                (*rank)[eIdx] = this->gridView_().comm().rank();
-
-                FVElementGeometry fvGeometry;
-                fvGeometry.update(this->gridView_(), element);
-
-                ElementVolumeVariables elemVolVars;
-                elemVolVars.update(this->problem_(),
-                                   element,
-                                   fvGeometry,
-                                   false /* oldSol? */);
-
-                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
-                {
-                    int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
-
-
-                    (*pw)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(wPhaseIdx);
-                    (*pn)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(nPhaseIdx);
-                    (*pc)[dofIdxGlobal] = elemVolVars[scvIdx].capillaryPressure();
-                    (*sw)[dofIdxGlobal] = elemVolVars[scvIdx].saturation(wPhaseIdx);
-                    (*sn)[dofIdxGlobal] = elemVolVars[scvIdx].saturation(nPhaseIdx);
-                    (*rhoW)[dofIdxGlobal] = elemVolVars[scvIdx].density(wPhaseIdx);
-                    (*rhoN)[dofIdxGlobal] = elemVolVars[scvIdx].density(nPhaseIdx);
-                    (*mobW)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(wPhaseIdx);
-                    (*mobN)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(nPhaseIdx);
-                    (*poro)[dofIdxGlobal] = elemVolVars[scvIdx].porosity();
-                    (*Te)[dofIdxGlobal] = elemVolVars[scvIdx].temperature();
-                    if (gravity)
-                        (*ph)[dofIdxGlobal] = elemVolVars[scvIdx].pressureHead(wPhaseIdx);
-                    (*wc)[dofIdxGlobal] = elemVolVars[scvIdx].waterContent(wPhaseIdx);
-
-                }
-
-                // velocity output
-                velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, element, /*phaseIdx=*/0);
+                int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
+
+
+                (*pw)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(wPhaseIdx);
+                (*pn)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(nPhaseIdx);
+                (*pc)[dofIdxGlobal] = elemVolVars[scvIdx].capillaryPressure();
+                (*sw)[dofIdxGlobal] = elemVolVars[scvIdx].saturation(wPhaseIdx);
+                (*sn)[dofIdxGlobal] = elemVolVars[scvIdx].saturation(nPhaseIdx);
+                (*rhoW)[dofIdxGlobal] = elemVolVars[scvIdx].density(wPhaseIdx);
+                (*rhoN)[dofIdxGlobal] = elemVolVars[scvIdx].density(nPhaseIdx);
+                (*mobW)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(wPhaseIdx);
+                (*mobN)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(nPhaseIdx);
+                (*poro)[dofIdxGlobal] = elemVolVars[scvIdx].porosity();
+                (*Te)[dofIdxGlobal] = elemVolVars[scvIdx].temperature();
+                if (gravity)
+                    (*ph)[dofIdxGlobal] = elemVolVars[scvIdx].pressureHead(wPhaseIdx);
+                (*wc)[dofIdxGlobal] = elemVolVars[scvIdx].waterContent(wPhaseIdx);
+
             }
+
+            // velocity output
+            velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, element, /*phaseIdx=*/0);
         }
 
         writer.attachDofData(*sn, "Sn", isBox);
diff --git a/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cniproblem.hh b/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cniproblem.hh
index c449cc59d62973ed43df2dc933e11e85da99611d..ad2860c0473e757b731ad063deb58382e4961895 100644
--- a/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cniproblem.hh
+++ b/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cniproblem.hh
@@ -174,12 +174,8 @@ public:
         mdGrid.startSubDomainMarking();
 
         // subdivide grid in two subdomains
-        for (const auto& element : elements(mdGrid.leafGridView()))
+        for (const auto& element : elements(mdGrid.leafGridView(), Dune::Partitions::interior))
         {
-            // this is required for parallelization, checks if element is within a partition
-            if (element.partitionType() != Dune::InteriorEntity)
-                continue;
-
             GlobalPosition globalPos = element.geometry().center();
 
             if (globalPos[1] > interfacePosY_)
diff --git a/test/multidomain/2cnizeroeq2p2cni/2cnizeroeq2p2cniproblem.hh b/test/multidomain/2cnizeroeq2p2cni/2cnizeroeq2p2cniproblem.hh
index d03d10fa0168bd6f643d4330ad4c60483087bdb8..afd1befda7a55e6eb52ffb8352ff37002b0d954d 100644
--- a/test/multidomain/2cnizeroeq2p2cni/2cnizeroeq2p2cniproblem.hh
+++ b/test/multidomain/2cnizeroeq2p2cni/2cnizeroeq2p2cniproblem.hh
@@ -160,12 +160,8 @@ public:
         mdGrid.startSubDomainMarking();
 
         // subdivide grid in two subdomains
-        for (const auto& element : elements(mdGrid.leafGridView()))
+        for (const auto& element : elements(mdGrid.leafGridView(), Dune::Partitions::interior))
         {
-            // this is required for parallelization checks if element is within a partition
-            if (element.partitionType() != Dune::InteriorEntity)
-                continue;
-
             GlobalPosition globalPos = element.geometry().center();
 
             if (globalPos[1] > interfacePosY_)
diff --git a/test/multidomain/2cstokes2p2c/2cstokes2p2cproblem.hh b/test/multidomain/2cstokes2p2c/2cstokes2p2cproblem.hh
index 88c4fc92ed84c5007f633769e6b06b59cd4b1d78..98949507864f35fbea0d7c5aaa5b337ad71ac829 100644
--- a/test/multidomain/2cstokes2p2c/2cstokes2p2cproblem.hh
+++ b/test/multidomain/2cstokes2p2c/2cstokes2p2cproblem.hh
@@ -173,12 +173,8 @@ public:
         mdGrid.startSubDomainMarking();
 
         // subdivide grid in two subdomains
-        for (const auto& element : elements(mdGrid.leafGridView()))
+        for (const auto& element : elements(mdGrid.leafGridView(), Dune::Partitions::interior))
         {
-            // this is required for parallelization, checks if element is within a partition
-            if (element.partitionType() != Dune::InteriorEntity)
-                continue;
-
             GlobalPosition globalPos = element.geometry().center();
 
             if (globalPos[1] > interfacePosY_)
diff --git a/test/multidomain/2czeroeq2p2c/2czeroeq2p2cproblem.hh b/test/multidomain/2czeroeq2p2c/2czeroeq2p2cproblem.hh
index 7ea688265c45f99dff78df2822a5db76db832e8f..2f1596456bbd4a1d4c9fae00ffc70174a5911984 100644
--- a/test/multidomain/2czeroeq2p2c/2czeroeq2p2cproblem.hh
+++ b/test/multidomain/2czeroeq2p2c/2czeroeq2p2cproblem.hh
@@ -157,13 +157,8 @@ public:
         mdGrid.startSubDomainMarking();
 
         // subdivide grid in two subdomains
-        for (const auto& element : elements(mdGrid.leafGridView()))
+        for (const auto& element : elements(mdGrid.leafGridView(), Dune::Partitions::interior))
         {
-            // this is required for parallelization
-            // checks if element is within a partition
-            if (element.partitionType() != Dune::InteriorEntity)
-                continue;
-
             GlobalPosition globalPos = element.geometry().center();
 
             if (globalPos[1] > interfacePosY_)