]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Finish documenting
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 28 Oct 2006 19:26:18 +0000 (19:26 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 28 Oct 2006 19:26:18 +0000 (19:26 +0000)
git-svn-id: https://svn.dealii.org/trunk@14119 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-21/step-21.cc

index 7fd03abb903ee499252a28ec287157d1c2d0318a..981b0694685831cc47a38302762d6eda7a38e16f 100644 (file)
 /*    to the file deal.II/doc/license.html for the  text  and     */
 /*    further information on this license.                        */
 
-                                // This program is an adaptation of step-20
-                                // and includes some technique of DG method
-                                // from step-12. A good part of the program
-                                // is therefore very similar to step-20 and
-                                // we will not comment again on these
-                                // parts. Only the new stuff will be
-                                // discussed in more detail.
-
-                                // @sect3{Include files}
-
-                                // All of these include files have been used
-                                // before:
+                                 // This program is an adaptation of step-20
+                                 // and includes some technique of DG method
+                                 // from step-12. A good part of the program
+                                 // is therefore very similar to step-20 and
+                                 // we will not comment again on these
+                                 // parts. Only the new stuff will be
+                                 // discussed in more detail.
+
+                                 // @sect3{Include files}
+
+                                 // All of these include files have been used
+                                 // before:
 #include <base/quadrature_lib.h>
 #include <base/logstream.h>
 #include <base/function.h>
 #include <fstream>
 #include <sstream>
 
-                                // In this program, we use a tensor-valued
-                                // coefficient. Since it may have a spatial
-                                // dependence, we consider it a tensor-valued
-                                // function. The following include file
-                                // provides the <code>TensorFunction</code>
-                                // class that offers such functionality:
+                                 // In this program, we use a tensor-valued
+                                 // coefficient. Since it may have a spatial
+                                 // dependence, we consider it a tensor-valued
+                                 // function. The following include file
+                                 // provides the <code>TensorFunction</code>
+                                 // class that offers such functionality:
 #include <base/tensor_function.h>
 
-                                // The last step is as in all
-                                // previous programs:
+                                 // The last step is as in all
+                                 // previous programs:
 using namespace dealii;
 
 
                                  // @sect3{The <code>TwoPhaseFlowProblem</code> class}
                                  
-                                // This is the main class of the program. It
-                                // is close to the one of step-20, but with a
-                                // few additional functions:
-                                //
-                                // <ul>
-                                //   <li><code>assemble_rhs_S</code> assembles the
-                                //   right hand side of the saturation
-                                //   equation. As explained in the
-                                //   introduction, this can't be integrated
-                                //   into <code>assemble_rhs</code> since it depends
-                                //   on the velocity that is computed in the
-                                //   first part of the time step.
-                                //
-                                //   <li><code>get_maximal_velocity</code> does as its
-                                //   name suggests. This function is used in
-                                //   the computation of the time step size.
-                                //
-                                //   <li><code>project_back_saturation</code> resets
-                                //   all saturation degrees of freedom with
-                                //   values less than zero to zero, and all
-                                //   those with saturations greater than one
-                                //   to one.
-                                // </ul>
-                                //
-                                // The rest of the class should be pretty
-                                // much obvious. The <code>viscosity</code> variable
-                                // stores the viscosity $\mu$ that enters
-                                // several of the formulas in the nonlinear
-                                // equations.
+                                 // This is the main class of the program. It
+                                 // is close to the one of step-20, but with a
+                                 // few additional functions:
+                                 //
+                                 // <ul>
+                                 //   <li><code>assemble_rhs_S</code> assembles the
+                                 //   right hand side of the saturation
+                                 //   equation. As explained in the
+                                 //   introduction, this can't be integrated
+                                 //   into <code>assemble_rhs</code> since it depends
+                                 //   on the velocity that is computed in the
+                                 //   first part of the time step.
+                                 //
+                                 //   <li><code>get_maximal_velocity</code> does as its
+                                 //   name suggests. This function is used in
+                                 //   the computation of the time step size.
+                                 //
+                                 //   <li><code>project_back_saturation</code> resets
+                                 //   all saturation degrees of freedom with
+                                 //   values less than zero to zero, and all
+                                 //   those with saturations greater than one
+                                 //   to one.
+                                 // </ul>
+                                 //
+                                 // The rest of the class should be pretty
+                                 // much obvious. The <code>viscosity</code> variable
+                                 // stores the viscosity $\mu$ that enters
+                                 // several of the formulas in the nonlinear
+                                 // equations.
 template <int dim>
 class TwoPhaseFlowProblem 
 {
@@ -138,14 +138,14 @@ class TwoPhaseFlowProblem
 };
 
 
-                                // @sect3{Equation data}
+                                 // @sect3{Equation data}
 
-                                // @sect4{Pressure right hand side}
-                                // At present, the right hand side of the
-                                // pressure equation is simply the zero
-                                // function. However, the rest of the program
-                                // is fully equipped to deal with anything
-                                // else, if this is desired:
+                                 // @sect4{Pressure right hand side}
+                                 // At present, the right hand side of the
+                                 // pressure equation is simply the zero
+                                 // function. However, the rest of the program
+                                 // is fully equipped to deal with anything
+                                 // else, if this is desired:
 template <int dim>
 class PressureRightHandSide : public Function<dim> 
 {
@@ -153,7 +153,7 @@ class PressureRightHandSide : public Function<dim>
     PressureRightHandSide () : Function<dim>(1) {};
     
     virtual double value (const Point<dim>   &p,
-                         const unsigned int  component = 0) const;
+                          const unsigned int  component = 0) const;
 };
 
 
@@ -161,16 +161,16 @@ class PressureRightHandSide : public Function<dim>
 template <int dim>
 double
 PressureRightHandSide<dim>::value (const Point<dim>  &/*p*/,
-                                  const unsigned int /*component*/) const 
+                                   const unsigned int /*component*/) const 
 {
   return 0;
 }
 
 
-                                // @sect4{Pressure boundary values}
-                                // The next are pressure boundary values. As
-                                // mentioned in the introduction, we choose a
-                                // linear pressure field:
+                                 // @sect4{Pressure boundary values}
+                                 // The next are pressure boundary values. As
+                                 // mentioned in the introduction, we choose a
+                                 // linear pressure field:
 template <int dim>
 class PressureBoundaryValues : public Function<dim> 
 {
@@ -178,29 +178,29 @@ class PressureBoundaryValues : public Function<dim>
     PressureBoundaryValues () : Function<dim>(1) {};
     
     virtual double value (const Point<dim>   &p,
-                         const unsigned int  component = 0) const;
+                          const unsigned int  component = 0) const;
 };
 
 
 template <int dim>
 double
 PressureBoundaryValues<dim>::value (const Point<dim>  &p,
-                                   const unsigned int /*component*/) const 
+                                    const unsigned int /*component*/) const 
 {
   return 1-p[0];
 }
 
 
-                                // @sect4{Saturation boundary values}
+                                 // @sect4{Saturation boundary values}
 
-                                // Then we also need boundary values on the
-                                // inflow portions of the boundary. The
-                                // question whether something is an inflow
-                                // part is decided when assembling the right
-                                // hand side, we only have to provide a
-                                // functional description of the boundary
-                                // values. This is as explained in the
-                                // introduction:
+                                 // Then we also need boundary values on the
+                                 // inflow portions of the boundary. The
+                                 // question whether something is an inflow
+                                 // part is decided when assembling the right
+                                 // hand side, we only have to provide a
+                                 // functional description of the boundary
+                                 // values. This is as explained in the
+                                 // introduction:
 template <int dim>
 class SaturationBoundaryValues : public Function<dim> 
 {
@@ -208,7 +208,7 @@ class SaturationBoundaryValues : public Function<dim>
     SaturationBoundaryValues () : Function<dim>(1) {};
     
     virtual double value (const Point<dim>   &p,
-                         const unsigned int  component = 0) const;
+                          const unsigned int  component = 0) const;
 };
 
 
@@ -216,7 +216,7 @@ class SaturationBoundaryValues : public Function<dim>
 template <int dim>
 double
 SaturationBoundaryValues<dim>::value (const Point<dim> &p,
-                                     const unsigned int /*component*/) const 
+                                      const unsigned int /*component*/) const 
 {
   if (p[0] == 0)
     return 1;
@@ -226,25 +226,25 @@ SaturationBoundaryValues<dim>::value (const Point<dim> &p,
 
 
 
-                                // @sect4{Initial data}
-
-                                // Finally, we need initial data. In reality,
-                                // we only need initial data for the
-                                // saturation, but we are lazy, so we will
-                                // later, before the first time step, simply
-                                // interpolate the entire solution for the
-                                // previous time step from a function that
-                                // contains all vector components.
-                                //
-                                // We therefore simply create a function that
-                                // returns zero in all components. We do that
-                                // by simply forward every function to the
-                                // ZeroFunction class. Why not use that right
-                                // away in the places of this program where
-                                // we presently use the <code>InitialValues</code>
-                                // class? Because this way it is simpler to
-                                // later go back and choose a different
-                                // function for initial values.
+                                 // @sect4{Initial data}
+
+                                 // Finally, we need initial data. In reality,
+                                 // we only need initial data for the
+                                 // saturation, but we are lazy, so we will
+                                 // later, before the first time step, simply
+                                 // interpolate the entire solution for the
+                                 // previous time step from a function that
+                                 // contains all vector components.
+                                 //
+                                 // We therefore simply create a function that
+                                 // returns zero in all components. We do that
+                                 // by simply forward every function to the
+                                 // ZeroFunction class. Why not use that right
+                                 // away in the places of this program where
+                                 // we presently use the <code>InitialValues</code>
+                                 // class? Because this way it is simpler to
+                                 // later go back and choose a different
+                                 // function for initial values.
 template <int dim>
 class InitialValues : public Function<dim> 
 {
@@ -252,10 +252,10 @@ class InitialValues : public Function<dim>
     InitialValues () : Function<dim>(dim+2) {};
     
     virtual double value (const Point<dim>   &p,
-                         const unsigned int  component = 0) const;
+                          const unsigned int  component = 0) const;
 
     virtual void vector_value (const Point<dim> &p, 
-                              Vector<double>   &value) const;
+                               Vector<double>   &value) const;
 
 };
 
@@ -263,7 +263,7 @@ class InitialValues : public Function<dim>
 template <int dim>
 double
 InitialValues<dim>::value (const Point<dim>  &p,
-                          const unsigned int component) const 
+                           const unsigned int component) const 
 {
   return ZeroFunction<dim>(dim+2).value (p, component);
 }
@@ -272,7 +272,7 @@ InitialValues<dim>::value (const Point<dim>  &p,
 template <int dim>
 void
 InitialValues<dim>::vector_value (const Point<dim> &p,
-                                 Vector<double>   &values) const 
+                                  Vector<double>   &values) const 
 {
   ZeroFunction<dim>(dim+2).vector_value (p, values);
 }
@@ -280,23 +280,23 @@ InitialValues<dim>::vector_value (const Point<dim> &p,
 
 
 
-                                // @sect3{The inverse permeability tensor}
+                                 // @sect3{The inverse permeability tensor}
 
-                                // As announced in the introduction, we
-                                // implement two different permeability
-                                // tensor fields. Each of them we put into a
-                                // namespace of its own, so that it will be
-                                // easy later to replace use of one by the
-                                // other in the code.
+                                 // As announced in the introduction, we
+                                 // implement two different permeability
+                                 // tensor fields. Each of them we put into a
+                                 // namespace of its own, so that it will be
+                                 // easy later to replace use of one by the
+                                 // other in the code.
 
-                                // @sect4{Single curving crack permeability}
+                                 // @sect4{Single curving crack permeability}
 
-                                // The first function for the permeability
-                                // was the one that models a single curving
-                                // crack. It was already used at the end of
-                                // step-20, and its functional form is given
-                                // in the introduction of the present
-                                // tutorial program:
+                                 // The first function for the permeability
+                                 // was the one that models a single curving
+                                 // crack. It was already used at the end of
+                                 // step-20, and its functional form is given
+                                 // in the introduction of the present
+                                 // tutorial program:
 namespace SingleCurvingCrack
 {
   template <int dim>
@@ -304,88 +304,88 @@ namespace SingleCurvingCrack
   {
     public:
       virtual void value_list (const std::vector<Point<dim> > &points,
-                              std::vector<Tensor<2,dim> >    &values) const;
+                               std::vector<Tensor<2,dim> >    &values) const;
   };
 
 
   template <int dim>
   void
   KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
-                            std::vector<Tensor<2,dim> >    &values) const
+                             std::vector<Tensor<2,dim> >    &values) const
   {
     Assert (points.size() == values.size(),
-           ExcDimensionMismatch (points.size(), values.size()));
+            ExcDimensionMismatch (points.size(), values.size()));
 
     for (unsigned int p=0; p<points.size(); ++p)
       {
-       values[p].clear ();
+        values[p].clear ();
 
-       const double distance_to_flowline
-         = std::fabs(points[p][1]-0.5-0.1*std::sin(10*points[p][0]));
+        const double distance_to_flowline
+          = std::fabs(points[p][1]-0.5-0.1*std::sin(10*points[p][0]));
       
-       const double permeability = std::max(std::exp(-(distance_to_flowline*
-                                                       distance_to_flowline)
-                                                     / (0.1 * 0.1)),
-                                            0.01);
+        const double permeability = std::max(std::exp(-(distance_to_flowline*
+                                                        distance_to_flowline)
+                                                      / (0.1 * 0.1)),
+                                             0.01);
       
-       for (unsigned int d=0; d<dim; ++d)
-         values[p][d][d] = 1./permeability;
+        for (unsigned int d=0; d<dim; ++d)
+          values[p][d][d] = 1./permeability;
       }
   }
 }
 
 
-                                // @sect4{Random medium permeability}
-
-                                // This function does as announced in the
-                                // introduction, i.e. it creates an overlay
-                                // of exponentials at random places. There is
-                                // one thing worth considering for this
-                                // class. The issue centers around the
-                                // problem that the class creates the centers
-                                // of the exponentials using a random
-                                // function. If we therefore created the
-                                // centers each time we create an object of
-                                // the present type, we would get a different
-                                // list of centers each time. That's not what
-                                // we expect from classes of this type: they
-                                // should reliably represent the same
-                                // function.
-                                //
-                                // The solution to this problem is to make
-                                // the list of centers a static member
-                                // variable of this class, i.e. there exists
-                                // exactly one such variable for the entire
-                                // program, rather than for each object of
-                                // this type. That's exactly what we are
-                                // going to do.
-                                //
-                                // The next problem, however, is that we need
-                                // a way to initialize this variable. Since
-                                // this variable is initialized at the
-                                // beginning of the program, we can't use a
-                                // regular member function for that since
-                                // there may not be an object of this type
-                                // around at the time. The C++ standard
-                                // therefore says that only non-member and
-                                // static member functions can be used to
-                                // initialize a static variable. We use the
-                                // latter possibility by defining a function
-                                // <code>get_centers</code> that computes the list of
-                                // center points when called.
-                                //
-                                // Note that this class works just fine in
-                                // both 2d and 3d, with the only difference
-                                // being that we use more points in 3d: by
-                                // experimenting we find that we need more
-                                // exponentials in 3d than in 2d (we have
-                                // more ground to cover, after all, if we
-                                // want to keep the distance between centers
-                                // roughly equal), so we choose 40 in 2d and
-                                // 100 in 3d. For any other dimension, the
-                                // function does presently not know what to
-                                // do so simply throws an exception
-                                // indicating exactly this.
+                                 // @sect4{Random medium permeability}
+
+                                 // This function does as announced in the
+                                 // introduction, i.e. it creates an overlay
+                                 // of exponentials at random places. There is
+                                 // one thing worth considering for this
+                                 // class. The issue centers around the
+                                 // problem that the class creates the centers
+                                 // of the exponentials using a random
+                                 // function. If we therefore created the
+                                 // centers each time we create an object of
+                                 // the present type, we would get a different
+                                 // list of centers each time. That's not what
+                                 // we expect from classes of this type: they
+                                 // should reliably represent the same
+                                 // function.
+                                 //
+                                 // The solution to this problem is to make
+                                 // the list of centers a static member
+                                 // variable of this class, i.e. there exists
+                                 // exactly one such variable for the entire
+                                 // program, rather than for each object of
+                                 // this type. That's exactly what we are
+                                 // going to do.
+                                 //
+                                 // The next problem, however, is that we need
+                                 // a way to initialize this variable. Since
+                                 // this variable is initialized at the
+                                 // beginning of the program, we can't use a
+                                 // regular member function for that since
+                                 // there may not be an object of this type
+                                 // around at the time. The C++ standard
+                                 // therefore says that only non-member and
+                                 // static member functions can be used to
+                                 // initialize a static variable. We use the
+                                 // latter possibility by defining a function
+                                 // <code>get_centers</code> that computes the list of
+                                 // center points when called.
+                                 //
+                                 // Note that this class works just fine in
+                                 // both 2d and 3d, with the only difference
+                                 // being that we use more points in 3d: by
+                                 // experimenting we find that we need more
+                                 // exponentials in 3d than in 2d (we have
+                                 // more ground to cover, after all, if we
+                                 // want to keep the distance between centers
+                                 // roughly equal), so we choose 40 in 2d and
+                                 // 100 in 3d. For any other dimension, the
+                                 // function does presently not know what to
+                                 // do so simply throws an exception
+                                 // indicating exactly this.
 namespace RandomMedium
 {
   template <int dim>
@@ -393,7 +393,7 @@ namespace RandomMedium
   {
     public:
       virtual void value_list (const std::vector<Point<dim> > &points,
-                              std::vector<Tensor<2,dim> >    &values) const;
+                               std::vector<Tensor<2,dim> >    &values) const;
 
     private:
       static std::vector<Point<dim> > centers;
@@ -413,15 +413,15 @@ namespace RandomMedium
   KInverse<dim>::get_centers ()
   {
     const unsigned int N = (dim == 2 ?
-                           40 :
-                           (dim == 3 ?
-                            100 :
-                            throw ExcNotImplemented()));
+                            40 :
+                            (dim == 3 ?
+                             100 :
+                             throw ExcNotImplemented()));
   
     std::vector<Point<dim> > centers_list (N);
     for (unsigned int i=0; i<N; ++i)
       for (unsigned int d=0; d<dim; ++d)
-       centers_list[i][d] = static_cast<double>(rand())/RAND_MAX;
+        centers_list[i][d] = static_cast<double>(rand())/RAND_MAX;
 
     return centers_list;
   }
@@ -431,46 +431,46 @@ namespace RandomMedium
   template <int dim>
   void
   KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
-                            std::vector<Tensor<2,dim> >    &values) const
+                             std::vector<Tensor<2,dim> >    &values) const
   {
     Assert (points.size() == values.size(),
-           ExcDimensionMismatch (points.size(), values.size()));
+            ExcDimensionMismatch (points.size(), values.size()));
 
     for (unsigned int p=0; p<points.size(); ++p)
       {
-       values[p].clear ();
+        values[p].clear ();
 
-       double permeability = 0;
-       for (unsigned int i=0; i<centers.size(); ++i)
-         permeability += std::exp(-(points[p]-centers[i]).square()
-                                  / (0.05 * 0.05));
+        double permeability = 0;
+        for (unsigned int i=0; i<centers.size(); ++i)
+          permeability += std::exp(-(points[p]-centers[i]).square()
+                                   / (0.05 * 0.05));
       
-       const double normalized_permeability
-         = std::min (std::max(permeability, 0.01), 4.);
+        const double normalized_permeability
+          = std::min (std::max(permeability, 0.01), 4.);
       
-       for (unsigned int d=0; d<dim; ++d)
-         values[p][d][d] = 1./normalized_permeability;
+        for (unsigned int d=0; d<dim; ++d)
+          values[p][d][d] = 1./normalized_permeability;
       }
   }
 }
 
 
 
-                                // @sect3{The inverse mobility and saturation functions}
+                                 // @sect3{The inverse mobility and saturation functions}
 
-                                // There are two more pieces of data that we
-                                // need to describe, namely the inverse
-                                // mobility function and the saturation
-                                // curve. Their form is also given in the
-                                // introduction:
+                                 // There are two more pieces of data that we
+                                 // need to describe, namely the inverse
+                                 // mobility function and the saturation
+                                 // curve. Their form is also given in the
+                                 // introduction:
 double mobility_inverse (const double S,
-                        const double viscosity)
+                         const double viscosity)
 {
   return 1.0 /(1.0/viscosity * S * S + (1-S) * (1-S));
 }
 
 double f_saturation (const double S,
-                    const double viscosity)
+                     const double viscosity)
 {   
   return S*S /( S * S +viscosity * (1-S) * (1-S));
 }
@@ -481,17 +481,17 @@ double f_saturation (const double S,
 
                                  // @sect3{extract_u and friends}
 
-                                // More tools: We need methods to extract the
-                                // velocity, pressure, and saturation
-                                // components of finite element shape
-                                // functions. These functions here are
-                                // completely analogous to the ones we have
-                                // already used in step-20:
+                                 // More tools: We need methods to extract the
+                                 // velocity, pressure, and saturation
+                                 // components of finite element shape
+                                 // functions. These functions here are
+                                 // completely analogous to the ones we have
+                                 // already used in step-20:
 template <int dim>
 Tensor<1,dim>
 extract_u (const FEValuesBase<dim> &fe_values,
-          const unsigned int i,
-          const unsigned int q)
+           const unsigned int i,
+           const unsigned int q)
 {
   Tensor<1,dim> tmp;
 
@@ -506,8 +506,8 @@ extract_u (const FEValuesBase<dim> &fe_values,
 template <int dim>
 double
 extract_div_u (const FEValuesBase<dim> &fe_values,
-              const unsigned int i,
-              const unsigned int q)
+               const unsigned int i,
+               const unsigned int q)
 {
   double divergence = 0;
   for (unsigned int d=0; d<dim; ++d)
@@ -541,8 +541,8 @@ double extract_s (const FEValuesBase<dim> &fe_values,
 template <int dim>
 Tensor<1,dim>
 extract_grad_s (const FEValuesBase<dim> &fe_values,
-               const unsigned int i,
-               const unsigned int q)
+                const unsigned int i,
+                const unsigned int q)
 {
   Tensor<1,dim> tmp;
   for (unsigned int d=0; d<dim; ++d)
@@ -555,10 +555,10 @@ extract_grad_s (const FEValuesBase<dim> &fe_values,
 
                                  // @sect3{Linear solvers and preconditioners}
 
-                                // The linear solvers we use are also
-                                // completely analogous to the ones used in
-                                // step-20. The following classes are
-                                // therefore copied verbatim from there.
+                                 // The linear solvers we use are also
+                                 // completely analogous to the ones used in
+                                 // step-20. The following classes are
+                                 // therefore copied verbatim from there.
 template <class Matrix>
 class InverseMatrix : public Subscriptor
 {
@@ -617,7 +617,7 @@ class SchurComplement : public Subscriptor
 
 SchurComplement::
 SchurComplement (const BlockSparseMatrix<double> &A,
-                const InverseMatrix<SparseMatrix<double> > &Minv)
+                 const InverseMatrix<SparseMatrix<double> > &Minv)
                 :
                 system_matrix (&A),
                 m_inverse (&Minv),
@@ -674,30 +674,30 @@ void ApproximateSchurComplement::vmult (Vector<double>       &dst,
 
                                  // @sect3{<code>TwoPhaseFlowProblem</code> class implementation}
 
-                                // Here now the implementation of the main
-                                // class. Much of it is actually copied from
-                                // step-20, so we won't comment on it in much
-                                // detail. You should try to get familiar
-                                // with that program first, then most of what
-                                // is happening here should be mostly clear.
+                                 // Here now the implementation of the main
+                                 // class. Much of it is actually copied from
+                                 // step-20, so we won't comment on it in much
+                                 // detail. You should try to get familiar
+                                 // with that program first, then most of what
+                                 // is happening here should be mostly clear.
 
                                  // @sect4{TwoPhaseFlowProblem::TwoPhaseFlowProblem}
                                  // First for the constructor. We use $RT_k
-                                 // \times DG_k \times DG_k$ spaces. The time
+                                 // \times DQ_k \times DQ_k$ spaces. The time
                                  // step is set to zero initially, but will be
                                  // computed before it is needed first, as
                                  // described in a subsection of the
                                  // introduction.
 template <int dim>
 TwoPhaseFlowProblem<dim>::TwoPhaseFlowProblem (const unsigned int degree)
-               :
-               degree (degree),
+                :
+                degree (degree),
                 fe (FE_RaviartThomas<dim>(degree), 1,
                     FE_DGQ<dim>(degree), 1,
-                   FE_DGQ<dim>(degree), 1),
-               dof_handler (triangulation),
-               n_refinement_steps (5),
-               time_step (0),
+                    FE_DGQ<dim>(degree), 1),
+                dof_handler (triangulation),
+                n_refinement_steps (5),
+                time_step (0),
                 viscosity (0.2)
 {}
 
@@ -724,16 +724,16 @@ void TwoPhaseFlowProblem<dim>::make_grid_and_dofs ()
   DoFTools::count_dofs_per_component (dof_handler, dofs_per_component);  
   const unsigned int n_u = dofs_per_component[0],
                      n_p = dofs_per_component[dim],
-                    n_s = dofs_per_component[dim+1];
+                     n_s = dofs_per_component[dim+1];
 
   std::cout << "Number of active cells: "
-           << triangulation.n_active_cells()
-           << std::endl
+            << triangulation.n_active_cells()
+            << std::endl
             << "Number of degrees of freedom: "
-           << dof_handler.n_dofs()
+            << dof_handler.n_dofs()
             << " (" << n_u << '+' << n_p << '+'<< n_s <<')'
-           << std::endl
-           << std::endl;
+            << std::endl
+            << std::endl;
   
   const unsigned int
     n_couplings = dof_handler.max_couplings_between_dofs();
@@ -805,11 +805,11 @@ void TwoPhaseFlowProblem<dim>::assemble_system ()
   QGauss<dim-1> face_quadrature_formula(degree+2);
 
   FEValues<dim> fe_values (fe, quadrature_formula, 
-                          update_values    | update_gradients |
+                           update_values    | update_gradients |
                            update_q_points  | update_JxW_values);
   FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula, 
-                                   update_values    | update_normal_vectors |
-                                   update_q_points  | update_JxW_values);
+                                    update_values    | update_normal_vectors |
+                                    update_q_points  | update_JxW_values);
 
   const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
   
@@ -842,67 +842,67 @@ void TwoPhaseFlowProblem<dim>::assemble_system ()
       local_matrix = 0;
       local_rhs = 0;
 
-                                      // Here's the first significant
-                                      // difference: We have to get the
-                                      // values of the saturation function of
-                                      // the previous time step at the
-                                      // quadrature points. To this end, we
-                                      // can use the
-                                      // FEValues::get_function_values
-                                      // (previously already used in step-9,
-                                      // step-14 and step-15), a function
-                                      // that takes a solution vector and
-                                      // returns a list of function values at
-                                      // the quadrature points of the present
-                                      // cell. In fact, it returns the
-                                      // complete vector-valued solution at
-                                      // each quadrature point, i.e. not only
-                                      // the saturation but also the
-                                      // velocities and pressure:
+                                       // Here's the first significant
+                                       // difference: We have to get the
+                                       // values of the saturation function of
+                                       // the previous time step at the
+                                       // quadrature points. To this end, we
+                                       // can use the
+                                       // FEValues::get_function_values
+                                       // (previously already used in step-9,
+                                       // step-14 and step-15), a function
+                                       // that takes a solution vector and
+                                       // returns a list of function values at
+                                       // the quadrature points of the present
+                                       // cell. In fact, it returns the
+                                       // complete vector-valued solution at
+                                       // each quadrature point, i.e. not only
+                                       // the saturation but also the
+                                       // velocities and pressure:
       fe_values.get_function_values (old_solution, old_solution_values);
 
-                                      // Then we also have to get the values
-                                      // of the pressure right hand side and
-                                      // of the inverse permeability tensor
-                                      // at the quadrature points:
+                                       // Then we also have to get the values
+                                       // of the pressure right hand side and
+                                       // of the inverse permeability tensor
+                                       // at the quadrature points:
       pressure_right_hand_side.value_list (fe_values.get_quadrature_points(),
-                                          pressure_rhs_values);
+                                           pressure_rhs_values);
       k_inverse.value_list (fe_values.get_quadrature_points(),
                             k_inverse_values);
 
-                                      // With all this, we can now loop over
-                                      // all the quadrature points and shape
-                                      // functions on this cell and assemble
-                                      // those parts of the matrix and right
-                                      // hand side that we deal with in this
-                                      // function. The individual terms in
-                                      // the contributions should be
-                                      // self-explanatory given the explicit
-                                      // form of the bilinear form stated in
-                                      // the introduction:
+                                       // With all this, we can now loop over
+                                       // all the quadrature points and shape
+                                       // functions on this cell and assemble
+                                       // those parts of the matrix and right
+                                       // hand side that we deal with in this
+                                       // function. The individual terms in
+                                       // the contributions should be
+                                       // self-explanatory given the explicit
+                                       // form of the bilinear form stated in
+                                       // the introduction:
       for (unsigned int q=0; q<n_q_points; ++q)            
         for (unsigned int i=0; i<dofs_per_cell; ++i)
           {
-           const double old_s = old_solution_values[q](dim+1);
+            const double old_s = old_solution_values[q](dim+1);
 
             const Tensor<1,dim> phi_i_u      = extract_u (fe_values, i, q);
-           const double        div_phi_i_u  = extract_div_u (fe_values, i, q);
+            const double        div_phi_i_u  = extract_div_u (fe_values, i, q);
             const double        phi_i_p      = extract_p (fe_values, i, q);
-           const double        phi_i_s      = extract_s (fe_values, i, q); 
-           const Tensor<1,dim> grad_phi_i_s = extract_grad_s(fe_values, i, q);
+            const double        phi_i_s      = extract_s (fe_values, i, q); 
+            const Tensor<1,dim> grad_phi_i_s = extract_grad_s(fe_values, i, q);
             
             for (unsigned int j=0; j<dofs_per_cell; ++j)
               {
                 const Tensor<1,dim> phi_j_u     = extract_u (fe_values, j, q);
-               const double        div_phi_j_u = extract_div_u (fe_values, j, q);
+                const double        div_phi_j_u = extract_div_u (fe_values, j, q);
                 const double        phi_j_p     = extract_p (fe_values, j, q);
                 const double        phi_j_s     = extract_s (fe_values, j, q);
-               
+                
                 local_matrix(i,j) += (phi_i_u * k_inverse_values[q] *
-                                     mobility_inverse(old_s,viscosity) * phi_j_u
+                                      mobility_inverse(old_s,viscosity) * phi_j_u
                                       - div_phi_i_u * phi_j_p
                                       - phi_i_p * div_phi_j_u
-                                     + phi_i_s * phi_j_s)
+                                      + phi_i_s * phi_j_s)
                                      * fe_values.JxW(q);     
               }
 
@@ -911,32 +911,32 @@ void TwoPhaseFlowProblem<dim>::assemble_system ()
           }
       
 
-                                      // Next, we also have to deal with the
-                                      // pressure boundary values. This,
-                                      // again is as in step-20:
+                                       // Next, we also have to deal with the
+                                       // pressure boundary values. This,
+                                       // again is as in step-20:
       for (unsigned int face_no=0;
-          face_no<GeometryInfo<dim>::faces_per_cell;
-          ++face_no)
-       if (cell->at_boundary(face_no))
-         {
-           fe_face_values.reinit (cell, face_no);
-           
-           pressure_boundary_values
-             .value_list (fe_face_values.get_quadrature_points(),
-                          boundary_values);
-
-           for (unsigned int q=0; q<n_face_q_points; ++q) 
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
-               {
-                 const Tensor<1,dim>
-                   phi_i_u = extract_u (fe_face_values, i, q);
-
-                 local_rhs(i) += -(phi_i_u *
-                                   fe_face_values.normal_vector(q) *
-                                   boundary_values[q] *
-                                   fe_face_values.JxW(q));
-               }
-         }
+           face_no<GeometryInfo<dim>::faces_per_cell;
+           ++face_no)
+        if (cell->at_boundary(face_no))
+          {
+            fe_face_values.reinit (cell, face_no);
+            
+            pressure_boundary_values
+              .value_list (fe_face_values.get_quadrature_points(),
+                           boundary_values);
+
+            for (unsigned int q=0; q<n_face_q_points; ++q) 
+              for (unsigned int i=0; i<dofs_per_cell; ++i)
+                {
+                  const Tensor<1,dim>
+                    phi_i_u = extract_u (fe_face_values, i, q);
+
+                  local_rhs(i) += -(phi_i_u *
+                                    fe_face_values.normal_vector(q) *
+                                    boundary_values[q] *
+                                    fe_face_values.JxW(q));
+                }
+          }
 
                                        // The final step in the loop
                                        // over all cells is to
@@ -947,10 +947,10 @@ void TwoPhaseFlowProblem<dim>::assemble_system ()
       for (unsigned int i=0; i<dofs_per_cell; ++i)
       
         for (unsigned int j=0; j<dofs_per_cell; ++j)
-         {    system_matrix.add (local_dof_indices[i],
-                                 local_dof_indices[j],
-                                 local_matrix(i,j));
-         }
+          {    system_matrix.add (local_dof_indices[i],
+                                  local_dof_indices[j],
+                                  local_matrix(i,j));
+          }
       
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         system_rhs(local_dof_indices[i]) += local_rhs(i);
@@ -958,33 +958,33 @@ void TwoPhaseFlowProblem<dim>::assemble_system ()
 }
 
 
-                                // So much for assembly of matrix and right
-                                // hand side. Note that we do not have to
-                                // interpolate and apply boundary values
-                                // since they have all been taken care of in
-                                // the weak form already.
+                                 // So much for assembly of matrix and right
+                                 // hand side. Note that we do not have to
+                                 // interpolate and apply boundary values
+                                 // since they have all been taken care of in
+                                 // the weak form already.
 
 
-                                // @sect4{TwoPhaseFlowProblem::assemble_rhs_S}
+                                 // @sect4{TwoPhaseFlowProblem::assemble_rhs_S}
 
-                                // As explained in the introduction, we can
-                                // only evaluate the right hand side of the
-                                // saturation equation once the velocity has
-                                // been computed. We therefore have this
-                                // separate function to this end.
+                                 // As explained in the introduction, we can
+                                 // only evaluate the right hand side of the
+                                 // saturation equation once the velocity has
+                                 // been computed. We therefore have this
+                                 // separate function to this end.
 template <int dim>
 void TwoPhaseFlowProblem<dim>::assemble_rhs_S () 
 {  
   QGauss<dim>   quadrature_formula(degree+2); 
   QGauss<dim-1> face_quadrature_formula(degree+2);  
   FEValues<dim> fe_values (fe, quadrature_formula, 
-                          update_values    | update_gradients |
-                          update_q_points  | update_JxW_values);
+                           update_values    | update_gradients |
+                           update_q_points  | update_JxW_values);
   FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula, 
-                                   update_values    | update_normal_vectors |
-                                   update_q_points  | update_JxW_values);
+                                    update_values    | update_normal_vectors |
+                                    update_q_points  | update_JxW_values);
   FEFaceValues<dim> fe_face_values_neighbor (fe, face_quadrature_formula, 
-                                            update_values);
+                                             update_values);
  
   const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
   const unsigned int   n_q_points      = quadrature_formula.n_quadrature_points;
@@ -1014,107 +1014,107 @@ void TwoPhaseFlowProblem<dim>::assemble_rhs_S ()
       fe_values.get_function_values (old_solution, old_solution_values);
       fe_values.get_function_values (solution, present_solution_values);
 
-                                      // First for the cell terms. These are,
-                                      // following the formulas in the
-                                      // introduction, $(S^n,\sigma)-(F(S^n)
-                                      // \mathbf{v}^{n+1},\nabla sigma)$,
-                                      // where $\sigma$ is the saturation
-                                      // component of the test function:
+                                       // First for the cell terms. These are,
+                                       // following the formulas in the
+                                       // introduction, $(S^n,\sigma)-(F(S^n)
+                                       // \mathbf{v}^{n+1},\nabla sigma)$,
+                                       // where $\sigma$ is the saturation
+                                       // component of the test function:
       for (unsigned int q=0; q<n_q_points; ++q) 
-       for (unsigned int i=0; i<dofs_per_cell; ++i)
-         {
-           const double old_s = old_solution_values[q](dim+1);
-           Tensor<1,dim> present_u;
-           for (unsigned int d=0; d<dim; ++d)
-             present_u[d] = present_solution_values[q](d);
-
-           const double        phi_i_s      = extract_s(fe_values, i, q);
-           const Tensor<1,dim> grad_phi_i_s = extract_grad_s(fe_values, i, q);
-                    
-           local_rhs(i) += (time_step *
-                            f_saturation(old_s,viscosity) *
-                            present_u *
-                            grad_phi_i_s
-                            +
-                            old_s * phi_i_s)
-                           *
-                           fe_values.JxW(q);
-         }
-
-                                      // Secondly, we have to deal with the
-                                      // flux parts on the face
-                                      // boundaries. This was a bit more
-                                      // involved because we first have to
-                                      // determine which are the influx and
-                                      // outflux parts of the cell
-                                      // boundary. If we have an influx
-                                      // boundary, we need to evaluate the
-                                      // saturation on the other side of the
-                                      // face (or the boundary values, if we
-                                      // are at the boundary of the domain).
-                                      //
-                                      // All this is a bit tricky, but has
-                                      // been explained in some detail
-                                      // already in step-9. Take a look there
-                                      // how this is supposed to work!
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          {
+            const double old_s = old_solution_values[q](dim+1);
+            Tensor<1,dim> present_u;
+            for (unsigned int d=0; d<dim; ++d)
+              present_u[d] = present_solution_values[q](d);
+
+            const double        phi_i_s      = extract_s(fe_values, i, q);
+            const Tensor<1,dim> grad_phi_i_s = extract_grad_s(fe_values, i, q);
+                     
+            local_rhs(i) += (time_step *
+                             f_saturation(old_s,viscosity) *
+                             present_u *
+                             grad_phi_i_s
+                             +
+                             old_s * phi_i_s)
+                            *
+                            fe_values.JxW(q);
+          }
+
+                                       // Secondly, we have to deal with the
+                                       // flux parts on the face
+                                       // boundaries. This was a bit more
+                                       // involved because we first have to
+                                       // determine which are the influx and
+                                       // outflux parts of the cell
+                                       // boundary. If we have an influx
+                                       // boundary, we need to evaluate the
+                                       // saturation on the other side of the
+                                       // face (or the boundary values, if we
+                                       // are at the boundary of the domain).
+                                       //
+                                       // All this is a bit tricky, but has
+                                       // been explained in some detail
+                                       // already in step-9. Take a look there
+                                       // how this is supposed to work!
       for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell;
-          ++face_no)
-       {
-         fe_face_values.reinit (cell, face_no);
-
-         fe_face_values.get_function_values (old_solution, old_solution_values_face);
-         fe_face_values.get_function_values (solution, present_solution_values_face);
-
-         if (cell->at_boundary(face_no))
-           saturation_boundary_values
-             .value_list (fe_face_values.get_quadrature_points(),
-                          neighbor_saturation);
-         else
-           {
-             const typename DoFHandler<dim>::active_cell_iterator
-               neighbor = cell->neighbor(face_no);
-             const unsigned int
-               neighbor_face = cell->neighbor_of_neighbor(face_no);
-
-             fe_face_values_neighbor.reinit (neighbor, neighbor_face);
-            
-             fe_face_values_neighbor
-               .get_function_values (old_solution,
-                                     old_solution_values_face_neighbor);
-            
-             for (unsigned int q=0; q<n_face_q_points; ++q)
-               neighbor_saturation[q] = old_solution_values_face_neighbor[q](dim+1);
-           }
+           ++face_no)
+        {
+          fe_face_values.reinit (cell, face_no);
+
+          fe_face_values.get_function_values (old_solution, old_solution_values_face);
+          fe_face_values.get_function_values (solution, present_solution_values_face);
+
+          if (cell->at_boundary(face_no))
+            saturation_boundary_values
+              .value_list (fe_face_values.get_quadrature_points(),
+                           neighbor_saturation);
+          else
+            {
+              const typename DoFHandler<dim>::active_cell_iterator
+                neighbor = cell->neighbor(face_no);
+              const unsigned int
+                neighbor_face = cell->neighbor_of_neighbor(face_no);
+
+              fe_face_values_neighbor.reinit (neighbor, neighbor_face);
+             
+              fe_face_values_neighbor
+                .get_function_values (old_solution,
+                                      old_solution_values_face_neighbor);
+             
+              for (unsigned int q=0; q<n_face_q_points; ++q)
+                neighbor_saturation[q] = old_solution_values_face_neighbor[q](dim+1);
+            }
           
 
-         for (unsigned int q=0; q<n_face_q_points; ++q)
-           {
-             Tensor<1,dim> present_u_face;
-             for (unsigned int d=0; d<dim; ++d)
-               present_u_face[d] = present_solution_values_face[q](d);
-
-             const double normal_flux = present_u_face *
-                                        fe_face_values.normal_vector(q);
-
-             const bool is_outflow_q_point = (normal_flux >= 0);
-                                    
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
-               local_rhs(i) -= time_step *
-                               normal_flux *
-                               f_saturation((is_outflow_q_point == true
-                                             ?
-                                             old_solution_values_face[q](dim+1)
-                                             :
-                                             neighbor_saturation[q]),
-                                            viscosity) *
-                               extract_s(fe_face_values,i,q) *
-                               fe_face_values.JxW(q);
-           }
-       }
+          for (unsigned int q=0; q<n_face_q_points; ++q)
+            {
+              Tensor<1,dim> present_u_face;
+              for (unsigned int d=0; d<dim; ++d)
+                present_u_face[d] = present_solution_values_face[q](d);
+
+              const double normal_flux = present_u_face *
+                                         fe_face_values.normal_vector(q);
+
+              const bool is_outflow_q_point = (normal_flux >= 0);
+                                     
+              for (unsigned int i=0; i<dofs_per_cell; ++i)
+                local_rhs(i) -= time_step *
+                                normal_flux *
+                                f_saturation((is_outflow_q_point == true
+                                              ?
+                                              old_solution_values_face[q](dim+1)
+                                              :
+                                              neighbor_saturation[q]),
+                                             viscosity) *
+                                extract_s(fe_face_values,i,q) *
+                                fe_face_values.JxW(q);
+            }
+        }
   
       cell->get_dof_indices (local_dof_indices);
       for (unsigned int i=0; i<dofs_per_cell; ++i)
-       system_rhs(local_dof_indices[i]) += local_rhs(i);       
+        system_rhs(local_dof_indices[i]) += local_rhs(i);       
     }
 } 
 
@@ -1137,9 +1137,9 @@ void TwoPhaseFlowProblem<dim>::solve ()
   Vector<double> tmp2 (solution.block(2).size());
   
 
-                                  // First the pressure, using the pressure
-                                  // Schur complement of the first two
-                                  // equations:
+                                   // First the pressure, using the pressure
+                                   // Schur complement of the first two
+                                   // equations:
   {
     m_inverse.vmult (tmp, system_rhs.block(0));
     system_matrix.block(1,0).vmult (schur_rhs, tmp);
@@ -1157,14 +1157,14 @@ void TwoPhaseFlowProblem<dim>::solve ()
 
     
     SolverControl solver_control (system_matrix.block(0,0).m(),
-                                 1e-12*schur_rhs.l2_norm());
+                                  1e-12*schur_rhs.l2_norm());
     SolverCG<>    cg (solver_control);
 
     cg.solve (schur_complement, solution.block(1), schur_rhs,
               preconditioner);
   
     std::cout << "   "
-             << solver_control.last_step()
+              << solver_control.last_step()
               << " CG Schur complement iterations for pressure."
               << std::endl;
   }
@@ -1178,43 +1178,43 @@ void TwoPhaseFlowProblem<dim>::solve ()
     m_inverse.vmult (solution.block(0), tmp);
   }
 
-                                  // Finally, we have to take care of the
-                                  // saturation equation. The first business
-                                  // we have here is to determine the time
-                                  // step using the formula in the
-                                  // introduction. Knowing the shape of our
-                                  // domain and that we created the mesh by
-                                  // regular subdivision of cells, we can
-                                  // compute the diameter of each of our
-                                  // cells quite easily (in fact we use the
-                                  // linear extensions in coordinate
-                                  // directions of the cells, not the
-                                  // diameter). The maximal velocity we
-                                  // compute using a helper function defined
-                                  // below:
+                                   // Finally, we have to take care of the
+                                   // saturation equation. The first business
+                                   // we have here is to determine the time
+                                   // step using the formula in the
+                                   // introduction. Knowing the shape of our
+                                   // domain and that we created the mesh by
+                                   // regular subdivision of cells, we can
+                                   // compute the diameter of each of our
+                                   // cells quite easily (in fact we use the
+                                   // linear extensions in coordinate
+                                   // directions of the cells, not the
+                                   // diameter). The maximal velocity we
+                                   // compute using a helper function defined
+                                   // below:
   time_step = std::pow(0.5, double(n_refinement_steps)) /
-             get_maximal_velocity();
+              get_maximal_velocity();
 
-                                  // The next step is to assemble the right
-                                  // hand side, and then to pass everything
-                                  // on for solution. At the end, we project
-                                  // back saturations onto the physically
-                                  // reasonable range:
+                                   // The next step is to assemble the right
+                                   // hand side, and then to pass everything
+                                   // on for solution. At the end, we project
+                                   // back saturations onto the physically
+                                   // reasonable range:
   assemble_rhs_S ();
   {
     
     SolverControl solver_control (system_matrix.block(2,2).m(),
-                                 1e-8*system_rhs.block(2).l2_norm());
+                                  1e-8*system_rhs.block(2).l2_norm());
     SolverCG<>   cg (solver_control);
     cg.solve (system_matrix.block(2,2), solution.block(2), system_rhs.block(2),
-             PreconditionIdentity());
-               
+              PreconditionIdentity());
+                
     project_back_saturation ();
-       
+        
     std::cout << "   "
-             << solver_control.last_step()
+              << solver_control.last_step()
               << " CG iterations for saturation."
-              << std::endl;            
+              << std::endl;             
   } 
 
    
@@ -1224,10 +1224,16 @@ void TwoPhaseFlowProblem<dim>::solve ()
 
                                  // @sect4{TwoPhaseFlowProblem::output_results}
 
-                                 // There is nothing surprising here:
+                                 // There is nothing surprising here. Since
+                                 // the program will do a lot of time steps,
+                                 // we create an output file only every fifth
+                                 // time step.
 template <int dim>
 void TwoPhaseFlowProblem<dim>::output_results ()  const
-{  
+{
+  if (timestep_number % 5 != 0)
+    return;
+  
   std::vector<std::string> solution_names;
   switch (dim)
     {
@@ -1235,7 +1241,7 @@ void TwoPhaseFlowProblem<dim>::output_results ()  const
             solution_names.push_back ("u");
             solution_names.push_back ("v");
             solution_names.push_back ("p");
-           solution_names.push_back ("S");
+            solution_names.push_back ("S");
             break;
             
       case 3:
@@ -1243,7 +1249,7 @@ void TwoPhaseFlowProblem<dim>::output_results ()  const
             solution_names.push_back ("v");
             solution_names.push_back ("w");
             solution_names.push_back ("p");
-           solution_names.push_back ("S");
+            solution_names.push_back ("S");
             break;
             
       default:
@@ -1266,27 +1272,53 @@ void TwoPhaseFlowProblem<dim>::output_results ()  const
 
 
 
+                                 // @sect4{TwoPhaseFlowProblem::project_back_saturation}
+
+                                 // In this function, we simply run over all
+                                 // saturation degrees of freedom and make
+                                 // sure that if they should have left the
+                                 // physically reasonable range, that they be
+                                 // reset to the interval $[0,1]$. To do this,
+                                 // we only have to loop over all saturation
+                                 // components of the solution vector; these
+                                 // are stored in the block 2 (block 0 are the
+                                 // velocities, block 1 are the pressures).
+                                 //
+                                 // It may be instructive to note that this
+                                 // function almost never triggers when the
+                                 // time step is chosen as mentioned in the
+                                 // introduction. However, if we choose the
+                                 // timestep only slightly larger, we get
+                                 // plenty of values outside the proper
+                                 // range. Strictly speaking, the function is
+                                 // therefore unnecessary if we choose the
+                                 // time step small enough. In a sense, the
+                                 // function is therefore only a safety device
+                                 // to avoid situations where our entire
+                                 // solution becomes unphysical because
+                                 // individual degrees of freedom have become
+                                 // unphysical a few time steps earlier.
 template <int dim>
 void
 TwoPhaseFlowProblem<dim>::project_back_saturation ()
 {
   for (unsigned int i=0; i<solution.block(2).size(); ++i)
     if (solution.block(2)(i) < 0)
-      {
-       std::cout << "xxx       " << solution.block(2)(i) << std::endl;
-       solution.block(2)(i) = 0;
-      }
-  
+      solution.block(2)(i) = 0;
     else
       if (solution.block(2)(i) > 1)
-       {
-         std::cout << "xxx       " << solution.block(2)(i) << std::endl;
-         solution.block(2)(i) = 1;
-       }
+        solution.block(2)(i) = 1;
 }
 
 
+                                 // @sect4{TwoPhaseFlowProblem::get_maximal_velocity}
 
+                                 // The following function is used in
+                                 // determining the maximal allowable time
+                                 // step. What it does is to loop over all
+                                 // quadrature points in the domain and find
+                                 // what the maximal magnitude of the velocity
+                                 // is.
 template <int dim>
 double
 TwoPhaseFlowProblem<dim>::get_maximal_velocity () const
@@ -1296,9 +1328,9 @@ TwoPhaseFlowProblem<dim>::get_maximal_velocity () const
     = quadrature_formula.n_quadrature_points;
 
   FEValues<dim> fe_values (fe, quadrature_formula, 
-                          update_values);
+                           update_values);
   std::vector<Vector<double> > solution_values(n_q_points,
-                                              Vector<double>(dim+2));
+                                               Vector<double>(dim+2));
   double max_velocity = 0;
   
   typename DoFHandler<dim>::active_cell_iterator
@@ -1310,14 +1342,14 @@ TwoPhaseFlowProblem<dim>::get_maximal_velocity () const
       fe_values.get_function_values (solution, solution_values);
 
       for (unsigned int q=0; q<n_q_points; ++q)
-       {
-         Tensor<1,dim> velocity;
-         for (unsigned int i=0; i<dim; ++i)
-           velocity[i] = solution_values[q](i);          
-         
-         max_velocity = std::max (max_velocity,
-                                  velocity.norm());
-       }
+        {
+          Tensor<1,dim> velocity;
+          for (unsigned int i=0; i<dim; ++i)
+            velocity[i] = solution_values[q](i);          
+          
+          max_velocity = std::max (max_velocity,
+                                   velocity.norm());
+        }
     }
 
   return max_velocity;
@@ -1326,23 +1358,51 @@ TwoPhaseFlowProblem<dim>::get_maximal_velocity () const
 
                                  // @sect4{TwoPhaseFlowProblem::run}
 
-                                 // This is the final function of our
-                                 // main class. It's only job is to
-                                 // call the other functions in their order:
+                                 // This is the final function of our main
+                                 // class. Its brevity speaks for
+                                 // itself. There are only two points worth
+                                 // noting: First, the function projects the
+                                 // initial values onto the finite element
+                                 // space at the beginning; the
+                                 // VectorTools::project function doing this
+                                 // requires an argument indicating the
+                                 // hanging node constraints. We have none in
+                                 // this program (we compute on a uniformly
+                                 // refined mesh), but the function requires
+                                 // the argument anyway, of course. So we have
+                                 // to create a constraint object. In its
+                                 // original state, constraint objects are
+                                 // unsorted, and have to be sorted (using the
+                                 // ConstraintMatrix::close function) before
+                                 // they can be used. This is what we do here,
+                                 // and which is why we can't simply call the
+                                 // VectorTools::project function with an
+                                 // anonymous temporary object
+                                 // <code>ConstraintMatrix()</code> as the
+                                 // second argument.
+                                 //
+                                 // The second point worth mentioning is that
+                                 // we only compute the length of the present
+                                 // time step in the middle of solving the
+                                 // linear system corresponding to each time
+                                 // step. We can therefore output the present
+                                 // end time of a time step only at the end of
+                                 // the time step.
 template <int dim>
 void TwoPhaseFlowProblem<dim>::run () 
 {
   make_grid_and_dofs();
   
-  ConstraintMatrix constraints;
-  constraints.close();
-
-  std::list<double> production_rate;
-  std::list<double> production_time;
-
-  Vector<double> tmp (old_solution.size());
-  VectorTools::project (dof_handler, constraints, QGauss<dim>(degree+2),InitialValues<dim>(),tmp);
-  std::copy (tmp.begin(), tmp.end(), old_solution.begin());
+  {
+    ConstraintMatrix constraints;
+    constraints.close();
+    
+    VectorTools::project (dof_handler,
+                          constraints,
+                          QGauss<dim>(degree+2),
+                          InitialValues<dim>(),
+                          old_solution);
+  }
   
   timestep_number = 1;
   double time = 0;
@@ -1350,7 +1410,7 @@ void TwoPhaseFlowProblem<dim>::run ()
   do
     { 
       std::cout << "Timestep " << timestep_number
-               << std::endl; 
+                << std::endl; 
 
       assemble_system ();
 
@@ -1361,9 +1421,9 @@ void TwoPhaseFlowProblem<dim>::run ()
       time += time_step;
       ++timestep_number;
       std::cout << "   Now at t=" << time
-               << ", dt=" << time_step << '.'
-               << std::endl
-               << std::endl;
+                << ", dt=" << time_step << '.'
+                << std::endl
+                << std::endl;
     }
   while (time <= 250);
 }
@@ -1371,12 +1431,13 @@ void TwoPhaseFlowProblem<dim>::run ()
     
                                  // @sect3{The <code>main</code> function}
 
-                                // In the main function, we pass the degree
-                                // of the finite element space to the
-                                // constructor of the TwoPhaseFlowProblem.
-                                // Here, we use zero-th degree elements,
-                                // i.e. $RT_0\times DQ_0 \times DQ_0$. The
-                                // rest is as in all the other programs.
+                                 // That's it. In the main function, we pass
+                                 // the degree of the finite element space to
+                                 // the constructor of the TwoPhaseFlowProblem
+                                 // object.  Here, we use zero-th degree
+                                 // elements, i.e. $RT_0\times DQ_0 \times
+                                 // DQ_0$. The rest is as in all the other
+                                 // programs.
 int main () 
 {
   try
@@ -1389,25 +1450,25 @@ int main ()
   catch (std::exception &exc)
     {
       std::cerr << std::endl << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
+                << "----------------------------------------------------"
+                << std::endl;
       std::cerr << "Exception on processing: " << std::endl
-               << exc.what() << std::endl
-               << "Aborting!" << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
       
       return 1;
     }
   catch (...) 
     {
       std::cerr << std::endl << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
+                << "----------------------------------------------------"
+                << std::endl;
       std::cerr << "Unknown exception!" << std::endl
-               << "Aborting!" << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
       return 1;
     }
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.