diff --git a/dumux/boxmodels/1p/1pfluxvariables.hh b/dumux/boxmodels/1p/1pfluxvariables.hh
index 67bc7d9e3878d31aad51a8ebaec10de44430effe..4b5b5be3756337225c5dc8e4cf0ae81fe690a1e3 100644
--- a/dumux/boxmodels/1p/1pfluxvariables.hh
+++ b/dumux/boxmodels/1p/1pfluxvariables.hh
@@ -148,15 +148,18 @@ private:
 
         // calculate potential gradient
         for (int idx = 0;
-             idx < fvGeometry_.numVertices;
+             idx < fvGeometry_.numFAP;
              idx++) // loop over adjacent vertices
         {
             // FE gradient at vertex idx
             const DimVector &feGrad = face().grad[idx];
 
+            // index for the element volume variables 
+            int volVarsIdx = face().fapIndices[idx];
+
             // the pressure gradient
             DimVector tmp(feGrad);
-            tmp *= elemVolVars[idx].pressure();
+            tmp *= elemVolVars[volVarsIdx].pressure();
             potentialGrad_ += tmp;
         }
 
diff --git a/dumux/boxmodels/1p/1plocalresidual.hh b/dumux/boxmodels/1p/1plocalresidual.hh
index 9d50acc139d39d165cf125a1e77620ed57e83f31..f8ad2ed03dfc368df1b31ac92407f269f97ee1ba 100644
--- a/dumux/boxmodels/1p/1plocalresidual.hh
+++ b/dumux/boxmodels/1p/1plocalresidual.hh
@@ -44,7 +44,7 @@ namespace Dumux
  *        using the one-phase box model.
  */
 template<class TypeTag>
-class OnePLocalResidual : public BoxLocalResidual<TypeTag>
+class OnePLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
 {
     typedef OnePLocalResidual<TypeTag> ThisType;
 
diff --git a/dumux/boxmodels/1p/1pmodel.hh b/dumux/boxmodels/1p/1pmodel.hh
index 6a2299fbbfd615443f13e58a93aaddc8d3108f0d..d1d2306f408be95d8c6a65074b347258429cb1a2 100644
--- a/dumux/boxmodels/1p/1pmodel.hh
+++ b/dumux/boxmodels/1p/1pmodel.hh
@@ -52,7 +52,7 @@ namespace Dumux
  * model supports compressible as well as incompressible fluids.
  */
 template<class TypeTag >
-class OnePBoxModel : public BoxModel<TypeTag>
+class OnePBoxModel : public GET_PROP_TYPE(TypeTag, BaseModel)
 {
     typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
     typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
diff --git a/dumux/boxmodels/1p2c/1p2cfluxvariables.hh b/dumux/boxmodels/1p2c/1p2cfluxvariables.hh
index c89d83adb7ab506f17cd527886767bb47af678d5..21c3dd7b260dae4e1719a5affbb94593690a580a 100644
--- a/dumux/boxmodels/1p2c/1p2cfluxvariables.hh
+++ b/dumux/boxmodels/1p2c/1p2cfluxvariables.hh
@@ -345,40 +345,43 @@ protected:
             // use finite-element gradients
             tmp = 0.0;
             for (int idx = 0;
-                    idx < fvGeometry_.numVertices;
+                    idx < fvGeometry_.numFAP;
                     idx++) // loop over adjacent vertices
             {
                 // FE gradient at vertex idx
                 const DimVector &feGrad = face().grad[idx];
 
+                // index for the element volume variables 
+                int volVarsIdx = face().fapIndices[idx];
+
                 // the pressure gradient
                 tmp = feGrad;
-                tmp *= elemVolVars[idx].pressure();
+                tmp *= elemVolVars[volVarsIdx].pressure();
                 potentialGrad_ += tmp;
 
                 // the concentration gradient [mol/m^3/m]
                 tmp = feGrad;
-                tmp *= elemVolVars[idx].molarity(comp1Idx);
+                tmp *= elemVolVars[volVarsIdx].molarity(comp1Idx);
                 concentrationGrad_ += tmp;
 
                 // the mole-fraction gradient
                 tmp = feGrad;
-                tmp *= elemVolVars[idx].moleFraction(comp1Idx);
+                tmp *= elemVolVars[volVarsIdx].moleFraction(comp1Idx);
                 moleFractionGrad_ += tmp;
 
                 // the mass-fraction gradient
                 tmp = feGrad;
-                tmp *= elemVolVars[idx].massFraction(comp1Idx);
+                tmp *= elemVolVars[volVarsIdx].massFraction(comp1Idx);
                 massFractionGrad_ += tmp;
 
                 // phase viscosity
-                viscosity_ += elemVolVars[idx].viscosity()*face().shapeValue[idx];
+                viscosity_ += elemVolVars[volVarsIdx].viscosity()*face().shapeValue[idx];
 
                 //phase molar density
-                molarDensity_ += elemVolVars[idx].molarDensity()*face().shapeValue[idx];
+                molarDensity_ += elemVolVars[volVarsIdx].molarDensity()*face().shapeValue[idx];
 
                 //phase density
-                density_ += elemVolVars[idx].density()*face().shapeValue[idx];
+                density_ += elemVolVars[volVarsIdx].density()*face().shapeValue[idx];
             }
         }
         else {
diff --git a/dumux/boxmodels/1p2c/1p2clocalresidual.hh b/dumux/boxmodels/1p2c/1p2clocalresidual.hh
index 9dbf104ec8f084edb3f4826717e6c0ce4d19a130..cbb43c182edd029f2f69966bc8eee03e74857fe5 100644
--- a/dumux/boxmodels/1p2c/1p2clocalresidual.hh
+++ b/dumux/boxmodels/1p2c/1p2clocalresidual.hh
@@ -55,7 +55,7 @@ namespace Dumux
  *  This class is used to fill the gaps in BoxLocalResidual for the 1p2c flow and transport.
  */
 template<class TypeTag>
-class OnePTwoCLocalResidual : public BoxLocalResidual<TypeTag>
+class OnePTwoCLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
 {
 protected:
     typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
diff --git a/dumux/boxmodels/1p2c/1p2cmodel.hh b/dumux/boxmodels/1p2c/1p2cmodel.hh
index 88ef8f1114993657d49d784dfa5ebaab5fa976c6..bd779f07f8dee2d0f535ea0d99cd724d450df1ff 100644
--- a/dumux/boxmodels/1p2c/1p2cmodel.hh
+++ b/dumux/boxmodels/1p2c/1p2cmodel.hh
@@ -75,7 +75,7 @@ namespace Dumux
  */
 
 template<class TypeTag >
-class OnePTwoCBoxModel : public BoxModel<TypeTag>
+class OnePTwoCBoxModel : public GET_PROP_TYPE(TypeTag, BaseModel)
 {
     typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
     typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
diff --git a/dumux/boxmodels/2p/2pfluxvariables.hh b/dumux/boxmodels/2p/2pfluxvariables.hh
index 4040e6a4b72fa46a0fb76f24a4eace48626e2b2f..72d4dbd38be598fc988ad70d7ea139c20c6f5654 100644
--- a/dumux/boxmodels/2p/2pfluxvariables.hh
+++ b/dumux/boxmodels/2p/2pfluxvariables.hh
@@ -172,8 +172,8 @@ private:
             // FE gradient at vertex idx
             const Vector &feGrad = face().grad[idx];
 
-	    // index for the element volume variables 
-	    int volVarsIdx = face().fapIndices[idx];
+            // index for the element volume variables 
+            int volVarsIdx = face().fapIndices[idx];
 	    
             // compute sum of pressure gradients for each phase
             for (int phase = 0; phase < numPhases; phase++)
diff --git a/dumux/boxmodels/2p2c/2p2cfluxvariables.hh b/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
index 07bcf16efb8f23e74f168acf8cacf27957911b4a..2387d0e7308a021e82b056ebff6fdfdff69b7a9c 100644
--- a/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
+++ b/dumux/boxmodels/2p2c/2p2cfluxvariables.hh
@@ -115,14 +115,17 @@ protected:
         // calculate densities at the integration points of the face
         DimVector tmp(0.0);
         for (int idx = 0;
-             idx < fvGeometry_.numVertices;
+             idx < fvGeometry_.numFAP;
              idx++) // loop over adjacent vertices
         {
+            // index for the element volume variables 
+            int volVarsIdx = face().fapIndices[idx];
+
             for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
             {
-                density_[phaseIdx] += elemVolVars[idx].density(phaseIdx)*
+                density_[phaseIdx] += elemVolVars[volVarsIdx].density(phaseIdx)*
                     face().shapeValue[idx];
-                molarDensity_[phaseIdx] += elemVolVars[idx].molarDensity(phaseIdx)*
+                molarDensity_[phaseIdx] += elemVolVars[volVarsIdx].molarDensity(phaseIdx)*
                     face().shapeValue[idx];
             }
         }
@@ -139,39 +142,42 @@ protected:
         // calculate gradients
         DimVector tmp(0.0);
         for (int idx = 0;
-             idx < fvGeometry_.numVertices;
+             idx < fvGeometry_.numFAP;
              idx++) // loop over adjacent vertices
         {
             // FE gradient at vertex idx
             const DimVector &feGrad = face().grad[idx];
 
+            // index for the element volume variables 
+            int volVarsIdx = face().fapIndices[idx];
+
             // compute sum of pressure gradients for each phase
             for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
             {
                 // the pressure gradient
                 tmp = feGrad;
-                tmp *= elemVolVars[idx].pressure(phaseIdx);
+                tmp *= elemVolVars[volVarsIdx].pressure(phaseIdx);
                 potentialGrad_[phaseIdx] += tmp;
             }
 
             // the concentration gradient of the non-wetting
             // component in the wetting phase
             tmp = feGrad;
-            tmp *= elemVolVars[idx].fluidState().massFraction(wPhaseIdx, nCompIdx);
+            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(wPhaseIdx, nCompIdx);
             massFractionGrad_[wPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemVolVars[idx].fluidState().moleFraction(wPhaseIdx, nCompIdx);
+            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(wPhaseIdx, nCompIdx);
             moleFractionGrad_[wPhaseIdx] += tmp;
 
             //            // the concentration gradient of the wetting component
             //            // in the non-wetting phase
             tmp = feGrad;
-            tmp *= elemVolVars[idx].fluidState().massFraction(nPhaseIdx, wCompIdx);
+            tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(nPhaseIdx, wCompIdx);
             massFractionGrad_[nPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemVolVars[idx].fluidState().moleFraction(nPhaseIdx, wCompIdx);
+            tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(nPhaseIdx, wCompIdx);
             moleFractionGrad_[nPhaseIdx] += tmp;
         }
 
diff --git a/dumux/boxmodels/2p2c/2p2clocalresidual.hh b/dumux/boxmodels/2p2c/2p2clocalresidual.hh
index 23cae4407587d6b3ff292c3148ac400448cd317e..959b25d50933677b5abc2877946a0e2937b491a1 100644
--- a/dumux/boxmodels/2p2c/2p2clocalresidual.hh
+++ b/dumux/boxmodels/2p2c/2p2clocalresidual.hh
@@ -121,7 +121,7 @@ class TwoPTwoCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
         ElementVolumeVariables elemVolVars;
         elemVolVars.update(this->problem_(), element, fvGeometry, false);
 
-        this->storageTerm_.resize(fvGeometry.numVertices);
+        this->storageTerm_.resize(fvGeometry.numSCV);
         this->storageTerm_ = 0;
 
         this->elemPtr_ = &element;
@@ -365,7 +365,7 @@ class TwoPTwoCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
     void evalPhaseStorage_(const int phaseIdx)
     {
         // evaluate the storage terms of a single phase
-        for (int i=0; i < this->fvGeometry_().numVertices; i++) {
+        for (int i=0; i < this->fvGeometry_().numSCV; i++) {
             PrimaryVariables &storage = this->storageTerm_[i];
             const ElementVolumeVariables &elemVolVars = this->curVolVars_();
             const VolumeVariables &volVars = elemVolVars[i];
diff --git a/dumux/boxmodels/2p2c/2p2cmodel.hh b/dumux/boxmodels/2p2c/2p2cmodel.hh
index 3bc1061b1e1af713588ccd9953a48558afa5629f..8f57fedb1a9d547d0250eb740979a21f28d0e384 100644
--- a/dumux/boxmodels/2p2c/2p2cmodel.hh
+++ b/dumux/boxmodels/2p2c/2p2cmodel.hh
@@ -87,9 +87,9 @@ namespace Dumux
  */
 
 template<class TypeTag>
-class TwoPTwoCModel: public BoxModel<TypeTag>
+class TwoPTwoCModel: public GET_PROP_TYPE(TypeTag, BaseModel)
 {
-    typedef BoxModel<TypeTag> ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, BaseModel) ParentType;
 
     typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
     typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
@@ -584,7 +584,7 @@ public:
         for (; eIt != eEndIt; ++eIt)
         {
             fvGeometry.update(this->gridView_(), *eIt);
-            for (int scvIdx = 0; scvIdx < fvGeometry.numVertices; ++scvIdx)
+            for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
             {
                 int globalIdx = this->vertexMapper().map(*eIt, scvIdx, dim);
 
diff --git a/dumux/boxmodels/2p2cni/2p2cnifluxvariables.hh b/dumux/boxmodels/2p2cni/2p2cnifluxvariables.hh
index abc1c5df8ff5f5b908e1c9efff818d8dc7fe59ba..144dc2064ca5bdcb3ae767606f71abf31039d6ac 100644
--- a/dumux/boxmodels/2p2cni/2p2cnifluxvariables.hh
+++ b/dumux/boxmodels/2p2cni/2p2cnifluxvariables.hh
@@ -108,10 +108,14 @@ protected:
         // gradients
         temperatureGrad_ = 0;
         DimVector tmp(0.0);
-        for (int idx = 0; idx < this->fvGeometry_.numVertices; idx++)
+        for (int idx = 0; idx < this->fvGeometry_.numFAP; idx++)
         {
             tmp = this->face().grad[idx];
-            tmp *= elemVolVars[idx].temperature();
+
+            // index for the element volume variables 
+            int volVarsIdx = this->face().fapIndices[idx];
+
+            tmp *= elemVolVars[volVarsIdx].temperature();
             temperatureGrad_ += tmp;
         }
 
diff --git a/dumux/boxmodels/2pni/2pnifluxvariables.hh b/dumux/boxmodels/2pni/2pnifluxvariables.hh
index 4a4f084e1f4b70f25e8aaddda405016bff83f04a..a45313b30ff1a135c0247cc81a2868d9dfe21d61 100644
--- a/dumux/boxmodels/2pni/2pnifluxvariables.hh
+++ b/dumux/boxmodels/2pni/2pnifluxvariables.hh
@@ -93,10 +93,14 @@ public:
         // gradients
         Vector temperatureGrad(0);
         Vector tmp(0.0);
-        for (int vertIdx = 0; vertIdx < elemGeom.numVertices; vertIdx++)
+        for (int vertIdx = 0; vertIdx < elemGeom.numFAP; vertIdx++)
         {
             tmp = this->face().grad[vertIdx];
-            tmp *= elemDat[vertIdx].temperature();
+
+            // index for the element volume variables 
+            int volVarsIdx = this->face().fapIndices[vertIdx];
+            
+            tmp *= elemDat[volVarsIdx].temperature();
             temperatureGrad += tmp;
         }
 
diff --git a/dumux/boxmodels/3p3c/3p3cfluxvariables.hh b/dumux/boxmodels/3p3c/3p3cfluxvariables.hh
index d17133f8fa599c2d02e253e65524f7320d9062ed..4b948d810eae3ab807462ca7ef8698d06453ccb8 100644
--- a/dumux/boxmodels/3p3c/3p3cfluxvariables.hh
+++ b/dumux/boxmodels/3p3c/3p3cfluxvariables.hh
@@ -135,95 +135,98 @@ private:
         // calculate gradients
         Vector tmp(0.0);
         for (int idx = 0;
-             idx < fvElemGeom_.numVertices;
+             idx < fvElemGeom_.numFAP;
              idx++) // loop over adjacent vertices
         {
             // FE gradient at vertex idx
             const Vector &feGrad = face().grad[idx];
 
+            // index for the element volume variables 
+            int volVarsIdx = face().fapIndices[idx];
+
             // compute sum of pressure gradients for each phase
             for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
             {
                 // the pressure gradient
                 tmp = feGrad;
-                tmp *= elemDat[idx].pressure(phaseIdx);
+                tmp *= elemDat[volVarsIdx].pressure(phaseIdx);
                 potentialGrad_[phaseIdx] += tmp;
             }
 
             // the concentration gradient of the components
             // component in the phases
             tmp = feGrad;
-            tmp *= elemDat[idx].fluidState().moleFraction(wPhaseIdx, wCompIdx);
+            tmp *= elemDat[volVarsIdx].fluidState().moleFraction(wPhaseIdx, wCompIdx);
             wConcentrationGrad_[wPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemDat[idx].fluidState().moleFraction(nPhaseIdx, wCompIdx);
+            tmp *= elemDat[volVarsIdx].fluidState().moleFraction(nPhaseIdx, wCompIdx);
             wConcentrationGrad_[nPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemDat[idx].fluidState().moleFraction(gPhaseIdx, wCompIdx);
+            tmp *= elemDat[volVarsIdx].fluidState().moleFraction(gPhaseIdx, wCompIdx);
             wConcentrationGrad_[gPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemDat[idx].fluidState().moleFraction(wPhaseIdx, cCompIdx);
+            tmp *= elemDat[volVarsIdx].fluidState().moleFraction(wPhaseIdx, cCompIdx);
             cConcentrationGrad_[wPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemDat[idx].fluidState().moleFraction(nPhaseIdx, cCompIdx);
+            tmp *= elemDat[volVarsIdx].fluidState().moleFraction(nPhaseIdx, cCompIdx);
             cConcentrationGrad_[nPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemDat[idx].fluidState().moleFraction(gPhaseIdx, cCompIdx);
+            tmp *= elemDat[volVarsIdx].fluidState().moleFraction(gPhaseIdx, cCompIdx);
             cConcentrationGrad_[gPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemDat[idx].fluidState().moleFraction(wPhaseIdx, aCompIdx);
+            tmp *= elemDat[volVarsIdx].fluidState().moleFraction(wPhaseIdx, aCompIdx);
             aConcentrationGrad_[wPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemDat[idx].fluidState().moleFraction(nPhaseIdx, aCompIdx);
+            tmp *= elemDat[volVarsIdx].fluidState().moleFraction(nPhaseIdx, aCompIdx);
             aConcentrationGrad_[nPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemDat[idx].fluidState().moleFraction(gPhaseIdx, aCompIdx);
+            tmp *= elemDat[volVarsIdx].fluidState().moleFraction(gPhaseIdx, aCompIdx);
             aConcentrationGrad_[gPhaseIdx] += tmp;
 
             // the molar concentration gradients of the components
             // in the phases
             tmp = feGrad;
-            tmp *= elemDat[idx].fluidState().moleFraction(wPhaseIdx, wCompIdx);
+            tmp *= elemDat[volVarsIdx].fluidState().moleFraction(wPhaseIdx, wCompIdx);
             molarWConcGrad_[wPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemDat[idx].fluidState().moleFraction(nPhaseIdx, wCompIdx);
+            tmp *= elemDat[volVarsIdx].fluidState().moleFraction(nPhaseIdx, wCompIdx);
             molarWConcGrad_[nPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemDat[idx].fluidState().moleFraction(gPhaseIdx, wCompIdx);
+            tmp *= elemDat[volVarsIdx].fluidState().moleFraction(gPhaseIdx, wCompIdx);
             molarWConcGrad_[gPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemDat[idx].fluidState().moleFraction(wPhaseIdx, cCompIdx);
+            tmp *= elemDat[volVarsIdx].fluidState().moleFraction(wPhaseIdx, cCompIdx);
             molarCConcGrad_[wPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemDat[idx].fluidState().moleFraction(nPhaseIdx, cCompIdx);
+            tmp *= elemDat[volVarsIdx].fluidState().moleFraction(nPhaseIdx, cCompIdx);
             molarCConcGrad_[nPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemDat[idx].fluidState().moleFraction(gPhaseIdx, cCompIdx);
+            tmp *= elemDat[volVarsIdx].fluidState().moleFraction(gPhaseIdx, cCompIdx);
             molarCConcGrad_[gPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemDat[idx].fluidState().moleFraction(wPhaseIdx, aCompIdx);
+            tmp *= elemDat[volVarsIdx].fluidState().moleFraction(wPhaseIdx, aCompIdx);
             molarAConcGrad_[wPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemDat[idx].fluidState().moleFraction(nPhaseIdx, aCompIdx);
+            tmp *= elemDat[volVarsIdx].fluidState().moleFraction(nPhaseIdx, aCompIdx);
             molarAConcGrad_[nPhaseIdx] += tmp;
 
             tmp = feGrad;
-            tmp *= elemDat[idx].fluidState().moleFraction(gPhaseIdx, aCompIdx);
+            tmp *= elemDat[volVarsIdx].fluidState().moleFraction(gPhaseIdx, aCompIdx);
             molarAConcGrad_[gPhaseIdx] += tmp;
         }
 
diff --git a/dumux/boxmodels/3p3c/3p3clocalresidual.hh b/dumux/boxmodels/3p3c/3p3clocalresidual.hh
index de85a64c31730c7f030d3bb9a4347620cb491b41..4cf1092bea2dc68f2fe7c2d19a0d8aaa0f84b805 100644
--- a/dumux/boxmodels/3p3c/3p3clocalresidual.hh
+++ b/dumux/boxmodels/3p3c/3p3clocalresidual.hh
@@ -55,7 +55,7 @@ namespace Dumux
  * This class is used to fill the gaps in BoxLocalResidual for the 3P3C flow.
  */
 template<class TypeTag>
-class ThreePThreeCLocalResidual: public BoxLocalResidual<TypeTag>
+class ThreePThreeCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
 {
 protected:
     typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
diff --git a/dumux/boxmodels/3p3c/3p3cmodel.hh b/dumux/boxmodels/3p3c/3p3cmodel.hh
index 81b059f81473a0108fc4fb16d6955a170c54397b..3b454ca9c6d9d29f297893978f8c4e75654404ec 100644
--- a/dumux/boxmodels/3p3c/3p3cmodel.hh
+++ b/dumux/boxmodels/3p3c/3p3cmodel.hh
@@ -95,9 +95,9 @@ namespace Dumux
  * </ul>
  */
 template<class TypeTag>
-class ThreePThreeCModel: public BoxModel<TypeTag>
+class ThreePThreeCModel: public GET_PROP_TYPE(TypeTag, BaseModel)
 {
-    typedef BoxModel<TypeTag> ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, BaseModel) ParentType;
 
     typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
     typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
@@ -445,7 +445,7 @@ public:
         for (; it != endit; ++it)
         {
             fvElemGeom.update(this->gridView_(), *it);
-            for (int i = 0; i < fvElemGeom.numVertices; ++i)
+            for (int i = 0; i < fvElemGeom.numSCV; ++i)
             {
                 int globalIdx = this->vertexMapper().map(*it, i, dim);
 
diff --git a/dumux/boxmodels/3p3cni/3p3cnifluxvariables.hh b/dumux/boxmodels/3p3cni/3p3cnifluxvariables.hh
index 7bd418cff65073736d9cc3fd9e854a564856675a..1750eca01efd9d3fd86d1abe0f842fa66cea64d0 100644
--- a/dumux/boxmodels/3p3cni/3p3cnifluxvariables.hh
+++ b/dumux/boxmodels/3p3cni/3p3cnifluxvariables.hh
@@ -96,10 +96,14 @@ public:
         // gradients
         Vector temperatureGrad(0);
         Vector tmp(0.0);
-        for (int vertIdx = 0; vertIdx < elemGeom.numVertices; vertIdx++)
+        for (int vertIdx = 0; vertIdx < elemGeom.numFAP; vertIdx++)
         {
             tmp = this->face().grad[vertIdx];
-            tmp *= elemDat[vertIdx].temperature();
+
+            // index for the element volume variables 
+            int volVarsIdx = this->face().fapIndices[vertIdx];
+
+            tmp *= elemDat[volVarsIdx].temperature();
             temperatureGrad += tmp;
         }
 
diff --git a/dumux/boxmodels/common/boxfvelementgeometry.hh b/dumux/boxmodels/common/boxfvelementgeometry.hh
index bd17b5509b29620072b9037c886bc33188b8214d..9e5ec2e9bbd3e9fb768f25cb2398ca1369ee0b46 100644
--- a/dumux/boxmodels/common/boxfvelementgeometry.hh
+++ b/dumux/boxmodels/common/boxfvelementgeometry.hh
@@ -588,6 +588,7 @@ public:
     int numVertices; //!< number of verts
     int numEdges; //!< number of edges
     int numFaces; //!< number of faces (0 in < 3D)
+    int numSCV; //!< number of subcontrol volumes
     int numFAP; //!< number of flux approximation points
 
     const LocalFiniteElementCache feCache_;
@@ -621,6 +622,7 @@ public:
         numVertices = referenceElement.size(dim);
         numEdges = referenceElement.size(dim-1);
         numFaces = (dim<3)?0:referenceElement.size(1);
+        numSCV = numVertices;
         numFAP = numVertices;
 
         // corners:
diff --git a/dumux/boxmodels/mpnc/energy/mpncfluxvariablesenergy.hh b/dumux/boxmodels/mpnc/energy/mpncfluxvariablesenergy.hh
index 5cc0882a3fbe803fbe20f821b04baabf24f895ef..3cbf49610ad4f5a1a71e82d850eb3f4a21fa4492 100644
--- a/dumux/boxmodels/mpnc/energy/mpncfluxvariablesenergy.hh
+++ b/dumux/boxmodels/mpnc/energy/mpncfluxvariablesenergy.hh
@@ -110,10 +110,14 @@ public:
         // gradients
         DimVector tmp(0.0);
         DimVector temperatureGradient(0.);
-        for (int scvIdx = 0; scvIdx < fvGeometry.numVertices; scvIdx++)
+        for (int scvIdx = 0; scvIdx < fvGeometry.numFAP; scvIdx++)
         {
             tmp = this->face().grad[scvIdx];
-            tmp *= elemVolVars[scvIdx].fluidState().temperature(/*phaseIdx=*/0);
+
+            // index for the element volume variables 
+            int volVarsIdx = this->face().fapIndices[scvIdx];
+            
+            tmp *= elemVolVars[volVarsIdx].fluidState().temperature(/*phaseIdx=*/0);
             temperatureGradient += tmp;
         }
 
diff --git a/dumux/boxmodels/mpnc/mpncfluxvariables.hh b/dumux/boxmodels/mpnc/mpncfluxvariables.hh
index 797a1cd54cbea7db055fd5b9ce889bf24486a684..0234c8b9b64a36886414f4937a6446e4e2db0407 100644
--- a/dumux/boxmodels/mpnc/mpncfluxvariables.hh
+++ b/dumux/boxmodels/mpnc/mpncfluxvariables.hh
@@ -248,12 +248,15 @@ private:
         // calculate pressure gradients using finite element gradients
         DimVector tmp(0.0);
         for (int idx = 0;
-             idx < fvGeometry_.numVertices;
+             idx < fvGeometry_.numFAP;
              idx++) // loop over adjacent vertices
         {
             // FE gradient at vertex idx
             const DimVector &feGrad = face().grad[idx];
 
+            // index for the element volume variables 
+            int volVarsIdx = face().fapIndices[idx];
+
             // TODO: only calculate the gradients for the present
             // phases.
             //
@@ -262,7 +265,7 @@ private:
             {
                 // the pressure gradient
                 tmp = feGrad;
-                tmp *= elemVolVars[idx].fluidState().pressure(phaseIdx);
+                tmp *= elemVolVars[volVarsIdx].fluidState().pressure(phaseIdx);
                 potentialGrad_[phaseIdx] += tmp;
             }
         }
diff --git a/dumux/boxmodels/mpnc/mpnclocalresidual.hh b/dumux/boxmodels/mpnc/mpnclocalresidual.hh
index 5b214216fb60f0f62d100134cd0fc6de6eb7809e..de8594843e380c5fb1d1b6f87b312f711792d779 100644
--- a/dumux/boxmodels/mpnc/mpnclocalresidual.hh
+++ b/dumux/boxmodels/mpnc/mpnclocalresidual.hh
@@ -42,13 +42,13 @@ namespace Dumux
  * two-phase, N-component twophase flow.
  */
 template<class TypeTag>
-class MPNCLocalResidual : public BoxLocalResidual<TypeTag>
+class MPNCLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
 {
     typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
 
 protected:
     typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
-    typedef BoxLocalResidual<TypeTag> ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, BaseLocalResidual) ParentType;
     typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
 
     enum {numPhases = GET_PROP_VALUE(TypeTag, NumPhases)};
@@ -120,7 +120,7 @@ public:
 
         // calculate the phase storage for all sub-control volumes
         for (int scvIdx=0;
-             scvIdx < fvGeometry.numVertices;
+             scvIdx < fvGeometry.numSCV;
              scvIdx++)
         {
             PrimaryVariables tmp(0.0);
@@ -224,7 +224,7 @@ public:
                          curVolVars,
                          bcType);
 
-        for (int i = 0; i < this->fvGeometry_().numVertices; ++i) {
+        for (int i = 0; i < this->fvGeometry_().numSCV; ++i) {
             // add the two auxiliary equations, make sure that the
             // dirichlet boundary condition is conserved
             for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
diff --git a/dumux/boxmodels/mpnc/mpncmodel.hh b/dumux/boxmodels/mpnc/mpncmodel.hh
index ef7e43194600d2fc1cc3b106003fe2edcf30653f..d5989ac261c2a71fe3cc655850a185050e7729ff 100644
--- a/dumux/boxmodels/mpnc/mpncmodel.hh
+++ b/dumux/boxmodels/mpnc/mpncmodel.hh
@@ -112,9 +112,9 @@ namespace Dumux
  * - Temperature \f$T\f$ if the energy equation is enabled
  */
 template<class TypeTag>
-class MPNCModel : public BoxModel<TypeTag>
+class MPNCModel : public GET_PROP_TYPE(TypeTag, BaseModel)
 {
-    typedef BoxModel<TypeTag> ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, BaseModel) ParentType;
     typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
     typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
     typedef typename GridView::template Codim<0>::Entity Element;
diff --git a/dumux/boxmodels/richards/richardsfluxvariables.hh b/dumux/boxmodels/richards/richardsfluxvariables.hh
index 6f5aac48ec40b793a90e44e04aced8bcabe844ac..fdfdce4931983fad5b9a86f986a8d170ccfe1087 100644
--- a/dumux/boxmodels/richards/richardsfluxvariables.hh
+++ b/dumux/boxmodels/richards/richardsfluxvariables.hh
@@ -152,15 +152,18 @@ protected:
         potentialGrad_ = 0.0;
         // calculate gradients
         for (int idx = 0;
-             idx < fvGeometry_.numVertices;
+             idx < fvGeometry_.numFAP;
              idx++) // loop over adjacent vertices
         {
             // FE gradient at vertex index
             const DimVector &feGrad = face().grad[idx];
 
+            // index for the element volume variables 
+            int volVarsIdx = face().fapIndices[idx];
+
             // the pressure gradient
             DimVector tmp(feGrad);
-            tmp *= elemVolVars[idx].pressure(wPhaseIdx);
+            tmp *= elemVolVars[volVarsIdx].pressure(wPhaseIdx);
 
             potentialGrad_ += tmp;
         }
diff --git a/dumux/boxmodels/richards/richardslocalresidual.hh b/dumux/boxmodels/richards/richardslocalresidual.hh
index 36347704139a11a763e89c0e10abc10420ace7c6..e1c8634dcfbfca993b79af902176a58f3d7950ef 100644
--- a/dumux/boxmodels/richards/richardslocalresidual.hh
+++ b/dumux/boxmodels/richards/richardslocalresidual.hh
@@ -42,7 +42,7 @@ namespace Dumux
  * \brief Element-wise calculation of the residual for the Richards box model.
  */
 template<class TypeTag>
-class RichardsLocalResidual : public BoxLocalResidual<TypeTag>
+class RichardsLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
 {
     typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
     typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
diff --git a/dumux/boxmodels/richards/richardsmodel.hh b/dumux/boxmodels/richards/richardsmodel.hh
index 36296581d56a682d2dfb3942420b41c91af4a89b..b1ce78b5ebe43f25836eaf160c81e84b3e9903f1 100644
--- a/dumux/boxmodels/richards/richardsmodel.hh
+++ b/dumux/boxmodels/richards/richardsmodel.hh
@@ -92,7 +92,7 @@ namespace Dumux
  * Richards model!
  */
 template<class TypeTag >
-class RichardsModel : public BoxModel<TypeTag>
+class RichardsModel : public GET_PROP_TYPE(TypeTag, BaseModel)
 {
     typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
     typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
diff --git a/dumux/boxmodels/richards/richardsnewtoncontroller.hh b/dumux/boxmodels/richards/richardsnewtoncontroller.hh
index ee2b7dddf5175b96a7f8e5b6f4479d6837848970..5bf5ed54672f3816df57cce2734c0c0344fcd52a 100644
--- a/dumux/boxmodels/richards/richardsnewtoncontroller.hh
+++ b/dumux/boxmodels/richards/richardsnewtoncontroller.hh
@@ -90,21 +90,21 @@ public:
                 return;
 
             // clamp saturation change to at most 20% per iteration
-            FVElementGeometry fvGeomtry;
+            FVElementGeometry fvGeometry;
             const GridView &gridView = this->problem_().gridView();
             ElementIterator elemIt = gridView.template begin<0>();
             const ElementIterator elemEndIt = gridView.template end<0>();
             for (; elemIt != elemEndIt; ++elemIt) {
-                fvGeomtry.update(gridView, *elemIt);
-                for (int i = 0; i < fvGeomtry.numVertices; ++i) {
+                fvGeometry.update(gridView, *elemIt);
+                for (int i = 0; i < fvGeometry.numSCV; ++i) {
                     int globI = this->problem_().vertexMapper().map(*elemIt, i, dim);
 
                     // calculate the old wetting phase saturation
                     const SpatialParams &spatialParams = this->problem_().spatialParams();
-                    const MaterialLawParams &mp = spatialParams.materialLawParams(*elemIt, fvGeomtry, i);
+                    const MaterialLawParams &mp = spatialParams.materialLawParams(*elemIt, fvGeometry, i);
                     Scalar pcMin = MaterialLaw::pC(mp, 1.0);
                     Scalar pW = uLastIter[globI][pwIdx];
-                    Scalar pN = std::max(this->problem_().referencePressure(*elemIt, fvGeomtry, i),
+                    Scalar pN = std::max(this->problem_().referencePressure(*elemIt, fvGeometry, i),
                                          pW + pcMin);
                     Scalar pcOld = pN - pW;
                     Scalar SwOld = std::max<Scalar>(0.0, MaterialLaw::Sw(mp, pcOld));