From 30fa9a5fc036b02b1db6846590fce0b060b34051 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Fri, 4 May 2012 08:50:36 +0000
Subject: [PATCH] Prepare box models for generalization to cell centered: -
 Inherit the model specific LocalResidual and Model from the properties  
 BaseLocalResidual and BaseModel instead of BoxLocalResidual and   BoxModel.

- Replace fvGeometry.numVertices by fvGeometry.numFAP for the flux
  calculations, or by fvGeometry.numSCV where appropriate.

Reviewed by Christoph.


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@8223 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
 dumux/boxmodels/1p/1pfluxvariables.hh         |  7 ++-
 dumux/boxmodels/1p/1plocalresidual.hh         |  2 +-
 dumux/boxmodels/1p/1pmodel.hh                 |  2 +-
 dumux/boxmodels/1p2c/1p2cfluxvariables.hh     | 19 ++++----
 dumux/boxmodels/1p2c/1p2clocalresidual.hh     |  2 +-
 dumux/boxmodels/1p2c/1p2cmodel.hh             |  2 +-
 dumux/boxmodels/2p/2pfluxvariables.hh         |  4 +-
 dumux/boxmodels/2p2c/2p2cfluxvariables.hh     | 24 +++++++----
 dumux/boxmodels/2p2c/2p2clocalresidual.hh     |  4 +-
 dumux/boxmodels/2p2c/2p2cmodel.hh             |  6 +--
 dumux/boxmodels/2p2cni/2p2cnifluxvariables.hh |  8 +++-
 dumux/boxmodels/2pni/2pnifluxvariables.hh     |  8 +++-
 dumux/boxmodels/3p3c/3p3cfluxvariables.hh     | 43 ++++++++++---------
 dumux/boxmodels/3p3c/3p3clocalresidual.hh     |  2 +-
 dumux/boxmodels/3p3c/3p3cmodel.hh             |  6 +--
 dumux/boxmodels/3p3cni/3p3cnifluxvariables.hh |  8 +++-
 .../boxmodels/common/boxfvelementgeometry.hh  |  2 +
 .../mpnc/energy/mpncfluxvariablesenergy.hh    |  8 +++-
 dumux/boxmodels/mpnc/mpncfluxvariables.hh     |  7 ++-
 dumux/boxmodels/mpnc/mpnclocalresidual.hh     |  8 ++--
 dumux/boxmodels/mpnc/mpncmodel.hh             |  4 +-
 .../richards/richardsfluxvariables.hh         |  7 ++-
 .../richards/richardslocalresidual.hh         |  2 +-
 dumux/boxmodels/richards/richardsmodel.hh     |  2 +-
 .../richards/richardsnewtoncontroller.hh      | 10 ++---
 25 files changed, 118 insertions(+), 79 deletions(-)

diff --git a/dumux/boxmodels/1p/1pfluxvariables.hh b/dumux/boxmodels/1p/1pfluxvariables.hh
index 67bc7d9e38..4b5b5be375 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 9d50acc139..f8ad2ed03d 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 6a2299fbbf..d1d2306f40 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 c89d83adb7..21c3dd7b26 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 9dbf104ec8..cbb43c182e 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 88ef8f1114..bd779f07f8 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 4040e6a4b7..72d4dbd38b 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 07bcf16efb..2387d0e730 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 23cae44075..959b25d509 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 3bc1061b1e..8f57fedb1a 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 abc1c5df8f..144dc2064c 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 4a4f084e1f..a45313b30f 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 d17133f8fa..4b948d810e 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 de85a64c31..4cf1092bea 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 81b059f814..3b454ca9c6 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 7bd418cff6..1750eca01e 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 bd17b5509b..9e5ec2e9bb 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 5cc0882a3f..3cbf49610a 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 797a1cd54c..0234c8b9b6 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 5b214216fb..de8594843e 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 ef7e431946..d5989ac261 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 6f5aac48ec..fdfdce4931 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 3634770413..e1c8634dcf 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 36296581d5..b1ce78b5eb 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 ee2b7dddf5..5bf5ed5467 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));
-- 
GitLab