]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Rewrite parts of the documentation.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 25 Jan 2006 05:04:36 +0000 (05:04 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 25 Jan 2006 05:04:36 +0000 (05:04 +0000)
git-svn-id: https://svn.dealii.org/trunk@12156 0785d39b-7218-0410-832d-ea1e28bc413d

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

index bfef563132d4cdf8c56d68eee764f997a9fec143..0eb5c4a5e32c684a842afa8d4b7a88ec18b63740 100644 (file)
@@ -4,26 +4,26 @@
 /*    $Id$       */
 /*    Version: $Name$                                          */
 /*                                                                */
-/*    Copyright (C) 1999, 2000, 2001, 2002, 2003 by the deal.II authors */
+/*    Copyright (C) 1999, 2000, 2001, 2002, 2003, 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.                        */
 
-                                // The following includes are just
-                                // like for the previous program, so
-                                // will not be commented further
+                                // The first few includes are just
+                                // like in the previous program, so
+                                // do not require additional comments:
 #include <grid/tria.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
 #include <grid/grid_generator.h>
 #include <grid/tria_boundary_lib.h>
 
-                                // We need this include file for the
-                                // association of degrees of freedom
-                                // ("DoF"s) to vertices, lines, and
-                                // cells.
+                                // However, the next file is new. We need
+                                // this include file for the association of
+                                // degrees of freedom ("DoF"s) to vertices,
+                                // lines, and cells:
 #include <dofs/dof_handler.h>
 
                                 // The following include contains the
                                 // but none on faces and none in the
                                 // interior of the cells.
                                 //
-                                // In fact, the file contains the
-                                // description of several more finite
-                                // elements as well, such as
-                                // biquadratic, bicubic and biquartic
-                                // elements, but not only for two
-                                // space dimensions, but also for one
-                                // and three dimensions.
+                                // (In fact, the file contains the
+                                // description of Lagrange elements in
+                                // general, i.e. also the quadratic, cubic,
+                                // etc versions, and not only for 2d but also
+                                // 1d and 3d.)
 #include <fe/fe_q.h>
                                 // In the following file, several
                                 // tools for manipulating degrees of
                                 // freedom. It is declared here:
 #include <dofs/dof_renumbering.h>
 
-                                // This is needed for C++ output:
+                                // And this is again needed for C++ output:
 #include <fstream>
 
 
+                                 // @sect3{Mesh generation}
 
-                                // This is the function that produced
-                                // the circular grid in the previous
-                                // example. The sole difference is
-                                // that it returns the grid it
-                                // produces via its argument.
+                                // This is the function that produced the
+                                // circular grid in the previous step-1
+                                // example program. The sole difference is
+                                // that it returns the grid it produces via
+                                // its argument.
                                 //
-                                // We won't comment on the internals
-                                // of this function, since this has
-                                // been done in the previous
-                                // example. If you don't understand
-                                // what is happening here, look
-                                // there.
+                                // The details of what the function does are
+                                // explained in step-1. The only thing we
+                                // would like to comment on is this:
+                                // 
+                                 // Since we want to export the triangulation
+                                 // through this function's parameter, we need
+                                 // to make sure that the boundary object
+                                 // lives at least as long as the
+                                 // triangulation does. However, in step-1,
+                                 // the boundary object is a local variable,
+                                 // and it would be deleted at the end of the
+                                 // function, which is too early. We avoid the
+                                 // problem by declaring it 'static' which
+                                 // makes sure that the object is initialized
+                                 // the first time control flow passes its
+                                 // point of declaration, but at the same time
+                                 // assures that it lives until the end of the
+                                 // program.
 void make_grid (Triangulation<2> &triangulation)
 {
   const Point<2> center (1,0);
   const double inner_radius = 0.5,
               outer_radius = 1.0;
   GridGenerator::hyper_shell (triangulation,
-                             center, inner_radius, outer_radius, 10);
+                             center, inner_radius, outer_radius,
+                              10);
 
-                                  // This is the single difference to
-                                  // the respective function in the
-                                  // previous program: since we want
-                                  // to export the triangulation
-                                  // through this function's
-                                  // parameter, we need to make sure
-                                  // that the boundary object lives
-                                  // at least as long as the
-                                  // triangulation does. However,
-                                  // since the boundary object is a
-                                  // local variable, it would be
-                                  // deleted at the end of this
-                                  // function, which is too early; by
-                                  // declaring it 'static', we can
-                                  // assure that it lives until the
-                                  // end of the program.
   static const HyperShellBoundary<2> boundary_description(center);
   triangulation.set_boundary (0, boundary_description);
   
   for (unsigned int step=0; step<5; ++step)
     {
-      Triangulation<2>::active_cell_iterator cell, endc;
-      cell = triangulation.begin_active();
-      endc = triangulation.end();
+      Triangulation<2>::active_cell_iterator
+       cell = triangulation.begin_active(),
+       endc = triangulation.end();
 
       for (; cell!=endc; ++cell)
        for (unsigned int vertex=0;
             vertex < GeometryInfo<2>::vertices_per_cell;
             ++vertex)
          {
-           const Point<2> vector_to_center
-             = (cell->vertex(vertex) - center);
-           const double distance_from_center
-             = std::sqrt(vector_to_center.square());
+            const double distance_from_center
+              = center.distance (cell->vertex(vertex));
            
            if (std::fabs(distance_from_center - inner_radius) < 1e-10)
              {
                cell->set_refine_flag ();
                break;
-             };
-         };
+             }
+         }
 
       triangulation.execute_coarsening_and_refinement ();
-    };
+    }
 }
 
+                                 // @sect3{Creation of a DoFHandler}
 
-                                // Up to now, we only have a grid,
-                                // i.e. some geometrical (the
-                                // position of the vertices and which
-                                // vertices make up which cell) and
-                                // some topological information
-                                // (neighborhoods of cells). To use
-                                // numerical algorithms, one needs
-                                // some logic information in addition
-                                // to that: we would like to
-                                // associate degree of freedom
-                                // numbers to each vertex (or line,
-                                // or cell, in case we were using
-                                // higher order elements) to later
-                                // generate matrices and vectors
-                                // which describe a finite element
+                                // Up to now, we only have a grid, i.e. some
+                                // geometrical (the position of the vertices)
+                                // and some topological information (how
+                                // vertices are connected to lines, and lines
+                                // to cells, as well as which cells neighbor
+                                // which other cells). To use numerical
+                                // algorithms, one needs some logic
+                                // information in addition to that: we would
+                                // like to associate degree of freedom
+                                // numbers to each vertex (or line, or cell,
+                                // in case we were using higher order
+                                // elements) to later generate matrices and
+                                // vectors which describe a finite element
                                 // field on the triangulation.
+                                 //
+                                 // This function shows how to do this. The
+                                 // object to consider is the ``DoFHandler''
+                                 // class template.  Before we do so, however,
+                                 // we first need something that describes how
+                                 // many degrees of freedom are to be
+                                 // associated to each of these objects. Since
+                                 // this is one aspect of the definition of a
+                                 // finite element space, the finite element
+                                 // base class stores this information. In the
+                                 // present context, we therefore create an
+                                 // object of the derived class ``FE_Q'' that
+                                 // describes Lagrange elements. Its
+                                 // constructor takes one argument that states
+                                 // the polynomial degree of the element,
+                                 // which here is one (indicating a bi-linear
+                                 // element); this then corresponds to one
+                                 // degree of freedom for each vertex, while
+                                 // there are none on lines and inside the
+                                 // quadrilateral. A value of, say, three
+                                 // given to the constructor would instead
+                                 // give us a bi-cubic element with one degree
+                                 // of freedom per vertex, two per line, and
+                                 // four inside the cell. In general, ``FE_Q''
+                                 // denotes the family of continuous elements
+                                 // with complete polynomials
+                                 // (i.e. tensor-product polynomials) up to
+                                 // the specified order.
+                                 //
+                                 // We first need to create an object of this
+                                 // class and then pass it on to the
+                                 // ``DoFHandler'' object to allocate storage
+                                 // for the degrees of freedom (in deal.II
+                                 // lingo: we ``distribute degrees of
+                                 // freedom''). Note that the DoFHandler
+                                 // object will store a reference to this
+                                 // finite element object, so we need have to
+                                 // make sure its lifetime is at least as long
+                                 // as that of the ``DoFHandler''; one way to
+                                 // make sure this is so is to make it static
+                                 // as well, in order to prevent its
+                                 // preemptive destruction. (However, the
+                                 // library would warn us if we forgot about
+                                 // this and abort the program if that
+                                 // occured. You can check this, if you want,
+                                 // by removing the 'static' declaration.)
 void distribute_dofs (DoFHandler<2> &dof_handler) 
 {
-                                  // In order to associate degrees of
-                                  // freedom with features of a
-                                  // triangulation (vertices, lines,
-                                  // quadrilaterals), we need an
-                                  // object which describes how many
-                                  // degrees of freedom are to be
-                                  // associated to each of these
-                                  // objects. This is done using the
-                                  // ``FE_Q'' class. Giving a
-                                   // constructor argument one
-                                  // states that one degree of
-                                  // freedom is to be assigned to
-                                  // each vertex, while there are
-                                  // none on lines and inside the
-                                  // quadrilateral. In fact, the
-                                  // argument denotes the polynomial
-                                  // degree, in this case, we get
-                                  // bilinear finite elements in two
-                                  // space dimensions; a value of,
-                                  // say, three would give us
-                                  // bi-cubic ones. In general, ``FE_Q''
-                                  // denotes the family of continuous
-                                  // elements with complete polynomials
-                                  // (i.e. tensor-product polynomials)
-                                  // up to the specified order 
-                                  //
-                                  // We first need to create an
-                                  // object of this class and use it
-                                  // to distribute the degrees of
-                                  // freedom. Note that the
-                                  // DoFHandler object will store a
-                                  // reference to this object, so we
-                                  // need to make it static as well,
-                                  // in order to prevent its
-                                  // preemptive
-                                  // destruction. (However, the
-                                  // library would warn us about this
-                                  // and exit the program if that
-                                  // occured. You can check this, if
-                                  // you want, by removing the
-                                  // 'static' declaration.)
+                                   // As described above, let us first create
+                                   // a finite element object, and then use it
+                                   // to allocate degrees of freedom on the
+                                   // triangulation with which the dof_handler
+                                   // object is associated:
   static const FE_Q<2> finite_element(1);
   dof_handler.distribute_dofs (finite_element);
 
-                                  // Now we have associated a number
-                                  // to each vertex, but how can we
-                                  // visualize this? Unfortunately,
-                                  // presently there is no way
-                                  // implemented to directly show the
-                                  // DoF number associated with each
-                                  // vertex. However, such
-                                  // information would hardly ever be
-                                  // truly important, since the
-                                  // numbering itself is more or less
-                                  // arbitrary. There are more
-                                  // important factors, of which we
-                                  // will visualize one in the
-                                  // following.
+                                  // Now that we have associated a degree of
+                                  // freedom with a global number to each
+                                  // vertex, we wonder how to visualize this?
+                                  // Unfortunately, presently no way is
+                                  // implemented to directly show the DoF
+                                  // number associated with each
+                                  // vertex. However, such information would
+                                  // hardly ever be truly important, since
+                                  // the numbering itself is more or less
+                                  // arbitrary. There are more important
+                                  // factors, of which we will visualize one
+                                  // in the following.
                                   //
-                                  // Associated with each vertex of
-                                  // the triangulation is a shape
-                                  // function. Assume we want to
-                                  // solve something like Laplace's
-                                  // equation, then the different
-                                  // matrix entries will be the
-                                  // integrals over the gradient of
-                                  // each two such shape
-                                  // functions. Obviously, since the
-                                  // shape functions are not equal to
-                                  // zero only on the cells adjacent
-                                  // to the vertex they are
-                                  // associated to, matrix entries
-                                  // will be nonzero only of the
-                                  // supports of the shape functions
-                                  // associated to the column and row
-                                  // numbers intersect. This is only
-                                  // the case for adjacent shape
-                                  // functions, and therefore only
-                                  // for adjacent vertices. Now,
-                                  // since the vertices are numbered
-                                  // more or less randomly be the
-                                  // above function
-                                  // (distribute_dofs), the pattern
-                                  // of nonzero entries in the matrix
-                                  // will be somewhat ragged, and we
-                                  // will take a look at it now.
+                                  // Associated with each vertex of the
+                                  // triangulation is a shape
+                                  // function. Assume we want to solve
+                                  // something like Laplace's equation, then
+                                  // the different matrix entries will be the
+                                  // integrals over the gradient of each pair
+                                  // of such shape functions. Obviously,
+                                  // since the shape functions are nonzero
+                                  // only on the cells adjacent to the vertex
+                                  // they are associated to, matrix entries
+                                  // will be nonzero only if the supports of
+                                  // the shape functions associated to that
+                                  // column and row numbers intersect. This
+                                  // is only the case for adjacent shape
+                                  // functions, and therefore only for
+                                  // adjacent vertices. Now, since the
+                                  // vertices are numbered more or less
+                                  // randomly by the above function
+                                  // (DoFHandler::distribute_dofs), the
+                                  // pattern of nonzero entries in the matrix
+                                  // will be somewhat ragged, and we will
+                                  // take a look at it now.
                                   //
                                   // First we have to create a
                                   // structure which we use to store
                                   // the places of nonzero
                                   // elements. We have to give it the
                                   // size of the matrix, which in our
-                                  // case will be square with that
+                                  // case will be square with as
                                   // many rows and columns as there
                                   // are degrees of freedom on the
                                   // grid:
   SparsityPattern sparsity_pattern (dof_handler.n_dofs(),
-                                   dof_handler.n_dofs());
-                                  // We fill it with the places where
-                                  // nonzero elements will be located
-                                  // given the present numbering of
+                                   dof_handler.n_dofs(),
+                                    20);
+                                   // The last argument to the constructor
+                                   // indicates the maximum number of entries
+                                   // we expect per row. If this were a
+                                   // uniformly refined square, then we know
+                                   // that each vertex degree of freedom would
+                                   // couple with itself and the eight degrees
+                                   // of freedom around it. However, our mesh
+                                   // is more complicated, and it may well be
+                                   // that more couplings will occur. The
+                                   // value 20 we use here is on the safe side
+                                   // of that, though it may actually be too
+                                   // large. In the step-3 tutorial program,
+                                   // we will see a way to compute a
+                                   // reasonable upper bound to the number of
+                                   // nonzero entries, and later programs will
+                                   // show several other methods to compute
+                                   // the numbers of zeros per row.
+
+                                  // We then fill this object with the
+                                  // places where nonzero elements will be
+                                  // located given the present numbering of
                                   // degrees of freedom:
   DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
                                   // Before further work can be done
@@ -256,8 +271,24 @@ void distribute_dofs (DoFHandler<2> &dof_handler)
                                   // for some internal
                                   // reorganization:
   sparsity_pattern.compress ();
+                                   // What actually happens in this call is
+                                   // the following: upon creation of a
+                                   // ``SparsityPattern'' object, memory is
+                                   // allocated for a maximum number of
+                                   // entries per row (20 in our case). The
+                                   // call to
+                                   // ``DoFTools::make_sparsity_pattern'' then
+                                   // actually allocates entries as necessary,
+                                   // but the number of nonzero entries in any
+                                   // given row may be less than the 20 we
+                                   // have allocated memory for. To save
+                                   // memory and to simplify some other
+                                   // operations, one then needs to
+                                   // ``compress'' the sparsity pattern before
+                                   // anything else.
 
-                                  // Now write the results to a file
+                                  // With this, we can now write the results
+                                  // to a file:
   std::ofstream out ("sparsity_pattern.1");
   sparsity_pattern.print_gnuplot (out);
                                   // The result is in GNUPLOT format,
@@ -266,23 +297,25 @@ void distribute_dofs (DoFHandler<2> &dof_handler)
                                   // nonzero entry are listed. The
                                   // output will be shown below.
                                   //
-                                  // If you look at it, you will note
-                                  // that the sparsity pattern is
-                                  // symmetric, which is quite often
-                                  // so, unless you have a rather
-                                  // special equation you want to
-                                  // solve. You will also note that
-                                  // it has several distinct region,
-                                  // which stem from the fact that
-                                  // the numbering starts from the
-                                  // coarsest cells and moves on to
-                                  // the finer ones; since they are
-                                  // all distributed symmetrically
-                                  // around the origin, this shows up
-                                  // again in the sparsity pattern.
+                                  // If you look at it, you will note that
+                                  // the sparsity pattern is symmetric. This
+                                  // should not come as a surprise, since we
+                                  // have not given the
+                                  // ``DoFTools::make_sparsity_pattern'' any
+                                  // information that would indicate that our
+                                  // bilinear form may couple shape functions
+                                  // in a non-symmetric way. You will also
+                                  // note that it has several distinct
+                                  // region, which stem from the fact that
+                                  // the numbering starts from the coarsest
+                                  // cells and moves on to the finer ones;
+                                  // since they are all distributed
+                                  // symmetrically around the origin, this
+                                  // shows up again in the sparsity pattern.
 }
 
 
+                                 // @sect3{Renumbering of DoFs}
 
                                 // In the sparsity pattern produced
                                 // above, the nonzero entries
@@ -320,65 +353,76 @@ void distribute_dofs (DoFHandler<2> &dof_handler)
                                 // close to the original one. Then,
                                 // their neighbors, if not yet
                                 // numbered, are numbered, and so
-                                // on. One such algorithm is the one
-                                // by Cuthill and McKee, which is a
-                                // little more complicated, but works
-                                // along the same lines. We will use
-                                // it to renumber the degrees of
-                                // freedom such that the resulting
-                                // sparsity pattern is more localized
-                                // around the diagonal.
+                                // on.
+                                //
+                                // One algorithm that adds a little bit of
+                                // sophistication along these lines is the
+                                // one by Cuthill and McKee. We will use it
+                                // in the following function to renumber the
+                                // degrees of freedom such that the resulting
+                                // sparsity pattern is more localized around
+                                // the diagonal. The only interesting part of
+                                // the function is the first call to
+                                // ``DoFRenumbering::Cuthill_McKee'', the
+                                // rest is essentially as before:
 void renumber_dofs (DoFHandler<2> &dof_handler) 
 {
-                                  // Renumber the degrees of freedom...
   DoFRenumbering::Cuthill_McKee (dof_handler);
-                                  // ...regenerate the sparsity pattern...
   SparsityPattern sparsity_pattern (dof_handler.n_dofs(),
                                    dof_handler.n_dofs());
+
   DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
   sparsity_pattern.compress ();
-                                  // ...and output the result:
+
   std::ofstream out ("sparsity_pattern.2");
   sparsity_pattern.print_gnuplot (out);
-                                  // Again, the output is shown
-                                  // below. Note that the nonzero
-                                  // entries are clustered far better
-                                  // around the diagonal than
-                                  // before. This effect is even more
-                                  // distinguished for larger
-                                  // matrices (the present one has
-                                  // 1260 rows and columns, but large
-                                  // matrices often have several
-                                  // 100,000s).
 }
 
+                                 // Again, the output is shown
+                                 // below. Note that the nonzero
+                                 // entries are clustered far better
+                                 // around the diagonal than
+                                 // before. This effect is even more
+                                 // distinguished for larger
+                                 // matrices (the present one has
+                                 // 1260 rows and columns, but large
+                                 // matrices often have several
+                                 // 100,000s).
+
+                                 // It is worth noting that the
+                                 // ``DoFRenumbering'' class offers a number
+                                 // of other algorithms as well to renumber
+                                 // degrees of freedom. For example, it would
+                                 // of course be ideal if all couplings were
+                                 // in the lower or upper triangular part of a
+                                 // matrix, since then solving the linear
+                                 // system would amoung to only forward or
+                                 // backward substitution. This is of course
+                                 // unachievable for symmetric sparsity
+                                 // patterns, but in some special situations
+                                 // involving transport equations, this is
+                                 // possible by enumerating degrees of freedom
+                                 // from the inflow boundary along streamlines
+                                 // to the outflow boundary. Not surprisingly,
+                                 // ``DoFRenumbering'' also has algorithms for
+                                 // this.
 
 
+                                 // @sect3{The main function}
 
-                                // This is the main program, which
-                                // only calls the other functions in
-                                // their respective order.
+                                // Finally, this is the main program. The
+                                // only thing it does is to allocate and
+                                // create the triangulation, then create a
+                                // ``DoFHandler'' object and associate it to
+                                // the triangulation, and finally call above
+                                // two functions on it:
 int main () 
 {
-                                  // Allocate space for a triangulation...
   Triangulation<2> triangulation;
-                                  // ...and create it
   make_grid (triangulation);
 
-                                  // A variable that will hold the
-                                  // information which vertex has
-                                  // which number. The geometric
-                                  // information is passed as
-                                  // parameter and a pointer to the
-                                  // triangulation will be stored
-                                  // inside the DoFHandler object.
   DoFHandler<2> dof_handler (triangulation);
-                                  // Associate vertices and degrees
-                                  // of freedom.
-  distribute_dofs (dof_handler);
 
-                                  // Show the effect of renumbering
-                                  // of degrees of freedom to the
-                                  // sparsity pattern of the matrix.
+  distribute_dofs (dof_handler);
   renumber_dofs (dof_handler);
 }

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.