]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Review the documentation.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 26 Jan 2006 06:52:55 +0000 (06:52 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 26 Jan 2006 06:52:55 +0000 (06:52 +0000)
git-svn-id: https://svn.dealii.org/trunk@12173 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 071dcc93508dc7c074827c06a6f51d524613f436..3d9991a8189446887a2e8032b051027c9358fa8c 100644 (file)
@@ -4,35 +4,42 @@
 /*    $Id$       */
 /*    Version: $Name$                                          */
 /*                                                                */
-/*    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors */
+/*    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2006 by the deal.II authors */
 /*                                                                */
 /*    This file is subject to QPL and may not be  distributed     */
 /*    without copyright and license information. Please refer     */
 /*    to the file deal.II/doc/license.html for the  text  and     */
 /*    further information on this license.                        */
 
+
+                                 // @sect3{Many new include files}
+
                                 // These include files are already
                                 // known to you. They declare the
                                 // classes which handle
-                                // triangulations and enumerate the
-                                // degrees of freedom.
+                                // triangulations and enumeration of
+                                // degrees of freedom:
 #include <grid/tria.h>
 #include <dofs/dof_handler.h>
                                 // And this is the file in which the
-                                // functions are declared which
-                                // create grids.
+                                // functions are declared that
+                                // create grids:
 #include <grid/grid_generator.h>
 
-                                // The next three files contain
-                                // classes which are needed for loops
-                                // over all cells and to get the
-                                // information from the cell objects.
+                                // The next three files contain classes which
+                                // are needed for loops over all cells and to
+                                // get the information from the cell
+                                // objects. The first two have been used
+                                // before to get geometric information from
+                                // cells; the last one is new and provides
+                                // information about the degrees of freedom
+                                // local to a cell:
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
 #include <dofs/dof_accessor.h>
 
-                                // In this file are the finite
-                                // element descriptions.
+                                // In this file contains the description of
+                                // the Lagrange interpolation finite element:
 #include <fe/fe_q.h>
 
                                 // And this file is needed for the
@@ -45,7 +52,7 @@
                                 // assembling the matrix using
                                 // quadrature on each cell. The
                                 // classes declared in them will be
-                                // explained below.
+                                // explained below:
 #include <fe/fe_values.h>
 #include <base/quadrature_lib.h>
 
 #include <numerics/vectors.h>
 #include <numerics/matrices.h>
 
-                                // These include files are for the
-                                // linear algebra which we employ to
-                                // solve the system of equations
-                                // arising from the finite element
-                                // discretization of the Laplace
-                                // equation. We will use vectors and
-                                // full matrices for assembling the
-                                // system of equations locally on
-                                // each cell, and transfer the
-                                // results into a sparse matrix. We
-                                // will then use a Conjugate Gradient
-                                // solver to solve the problem, for
-                                // which we need a preconditioner (in
-                                // this program, we use the identity
-                                // preconditioner which does nothing,
-                                // but we need to include the file
-                                // anyway), and a class which
-                                // provides the solver with some
-                                // memory for temporary vectors.
+                                // We're now almost to the end. The second to
+                                // last group of include files is for the
+                                // linear algebra which we employ to solve
+                                // the system of equations arising from the
+                                // finite element discretization of the
+                                // Laplace equation. We will use vectors and
+                                // full matrices for assembling the system of
+                                // equations locally on each cell, and
+                                // transfer the results into a sparse
+                                // matrix. We will then use a Conjugate
+                                // Gradient solver to solve the problem, for
+                                // which we need a preconditioner (in this
+                                // program, we use the identity
+                                // preconditioner which does nothing, but we
+                                // need to include the file anyway):
 #include <lac/vector.h>
 #include <lac/full_matrix.h>
 #include <lac/sparse_matrix.h>
 #include <iostream>
 
 
-                                // Instead of the procedural
-                                // programming of previous examples,
-                                // we encapsulate everything into a
-                                // class for this program. The class
-                                // consists of functions which do
-                                // certain aspects of a finite
-                                // element program, a `main' function
-                                // which controls what is done first
-                                // and what is done next, and a list
-                                // of member variables.
+                                 // @sect3{The ``LaplaceProblem'' class}
+
+                                // Instead of the procedural programming of
+                                // previous examples, we encapsulate
+                                // everything into a class for this
+                                // program. The class consists of functions
+                                // which each perform certain aspects of a
+                                // finite element program, a `main' function
+                                // which controls what is done first and what
+                                // is done next, and a list of member
+                                // variables.
+
+                                 // The public part of the class is rather
+                                 // short: it has a constructor and a function
+                                 // `run' that is called from the outside and
+                                 // acts as something like the `main'
+                                 // function: it coordinates which operations
+                                 // of this class shall be run in which
+                                 // order. Everything else in the class,
+                                 // i.e. all the functions that actually do
+                                 // anything, are in the private section of
+                                 // the class:
 class LaplaceProblem 
 {
   public:
-                                    // This is the constructor:
     LaplaceProblem ();
 
-                                    // And the top-level function,
-                                    // which is called from the
-                                    // outside to start the whole
-                                    // program (see the `main'
-                                    // function at the bottom of this
-                                    // file):
     void run ();
     
-                                    // Then there are some member
-                                    // functions that mostly do what
-                                    // their names suggest. Since
-                                    // they do not need to be called
-                                    // from outside, they are made
+                                    // Then there are the member functions
+                                    // that mostly do what their names
+                                    // suggest. Since they do not need to be
+                                    // called from outside, they are made
                                     // private to this class.
   private:
     void make_grid_and_dofs ();
@@ -124,10 +131,10 @@ class LaplaceProblem
     void solve ();
     void output_results () const;
 
-                                    // And then we have the member
+                                    // And finally we have some member
                                     // variables. There are variables
                                     // describing the triangulation
-                                    // and the numbering of the
+                                    // and the global numbering of the
                                     // degrees of freedom (we will
                                     // specify the exact polynomial
                                     // degree of the finite element
@@ -152,27 +159,35 @@ class LaplaceProblem
     Vector<double>       system_rhs;
 };
 
-
-                                // Here comes the constructor. It
-                                // does not much more than first to
-                                // specify that we want bi-linear
-                                // elements (denoted by the parameter
-                                // to the finite element object,
-                                // which specifies the polynomial
-                                // degree), and to associate the
-                                // dof_handler variable to the
-                                // triangulation we use. All the
-                                // other member variables of the
-                                // LaplaceProblem class have a
-                                // default constructor which does all
-                                // we want.
+                                 // @sect4{LaplaceProblem::LaplaceProblem}
+
+                                // Here comes the constructor. It does not
+                                // much more than first to specify that we
+                                // want bi-linear elements (denoted by the
+                                // parameter to the finite element object,
+                                // which indicates the polynomial degree),
+                                // and to associate the dof_handler variable
+                                // to the triangulation we use. (Note that
+                                // the triangulation isn't set up with a mesh
+                                // at all at the present time, but the
+                                // DoFHandler doesn't care: it only wants to
+                                // know which triangulation it will be
+                                // associated with, and it only starts to
+                                // care about an actual mesh once you try to
+                                // distribute degree of freedom on the mesh
+                                // using the distribute_dofs() function.) All
+                                // the other member variables of the
+                                // LaplaceProblem class have a default
+                                // constructor which does all we want.
 LaplaceProblem::LaplaceProblem () :
                 fe (1),
                dof_handler (triangulation)
 {}
 
 
-                                // Now, the first thing we've got to
+                                 // @sect4{LaplaceProblem::make_grid_and_dofs}
+
+                                 // Now, the first thing we've got to
                                 // do is to generate the
                                 // triangulation on which we would
                                 // like to do our computation and
@@ -194,48 +209,43 @@ void LaplaceProblem::make_grid_and_dofs ()
                                   // total of 1024.
   GridGenerator::hyper_cube (triangulation, -1, 1);
   triangulation.refine_global (5);
-                                  // Unsure that 1024 is the correct
-                                  // number? Let's see:
-                                  // n_active_cells return the number
-                                  // of terminal cells. By terminal
-                                  // we mean the cells on the finest
-                                  // grid.
+                                  // Unsure that 1024 is the correct number?
+                                  // Let's see: n_active_cells return the
+                                  // number of active cells:
   std::cout << "Number of active cells: "
            << triangulation.n_active_cells()
            << std::endl;
-                                  // We stress the adjective
-                                  // `terminal' or `active', since
-                                  // there are more cells, namely the
-                                  // parent cells of the finest
-                                  // cells, their parents, etc, up to
-                                  // the one cell which made up the
-                                  // initial grid. Of course, on the
-                                  // next coarser level, the number
-                                  // of cells is one quarter that of
-                                  // the cells on the finest level,
-                                  // i.e. 256, then 64, 16, 4, and
-                                  // 1. We can get the total number
-                                  // of cells like this:
+                                   // Here, by active we mean the cells on the
+                                  // finest level, i.e. cells that aren't
+                                  // refined any further.  We stress the
+                                  // adjective `active', since there are more
+                                  // cells, namely the parent cells of the
+                                  // finest cells, their parents, etc, up to
+                                  // the one cell which made up the initial
+                                  // grid. Of course, on the next coarser
+                                  // level, the number of cells is one
+                                  // quarter that of the cells on the finest
+                                  // level, i.e. 256, then 64, 16, 4, and
+                                  // 1. We can get the total number of cells
+                                  // like this:
   std::cout << "Total number of cells: "
            << triangulation.n_cells()
            << std::endl;
                                   // Note the distinction between
                                   // n_active_cells() and n_cells().
   
-                                  // Next we enumerate all the
-                                  // degrees of freedom. This is done
-                                  // by using the distribute_dofs
-                                  // function, as we have seen in
-                                  // previous examples. Since we use
-                                  // the FEQ1 class, i.e. bilinear
-                                  // elements, this associates one
-                                  // degree of freedom with each
-                                  // vertex.
+                                  // Next we enumerate all the degrees of
+                                  // freedom. This is done by using the
+                                  // distribute_dofs function, as we have
+                                  // seen in the step-2 example. Since we use
+                                  // the ``FE_Q'' class with a polynomial
+                                  // degree of 1, i.e. bilinear elements,
+                                  // this associates one degree of freedom
+                                  // with each vertex. While we're at
+                                  // generating output, let us also take a
+                                  // look at how many degrees of freedom are
+                                  // generated:
   dof_handler.distribute_dofs (fe);
-
-                                  // Now that we have the degrees of
-                                  // freedom, we can take a look at
-                                  // how many there are:
   std::cout << "Number of degrees of freedom: "
            << dof_handler.n_dofs()
            << std::endl;
@@ -244,13 +254,16 @@ void LaplaceProblem::make_grid_and_dofs ()
                                   // 32 grid, the number of DoFs
                                   // should be 33 times 33, or 1089.
 
-                                  // As we have seen in the previous
-                                  // example, we set up a sparse
-                                  // matrix for the system matrix and
-                                  // tag those entries that might be
-                                  // nonzero. Since that has already
-                                  // been done, we won't discuss the
-                                  // next few lines:
+                                  // As we have seen in the previous example,
+                                  // we set up a sparsity pattern for the
+                                  // system matrix and tag those entries that
+                                  // might be nonzero. Compared to what we
+                                  // did in step-2, the only change is that
+                                  // instead of giving a magically obtained
+                                  // maximal number of nonzero entries per
+                                  // row, we now use a function in the
+                                  // ``DoFHandler'' class that can compute
+                                  // this number for us:
   sparsity_pattern.reinit (dof_handler.n_dofs(),
                           dof_handler.n_dofs(),
                           dof_handler.max_couplings_between_dofs());
@@ -262,7 +275,7 @@ void LaplaceProblem::make_grid_and_dofs ()
                                   // `compress' has been called, you
                                   // can't add nonzero entries
                                   // anymore; the sparsity pattern is
-                                  // `sealed', so to say), and we can
+                                  // `sealed', so to say), we can
                                   // initialize the matrix itself
                                   // with it. Note that the
                                   // SparsityPattern object does
@@ -296,6 +309,8 @@ void LaplaceProblem::make_grid_and_dofs ()
   system_rhs.reinit (dof_handler.n_dofs());
 }
 
+                                 // @sect4{LaplaceProblem::assemble_system}
+
 
                                 // Now comes the difficult part:
                                 // assembling matrices and
@@ -306,132 +321,178 @@ void LaplaceProblem::make_grid_and_dofs ()
                                 // other things in the functions
                                 // above and below.
                                 //
-                                // The general way to assemble
-                                // matrices and vectors is to loop
-                                // over all cells, and on each cell
-                                // compute the contribution of that
-                                // cell to the global matrix and
-                                // right hand side by quadrature. The
-                                // idea now is that since we only
-                                // need the finite element shape
-                                // functions on the quadrature points
-                                // of each cell, we don't need the
-                                // shape functions of the finite
-                                // element themselves any
-                                // more. Therefore, we won't deal
-                                // with the finite element object
-                                // `fe' (which was of type FEQ1), but
-                                // with another object which only
-                                // provides us with the values,
-                                // gradients, etc of the shape
-                                // functions at the quadrature
-                                // points. The objects which do this
-                                // are of type FEValues.
+                                // The general way to assemble matrices and
+                                // vectors is to loop over all cells, and on
+                                // each cell compute the contribution of that
+                                // cell to the global matrix and right hand
+                                // side by quadrature. The point to realize
+                                // now is that we need the values of the
+                                // shape functions at the locations of
+                                // quadrature points on the real
+                                // cell. However, both the finite element
+                                // shape functions as well as the quadrature
+                                // points are only defined on the unit
+                                // cell. They are therefore of little help to
+                                // us, and we will in fact hardly ever query
+                                // information about finite element shape
+                                // functions or quadrature points from these
+                                // objects directly.
+                                //
+                                // Rather, what is required is a way to map
+                                // this data from the unit cell to the real
+                                // cell. Classes that can do that are derived
+                                // from the Mapping class, though one again
+                                // often does not have to deal with them
+                                // directly: many functions in the library
+                                // can take a mapping object as argument, but
+                                // when it is omitted they simply resort to
+                                // the standard bilinear Q1 mapping. We will
+                                // go this route, and not bother with it for
+                                // the moment (we come back to this in
+                                // step-10, step-11, and step-12).
+                                //
+                                // So what we now have is a collection of
+                                // three classes to deal with: finite
+                                // element, quadrature, and mapping
+                                // objects. That's too much, so there is one
+                                // type of class that orchestrates
+                                // information exchange between these three:
+                                // the ``FEValues'' class. If given one
+                                // instance of each three of these objects,
+                                // it will be able to provide you with
+                                // information about values and gradients of
+                                // shape functions at quadrature points on a
+                                // real cell.
+                                 //
+                                 // Using all this, we will assemble the
+                                 // linear system for this problem in the
+                                 // following function:
 void LaplaceProblem::assemble_system () 
 {
-                                  // Ok, let's start: we need a
-                                  // quadrature formula for the
-                                  // evaluation of the integrals on
-                                  // each cell. Let's take a Gauss
-                                  // formula with two quadrature
-                                  // points in each direction, i.e. a
-                                  // total of four points since we
-                                  // are in 2D. This quadrature
-                                  // formula integrates polynomials
-                                  // of degrees up to three exactly
-                                  // (in 1D). Since the integrands in
-                                  // the matrix entries are quadratic
-                                  // (in 1D), this is sufficient. The
-                                  // same holds for 2D.
+                                  // Ok, let's start: we need a quadrature
+                                  // formula for the evaluation of the
+                                  // integrals on each cell. Let's take a
+                                  // Gauss formula with two quadrature points
+                                  // in each direction, i.e. a total of four
+                                  // points since we are in 2D. This
+                                  // quadrature formula integrates
+                                  // polynomials of degrees up to three
+                                  // exactly (in 1D). It is easy to check
+                                  // that this is sufficient for the present
+                                  // problem:
   QGauss<2>  quadrature_formula(2);
-                                  // And we initialize the object
-                                  // which we have briefly talked
-                                  // about above. It needs to be told
-                                  // which the finite element is that
-                                  // we want to use, the quadrature
-                                  // points and their
-                                  // weights. Finally, we have to
-                                  // tell it what we want it to
-                                  // compute on each cell: we need
-                                  // the values of the shape
-                                  // functions at the quadrature
-                                  // points, their gradients, and
-                                  // also the weights of the
-                                  // quadrature points and the
+                                  // And we initialize the object which we
+                                  // have briefly talked about above. It
+                                  // needs to be told which finite element we
+                                  // want to use, and the quadrature points
+                                  // and their weights (jointly described by
+                                  // a Quadrature object). As mentioned, we
+                                  // use the implied Q1 mapping, rather than
+                                  // specifying one ourselves
+                                  // explicitly. Finally, we have to tell it
+                                  // what we want it to compute on each cell:
+                                  // we need the values of the shape
+                                  // functions at the quadrature points (for
+                                  // the right hand side (f,phi)), their
+                                  // gradients (for the matrix entries (grad
+                                  // phi_i, grad phi_j)), and also the
+                                  // weights of the quadrature points and the
                                   // determinants of the Jacobian
-                                  // transformations from the unit
-                                  // cell to the real cells. The
-                                  // values of the shape functions
-                                  // computed by specifying
-                                  // update_values; the gradients are
-                                  // done alike, using
-                                  // update_gradients. The
-                                  // determinants of the Jacobians
-                                  // and the weights are always used
-                                  // together, so only the products
-                                  // (Jacobians times weights, or
-                                  // short JxW) are computed; since
-                                  // we also need them, we have to
-                                  // list them as well. The advantage
-                                  // of this proceeding is that we
-                                  // calculate only what we
-                                  // need. This optimatizes the
-                                  // process of solving:
+                                  // transformations from the unit cell to
+                                  // the real cells.
+                                  //
+                                  // This list of what kind of information we
+                                  // actually need is given as a bitwise
+                                  // connection of flags as the third
+                                  // argument to the constructor of
+                                  // ``FEValues''. Since these values have to
+                                  // be recomputed, or updated, every time we
+                                  // go to a new cell, all of these flags
+                                  // start with the prefix ``update_'' and
+                                  // then indicate what it actually is that
+                                  // we want updated. The flag to give if we
+                                  // want the values of the shape functions
+                                  // computed is ``update_values''; for the
+                                  // gradients it is
+                                  // ``update_gradients''. The determinants
+                                  // of the Jacobians and the quadrature
+                                  // weights are always used together, so
+                                  // only the products (Jacobians times
+                                  // weights, or short ``JxW'') are computed;
+                                  // since we need them, we have to list
+                                  // ``update_JxW_values'' as well:
   FEValues<2> fe_values (fe, quadrature_formula, 
-                        UpdateFlags(update_values    |
-                                    update_gradients |
-                                    update_JxW_values));
-
-                                  // For use further down below, we
-                                  // define two short cuts for the
-                                  // number of degrees of freedom on
-                                  // each cell (since we are in 2D
-                                  // and degrees of freedom are
-                                  // associated with vertices only,
-                                  // this number is four). We also
-                                  // define an abbreviation for the
-                                  // number of quadrature points
-                                  // (here that should be nine). In
-                                  // general, it is a good idea to
-                                  // use their symbolic names instead
-                                  // of hard-coding these number even
-                                  // if you know them, since you may
-                                  // want to change the quadrature
-                                  // formula and/or finite element at
-                                  // some time; the program will just
-                                  // work with these changes, without
-                                  // the need to change the matrix
-                                  // assemblage.
+                        update_values | update_gradients | update_JxW_values);
+                                   // The advantage of this proceeding is that
+                                  // we can specify what kind of information
+                                  // we actually need on each cell. It is
+                                  // easily understandable that this approach
+                                  // can significant speed up finite element
+                                  // computations, compared to approaches
+                                  // where everything, including second
+                                  // derivatives, normal vectors to cells,
+                                  // etc are computed on each cell,
+                                  // regardless whether they are needed or
+                                  // not.
+  
+                                  // For use further down below, we define
+                                  // two short cuts for values that will be
+                                  // used very frequently. First, an
+                                  // abbreviation for the number of degrees
+                                  // of freedom on each cell (since we are in
+                                  // 2D and degrees of freedom are associated
+                                  // with vertices only, this number is four,
+                                  // but we rather want to write the
+                                  // definition of this variable in a way
+                                  // that does not preclude us from later
+                                  // choosing a different finite element that
+                                  // has a different number of degrees of
+                                  // freedom per cell, or work in a different
+                                  // space dimension).
+                                  //
+                                  // Secondly, we also define an abbreviation
+                                  // for the number of quadrature points
+                                  // (here that should be four). In general,
+                                  // it is a good idea to use their symbolic
+                                  // names instead of hard-coding these
+                                  // number even if you know them, since you
+                                  // may want to change the quadrature
+                                  // formula and/or finite element at some
+                                  // time; the program will just work with
+                                  // these changes, without the need to
+                                  // change anything in this function.
                                   //
-                                  // The shortcuts, finally, are only
-                                  // defined to make the following
-                                  // loops a bit more readable. You
-                                  // will see them in many places in
-                                  // larger programs, and
-                                  // `dofs_per_cell' and `n_q_points'
-                                  // are more or less standard names
-                                  // for these purposes.
+                                  // The shortcuts, finally, are only defined
+                                  // to make the following loops a bit more
+                                  // readable. You will see them in many
+                                  // places in larger programs, and
+                                  // `dofs_per_cell' and `n_q_points' are
+                                  // more or less by convention the standard
+                                  // names for these purposes:
   const unsigned int   dofs_per_cell = fe.dofs_per_cell;
   const unsigned int   n_q_points    = quadrature_formula.n_quadrature_points;
 
-                                  // Now, we said that we wanted to
-                                  // assemble the global matrix and
-                                  // vector cell-by-cell. We could
-                                  // write the results directly into
-                                  // the global matrix, but this is
-                                  // not very efficient since access
-                                  // to the elements of a sparse
-                                  // matrix is slow. Rather, we first
-                                  // compute the contribution of each
-                                  // ell in a small matrix with the
-                                  // degrees of freedom on the
-                                  // present cell, and only transfer
-                                  // them to the global matrix when
-                                  // the computations are finished
-                                  // for this cell. We do the same
-                                  // for the right hand side vector,
-                                  // although access times are not so
-                                  // problematic for them.
+                                  // Now, we said that we wanted to assemble
+                                  // the global matrix and vector
+                                  // cell-by-cell. We could write the results
+                                  // directly into the global matrix, but
+                                  // this is not very efficient since access
+                                  // to the elements of a sparse matrix is
+                                  // slow. Rather, we first compute the
+                                  // contribution of each cell in a small
+                                  // matrix with the degrees of freedom on
+                                  // the present cell, and only transfer them
+                                  // to the global matrix when the
+                                  // computations are finished for this
+                                  // cell. We do the same for the right hand
+                                  // side vector. So let's first allocate
+                                  // these objects (these being local
+                                  // objects, all degrees of freedom are
+                                  // coupling with all others, and we should
+                                  // use a full matrix object rather than a
+                                  // sparse one for the local operations;
+                                  // everything will be transferred to a
+                                  // global sparse matrix later on):
   FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
   Vector<double>       cell_rhs (dofs_per_cell);
 
@@ -445,72 +506,70 @@ void LaplaceProblem::assemble_system ()
                                   // we transfer the result into the
                                   // global matrix, we have to know
                                   // the global numbers of the
-                                  // degrees of freedom. When we get
+                                  // degrees of freedom. When we query
                                   // them, we need a scratch
                                   // (temporary) array for these
                                   // numbers:
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
-                                  // Now for the loop over all
-                                  // cells. You have seen before how
-                                  // this works, so this should be
-                                  // familiar to you:
-  DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(),
-                                     endc = dof_handler.end();
+                                  // Now for the loop over all cells. We have
+                                  // seen before how this works, so this
+                                  // should be familiar including the
+                                  // conventional names for these variables:
+  DoFHandler<2>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
   for (; cell!=endc; ++cell)
     {
-                                      // We are on one cell, and we
-                                      // would like the values and
-                                      // gradients of the shape
-                                      // functions be computed, as
-                                      // well as the determinants of
-                                      // the Jacobian matrices of the
-                                      // mapping between unit cell
-                                      // and true cell, at the
-                                      // quadrature points. Since all
-                                      // these values depend on the
-                                      // geometry of the cell, we
-                                      // have to have the FEValues
-                                      // object re-compute them on
-                                      // each cell:
+                                      // We are now sitting on one cell, and
+                                      // we would like the values and
+                                      // gradients of the shape functions be
+                                      // computed, as well as the
+                                      // determinants of the Jacobian
+                                      // matrices of the mapping between unit
+                                      // cell and true cell, at the
+                                      // quadrature points. Since all these
+                                      // values depend on the geometry of the
+                                      // cell, we have to have the FEValues
+                                      // object re-compute them on each cell:
       fe_values.reinit (cell);
 
-                                      // Reset the values of the
-                                      // contributions of this cell
-                                      // to global matrix and global
-                                      // right hand side to zero,
-                                      // before we fill them.
+                                      // Next, reset the local cell's
+                                      // contributions contributions to
+                                      // global matrix and global right hand
+                                      // side to zero, before we fill them:
       cell_matrix = 0;
       cell_rhs = 0;
 
-                                      // Assemble the matrix: For the
-                                      // Laplace problem, the matrix
-                                      // on each cell is the integral
-                                      // over the gradients of shape
-                                      // function i and j. Since we
-                                      // do not integrate, but rather
-                                      // use quadrature, this is the
-                                      // sum over all quadrature
-                                      // points of the integrands
-                                      // times the determinant of the
-                                      // Jacobian matrix at the
-                                      // quadrature point times the
-                                      // weight of this quadrature
-                                      // point. You can get the
-                                      // gradient of shape function i
-                                      // at quadrature point q_point
-                                      // by using
+                                      // Then finally assemble the matrix:
+                                      // For the Laplace problem, the matrix
+                                      // on each cell is the integral over
+                                      // the gradients of shape function i
+                                      // and j. Since we do not integrate,
+                                      // but rather use quadrature, this is
+                                      // the sum over all quadrature points
+                                      // of the integrands times the
+                                      // determinant of the Jacobian matrix
+                                      // at the quadrature point times the
+                                      // weight of this quadrature point. You
+                                      // can get the gradient of shape
+                                      // function i at quadrature point
+                                      // q_point by using
                                       // fe_values.shape_grad(i,q_point);
-                                      // this gradient is a
-                                      // 2-dimensional vector (in
-                                      // fact it is of type
-                                      // Tensor<1,dim>, with here
-                                      // dim=2) and the product of
-                                      // two such vectors is the
-                                      // scalar product, i.e. the
-                                      // product of the two
-                                      // shape_grad function calls is
-                                      // the dot product.
+                                      // this gradient is a 2-dimensional
+                                      // vector (in fact it is of type
+                                      // Tensor<1,dim>, with here dim=2) and
+                                      // the product of two such vectors is
+                                      // the scalar product, i.e. the product
+                                      // of the two shape_grad function calls
+                                      // is the dot product. This is in turn
+                                      // multiplied by the Jacobian
+                                      // determinant and the quadrature point
+                                      // weight (that one gets together by
+                                      // the call to
+                                      // ``fe_values.JxW''). Finally, this is
+                                      // repeated for all shape functions
+                                      // phi_i and phi_j:
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        for (unsigned int j=0; j<dofs_per_cell; ++j)
          for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
@@ -575,86 +634,73 @@ void LaplaceProblem::assemble_system ()
                                       // side vector.
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        system_rhs(local_dof_indices[i]) += cell_rhs(i);
-    };
-
-
-                                  // Now almost everything is set up
-                                  // for the solution of the discrete
-                                  // system. However, we have not yet
-                                  // taken care of boundary values
-                                  // (in fact, Laplace's equation
-                                  // without Dirichlet boundary
-                                  // values is not even uniquely
-                                  // solvable, since you can add an
-                                  // arbitrary constant to the
-                                  // discrete solution). We therefore
-                                  // have to take into account
-                                  // boundary values.
+    }
+
+
+                                  // Now almost everything is set up for the
+                                  // solution of the discrete
+                                  // system. However, we have not yet taken
+                                  // care of boundary values (in fact,
+                                  // Laplace's equation without Dirichlet
+                                  // boundary values is not even uniquely
+                                  // solvable, since you can add an arbitrary
+                                  // constant to the discrete solution). We
+                                  // therefore have to do something about the
+                                  // situation.
                                   //
-                                  // For this, we first obtain a list
-                                  // of the degrees of freedom on the
-                                  // boundary and the value the shape
-                                  // function shall have there. For
-                                  // simplicity, we only interpolate
-                                  // the boundary value function,
-                                  // rather than projecting them onto
-                                  // the boundary. There is a
-                                  // function in the library which
-                                  // does exactly this:
-                                  // interpolate_boundary_values. Its
-                                  // parameters are (omitting
-                                  // parameters for which default
-                                  // values exist which are
-                                  // sufficient here): the DoFHandler
-                                  // object to get the global numbers
-                                  // of the degrees of freedom on the
-                                  // boundary; the component of the
-                                  // boundary where the boundary
-                                  // values shall be interpolated;
-                                  // the boundary value function
-                                  // itself; and the output object.
+                                  // For this, we first obtain a list of the
+                                  // degrees of freedom on the boundary and
+                                  // the value the shape function shall have
+                                  // there. For simplicity, we only
+                                  // interpolate the boundary value function,
+                                  // rather than projecting it onto the
+                                  // boundary. There is a function in the
+                                  // library which does exactly this:
+                                  // ``VectorTools::interpolate_boundary_values''. Its
+                                  // parameters are (omitting parameters for
+                                  // which default values exist and that we
+                                  // don't care about): the DoFHandler object
+                                  // to get the global numbers of the degrees
+                                  // of freedom on the boundary; the
+                                  // component of the boundary where the
+                                  // boundary values shall be interpolated;
+                                  // the boundary value function itself; and
+                                  // the output object.
                                   //
-                                  // The component of the boundary is
-                                  // meant as follows: in many cases,
-                                  // you may want to impose certain
-                                  // boundary values only on parts of
-                                  // the boundary. For example, you
-                                  // may have inflow and outflow
-                                  // boundaries in fluid dynamics,
-                                  // are clamped and free parts of
-                                  // bodies in deformation
-                                  // computations of bodies. Then you
-                                  // will want to denote these
-                                  // different parts of the boundary
-                                  // by different numbers and tell
-                                  // the interpolate_boundary_values
-                                  // function to only compute the
-                                  // boundary values on a certain
-                                  // part of the boundary (e.g. the
-                                  // clamped part, or the inflow
-                                  // boundary). By default, all
-                                  // boundaries have the number `0',
-                                  // and since we have not changed
-                                  // that, this is still so;
-                                  // therefore, if we give `0' as the
-                                  // desired portion of the boundary,
-                                  // this means we get the whole
-                                  // boundary. If you got different
-                                  // boundaries, you have to number
-                                  // them differently and have to be
-                                  // attentive when using the method
-                                  // above.
+                                  // The component of the boundary is meant
+                                  // as follows: in many cases, you may want
+                                  // to impose certain boundary values only
+                                  // on parts of the boundary. For example,
+                                  // you may have inflow and outflow
+                                  // boundaries in fluid dynamics, or clamped
+                                  // and free parts of bodies in deformation
+                                  // computations of bodies. Then you will
+                                  // want to denote these different parts of
+                                  // the boundary by different numbers and
+                                  // tell the interpolate_boundary_values
+                                  // function to only compute the boundary
+                                  // values on a certain part of the boundary
+                                  // (e.g. the clamped part, or the inflow
+                                  // boundary). By default, all boundaries
+                                  // have the number `0', and since we have
+                                  // not changed that, this is still so;
+                                  // therefore, if we give `0' as the desired
+                                  // portion of the boundary, this means we
+                                  // get the whole boundary. If you have
+                                  // boundaries with kinds of boundaries, you
+                                  // have to number them differently. The
+                                  // function call below will then only
+                                  // determine boundary values for parts of
+                                  // the boundary.
                                   //
-                                  // The function describing the
-                                  // boundary values is an object of
-                                  // type `Function' or of a derived
-                                  // class. One of the derived
-                                  // classes is ZeroFunction, which
-                                  // described a function which is
-                                  // zero everywhere. We create such
-                                  // an object in-place and pass it
-                                  // to the
-                                  // interpolate_boundary_values
+                                  // The function describing the boundary
+                                  // values is an object of type ``Function''
+                                  // or of a derived class. One of the
+                                  // derived classes is ``ZeroFunction'',
+                                  // which describes (not unexpectedly) a
+                                  // function which is zero everywhere. We
+                                  // create such an object in-place and pass
+                                  // it to the interpolate_boundary_values
                                   // function.
                                   //
                                   // Finally, the output object is a
@@ -666,7 +712,7 @@ void LaplaceProblem::assemble_system ()
                                   // here for all entries). This
                                   // mapping of DoF numbers to
                                   // boundary values is done by the
-                                  // `map' class.
+                                  // ``std::map'' class.
   std::map<unsigned int,double> boundary_values;
   VectorTools::interpolate_boundary_values (dof_handler,
                                            0,
@@ -686,7 +732,9 @@ void LaplaceProblem::assemble_system ()
 }
 
 
-                                // The following function simply
+                                 // @sect4{LaplaceProblem::solve}
+
+                                 // The following function simply
                                 // solves the discretized
                                 // equation. As the system is quite a
                                 // large one for direct solvers such
@@ -697,42 +745,42 @@ void LaplaceProblem::assemble_system ()
                                 // variables here (only 1089) is a
                                 // very small number for finite
                                 // element computations, where
-                                // 100.000 is a more usual number;
-                                // for this number of variables,
+                                // 100.000 is a more usual number.
+                                // For this number of variables,
                                 // direct methods are no longer
                                 // usable and you are forced to use
                                 // methods like CG.
 void LaplaceProblem::solve () 
 {
-                                  // We need to tell the algorithm
-                                  // where to stop. This is done by
-                                  // using a SolverControl object,
-                                  // and as stopping criterion we
-                                  // say: maximally 1000 iterations
-                                  // (which is far more than is
-                                  // needed for 1089 variables; see
-                                  // the results section to find out
-                                  // how many were really used), and
-                                  // stop if the norm of the residual
-                                  // is below 1e-12. In practice, the
-                                  // latter criterion will be the one
-                                  // which stops the iteration.
+                                  // First, we need to have an object that
+                                  // knows how to tell the CG algorithm when
+                                  // to stop. This is done by using a
+                                  // ``SolverControl'' object, and as
+                                  // stopping criterion we say: stop after a
+                                  // maximum of 1000 iterations (which is far
+                                  // more than is needed for 1089 variables;
+                                  // see the results section to find out how
+                                  // many were really used), and stop if the
+                                  // norm of the residual is below 1e-12. In
+                                  // practice, the latter criterion will be
+                                  // the one which stops the iteration:
   SolverControl           solver_control (1000, 1e-12);
-                                  // Then we need the solver
-                                  // itself. The template parameters
-                                  // here are the matrix type and the
-                                  // type of the vectors, but the
-                                  // empty angle brackets indicate
-                                  // that we simply take the default
-                                  // arguments.
+                                  // Then we need the solver itself. The
+                                  // template parameters to the ``SolverCG''
+                                  // class are the matrix type and the type
+                                  // of the vectors, but the empty angle
+                                  // brackets indicate that we simply take
+                                  // the default arguments (which are
+                                  // ``SparseMatrix<double>'' and
+                                  // ``Vector<double>''):
   SolverCG<>              cg (solver_control);
 
-                                  // Now solve the system of
-                                  // equations. The CG solver takes a
-                                  // preconditioner, but we don't
-                                  // want to use one, so we tell it
-                                  // to use the identity operation as
-                                  // preconditioner.
+                                  // Now solve the system of equations. The
+                                  // CG solver takes a preconditioner as its
+                                  // fourth argument. We don't feel ready to
+                                  // delve into this yet, so we tell it to
+                                  // use the identity operation as
+                                  // preconditioner:
   cg.solve (system_matrix, solution, system_rhs,
            PreconditionIdentity());
                                   // Now that the solver has done its
@@ -742,6 +790,8 @@ void LaplaceProblem::solve ()
 }
 
 
+                                 // @sect4{LaplaceProblem::output_results}
+
                                 // The last part of a typical finite
                                 // element program is to output the
                                 // results and maybe do some
@@ -757,23 +807,20 @@ void LaplaceProblem::output_results () const
                                   // To write the output to a file,
                                   // we need an object which knows
                                   // about output formats and the
-                                  // like. This is the DataOut class,
+                                  // like. This is the ``DataOut'' class,
                                   // and we need an object of that
                                   // type:
   DataOut<2> data_out;
-                                  // Now we have to tell it where to
-                                  // take the values from which it
-                                  // shall write. We tell it which
-                                  // DoFHandler object to use, and we
-                                  // add the solution vector (and the
-                                  // name by which it shall be
-                                  // written to disk) to the list of
-                                  // data that is to be written. If
-                                  // we had more than one vector
-                                  // which we would like to look at
-                                  // in the output (for example right
-                                  // hand sides, errors per cell,
-                                  // etc) we would add them as well:
+                                  // Now we have to tell it where to take the
+                                  // values from which it shall write. We
+                                  // tell it which ``DoFHandler'' object to
+                                  // use, and the solution vector (and
+                                  // the name by which the solution variable
+                                  // shall appear in the output file). If
+                                  // we had more than one vector which we
+                                  // would like to look at in the output (for
+                                  // example right hand sides, errors per
+                                  // cell, etc) we would add them as well:
   data_out.attach_dof_handler (dof_handler);
   data_out.add_data_vector (solution, "solution");
                                   // After the DataOut object knows
@@ -783,9 +830,9 @@ void LaplaceProblem::output_results () const
                                   // handle. The reason is that we
                                   // have separated the frontend
                                   // (which knows about how to treat
-                                  // DoFHandler objects and data
+                                  // ``DoFHandler'' objects and data
                                   // vectors) from the back end (which
-                                  // knows several output formats)
+                                  // knows many different output formats)
                                   // and use an intermediate data
                                   // format to transfer data from the
                                   // front- to the backend. The data
@@ -806,15 +853,16 @@ void LaplaceProblem::output_results () const
 }
 
 
-                                // The following function is the main
-                                // function which calls all the other
-                                // functions of the LaplaceProblem
-                                // class. The order in which this is
-                                // done resembles the order in which
-                                // most finite element programs
-                                // work. Since the names are mostly
-                                // self-explanatory, there is not
-                                // much to comment about:
+                                 // @sect4{LaplaceProblem::run}
+
+                                // Finally, the last function of this class
+                                // is the main function which calls all the
+                                // other functions of the ``LaplaceProblem''
+                                // class. The order in which this is done
+                                // resembles the order in which most finite
+                                // element programs work. Since the names are
+                                // mostly self-explanatory, there is not much
+                                // to comment about:
 void LaplaceProblem::run () 
 {
   make_grid_and_dofs ();
@@ -823,7 +871,8 @@ void LaplaceProblem::run ()
   output_results ();
 }
 
-    
+
+                                 // @sect3{The ``main'' function}
 
                                 // This is the main function of the
                                 // program. Since the concept of a
@@ -832,11 +881,12 @@ void LaplaceProblem::run ()
                                 // programming, it often does not
                                 // much more than creating an object
                                 // of the top-level class and calling
-                                // it principle function. This is
-                                // what is done here as well.
+                                // its principle function. This is
+                                // what is done here as well:
 int main () 
 {
   LaplaceProblem laplace_problem;
   laplace_problem.run ();
+
   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.