From fa63c17b6f22048d319bfdab3cae85b8c84be8ec Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Thu, 9 Dec 2010 18:55:48 +0000
Subject: [PATCH] extend/correct a few comments, remove artifacts of renamed
 concepts

(phaseState -> fluidState, material context -> material parameters, etc)

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4834 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
 .../msm/1p2cvs2p/lensspatialparameters2p.hh   |  2 +-
 .../buckleyleverett_spatialparams.hh          |  2 +-
 .../msm/mcwhorter/mcwhorter_spatialparams.hh  |  4 +-
 dumux/boxmodels/2p/2plocalresidual.hh         | 10 ++-
 dumux/boxmodels/2p2c/2p2cproperties.hh        |  2 +-
 dumux/boxmodels/common/pdelabboxassembler.hh  | 82 +++++++++++--------
 .../boxmodels/richards/richardsproperties.hh  |  2 +-
 dumux/material/fluidsystems/2p_system.hh      | 24 +++---
 test/boxmodels/2p/lensspatialparameters.hh    |  2 +-
 .../2p2c/injectionspatialparameters.hh        |  2 +-
 .../2p2cni/waterairspatialparameters.hh       |  2 +-
 .../1p/test_diffusion_spatialparams.hh        |  2 +-
 test/decoupled/2p/test_impes_spatialparams.hh |  2 +-
 .../2p/test_transport_spatialparams.hh        |  2 +-
 .../2p2c/test_dec2p2c_spatialparams.hh        |  2 +-
 tutorial/tutorialspatialparameters_coupled.hh | 11 +--
 .../tutorialspatialparameters_decoupled.hh    |  9 +-
 17 files changed, 92 insertions(+), 70 deletions(-)

diff --git a/appl/lecture/msm/1p2cvs2p/lensspatialparameters2p.hh b/appl/lecture/msm/1p2cvs2p/lensspatialparameters2p.hh
index 50bb1459f8..c0cab07b07 100644
--- a/appl/lecture/msm/1p2cvs2p/lensspatialparameters2p.hh
+++ b/appl/lecture/msm/1p2cvs2p/lensspatialparameters2p.hh
@@ -144,7 +144,7 @@ public:
         return outerPorosity_;
     }
 
-    // return the brooks-corey context depending on the position
+    // return the parameter object for the Brooks-Corey material law which depends on the position
     const MaterialLawParams& materialLawParams(const Element &element,
                                                 const FVElementGeometry &fvElemGeom,
                                                 int scvIdx) const
diff --git a/appl/lecture/msm/buckleyleverett/buckleyleverett_spatialparams.hh b/appl/lecture/msm/buckleyleverett/buckleyleverett_spatialparams.hh
index 150d22ece9..7a0d85f01c 100644
--- a/appl/lecture/msm/buckleyleverett/buckleyleverett_spatialparams.hh
+++ b/appl/lecture/msm/buckleyleverett/buckleyleverett_spatialparams.hh
@@ -65,7 +65,7 @@ public:
     }
 
 
-    // return the brooks-corey context depending on the position
+    // return the parameter object for the Brooks-Corey material law which depends on the position
     const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
     {
             return materialLawParams_;
diff --git a/appl/lecture/msm/mcwhorter/mcwhorter_spatialparams.hh b/appl/lecture/msm/mcwhorter/mcwhorter_spatialparams.hh
index 868caa590d..1e6be79da4 100644
--- a/appl/lecture/msm/mcwhorter/mcwhorter_spatialparams.hh
+++ b/appl/lecture/msm/mcwhorter/mcwhorter_spatialparams.hh
@@ -65,10 +65,10 @@ public:
     }
 
 
-    // return the brooks-corey context depending on the position
+    // return the brooks-corey material parameter object which depends on the position
     const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
     {
-            return materialLawParams_;
+        return materialLawParams_;
     }
 
 
diff --git a/dumux/boxmodels/2p/2plocalresidual.hh b/dumux/boxmodels/2p/2plocalresidual.hh
index ac135ff374..953a657512 100644
--- a/dumux/boxmodels/2p/2plocalresidual.hh
+++ b/dumux/boxmodels/2p/2plocalresidual.hh
@@ -152,9 +152,13 @@ public:
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
         {
             // calculate the flux in the normal direction of the
-            // current sub control volume face
-            fluxVars.intrinsicPermeability().mv(fluxVars.potentialGrad(phaseIdx),
-                                                tmpVec);
+            // current sub control volume face:
+            //
+            // v = - (K grad p) * n
+            //
+            // (the minus comes from the Darcy law which states that
+            // the flux is from high to low pressure potentials.)
+            fluxVars.intrinsicPermeability().mv(fluxVars.potentialGrad(phaseIdx), tmpVec);
             Scalar normalFlux = - (tmpVec*fluxVars.face().normal);
 
             // data attached to upstream and the downstream vertices
diff --git a/dumux/boxmodels/2p2c/2p2cproperties.hh b/dumux/boxmodels/2p2c/2p2cproperties.hh
index dfd1c9e183..9ebcfa6fed 100644
--- a/dumux/boxmodels/2p2c/2p2cproperties.hh
+++ b/dumux/boxmodels/2p2c/2p2cproperties.hh
@@ -53,7 +53,7 @@ NEW_PROP_TAG(SpatialParameters); //!< The type of the spatial parameters
 NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations
 
 NEW_PROP_TAG(MaterialLaw);   //!< The material law which ought to be used (extracted from the spatial parameters)
-NEW_PROP_TAG(MaterialLawParams); //!< The context material law (extracted from the spatial parameters)
+NEW_PROP_TAG(MaterialLawParams); //!< The parameters of the material law (extracted from the spatial parameters)
 
 NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem
 NEW_PROP_TAG(MobilityUpwindAlpha); //!< The value of the upwind parameter for the mobility
diff --git a/dumux/boxmodels/common/pdelabboxassembler.hh b/dumux/boxmodels/common/pdelabboxassembler.hh
index 2f3cce1dbf..fbac1a1985 100644
--- a/dumux/boxmodels/common/pdelabboxassembler.hh
+++ b/dumux/boxmodels/common/pdelabboxassembler.hh
@@ -337,9 +337,9 @@ public:
      *
      * \param relTol The relative error below which a vertex won't be
      *               reassembled. Note that this specifies the
-     *               relative error between the current evaluation
-     *               point and the current solution and _not_ the
-     *               delta vector of the Newton iteration!
+     *               worst-case relative error between the last
+     *               linearization point and the current solution and
+     *               _not_ the delta vector of the Newton iteration!
      */
     void computeColors(Scalar relTol)
     { 
@@ -349,10 +349,9 @@ public:
         ElementIterator elemIt = gridView_().template begin<0>();
         ElementIterator elemEndIt = gridView_().template end<0>();
 
-        greenElems_ = 0;
+        // mark the red vertices and update the tolerance of the
+        // linearization which actually will get achieved
         nextReassembleTolerance_ = 0;
-
-        // mark the red vertices
         for (int i = 0; i < vertexColor_.size(); ++i) {
             vertexColor_[i] = Green;
             if (vertexDelta_[i] > relTol) {
@@ -366,26 +365,33 @@ public:
 
         // Mark all red elements
         for (; elemIt != elemEndIt; ++elemIt) {
-            bool needReassemble = false;
+            // find out whether the current element features a red
+            // vertex
+            bool isRed = false;
             int numVerts = elemIt->template count<dim>();
             for (int i=0; i < numVerts; ++i) {
                 int globalI = vertexMapper_().map(*elemIt, i, dim);
                 if (vertexColor_[globalI] == Red) {
-                    needReassemble = true;
+                    isRed = true;
                     break;
                 }
             };
             
+            // if yes, the element color is also red, else it is not
+            // red, i.e. green for the mean time
             int globalElemIdx = elementMapper_().map(*elemIt);
-            elementColor_[globalElemIdx] = needReassemble?Red:Green;
+            if (isRed)
+                elementColor_[globalElemIdx] = Red;
+            else
+                elementColor_[globalElemIdx] = Green;
         }
         
-        // Mark all yellow vertices
+        // Mark yellow vertices (as orange for the mean time)
         elemIt = gridView_().template begin<0>();
         for (; elemIt != elemEndIt; ++elemIt) {
             int elemIdx = this->elementMapper_().map(*elemIt);
-            if (elementColor_[elemIdx] == Green)
-                continue; // green elements do not tint vertices
+            if (elementColor_[elemIdx] != Red)
+                continue; // non-red elements do not tint vertices
                           // yellow!
             
             int numVerts = elemIt->template count<dim>();
@@ -398,8 +404,7 @@ public:
             };
         }
 
-        // Mark all yellow elements
-        int numGreen = 0;
+        // Mark yellow elements
         elemIt = gridView_().template begin<0>();
         for (; elemIt != elemEndIt; ++elemIt) {
             int elemIdx = this->elementMapper_().map(*elemIt);
@@ -407,6 +412,8 @@ public:
                 continue; // element is red already!
             }
 
+            // check whether the element features a yellow
+            // (resp. orange at this point) vertex
             bool isYellow = false;
             int numVerts = elemIt->template count<dim>();
             for (int i=0; i < numVerts; ++i) {
@@ -417,15 +424,14 @@ public:
                 }
             };
             
-            if (isYellow) {
+            if (isYellow)
                 elementColor_[elemIdx] = Yellow;
-            }
-            else // elementColor_[elemIdx] == Green
-                ++ numGreen;
         }
 
         // Demote orange vertices to yellow ones if it has at least
-        // one green element as a neighbor
+        // one green element as a neighbor. also count the total
+        // number of green elements for statistical purposes.
+        greenElems_ = 0;
         elemIt = gridView_().template begin<0>();
         for (; elemIt != elemEndIt; ++elemIt) {
             int elemIdx = this->elementMapper_().map(*elemIt);
@@ -444,17 +450,19 @@ public:
             };
         }
 
-        // promote the remaining orange vertices to red ones
+        // promote the remaining orange vertices to red
         for (int i=0; i < vertexColor_.size(); ++i) {
-            // if a vertex is green or yellow don't recolor it!
+            // if a vertex is green or yellow don't do anything!
             if (vertexColor_[i] == Green || vertexColor_[i] == Yellow)
                 continue;
+
+            // make sure the vertex is red (this is a no-op vertices
+            // which are already red!)
+            vertexColor_[i] = Red;
             
-            // set discrepancy for this vertex to 0 because the
-            // system will be consistently linearized for this
-            // vertex
+            // set the error of this vertex to 0 because the system
+            // will be consistently linearized at this vertex
             vertexDelta_[i] = 0.0;
-            vertexColor_[i] = Red;
         };
 
         greenElems_ = gridView_().comm().sum(greenElems_);
@@ -555,9 +563,10 @@ private:
         // allocate raw matrix
         matrix_ = new Matrix(nVerts, nVerts, Matrix::random);
         
-        // find out how many neighbors each vertex has
+        // find out the global indices of the neighboring vertices of
+        // each vertex
         typedef std::set<int> NeighborSet;
-        std::vector<NeighborSet > neighbors(nVerts);
+        std::vector<NeighborSet> neighbors(nVerts);
         ElementIterator eIt = gridView_().template begin<0>();
         const ElementIterator eEndIt = gridView_().template end<0>();
         for (; eIt != eEndIt; ++eIt) {
@@ -565,25 +574,32 @@ private:
 
             // loop over all element vertices
             int n = elem.template count<dim>();
-            for (int i = 0; i < n; ++i) {
+            for (int i = 0; i < n - 1; ++i) {
                 int globalI = vertexMapper_().map(*eIt, i, dim);
-                for (int j = 0; j < n; ++j) {
+                for (int j = i + 1; j < n; ++j) {
                     int globalJ = vertexMapper_().map(*eIt, j, dim);
-                    // insert into BCRS matrix
+                    // make sure that vertex j is in the neighbor set
+                    // of vertex i and vice-versa
                     neighbors[globalI].insert(globalJ);
+                    neighbors[globalJ].insert(globalI);
                 }
             }
         };
         
-        // allocate space for the rows
+        // make vertices neighbors to themselfs
+        for (int i = 0; i < nVerts; ++i)
+            neighbors[i].insert(i);
+        
+        // allocate space for the rows of the matrix
         for (int i = 0; i < nVerts; ++i) {
             matrix_->setrowsize(i, neighbors[i].size());
         }
         matrix_->endrowsizes();
 
-        // allocate space for the rows
+        // fill the rows with indices. each vertex talks to all of its
+        // neighbors. (it also talks to itself since vertices are
+        // sometimes quite egocentric.)
         for (int i = 0; i < nVerts; ++i) {
-            // off-diagonal entries
             typename NeighborSet::iterator nIt = neighbors[i].begin();
             typename NeighborSet::iterator nEndIt = neighbors[i].end();
             for (; nIt != nEndIt; ++nIt) {
diff --git a/dumux/boxmodels/richards/richardsproperties.hh b/dumux/boxmodels/richards/richardsproperties.hh
index 52f64170aa..b76b4a48a3 100644
--- a/dumux/boxmodels/richards/richardsproperties.hh
+++ b/dumux/boxmodels/richards/richardsproperties.hh
@@ -50,7 +50,7 @@ NEW_PROP_TAG(NumPhases);   //!< Number of fluid phases in the system
 NEW_PROP_TAG(RichardsIndices); //!< Enumerations used by the Richards models
 NEW_PROP_TAG(SpatialParameters); //!< The type of the spatial parameters object
 NEW_PROP_TAG(MaterialLaw);   //!< The material law which ought to be used (by default extracted from the spatial parameters)
-NEW_PROP_TAG(MaterialLawParams); //!< The context material law (by default extracted from the spatial parameters)
+NEW_PROP_TAG(MaterialLawParams); //!< The type of the parameter object for the material law (by default extracted from the spatial parameters)
 NEW_PROP_TAG(FluidSystem); //!< The fluid system to be used for the Richards model
 NEW_PROP_TAG(FluidState); //!< The fluid state to be used for the Richards model
 NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem
diff --git a/dumux/material/fluidsystems/2p_system.hh b/dumux/material/fluidsystems/2p_system.hh
index ad65252308..583f121575 100644
--- a/dumux/material/fluidsystems/2p_system.hh
+++ b/dumux/material/fluidsystems/2p_system.hh
@@ -135,7 +135,7 @@ public:
      * \param phaseIdx index of the phase
      * \param temperature phase temperature in [K]
      * \param pressure phase pressure in [Pa]
-     * \param phaseState The fluid state of the two-phase model
+     * \param fluidState The fluid state of the two-phase model
      * \tparam FluidState the fluid state class of the two-phase model
      * \return returns the density of the phase
      */
@@ -143,7 +143,7 @@ public:
     static Scalar phaseDensity(int phaseIdx,
                                Scalar temperature,
                                Scalar pressure,
-                               const FluidState &phaseState)
+                               const FluidState &fluidState)
     {
         switch (phaseIdx) {
         case wPhaseIdx:
@@ -160,7 +160,7 @@ public:
      * \param phaseIdx index of the phase
      * \param temperature phase temperature in [K]
      * \param pressure phase pressure in [Pa]
-     * \param phaseState The fluid state of the two-phase model
+     * \param fluidState The fluid state of the two-phase model
      * \tparam FluidState the fluid state class of the two-phase model
      * \return returns the viscosity of the phase [Pa*s]
      */
@@ -168,7 +168,7 @@ public:
     static Scalar phaseViscosity(int phaseIdx,
                                  Scalar temperature,
                                  Scalar pressure,
-                                 const FluidState &phaseState)
+                                 const FluidState &fluidState)
     {
         switch (phaseIdx) {
         case wPhaseIdx:
@@ -185,15 +185,15 @@ public:
      * \param phaseIdx index of the phase
      * \param temperature phase temperature in [K]
      * \param pressure phase pressure in [Pa]
-     * \param phaseState The fluid state of the two-phase model
+     * \param fluidState The fluid state of the two-phase model
      * \tparam FluidState the fluid state class of the two-phase model
      * \return returns the specific enthalpy of the phase [J/kg]
      */
     template <class FluidState>
     static Scalar phaseEnthalpy(int phaseIdx,
-                           Scalar temperature,
-                           Scalar pressure,
-                           const FluidState &phaseState)
+                                Scalar temperature,
+                                Scalar pressure,
+                                const FluidState &fluidState)
     {
         switch (phaseIdx) {
         case wPhaseIdx:
@@ -210,15 +210,15 @@ public:
      * \param phaseIdx index of the phase
      * \param temperature phase temperature in [K]
      * \param pressure phase pressure in [Pa]
-     * \param phaseState The fluid state of the two-phase model
+     * \param fluidState The fluid state of the two-phase model
      * \tparam FluidState the fluid state class of the two-phase model
      * \return returns the specific internal energy of the phase [J/kg]
      */
     template <class FluidState>
     static Scalar phaseInternalEnergy(int phaseIdx,
-                                 Scalar temperature,
-                                 Scalar pressure,
-                                 const FluidState &phaseState)
+                                      Scalar temperature,
+                                      Scalar pressure,
+                                      const FluidState &fluidState)
     {
         switch (phaseIdx) {
         case wPhaseIdx:
diff --git a/test/boxmodels/2p/lensspatialparameters.hh b/test/boxmodels/2p/lensspatialparameters.hh
index e74ad8493b..7812b97ea4 100644
--- a/test/boxmodels/2p/lensspatialparameters.hh
+++ b/test/boxmodels/2p/lensspatialparameters.hh
@@ -125,7 +125,7 @@ public:
                     int scvIdx) const
     { return 0.4; }
 
-    // return the brooks-corey context depending on the position
+    // return the parameter object for the Brooks-Corey material law which depends on the position
     const MaterialLawParams& materialLawParams(const Element &element,
                                                 const FVElementGeometry &fvElemGeom,
                                                 int scvIdx) const
diff --git a/test/boxmodels/2p2c/injectionspatialparameters.hh b/test/boxmodels/2p2c/injectionspatialparameters.hh
index 392df67e5a..38ecf667f5 100644
--- a/test/boxmodels/2p2c/injectionspatialparameters.hh
+++ b/test/boxmodels/2p2c/injectionspatialparameters.hh
@@ -160,7 +160,7 @@ public:
 
 
     /*!
-     * \brief return the brooks-corey context depending on the position
+     * \brief return the parameter object for the Brooks-Corey material law which depends on the position
      *
     * \param element The current finite element
     * \param fvElemGeom The current finite volume geometry of the element
diff --git a/test/boxmodels/2p2cni/waterairspatialparameters.hh b/test/boxmodels/2p2cni/waterairspatialparameters.hh
index 2e9dcde0a4..a3d4170313 100644
--- a/test/boxmodels/2p2cni/waterairspatialparameters.hh
+++ b/test/boxmodels/2p2cni/waterairspatialparameters.hh
@@ -160,7 +160,7 @@ public:
 
 
     /*!
-     * \brief return the brooks-corey context depending on the position
+     * \brief return the parameter object for the Brooks-Corey material law which depends on the position
      *
     * \param element The current finite element
     * \param fvElemGeom The current finite volume geometry of the element
diff --git a/test/decoupled/1p/test_diffusion_spatialparams.hh b/test/decoupled/1p/test_diffusion_spatialparams.hh
index fd39217441..b18507c837 100644
--- a/test/decoupled/1p/test_diffusion_spatialparams.hh
+++ b/test/decoupled/1p/test_diffusion_spatialparams.hh
@@ -73,7 +73,7 @@ public:
     }
 
 
-    // return the brooks-corey context depending on the position
+    // return the parameter object for the Brooks-Corey material law which depends on the position
     const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
     {
             return materialLawParams_;
diff --git a/test/decoupled/2p/test_impes_spatialparams.hh b/test/decoupled/2p/test_impes_spatialparams.hh
index cab3c4125b..9049f81a58 100644
--- a/test/decoupled/2p/test_impes_spatialparams.hh
+++ b/test/decoupled/2p/test_impes_spatialparams.hh
@@ -70,7 +70,7 @@ public:
     }
 
 
-    // return the brooks-corey context depending on the position
+    // return the parameter object for the Brooks-Corey material law which depends on the position
     const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
     {
             return materialLawParams_;
diff --git a/test/decoupled/2p/test_transport_spatialparams.hh b/test/decoupled/2p/test_transport_spatialparams.hh
index 14e3fd6d0e..8c17642e68 100644
--- a/test/decoupled/2p/test_transport_spatialparams.hh
+++ b/test/decoupled/2p/test_transport_spatialparams.hh
@@ -68,7 +68,7 @@ public:
     }
 
 
-    // return the brooks-corey context depending on the position
+    // return the parameter object for the Brooks-Corey material law which depends on the position
     const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
     {
             return materialLawParams_;
diff --git a/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh b/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh
index 4c78fccda3..088e31c295 100644
--- a/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh
+++ b/test/decoupled/2p2c/test_dec2p2c_spatialparams.hh
@@ -70,7 +70,7 @@ public:
     }
 
 
-    // return the brooks-corey context depending on the position
+    // return the parameter object for the Brooks-Corey material law which depends on the position
     const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
     {
             return materialLawParams_;
diff --git a/tutorial/tutorialspatialparameters_coupled.hh b/tutorial/tutorialspatialparameters_coupled.hh
index b6f62f29aa..f595b6c5bd 100644
--- a/tutorial/tutorialspatialparameters_coupled.hh
+++ b/tutorial/tutorialspatialparameters_coupled.hh
@@ -56,11 +56,11 @@ class TutorialSpatialParametersCoupled: public BoxSpatialParameters<TypeTag> /*@
     typedef typename Grid::Traits::template Codim<0>::Entity Element;
 
     // select material law to be used
-    typedef RegularizedBrooksCorey<Scalar>        RawMaterialLaw;     /*@\label{tutorial-coupled:rawlaw}@*/
+    typedef RegularizedBrooksCorey<Scalar> EffectiveMaterialLaw;     /*@\label{tutorial-coupled:rawlaw}@*/
 
 public:
     // adapter for absolute law
-    typedef EffToAbsLaw<RawMaterialLaw> MaterialLaw;        /*@\label{tutorial-coupled:eff2abs}@*/
+    typedef EffToAbsLaw<EffectiveMaterialLaw> MaterialLaw;        /*@\label{tutorial-coupled:eff2abs}@*/
     // determine appropriate parameters depening on selected materialLaw
     typedef typename MaterialLaw::Params MaterialLawParams;    /*@\label{tutorial-coupled:matLawObjectType}@*/
 
@@ -83,10 +83,11 @@ public:
         return 0.2;
     }
 
-    // return the materialLaw context (i.e. BC, regularizedVG, etc) depending on the position
+    // return the parameter object for the material law (i.e. Brooks-Corey)
+    // which may vary with the spatial position
     const MaterialLawParams& materialLawParams(const Element &element,            /*@\label{tutorial-coupled:matLawParams}@*/
-                                            const FVElementGeometry &fvElemGeom,
-                                            int scvIdx) const
+                                               const FVElementGeometry &fvElemGeom,
+                                               int scvIdx) const
     {
         return materialParams_;
     }
diff --git a/tutorial/tutorialspatialparameters_decoupled.hh b/tutorial/tutorialspatialparameters_decoupled.hh
index 22479f852b..bc021ae22e 100644
--- a/tutorial/tutorialspatialparameters_decoupled.hh
+++ b/tutorial/tutorialspatialparameters_decoupled.hh
@@ -47,10 +47,10 @@ class TutorialSpatialParametersDecoupled
     typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix;
 
     // material law typedefs
-    typedef RegularizedBrooksCorey<Scalar>                RawMaterialLaw;
-//    typedef LinearMaterial<Scalar>                        RawMaterialLaw;
+    typedef RegularizedBrooksCorey<Scalar> EffectiveMaterialLaw;
+//    typedef LinearMaterial<Scalar> EffectiveMaterialLaw;
 public:
-    typedef EffToAbsLaw<RawMaterialLaw>               MaterialLaw;
+    typedef EffToAbsLaw<EffectiveMaterialLaw> MaterialLaw;
     typedef typename MaterialLaw::Params MaterialLawParams;
 
     //! Update the spatial parameters with the flow solution after a timestep.
@@ -73,7 +73,8 @@ public:
         return 0.2;
     }
 
-    //! return the material law context (i.e. BC, regularizedVG, etc) depending on the position
+    //! return the parameter object for the material law (i.e. Brooks-Corey)
+    //! which may vary with the spatial position
     const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
     {
             return materialLawParams_;
-- 
GitLab