diff --git a/dumux/geomechanics/el1p2c/el1p2cmodel.hh b/dumux/geomechanics/el1p2c/el1p2cmodel.hh
index 1db85bdf76e4af4fc6db3dbc0c5e51143a742fb3..8be6716d7fff353cfc7e1a69f0cc0563c04c6e5b 100644
--- a/dumux/geomechanics/el1p2c/el1p2cmodel.hh
+++ b/dumux/geomechanics/el1p2c/el1p2cmodel.hh
@@ -253,130 +253,133 @@ public:
         // loop over all elements (cells)
         for (; eIt != eEndIt; ++eIt)
         {
-            unsigned int eIdx = this->problem_().model().elementMapper().map(*eIt);
-            rank[eIdx] = this->gridView_().comm().rank();
+            if(eIt->partitionType() == Dune::InteriorEntity)
+            {
+                unsigned int eIdx = this->problem_().model().elementMapper().map(*eIt);
+                rank[eIdx] = this->gridView_().comm().rank();
 
-            fvGeometry.update(this->gridView_(), *eIt);
-            elemBcTypes.update(this->problem_(), *eIt, fvGeometry);
-            elemVolVars.update(this->problem_(), *eIt, fvGeometry, false);
+                fvGeometry.update(this->gridView_(), *eIt);
+                elemBcTypes.update(this->problem_(), *eIt, fvGeometry);
+                elemVolVars.update(this->problem_(), *eIt, fvGeometry, false);
 
-            // loop over all local vertices of the cell
+                // loop over all local vertices of the cell
 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
-            int numScv = eIt->subEntities(dim);
+                int numScv = eIt->subEntities(dim);
 #else
-            int numScv = eIt->template count<dim>();
+                int numScv = eIt->template count<dim>();
 #endif
 
-            for (int scvIdx = 0; scvIdx < numScv; ++scvIdx)
-            {
-                unsigned int globalIdx = this->dofMapper().map(*eIt, scvIdx, dim);
-
-                pressure[globalIdx] = elemVolVars[scvIdx].pressure();
-                moleFraction0[globalIdx] = elemVolVars[scvIdx].moleFraction(0);
-                moleFraction1[globalIdx] = elemVolVars[scvIdx].moleFraction(1);
-                massFraction0[globalIdx] = elemVolVars[scvIdx].massFraction(0);
-                massFraction1[globalIdx] = 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[globalIdx] = tmpDispl;
-                    }
+                for (int scvIdx = 0; scvIdx < numScv; ++scvIdx)
+                {
+                    unsigned int globalIdx = this->dofMapper().map(*eIt, scvIdx, dim);
+
+                    pressure[globalIdx] = elemVolVars[scvIdx].pressure();
+                    moleFraction0[globalIdx] = elemVolVars[scvIdx].moleFraction(0);
+                    moleFraction1[globalIdx] = elemVolVars[scvIdx].moleFraction(1);
+                    massFraction0[globalIdx] = elemVolVars[scvIdx].massFraction(0);
+                    massFraction1[globalIdx] = 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[globalIdx] = tmpDispl;
+                        }
 
-                else
-                    displacement[globalIdx] = elemVolVars[scvIdx].displacement();
-
-                density[globalIdx] = elemVolVars[scvIdx].density();
-                viscosity[globalIdx] = elemVolVars[scvIdx].viscosity();
-                porosity[globalIdx] = elemVolVars[scvIdx].porosity();
-                Kx[globalIdx] =    this->problem_().spatialParams().intrinsicPermeability(
-                                *eIt, 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(
-                                *eIt, 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 faceIdx = 0; faceIdx < fvGeometry.numScvf; faceIdx++) {
-
-                //prepare the flux calculations (set up and prepare geometry, FE gradients)
-                FluxVariables fluxVars(this->problem_(),
-                                    *eIt, fvGeometry,
-                                    faceIdx,
-                                    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;
+                    else
+                        displacement[globalIdx] = elemVolVars[scvIdx].displacement();
+
+                    density[globalIdx] = elemVolVars[scvIdx].density();
+                    viscosity[globalIdx] = elemVolVars[scvIdx].viscosity();
+                    porosity[globalIdx] = elemVolVars[scvIdx].porosity();
+                    Kx[globalIdx] =    this->problem_().spatialParams().intrinsicPermeability(
+                                    *eIt, 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(
+                                    *eIt, 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 faceIdx = 0; faceIdx < fvGeometry.numScvf; faceIdx++) {
+
+                    //prepare the flux calculations (set up and prepare geometry, FE gradients)
+                    FluxVariables fluxVars(this->problem_(),
+                                        *eIt, fvGeometry,
+                                        faceIdx,
+                                        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];
+                        }
+                    }
+                }
 
+                // calculate total stresses
                 // in case of rock mechanics sign convention compressive stresses
-                // are defined to be positive
+                // are defined to be positive and total stress is calculated by adding the pore pressure
                 if(rockMechanicsSignConvention_){
-                    effStressX[eIdx] -= tmpEffStress[0];
+                    totalStressX[eIdx][0] = effStressX[eIdx][0]    + cellPressure[eIdx];
+                    totalStressX[eIdx][1] = effStressX[eIdx][1];
+                    totalStressX[eIdx][2] = effStressX[eIdx][2];
                     if (dim >= 2) {
-                        effStressY[eIdx] -= tmpEffStress[1];
+                        totalStressY[eIdx][0] = effStressY[eIdx][0];
+                        totalStressY[eIdx][1] = effStressY[eIdx][1]    + cellPressure[eIdx];
+                        totalStressY[eIdx][2] = effStressY[eIdx][2];
                     }
                     if (dim >= 3) {
-                        effStressZ[eIdx] -= tmpEffStress[2];
+                        totalStressZ[eIdx][0] = effStressZ[eIdx][0];
+                        totalStressZ[eIdx][1] = effStressZ[eIdx][1];
+                        totalStressZ[eIdx][2] = effStressZ[eIdx][2]    + cellPressure[eIdx];
                     }
                 }
                 else{
-                    effStressX[eIdx] += tmpEffStress[0];
+                    totalStressX[eIdx][0] = effStressX[eIdx][0]    - cellPressure[eIdx];
+                    totalStressX[eIdx][1] = effStressX[eIdx][1];
+                    totalStressX[eIdx][2] = effStressX[eIdx][2];
                     if (dim >= 2) {
-                        effStressY[eIdx] += tmpEffStress[1];
+                        totalStressY[eIdx][0] = effStressY[eIdx][0];
+                        totalStressY[eIdx][1] = effStressY[eIdx][1]    - cellPressure[eIdx];
+                        totalStressY[eIdx][2] = effStressY[eIdx][2];
                     }
                     if (dim >= 3) {
-                        effStressZ[eIdx] += tmpEffStress[2];
+                        totalStressZ[eIdx][0] = effStressZ[eIdx][0];
+                        totalStressZ[eIdx][1] = effStressZ[eIdx][1];
+                        totalStressZ[eIdx][2] = effStressZ[eIdx][2]    - cellPressure[eIdx];
                     }
                 }
             }
-
-            // 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/el2pmodel.hh b/dumux/geomechanics/el2p/el2pmodel.hh
index f0d2bf48ecf337c2761e49c44f1e94572ebf6bf5..0520132e150f25500eb8ef4d92cd5828d3591703 100644
--- a/dumux/geomechanics/el2p/el2pmodel.hh
+++ b/dumux/geomechanics/el2p/el2pmodel.hh
@@ -284,6 +284,8 @@ public:
 
 
         for (unsigned int eIdx = 0; eIdx < numElements; ++eIdx) {
+            if(eIt->partitionType() == Dune::InteriorEntity)
+            {
             deltaEffStressX[eIdx] = Scalar(0.0);
             if (dim >= 2)
                 deltaEffStressY[eIdx] = Scalar(0.0);
@@ -630,7 +632,7 @@ public:
             // Pressure margins according to J. Rutqvist et al. / International Journal of Rock Mecahnics & Mining Sciences 45 (2008), 132-143
             Pcrtens[eIdx] = Peff - principalStress3[eIdx];
             Pcrshe[eIdx] = Peff - Psc;
-
+        }
         }
 
         writer.attachVertexData(Te, "T");
diff --git a/dumux/geomechanics/elastic/elasticmodel.hh b/dumux/geomechanics/elastic/elasticmodel.hh
index 83652bf9bc40e8400a6a2f6ee0de6e81fb4230b9..d67539221bc686cd84cebfa72a87ef74cbca8d9d 100644
--- a/dumux/geomechanics/elastic/elasticmodel.hh
+++ b/dumux/geomechanics/elastic/elasticmodel.hh
@@ -124,6 +124,8 @@ public:
         ElementIterator eEndIt = this->gridView_().template end<0>();
         for (; eIt != eEndIt; ++eIt)
         {
+            if(eIt->partitionType() == Dune::InteriorEntity)
+            {
             int eIdx = this->problem_().model().elementMapper().map(*eIt);
             rank[eIdx] = this->gridView_().comm().rank();
 
@@ -184,6 +186,7 @@ public:
                     sigmaz[eIdx] += stress[2];
                 }
             }
+            }
         }
 
 
diff --git a/dumux/implicit/1p/1pmodel.hh b/dumux/implicit/1p/1pmodel.hh
index 0148ddc7ce388b702f18c699a4b2045b73eb1b64..009b292da1717ebde4909a66df43ce743d7ec684 100644
--- a/dumux/implicit/1p/1pmodel.hh
+++ b/dumux/implicit/1p/1pmodel.hh
@@ -105,32 +105,35 @@ public:
         ElementIterator eEndIt = this->gridView_().template end<0>();
         for (; eIt != eEndIt; ++eIt)
         {
-            int eIdx = this->problem_().model().elementMapper().map(*eIt);
-            (*rank)[eIdx] = this->gridView_().comm().rank();
+            if(eIt->partitionType() == Dune::InteriorEntity)
+            {
+                int eIdx = this->problem_().model().elementMapper().map(*eIt);
+                (*rank)[eIdx] = this->gridView_().comm().rank();
 
-            FVElementGeometry fvGeometry;
-            fvGeometry.update(this->gridView_(), *eIt);
+                FVElementGeometry fvGeometry;
+                fvGeometry.update(this->gridView_(), *eIt);
 
-            ElementVolumeVariables elemVolVars;
-            elemVolVars.update(this->problem_(),
-                               *eIt,
-                               fvGeometry,
-                               false /* oldSol? */);
+                ElementVolumeVariables elemVolVars;
+                elemVolVars.update(this->problem_(),
+                                   *eIt,
+                                   fvGeometry,
+                                   false /* oldSol? */);
 
-            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
-            {
-                int globalIdx = this->dofMapper().map(*eIt, scvIdx, dofCodim);
+                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+                {
+                    int globalIdx = this->dofMapper().map(*eIt, scvIdx, dofCodim);
 
-                const SpatialParams &spatialParams = this->problem_().spatialParams();
+                    const SpatialParams &spatialParams = this->problem_().spatialParams();
 
-                (*p)[globalIdx] = elemVolVars[scvIdx].pressure();
-                (*K)[globalIdx] = spatialParams.intrinsicPermeability(*eIt,
-                                                                     fvGeometry,
-                                                                     scvIdx);
-            }
+                    (*p)[globalIdx] = elemVolVars[scvIdx].pressure();
+                    (*K)[globalIdx] = spatialParams.intrinsicPermeability(*eIt,
+                                                                         fvGeometry,
+                                                                         scvIdx);
+                }
 
-            // velocity output
-            velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, *eIt, /*phaseIdx=*/0);
+                // velocity output
+                velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, *eIt, /*phaseIdx=*/0);
+            }
         }
 
         writer.attachDofData(*p, "p", isBox);
diff --git a/dumux/implicit/1p2c/1p2cmodel.hh b/dumux/implicit/1p2c/1p2cmodel.hh
index 85ac63fd3d0820f6266efd0263b47492a4ed3817..4abe9faedfd1c50fe6576c44014c34f3db46e2bc 100644
--- a/dumux/implicit/1p2c/1p2cmodel.hh
+++ b/dumux/implicit/1p2c/1p2cmodel.hh
@@ -135,34 +135,37 @@ public:
         ElementIterator eEndIt = this->gridView_().template end<0>();
         for (; eIt != eEndIt; ++eIt)
         {
-            int eIdx = this->problem_().model().elementMapper().map(*eIt);
-            rank[eIdx] = this->gridView_().comm().rank();
-
-            FVElementGeometry fvGeometry;
-            fvGeometry.update(this->gridView_(), *eIt);
-
-            ElementVolumeVariables elemVolVars;
-            elemVolVars.update(this->problem_(),
-                               *eIt,
-                               fvGeometry,
-                               false /* oldSol? */);
-
-            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+            if(eIt->partitionType() == Dune::InteriorEntity)
             {
-                int globalIdx = this->dofMapper().map(*eIt, scvIdx, dofCodim);
-
-                pressure[globalIdx] = elemVolVars[scvIdx].pressure();
-                delp[globalIdx] = elemVolVars[scvIdx].pressure() - 1e5;
-                moleFraction0[globalIdx] = elemVolVars[scvIdx].moleFraction(0);
-                moleFraction1[globalIdx] = elemVolVars[scvIdx].moleFraction(1);
-                massFraction0[globalIdx] = elemVolVars[scvIdx].massFraction(0);
-                massFraction1[globalIdx] = elemVolVars[scvIdx].massFraction(1);
-                rho[globalIdx] = elemVolVars[scvIdx].density();
-                mu[globalIdx] = elemVolVars[scvIdx].viscosity();
+                int eIdx = this->problem_().model().elementMapper().map(*eIt);
+                rank[eIdx] = this->gridView_().comm().rank();
+
+                FVElementGeometry fvGeometry;
+                fvGeometry.update(this->gridView_(), *eIt);
+
+                ElementVolumeVariables elemVolVars;
+                elemVolVars.update(this->problem_(),
+                                   *eIt,
+                                   fvGeometry,
+                                   false /* oldSol? */);
+
+                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+                {
+                    int globalIdx = this->dofMapper().map(*eIt, scvIdx, dofCodim);
+
+                    pressure[globalIdx] = elemVolVars[scvIdx].pressure();
+                    delp[globalIdx] = elemVolVars[scvIdx].pressure() - 1e5;
+                    moleFraction0[globalIdx] = elemVolVars[scvIdx].moleFraction(0);
+                    moleFraction1[globalIdx] = elemVolVars[scvIdx].moleFraction(1);
+                    massFraction0[globalIdx] = elemVolVars[scvIdx].massFraction(0);
+                    massFraction1[globalIdx] = elemVolVars[scvIdx].massFraction(1);
+                    rho[globalIdx] = elemVolVars[scvIdx].density();
+                    mu[globalIdx] = elemVolVars[scvIdx].viscosity();
+                }
+
+                // velocity output
+                velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, *eIt, phaseIdx);
             }
-
-            // velocity output
-            velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, *eIt, phaseIdx);
         }
 
         writer.attachDofData(pressure, "P", isBox);
diff --git a/dumux/implicit/2p/2pmodel.hh b/dumux/implicit/2p/2pmodel.hh
index ed79bd4af056a622a953d466db3bbe8fbaafd71d..ca3d004fa0335d7e5febc5f7d2e02c668d8ce077 100644
--- a/dumux/implicit/2p/2pmodel.hh
+++ b/dumux/implicit/2p/2pmodel.hh
@@ -144,38 +144,41 @@ public:
         ElementIterator eEndIt = this->gridView_().template end<0>();
         for (; eIt != eEndIt; ++eIt)
         {
-            int eIdx = this->elementMapper().map(*eIt);
-            (*rank)[eIdx] = this->gridView_().comm().rank();
-
-            FVElementGeometry fvGeometry;
-            fvGeometry.update(this->gridView_(), *eIt);
-
-            ElementVolumeVariables elemVolVars;
-            elemVolVars.update(this->problem_(),
-                               *eIt,
-                               fvGeometry,
-                               false /* oldSol? */);
-
-            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+            if(eIt->partitionType() == Dune::InteriorEntity)
             {
-                int globalIdx = this->dofMapper().map(*eIt, scvIdx, dofCodim);
-
-                (*pw)[globalIdx] = elemVolVars[scvIdx].pressure(wPhaseIdx);
-                (*pn)[globalIdx] = elemVolVars[scvIdx].pressure(nPhaseIdx);
-                (*pc)[globalIdx] = elemVolVars[scvIdx].capillaryPressure();
-                (*sw)[globalIdx] = elemVolVars[scvIdx].saturation(wPhaseIdx);
-                (*sn)[globalIdx] = elemVolVars[scvIdx].saturation(nPhaseIdx);
-                (*rhoW)[globalIdx] = elemVolVars[scvIdx].density(wPhaseIdx);
-                (*rhoN)[globalIdx] = elemVolVars[scvIdx].density(nPhaseIdx);
-                (*mobW)[globalIdx] = elemVolVars[scvIdx].mobility(wPhaseIdx);
-                (*mobN)[globalIdx] = elemVolVars[scvIdx].mobility(nPhaseIdx);
-                (*poro)[globalIdx] = elemVolVars[scvIdx].porosity();
-                (*Te)[globalIdx] = elemVolVars[scvIdx].temperature();
+                int eIdx = this->elementMapper().map(*eIt);
+                (*rank)[eIdx] = this->gridView_().comm().rank();
+
+                FVElementGeometry fvGeometry;
+                fvGeometry.update(this->gridView_(), *eIt);
+
+                ElementVolumeVariables elemVolVars;
+                elemVolVars.update(this->problem_(),
+                                   *eIt,
+                                   fvGeometry,
+                                   false /* oldSol? */);
+
+                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+                {
+                    int globalIdx = this->dofMapper().map(*eIt, scvIdx, dofCodim);
+
+                    (*pw)[globalIdx] = elemVolVars[scvIdx].pressure(wPhaseIdx);
+                    (*pn)[globalIdx] = elemVolVars[scvIdx].pressure(nPhaseIdx);
+                    (*pc)[globalIdx] = elemVolVars[scvIdx].capillaryPressure();
+                    (*sw)[globalIdx] = elemVolVars[scvIdx].saturation(wPhaseIdx);
+                    (*sn)[globalIdx] = elemVolVars[scvIdx].saturation(nPhaseIdx);
+                    (*rhoW)[globalIdx] = elemVolVars[scvIdx].density(wPhaseIdx);
+                    (*rhoN)[globalIdx] = elemVolVars[scvIdx].density(nPhaseIdx);
+                    (*mobW)[globalIdx] = elemVolVars[scvIdx].mobility(wPhaseIdx);
+                    (*mobN)[globalIdx] = elemVolVars[scvIdx].mobility(nPhaseIdx);
+                    (*poro)[globalIdx] = elemVolVars[scvIdx].porosity();
+                    (*Te)[globalIdx] = elemVolVars[scvIdx].temperature();
+                }
+
+                // velocity output
+                velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, *eIt, wPhaseIdx);
+                velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, nPhaseIdx);
             }
-
-            // velocity output
-            velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, *eIt, wPhaseIdx);
-            velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, nPhaseIdx);
         }
 
         writer.attachDofData(*sn, "Sn", isBox);
diff --git a/dumux/implicit/2p2c/2p2cmodel.hh b/dumux/implicit/2p2c/2p2cmodel.hh
index fbe9e9f104f14577832ae8d15c44a6838198c72e..5c488079324aba06c208105a796607f9bca82d1a 100644
--- a/dumux/implicit/2p2c/2p2cmodel.hh
+++ b/dumux/implicit/2p2c/2p2cmodel.hh
@@ -333,57 +333,60 @@ public:
         ElementIterator eEndIt = this->gridView_().template end<0>();
         for (; eIt != eEndIt; ++eIt)
         {
-            int eIdx = this->elementMapper().map(*eIt);
-            (*rank)[eIdx] = this->gridView_().comm().rank();
+            if(eIt->partitionType() == Dune::InteriorEntity)
+            {
+                int eIdx = this->elementMapper().map(*eIt);
+                (*rank)[eIdx] = this->gridView_().comm().rank();
 
-            FVElementGeometry fvGeometry;
-            fvGeometry.update(this->gridView_(), *eIt);
+                FVElementGeometry fvGeometry;
+                fvGeometry.update(this->gridView_(), *eIt);
 
-            ElementVolumeVariables elemVolVars;
-            elemVolVars.update(this->problem_(),
-                               *eIt,
-                               fvGeometry,
-                               false /* oldSol? */);
+                ElementVolumeVariables elemVolVars;
+                elemVolVars.update(this->problem_(),
+                                   *eIt,
+                                   fvGeometry,
+                                   false /* oldSol? */);
 
-            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
-            {
-                int globalIdx = this->dofMapper().map(*eIt, scvIdx, dofCodim);
-
-                (*sN)[globalIdx]    = elemVolVars[scvIdx].saturation(nPhaseIdx);
-                (*sW)[globalIdx]    = elemVolVars[scvIdx].saturation(wPhaseIdx);
-                (*pn)[globalIdx]    = elemVolVars[scvIdx].pressure(nPhaseIdx);
-                (*pw)[globalIdx]    = elemVolVars[scvIdx].pressure(wPhaseIdx);
-                (*pc)[globalIdx]    = elemVolVars[scvIdx].capillaryPressure();
-                (*rhoW)[globalIdx]  = elemVolVars[scvIdx].density(wPhaseIdx);
-                (*rhoN)[globalIdx]  = elemVolVars[scvIdx].density(nPhaseIdx);
-                (*mobW)[globalIdx]  = elemVolVars[scvIdx].mobility(wPhaseIdx);
-                (*mobN)[globalIdx]  = elemVolVars[scvIdx].mobility(nPhaseIdx);
-                for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-                    for (int compIdx = 0; compIdx < numComponents; ++compIdx)
-                    {
-                        (*massFrac[phaseIdx][compIdx])[globalIdx]
-                            = elemVolVars[scvIdx].massFraction(phaseIdx, compIdx);
+                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+                {
+                    int globalIdx = this->dofMapper().map(*eIt, scvIdx, dofCodim);
+
+                    (*sN)[globalIdx]    = elemVolVars[scvIdx].saturation(nPhaseIdx);
+                    (*sW)[globalIdx]    = elemVolVars[scvIdx].saturation(wPhaseIdx);
+                    (*pn)[globalIdx]    = elemVolVars[scvIdx].pressure(nPhaseIdx);
+                    (*pw)[globalIdx]    = elemVolVars[scvIdx].pressure(wPhaseIdx);
+                    (*pc)[globalIdx]    = elemVolVars[scvIdx].capillaryPressure();
+                    (*rhoW)[globalIdx]  = elemVolVars[scvIdx].density(wPhaseIdx);
+                    (*rhoN)[globalIdx]  = elemVolVars[scvIdx].density(nPhaseIdx);
+                    (*mobW)[globalIdx]  = elemVolVars[scvIdx].mobility(wPhaseIdx);
+                    (*mobN)[globalIdx]  = elemVolVars[scvIdx].mobility(nPhaseIdx);
+                    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                        for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+                        {
+                            (*massFrac[phaseIdx][compIdx])[globalIdx]
+                                = elemVolVars[scvIdx].massFraction(phaseIdx, compIdx);
+
+                            Valgrind::CheckDefined((*massFrac[phaseIdx][compIdx])[globalIdx][0]);
+                        }
+                    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                        for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+                        {
+                            (*moleFrac[phaseIdx][compIdx])[globalIdx]
+                                = elemVolVars[scvIdx].moleFraction(phaseIdx, compIdx);
+
+                            Valgrind::CheckDefined((*moleFrac[phaseIdx][compIdx])[globalIdx][0]);
+                        }
+                    (*poro)[globalIdx]  = elemVolVars[scvIdx].porosity();
+                    (*temperature)[globalIdx] = elemVolVars[scvIdx].temperature();
+                    (*phasePresence)[globalIdx]
+                        = staticDat_[globalIdx].phasePresence;
+                }
 
-                        Valgrind::CheckDefined((*massFrac[phaseIdx][compIdx])[globalIdx][0]);
-                    }
-                for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-                    for (int compIdx = 0; compIdx < numComponents; ++compIdx)
-                    {
-                        (*moleFrac[phaseIdx][compIdx])[globalIdx]
-                            = elemVolVars[scvIdx].moleFraction(phaseIdx, compIdx);
-
-                        Valgrind::CheckDefined((*moleFrac[phaseIdx][compIdx])[globalIdx][0]);
-                    }
-                (*poro)[globalIdx]  = elemVolVars[scvIdx].porosity();
-                (*temperature)[globalIdx] = elemVolVars[scvIdx].temperature();
-                (*phasePresence)[globalIdx]
-                    = staticDat_[globalIdx].phasePresence;
+                // velocity output
+                velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, *eIt, wPhaseIdx);
+                velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, nPhaseIdx);
             }
 
-            // velocity output
-            velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, *eIt, wPhaseIdx);
-            velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, nPhaseIdx);
-
         } // loop over elements
 
         writer.attachDofData(*sN,     "Sn", isBox);
diff --git a/dumux/implicit/3p/3pmodel.hh b/dumux/implicit/3p/3pmodel.hh
index 25bdacce9d6a755743847233b8ddee96e06f8196..18fa8c24910472c4de45ff7d6e2ede4f57ff350c 100644
--- a/dumux/implicit/3p/3pmodel.hh
+++ b/dumux/implicit/3p/3pmodel.hh
@@ -140,39 +140,41 @@ public:
         ElementIterator eEndIt = this->gridView_().template end<0>();
         for (; eIt != eEndIt; ++eIt)
         {
-            int eIdx = this->problem_().elementMapper().map(*eIt);
-            (*rank)[eIdx] = this->gridView_().comm().rank();
+            if(eIt->partitionType() == Dune::InteriorEntity)
+            {
+                int eIdx = this->problem_().elementMapper().map(*eIt);
+                (*rank)[eIdx] = this->gridView_().comm().rank();
 
-            FVElementGeometry fvGeometry;
-            fvGeometry.update(this->gridView_(), *eIt);
+                FVElementGeometry fvGeometry;
+                fvGeometry.update(this->gridView_(), *eIt);
 
 
-            ElementVolumeVariables elemVolVars;
-            elemVolVars.update(this->problem_(),
-                               *eIt,
-                               fvGeometry,
-                               false /* oldSol? */);
+                ElementVolumeVariables elemVolVars;
+                elemVolVars.update(this->problem_(),
+                                   *eIt,
+                                   fvGeometry,
+                                   false /* oldSol? */);
 
-            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
-            {
-                int globalIdx = this->dofMapper().map(*eIt, scvIdx, dofCodim);
+                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+                {
+                    int globalIdx = this->dofMapper().map(*eIt, scvIdx, dofCodim);
+
+                    for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+                        (*saturation[phaseIdx])[globalIdx] = elemVolVars[scvIdx].saturation(phaseIdx);
+                        (*pressure[phaseIdx])[globalIdx] = elemVolVars[scvIdx].pressure(phaseIdx);
+                        (*density[phaseIdx])[globalIdx] = elemVolVars[scvIdx].density(phaseIdx);
+                    }
 
-                for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
-                    (*saturation[phaseIdx])[globalIdx] = elemVolVars[scvIdx].saturation(phaseIdx);
-                    (*pressure[phaseIdx])[globalIdx] = elemVolVars[scvIdx].pressure(phaseIdx);
-                    (*density[phaseIdx])[globalIdx] = elemVolVars[scvIdx].density(phaseIdx);
+                    (*poro)[globalIdx] = elemVolVars[scvIdx].porosity();
+                    (*perm)[globalIdx] = elemVolVars[scvIdx].permeability();
+                    (*temperature)[globalIdx] = elemVolVars[scvIdx].temperature();
                 }
 
-                (*poro)[globalIdx] = elemVolVars[scvIdx].porosity();
-                (*perm)[globalIdx] = elemVolVars[scvIdx].permeability();
-                (*temperature)[globalIdx] = elemVolVars[scvIdx].temperature();
+                // velocity output
+                velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, *eIt, wPhaseIdx);
+                velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, nPhaseIdx);
+                velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, gPhaseIdx);
             }
-
-            // velocity output
-            velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, *eIt, wPhaseIdx);
-            velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, nPhaseIdx);
-            velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, gPhaseIdx);
-
         }
 
         writer.attachDofData(*saturation[wPhaseIdx], "sw", isBox);
diff --git a/dumux/implicit/3p3c/3p3cmodel.hh b/dumux/implicit/3p3c/3p3cmodel.hh
index 0ef316f2bbe5f95a16b80e900d092ce856fc1fe9..0c8c10fc269615870f8008a7d03f2974e64df1b0 100644
--- a/dumux/implicit/3p3c/3p3cmodel.hh
+++ b/dumux/implicit/3p3c/3p3cmodel.hh
@@ -335,50 +335,52 @@ public:
         ElementIterator eEndIt = this->gridView_().template end<0>();
         for (; eIt != eEndIt; ++eIt)
         {
-            int eIdx = this->problem_().elementMapper().map(*eIt);
-            (*rank)[eIdx] = this->gridView_().comm().rank();
+            if(eIt->partitionType() == Dune::InteriorEntity)
+            {
+                int eIdx = this->problem_().elementMapper().map(*eIt);
+                (*rank)[eIdx] = this->gridView_().comm().rank();
 
-            FVElementGeometry fvGeometry;
-            fvGeometry.update(this->gridView_(), *eIt);
+                FVElementGeometry fvGeometry;
+                fvGeometry.update(this->gridView_(), *eIt);
 
 
-            ElementVolumeVariables elemVolVars;
-            elemVolVars.update(this->problem_(),
-                               *eIt,
-                               fvGeometry,
-                               false /* oldSol? */);
+                ElementVolumeVariables elemVolVars;
+                elemVolVars.update(this->problem_(),
+                                   *eIt,
+                                   fvGeometry,
+                                   false /* oldSol? */);
 
-            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
-            {
-                int globalIdx = this->dofMapper().map(*eIt, scvIdx, dofCodim);
+                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+                {
+                    int globalIdx = this->dofMapper().map(*eIt, scvIdx, dofCodim);
 
-                for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
-                    (*saturation[phaseIdx])[globalIdx] = elemVolVars[scvIdx].saturation(phaseIdx);
-                    (*pressure[phaseIdx])[globalIdx] = elemVolVars[scvIdx].pressure(phaseIdx);
-                    (*density[phaseIdx])[globalIdx] = elemVolVars[scvIdx].density(phaseIdx);
-                }
+                    for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+                        (*saturation[phaseIdx])[globalIdx] = elemVolVars[scvIdx].saturation(phaseIdx);
+                        (*pressure[phaseIdx])[globalIdx] = elemVolVars[scvIdx].pressure(phaseIdx);
+                        (*density[phaseIdx])[globalIdx] = elemVolVars[scvIdx].density(phaseIdx);
+                    }
 
-                for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-                    for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
-                        (*moleFraction[phaseIdx][compIdx])[globalIdx] =
-                            elemVolVars[scvIdx].moleFraction(phaseIdx,
-                                                              compIdx);
+                    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+                        for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
+                            (*moleFraction[phaseIdx][compIdx])[globalIdx] =
+                                elemVolVars[scvIdx].moleFraction(phaseIdx,
+                                                                  compIdx);
 
-                        Valgrind::CheckDefined((*moleFraction[phaseIdx][compIdx])[globalIdx]);
+                            Valgrind::CheckDefined((*moleFraction[phaseIdx][compIdx])[globalIdx]);
+                        }
                     }
+
+                    (*poro)[globalIdx] = elemVolVars[scvIdx].porosity();
+                    (*perm)[globalIdx] = elemVolVars[scvIdx].permeability();
+                    (*temperature)[globalIdx] = elemVolVars[scvIdx].temperature();
+                    (*phasePresence)[globalIdx] = staticDat_[globalIdx].phasePresence;
                 }
 
-                (*poro)[globalIdx] = elemVolVars[scvIdx].porosity();
-                (*perm)[globalIdx] = elemVolVars[scvIdx].permeability();
-                (*temperature)[globalIdx] = elemVolVars[scvIdx].temperature();
-                (*phasePresence)[globalIdx] = staticDat_[globalIdx].phasePresence;
+                // velocity output
+                velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, *eIt, wPhaseIdx);
+                velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, nPhaseIdx);
+                velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, gPhaseIdx);
             }
-
-            // velocity output
-            velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, *eIt, wPhaseIdx);
-            velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, nPhaseIdx);
-            velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, gPhaseIdx);
-
         }
 
         writer.attachDofData(*saturation[wPhaseIdx], "Sw", isBox);
diff --git a/dumux/implicit/mpnc/mpncvtkwriter.hh b/dumux/implicit/mpnc/mpncvtkwriter.hh
index 2e7163678f3d4b71e19c2ec02834fb0f0cd37512..5c58049da3afa81817a8d25d4bcb7d4fe1c59719 100644
--- a/dumux/implicit/mpnc/mpncvtkwriter.hh
+++ b/dumux/implicit/mpnc/mpncvtkwriter.hh
@@ -80,33 +80,36 @@ public:
         ElementIterator eEndIt = problem_.gridView().template end<0>();
         for (; eIt != eEndIt; ++eIt)
         {
-            fvGeometry.update(problem_.gridView(), *eIt);
-            elemBcTypes.update(problem_, *eIt);
-            this->problem_.model().setHints(*eIt, elemVolVars);
-            elemVolVars.update(problem_,
-                               *eIt,
-                               fvGeometry,
-                               false);
-            this->problem_.model().updateCurHints(*eIt, elemVolVars);
+            if(eIt->partitionType() == Dune::InteriorEntity)
+            {
+                fvGeometry.update(problem_.gridView(), *eIt);
+                elemBcTypes.update(problem_, *eIt);
+                this->problem_.model().setHints(*eIt, elemVolVars);
+                elemVolVars.update(problem_,
+                                   *eIt,
+                                   fvGeometry,
+                                   false);
+                this->problem_.model().updateCurHints(*eIt, elemVolVars);
 
-            // tell the sub-writers to do what ever they need to with
-            // their internal buffers when a given element is seen.
-            commonWriter_.processElement(*eIt,
-                                         fvGeometry,
-                                         elemVolVars,
-                                         elemBcTypes);
-            massWriter_.processElement(*eIt,
-                                       fvGeometry,
-                                       elemVolVars,
-                                       elemBcTypes);
-            energyWriter_.processElement(*eIt,
-                                         fvGeometry,
-                                         elemVolVars,
-                                         elemBcTypes);
-            customWriter_.processElement(*eIt,
-                                         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(*eIt,
+                                             fvGeometry,
+                                             elemVolVars,
+                                             elemBcTypes);
+                massWriter_.processElement(*eIt,
+                                           fvGeometry,
+                                           elemVolVars,
+                                           elemBcTypes);
+                energyWriter_.processElement(*eIt,
+                                             fvGeometry,
+                                             elemVolVars,
+                                             elemBcTypes);
+                customWriter_.processElement(*eIt,
+                                             fvGeometry,
+                                             elemVolVars,
+                                             elemBcTypes);
+            }
         }
 
         // write everything to the output file
diff --git a/dumux/implicit/nonisothermal/nimodel.hh b/dumux/implicit/nonisothermal/nimodel.hh
index b2ff963493eb22ec458f17804489c6ec43086d4d..5bba98d120f176d8a23f1410ac61ea440e2275ff 100644
--- a/dumux/implicit/nonisothermal/nimodel.hh
+++ b/dumux/implicit/nonisothermal/nimodel.hh
@@ -74,19 +74,22 @@ public:
         ElementIterator eEndIt = this->gridView().template end<0>();
         for (; eIt != eEndIt; ++eIt)
         {
-            FVElementGeometry fvGeometry;
-            fvGeometry.update(this->gridView_(), *eIt);
+            if(eIt->partitionType() == Dune::InteriorEntity)
+            {
+                FVElementGeometry fvGeometry;
+                fvGeometry.update(this->gridView_(), *eIt);
 
-            ElementVolumeVariables elemVolVars;
-            elemVolVars.update(this->problem_(),
-                               *eIt,
-                               fvGeometry,
-                               false /* oldSol? */);
+                ElementVolumeVariables elemVolVars;
+                elemVolVars.update(this->problem_(),
+                                   *eIt,
+                                   fvGeometry,
+                                   false /* oldSol? */);
 
-            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
-            {
-                int globalIdx = this->dofMapper().map(*eIt, scvIdx, dofCodim);
-                temperature[globalIdx] = elemVolVars[scvIdx].temperature();
+                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+                {
+                    int globalIdx = this->dofMapper().map(*eIt, scvIdx, dofCodim);
+                    temperature[globalIdx] = elemVolVars[scvIdx].temperature();
+                }
             }
         }
 
diff --git a/dumux/implicit/richards/richardsmodel.hh b/dumux/implicit/richards/richardsmodel.hh
index d9c4ffe2a440f989316a0cda35a5d8f3dc688eea..a44836c5741c1b25fbc932e0a8c93914c4f1775e 100644
--- a/dumux/implicit/richards/richardsmodel.hh
+++ b/dumux/implicit/richards/richardsmodel.hh
@@ -172,37 +172,40 @@ public:
         ElementIterator eEndIt = this->gridView_().template end<0>();
         for (; eIt != eEndIt; ++eIt)
         {
-            int eIdx = this->problem_().model().elementMapper().map(*eIt);
-            (*rank)[eIdx] = this->gridView_().comm().rank();
-
-            FVElementGeometry fvGeometry;
-            fvGeometry.update(this->gridView_(), *eIt);
-
-            ElementVolumeVariables elemVolVars;
-            elemVolVars.update(this->problem_(),
-                               *eIt,
-                               fvGeometry,
-                               false /* oldSol? */);
-
-            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+            if(eIt->partitionType() == Dune::InteriorEntity)
             {
-                int globalIdx = this->dofMapper().map(*eIt, scvIdx, dofCodim);
-
-                (*pw)[globalIdx] = elemVolVars[scvIdx].pressure(wPhaseIdx);
-                (*pn)[globalIdx] = elemVolVars[scvIdx].pressure(nPhaseIdx);
-                (*pc)[globalIdx] = elemVolVars[scvIdx].capillaryPressure();
-                (*sw)[globalIdx] = elemVolVars[scvIdx].saturation(wPhaseIdx);
-                (*sn)[globalIdx] = elemVolVars[scvIdx].saturation(nPhaseIdx);
-                (*rhoW)[globalIdx] = elemVolVars[scvIdx].density(wPhaseIdx);
-                (*rhoN)[globalIdx] = elemVolVars[scvIdx].density(nPhaseIdx);
-                (*mobW)[globalIdx] = elemVolVars[scvIdx].mobility(wPhaseIdx);
-                (*mobN)[globalIdx] = elemVolVars[scvIdx].mobility(nPhaseIdx);
-                (*poro)[globalIdx] = elemVolVars[scvIdx].porosity();
-                (*Te)[globalIdx] = elemVolVars[scvIdx].temperature();
+                int eIdx = this->problem_().model().elementMapper().map(*eIt);
+                (*rank)[eIdx] = this->gridView_().comm().rank();
+
+                FVElementGeometry fvGeometry;
+                fvGeometry.update(this->gridView_(), *eIt);
+
+                ElementVolumeVariables elemVolVars;
+                elemVolVars.update(this->problem_(),
+                                   *eIt,
+                                   fvGeometry,
+                                   false /* oldSol? */);
+
+                for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+                {
+                    int globalIdx = this->dofMapper().map(*eIt, scvIdx, dofCodim);
+
+                    (*pw)[globalIdx] = elemVolVars[scvIdx].pressure(wPhaseIdx);
+                    (*pn)[globalIdx] = elemVolVars[scvIdx].pressure(nPhaseIdx);
+                    (*pc)[globalIdx] = elemVolVars[scvIdx].capillaryPressure();
+                    (*sw)[globalIdx] = elemVolVars[scvIdx].saturation(wPhaseIdx);
+                    (*sn)[globalIdx] = elemVolVars[scvIdx].saturation(nPhaseIdx);
+                    (*rhoW)[globalIdx] = elemVolVars[scvIdx].density(wPhaseIdx);
+                    (*rhoN)[globalIdx] = elemVolVars[scvIdx].density(nPhaseIdx);
+                    (*mobW)[globalIdx] = elemVolVars[scvIdx].mobility(wPhaseIdx);
+                    (*mobN)[globalIdx] = elemVolVars[scvIdx].mobility(nPhaseIdx);
+                    (*poro)[globalIdx] = elemVolVars[scvIdx].porosity();
+                    (*Te)[globalIdx] = elemVolVars[scvIdx].temperature();
+                }
+
+                // velocity output
+                velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, *eIt, /*phaseIdx=*/0);
             }
-
-            // velocity output
-            velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, *eIt, /*phaseIdx=*/0);
         }
 
         writer.attachDofData(*sn, "Sn", isBox);