]> https://gitweb.dealii.org/ - dealii.git/commitdiff
A couple of changes and corrections.
authorRalf Hartmann <Ralf.Hartmann@dlr.de>
Thu, 22 Nov 2001 15:37:47 +0000 (15:37 +0000)
committerRalf Hartmann <Ralf.Hartmann@dlr.de>
Thu, 22 Nov 2001 15:37:47 +0000 (15:37 +0000)
git-svn-id: https://svn.dealii.org/trunk@5241 0785d39b-7218-0410-832d-ea1e28bc413d

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

index f3dddad75d3088640aff0ec407a3af07f572cda3..8d568d396ab91948d349b34f793e49e2f76d2836 100644 (file)
@@ -1,5 +1,5 @@
 /* $Id$ */
-/* Author: Ralf Hartmann, University of Heidelberg, 2000 */
+/* Author: Ralf Hartmann, University of Heidelberg, 2001 */
 
                                 // The first few files have already
                                 // been covered in previous examples
@@ -38,7 +38,7 @@
                                 // used as all other finite elements.
 #include <fe/fe_dgq.h>
                                 // We are going to use the simplest
-                                // possible solver, called richardson
+                                // possible solver, called Richardson
                                 // iteration, that represents a simple
                                 // defect correction. This, in
                                 // combination with a block SSOR
 #include <fstream>
 
 
-                                // First we define the class
+                                // @sect3{Equation data}
+                                //
+                                // First we define the classes
                                 // representing the equation-specific
                                 // functions. Both classes, ``RHS''
                                 // and ``BoundaryValues'', are
-                                // derived from the Function
+                                // derived from the ``Function''
                                 // class. Only the ``value_list''
                                 // function are implemented because
                                 // only lists of function values are
@@ -74,8 +76,6 @@ template <int dim>
 class RHS:  public Function<dim>
 {
   public:
-    RHS() {};
-    
     virtual void value_list (const std::vector<Point<dim> > &points,
                             std::vector<double> &values,
                             const unsigned int component=0) const;
@@ -86,18 +86,16 @@ template <int dim>
 class BoundaryValues:  public Function<dim>
 {
   public:
-    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
+                                // The class ``Beta'' represents the
+                                // vector valued flow field of the
+                                // linear transport equation and is
+                                // not derived from the ``Function''
                                 // class as we prefer to get function
                                 // values of type ``Point'' rather
                                 // than of type
@@ -112,8 +110,6 @@ template <int dim>
 class Beta
 {
   public:
-    Beta () {};
-
     void value_list (const std::vector<Point<dim> > &points,
                     std::vector<Point<dim> > &values) const;
 };
@@ -123,35 +119,41 @@ class Beta
                                 // ``value_list'' functions of these
                                 // classes are rather simple.  For
                                 // simplicity the right hand side is
-                                // set to be zero.
+                                // set to be zero but will be
+                                // assembled anyway.
 template <int dim>
-void RHS<dim>::value_list(const std::vector<Point<dim> > &,
+void RHS<dim>::value_list(const std::vector<Point<dim> > &points,
                          std::vector<double> &values,
                          const unsigned int) const
 {
+                                  // Usually we check whether input
+                                  // parameter have the right sizes.
+  Assert(values.size()==points.size(),
+        ExcDimensionMismatch(values.size(),points.size()));
+
   for (unsigned int i=0; i<values.size(); ++i)
     values[i]=0;
 }
 
 
                                 // The flow field is chosen to be
-                                // circular, anticlockwise, and with
+                                // circular, counterclockwise, and with
                                 // the origin as midpoint.
-template <>
-void Beta<2>::value_list(const std::vector<Point<2> > &points,
-                        std::vector<Point<2> > &values) const
+template <int dim>
+void Beta<dim>::value_list(const std::vector<Point<dim> > &points,
+                          std::vector<Point<dim> > &values) const
 {
   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];
-      Point<2> &beta=values[i];
+      const Point<dim> &p=points[i];
+      Point<dim> &beta=values[i];
 
       beta(0)=-p(1);
       beta(1)=p(0);
-      beta/=sqrt(beta*beta);
+      beta/=sqrt(beta.square());
     }
 }
 
@@ -185,9 +187,10 @@ void BoundaryValues<dim>::value_list(const std::vector<Point<dim> > &points,
 
                                 // @sect3{Class: DGTransportEquation}
                                 //
-                                // Next we define the equation-
-                                // dependent and DG-method-dependent
-                                // class ``DGTransportEquation''. Its
+                                // 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
@@ -204,27 +207,27 @@ class DGTransportEquation
 
     void assemble_cell_term(const FEValues<dim>& fe_v,
                            FullMatrix<double> &u_v_matrix,
-                           Vector<double> &cell_vector);
+                           Vector<double> &cell_vector) const;
     
     void assemble_boundary_term(const FEFaceValues<dim>& fe_v,
                                FullMatrix<double> &u_v_matrix,
-                               Vector<double> &cell_vector);
+                               Vector<double> &cell_vector) const;
     
     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);
+                            FullMatrix<double> &un_v_matrix) const;
 
     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);
+                            FullMatrix<double> &un_vn_matrix) const;
   private:
-    Beta<dim> beta_function;
-    RHS<dim> rhs_function;
-    BoundaryValues<dim> boundary_function;
+    const Beta<dim> beta_function;
+    const RHS<dim> rhs_function;
+    const BoundaryValues<dim> boundary_function;
 };
 
 
@@ -238,27 +241,27 @@ class DGTransportEquation
                                 // 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
+                                // ``fe_v'' is already reinit'ed with the
                                 // current cell before and includes
                                 // all shape values needed.
 template <int dim>
 void DGTransportEquation<dim>::assemble_cell_term(
-  const FEValues<dim>fe_v,
+  const FEValues<dim> &fe_v,
   FullMatrix<double> &u_v_matrix,
-  Vector<double> &cell_vector)
+  Vector<double> &cell_vector) const
 {
                                   // First we ask ``fe_v'' for the
                                   // shape gradients, shape values and
                                   // quadrature weights,
-  const vector<vector<Tensor<1,2> > > &grad_v = fe_v.get_shape_grads ();
+  const std::vector<std::vector<Tensor<1,dim> > > &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 ();
+  const std::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);
+  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
+  std::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);
@@ -286,7 +289,7 @@ void DGTransportEquation<dim>::assemble_cell_term(
                                 // function assembles the face terms
                                 // at boundary faces.  When this
                                 // function is invoked, ``fe_v'' is
-                                // already reinited with the current
+                                // already reinit'ed with the current
                                 // cell and current face. Hence it
                                 // provides the shape values on that
                                 // boundary face.
@@ -294,7 +297,7 @@ template <int dim>
 void DGTransportEquation<dim>::assemble_boundary_term(
   const FEFaceValues<dim>& fe_v,    
   FullMatrix<double> &u_v_matrix,
-  Vector<double> &cell_vector)
+  Vector<double> &cell_vector) const
 {
                                   // First we check whether the
                                   // current face is really at the
@@ -302,19 +305,19 @@ void DGTransportEquation<dim>::assemble_boundary_term(
   Assert(fe_v.get_face()->at_boundary(), ExcInternalError());
   
                                   // Again, as in the previous
-                                  // function, we ask the FEValues
+                                  // function, we ask the ``FEValues''
                                   // object for the shape values and
                                   // the quadrature weights
   const FullMatrix<double> &v = fe_v.get_shape_values ();
-  const vector<double> &JxW = fe_v.get_JxW_values ();
+  const std::vector<double> &JxW = fe_v.get_JxW_values ();
                                   // but here also for the normals.
-  const vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
+  const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
 
                                   // We evaluate the flow field
                                   // and the boundary values at the
                                   // quadrature points.
-  vector<Point<dim> > beta (fe_v.n_quadrature_points);
-  vector<double> g(fe_v.n_quadrature_points);
+  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
+  std::vector<double> g(fe_v.n_quadrature_points);
   
   beta_function.value_list (fe_v.get_quadrature_points(), beta);
   boundary_function.value_list (fe_v.get_quadrature_points(), g);
@@ -325,11 +328,11 @@ void DGTransportEquation<dim>::assemble_boundary_term(
                                   // introduction.
   for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
     {
-      double beta_n=beta[point] * normals[point];      
-      if (beta_n>0)
+      const double beta_n=beta[point] * normals[point];      
                                         // We assemble the term
                                         // $(\beta\cdot n
                                         // u,v)_{\partial K_+}$,
+      if (beta_n>0)
        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 *
@@ -361,7 +364,7 @@ void DGTransportEquation<dim>::assemble_boundary_term(
                                 //
                                 // When this function is invoked,
                                 // ``fe_v'' and ``fe_v_neighbor'' are
-                                // already reinited with the current
+                                // already reinit'ed with the current
                                 // cell and the neighoring cell,
                                 // respectively, as well as with the
                                 // current face. Hence they provide
@@ -375,7 +378,7 @@ void DGTransportEquation<dim>::assemble_boundary_term(
                                 // contributions to the system matrix
                                 // that are based on outer values of
                                 // u, see $\hat u_h$ in the
-                                // Introduction, and inner values of
+                                // introduction, and inner values of
                                 // v, see $v_h$. Here we note that
                                 // ``un'' is the short notation for
                                 // ``u_neighbor'' and represents
@@ -385,7 +388,7 @@ 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)
+  FullMatrix<double> &un_v_matrix) const
 {
                                   // First we check that the current
                                   // face is not at the boundary by
@@ -399,12 +402,12 @@ void DGTransportEquation<dim>::assemble_face_term1(
                                   // normals
   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 ();
+  const std::vector<double> &JxW = fe_v.get_JxW_values ();
+  const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
 
                                   // and we evaluate the flow field
                                   // at the quadrature points.
-  vector<Point<dim> > beta (fe_v.n_quadrature_points);
+  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
   beta_function.value_list (fe_v.get_quadrature_points(), beta);
 
                                   // Then we assemble the cell
@@ -413,11 +416,11 @@ void DGTransportEquation<dim>::assemble_face_term1(
                                   // introduction.
   for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
     {
-      double beta_n=beta[point] * normals[point];
-      if (beta_n>0)
+      const double beta_n=beta[point] * normals[point];
                                         // We assemble the term
                                         // $(\beta\cdot n
                                         // u,v)_{\partial K_+}$,
+      if (beta_n>0)
        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 *
@@ -462,21 +465,21 @@ void DGTransportEquation<dim>::assemble_face_term2(
   FullMatrix<double> &u_v_matrix,
   FullMatrix<double> &un_v_matrix,
   FullMatrix<double> &u_vn_matrix,
-  FullMatrix<double> &un_vn_matrix)
+  FullMatrix<double> &un_vn_matrix) const
 {
                                   // 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 ();
+  const std::vector<double> &JxW = fe_v.get_JxW_values ();
+  const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
 
-  vector<Point<dim> > beta (fe_v.n_quadrature_points);
+  std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
   
   beta_function.value_list (fe_v.get_quadrature_points(), beta);
 
   for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
     {
-      double beta_n=beta[point] * normals[point];
+      const double beta_n=beta[point] * normals[point];
       if (beta_n>0)
        {
                                           // This terms we've already seen.
@@ -533,8 +536,8 @@ void DGTransportEquation<dim>::assemble_face_term2(
                                 // 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.
+                                // because there are no hanging node
+                                // constraints in DG discretizations.
 template <int dim>
 class DGMethod
 {
@@ -553,12 +556,17 @@ class DGMethod
     void output_results (const unsigned int cycle) const;
     
     Triangulation<dim>   triangulation;
-    MappingQ1<dim>       mapping;
+    const MappingQ1<dim> mapping;
     
-                                    // Furthermore we want to
-                                    // use DG elements of degree 1
-                                    // (but this is only specified in
-                                    // the constructor):
+                                    // Furthermore we want to use DG
+                                    // elements of degree 1 (but this
+                                    // is only specified in the
+                                    // constructor). If you want to
+                                    // use a DG method of a different
+                                    // degree the whole program stays
+                                    // the same, only replace 1 in
+                                    // the constructor by the wanted
+                                    // degree.
     FE_DGQ<dim>          fe;
     DoFHandler<dim>      dof_handler;
 
@@ -569,8 +577,8 @@ class DGMethod
                                     // formulae for the cell and the
                                     // face terms of the
                                     // discretization.
-    QGauss4<dim>   quadrature;
-    QGauss4<dim-1> face_quadrature;
+    const QGauss4<dim>   quadrature;
+    const QGauss4<dim-1> face_quadrature;
     
                                     // And there are two solution
                                     // vectors, that store the
@@ -587,18 +595,15 @@ class DGMethod
                                     // object of the
                                     // DGTransportEquations class
                                     // described above.
-    DGTransportEquation<dim>     dg;
+    const 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>
 DGMethod<dim>::DGMethod () :
+                                                // Change here for DG
+                                                // methods of
+                                                // different degrees.
                 fe (1),
                dof_handler (triangulation)
 {}
@@ -621,10 +626,11 @@ void DGMethod<dim>::setup_system ()
                                   // 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.
+                                  // number of matrix entries per row
+                                  // 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(),
                           (GeometryInfo<dim>::faces_per_cell
@@ -657,67 +663,75 @@ void DGMethod<dim>::setup_system ()
                                 // ``assemble_cell_term'',
                                 // ``assemble_boundary_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:
+                                // of the ``DGTransportEquation''
+                                // object.  The
+                                // ``assemble_boundary_term'' covers
+                                // the first case mentioned in the
+                                // introduction.
+                                //
+                                // 1. face is at boundary
                                 //
-                                // 1. face is at boundary (current
-                                // cell: FEFaceValues);
+                                // This function takes a
+                                // ``FEFaceValues'' object as
+                                // argument.  In contrast to that
+                                // ``assemble_face_term1''
+                                // 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
+                                // the remaining cases mentioned
+                                // in the introduction:
                                 //
                                 // 2. neighboring cell is finer
-                                // (current cell: FESubfaceValues,
-                                // neighboring cell: FEFaceValues);
+                                // (current cell: ``FESubfaceValues'',
+                                // neighboring cell: ``FEFaceValues'');
                                 //
                                 // 3. neighboring cell is of the same
                                 // refinement level (both, current
                                 // and neighboring cell:
-                                // FEFaceValues);
+                                // ``FEFaceValues'');
                                 //
                                 // 4. neighboring cell is coarser
-                                // (current cell: FEFaceValues,
+                                // (current cell: ``FEFaceValues'',
                                 // neighboring cell:
-                                // FESubfaceValues).
+                                // ``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.
+                                // meshes then only case 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 DGMethod<dim>::assemble_system1 () 
 {
   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);
-  
+  std::vector<unsigned int> dofs (dofs_per_cell);
+  std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
+
+                                  // First we create the
+                                  // ``UpdateFlags'' for the
+                                  // ``FEValues'' and the
+                                  // ``FEFaceValues'' objects.
+  UpdateFlags update_flags = 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);
+  UpdateFlags face_update_flags = update_values
+                                 | update_q_points
+                                 | update_JxW_values
+                                 | update_normal_vectors;
   
                                   // On the neighboring cell we only
                                   // need the shape values. Given a
@@ -727,34 +741,36 @@ void DGMethod<dim>::assemble_system1 ()
                                   // 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);
+  UpdateFlags neighbor_face_update_flags = update_values;
    
-                                  // Then we create the FEValues
+                                  // 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
+                                  // 3.2.0 of deal.II 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.
+                                  // 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'' 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 (
@@ -774,8 +790,8 @@ void DGMethod<dim>::assemble_system1 ()
                                   // 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'
+                                  // the convention that `un' is
+                                  // the shortcut for `u_neighbor'
                                   // and represents the $u_hat$, see
                                   // introduction.
   FullMatrix<double> u_v_matrix (dofs_per_cell, dofs_per_cell);
@@ -784,34 +800,45 @@ void DGMethod<dim>::assemble_system1 ()
   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>::active_cell_iterator neighbor_child;
+                                  // iterators.
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
 
                                   // Now we start the loop over all
-                                  // active cells
+                                  // active cells.
   for (;cell!=endc; ++cell) 
     {
-                                      // and reinit the FEValues
-                                      // object for the current cell,
+                                      // In the
+                                      // ``assemble_face_term1''
+                                      // function contributions to
+                                      // the cell matrices and the
+                                      // cell vector are only
+                                      // ADDED. Therefore on each
+                                      // cell we need to reset the
+                                      // ``u_v_matrix'' and
+                                      // ``cell_vector'' to zero,
+                                      // before assembling the cell terms.
+      u_v_matrix.clear ();
+      cell_vector.clear ();
+      
+                                      // Now we reinit the ``FEValues''
+                                      // object for the current cell
       fe_v.reinit (cell);
 
-                                      // Call the function that
-                                      // assembles the cell
+                                      // and call the function
+                                      // that assembles the cell
                                       // terms. The first argument is
-                                      // the FEValues that was
-                                      // already reinited on current
-                                      // the cell.
+                                      // the ``FEValues'' that was
+                                      // previously reinit'ed on the
+                                      // current cell.
       dg.assemble_cell_term(fe_v,
                            u_v_matrix,
                            cell_vector);
 
-                                      // As in previous example steps
-                                      // the vector `dofs' includes
-                                      // the dof_indices.
+                                      // As in previous examples the
+                                      // vector `dofs' includes the
+                                      // dof_indices.
       cell->get_dof_indices (dofs);
 
                                       // This is the start of the
@@ -819,9 +846,14 @@ void DGMethod<dim>::assemble_system1 ()
       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);
-
+                                          // iterator
+         typename DoFHandler<dim>::face_iterator face=cell->face(face_no);
+         
+                                          // and clear the
+                                          // ``un_v_matrix'' on each
+                                          // face.
+         un_v_matrix.clear();
+                 
                                           // Now we distinguish the
                                           // four different cases in
                                           // the ordering mentioned
@@ -831,8 +863,9 @@ void DGMethod<dim>::assemble_system1 ()
          if (face->at_boundary())
            {
                                               // We reinit the
-                                              // FEFaceValues object
-                                              // to the current face
+                                              // ``FEFaceValues''
+                                              // object to the
+                                              // current face
              fe_v_face.reinit (cell, face_no);
 
                                               // and assemble the
@@ -849,7 +882,8 @@ void DGMethod<dim>::assemble_system1 ()
                                               // domain, therefore
                                               // there must exist a
                                               // neighboring cell.
-             neighbor = cell->neighbor(face_no);
+             typename DoFHandler<dim>::cell_iterator neighbor=
+               cell->neighbor(face_no);;
 
                                               // We proceed with the
                                               // second and most
@@ -864,17 +898,21 @@ void DGMethod<dim>::assemble_system1 ()
                                               // difference of not
                                               // more than one, the
                                               // neighboring cell is
-                                              // known to be only
+                                              // known to be at most
                                               // ONCE more refined
                                               // than the current
                                               // cell. Furthermore
                                               // also the face is
-                                              // once more refined,
+                                              // more refined,
                                               // i.e. it has
-                                              // children.
+                                              // children. Here we
+                                              // note that the
+                                              // following part of
+                                              // code will not work
+                                              // for ``dim==1''.
              if (face->has_children())
                {
-                                                  // first we store
+                                                  // First we store
                                                   // which number the
                                                   // current cell has
                                                   // in the list of
@@ -904,34 +942,56 @@ void DGMethod<dim>::assemble_system1 ()
                                                       // `behind' the
                                                       // current
                                                       // subface.
-                     neighbor_child = neighbor->child(GeometryInfo<dim>::
-                                                      child_cell_on_face(neighbor2,subface_no));
-
+                     typename DoFHandler<dim>::active_cell_iterator neighbor_child=
+                       neighbor->child(GeometryInfo<dim>::
+                                       child_cell_on_face(neighbor2,subface_no));
+                     
                                                       // As these are
                                                       // quite
                                                       // complicated
                                                       // indirections
-                                                      // we check for
+                                                      // which one
+                                                      // does not
+                                                      // usually get
+                                                      // right at
+                                                      // first
+                                                      // attempt we
+                                                      // check for
                                                       // the internal
                                                       // consistency.
                      Assert (neighbor_child->face(neighbor2) == face->child(subface_no),
                              ExcInternalError());
                      Assert (!neighbor_child->has_children(), ExcInternalError());
 
+                                                      // We need to
+                                                      // reset the
+                                                      // ``un_v_matrix''
+                                                      // on each
+                                                      // subface
+                                                      // because on
+                                                      // each subface
+                                                      // the ``un''
+                                                      // belong to
+                                                      // different
+                                                      // neighboring
+                                                      // cells.
+                     un_v_matrix.clear();
+                     
                                                       // As already
                                                       // mentioned
                                                       // above for
-                                                      // this case
-                                                      // (case 2) we
-                                                      // employ the
-                                                      // FESubfaceValues
+                                                      // the current
+                                                      // case (case
+                                                      // 2) we employ
+                                                      // the
+                                                      // ``FESubfaceValues''
                                                       // of the
                                                       // current
-                                                      // cellhere
+                                                      // cell (here
                                                       // reinited for
                                                       // the current
                                                       // cell, face
-                                                      // and subface,
+                                                      // and subface)
                                                       // and we
                                                       // employ the
                                                       // FEFaceValues
@@ -946,38 +1006,27 @@ void DGMethod<dim>::assemble_system1 ()
                                             u_v_matrix,
                                             un_v_matrix);
                      
-                                                      // get dof
+                                                      // Then we get
+                                                      // the dof
                                                       // indices of
                                                       // the
                                                       // neighbor_child
                                                       // cell
                      neighbor_child->get_dof_indices (dofs_neighbor);
                                                                
+                                                      // and
                                                       // distribute
-                                                      // cell matrix
+                                                      // ``un_v_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],
                                            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();
                    }
+                                                  // End of ``if
+                                                  // (face->has_children())''
                }
-                                              // End of ``if
-                                              // (face->has_children())''
              else
                {
                                                   // We proceed with
@@ -989,19 +1038,21 @@ void DGMethod<dim>::assemble_system1 ()
                                                   // current cell.
                  if (neighbor->level() == cell->level()) 
                    {
-                                                  // Like before we
-                                                  // store which
-                                                  // number the
-                                                  // current cell has
-                                                  // in the list of
-                                                  // neighbors of the
-                                                  // neighboring
-                                                  // cell.
+                                                      // 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
+                                                      // ``FEFaceValues''
                                                       // of the
                                                       // current and
                                                       // neighboring
@@ -1061,7 +1112,7 @@ void DGMethod<dim>::assemble_system1 ()
 
                                                       // Reinit the
                                                       // appropriate
-                                                      // FEFaceValues
+                                                      // ``FEFaceValues''
                                                       // and assemble
                                                       // the face
                                                       // terms.
@@ -1075,23 +1126,19 @@ void DGMethod<dim>::assemble_system1 ()
                                             un_v_matrix);
                    }
 
-                                                  // Get dof indices
-                                                  // of the
-                                                  // neighbor_child
+                                                  // Now we get the
+                                                  // dof indices of
+                                                  // the
+                                                  // ``neighbor_child''
                                                   // cell,
                  neighbor->get_dof_indices (dofs_neighbor);
                                                                
-                                                  // distribute the
-                                                  // un_v_matrix,
+                                                  // and 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],
                                        un_v_matrix(i,k));
-                 
-                                                  // and clear the
-                                                  // ``un_v_matrix''
-                                                  // on each face.
-                 un_v_matrix.clear();
                }
                                               // End of ``face not at boundary'':
            }
@@ -1099,18 +1146,14 @@ void DGMethod<dim>::assemble_system1 ()
        }
       
                                       // Finally we distribute the
-                                      // u_v_matrix,
+                                      // ``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], u_v_matrix(i,j));
       
-                                      // the cell vector
+                                      // and the cell vector.
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        right_hand_side(dofs[i]) += cell_vector(i);
-
-                                      // and clear them both.
-      u_v_matrix.clear ();
-      cell_vector.clear ();
     }
 };
 
@@ -1155,20 +1198,20 @@ template <int dim>
 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);
+  std::vector<unsigned int> dofs (dofs_per_cell);
+  std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
 
-  UpdateFlags update_flags = UpdateFlags(update_values
-                                        | update_gradients
-                                        | update_q_points
-                                        | update_JxW_values);
+  UpdateFlags update_flags = 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);
+  UpdateFlags face_update_flags = update_values
+                                 | update_q_points
+                                 | update_JxW_values
+                                 | update_normal_vectors;
+  
+  UpdateFlags neighbor_face_update_flags = update_values;
 
                                   // Here we do not need
                                   // ``fe_v_face_neighbor'' as case 4
@@ -1186,15 +1229,16 @@ void DGMethod<dim>::assemble_system2 ()
   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'.
+                                  // Additionally we need the
+                                  // 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);
   
@@ -1203,14 +1247,14 @@ void DGMethod<dim>::assemble_system2 ()
                                   // 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;
-
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
   for (;cell!=endc; ++cell) 
     {
+      u_v_matrix.clear ();
+      cell_vector.clear ();
+
       fe_v.reinit (cell);
 
       dg.assemble_cell_term(fe_v,
@@ -1221,7 +1265,8 @@ void DGMethod<dim>::assemble_system2 ()
 
       for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
        {
-         face = cell->face(face_no);
+         typename DoFHandler<dim>::face_iterator face=
+           cell->face(face_no);
 
                                           // Case 1:
          if (face->at_boundary())
@@ -1235,8 +1280,8 @@ void DGMethod<dim>::assemble_system2 ()
          else
            {
              Assert (cell->neighbor(face_no).state() == valid, ExcInternalError());
-             neighbor = cell->neighbor(face_no);
-
+             typename DoFHandler<dim>::cell_iterator neighbor=
+               cell->neighbor(face_no);
                                               // Case 2:
              if (face->has_children())
                {
@@ -1247,12 +1292,17 @@ void DGMethod<dim>::assemble_system2 ()
                       subface_no<GeometryInfo<dim>::subfaces_per_face;
                       ++subface_no)
                    {
-                     neighbor_child = neighbor->child(
-                       GeometryInfo<dim>::child_cell_on_face(neighbor2,subface_no));
+                     typename DoFHandler<dim>::cell_iterator 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());
-
+                     
+                     un_v_matrix.clear();
+                     u_vn_matrix.clear();
+                     un_vn_matrix.clear();
+                     
                      fe_v_subface.reinit (cell, face_no, subface_no);
                      fe_v_face_neighbor.reinit (neighbor_child, neighbor2);
 
@@ -1275,10 +1325,6 @@ void DGMethod<dim>::assemble_system2 ()
                            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
@@ -1290,7 +1336,11 @@ void DGMethod<dim>::assemble_system2 ()
                      neighbor->index() > cell->index()) 
                    {
                      const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
-
+                     
+                     un_v_matrix.clear();
+                     u_vn_matrix.clear();
+                     un_vn_matrix.clear();
+                     
                      fe_v_face.reinit (cell, face_no);
                      fe_v_face_neighbor.reinit (neighbor, neighbor2);
                      
@@ -1313,10 +1363,6 @@ void DGMethod<dim>::assemble_system2 ()
                            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)
@@ -1333,15 +1379,15 @@ void DGMethod<dim>::assemble_system2 ()
       
       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
+
+                                // @sect3{All the rest}
+                                //
+                                // For this simple problem we use the
                                 // simplest possible solver, called
-                                // richardson iteration, that
+                                // Richardson iteration, that
                                 // represents a simple defect
                                 // correction. This, in combination
                                 // with a block SSOR preconditioner,
@@ -1349,8 +1395,8 @@ void DGMethod<dim>::assemble_system2 ()
                                 // 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
+                                // 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
@@ -1380,8 +1426,9 @@ void DGMethod<dim>::solve (Vector<double> &solution)
                                   // 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
+                                  // solver. Otherwise, it takes less
+                                  // memory, but the diagonal blocks
+                                  // are inverted in each
                                   // preconditioner step,
                                   // significantly slowing down the
                                   // linear solving process.
@@ -1396,42 +1443,47 @@ void DGMethod<dim>::solve (Vector<double> &solution)
 
                                 // 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
+                                // namely an approximation to the
+                                // gradient 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
+                                // ``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
+                                // ``GradientEstimation'' described
+                                // in Step 9 of this tutorial. In
+                                // fact, the
+                                // ``DerivativeApproximation'' class
+                                // was developed following the
+                                // ``GradientEstimation'' class of
+                                // Step 9. Relating to the
+                                // discussion 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.
+                                // 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 DGMethod<dim>::refine_grid ()
 {
-                                  // The DerivativeApproximation
+                                  // The ``DerivativeApproximation''
                                   // class computes the gradients to
                                   // float precision. This is
                                   // sufficient as they are
@@ -1492,7 +1544,7 @@ void DGMethod<dim>::output_results (const unsigned int cycle) const
   Assert (cycle < 10, ExcInternalError());
   
   filename += ".gnuplot";
-  cout << "Writing solution to <" << filename << ">..." << endl;
+  cout << "Writing solution to <" << filename << ">..." << endl << endl;
   std::ofstream gnuplot_output (filename.c_str());
   
   DataOut<dim> data_out;
@@ -1583,8 +1635,8 @@ void DGMethod<dim>::run ()
                                       // two solutions for equality.
       solution1-=solution2;
       const double difference=solution1.linfty_norm();
-      if (difference<1e-13)
-       cout << "solution1 and solution2 do not differ." << endl;
+      if (difference>1e-13)
+       cout << "solution1 and solution2 differ!!" << endl;
 
                                       // Finally we perform the
                                       // output.
@@ -1592,13 +1644,40 @@ void DGMethod<dim>::run ()
     }
 }
 
-
-
+                                // The following ``main'' function is
+                                // similar to previous examples and
+                                // need not to be commented on.
 int main () 
 {
-  DGMethod<2> dgmethod_2d;
-  dgmethod_2d.run ();
-
+  try
+    {
+      DGMethod<2> dgmethod;
+      dgmethod.run ();
+    }
+  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;
+      return 1;
+    }
+  catch (...) 
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    };
+  
   return 0;
 };
 

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.