From ca451d24f8be4fa6c3703311eec413d6bc1bbfba Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Wed, 28 Nov 2012 15:49:35 +0000
Subject: [PATCH] implicit branch: further unify the treatment of
 intrinsicPermeability in the problem-specific spatial parameters. This comes
 at the expense of adding an additional method to the BoxSpatialParams1P base
 class and a separated treatment in the flux variables.

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/branches/implicit@9703 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
 dumux/implicit/1p2c/1p2cfluxvariables.hh      | 24 ++++++---
 dumux/implicit/box/boxdarcyfluxvariables.hh   | 49 ++++++++++++-------
 .../box/boxforchheimerfluxvariables.hh        | 25 +++++++---
 dumux/implicit/box/boxfvelementgeometry.hh    |  5 +-
 .../cellcentered/ccelementvolumevariables.hh  |  2 +-
 .../cellcentered/ccfvelementgeometry.hh       |  4 +-
 .../implicit/cellcentered/cclocaljacobian.hh  |  2 +-
 .../spatialparams/boxspatialparams1p.hh       |  7 +++
 test/implicit/2p/lensspatialparams.hh         |  8 +--
 9 files changed, 79 insertions(+), 47 deletions(-)

diff --git a/dumux/implicit/1p2c/1p2cfluxvariables.hh b/dumux/implicit/1p2c/1p2cfluxvariables.hh
index 81a09a6098..9dd73c1563 100644
--- a/dumux/implicit/1p2c/1p2cfluxvariables.hh
+++ b/dumux/implicit/1p2c/1p2cfluxvariables.hh
@@ -351,14 +351,22 @@ protected:
                      const ElementVolumeVariables &elemVolVars)
     {
         const SpatialParams &sp = problem.spatialParams();
-        sp.meanK(K_,
-                 sp.intrinsicPermeability(element,
-                                          fvGeometry_,
-                                          face().i),
-                 sp.intrinsicPermeability(element,
-                                          fvGeometry_,
-                                          face().j));
-
+        if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
+        {
+            sp.meanK(K_,
+                     sp.intrinsicPermeability(element,
+                                              fvGeometry_,
+                                              face().i),
+                     sp.intrinsicPermeability(element,
+                                              fvGeometry_,
+                                              face().j));
+        }
+        else
+        {
+            sp.meanK(K_,
+                     sp.elementIntrinsicPermeability(*fvGeometry_.neighbors[face().i]),
+                     sp.elementIntrinsicPermeability(*fvGeometry_.neighbors[face().j]));
+        }
     }
 
     /*!
diff --git a/dumux/implicit/box/boxdarcyfluxvariables.hh b/dumux/implicit/box/boxdarcyfluxvariables.hh
index 6e5d3df449..e950b5323e 100644
--- a/dumux/implicit/box/boxdarcyfluxvariables.hh
+++ b/dumux/implicit/box/boxdarcyfluxvariables.hh
@@ -91,7 +91,7 @@ public:
                  const int faceIdx,
                  const ElementVolumeVariables &elemVolVars,
                  const bool onBoundary = false)
-        : fvGeometry_(fvGeometry), faceIdx_(faceIdx), onBoundary_(onBoundary)
+    : fvGeometry_(fvGeometry), faceIdx_(faceIdx), onBoundary_(onBoundary)
     {
         mobilityUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MobilityUpwindWeight);
         calculateGradients_(problem, element, elemVolVars);
@@ -165,17 +165,7 @@ public:
     }
 
 protected:
-    const FVElementGeometry &fvGeometry_;   //!< Information about the geometry of discretization
-    const unsigned int faceIdx_;            //!< The index of the sub control volume face
-    const bool      onBoundary_;                //!< Specifying whether we are currently on the boundary of the simulation domain
-    unsigned int    upstreamIdx_[numPhases] , downstreamIdx_[numPhases]; //!< local index of the upstream / downstream vertex
-    Scalar          volumeFlux_[numPhases] ;    //!< Velocity multiplied with normal (magnitude=area)
-    DimVector       velocity_[numPhases] ;      //!< The velocity as determined by Darcy's law or by the Forchheimer relation
-    Scalar          kGradPNormal_[numPhases] ;  //!< Permeability multiplied with gradient in potential, multiplied with normal (magnitude=area)
-    DimVector       kGradP_[numPhases] ; //!< Permeability multiplied with gradient in potential
-    DimVector       gradPotential_[numPhases] ; //!< Gradient of potential, which drives flow
-    Scalar          mobilityUpwindWeight_;      //!< Upwind weight for mobility. Set to one for full upstream weighting
-
+    
     /*
      * \brief Calculation of the potential gradients
      *
@@ -255,13 +245,23 @@ protected:
         // calculate the mean intrinsic permeability
         const SpatialParams &spatialParams = problem.spatialParams();
         Tensor K;
-        spatialParams.meanK(K,
-                            spatialParams.intrinsicPermeability(element,
-                                                                fvGeometry_,
-                                                                face().i),
-                            spatialParams.intrinsicPermeability(element,
-                                                                fvGeometry_,
-                                                                face().j));
+        if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
+        {
+            spatialParams.meanK(K,
+                                spatialParams.intrinsicPermeability(element,
+                                                                    fvGeometry_,
+                                                                    face().i),
+                                spatialParams.intrinsicPermeability(element,
+                                                                    fvGeometry_,
+                                                                    face().j));
+        }
+        else
+        {
+            spatialParams.meanK(K,
+                                spatialParams.elementIntrinsicPermeability(*fvGeometry_.neighbors[face().i]),
+                                spatialParams.elementIntrinsicPermeability(*fvGeometry_.neighbors[face().j]));
+        }
+        
         // loop over all phases
         for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
         {
@@ -306,6 +306,17 @@ protected:
             volumeFlux_[phaseIdx] = velocity_[phaseIdx] * face().normal ;
         }// loop all phases
     }
+
+    const FVElementGeometry &fvGeometry_;   //!< Information about the geometry of discretization
+    const unsigned int faceIdx_;            //!< The index of the sub control volume face
+    const bool      onBoundary_;                //!< Specifying whether we are currently on the boundary of the simulation domain
+    unsigned int    upstreamIdx_[numPhases] , downstreamIdx_[numPhases]; //!< local index of the upstream / downstream vertex
+    Scalar          volumeFlux_[numPhases] ;    //!< Velocity multiplied with normal (magnitude=area)
+    DimVector       velocity_[numPhases] ;      //!< The velocity as determined by Darcy's law or by the Forchheimer relation
+    Scalar          kGradPNormal_[numPhases] ;  //!< Permeability multiplied with gradient in potential, multiplied with normal (magnitude=area)
+    DimVector       kGradP_[numPhases] ; //!< Permeability multiplied with gradient in potential
+    DimVector       gradPotential_[numPhases] ; //!< Gradient of potential, which drives flow
+    Scalar          mobilityUpwindWeight_;      //!< Upwind weight for mobility. Set to one for full upstream weighting
 };
 
 } // end namespace
diff --git a/dumux/implicit/box/boxforchheimerfluxvariables.hh b/dumux/implicit/box/boxforchheimerfluxvariables.hh
index ec3117cb08..78d008f203 100644
--- a/dumux/implicit/box/boxforchheimerfluxvariables.hh
+++ b/dumux/implicit/box/boxforchheimerfluxvariables.hh
@@ -133,14 +133,23 @@ protected:
         // calculate the mean intrinsic permeability
         const SpatialParams &spatialParams = problem.spatialParams();
         Tensor K;
-        spatialParams.meanK(K,
-                            spatialParams.intrinsicPermeability(element,
-                                                                this->fvGeometry_,
-                                                                this->face().i),
-                            spatialParams.intrinsicPermeability(element,
-                                                                this->fvGeometry_,
-                                                                this->face().j));
-
+        if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
+        {
+            spatialParams.meanK(K,
+                                spatialParams.intrinsicPermeability(element,
+                                                                    fvGeometry_,
+                                                                    face().i),
+                                spatialParams.intrinsicPermeability(element,
+                                                                    fvGeometry_,
+                                                                    face().j));
+        }
+        else
+        {
+            spatialParams.meanK(K,
+                                spatialParams.elementIntrinsicPermeability(*fvGeometry_.neighbors[face().i]),
+                                spatialParams.elementIntrinsicPermeability(*fvGeometry_.neighbors[face().j]));
+        }
+        
         // obtain the Forchheimer coefficient from the spatial parameters
         const Scalar forchCoeff = spatialParams.forchCoeff(element,
                                                           this->fvGeometry_,
diff --git a/dumux/implicit/box/boxfvelementgeometry.hh b/dumux/implicit/box/boxfvelementgeometry.hh
index db67464594..2a8dcf4c97 100644
--- a/dumux/implicit/box/boxfvelementgeometry.hh
+++ b/dumux/implicit/box/boxfvelementgeometry.hh
@@ -591,7 +591,10 @@ public:
     std::vector<ElementPointer> neighbors; //!< needed for compatibility with cc models
     
     const LocalFiniteElementCache feCache_;
-
+    
+    void updateInner(const Element& element) //!< needed for compatibility with cc models
+    {}
+    
     void update(const GridView& gridView, const Element& element)
     {
         const Geometry& geometry = element.geometry();
diff --git a/dumux/implicit/cellcentered/ccelementvolumevariables.hh b/dumux/implicit/cellcentered/ccelementvolumevariables.hh
index 6700444b86..8a34634c1f 100644
--- a/dumux/implicit/cellcentered/ccelementvolumevariables.hh
+++ b/dumux/implicit/cellcentered/ccelementvolumevariables.hh
@@ -90,7 +90,7 @@ public:
                     = globalSol[problem.elementMapper().map(neighbor)];
 
             FVElementGeometry neighborFVGeom;
-            neighborFVGeom.updateInner(problem.gridView(), neighbor);
+            neighborFVGeom.updateInner(neighbor);
             
             (*this)[i].update(solI,
                               problem,
diff --git a/dumux/implicit/cellcentered/ccfvelementgeometry.hh b/dumux/implicit/cellcentered/ccfvelementgeometry.hh
index 9c9020db23..79e88ddee0 100644
--- a/dumux/implicit/cellcentered/ccfvelementgeometry.hh
+++ b/dumux/implicit/cellcentered/ccfvelementgeometry.hh
@@ -114,7 +114,7 @@ public:
     int numFAP; //!< number of flux approximation points
     std::vector<ElementPointer> neighbors; //!< stores pointers for the neighboring elements
     
-    void updateInner(const GridView& gridView, const Element& e)
+    void updateInner(const Element& e)
     {
         const Geometry& geometry = e.geometry();
 
@@ -143,7 +143,7 @@ public:
     
     void update(const GridView& gridView, const Element& e)
     {
-        updateInner(gridView, e);
+        updateInner(e);
         
         const Geometry& geometry = e.geometry();
 
diff --git a/dumux/implicit/cellcentered/cclocaljacobian.hh b/dumux/implicit/cellcentered/cclocaljacobian.hh
index 2da777b72a..1aa4d7ff50 100644
--- a/dumux/implicit/cellcentered/cclocaljacobian.hh
+++ b/dumux/implicit/cellcentered/cclocaljacobian.hh
@@ -390,7 +390,7 @@ protected:
     {
         const Element& neighbor = *(fvElemGeom_.neighbors[neighborIdx]);
 	FVElementGeometry neighborFVGeom;
-	neighborFVGeom.updateInner(problemPtr_->gridView(), neighbor);
+	neighborFVGeom.updateInner(neighbor);
 	
         int globalIdx = problemPtr_->elementMapper().map(neighbor);
 
diff --git a/dumux/material/spatialparams/boxspatialparams1p.hh b/dumux/material/spatialparams/boxspatialparams1p.hh
index bb73bcfd7b..b1c6e0e706 100644
--- a/dumux/material/spatialparams/boxspatialparams1p.hh
+++ b/dumux/material/spatialparams/boxspatialparams1p.hh
@@ -107,6 +107,13 @@ public:
                 result[i][j] = harmonicMean(K1[i][j], K2[i][j]);
     }
 
+    Scalar elementIntrinsicPermeability (const Element &element) const
+    {
+        FVElementGeometry fvGeometry;
+        fvGeometry.subContVol[0].global = element.geometry().center();
+        return asImp_().intrinsicPermeability(element, fvGeometry, /*scvIdx=*/0);
+    }
+    
     /*!
      * \brief Function for defining the intrinsic (absolute) permeability.
      *
diff --git a/test/implicit/2p/lensspatialparams.hh b/test/implicit/2p/lensspatialparams.hh
index 50d24f9783..63b909c947 100644
--- a/test/implicit/2p/lensspatialparams.hh
+++ b/test/implicit/2p/lensspatialparams.hh
@@ -146,16 +146,13 @@ public:
                                  const FVElementGeometry &fvElemGeom,
                                  int scvIdx) const
     {
-        if (!isBox_ && scvIdx > 0)
-            return intrinsicPermeability(*(fvElemGeom.neighbors[scvIdx]), fvElemGeom, 0);
-        
         const GlobalPosition& globalPos = fvElemGeom.subContVol[scvIdx].global;
         
         if (isInLens_(globalPos))
             return lensK_;
         return outerK_;
     }
-
+    
     /*!
      * \brief Porosity
      *
@@ -181,9 +178,6 @@ public:
                                                 const FVElementGeometry &fvElemGeom,
                                                 int scvIdx) const
     {
-        if (!isBox_ && scvIdx > 0)
-            return materialLawParams(*(fvElemGeom.neighbors[scvIdx]), fvElemGeom, 0);
-        
         const GlobalPosition& globalPos = fvElemGeom.subContVol[scvIdx].global;
         
         if (isInLens_(globalPos))
-- 
GitLab