]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Small edits.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 8 Jan 2010 09:12:39 +0000 (09:12 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 8 Jan 2010 09:12:39 +0000 (09:12 +0000)
git-svn-id: https://svn.dealii.org/trunk@20331 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-38/doc/intro.dox
deal.II/examples/step-38/step-38.cc

index fdc96eed63b0988abb0e8a5a43063fd7bf9e114c..7eca8fa7bcd8a70a9dd34ee9679bf3a3cb64af2c 100644 (file)
@@ -1,21 +1,23 @@
 <a name="Intro"></a>
-<h1>Example of the MeshWorker framework with an advection problem</h1>
+<h1>An example of the MeshWorker framework with an advection problem</h1>
 
 <h3>Overview</h3>
 
-This example is devoted to the MeshWorker framework and the
-<em>discontinuous Galerkin method</em>, or in short: DG method. It
-solves the same problem as @ref step_12 "step-12" (see there for a
-description of the problem and discretization), but here we use the
-MeshWorker framework in order to save reprogramming the cell/face
-loops. We have tried to strip this example of peripheral information
-such that the structure becomes more clear.
+This example is devoted to the MeshWorker framework and the <em>discontinuous
+Galerkin method</em>, or in short: DG method. It solves the same problem as
+@ref step_12 "step-12" (see there for a description of the problem and
+discretization), but here we use the MeshWorker framework in order to save
+programming the cell/face loops that are often rather. The aim of the
+MeshWorker framework is to simplify this process, by putting the majority of
+the boring setup into a framework class and leaving to user code only things
+that are specific to the application.  We have tried to strip this example of
+peripheral information such that the structure becomes more clear.
 
-In particular the loops of DG methods turn out to be complex, because
-for the face terms, we have to distinguish the cases of boundary,
-regular interior faces and interior faces with hanging nodes,
-respectively. The MeshWorker framework implements the standard loop
-over all cells and faces in MeshWorker::loop() and takes care of
+The particular concern of this program are the loops of DG methods. These turn
+out to be especially complex, primarily because for the face terms, we have to
+distinguish the cases of boundary, regular interior faces and interior faces
+with hanging nodes, respectively. The MeshWorker framework implements the
+standard loop over all cells and faces in MeshWorker::loop() and takes care of
 distinguishing between all the different faces.
 
 There are two things left to do if you use MeshWorker: first, you need
index 090e4edb867d9f54d2d86c07f25e0fcec0eac2c4..383a724020fca4dfbc8c7d116dbce58dd8da3bc2 100644 (file)
@@ -3,7 +3,7 @@
 
 /*    $Id$       */
 /*                                                                */
-/*    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors */
+/*    Copyright (C) 2010 by the deal.II authors */
 /*                                                                */
 /*    This file is subject to QPL and may not be  distributed     */
 /*    without copyright and license information. Please refer     */
@@ -13,7 +13,7 @@
                                 // The first few files have already
                                 // been covered in step-12
                                 // and will thus not be further
-                                // commented on.
+                                // commented on:
 #include <base/quadrature_lib.h>
 #include <base/function.h>
 #include <lac/vector.h>
@@ -37,7 +37,7 @@
 #include <base/timer.h>
 
                                 // Here come the new include files
-                                // for using the MeshWorker framework
+                                // for using the MeshWorker framework:
 #include <numerics/mesh_worker.h>
 #include <numerics/mesh_worker_loop.h>
 
@@ -50,6 +50,10 @@ using namespace dealii;
 
                                 // @sect3{Equation data}
                                 //
+                                // First, we need to describe the
+                                // coefficients in the equation. Here, this
+                                // concerns the boundary values, which we
+                                // choose in the same way as for step-12:
 template <int dim>
 class BoundaryValues:  public Function<dim>
 {
@@ -60,15 +64,14 @@ class BoundaryValues:  public Function<dim>
                             const unsigned int component=0) const;
 };
 
-                                // 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.
+                                // Given the flow direction, 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 BoundaryValues<dim>::value_list(const std::vector<Point<dim> > &points,
                                       std::vector<double> &values,
@@ -93,7 +96,7 @@ void BoundaryValues<dim>::value_list(const std::vector<Point<dim> > &points,
                                 // the MeshWorker framework. Since it
                                 // will be used by
                                 // MeshWorker::AssemblingIntegrator,
-                                // it needs the functions for cell,
+                                // it needs functions for cell,
                                 // boundary and interior face
                                 // integration specified exactly as
                                 // below.
@@ -174,25 +177,24 @@ void DGIntegrator<dim>::cell(CellInfo& info) const
       beta(1) = fe_v.quadrature_point(point)(0);
       beta /= beta.norm();
       
+                                      // We solve a homogeneous
+                                      // equation, thus no right
+                                      // hand side shows up in
+                                      // the cell term.
+                                      // What's left is
+                                      // integrating the matrix entries.
       for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) 
-       {
-                                          // We solve a homogeneous
-                                          // equation, thus no right
-                                          // hand side shows up in
-                                          // the cell term.
-                                          // What's left is
-                                          // integrating the matrix entries.
-         for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
-           local_matrix(i,j) -= beta*fe_v.shape_grad(i,point)*
-                                fe_v.shape_value(j,point) *
-                                JxW[point];
-       }
+       for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
+         local_matrix(i,j) -= beta*fe_v.shape_grad(i,point)*
+                              fe_v.shape_value(j,point) *
+                              JxW[point];
     }
 }
 
-                                // Now the same for the boundary
-                                // terms. Note that now we use
-                                // FEFaceValuesBase in order to get
+                                // Now the same for the boundary terms. Note
+                                // that now we use FEFaceValuesBase, the base
+                                // class for both FEFaceValues and
+                                // FESubfaceValues, in order to get access to
                                 // normal vectors.
 template <int dim>
 void DGIntegrator<dim>::bdry(FaceInfo& info) const
@@ -257,7 +259,7 @@ void DGIntegrator<dim>::face(FaceInfo& info1, FaceInfo& info2) const
                                   // four local matrices. The letters
                                   // u and v refer to trial and test
                                   // functions, respectively. The
-                                  // numbers indicate the cells
+                                  // %numbers indicate the cells
                                   // provided by info1 and info2. By
                                   // convention, the two matrices in
                                   // each info object refer to the
@@ -277,7 +279,7 @@ void DGIntegrator<dim>::face(FaceInfo& info1, FaceInfo& info2) const
                                   // vectors. Fortunately, the
                                   // interface terms only involve the
                                   // solution and the right hand side
-                                  // does not obtain a contribution.
+                                  // does not receive any contributions.
   
   const std::vector<double> &JxW = fe_v.get_JxW_values ();
   const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
@@ -293,7 +295,7 @@ void DGIntegrator<dim>::face(FaceInfo& info1, FaceInfo& info2) const
       if (beta_n>0)
        {
                                           // This term we've already
-                                          // seen.
+                                          // 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)
              u1_v1_matrix(i,j) += beta_n *
@@ -315,7 +317,7 @@ void DGIntegrator<dim>::face(FaceInfo& info1, FaceInfo& info2) const
       else
        {
                                           // This one we've already
-                                          // seen, too.
+                                          // seen, too:
          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)
              u2_v1_matrix(i,l) += beta_n *
@@ -326,7 +328,7 @@ void DGIntegrator<dim>::face(FaceInfo& info1, FaceInfo& info2) const
                                           // And this is another new
                                           // one: $(\beta\cdot n \hat
                                           // u,\hat v)_{\partial
-                                          // \kappa_-}$.
+                                          // \kappa_-}$:
          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)
              u2_v2_matrix(k,l) -= beta_n *
@@ -375,32 +377,41 @@ class DGMethod
                                     // use a DG method of a different
                                     // degree the whole program stays
                                     // the same, only replace 1 in
-                                    // the constructor by the wanted
-                                    // degree.
+                                    // the constructor by the desired
+                                    // polynomial degree.
     FE_DGQ<dim>          fe;
     DoFHandler<dim>      dof_handler;
 
     SparsityPattern      sparsity_pattern;
     SparseMatrix<double> system_matrix;
 
-                                    // And there are two solution
-                                    // vectors, that store the
-                                    // solutions to the problems
-                                    // corresponding to the two
+                                    // In step-12 we had two solution vectors
+                                    // that stored the solutions to the
+                                    // problems corresponding to the two
                                     // different assembling routines
                                     // <code>assemble_system1</code> and
-                                    // <code>assemble_system2</code>;
+                                    // <code>assemble_system2</code>. In this
+                                    // program, the goal is only to show the
+                                    // MeshWorker framework, so we only
+                                    // assemble the system in one of the two
+                                    // ways, and consequently we have only
+                                    // one solution vector along with the
+                                    // single <code>assemble_system</code>
+                                    // function declared above:
     Vector<double>       solution;
     Vector<double>       right_hand_side;    
 };
 
 
+                                                // We start with the
+                                                // constructor. This is the
+                                                // place to change the
+                                                // polynomial degree of the
+                                                // finite element shape
+                                                // functions.
 template <int dim>
 DGMethod<dim>::DGMethod ()
                :
-                                                // Change here for DG
-                                                // methods of
-                                                // different degrees.
                 fe (1),
                dof_handler (triangulation)
 {}
@@ -413,29 +424,35 @@ DGMethod<dim>::~DGMethod ()
 }
 
 
+                                // In the function that sets up the usual
+                                // finite element data structures, we first
+                                // need to distribute the DoFs.
 template <int dim>
 void DGMethod<dim>::setup_system ()
 {
-                                  // 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 per row
-                                  // is needed when all neighbors of
-                                  // a cell are once more refined
-                                  // than the cell under
+                                  // The DoFs of a cell are coupled with all
+                                  // DoFs of all neighboring cells, along
+                                  // with all of its siblings on the current
+                                  // cell.  Therefore the maximum 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
-                           *GeometryInfo<dim>::max_children_per_face+1)*fe.dofs_per_cell);
+                          (GeometryInfo<dim>::faces_per_cell *
+                           GeometryInfo<dim>::max_children_per_face
+                           +
+                           1)*fe.dofs_per_cell);
   
-                                  // For DG discretizations we call
-                                  // the function analogue to
-                                  // DoFTools::make_sparsity_pattern.
+                                  // To build the sparsity pattern for DG
+                                  // discretizations, we can call the
+                                  // function analogue to
+                                  // DoFTools::make_sparsity_pattern, which
+                                  // is called
+                                  // DoFTools::make_flux_sparsity_pattern:
   DoFTools::make_flux_sparsity_pattern (dof_handler, sparsity_pattern);
   
                                   // All following function calls are
@@ -449,15 +466,15 @@ void DGMethod<dim>::setup_system ()
 }
 
                                 // @sect4{Function: assemble_system}
-                                // Here we see the major difference
-                                // to assembling by hand. Instead of
-                                // writing loops over cells and
-                                // faces, we leave all this to the
-                                // MeshWorker framework. In order to
-                                // do so, we just have to define
-                                // local integration objects and use
-                                // one of the classes in Assembler to
-                                // build the global system. 
+
+                                // Here we see the major difference to
+                                // assembling by hand. Instead of writing
+                                // loops over cells and faces, we leave all
+                                // this to the MeshWorker framework. In order
+                                // to do so, we just have to define local
+                                // integration objects and use one of the
+                                // classes in namespace MeshWorker::Assembler
+                                // to build the global system.
 template <int dim>
 void DGMethod<dim>::assemble_system () 
 {
@@ -485,18 +502,21 @@ void DGMethod<dim>::assemble_system ()
                                   // the global sparse matrix and the
                                   // right hand side vector.
                                   //
-                                  // It turns out,
                                   // MeshWorker::AssemblingIntegrator
-                                  // itself is not that clever at
-                                  // all, but all its capabilities
-                                  // are provided by the two later
+                                  // is not all that clever by itself,
+                                  // but its capabilities
+                                  // are provided by its two latter
                                   // template arguments. By
                                   // exchanging
                                   // MeshWorker::Assembler::SystemSimple,
                                   // we could for instance assemble a
                                   // BlockMatrix or just a Vector
                                   // instead.
-  MeshWorker::AssemblingIntegrator<dim, MeshWorker::Assembler::SystemSimple<SparseMatrix<double>, Vector<double> >, DGIntegrator<dim> >
+  MeshWorker::AssemblingIntegrator
+    <dim,
+     MeshWorker::Assembler::SystemSimple<SparseMatrix<double>,
+                                         Vector<double> >,
+     DGIntegrator<dim> >
     integrator(dg);
   
                                   // First, we initialize the
@@ -512,7 +532,9 @@ void DGMethod<dim>::assemble_system ()
                                   // independently, we have to hand
                                   // over this value three times.
   const unsigned int n_gauss_points = dof_handler.get_fe().degree+1;
-  integrator.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points);
+  integrator.initialize_gauss_quadrature(n_gauss_points,
+                                        n_gauss_points,
+                                        n_gauss_points);
 
                                   // These are the types of values we
                                   // need for integrating our
@@ -521,19 +543,21 @@ void DGMethod<dim>::assemble_system ()
                                   // and interior faces, as well as
                                   // interior neighbor faces, which is
                                   // forced by the four @p true values.
-  UpdateFlags update_flags = update_quadrature_points | update_values | update_gradients;
+  UpdateFlags update_flags = update_quadrature_points |
+                            update_values            |
+                            update_gradients;
   integrator.add_update_flags(update_flags, true, true, true, true);
 
                                   // Finally, we have to tell the
-                                  // assembler base class, where to
+                                  // assembler base class where to
                                   // put the local data. These will
                                   // be our system matrix and the
                                   // right hand side.
   integrator.initialize(system_matrix, right_hand_side);
 
-                                  // Finally, we get to the
+                                  // We are now ready to get to the
                                   // integration loop. @p info_box is
-                                  // an object, that generates the
+                                  // an object that generates the
                                   // extended iterators for cells and
                                   // faces of type
                                   // MeshWorker::IntegrationInfo. Since
@@ -556,25 +580,21 @@ void DGMethod<dim>::assemble_system ()
                                 //
                                 // For this simple problem 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
-                                // structure 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
+                                // Richardson iteration, that represents a
+                                // simple defect correction. This, in
+                                // combination with a block SSOR
+                                // preconditioner, that uses the special
+                                // block matrix structure 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 in the downstream direction
+                                // of the flow, then a block Gauss-Seidel
                                 // preconditioner (see the
                                 // PreconditionBlockSOR class with
-                                // relaxation=1) makes a much better
-                                // job.
+                                // relaxation=1) does a much better job.
 template <int dim>
 void DGMethod<dim>::solve (Vector<double> &solution) 
 {
@@ -585,8 +605,8 @@ void DGMethod<dim>::solve (Vector<double> &solution)
                                   // preconditioner,
   PreconditionBlockSSOR<SparseMatrix<double> > preconditioner;
 
-                                  // we assigned the matrix to it and
-                                  // set the right block size.
+                                  // then assign the matrix to it and
+                                  // set the right block size:
   preconditioner.initialize(system_matrix, fe.dofs_per_cell);
 
                                   // After these preparations we are
@@ -715,18 +735,7 @@ void DGMethod<dim>::output_results (const unsigned int cycle) const
 
 
                                 // The following <code>run</code> function is
-                                // similar to previous examples. The
-                                // only difference is that the
-                                // problem is assembled and solved
-                                // twice on each refinement step;
-                                // first by <code>assemble_system1</code> that
-                                // implements the first version and
-                                // then by <code>assemble_system2</code> that
-                                // implements the second version of
-                                // writing the DG
-                                // discretization. Furthermore the
-                                // time needed by each of the two
-                                // assembling routines is measured.
+                                // similar to previous examples.
 template <int dim>
 void DGMethod<dim>::run () 
 {
@@ -757,15 +766,13 @@ void DGMethod<dim>::run ()
       assemble_system ();
       solve (solution);
 
-                                      // Finally we perform the
-                                      // output.
       output_results (cycle);
     }
 }
 
                                 // The following <code>main</code> function is
-                                // similar to previous examples and
-                                // need not to be commented on.
+                                // similar to previous examples as well, and
+                                // need not be commented on.
 int main () 
 {
   try

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.