From ebb63b7332467b15f5cb45d8b0e88e4f3c6ff49f Mon Sep 17 00:00:00 2001
From: hartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Tue, 20 Nov 2001 16:23:23 +0000
Subject: [PATCH] Complete program, and major part of documentation.

git-svn-id: https://svn.dealii.org/trunk@5220 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/examples/step-12/step-12.cc | 1568 +++++++++++++++++++++------
 1 file changed, 1220 insertions(+), 348 deletions(-)

diff --git a/deal.II/examples/step-12/step-12.cc b/deal.II/examples/step-12/step-12.cc
index 741eb6c9ab..5bea1272c8 100644
--- a/deal.II/examples/step-12/step-12.cc
+++ b/deal.II/examples/step-12/step-12.cc
@@ -21,28 +21,60 @@
 #include <numerics/data_out.h>
 #include <grid/grid_out.h>
 #include <grid/grid_refinement.h>
-#include <numerics/error_estimator.h>
 
+				 // This is the first new file. It
+				 // declares the MappingQ1 class that
+				 // gives the standard bilinear
+				 // mapping. For bilinear mappings use
+				 // an object of this class rather
+				 // than an object of the MappingQ(1)
+				 // class, as the MappingQ1 class is
+				 // optimized due to the
+				 // pre-knowledge of the actual
+				 // polynomial degree 1.
 #include <fe/mapping_q1.h>
+
+				 // Here the discontinuous finite
+				 // elements are defined. They are
+				 // used as all other finite elements.
 #include <fe/fe_dgq.h>
-#include <lac/precondition_block.h>
-#include <lac/solver_richardson.h>
 
+				 // We are going to use the simplest
+				 // possible solver, called richardson
+				 // iteration, that represents a simple
+				 // defect correction. This, in
+				 // combination with a block SSOR
+				 // preconditioner (defined in
+				 // precondition_block.h), that uses
+				 // the special block matrix structur
+				 // of system matrices arising from DG
+				 // discretizations.
+#include <lac/solver_richardson.h>
+#include <lac/precondition_block.h>
 
-#include <fstream>
+				 // We are going to use gradients as
+				 // refinement indicator.
+#include <numerics/derivative_approximation.h>
 
 
-template <int dim>
-class Beta
-{
-  public:
-    Beta () {};
+				 // Finally we do some time comparison
+				 // using the ``Timer'' class.
+#include <base/timer.h>
 
-    void value_list (const std::vector<Point<dim> > &points,
-		     std::vector<Point<dim> > &values) const;
-};
+				 // And this again is C++:
+#include <fstream>
 
 
+				 // First we define the class
+				 // representing the equation-specific
+				 // functions. Both classes, ``RHS''
+				 // and ``BoundaryValues'', are
+				 // derived from the Function
+				 // class. Only the ``value_list''
+				 // function are implemented because
+				 // only lists of function values are
+				 // computed rather than single
+				 // values.
 template <int dim>
 class RHS:  public Function<dim>
 {
@@ -56,90 +88,66 @@ class RHS:  public Function<dim>
 
 
 template <int dim>
-class BoundaryFunction:  public Function<dim>
+class BoundaryValues:  public Function<dim>
 {
   public:
-    BoundaryFunction() {};
+    BoundaryValues() {};
     
     virtual void value_list (const std::vector<Point<dim> > &points,
 			     std::vector<double> &values,
 			     const unsigned int component=0) const;
 };
 
+
+				 // The class ``Beta'' that represents
+				 // the vector valued flow field of
+				 // the linear transport equation is
+				 // not derived from the Function
+				 // class as we prefer to get function
+				 // values of type ``Point'' rather
+				 // than of type
+				 // ``Vector<double>''. This, because
+				 // there exist scalar products
+				 // between ``Point'' and ``Point'' as
+				 // well as between ``Point'' and
+				 // ``Tensor'', simplifying terms like
+				 // $\beta\cdot n$ and
+				 // $\beta\cdot\nabla v$.
 template <int dim>
-class DGAssembler
+class Beta
 {
   public:
-    DGAssembler() {};
+    Beta () {};
 
-    void assemble_cell_term(const FEValuesBase<dim>& fe_v,
-			    FullMatrix<double> &cell_matrix,
-			    Vector<double> &cell_vector);
-    
-    void assemble_face_term(const FEFaceValuesBase<dim>& fe_v,
-			    const FEFaceValuesBase<dim>& fe_v_neighbor,
-			    FullMatrix<double> &cell_matrix,
-			    FullMatrix<double> &cell_inflow_matrix,
-			    Vector<double> &cell_vector);
-    
-  private:
-    Beta<dim> beta_function;
-    RHS<dim> rhs_function;
-    BoundaryFunction<dim> boundary_function;
+    void value_list (const std::vector<Point<dim> > &points,
+		     std::vector<Point<dim> > &values) const;
 };
 
-				 // The main class is again almost
-				 // unchanged. Two additions, however,
-				 // are made: we have added the
-				 // ``refine'' function, which is used
-				 // to adaptively refine the grid
-				 // (instead of the global refinement
-				 // in the previous examples), and a
-				 // variable which will hold the
-				 // constraints associated to the
-				 // hanging nodes.
+
+				 // The implementation of the
+				 // ``value_list'' functions of these
+				 // classes are rather simple.  For
+				 // simplicity the right hand side is
+				 // set to be zero.
 template <int dim>
-class TransportProblem
+void RHS<dim>::value_list(const std::vector<Point<dim> > &,
+			  std::vector<double> &values,
+			  const unsigned int) const
 {
-  public:
-    TransportProblem ();
-    ~TransportProblem ();
-
-    void run ();
-    
-  private:
-    void setup_system ();
-    void assemble_system ();
-    void solve ();
-    void refine_grid ();
-    void output_results (const unsigned int cycle) const;
-
-    Triangulation<dim>   triangulation;
-    MappingQ1<dim>       mapping;
-    
-				     // We need a finite element
-				     // again. This time, we will want
-				     // to use quadratic polynomials
-				     // (but this is only specified in
-				     // the constructor):
-    FE_DGQ<dim>          fe;
-    DoFHandler<dim>      dof_handler;
-
-    SparsityPattern      sparsity_pattern;
-    SparseMatrix<double> system_matrix;
-    
-    Vector<double>       solution;
-    Vector<double>       right_hand_side;
-
-    DGAssembler<dim>     dg_assembler;
-};
-
+  for (unsigned int i=0; i<values.size(); ++i)
+    values[i]=0;
+}
 
+				 // The flow field is chosen to be
+				 // circular, anticlockwise, and with
+				 // the origin as midpoint.
 template <>
 void Beta<2>::value_list(const std::vector<Point<2> > &points,
 			 std::vector<Point<2> > &values) const
 {
-  Assert(values.size()==points.size(), ExcDimensionMismatch(values.size(),points.size()));
+  Assert(values.size()==points.size(),
+	 ExcDimensionMismatch(values.size(),points.size()));
+
   for (unsigned int i=0; i<points.size(); ++i)
     {
       const Point<2> &p=points[i];
@@ -151,27 +159,23 @@ void Beta<2>::value_list(const std::vector<Point<2> > &points,
     }
 }
 
-
-
-
-template <int dim>
-void RHS<dim>::value_list(const std::vector<Point<dim> > &,
-			  std::vector<double> &values,
-			  const unsigned int) const
-{
-  for (unsigned int i=0; i<values.size(); ++i)
-    values[i]=0;
-}
-
-
-
-
+				 // Hence the inflow boundary of the
+				 // unit square [0,1]^2 are the right
+				 // and the lower boundaries. We
+				 // prescribe discontinuous boundary
+				 // values 1 and 0 on the x-axis and
+				 // value 0 on the right boundary. The
+				 // values of this function on the
+				 // outflow boundaries will not be
+				 // used within the DG scheme.
 template <int dim>
-void BoundaryFunction<dim>::value_list(const std::vector<Point<dim> > &points,
+void BoundaryValues<dim>::value_list(const std::vector<Point<dim> > &points,
 				       std::vector<double> &values,
 				       const unsigned int) const
 {
-  Assert(values.size()==points.size(), ExcDimensionMismatch(values.size(),points.size()));
+  Assert(values.size()==points.size(),
+	 ExcDimensionMismatch(values.size(),points.size()));
+
   for (unsigned int i=0; i<values.size(); ++i)
     {
       if (points[i](0)<0.5)
@@ -181,29 +185,87 @@ void BoundaryFunction<dim>::value_list(const std::vector<Point<dim> > &points,
     }
 }
 
+				 // Next we define the equation-
+				 // dependent and DG-method-dependent
+				 // class ``DGTransportEquation''. Its
+				 // member functions were already
+				 // mentioned in the Introduction and
+				 // will be explained
+				 // below. Furthermore it includes
+				 // objects of the previously defined
+				 // ``Beta'', ``RHS'' and
+				 // ``BoundaryValues'' function
+				 // classes.
+template <int dim>
+class DGTransportEquation
+{
+  public:
+    DGTransportEquation() {};
 
+    void assemble_cell_term(const FEValues<dim>& fe_v,
+			    FullMatrix<double> &u_v_matrix,
+			    Vector<double> &cell_vector);
+    
+    void assemble_face_term1(const FEFaceValuesBase<dim>& fe_v,
+			     const FEFaceValuesBase<dim>& fe_v_neighbor,
+			     FullMatrix<double> &u_v_matrix,
+			     FullMatrix<double> &un_v_matrix,
+			     Vector<double> &cell_vector);
+
+    void assemble_face_term2(const FEFaceValuesBase<dim>& fe_v,
+			     const FEFaceValuesBase<dim>& fe_v_neighbor,
+			     FullMatrix<double> &u_v_matrix,
+			     FullMatrix<double> &un_v_matrix,
+			     FullMatrix<double> &u_vn_matrix,
+			     FullMatrix<double> &un_vn_matrix,
+			     Vector<double> &cell_vector);
+  private:
+    Beta<dim> beta_function;
+    RHS<dim> rhs_function;
+    BoundaryValues<dim> boundary_function;
+};
 
-
+				 // ``u_v_matrix'' is a cell matrix,
+				 // i.e. for a DG method of degree 1,
+				 // it is of size 4 times 4, and
+				 // ``cell_vector'' is of size 4.
+				 // When this function is invoked,
+				 // ``fe_v'' was reinited with the
+				 // current cell before and includes
+				 // all shape values needed.
 template <int dim>
-void DGAssembler<dim>::assemble_cell_term(const FEValuesBase<dim>& fe_v,
-					  FullMatrix<double> &cell_matrix,
-					  Vector<double> &cell_vector)
+void DGTransportEquation<dim>::assemble_cell_term(
+  const FEValues<dim>& fe_v,
+  FullMatrix<double> &u_v_matrix,
+  Vector<double> &cell_vector)
 {
+				   // First we ask ``fe_v'' for the
+				   // shape grads, shape values and
+				   // quadrature weights,
   const vector<vector<Tensor<1,2> > > &grad_v = fe_v.get_shape_grads ();
   const FullMatrix<double> &v = fe_v.get_shape_values ();
   const vector<double> &JxW = fe_v.get_JxW_values ();
 
+				   // Then the flow field beta and the
+				   // ``rhs_function'' are evaluated at
+				   // the quadrature points,
   vector<Point<dim> > beta (fe_v.n_quadrature_points);
   vector<double> rhs (fe_v.n_quadrature_points);
   
   beta_function.value_list (fe_v.get_quadrature_points(), beta);
   rhs_function.value_list (fe_v.get_quadrature_points(), rhs);
   
+				   // and the cell matrix and cell
+				   // vector are assembled as in
+				   // previous tutorial steps.  Here,
+				   // the terms $-(u,\beta\cdot\nabla
+				   // v)_K$ and $(f,v)_K$ are
+				   // assembled.
   for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
     for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) 
       {
 	for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
-	  cell_matrix(i,j) -= beta[point]*grad_v[i][point]*
+	  u_v_matrix(i,j) -= beta[point]*grad_v[i][point]*
 			      v(j,point) *
 			      JxW[point];
 	
@@ -212,40 +274,89 @@ void DGAssembler<dim>::assemble_cell_term(const FEValuesBase<dim>& fe_v,
 }
 
 
+				 // The ``assemble_face_term1''
+				 // function assembles the face terms
+				 // corresponding to the first version
+				 // of the DG method, cf. above. Then,
+				 // the face terms are given as a sum
+				 // of integrals over all cell
+				 // boundaries.
+				 //
+				 // When this function is invoked,
+				 // ``fe_v'' and ``fe_v_neighbor'' are
+				 // already reinited with the current
+				 // cell and the neighoring cell,
+				 // respectively, as well as with the
+				 // current face. Hence they provide
+				 // the inner and outer shape values
+				 // on the face.
+				 //
+				 // In addition to the cell matrix
+				 // ``u_v_matrix'' and the
+				 // ``cell_vector'' this function has
+				 // got a new argument
+				 // ``un_v_matrix'', that stores
+				 // contributions to the system matrix
+				 // that are based on outer values of
+				 // u, see $\hat u_h$ in the
+				 // Introduction, and inner values of
+				 // v, see $v_h$. Here we note that
+				 // ``un'' is the short notation for
+				 // ``u_neighbor'' and represents
+				 // $\hat u_h$.
 template <int dim>
-void DGAssembler<dim>::assemble_face_term(const FEFaceValuesBase<dim>& fe_v,
-					  const FEFaceValuesBase<dim>& fe_v_neighbor,      
-					  FullMatrix<double> &cell_matrix,
-					  FullMatrix<double> &cell_inflow_matrix,
-					  Vector<double> &cell_vector)
+void DGTransportEquation<dim>::assemble_face_term1(
+  const FEFaceValuesBase<dim>& fe_v,
+  const FEFaceValuesBase<dim>& fe_v_neighbor,      
+  FullMatrix<double> &u_v_matrix,
+  FullMatrix<double> &un_v_matrix,
+  Vector<double> &cell_vector)
 {
-  DoFHandler<dim>::face_iterator face=fe_v.get_face();
-  
+				   // Again, we ask the FEValues
+				   // objects for the shape values and
+				   // the quadrature weights
   const FullMatrix<double> &v = fe_v.get_shape_values ();
   const FullMatrix<double> &v_neighbor = fe_v_neighbor.get_shape_values ();  
   const vector<double> &JxW = fe_v.get_JxW_values ();
+				   // but also for the normals.
   const vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
 
+				   // We also evaluate the flow field
+				   // at the quadrature points
   vector<Point<dim> > beta (fe_v.n_quadrature_points);
-  vector<double> g(fe_v.n_quadrature_points);
   
   beta_function.value_list (fe_v.get_quadrature_points(), beta);
 
+				   // and the boundary values if the
+				   // current face belongs to the
+				   // boundary.
+  vector<double> g(fe_v.n_quadrature_points);
+  DoFHandler<dim>::face_iterator face=fe_v.get_face();
   if (face->at_boundary())
     boundary_function.value_list (fe_v.get_quadrature_points(), g);
 
+				   // Then we assemble the cell matrix
+				   // and cell vector according to the
+				   // DG method given in the
+				   // introduction.
   for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
     {
       double beta_n=beta[point] * normals[point];
       if (beta_n>0)
+					 // The term $(\beta\cdot n
+					 // u,v)_{\partial K_+}$,
 	for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
 	  for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
-	    cell_matrix(i,j) += beta_n *
-				v(j,point) *
-				v(i,point) *
-				JxW[point];
+	    u_v_matrix(i,j) += beta_n *
+			       v(j,point) *
+			       v(i,point) *
+			       JxW[point];
       else
 	{
+					   // at the boundary the term
+					   // $(\beta\cdot n
+					   // g,v)_{\partial
+					   // K_-\cap\partial\Omega}$,
 	  if (face->at_boundary())
 	    for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
 	      cell_vector(i) -= beta_n *
@@ -253,26 +364,200 @@ void DGAssembler<dim>::assemble_face_term(const FEFaceValuesBase<dim>& fe_v,
 				v(i,point) *
 				JxW[point];
 	  else
+					     // and on inner faces the
+					     // term $(\beta\cdot n
+					     // \hat u,v)_{\partial
+					     // K_-}$
 	    for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
 	      for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
-		cell_inflow_matrix(i,k) += beta_n *
-					   v_neighbor(k,point) *
-					   v(i,point) *
-					   JxW[point];
+		un_v_matrix(i,k) += beta_n *
+				    v_neighbor(k,point) *
+				    v(i,point) *
+				    JxW[point];
 	}
     }
 }
 
+				 // Now we look at the assembling
+				 // function that assembles the face
+				 // terms corresponding to the second
+				 // version of the DG method,
+				 // cf. above. Then, the face terms
+				 // are given as a sum of integrals
+				 // over all faces.  Here we need two
+				 // additional cell matrices
+				 // ``u_vn_matrix'' and
+				 // ``un_vn_matrix'' that will store
+				 // contributions due to terms
+				 // involving u and vn as well as un
+				 // and vn.
+template <int dim>
+void DGTransportEquation<dim>::assemble_face_term2(
+  const FEFaceValuesBase<dim>& fe_v,
+  const FEFaceValuesBase<dim>& fe_v_neighbor,      
+  FullMatrix<double> &u_v_matrix,
+  FullMatrix<double> &un_v_matrix,
+  FullMatrix<double> &u_vn_matrix,
+  FullMatrix<double> &un_vn_matrix,
+  Vector<double> &cell_vector)
+{
+				   // the first few lines are the same
+  const FullMatrix<double> &v = fe_v.get_shape_values ();
+  const FullMatrix<double> &v_neighbor = fe_v_neighbor.get_shape_values ();  
+  const vector<double> &JxW = fe_v.get_JxW_values ();
+  const vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+
+  vector<Point<dim> > beta (fe_v.n_quadrature_points);
   
+  beta_function.value_list (fe_v.get_quadrature_points(), beta);
+
+  vector<double> g(fe_v.n_quadrature_points);
+  DoFHandler<dim>::face_iterator face=fe_v.get_face();
+  if (face->at_boundary())
+    boundary_function.value_list (fe_v.get_quadrature_points(), g);
+
+  for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
+    {
+      double beta_n=beta[point] * normals[point];
+      if (beta_n>0)
+	{
+					   // This terms we've already seen,
+	  for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
+	    for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
+	      u_v_matrix(i,j) += beta_n *
+				 v(j,point) *
+				 v(i,point) *
+				 JxW[point];
+
+					   // on inner faces we
+					   // additionally have the
+					   // term $(\beta\cdot n
+					   // u,\hat v)_{\partial K_+},
+	  if (!face->at_boundary())
+	    for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
+	      for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
+		u_vn_matrix(k,j) -= beta_n *
+				    v(j,point) *
+				    v_neighbor(k,point) *
+				    JxW[point];
+	}
+      else
+	{
+					   // this one we already know,
+	  if (face->at_boundary())
+	    for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
+	      cell_vector(i) -= beta_n *
+				g[point] *
+				v(i,point) *
+				JxW[point];
+	  else
+	    {
+					       // this one also,
+	      for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
+		for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l)
+		  un_v_matrix(i,l) += beta_n *
+				      v_neighbor(l,point) *
+				      v(i,point) *
+				      JxW[point];
+
+					       // and this is another
+					       // new one:
+					       // $(\beta\cdot n \hat
+					       // u,\hat v)_{\partial
+					       // K_-}$.
+	      for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
+		for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l)
+		  un_vn_matrix(k,l) -= beta_n *
+				       v_neighbor(l,point) *
+				       v_neighbor(k,point) *
+				       JxW[point];
+	    }
+	}
+    }
+}
+
+
+				 // After these preparations, we
+				 // proceed with the main part of this
+				 // program. The main class, here
+				 // called ``DGMethod'' is basically
+				 // the main class of step 6. One of
+				 // the differences is that there's no
+				 // ConstraintMatrix object. This is,
+				 // because there are no hanging nodes
+				 // in DG discretizations.
+template <int dim>
+class DGMethod
+{
+  public:
+    DGMethod ();
+    ~DGMethod ();
+
+    void run ();
+    
+  private:
+    void setup_system ();
+    void assemble_system1 ();
+    void assemble_system2 ();
+    void solve (Vector<double> &solution);
+    void refine_grid ();
+    void output_results (const unsigned int cycle) const;
+    
+    Triangulation<dim>   triangulation;
+    MappingQ1<dim>       mapping;
+    
+				     // Furthermore we want to
+				     // use DG elements of degree 1
+				     // (but this is only specified in
+				     // the constructor):
+    FE_DGQ<dim>          fe;
+    DoFHandler<dim>      dof_handler;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+				     // We define the quadrature
+				     // formulae for the cell and the
+				     // face terms of the
+				     // discretization.
+    QGauss4<dim>   quadrature;
+    QGauss4<dim-1> face_quadrature;
+    
+				     // And there are two solution
+				     // vectors, that store the
+				     // solutions to the problems
+				     // corresponding to the two
+				     // different assembling routines
+				     // ``assemble_system1'' and
+				     // ``assemble_system2'';
+    Vector<double>       solution1;
+    Vector<double>       solution2;
+    Vector<double>       right_hand_side;
+    
+				     // Finally this class includes an
+				     // object of the
+				     // DGTransportEquations class
+				     // described above.
+    DGTransportEquation<dim>     dg;
+};
+
+
+
+				 // Now for the implementation of the
+				 // main class. Constructor and
+				 // destructor follow the same
+				 // pattern that was used previously,
+				 // so we need not comment on these
+				 // two functions:  
 template <int dim>
-TransportProblem<dim>::TransportProblem () :
+DGMethod<dim>::DGMethod () :
                 fe (1),
 		dof_handler (triangulation)
 {}
 
 
 template <int dim>
-TransportProblem<dim>::~TransportProblem () 
+DGMethod<dim>::~DGMethod () 
 {
   dof_handler.clear ();
 };
@@ -280,172 +565,460 @@ TransportProblem<dim>::~TransportProblem ()
 
 
 template <int dim>
-void TransportProblem<dim>::setup_system ()
+void DGMethod<dim>::setup_system ()
 {
-				   // To distribute degrees of
-				   // freedom, the ``dof_handler''
-				   // variable takes only the finite
-				   // element object. In this case, it
-				   // will distribute four degrees of
-				   // freedom per cell.
+				   // First we need to distribute the
+				   // DoFs.
   dof_handler.distribute_dofs (fe);
 
+				   // The DoFs of a cell are coupled
+				   // with all DoFs of all neighboring
+				   // cells.  Therefore the maximum
+				   // number of matrix entries is
+				   // needed when all neighbors of a
+				   // cell are once more refined than
+				   // the cell under consideration.
   sparsity_pattern.reinit (dof_handler.n_dofs(),
 			   dof_handler.n_dofs(),
-			   dof_handler.max_couplings_between_dofs());
+			   (GeometryInfo<dim>::faces_per_cell
+			    *GeometryInfo<dim>::subfaces_per_face+1)*fe.dofs_per_cell);
+  
+				   // For DG discretizations we call
+				   // the function analogue to
+				   // DoFTools::make_sparsity_pattern.
   DoFTools::make_flux_sparsity_pattern (dof_handler, sparsity_pattern);
+  
+				   // All following function calls are
+				   // already known.
   sparsity_pattern.compress();
   
   system_matrix.reinit (sparsity_pattern);
 
-  solution.reinit (dof_handler.n_dofs());
+  solution1.reinit (dof_handler.n_dofs());
+  solution2.reinit (dof_handler.n_dofs());
   right_hand_side.reinit (dof_handler.n_dofs());
 };
 
 
-
+				 // We proceed with the
+				 // ``assemble_system1'' function that
+				 // implements the DG discretization
+				 // in its first version. This
+				 // function repeatedly calls the
+				 // ``assemble_cell_term'' and
+				 // ``assemble_face_term1'' functions
+				 // of the DGTransportEquation object.
+				 // The ``assemble_face_term1''
+				 // function takes two
+				 // FEFaceValuesBase objects; one for
+				 // the shape functions on the current
+				 // cell and the other for shape
+				 // functions on the neighboring cell
+				 // under consideration. Both objects
+				 // are either of class FEFaceValues
+				 // or of class FESubfaceValues (both
+				 // derived from FEFaceValuesBase)
+				 // according to following cases
+				 // already mentioned in the
+				 // introduction:
+				 //
+				 // 1. face is at boundary (current
+				 // cell: FEFaceValues, neighboring
+				 // cell does not exist);
+				 //
+				 // 2. neighboring cell is finer
+				 // (current cell: FESubfaceValues,
+				 // neighboring cell: FEFaceValues);
+				 //
+				 // 3. neighboring cell is of the same
+				 // refinement level (both, current
+				 // and neighboring cell:
+				 // FEFaceValues);
+				 //
+				 // 4. neighboring cell is coarser
+				 // (current cell: FEFaceValues,
+				 // neighboring cell:
+				 // FESubfaceValues).
+				 //
+				 // If we considered globally refined
+				 // meshes then only cases 1 and 3
+				 // would occur. But as we consider
+				 // also locally refined meshes we
+				 // need to distinguish all four cases
+				 // making the following assembling
+				 // function a bit longish.
 template <int dim>
-void TransportProblem<dim>::assemble_system () 
+void DGMethod<dim>::assemble_system1 () 
 {
-				   // See Cockburn paper for the proper quadrature.
-  QGauss4<dim>  quadrature;
-  QGauss4<dim-1>  face_quadrature;
-  
   const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
   vector<unsigned int> dofs (dofs_per_cell);
   vector<unsigned int> dofs_neighbor (dofs_per_cell);
 
+				   // First we create the Update flags
+				   // for the FEValues and the
+				   // FEFaceValues objects.
   UpdateFlags update_flags = UpdateFlags(update_values
 					 | update_gradients
 					 | update_q_points
 					 | update_JxW_values);
   
+				   // Note, that on faces we do not
+				   // need gradients but we need
+				   // normal vectors.
   UpdateFlags face_update_flags = UpdateFlags(update_values
 					      | update_q_points
 					      | update_JxW_values
 					      | update_normal_vectors);
+  
+				   // On the neighboring cell we only
+				   // need the shape values. Given a
+				   // specific face, the quadrature
+				   // points and `JxW values' are the
+				   // same as for the current cells,
+				   // the normal vectors are known to
+				   // be the negative of the normal
+				   // vectors of the current cell.
+  UpdateFlags neighbor_face_update_flags = UpdateFlags(update_values);
    
-
+				   // Then we create the FEValues
+				   // object. Note, that since version
+				   // 3.2.0 the constructor of this
+				   // class takes a Mapping object as
+				   // first argument. Although the
+				   // constructor without Mapping
+				   // argument is still supported it
+				   // is recommended to use the new
+				   // constructor. This reduces the
+				   // effect of `hidden magic' (the
+				   // old constructor implicitely
+				   // assumes a MappingQ1 mapping) and
+				   // makes it easier to change the
+				   // Mapping object later.
   FEValues<dim> fe_v (
     mapping, fe, quadrature, update_flags);
+  
+				   // Similarly we create the
+				   // FEFaceValues and FESubfaceValues
+				   // objects for both, the current
+				   // and the neighboring cell. Within
+				   // the following nested loop over
+				   // all cells and all faces of the
+				   // cell they will be reinited to
+				   // the current cell and the face
+				   // (and subface) number.
   FEFaceValues<dim> fe_v_face (
     mapping, fe, face_quadrature, face_update_flags);
   FESubfaceValues<dim> fe_v_subface (
     mapping, fe, face_quadrature, face_update_flags);
   FEFaceValues<dim> fe_v_face_neighbor (
-    mapping, fe, face_quadrature, UpdateFlags(update_values | update_default));
+    mapping, fe, face_quadrature, neighbor_face_update_flags);
   FESubfaceValues<dim> fe_v_subface_neighbor (
-    mapping, fe, face_quadrature, UpdateFlags(update_values | update_default));
-
-				   // includes the u and v terms
-  FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
-				   // includes u_hat and v terms
-  FullMatrix<double> cell_inflow_matrix (dofs_per_cell, dofs_per_cell);
+    mapping, fe, face_quadrature, neighbor_face_update_flags);
+
+				   // Now we create the cell matrices
+				   // and vectors. Here we need two
+				   // cell matrices, both for face
+				   // terms that include test
+				   // functions ``v'' (shape functions
+				   // of the current cell). To be more
+				   // precise, the first matrix will
+				   // include the `u and v terms' and
+				   // the second that will include the
+				   // `un and v terms'. Here we recall
+				   // our the convention that `un' is
+				   // the short cut for `u_neighbor'
+				   // and represents the $u_hat$, see
+				   // introduction.
+  FullMatrix<double> u_v_matrix (dofs_per_cell, dofs_per_cell);
+  FullMatrix<double> un_v_matrix (dofs_per_cell, dofs_per_cell);
 
   Vector<double>  cell_vector (dofs_per_cell);
 
+				   // Furthermore we need some cell
+				   // and face iterators
   DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
 					endc = dof_handler.end();
   DoFHandler<dim>::face_iterator face;
   DoFHandler<dim>::cell_iterator neighbor;
-  DoFHandler<dim>::cell_iterator neighbor_child;
+  DoFHandler<dim>::active_cell_iterator neighbor_child;
 
+				   // Now we start the loop over all
+				   // active cells
   for (;cell!=endc; ++cell) 
     {
-				       // re-init fe values for this cell
+				       // and reinit the FEValues
+				       // object for the current cell,
       fe_v.reinit (cell);
 
-      cell_matrix.clear ();
-      cell_vector.clear ();
-
-      dg_assembler.assemble_cell_term(fe_v,
-				      cell_matrix,
-				      cell_vector);
-      
+				       // Call the function that
+				       // assembles the cell
+				       // terms. The first argument is
+				       // the FEValues that was
+				       // already reinited on current
+				       // the cell.
+      dg.assemble_cell_term(fe_v,
+			    u_v_matrix,
+			    cell_vector);
+
+				       // As in previous example steps
+				       // the vector `dofs' includes
+				       // the dof_indices.
       cell->get_dof_indices (dofs);
 
+				       // This is the start of the
+				       // nested loop over all faces.
       for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
 	{
+					   // First we set the face
+					   // iterator.
 	  face = cell->face(face_no);
-	  
-	  cell_inflow_matrix.clear();
 
+					   // Now we distinguish the
+					   // four different cases in
+					   // the ordering mentioned
+					   // above. We start with
+					   // faces belonging to the
+					   // boundary of the domain.
 	  if (face->at_boundary())
 	    {
+					       // We reinit the
+					       // FEFaceValues object
+					       // to the current face
 	      fe_v_face.reinit (cell, face_no);
 
-	      dg_assembler.assemble_face_term(fe_v_face,
-					      fe_v_face,
-					      cell_matrix,
-					      cell_inflow_matrix,
-					      cell_vector);
+					       // and assemble the
+					       // corresponding face
+					       // terms. Here, the
+					       // second and fourth
+					       // arguments are only
+					       // dummy arguments. On
+					       // the boundary of the
+					       // domain the
+					       // ``assemble_face_term1''
+					       // function will not
+					       // access to shape
+					       // values on the
+					       // non-existent
+					       // neighboring
+					       // cell. Also,
+					       // ``un_v_matrix'' will
+					       // be unchanged.
+	      dg.assemble_face_term1(fe_v_face,
+				     fe_v_face,
+				     u_v_matrix,
+				     un_v_matrix,
+				     cell_vector);
 	    }
-	  else // if (!face->at_boundary())
+	  else
 	    {
-	      Assert (cell->neighbor(face_no).state() == valid, ExcInternalError());
+					       // When we are not on the
+					       // boundary of the
+					       // domain then there
+					       // must exist a
+					       // neighboring cell.
 	      neighbor = cell->neighbor(face_no);
-	      
-	      if (face->has_children())  // i.e. neighbor is one level more refined than cell
+
+					       // We proceed with the
+					       // second and most
+					       // complicated case:
+					       // the neighboring cell
+					       // is more refined than
+					       // the current cell. As
+					       // in deal.II
+					       // neighboring cells
+					       // are restricted to
+					       // have a level
+					       // difference of not
+					       // more than one, the
+					       // neighboring cell is
+					       // known to be only
+					       // ONCE more refined
+					       // than the current
+					       // cell. Furthermore
+					       // also the face is
+					       // once more refined,
+					       // i.e. it has
+					       // children.
+	      if (face->has_children())
 		{
-						   // store which number #cell# has in the
-						   // list of neighbors of #neighbor#
-		  const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
+						   // first we store
+						   // which number the
+						   // current cell has
+						   // in the list of
+						   // neighbors of the
+						   // neighboring
+						   // cell. Hence,
+						   // neighbor->neighbor(neighbor2)
+						   // equals the
+						   // current cell
+						   // ``cell''.
+		  const unsigned int neighbor2=
+		    cell->neighbor_of_neighbor(face_no);
 		  
 		  
-						   // loop over all subfaces
-		  for (unsigned int subface_no=0; subface_no<GeometryInfo<dim>::subfaces_per_face;
+						   // We loop over
+						   // subfaces
+		  for (unsigned int subface_no=0;
+		       subface_no<GeometryInfo<dim>::subfaces_per_face;
 		       ++subface_no)
 		    {
-						       // get an iterator pointing to the
-						       // cell behind the present subface
+						       // and set the
+						       // cell
+						       // iterator
+						       // ``neighbor_child''
+						       // to the cell
+						       // placed
+						       // `behind' the
+						       // current
+						       // subface.
 		      neighbor_child = neighbor->child(GeometryInfo<dim>::
 						       child_cell_on_face(neighbor2,subface_no));
+
+						       // As these are
+						       // quite
+						       // complicated
+						       // indirections
+						       // we check for
+						       // the internal
+						       // consistency.
 		      Assert (neighbor_child->face(neighbor2) == face->child(subface_no),
 			      ExcInternalError());
 		      Assert (!neighbor_child->has_children(), ExcInternalError());
 
+						       // As already
+						       // mentioned
+						       // above for
+						       // this case
+						       // (case 2) we
+						       // employ the
+						       // FESubfaceValues
+						       // of the
+						       // current
+						       // cell, here
+						       // reinited for
+						       // the current
+						       // cell, face
+						       // and subface,
+						       // and we
+						       // employ the
+						       // FEFaceValues
+						       // of the
+						       // neighboring
+						       // child cell.
 		      fe_v_subface.reinit (cell, face_no, subface_no);
 		      fe_v_face_neighbor.reinit (neighbor_child, neighbor2);
-	
-		      cell_inflow_matrix.clear();
-	    
-		      dg_assembler.assemble_face_term(fe_v_subface,
-						      fe_v_face_neighbor,
-						      cell_matrix,
-						      cell_inflow_matrix,
-						      cell_vector);
-		  
-						       // get indices of dofs of neighbor_child cell
+
+		      dg.assemble_face_term1(fe_v_subface,
+					     fe_v_face_neighbor,
+					     u_v_matrix,
+					     un_v_matrix,
+					     cell_vector);
+		      
+						       // get dof
+						       // indices of
+						       // the
+						       // neighbor_child
+						       // cell
 		      neighbor_child->get_dof_indices (dofs_neighbor);
 		      						
-						       // distribute cell matrix
+						       // distribute
+						       // cell matrix
+						       // to the
+						       // system_matrix
 		      for (unsigned int i=0; i<dofs_per_cell; ++i)
 			for (unsigned int k=0; k<dofs_per_cell; ++k)
 			  system_matrix.add(dofs[i], dofs_neighbor[k],
-					    cell_inflow_matrix(i,k));
+					    un_v_matrix(i,k));
+
+		      				       // In the
+						       // ``assemble_face_term1''
+						       // function contributions to
+						       // the cell matrices and the
+						       // cell vector are only
+						       // ADDED. Therefore on each
+						       // subface we need to reset the
+						       // un_v_matrix
+						       // to zero, before assembling
+						       // the face terms corresponding
+						       // to the following neighbor_child cell.
+		      un_v_matrix.clear();
 		    }
 		}
-	      else // if (!face->has_children())
+					       // End of ``if
+					       // (face->has_children())''
+	      else
 		{
+						   // We proceed with
+						   // case 3,
+						   // i.e. neighboring
+						   // cell is of the
+						   // same refinement
+						   // level as the
+						   // current cell.
 		  if (neighbor->level() == cell->level()) 
 		    {
-						       // store which number #cell# has in the
-						       // list of neighbors of #neighbor#
+						   // Like before we
+						   // store which
+						   // number the
+						   // current cell has
+						   // in the list of
+						   // neighbors of the
+						   // neighboring
+						   // cell.
 		      const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
 
+						       // We reinit
+						       // the
+						       // FEFaceValues
+						       // of the
+						       // current and
+						       // neighboring
+						       // cell to the
+						       // current face
+						       // and assemble
+						       // the
+						       // corresponding
+						       // face terms.
 		      fe_v_face.reinit (cell, face_no);
 		      fe_v_face_neighbor.reinit (neighbor, neighbor2);
 		      
-		      dg_assembler.assemble_face_term(fe_v_face,
-						      fe_v_face_neighbor,
-						      cell_matrix,
-						      cell_inflow_matrix,
-						      cell_vector);
+		      dg.assemble_face_term1(fe_v_face,
+					     fe_v_face_neighbor,
+					     u_v_matrix,
+					     un_v_matrix,
+					     cell_vector);
+						       // End of ``if
+						       // (neighbor->level()
+						       // ==
+						       // cell->level())''
 		    }
-		  else // if (neighbor->level() < cell->level()) i.e. neighbor is one level coarser than cell
+		  else
 		    {
+						       // Finally we
+						       // consider
+						       // case 4. When
+						       // the
+						       // neighboring
+						       // cell is not
+						       // finer and
+						       // not of the
+						       // same
+						       // refinement
+						       // level as the
+						       // current cell
+						       // it must be
+						       // coarser.
 		      Assert(neighbor->level() < cell->level(), ExcInternalError());
 
+						       // Find out the
+						       // how many'th
+						       // face_no and
+						       // subface_no
+						       // the current
+						       // face is
+						       // w.r.t. the
+						       // neighboring
+						       // cell.
 		      const std::pair<unsigned int, unsigned int> faceno_subfaceno=
 			cell->neighbor_of_coarser_neighbor(face_no);
 		      const unsigned int neighbor_face_no=faceno_subfaceno.first,
@@ -454,127 +1027,460 @@ void TransportProblem<dim>::assemble_system ()
 		      Assert (neighbor->neighbor(neighbor_face_no)
 			      ->child(GeometryInfo<dim>::child_cell_on_face(
 				face_no,neighbor_subface_no)) == cell, ExcInternalError());
-			
-						       // now 'neighbor_face_no' stores the number
-						       // of a face in the list of faces of 'neighbor'.
-						       // This face has got a subface that is 
-						       // between 'cell' and 'neighbor'.
-						       // 'neighbor_subface_no' stores the number
-						       // of this subface in the list of subfaces of this
-						       // face 'neighbor->face(neighbor_face_no)'
-						       // that is between 'cell' and 'neighbor'
+
+						       // Reinit the
+						       // appropriate
+						       // FEFaceValues
+						       // and assemble
+						       // the face
+						       // terms.
 		      fe_v_face.reinit (cell, face_no);
 		      fe_v_subface_neighbor.reinit (neighbor, neighbor_face_no,
 						    neighbor_subface_no);
 		      
-		      dg_assembler.assemble_face_term(fe_v_face,
-						      fe_v_subface_neighbor,
-						      cell_matrix,
-						      cell_inflow_matrix,
-						      cell_vector);
-		    } // else // if (neighbor->level() < cell->level())
-
-						   // get indices of dofs of neighbor_child cell
+		      dg.assemble_face_term1(fe_v_face,
+					     fe_v_subface_neighbor,
+					     u_v_matrix,
+					     un_v_matrix,
+					     cell_vector);
+		    }
+
+						   // Get dof indices
+						   // of the
+						   // neighbor_child
+						   // cell,
 		  neighbor->get_dof_indices (dofs_neighbor);
 		 						
-						   // distribute cell_inflow_matrix
+						   // distribute the
+						   // un_v_matrix,
 		  for (unsigned int i=0; i<dofs_per_cell; ++i)
 		    for (unsigned int k=0; k<dofs_per_cell; ++k)
 		      system_matrix.add(dofs[i], dofs_neighbor[k],
-					cell_inflow_matrix(i,k));
-		} // else // if (!face->has_children())
-	    }  // else // if (!face->at_boundary())
-	} //for (face_no...)
+					un_v_matrix(i,k));
+		  
+						   // and clear the
+						   // ``un_v_matrix''
+						   // on each face.
+		  un_v_matrix.clear();
+		}
+					       // End of ``face not at boundary'':
+	    }
+					   // End of loop over all faces:
+	}
       
-				       // distribute cell matrix
+				       // Finally we distribute the
+				       // u_v_matrix,
       for (unsigned int i=0; i<dofs_per_cell; ++i)
 	for (unsigned int j=0; j<dofs_per_cell; ++j)
-	  system_matrix.add(dofs[i], dofs[j], cell_matrix(i,j));
+	  system_matrix.add(dofs[i], dofs[j], u_v_matrix(i,j));
       
-				       // distribute cell vector
+				       // the cell vector
       for (unsigned int i=0; i<dofs_per_cell; ++i)
 	right_hand_side(dofs[i]) += cell_vector(i);
-    }  // for (cell...)
+
+				       // and clear them both.
+      u_v_matrix.clear ();
+      cell_vector.clear ();
+    }
 };
 
 
 
+				 // We proceed with the
+				 // ``assemble_system2'' function that
+				 // implements the DG discretization
+				 // in its second version. This
+				 // function is very similar to the
+				 // ``assemble_system1''
+				 // function. Therefore, here we only
+				 // discuss the differences between
+				 // the two functions. This function
+				 // repeatedly calls the
+				 // ``assemble_face_term2'' function
+				 // of the DGTransportEquation object,
+				 // that assembles the face terms
+				 // written as a sum of integrals over
+				 // all faces. Therefore, we need to
+				 // make sure that each face is
+				 // treated only once. This is achieved
+				 // by introducing the rule:
+				 // 
+				 // a) If the current and the
+				 // neighboring cells are of the same
+				 // refinement level we access and
+				 // treat the face from the cell with
+				 // lower index.
+				 //
+				 // b) If the two cells are of
+				 // different refinement levels we
+				 // access and treat the face from the
+				 // coarser cell.
+				 //
+				 // Due to rule b) we do not need to
+				 // consider case 4 (neighboring cell
+				 // is coarser) any more.
+
 template <int dim>
-void TransportProblem<dim>::solve () 
-{  
-  SolverControl           solver_control (1000, 1e-12);
+void DGMethod<dim>::assemble_system2 () 
+{
+  const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+  vector<unsigned int> dofs (dofs_per_cell);
+  vector<unsigned int> dofs_neighbor (dofs_per_cell);
+
+  UpdateFlags update_flags = UpdateFlags(update_values
+					 | update_gradients
+					 | update_q_points
+					 | update_JxW_values);
+  
+  UpdateFlags face_update_flags = UpdateFlags(update_values
+					      | update_q_points
+					      | update_JxW_values
+					      | update_normal_vectors);
+   
+  UpdateFlags neighbor_face_update_flags = UpdateFlags(update_values);
+
+				   // Here we do not need
+				   // ``fe_v_face_neighbor'' as case 4
+				   // does not occur.
+  FEValues<dim> fe_v (
+    mapping, fe, quadrature, update_flags);
+  FEFaceValues<dim> fe_v_face (
+    mapping, fe, face_quadrature, face_update_flags);
+  FESubfaceValues<dim> fe_v_subface (
+    mapping, fe, face_quadrature, face_update_flags);
+  FEFaceValues<dim> fe_v_face_neighbor (
+    mapping, fe, face_quadrature, neighbor_face_update_flags);
+
+
+  FullMatrix<double> u_v_matrix (dofs_per_cell, dofs_per_cell);
+  FullMatrix<double> un_v_matrix (dofs_per_cell, dofs_per_cell);
+  
+				   // Additionally we need following
+				   // two cell matrices, both for face
+				   // term that include test function
+				   // ``vn'' (shape functions of the
+				   // neighboring cell). To be more
+				   // precise, the first matrix will
+				   // include the `u and vn terms' and
+				   // the second that will include the
+				   // `un and vn terms'.
+  FullMatrix<double> u_vn_matrix (dofs_per_cell, dofs_per_cell);
+  FullMatrix<double> un_vn_matrix (dofs_per_cell, dofs_per_cell);
+  
+  Vector<double>  cell_vector (dofs_per_cell);
+
+				   // Furthermore, here we define a
+				   // dummy matrix and rhs to
+				   // emphasize when arguments of the
+				   // ``assemble_face_term2''
+				   // functions will not be access.
+  FullMatrix<double> dummy_matrix;
+  Vector<double>     dummy_rhs;
+
+				   // The following lines are roughly
+				   // the same as in the previous
+				   // function.
+  DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+					endc = dof_handler.end();
+  DoFHandler<dim>::face_iterator face;
+  DoFHandler<dim>::cell_iterator neighbor;
+  DoFHandler<dim>::cell_iterator neighbor_child;
+
+  for (;cell!=endc; ++cell) 
+    {
+      fe_v.reinit (cell);
+
+      dg.assemble_cell_term(fe_v,
+			    u_v_matrix,
+			    cell_vector);
+      
+      cell->get_dof_indices (dofs);
+
+      for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+	{
+	  face = cell->face(face_no);
+
+					   // Case 1:
+	  if (face->at_boundary())
+	    {
+	      fe_v_face.reinit (cell, face_no);
+
+	      dg.assemble_face_term2(fe_v_face,
+				     fe_v_face,
+				     u_v_matrix,
+				     dummy_matrix,
+				     dummy_matrix,
+				     dummy_matrix,
+				     cell_vector);
+	    }
+	  else
+	    {
+	      Assert (cell->neighbor(face_no).state() == valid, ExcInternalError());
+	      neighbor = cell->neighbor(face_no);
+
+					       // Case 2:
+	      if (face->has_children())
+		{
+		  const unsigned int neighbor2=
+		    cell->neighbor_of_neighbor(face_no);
+		  
+		  for (unsigned int subface_no=0;
+		       subface_no<GeometryInfo<dim>::subfaces_per_face;
+		       ++subface_no)
+		    {
+		      neighbor_child = neighbor->child(
+			GeometryInfo<dim>::child_cell_on_face(neighbor2,subface_no));
+		      Assert (neighbor_child->face(neighbor2) == face->child(subface_no),
+			      ExcInternalError());
+		      Assert (!neighbor_child->has_children(), ExcInternalError());
+
+		      fe_v_subface.reinit (cell, face_no, subface_no);
+		      fe_v_face_neighbor.reinit (neighbor_child, neighbor2);
+
+		      dg.assemble_face_term2(fe_v_subface,
+					     fe_v_face_neighbor,
+					     u_v_matrix,
+					     un_v_matrix,
+					     u_vn_matrix,
+					     un_vn_matrix,
+					     dummy_rhs);
+		  
+		      neighbor_child->get_dof_indices (dofs_neighbor);
+		      						
+		      for (unsigned int i=0; i<dofs_per_cell; ++i)
+			for (unsigned int j=0; j<dofs_per_cell; ++j)
+			  {
+			    system_matrix.add(dofs[i], dofs_neighbor[j],
+					      un_v_matrix(i,j));
+			    system_matrix.add(dofs_neighbor[i], dofs[j],
+					      u_vn_matrix(i,j));
+			    system_matrix.add(dofs_neighbor[i], dofs_neighbor[j],
+					      un_vn_matrix(i,j));
+			  }
+		      
+		      un_v_matrix.clear();
+		      u_vn_matrix.clear();
+		      un_vn_matrix.clear();
+		    }
+		}
+	      else
+		{
+						   // Case 3, with the
+						   // additional rule
+						   // a)
+		  if (neighbor->level() == cell->level() &&
+		      neighbor->index() > cell->index()) 
+		    {
+		      const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
+
+		      fe_v_face.reinit (cell, face_no);
+		      fe_v_face_neighbor.reinit (neighbor, neighbor2);
+		      
+		      dg.assemble_face_term2(fe_v_face,
+					     fe_v_face_neighbor,
+					     u_v_matrix,
+					     un_v_matrix,
+					     u_vn_matrix,
+					     un_vn_matrix,
+					     dummy_rhs);
+
+		      neighbor->get_dof_indices (dofs_neighbor);
+
+		      for (unsigned int i=0; i<dofs_per_cell; ++i)
+			for (unsigned int j=0; j<dofs_per_cell; ++j)
+			  {
+			    system_matrix.add(dofs[i], dofs_neighbor[j],
+					      un_v_matrix(i,j));
+			    system_matrix.add(dofs_neighbor[i], dofs[j],
+					      u_vn_matrix(i,j));
+			    system_matrix.add(dofs_neighbor[i], dofs_neighbor[j],
+					      un_vn_matrix(i,j));
+			  }
+		      
+		      un_v_matrix.clear();
+		      u_vn_matrix.clear();
+		      un_vn_matrix.clear();
+		    }
+
+						   // Due to rule b)
+						   // we do not need
+						   // to consider case
+						   // 4.
+		}
+	    }
+	}
+      
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+	for (unsigned int j=0; j<dofs_per_cell; ++j)
+	  system_matrix.add(dofs[i], dofs[j], u_v_matrix(i,j));
+      
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+	right_hand_side(dofs[i]) += cell_vector(i);
+      
+      u_v_matrix.clear ();
+      cell_vector.clear ();
+    }
+};
+
+				 // For this simple solver we use the
+				 // simplest possible solver, called
+				 // richardson iteration, that
+				 // represents a simple defect
+				 // correction. This, in combination
+				 // with a block SSOR preconditioner,
+				 // that uses the special block matrix
+				 // structur of system matrices
+				 // arising from DG
+				 // discretizations. The size of these
+				 // blocks are the number of DoFs
+				 // per cell. Here, we use a SSOR
+				 // preconditioning as we have not
+				 // renumbered the DoFs according to
+				 // the flow field. If the DoFs are
+				 // renumbered downstream the flow,
+				 // then a block Gauss-Seidel
+				 // preconditioner (see the
+				 // PreconditionBlockSOR class with
+				 // relaxation=1) makes a much better
+				 // job.
+template <int dim>
+void DGMethod<dim>::solve (Vector<double> &solution) 
+{
+  SolverControl           solver_control (1000, 1e-12, false, false);
   PrimitiveVectorMemory<> vector_memory;
   SolverRichardson<>      solver (solver_control, vector_memory);
 
+				   // Here we create the
+				   // preconditioner,
   PreconditionBlockSSOR<double> preconditioner;
+
+				   // we asigned the matrix to it and
+				   // set the right block size.
   preconditioner.initialize(system_matrix, fe.dofs_per_cell);
+
+				   // As the inverses of the diagonal
+				   // blocks are needed in each
+				   // preconditioner step, it is wise
+				   // to invert the diagonal blocks of
+				   // the matrix before starting the
+				   // solver. Otherwise, the diagonal
+				   // blocks are inverted in each
+				   // preconditioner step,
+				   // significantly slowing down the
+				   // linear solving process.
   preconditioner.invert_diagblocks();
-  
+
+				   // After these preparations we are
+				   // ready to start the linear solver.
   solver.solve (system_matrix, solution, right_hand_side,
 		preconditioner);
 };
 
 
+				 // We refine the grid according to a
+				 // very simple refinement criterion,
+				 // namely the gradients of the
+				 // solution. As here we consider the
+				 // DG(1) method (i.e. we use
+				 // piecewise bilinear shape
+				 // functions) we could simply compute
+				 // the gradients on each cell. But we
+				 // do not want to base our refinement
+				 // indicator on the gradients on each
+				 // cell only, but want to base them
+				 // also on jumps of the discontinuous
+				 // solution function over faces
+				 // between neighboring cells. The
+				 // simpliest way of doing that is to
+				 // compute approximative gradients by
+				 // difference quotients including the
+				 // cell under consideration and its
+				 // neighbors. This is done by the
+				 // DerivativeApproximation class that
+				 // computes the approximate
+				 // gradients in a way similar to the
+				 // GradientEstimation described in
+				 // Step 9 of this tutorial. According
+				 // to the argumentation in Step 9,
+				 // here we consider
+				 // $h^{1+d/2}|\nabla_h
+				 // u_h|$. Futhermore we note that we
+				 // do not consider approximate
+				 // second derivatives because
+				 // solutions to the linear advection
+				 // equation are in general not in H^2
+				 // but in H^1 (to be more precise, in
+				 // H^1_\beta) only.
 template <int dim>
-void TransportProblem<dim>::refine_grid ()
+void DGMethod<dim>::refine_grid ()
 {
-  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
-
-  FunctionMap<dim>::type neumann_boundary;
-
-  KellyErrorEstimator<dim>::estimate (dof_handler,
-				      QGauss3<dim-1>(),
-				      neumann_boundary,
-				      solution,
-				      estimated_error_per_cell);
+				   // The DerivativeApproximation
+				   // class computes the gradients to
+				   // float precision. This is
+				   // sufficient as they are
+				   // approximate and serve as
+				   // refinement indicators only.
+  Vector<float> gradient_indicator (triangulation.n_active_cells());
+
+				   // Now the approximate gradients
+				   // are computed
+  DerivativeApproximation::approximate_gradient (mapping,
+						 dof_handler,
+						 solution2,
+						 gradient_indicator);
+
+				   // and they are cell-wise scaled by
+				   // the factor $h^{1+d/2}$
+  DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+					endc = dof_handler.end();
+  for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no)
+    gradient_indicator(cell_no)*=std::pow(cell->diameter(), 1+1.0*dim/2);
 
+				   // Finally they serve as refinement
+				   // indicator.
   GridRefinement::refine_and_coarsen_fixed_number (triangulation,
-						   estimated_error_per_cell,
-						   0.3, 0.03);
+						   gradient_indicator,
+						   0.3, 0.1);
 
   triangulation.execute_coarsening_and_refinement ();
-};
-
+}
 
 
+				 // The output of this program
+				 // consists of eps-files of the
+				 // adaptively refined grids and the
+				 // numerical solutions given in
+				 // gnuplot format. This was covered
+				 // in previous examples and will not
+				 // be further commented on.
 template <int dim>
-void TransportProblem<dim>::output_results (const unsigned int cycle) const
+void DGMethod<dim>::output_results (const unsigned int cycle) const
 {
-				   // We want to write the grid in
-				   // each cycle. Here is another way
-				   // to quickly produce a filename
-				   // based on the cycle number. It
-				   // assumes that the numbers `0'
-				   // through `9' are represented
-				   // consecutively in the character
-				   // set (which is the case in all
-				   // known character sets). However,
-				   // this will only work if the cycle
-				   // number is less than ten, which
-				   // we check by an assertion.
+				   // Write the grid in eps format.
   std::string filename = "grid-";
   filename += ('0' + cycle);
   Assert (cycle < 10, ExcInternalError());
   
   filename += ".eps";
+  cout << "Writing grid to <" << filename << ">..." << endl;
   std::ofstream eps_output (filename.c_str());
 
-				   // Using this filename, we write
-				   // each grid as a postscript file.
   GridOut grid_out;
   grid_out.write_eps (triangulation, eps_output);
-
-				   // output of the solution
+  
+				   // Output of the solution in
+				   // gnuplot format.
   filename = "sol-";
   filename += ('0' + cycle);
   Assert (cycle < 10, ExcInternalError());
   
   filename += ".gnuplot";
+  cout << "Writing solution to <" << filename << ">..." << endl;
   std::ofstream gnuplot_output (filename.c_str());
   
   DataOut<dim> data_out;
   data_out.attach_dof_handler (dof_handler);
-  data_out.add_data_vector (solution, "u");
+  data_out.add_data_vector (solution2, "u");
 
   data_out.build_patches ();
   
@@ -582,11 +1488,23 @@ void TransportProblem<dim>::output_results (const unsigned int cycle) const
 };
 
 
-
+				 // The following ``run'' function is
+				 // similar to previous examples. The
+				 // only difference is that the
+				 // problem is assembled and solved
+				 // twice on each refinement step;
+				 // first by ``assemble_system1'' that
+				 // implements the first version and
+				 // then by ``assemble_system2'' that
+				 // implements the second version of
+				 // writing the DG
+				 // discretization. Furthermore the
+				 // time needed by each of the two
+				 // assembling routines is measured.
 template <int dim>
-void TransportProblem<dim>::run () 
+void DGMethod<dim>::run () 
 {
-  for (unsigned int cycle=0; cycle<3; ++cycle)
+  for (unsigned int cycle=0; cycle<6; ++cycle)
     {
       std::cout << "Cycle " << cycle << ':' << std::endl;
 
@@ -597,18 +1515,7 @@ void TransportProblem<dim>::run ()
 	  triangulation.refine_global (3);
 	}
       else
-					 // In case this is not the
-					 // first cycle, we want to
-					 // refine the grid. Unlike
-					 // the global refinement
-					 // employed in the last
-					 // example, we now use the
-					 // adaptive procedure
-					 // described in the function
-					 // which we now call:
-	{
-	  refine_grid ();
-	};
+	refine_grid ();
       
 
       std::cout << "   Number of active cells:       "
@@ -620,97 +1527,62 @@ void TransportProblem<dim>::run ()
       std::cout << "   Number of degrees of freedom: "
 		<< dof_handler.n_dofs()
 		<< std::endl;
-      
-      assemble_system ();
-      solve ();
+
+				       // The constructor of the Timer
+				       // class automatically starts
+				       // the time measurement.
+      Timer assemble_timer;
+				       // First assembling routine.
+      assemble_system1 ();
+				       // The operator () accesses the
+				       // current time without
+				       // disturbing the time
+				       // measurement.
+      cout << "Time of assemble_system1: " << assemble_timer() << endl;
+      solve (solution1);
+
+				       // As preparation for the
+				       // second assembling routine we
+				       // reinit the system matrix, the
+				       // right hand side vector and
+				       // the Timer object.
+      system_matrix.reinit();
+      right_hand_side.clear();
+      assemble_timer.reset();
+
+				       // We start the Timer,
+      assemble_timer.start();
+				       // call the second assembling routine
+      assemble_system2 ();
+				       // and access the current time.
+      cout << "Time of assemble_system2: " << assemble_timer() << endl;
+      solve (solution2);
+
+				       // To make sure that both
+				       // versions of the DG method
+				       // yield the same
+				       // discretization and hence the
+				       // same solution we check the
+				       // two solutions for equality.
+      solution1-=solution2;
+      const double difference=solution1.linfty_norm();
+      if (difference<1e-13)
+	cout << "solution1 and solution2 do not differ." << endl;
+
+				       // Finally we perform the
+				       // output.
       output_results (cycle);
     }
 }
 
+
+
 int main () 
 {
+  DGMethod<2> dgmethod_2d;
+  dgmethod_2d.run ();
 
-				   // The general idea behind the
-				   // layout of this function is as
-				   // follows: let's try to run the
-				   // program as we did before...
-  try
-    {
-      TransportProblem<2> Transport_problem_2d;
-      Transport_problem_2d.run ();
-    }
-				   // ...and if this should fail, try
-				   // to gather as much information as
-				   // possible. Specifically, if the
-				   // exception that was thrown is an
-				   // object of a class that is
-				   // derived from the C++ standard
-				   // class ``exception'', then we can
-				   // use the ``what'' member function
-				   // to get a string which describes
-				   // the reason why the exception was
-				   // thrown. 
-				   //
-				   // The deal.II exception classes
-				   // are all derived from the
-				   // standard class, and in
-				   // particular, the ``exc.what()''
-				   // function will return
-				   // approximately the same string as
-				   // would be generated if the
-				   // exception was thrown using the
-				   // ``Assert'' macro. You have seen
-				   // the output of such an exception
-				   // in the previous example, and you
-				   // then know that it contains the
-				   // file and line number of where
-				   // the exception occured, and some
-				   // other information. This is also
-				   // what would be printed in the
-				   // following.
-  catch (std::exception &exc)
-    {
-      std::cerr << std::endl << std::endl
-		<< "----------------------------------------------------"
-		<< std::endl;
-      std::cerr << "Exception on processing: " << std::endl
-		<< exc.what() << std::endl
-		<< "Aborting!" << std::endl
-		<< "----------------------------------------------------"
-		<< std::endl;
-				       // We can't do much more than
-				       // printing as much information
-				       // as we can get to, so abort
-				       // with error:
-      return 1;
-    }
-				   // If the exception that was thrown
-				   // somewhere was not an object of a
-				   // class derived from the standard
-				   // ``exception'' class, then we
-				   // can't do anything at all. We
-				   // then simply print an error
-				   // message and exit.
-  catch (...) 
-    {
-      std::cerr << std::endl << std::endl
-		<< "----------------------------------------------------"
-		<< std::endl;
-      std::cerr << "Unknown exception!" << std::endl
-		<< "Aborting!" << std::endl
-		<< "----------------------------------------------------"
-		<< std::endl;
-      return 1;
-    };
-
-				   // If we got to this point, there
-				   // was no exception which
-				   // propagated up to the main
-				   // function (maybe there were some,
-				   // but they were caught somewhere
-				   // in the program or the
-				   // library). Therefore, the program
-				   // performed as was expected and we
-				   // can return without error.
   return 0;
 };
+
+
-- 
2.39.5