]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
documentation
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 17 Jun 2010 04:01:18 +0000 (04:01 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 17 Jun 2010 04:01:18 +0000 (04:01 +0000)
git-svn-id: https://svn.dealii.org/trunk@21220 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/multigrid/multigrid.h
deal.II/examples/step-39/doc/results.dox
deal.II/examples/step-39/step-39.cc

index 9ffeddc4b0af45b5616900ab7ceaa1c8db6cc011..f48b57a0e14d441d5eeb795192b5bcc90e47abfd 100644 (file)
@@ -254,7 +254,7 @@ class Multigrid : public Subscriptor
                                      * matrices originate from
                                      * discontinuous Galerkin methods
                                      * (see FE_DGQ etc.), where they
-                                     * correspond tu the edge fluxes
+                                     * correspond to the edge fluxes
                                      * at the refinement edge between
                                      * two levels.
                                      *
index 91c515a668c83532ac4dd6f9c2f26d189edcf7f3..20b3dc462d52bf894cc6dcaeaffaada1340ad71e 100644 (file)
@@ -14,7 +14,7 @@ DEAL::Assemble right hand side
 DEAL::Solve
 DEAL:cg::Starting value 27.1275
 DEAL:cg::Convergence step 1 value 1.97998e-14
-DEAL::Error    0.161536
+DEAL::Error    0.161172
 DEAL::Estimate 1.35839
 DEAL::Writing solution to <sol-00.gnuplot>...
 DEAL::
@@ -26,8 +26,8 @@ DEAL::Assemble multilevel matrix
 DEAL::Assemble right hand side
 DEAL::Solve
 DEAL:cg::Starting value 35.5356
-DEAL:cg::Convergence step 10 value 1.24790e-13
-DEAL::Error    0.164997
+DEAL:cg::Convergence step 14 value 3.21479e-13
+DEAL::Error    0.164760
 DEAL::Estimate 1.08528
 DEAL::Writing solution to <sol-01.gnuplot>...
 DEAL::
@@ -39,8 +39,8 @@ DEAL::Assemble multilevel matrix
 DEAL::Assemble right hand side
 DEAL::Solve
 DEAL:cg::Starting value 37.0552
-DEAL:cg::Convergence step 10 value 4.74024e-14
-DEAL::Error    0.113596
+DEAL:cg::Convergence step 14 value 6.05416e-13
+DEAL::Error    0.113503
 DEAL::Estimate 0.990460
 DEAL::Writing solution to <sol-02.gnuplot>...
 
@@ -54,8 +54,8 @@ DEAL::Assemble multilevel matrix
 DEAL::Assemble right hand side
 DEAL::Solve
 DEAL:cg::Starting value 38.5798
-DEAL:cg::Convergence step 12 value 7.56337e-14
-DEAL::Error    0.0101279
+DEAL:cg::Convergence step 17 value 2.64999e-13
+DEAL::Error    0.0101278
 DEAL::Estimate 0.0957571
 DEAL::Writing solution to <sol-10.gnuplot>...
 DEAL::
@@ -67,21 +67,22 @@ DEAL::Assemble multilevel matrix
 DEAL::Assemble right hand side
 DEAL::Solve
 DEAL:cg::Starting value 44.1721
-DEAL:cg::Convergence step 12 value 9.49720e-14
-DEAL::Error    0.00716966
+DEAL:cg::Convergence step 17 value 3.18657e-13
+DEAL::Error    0.00716962
 DEAL::Estimate 0.0681646
 DEAL::Writing solution to <sol-11.gnuplot>...
+DEAL::
 @endcode
 
 This log for instance shows that the number of conjugate gradient
-iteration steps is constant at approximately 12.
+iteration steps is constant at approximately 17.
 
 <h2>Postprocessing of the logfile</h2>
 
 @image html "step-39-convergence.png"
 Using the perl script <tt>postprocess.pl</tt>, we extract relevant
 data into <tt>output.dat</tt>, which can be used to plot graphs with
-<tt>gnuplot</tt>. The graph on the right for instance was produced with
+<tt>gnuplot</tt>. The graph above for instance was produced with
 
 @code
 set style data linespoints
index 8563ea41019be958616334b9c483f3e02cc8697a..2028ed867248368525f68608284c0ecd548f210c 100644 (file)
@@ -88,9 +88,16 @@ Functions::SlitSingularityFunction<2> exact_solution;
                                 // the information needed for the
                                 // local integration is provided by
                                 // MeshWorker::IntegrationInfo<dim>. Note
-                                // that this public interface cannot
+                                // that the signature of the functions cannot
                                 // be changed, because it is expected
                                 // by MeshWorker::integration_loop().
+
+                                // The first class defining local
+                                // integrators is responsible for
+                                // computing cell and face
+                                // matrices. It is used to assemble
+                                // the global matrix as well as the
+                                // level matrices.
 template <int dim>
 class MatrixIntegrator : public Subscriptor
 {
@@ -150,7 +157,8 @@ void MatrixIntegrator<dim>::bdry(
                             * fe.JxW(k);
 }
 
-
+                                // Interior faces use the interior
+                                // penalty method
 template <int dim>
 void MatrixIntegrator<dim>::face(
   MeshWorker::DoFInfo<dim>& dinfo1,
@@ -200,7 +208,12 @@ const double penalty = penalty1 + penalty2;
        }
 }
 
-
+                                // The second local integrator builds
+                                // the right hand side. In our
+                                // example, the right hand side
+                                // function is zero, such that only
+                                // the boundary condition is set here
+                                // in weak form.
 template <int dim>
 class RHSIntegrator : public Subscriptor
 {
@@ -247,6 +260,11 @@ void RHSIntegrator<dim>::face(MeshWorker::DoFInfo<dim>&,
 {}
 
 
+                                // The third local integrator is
+                                // responsible for the contributions
+                                // to the error estimate. This is the
+                                // standard energy estimator due to
+                                // Karakashian and Pascal (2003).
 template <int dim>
 class Estimator : public Subscriptor
 {
@@ -260,6 +278,10 @@ class Estimator : public Subscriptor
 };
 
 
+                                // The cell contribution is the
+                                // Laplacian of the discrete
+                                // solution, since the right hand
+                                // side is zero.
 template <int dim>
 void Estimator<dim>::cell(MeshWorker::DoFInfo<dim>& dinfo, typename MeshWorker::IntegrationInfo<dim>& info)
 {
@@ -274,7 +296,12 @@ void Estimator<dim>::cell(MeshWorker::DoFInfo<dim>& dinfo, typename MeshWorker::
   dinfo.value(0) = std::sqrt(dinfo.value(0));
 }
 
-
+                                // At the boundary, we use simply a
+                                // weighted form of the boundary
+                                // residual, namely the norm of the
+                                // difference between the finite
+                                // element solution and the correct
+                                // boundary condition.
 template <int dim>
 void Estimator<dim>::bdry(MeshWorker::DoFInfo<dim>& dinfo, typename MeshWorker::IntegrationInfo<dim>& info)
 {
@@ -295,6 +322,10 @@ void Estimator<dim>::bdry(MeshWorker::DoFInfo<dim>& dinfo, typename MeshWorker::
 }
 
 
+                                // Finally, on interior faces, the
+                                // estimator consists of the jumps of
+                                // the solution and its normal
+                                // derivative, weighted appropriately.
 template <int dim>
 void Estimator<dim>::face(MeshWorker::DoFInfo<dim>& dinfo1,
                          MeshWorker::DoFInfo<dim>& dinfo2,
@@ -326,11 +357,17 @@ void Estimator<dim>::face(MeshWorker::DoFInfo<dim>& dinfo1,
 
 
                                 // @sect3{The main class}
+
+                                // This class does the main job, like
+                                // in previous examples. For a
+                                // description of the functions
+                                // declared here, please refer to
+                                // the implementation below.
 template <int dim>
 class Step39
 {
   public:
-    typedef typename MeshWorker::IntegrationInfo<dim> CellInfo;
+    typedef MeshWorker::IntegrationInfo<dim> CellInfo;
     
     Step39(const FiniteElement<dim>& fe);
 
@@ -346,27 +383,66 @@ class Step39
     void solve ();
     void refine_grid ();
     void output_results (const unsigned int cycle) const;
-    
+
+                                    // The member objects related to
+                                    // the discretization are here.
     Triangulation<dim>        triangulation;
     const MappingQ1<dim>      mapping;
     const FiniteElement<dim>& fe;
     MGDoFHandler<dim>         mg_dof_handler;
     DoFHandler<dim>&          dof_handler;
 
+                                    // Then, we have the matrices and
+                                    // vectors related to the global
+                                    // discrete system.
     SparsityPattern      sparsity;
     SparseMatrix<double> matrix;
     Vector<double>       solution;
     Vector<double>       right_hand_side;
     BlockVector<double>  estimates;
-    
+
+                                    // Finally, we have a group of
+                                    // sparsity patterns and sparse
+                                    // matrices related to the
+                                    // multilevel preconditioner.
+                                    // First, we have a level matrix
+                                    // and its sparsity pattern.
     MGLevelObject<SparsityPattern> mg_sparsity;
-    MGLevelObject<SparsityPattern> mg_sparsity_dg_interface;
     MGLevelObject<SparseMatrix<double> > mg_matrix;
-    MGLevelObject<SparseMatrix<double> > mg_matrix_dg_up;
+    
+                                    // When we perform multigrid with
+                                    // local smoothing on locally
+                                    // refined meshes, additional
+                                    // matrices are required; see
+                                    // Kanschat (2004). Here is the
+                                    // sparsity pattern for these
+                                    // edge matrices. We only need
+                                    // one, because the pattern of
+                                    // the up matrix is the
+                                    // transpose of that of the down
+                                    // matrix. Actually, we do not
+                                    // care too much about these
+                                    // details, since the MeshWorker
+                                    // is filling these matrices.
+    MGLevelObject<SparsityPattern> mg_sparsity_dg_interface;
+                                    // The flux matrix at the
+                                    // refinement edge, coupling fine
+                                    // level degrees of freedom to
+                                    // coarse level.
     MGLevelObject<SparseMatrix<double> > mg_matrix_dg_down;
+                                    // The transpose of the flux
+                                    // matrix at the refinement edge,
+                                    // coupling coarse level degrees
+                                    // of freedom to fine level.
+    MGLevelObject<SparseMatrix<double> > mg_matrix_dg_up;
 };
 
 
+                                // The constructor simply sets up the
+                                // coarse grid and the
+                                // DoFHandler. The FiniteElement is
+                                // provided as a parameter to allow
+                                // flexibility.
 template <int dim>
 Step39<dim>::Step39(const FiniteElement<dim>& fe)
                :
@@ -379,49 +455,93 @@ Step39<dim>::Step39(const FiniteElement<dim>& fe)
 }
 
 
+                                // In this function, we set up the
+                                // dimension of the linear system and
+                                // the sparsity patterns for the
+                                // global matrix as well as the level
+                                // matrices.
 template <int dim>
 void
 Step39<dim>::setup_system()
 {
+                                  // First, we use the finite element
+                                  // to distribute degrees of
+                                  // freedom over the mesh and number
+                                  // them.
   dof_handler.distribute_dofs(fe);
   unsigned int n_dofs = dof_handler.n_dofs();
-  
+                                  // Then, we already know the size
+                                  // of the vectors representing
+                                  // finite element functions.
+  solution.reinit(n_dofs);
+  right_hand_side.reinit(n_dofs);
+
+                                  // Next, we set up the sparsity
+                                  // pattern for the global
+                                  // matrix. Since we do not know the
+                                  // row sizes in advance, we first
+                                  // fill a temporary
+                                  // CompressedSparsityPattern object
+                                  // and copy it to the regular
+                                  // SparsityPattern once it is
+                                  // complete.
   CompressedSparsityPattern c_sparsity(n_dofs);
   DoFTools::make_flux_sparsity_pattern(dof_handler, c_sparsity);
   sparsity.copy_from(c_sparsity);
   matrix.reinit(sparsity);
-  
-  solution.reinit(n_dofs);
-  right_hand_side.reinit(n_dofs);
 
   const unsigned int n_levels = triangulation.n_levels();
-  
+                                  // The global system is set up, now
+                                  // we attend to the level
+                                  // matrices. We resize all matrix
+                                  // objects to hold one matrix per level.
   mg_matrix.resize(0, n_levels-1);
   mg_matrix.clear();
   mg_matrix_dg_up.resize(0, n_levels-1);
   mg_matrix_dg_up.clear();
   mg_matrix_dg_down.resize(0, n_levels-1);
   mg_matrix_dg_down.clear();
-
+                                  // It is important to update the
+                                  // sparsity patterns after
+                                  // <tt>clear()</tt> was called for
+                                  // the level matrices, since the
+                                  // matrices lock the sparsity
+                                  // pattern through the Smartpointer
+                                  // ans Subscriptor mechanism.
   mg_sparsity.resize(0, n_levels-1);
   mg_sparsity_dg_interface.resize(0, n_levels-1);
-  
+
+                                  // Now all objects are prepared to
+                                  // hold one sparsity pattern or
+                                  // matrix per level. What's left is
+                                  // setting up the sparsity patterns
+                                  // on each level.
   for (unsigned int level=mg_sparsity.get_minlevel();
        level<=mg_sparsity.get_maxlevel();++level)
     {
-      CompressedSparsityPattern c_sparsity(mg_dof_handler.n_dofs(level));
-      CompressedSparsityPattern ci_sparsity;
-      if (level>0)
-       ci_sparsity.reinit(mg_dof_handler.n_dofs(level-1), mg_dof_handler.n_dofs(level));
-      
+                                      // These are roughly the same
+                                      // lines as above for the
+                                      // global matrix, now for each
+                                      // level.
+      CompressedSparsityPattern c_sparsity(mg_dof_handler.n_dofs(level));      
       MGTools::make_flux_sparsity_pattern(mg_dof_handler, c_sparsity, level);
-      if (level>0)
-       MGTools::make_flux_sparsity_pattern_edge(mg_dof_handler, ci_sparsity, level);
-      
       mg_sparsity[level].copy_from(c_sparsity);
       mg_matrix[level].reinit(mg_sparsity[level]);
+
+                                      // Additionally, we need to
+                                      // initialize the transfer
+                                      // matrices at the refinement
+                                      // edge between levels. They
+                                      // are stored at the index
+                                      // referring to the finer of
+                                      // the two indices, thus there
+                                      // is no such object on level
+                                      // 0.
       if (level>0)
        {
+         CompressedSparsityPattern ci_sparsity;
+         ci_sparsity.reinit(mg_dof_handler.n_dofs(level-1), mg_dof_handler.n_dofs(level));
+         MGTools::make_flux_sparsity_pattern_edge(mg_dof_handler, ci_sparsity, level);
          mg_sparsity_dg_interface[level].copy_from(ci_sparsity);
          mg_matrix_dg_up[level].reinit(mg_sparsity_dg_interface[level]);
          mg_matrix_dg_down[level].reinit(mg_sparsity_dg_interface[level]);
@@ -430,23 +550,61 @@ Step39<dim>::setup_system()
 }
 
 
+                                // In this function, we assemble the
+                                // global system matrix, where by
+                                // global we indicate that this is
+                                // the matrix of the discrete system
+                                // we solve and it is covering the
+                                // whole mesh.
 template <int dim>
 void
 Step39<dim>::assemble_matrix()
 {
+                                  // First, we need t set up the
+                                  // object providing the values we
+                                  // integrate. This object contains
+                                  // all FEValues and FEFaceValues
+                                  // objects needed and also
+                                  // maintains them automatically
+                                  // such that they always point to
+                                  // the current cell. To this end,
+                                  // we need to tell it first, where
+                                  // and what to compute,
   MeshWorker::IntegrationInfoBox<dim> info_box;
+                                  // namely, which quadrature
+                                  // formulas to use and
   const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1;
   info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points);
+                                  // which values to update in these
+                                  // points. We call
+                                  // <tt>initialize_update_flags</tt>
+                                  // first in order to set default
+                                  // values. Then, we add what we
+                                  // need additionally.
   info_box.initialize_update_flags();
   UpdateFlags update_flags = update_values | update_gradients;
   info_box.add_update_flags(update_flags, true, true, true, true);
   info_box.initialize(fe, mapping);
-  
+
+                                  // This is the object into which we
+                                  // integrate local data.
   MeshWorker::DoFInfo<dim> dof_info(dof_handler);
-  
+
+                                  // Finally, we need an object that
+                                  // assembles the local matrix into
+                                  // the global matrix.
   MeshWorker::Assembler::MatrixSimple<SparseMatrix<double> > assembler;
   assembler.initialize(matrix);
-  
+
+                                  // Now, we throw everything into a
+                                  // MeshWorker::loop(), which here
+                                  // traverses all active cells of
+                                  // the mesh, computes cell and face
+                                  // matrices and assembles them into
+                                  // the global matrix. We use the
+                                  // variable <tt>dof_handler</tt>
+                                  // here in order to use the global
+                                  // numbering of degrees of freedom.
   MeshWorker::integration_loop<dim, dim>(
     dof_handler.begin_active(), dof_handler.end(),
     dof_info, info_box,
@@ -457,6 +615,11 @@ Step39<dim>::assemble_matrix()
 }
 
 
+                                // Now, we do the same for the level
+                                // matrices. Not too surprisingly,
+                                // this function looks like a twin of
+                                // the previous one. Indeed, there
+                                // are only two minor differences.
 template <int dim>
 void
 Step39<dim>::assemble_mg_matrix()
@@ -471,10 +634,23 @@ Step39<dim>::assemble_mg_matrix()
 
   MeshWorker::DoFInfo<dim> dof_info(mg_dof_handler);
 
+                                  // Obviously, the assembler needs
+                                  // to be replaced by one filling
+                                  // level matrices. Note that it
+                                  // automatically fills the edge
+                                  // matrices as well.
   MeshWorker::Assembler::MGMatrixSimple<SparseMatrix<double> > assembler;
   assembler.initialize(mg_matrix);
   assembler.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down);
 
+                                  // Here is the other difference to
+                                  // the previous function: we run
+                                  // over all cells, not only the
+                                  // active ones. And we use
+                                  // <tt>mg_dof_handler</tt>, since
+                                  // we need the degrees of freedom
+                                  // on each level, not the global
+                                  // numbering.
   MeshWorker::integration_loop<dim, dim> (
     mg_dof_handler.begin(), mg_dof_handler.end(),
     dof_info, info_box,
@@ -485,6 +661,11 @@ Step39<dim>::assemble_mg_matrix()
 }
 
 
+                                // Here we have another clone of the
+                                // assemble function. The difference
+                                // to assembling the system matrix
+                                // consists in that we assemble a
+                                // vector here.
 template <int dim>
 void
 Step39<dim>::assemble_right_hand_side()
@@ -498,7 +679,17 @@ Step39<dim>::assemble_right_hand_side()
   info_box.initialize(fe, mapping);
   
   MeshWorker::DoFInfo<dim> dof_info(dof_handler);
-  
+
+                                  // Since this assembler alows us to
+                                  // fill several vectors, the
+                                  // interface is a little more
+                                  // complicated as above. The
+                                  // pointers to the vectors have to
+                                  // be stored in a NamedData
+                                  // object. While this seems to
+                                  // cause two extra lines of code
+                                  // here, it actually comes handy in
+                                  // more complex applications.
   MeshWorker::Assembler::ResidualSimple<Vector<double> > assembler;  
   NamedData<Vector<double>* > data;
   Vector<double>* rhs = &right_hand_side;
@@ -517,27 +708,52 @@ Step39<dim>::assemble_right_hand_side()
 }
 
 
+                                // Now that we have coded all
+                                // functions building the discrete
+                                // linear system, it is about time
+                                // that we actually solve it.
 template <int dim>
 void
 Step39<dim>::solve()
 {
+                                  // The solver of choice is
+                                  // conjugate gradient.
   SolverControl control(1000, 1.e-12);
   SolverCG<Vector<double> > cg(control);
 
-  GrowingVectorMemory<Vector<double> > mem;
+                                  // Now we are setting up the
+                                  // components of the multilevel
+                                  // preconditioner. First, we need
+                                  // transfer between grid
+                                  // levels. The object we are using
+                                  // here generates sparse matrices
+                                  // for these transfers.
   MGTransferPrebuilt<Vector<double> > mg_transfer;
   mg_transfer.build_matrices(mg_dof_handler);
+
+                                  // Then, we need an exact solver
+                                  // for the matrix on the coarsest
+                                  // level.
   FullMatrix<double> coarse_matrix;
   coarse_matrix.copy_from (mg_matrix[0]);
   MGCoarseGridHouseholder<double, Vector<double> > mg_coarse;
   mg_coarse.initialize(coarse_matrix);
-  typedef PreconditionSSOR<SparseMatrix<double> > RELAXATION;
+
+                                  // While transfer and coarse grid
+                                  // solver are pretty much generic,
+                                  // more flexibility is offered for
+                                  // the smoother. First, we choose
+                                  // Gauss-Seidel as our smoothing
+                                  // method.
+  GrowingVectorMemory<Vector<double> > mem;
+  typedef PreconditionSOR<SparseMatrix<double> > RELAXATION;
   MGSmootherRelaxation<SparseMatrix<double>, RELAXATION, Vector<double> >
     mg_smoother(mem);
   RELAXATION::AdditionalData smoother_data(1.);
   mg_smoother.initialize(mg_matrix, smoother_data);
   
-                                  // Do two smoothing steps per level
+                                  // Do two smoothing steps on each
+                                  // level.
   mg_smoother.set_steps(2);
                                   // Since the SOR method is not
                                   // symmetric, but we use conjugate
@@ -547,15 +763,18 @@ Step39<dim>::solve()
                                   // symmetric operator even for
                                   // nonsymmetric smoothers.
   mg_smoother.set_symmetric(true);
+                                  // The smoother class optionally
+                                  // implements the variable V-cycle,
+                                  // which we do not want here.
   mg_smoother.set_variable(false);
 
-                                  // We must wrap our matrices in an
-                                  // object having the required
-                                  // multiplication functions.
+                                  // Finally, we must wrap our
+                                  // matrices in an object having the
+                                  // required multiplication
+                                  // functions.
   MGMatrix<SparseMatrix<double>, Vector<double> > mgmatrix(&mg_matrix);
   MGMatrix<SparseMatrix<double>, Vector<double> > mgdown(&mg_matrix_dg_down);
   MGMatrix<SparseMatrix<double>, Vector<double> > mgup(&mg_matrix_dg_up);
-
   
                                   // Now, we are ready to set up the
                                   // V-cycle operator and the
@@ -563,18 +782,26 @@ Step39<dim>::solve()
   Multigrid<Vector<double> > mg(mg_dof_handler, mgmatrix,
                                mg_coarse, mg_transfer,
                                mg_smoother, mg_smoother);
+                                  // Let us not forget the edge
+                                  // matrices needed because of the
+                                  // adaptive refinement.
   mg.set_edge_flux_matrices(mgdown, mgup);
-  mg.set_debug(0);
-  mg_smoother.set_debug(0);
-  
+
+                                  // After all preparations, wrap the
+                                  // Multigrid object into another
+                                  // object, which can be used as a
+                                  // regular preconditioner,
   PreconditionMG<dim, Vector<double>,
     MGTransferPrebuilt<Vector<double> > >
     preconditioner(mg_dof_handler, mg, mg_transfer);
-  
+                                  // and use it to solve the system.
   cg.solve(matrix, solution, right_hand_side, preconditioner);
 }
 
-
+                                // Here we compare our finite element
+                                // solution with the (known) exact
+                                // solution and compute the mean
+                                // quadratic error of the gradient.
 template <int dim>
 void
 Step39<dim>::error()
@@ -584,42 +811,71 @@ Step39<dim>::error()
   
   QGauss<dim> quadrature(n_gauss_points);
   VectorTools::integrate_difference(mapping, dof_handler, solution, exact_solution,
-                                   cell_errors, quadrature, VectorTools::H1_norm);
+                                   cell_errors, quadrature, VectorTools::H1_seminorm);
   deallog << "Error    " << cell_errors.l2_norm() << std::endl;
 }
 
 
+                                // Another clone of the assemble
+                                // function. The big difference to
+                                // the previous ones is here that we
+                                // also have an input vector.
 template <int dim>
 double
 Step39<dim>::estimate()
 {
+                                  // The results of the estimator are
+                                  // stored in a vector with one
+                                  // entry per cell. Since cells in
+                                  // deal.II are not numbered, we
+                                  // have to create our own numbering
+                                  // in order to use this vector.
   estimates.block(0).reinit(triangulation.n_active_cells());
   unsigned int i=0;
   for (typename Triangulation<dim>::active_cell_iterator cell = triangulation.begin_active();
        cell != triangulation.end();++cell,++i)
     cell->set_user_index(i);
-  
+
+                                  // This starts like before,
   MeshWorker::IntegrationInfoBox<dim> info_box;
   const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1;
   info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points);
-  
+
+                                  // but now we need to notify the
+                                  // info box of the finite element
+                                  // functio we want to evaluate in
+                                  // the quadrature points. First, we
+                                  // create a NamedData object with
+                                  // this vector, which is the
+                                  // solution we just computed.
   NamedData<Vector<double>* > solution_data;
   solution_data.add(&solution, "solution");
-  
-  MeshWorker::VectorSelector cs;
-  MeshWorker::VectorSelector fs;
-  cs.add("solution", true, true, true);
-  fs.add("solution", true, true, false);
-  
-  info_box.cell_selector = cs;
-  info_box.boundary_selector = fs;
-  info_box.face_selector = fs;
+                                  // Then, we tell the Meshworker::VectorSelector
+                                  // for cells, that we need the
+                                  // second derivatives of this
+                                  // solution (to compute the Laplacian).
+  info_box.cell_selector.add("solution", false, false, true);
+                                  // On interior and boundary faces,
+                                  // we need the function values and
+                                  // the first derivatives.
+  info_box.boundary_selector.add("solution", true, true, false);
+  info_box.face_selector.add("solution", true, true, false);
+
+                                  // And we continue as before, with
+                                  // the exception that the default
+                                  // update flags are already
+                                  // adjusted to the values and
+                                  // derivatives we requested above.
   info_box.initialize_update_flags();
   info_box.add_update_flags(update_quadrature_points, false, true, false, false);
   info_box.initialize(fe, mapping, solution_data);
   
   MeshWorker::DoFInfo<dim> dof_info(dof_handler);
-  
+
+                                  // The assembler stores one number
+                                  // per cell, but else this is the
+                                  // same as in the computation of
+                                  // the right hand side.
   MeshWorker::Assembler::CellsAndFaces<double> assembler;  
   NamedData<BlockVector<double>* > out_data;
   BlockVector<double>* est = &estimates;
@@ -637,6 +893,7 @@ Step39<dim>::estimate()
 }
 
 
+                                // Some graphical output
 template <int dim>
 void Step39<dim>::output_results (const unsigned int cycle) const
 {
@@ -661,7 +918,9 @@ void Step39<dim>::output_results (const unsigned int cycle) const
   data_out.write_gnuplot(gnuplot_output);
 }
 
-
+                                // And finally the adaptive loop,
+                                // more or less like in previous
+                                // examples.
 template <int dim>
 void
 Step39<dim>::run(unsigned int n_steps)

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.