]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Re-indent the entire file. Fix up a few places with formulas and other formatting.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 23 Feb 2008 19:26:27 +0000 (19:26 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 23 Feb 2008 19:26:27 +0000 (19:26 +0000)
git-svn-id: https://svn.dealii.org/trunk@15760 0785d39b-7218-0410-832d-ea1e28bc413d

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

index a726cb4ed3c11b5cefb9f273e5f563fdd7b9ad21..516dfbb38b8bfd47f664da04ff1be44bc71af6db 100644 (file)
 /*    further information on this license.                        */
 
 
-                        // @sect3{Include files}
+                                // @sect3{Include files}
                         
-                    // As usual, we start by including 
-                    // some well-known files. 
+                                // As usual, we start by including 
+                                // some well-known files. 
 #include <base/quadrature_lib.h>
 #include <base/logstream.h>
 #include <base/function.h>
 #include <numerics/data_out.h>
 #include <numerics/error_estimator.h>
 
-                    // As in 
-                    // step-29, we include the libary
-                    // for the sparse direct solver
-                    // UMFPACK.
+                                // As in step-29, we include the libary for
+                                // the sparse direct solver UMFPACK.
 #include <lac/sparse_direct.h>
 
-                    // This includes the libary for the 
-                    // incomplete LU factorization that will 
-                    // be used as a preconditioner in 3D.
+                                // This includes the libary for the 
+                                // incomplete LU factorization that will 
+                                // be used as a preconditioner in 3D.
 #include <lac/sparse_ilu.h>
 
-                    // This is C++:
+                                // This is C++:
 #include <fstream>
 #include <sstream>
 
-                    // As in all programs, the namespace is set
-                    // to dealii.
+                                // As in all programs, the namespace is set
+                                // to dealii.
 using namespace dealii;
 
-                        // @sect3{Defining the inner preconditioner type}
+                                // @sect3{Defining the inner preconditioner type}
                         
-                    // As explained in the introduction, we
-                    // are going to use different preconditioners
-                    // for two and three space dimensions, 
-                    // respectively. We differentiate between
-                    // them by the use of the spatial dimension
-                    // as a template parameter. See step-4 for
-                    // details on templates.
-                    // We are not going to create any preconditioner
-                    // object here, all we do is to create a 
-                    // data structure that holds the information
-                    // on it so we can write our program in a
-                    // dimension-independent way.
+                                // As explained in the introduction, we are
+                                // going to use different preconditioners for
+                                // two and three space dimensions,
+                                // respectively. We differentiate between
+                                // them by the use of the spatial dimension
+                                // as a template parameter. See step-4 for
+                                // details on templates.  We are not going to
+                                // create any preconditioner object here, all
+                                // we do is to create a data structure that
+                                // holds the information on it so we can
+                                // write our program in a
+                                // dimension-independent way.
 template <int dim>
 struct InnerPreconditioner;
 
-                    // In 2D, we are going to use a sparse direct
-                    // solve as preconditioner. The syntax is 
-                    // known from step-29.
+                                // In 2D, we are going to use a sparse direct
+                                // solve as preconditioner. The syntax is 
+                                // known from step-29.
 template <>
 struct InnerPreconditioner<2> 
 {
     typedef SparseDirectUMFPACK type;
 };
 
-                    // And the ILU preconditioning in 3D, called
-                    // by <code>SparseILU@<double></code>.
+                                // And the ILU preconditioning in 3D, called
+                                // by <code>SparseILU@<double></code>.
 template <>
 struct InnerPreconditioner<3> 
 {
@@ -105,14 +103,13 @@ struct InnerPreconditioner<3>
 };
 
 
-                        // @sect3{The <code>StokesProblem</code> class template}
+                                // @sect3{The <code>StokesProblem</code> class template}
                     
-                    // This is an adaptation of step-20, 
-                    // so the main class and the data types
-                    // are the same as used there. In this
-                    // example we also use adaptive grid
-                    // refinement, which is handled in complete
-                    // analogy to step-6.
+                                // This is an adaptation of step-20, so the
+                                // main class and the data types are the same
+                                // as used there. In this example we also use
+                                // adaptive grid refinement, which is handled
+                                // in complete analogy to step-6.
 template <int dim>
 class StokesProblem 
 {
@@ -141,38 +138,39 @@ class StokesProblem
     BlockVector<double> solution;
     BlockVector<double> system_rhs;
 
-                    // This one is new: We shall use a so-called
-                    // shared pointer structure to access
-                    // the preconditioner. This provides
-                    // flexibility when using the object 
-                    // that the pointer refers to, as e.g. 
-                    // the reset option.
+                                    // This one is new: We shall use a
+                                    // so-called shared pointer structure to
+                                    // access the preconditioner. This
+                                    // provides flexibility when using the
+                                    // object that the pointer refers to, as
+                                    // e.g. the reset option.
     boost::shared_ptr<typename InnerPreconditioner<dim>::type> A_preconditioner;
 };
 
-                                       // @sect3{Boundary values and right hand side}
-
-                    // As in step-20 and most other example
-                    // programs, the next task is to define
-                    // the parameter functions for the PDE:
-                    // For the Stokes problem, we are going to
-                    // use pressure boundary values at some portion
-                    // of the boundary (Neumann-type), and 
-                    // boundary conditions on the velocity 
-                    // (Dirichlet type) on the rest of the boundary.
-                    // The pressure boundary condition is
-                    // scalar, and so is the respective function, 
-                    // whereas the Dirichlet (velocity) 
-                    // condition is vector-valued. Due to the
-                    // structure of deal.II's libraries, we have to
-                    // define the function on the (u,p)-space, but
-                    // we are going to filter out the pressure
-                    // component when condensating the Dirichlet
-                    // data in <code>assemble_system</code>.
+                                // @sect3{Boundary values and right hand side}
+
+                                // As in step-20 and most other example
+                                // programs, the next task is to define the
+                                // parameter functions for the PDE: For the
+                                // Stokes problem, we are going to use
+                                // pressure boundary values at some portion
+                                // of the boundary (Neumann-type), and
+                                // boundary conditions on the velocity
+                                // (Dirichlet type) on the rest of the
+                                // boundary.  The pressure boundary condition
+                                // is scalar, and so is the respective
+                                // function, whereas the Dirichlet (velocity)
+                                // condition is vector-valued. Due to the
+                                // structure of deal.II's libraries, we have
+                                // to define the function on the (u,p)-space,
+                                // but we are going to filter out the
+                                // pressure component when condensating the
+                                // Dirichlet data in
+                                // <code>assemble_system</code>.
                     
-                    // Given the problem described in the 
-                    // introduction, we know which values to 
-                    // set for the respective functions.
+                                // Given the problem described in the 
+                                // introduction, we know which values to 
+                                // set for the respective functions.
 template <int dim>
 class PressureBoundaryValues : public Function<dim> 
 {
@@ -211,7 +209,7 @@ class BoundaryValues : public Function<dim>
 template <int dim>
 double
 BoundaryValues<dim>::value (const Point<dim>  &p,
-                           const unsigned int component) const 
+                           const unsigned int component) const 
 {
   if (component == 0)
     return (p[0] < 0 ? -1 : (p[0] > 0 ? 1 : 0));
@@ -222,7 +220,7 @@ BoundaryValues<dim>::value (const Point<dim>  &p,
 template <int dim>
 void
 BoundaryValues<dim>::vector_value (const Point<dim> &p,
-                                  Vector<double>   &values) const 
+                                  Vector<double>   &values) const 
 {
   for (unsigned int c=0; c<this->n_components; ++c)
     values(c) = BoundaryValues<dim>::value (p, c);
@@ -230,8 +228,8 @@ BoundaryValues<dim>::vector_value (const Point<dim> &p,
 
 
 
-                    // We implement similar functions
-                    // for the right hand side.
+                                // We implement similar functions
+                                // for the right hand side.
 template <int dim>
 class RightHandSide : public Function<dim> 
 {
@@ -266,27 +264,29 @@ RightHandSide<dim>::vector_value (const Point<dim> &p,
 }
 
 
-                        // @sect3{Linear solvers and preconditioners}
+                                // @sect3{Linear solvers and preconditioners}
                         
-                    // The linear solvers and preconditioners are 
-                    // discussed extensively in the introduction. Here,
-                    // we create the respective objects that will be used.
+                                // The linear solvers and preconditioners are
+                                // discussed extensively in the
+                                // introduction. Here, we create the
+                                // respective objects that will be used.
                         
-                        // @sect4{The <code>InverseMatrix</code> class template}
+                                // @sect4{The <code>InverseMatrix</code> class template}
                         
-                    // This is going to represent the data
-                    // structure for an inverse matrix. This class
-                    // is derived from the one in step-20. The
-                    // only difference is that we now
-                    // do include a preconditioner to the matrix.
-                    // This is going to happen via a template parameter
-                    // <code>class Preconditioner</code>, so 
-                    // the preconditioner type will be set when
-                    // an <code>InverseMatrix</code> object is 
-                    // created. The member function 
-                    // <code>vmult</code> is, as in
-                    // step-20, a multiplication with a vector,
-                    // obtained by solving a linear system.
+                                // This is going to represent the data
+                                // structure for an inverse matrix. This
+                                // class is derived from the one in
+                                // step-20. The only difference is that we
+                                // now do include a preconditioner to the
+                                // matrix.  This is going to happen via a
+                                // template parameter <code>class
+                                // Preconditioner</code>, so the
+                                // preconditioner type will be set when an
+                                // <code>InverseMatrix</code> object is
+                                // created. The member function
+                                // <code>vmult</code> is, as in step-20, a
+                                // multiplication with a vector, obtained by
+                                // solving a linear system.
 template <class Matrix, class Preconditioner>
 class InverseMatrix : public Subscriptor
 {
@@ -307,62 +307,57 @@ class InverseMatrix : public Subscriptor
 
 template <class Matrix, class Preconditioner>
 InverseMatrix<Matrix,Preconditioner>::InverseMatrix (const Matrix &m,
-                                const Preconditioner &preconditioner)
-                 :
-                 matrix (&m),
-                 preconditioner (preconditioner)
+                                                    const Preconditioner &preconditioner)
+               :
+               matrix (&m),
+               preconditioner (preconditioner)
 {}
 
 
-                    // This is the implementation of the 
-                    // <code>vmult</code> function. We note 
-                    // two things: 
+                                // This is the implementation of the 
+                                // <code>vmult</code> function. We note 
+                                // two things: 
                     
-                    // Firstly, we use
-                    // a rather large tolerance for the
-                    // solver control. The reason for this
-                    // is that the function is used very 
-                    // frequently, and hence, any additional
-                    // effort to make the residual in
-                    // the CG solve smaller makes the
-                    // solution more expensive. Note that
-                    // we do not only use this class as a
-                    // preconditioner for the Schur complement,
-                    // but also when forming the inverse of
-                    // the Laplace matrix - which has to
-                    // be accurate in order to obtain a
-                    // solution to the right problem.
+                                // Firstly, we use a rather large tolerance
+                                // for the solver control. The reason for
+                                // this is that the function is used very
+                                // frequently, and hence, any additional
+                                // effort to make the residual in the CG
+                                // solve smaller makes the solution more
+                                // expensive. Note that we do not only use
+                                // this class as a preconditioner for the
+                                // Schur complement, but also when forming
+                                // the inverse of the Laplace matrix - which
+                                // has to be accurate in order to obtain a
+                                // solution to the right problem.
                     
-                    // Secondly, we catch exceptions from
-                    // the solver at this stage. While this
-                    // is not of greater interest our
-                    // general setting with the requirement
-                    // of accurate inverses (and we indeed 
-                    // abort the program when any exception
-                    // occurs), the situation would
-                    // change if an object of the class
-                    // <code>InverseMatrix</code> is only 
-                    // used for preconditioning. In such a 
-                    // setting, one could imagine to use 
-                    // a few CG sweeps as a preconditioner - 
-                    // which is done e.g. for mass
-                    // matrices, see the results section
-                    // below. Using
-                    // <code>catch (SolverControl::NoConvergence) {}</code>
-                    // in conjunction with only a few iterations,
-                    // say 10, would result in that effect - 
-                    // the program would continue to run
-                    // even though the solver has not converged.
-                    // Note, though, that applying the CG method 
-                    // is not a linear operation (see the
-                    // actual CG algorithm for details
-                    // on that), so unconverged
-                    // preconditioners are to be used with
-                    // care in order to not yield a wrong
-                    // solution.
+                                // Secondly, we catch exceptions from the
+                                // solver at this stage. While this is not of
+                                // greater interest our general setting with
+                                // the requirement of accurate inverses (and
+                                // we indeed abort the program when any
+                                // exception occurs), the situation would
+                                // change if an object of the class
+                                // <code>InverseMatrix</code> is only used
+                                // for preconditioning. In such a setting,
+                                // one could imagine to use a few CG sweeps
+                                // as a preconditioner - which is done
+                                // e.g. for mass matrices, see the results
+                                // section below. Using <code>catch
+                                // (SolverControl::NoConvergence) {}</code>
+                                // in conjunction with only a few iterations,
+                                // say 10, would result in that effect - the
+                                // program would continue to run even though
+                                // the solver has not converged.  Note,
+                                // though, that applying the CG method is not
+                                // a linear operation (see the actual CG
+                                // algorithm for details on that), so
+                                // unconverged preconditioners are to be used
+                                // with care in order to not yield a wrong
+                                // solution.
 template <class Matrix, class Preconditioner>
 void InverseMatrix<Matrix,Preconditioner>::vmult (Vector<double>       &dst,
-                          const Vector<double> &src) const
+                                                 const Vector<double> &src) const
 {
   SolverControl solver_control (src.size(), 1e-6*src.l2_norm());
   SolverCG<> cg (solver_control, vector_memory);
@@ -380,31 +375,30 @@ void InverseMatrix<Matrix,Preconditioner>::vmult (Vector<double>       &dst,
 }
 
 
-                        // @sect4{The <code>SchurComplement</code> class template}
-
-                    // This class implements the Schur 
-                    // complement discussed in the introduction.
-                    // It is in analogy to step-20.
-                    // Though, we now call it with a template
-                    // parameter <code>Preconditioner</code>
-                    // in order to access that when specifying
-                    // the respective type of the inverse
-                    // matrix class. As a consequence of the
-                    // definition above, the declaration
-                    // <code>InverseMatrix</code> 
-                    // now contains the second template
-                    // parameter from preconditioning as above,
-                    // which effects the <code>SmartPointer</code>
-                    // object <code>m_inverse</code> as well.
+                                // @sect4{The <code>SchurComplement</code> class template}
+
+                                // This class implements the Schur complement
+                                // discussed in the introduction.  It is in
+                                // analogy to step-20.  Though, we now call
+                                // it with a template parameter
+                                // <code>Preconditioner</code> in order to
+                                // access that when specifying the respective
+                                // type of the inverse matrix class. As a
+                                // consequence of the definition above, the
+                                // declaration <code>InverseMatrix</code> now
+                                // contains the second template parameter
+                                // from preconditioning as above, which
+                                // effects the <code>SmartPointer</code>
+                                // object <code>m_inverse</code> as well.
 template <class Preconditioner>
 class SchurComplement : public Subscriptor
 {
   public:
     SchurComplement (const BlockSparseMatrix<double> &A,
-             const InverseMatrix<SparseMatrix<double>, Preconditioner> &Minv);
+                    const InverseMatrix<SparseMatrix<double>, Preconditioner> &Minv);
 
     void vmult (Vector<double>       &dst,
-        const Vector<double> &src) const;
+               const Vector<double> &src) const;
 
   private:
     const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
@@ -418,18 +412,18 @@ class SchurComplement : public Subscriptor
 template <class Preconditioner>
 SchurComplement<Preconditioner>::
 SchurComplement (const BlockSparseMatrix<double> &A,
-         const InverseMatrix<SparseMatrix<double>,Preconditioner> &Minv)
-        :
-        system_matrix (&A),
-        m_inverse (&Minv),
-        tmp1 (A.block(0,0).m()),
-        tmp2 (A.block(0,0).m())
+                const InverseMatrix<SparseMatrix<double>,Preconditioner> &Minv)
+               :
+               system_matrix (&A),
+               m_inverse (&Minv),
+               tmp1 (A.block(0,0).m()),
+               tmp2 (A.block(0,0).m())
 {}
 
 
 template <class Preconditioner>
 void SchurComplement<Preconditioner>::vmult (Vector<double>       &dst,
-                         const Vector<double> &src) const
+                                            const Vector<double> &src) const
 {
   system_matrix->block(0,1).vmult (tmp1, src);
   m_inverse->vmult (tmp2, tmp1);
@@ -437,30 +431,32 @@ void SchurComplement<Preconditioner>::vmult (Vector<double>       &dst,
 }
 
 
-                        // @sect3{StokesProblem class implementation}
+                                // @sect3{StokesProblem class implementation}
                         
-                        // @sect4{StokesProblem::StokesProblem}
-
-                    // The constructor of this class looks very
-                    // similar to the one of step-20. The constructor
-                    // initializes the variables for the polynomial
-                    // degree, triangulation, finite element system
-                    // and the dof handler. The underlying polynomial
-                    // functions are of order <code>degree+1</code> for
-                    // the vector-valued velocity components and 
-                    // of order <code>degree</code> in pressure.
-                    // This gives the LBB-stable element pair 
-                    // Q(degree+1)Q(degree).
-                    
-                    // Note that we initialize the triangulation
-                    // with a MeshSmoothing argument, which ensures
-                    // that the refinement of cells is done
-                    // in a way that the approximation of the
-                    // PDE solution remains well-behaved (problems
-                    // arise if grids are too unstructered),
-                    // see the documentation of
-                    // <code>Triangulation::MeshSmoothing</code>
-                    // for details.
+                                // @sect4{StokesProblem::StokesProblem}
+
+                                // The constructor of this class looks very
+                                // similar to the one of step-20. The
+                                // constructor initializes the variables for
+                                // the polynomial degree, triangulation,
+                                // finite element system and the dof
+                                // handler. The underlying polynomial
+                                // functions are of order
+                                // <code>degree+1</code> for the
+                                // vector-valued velocity components and of
+                                // order <code>degree</code> in pressure.
+                                // This gives the LBB-stable element pair
+                                // Q(degree+1)Q(degree).
+                                //                    
+                                // Note that we initialize the triangulation
+                                // with a MeshSmoothing argument, which
+                                // ensures that the refinement of cells is
+                                // done in a way that the approximation of
+                                // the PDE solution remains well-behaved
+                                // (problems arise if grids are too
+                                // unstructered), see the documentation of
+                                // <code>Triangulation::MeshSmoothing</code>
+                                // for details.
 template <int dim>
 StokesProblem<dim>::StokesProblem (const unsigned int degree)
                 :
@@ -472,72 +468,72 @@ StokesProblem<dim>::StokesProblem (const unsigned int degree)
 {}
 
 
-                        // @sect4{StokesProblem::setup_dofs}
+                                // @sect4{StokesProblem::setup_dofs}
                         
-                    // Given a mesh, this function associates
-                    // the degrees of freedom with it and
-                    // creates the corresponding matrices and
-                    // vectors.
+                                // Given a mesh, this function associates
+                                // the degrees of freedom with it and
+                                // creates the corresponding matrices and
+                                // vectors.
 template <int dim>
 void StokesProblem<dim>::setup_dofs ()
 {
-                   // Release preconditioner from
-                   // previous steps since it
-                   // will definitely not be needed
-                   // any more after this point.
+                                  // Release preconditioner from
+                                  // previous steps since it
+                                  // will definitely not be needed
+                                  // any more after this point.
   A_preconditioner.reset ();
   
   dof_handler.distribute_dofs (fe); 
   
-                    // In order to make the ILU preconditioner
-                    // (in 3D) to work efficiently, the dofs
-                    // are renumbered using the Cuthill-McKee
-                    // algorithm. Though, the block structure 
-                    // of velocity and pressure shall be as in
-                    // step-20. This is done in two steps. First,
-                    // all dofs are renumbered by 
-                    // <code>DoFRenumbering::Cuthill_McKee</code>,
-                    // and then we renumber once again by 
-                    // components. Since 
-                    // <code>DoFRenumbering::component_wise</code>
-                    // does not touch the renumbering within 
-                    // the individual blocks, the basic
-                    // renumbering from Cuthill-McKee remains.
+                                  // In order to make the ILU preconditioner
+                                  // (in 3D) to work efficiently, the dofs
+                                  // are renumbered using the Cuthill-McKee
+                                  // algorithm. Though, the block structure
+                                  // of velocity and pressure shall be as in
+                                  // step-20. This is done in two
+                                  // steps. First, all dofs are renumbered by
+                                  // <code>DoFRenumbering::Cuthill_McKee</code>,
+                                  // and then we renumber once again by
+                                  // components. Since
+                                  // <code>DoFRenumbering::component_wise</code>
+                                  // does not touch the renumbering within
+                                  // the individual blocks, the basic
+                                  // renumbering from Cuthill-McKee remains.
   DoFRenumbering::Cuthill_McKee (dof_handler);
 
-                    // There is one more change: There
-                    // is no reason in creating <code>dim</code>
-                    // blocks for the velocity components,
-                    // so they can all be grouped in only
-                    // one block. The vector 
-                    // <code>block_component</code> does precisely
-                    // this: velocity values correspond to block
-                    // 0, and pressure values will sit in block
-                    // 1.
+                                  // There is one more change: There is no
+                                  // reason in creating <code>dim</code>
+                                  // blocks for the velocity components, so
+                                  // they can all be grouped in only one
+                                  // block. The vector
+                                  // <code>block_component</code> does
+                                  // precisely this: velocity values
+                                  // correspond to block 0, and pressure
+                                  // values will sit in block 1.
   std::vector<unsigned int> block_component (dim+1,0);
   block_component[dim] = 1;
   DoFRenumbering::component_wise (dof_handler, block_component);
 
-                    // Since we use adaptively refined grids
-                    // the constraint matrix for hanging
-                    // node constraints is generated from
-                    // the dof handler.
+                                  // Since we use adaptively refined grids
+                                  // the constraint matrix for hanging node
+                                  // constraints is generated from the dof
+                                  // handler.
   hanging_node_constraints.clear ();
   DoFTools::make_hanging_node_constraints (dof_handler,
                                            hanging_node_constraints);
   hanging_node_constraints.close ();
 
-                    // In analogy to step-20, we count
-                    // the dofs in the individual components.
-                    // We could do this in the same way as
-                    // there, but we want to operate on
-                    // the block structure we used already for 
-                    // the renumbering: The function 
-                    // <code>DoFTools::count_dofs_per_block</code>
-                    // does the same as 
-                    // <code>DoFTools::count_dofs_per_component</code>,
-                    // but now grouped as velocity and 
-                    // pressure block via <code>block_component</code>.
+                                  // In analogy to step-20, we count the dofs
+                                  // in the individual components.  We could
+                                  // do this in the same way as there, but we
+                                  // want to operate on the block structure
+                                  // we used already for the renumbering: The
+                                  // function
+                                  // <code>DoFTools::count_dofs_per_block</code>
+                                  // does the same as
+                                  // <code>DoFTools::count_dofs_per_component</code>,
+                                  // but now grouped as velocity and pressure
+                                  // block via <code>block_component</code>.
   std::vector<unsigned int> dofs_per_block (2);
   DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);  
   const unsigned int n_u = dofs_per_block[0],
@@ -551,40 +547,37 @@ void StokesProblem<dim>::setup_dofs ()
             << " (" << n_u << '+' << n_p << ')'
             << std::endl;
 
-                    // Release the memory previously attached
-                    // to the system matrix and untie it
-                    // from the old sparsity pattern prior to
-                    // generating the current data structure.
+                                  // Release the memory previously attached
+                                  // to the system matrix and untie it from
+                                  // the old sparsity pattern prior to
+                                  // generating the current data structure.
   system_matrix.clear ();
       
-                    // The next task is to allocate a
-                    // sparsity pattern for the system
-                    // matrix we will create. We could do
-                    // this in the same way as in step-20, 
-                    // though, there is a major reason
-                    // not to do so. In 3D, the function
-                    // <code>DoFTools::max_couplings_between_dofs</code>
-                    // yields a very large number for the
-                    // coupling between the individual dofs,
-                    // so that the memory initially provided
-                    // for the creation of the sparsity pattern
-                    // of the matrix is far too much - so 
-                    // much actually that it won't even fit
-                    // into the physical memory of most
-                    // systems already for moderately-sized 3D
-                    // problems. See also the discussion in
-                    // step-18.
-                    // Instead, we use a temporary object of
-                    // the class
-                    // BlockCompressedSparsityPattern,
-                    // which is a block version of the
-                    // compressed sparsity patterns from
-                    // step-11 and step-18. All this is done
-                    // inside a new scope, which means that
-                    // the memory of <code>csp</code> will be
-                    // released once the information has been
-                    // copied to 
-                    // <code>sparsity_pattern</code>.
+                                  // The next task is to allocate a sparsity
+                                  // pattern for the system matrix we will
+                                  // create. We could do this in the same way
+                                  // as in step-20, though, there is a major
+                                  // reason not to do so. In 3D, the function
+                                  // <code>DoFTools::max_couplings_between_dofs</code>
+                                  // yields a very large number for the
+                                  // coupling between the individual dofs, so
+                                  // that the memory initially provided for
+                                  // the creation of the sparsity pattern of
+                                  // the matrix is far too much - so much
+                                  // actually that it won't even fit into the
+                                  // physical memory of most systems already
+                                  // for moderately-sized 3D problems. See
+                                  // also the discussion in step-18.
+                                  // Instead, we use a temporary object of
+                                  // the class
+                                  // BlockCompressedSparsityPattern, which is
+                                  // a block version of the compressed
+                                  // sparsity patterns from step-11 and
+                                  // step-18. All this is done inside a new
+                                  // scope, which means that the memory of
+                                  // <code>csp</code> will be released once
+                                  // the information has been copied to
+                                  // <code>sparsity_pattern</code>.
   {
     BlockCompressedSparsityPattern csp;
 
@@ -604,10 +597,10 @@ void StokesProblem<dim>::setup_dofs ()
   std::ofstream out ("sparsity_pattern.gpl");
   sparsity_pattern.block(0,0).print_gnuplot(out);
 
-                    // Finally, the system matrix,
-                    // solution and right hand side are 
-                    // created from the block
-                    // structure as in step-20.
+                                  // Finally, the system matrix,
+                                  // solution and right hand side are 
+                                  // created from the block
+                                  // structure as in step-20.
   system_matrix.reinit (sparsity_pattern);
                                    
   solution.reinit (2);
@@ -622,16 +615,16 @@ void StokesProblem<dim>::setup_dofs ()
 }
 
 
-                        // @sect4{StokesProblem::assemble_system}
+                                // @sect4{StokesProblem::assemble_system}
                         
-                    // The assembly process follows the
-                    // discussion in step-20 and in the
-                    // introduction. We use the well-known
-                    // abbreviations for the data structures
-                    // that hold the local matrix, right
-                    // hand side, and global
-                    // numbers of the degrees of freedom
-                    // for the present cell.
+                                // The assembly process follows the
+                                // discussion in step-20 and in the
+                                // introduction. We use the well-known
+                                // abbreviations for the data structures
+                                // that hold the local matrix, right
+                                // hand side, and global
+                                // numbers of the degrees of freedom
+                                // for the present cell.
 template <int dim>
 void StokesProblem<dim>::assemble_system () 
 {
@@ -662,13 +655,13 @@ void StokesProblem<dim>::assemble_system ()
 
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
   
-                    // As usual, we create objects that
-                    // hold the functions for the right 
-                    // hand side and Neumann boundary
-                    // function, and, additionally,
-                    // an array that holds the respective
-                    // function values at the quadrature
-                    // points.
+                                  // As usual, we create objects that
+                                  // hold the functions for the right 
+                                  // hand side and Neumann boundary
+                                  // function, and, additionally,
+                                  // an array that holds the respective
+                                  // function values at the quadrature
+                                  // points.
   const PressureBoundaryValues<dim> pressure_boundary_values;
   
   std::vector<double>               boundary_values (n_face_q_points);
@@ -684,11 +677,11 @@ void StokesProblem<dim>::assemble_system ()
   const FEValuesExtractors::Vector velocities (0);
   const FEValuesExtractors::Scalar pressure (dim);
 
-                    // This starts the loop over all
-                    // cells. With the abbreviations
-                    // <code>extract_u</code> etc. 
-                    // introduced above, it is 
-                    // evident what is going on.
+                                  // This starts the loop over all
+                                  // cells. With the abbreviations
+                                  // <code>extract_u</code> etc. 
+                                  // introduced above, it is 
+                                  // evident what is going on.
   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
@@ -711,25 +704,29 @@ void StokesProblem<dim>::assemble_system ()
              const double        div_phi_i_u   = fe_values[velocities].divergence (i, q);
              const double        phi_i_p       = fe_values[pressure].value (i, q);
           
-          for (unsigned int j=0; j<dofs_per_cell; ++j)
-          {
+             for (unsigned int j=0; j<dofs_per_cell; ++j)
+               {
                  const SymmetricTensor<2,dim>
                    phi_j_grads_u = fe_values[velocities].symmetric_gradient (j, q);
                  const double        div_phi_j_u   = fe_values[velocities].divergence (j, q);
                  const double        phi_j_p       = fe_values[pressure].value (j, q);
-                    // Note the way we write the 
-                    // contributions
-                    // <code> phi_i_p * phi_j_p </code>,
-                    // yielding a pressure mass matrix,
-                    // into the same data structure as
-                    // the terms for the actual 
-                    // Stokes system - in accordance with
-                    // the description in the introduction.
-                    // They won't be mixed up, since
-                    // <code>phi_i_p * phi_j_p</code>
-                    // is only non-zero when all the
-                    // other terms vanish (and the other
-                    // way around).
+                                                  // Note the way we write
+                                                  // the contributions <code>
+                                                  // phi_i_p * phi_j_p
+                                                  // </code>, yielding a
+                                                  // pressure mass matrix,
+                                                  // into the same data
+                                                  // structure as the terms
+                                                  // for the actual Stokes
+                                                  // system - in accordance
+                                                  // with the description in
+                                                  // the introduction.  They
+                                                  // won't be mixed up, since
+                                                  // <code>phi_i_p *
+                                                  // phi_j_p</code> is only
+                                                  // non-zero when all the
+                                                  // other terms vanish (and
+                                                  // the other way around).
                                                   //
                                                   // Note also that operator*
                                                   // is overloaded for
@@ -738,27 +735,27 @@ void StokesProblem<dim>::assemble_system ()
                                                   // product between the two
                                                   // tensors in the first
                                                   // line:
-            local_matrix(i,j) += (phi_i_grads_u * phi_j_grads_u
-                      - div_phi_i_u * phi_j_p
-                      - phi_i_p * div_phi_j_u
-                      + phi_i_p * phi_j_p)
-                         * fe_values.JxW(q);     
-
-          }
-          const unsigned int component_i =
-            fe.system_to_component_index(i).first;
-          local_rhs(i) += fe_values.shape_value(i,q) * 
-                          rhs_values[q](component_i) *
-                          fe_values.JxW(q);
-        }
-      }
+                 local_matrix(i,j) += (phi_i_grads_u * phi_j_grads_u
+                                       - div_phi_i_u * phi_j_p
+                                       - phi_i_p * div_phi_j_u
+                                       + phi_i_p * phi_j_p)
+                                      * fe_values.JxW(q);     
+
+               }
+             const unsigned int component_i =
+               fe.system_to_component_index(i).first;
+             local_rhs(i) += fe_values.shape_value(i,q) * 
+                             rhs_values[q](component_i) *
+                             fe_values.JxW(q);
+           }
+       }
       
-                    // Here we add the contributions from
-                    // Neumann (pressure) boundary conditions
-                    // at faces on the domain boundary that
-                    // have the boundary flag "0", i.e. those
-                    // that are not subject to Dirichlet
-                    // conditions.
+                                      // Here we add the contributions from
+                                      // Neumann (pressure) boundary
+                                      // conditions at faces on the domain
+                                      // boundary that have the boundary flag
+                                      // "0", i.e. those that are not subject
+                                      // to Dirichlet conditions.
       for (unsigned int face_no=0;
            face_no<GeometryInfo<dim>::faces_per_cell;
            ++face_no)
@@ -784,16 +781,15 @@ void StokesProblem<dim>::assemble_system ()
                 }
           }
 
-                    // The final step is, as usual,
-                    // the transfer of the local
-                    // contributions to the global
-                    // system matrix. This works
-                    // also in the case of block
-                    // vectors and matrices, and
-                    // also the terms constituting
-                    // the pressure mass matrix are
-                    // written into the correct position
-                    // without any further interaction.
+                                      // The final step is, as usual, the
+                                      // transfer of the local contributions
+                                      // to the global system matrix. This
+                                      // works also in the case of block
+                                      // vectors and matrices, and also the
+                                      // terms constituting the pressure mass
+                                      // matrix are written into the correct
+                                      // position without any further
+                                      // interaction.
       cell->get_dof_indices (local_dof_indices);
 
       for (unsigned int i=0; i<dofs_per_cell; ++i)
@@ -806,30 +802,24 @@ void StokesProblem<dim>::assemble_system ()
         system_rhs(local_dof_indices[i]) += local_rhs(i);
     }
 
-                    // After the addition of the local
-                    // contributions, we have to
-                    // condense the hanging node
-                    // constraints and interpolate
-                    // Dirichlet boundary conditions.
-                    // Note that Dirichlet boundary
-                    // conditions are only condensed
-                    // in boundary points that are 
-                    // labeled with "1", indicating
-                    // that Dirichlet data is to be
-                    // set.
-                    // There is one more thing, though.
-                    // The function describing the
-                    // Dirichlet conditions was 
-                    // defined for all components, both
-                    // velocity and pressure. However, 
-                    // the Dirichlet conditions are to
-                    // be set for the velocity only.
-                    // To this end, we use a 
-                    // <code>component_mask</code> that
-                    // filters away the pressure 
-                    // component, so that the condensation
-                    // is performed on
-                    // velocity dofs.
+                                  // After the addition of the local
+                                  // contributions, we have to condense the
+                                  // hanging node constraints and interpolate
+                                  // Dirichlet boundary conditions.  Note
+                                  // that Dirichlet boundary conditions are
+                                  // only condensed in boundary points that
+                                  // are labeled with "1", indicating that
+                                  // Dirichlet data is to be set.  There is
+                                  // one more thing, though.  The function
+                                  // describing the Dirichlet conditions was
+                                  // defined for all components, both
+                                  // velocity and pressure. However, the
+                                  // Dirichlet conditions are to be set for
+                                  // the velocity only.  To this end, we use
+                                  // a <code>component_mask</code> that
+                                  // filters away the pressure component, so
+                                  // that the condensation is performed on
+                                  // velocity dofs.
   hanging_node_constraints.condense (system_matrix);
   hanging_node_constraints.condense (system_rhs);  
 
@@ -849,48 +839,43 @@ void StokesProblem<dim>::assemble_system ()
                                         system_rhs);
   }
   
-                    // Before we're going to solve 
-                    // this linear system, we generate
-                    // a preconditioner for the
-                    // velocity-velocity matrix, 
-                    // i.e., <code>block(0,0)</code>
-                    // in the system matrix. As mentioned
-                    // above, this depends on the
-                    // spatial dimension. Since this
-                    // handled automatically by the
-                    // template <code><dim></code>
-                    // in <code>InnerPreconditioner</code>,
-                    // we don't have to manually
-                    // intervene at this point any 
-                    // further.
+                                  // Before we're going to solve this linear
+                                  // system, we generate a preconditioner for
+                                  // the velocity-velocity matrix, i.e.,
+                                  // <code>block(0,0)</code> in the system
+                                  // matrix. As mentioned above, this depends
+                                  // on the spatial dimension. Since this
+                                  // handled automatically by the template
+                                  // <code><dim></code> in
+                                  // <code>InnerPreconditioner</code>, we
+                                  // don't have to manually intervene at this
+                                  // point any further.
   std::cout << "   Computing preconditioner..." << std::endl << std::flush;
       
   A_preconditioner
     = boost::shared_ptr<typename InnerPreconditioner<dim>::type>(new typename InnerPreconditioner<dim>::type());
   A_preconditioner->initialize (system_matrix.block(0,0),
-                typename InnerPreconditioner<dim>::type::AdditionalData());
+                               typename InnerPreconditioner<dim>::type::AdditionalData());
 
 }
 
 
 
-                        // @sect4{StokesProblem::solve}
+                                // @sect4{StokesProblem::solve}
                         
-                    // After the discussion in the
-                    // introduction and the definition
-                    // of the respective classes above,
-                    // the implementation of the 
-                    // <code>solve</code> function is
-                    // rather straigt-forward and done in
-                    // a similar way as in step-20. To 
-                    // start with, we need an object of
-                    // the <code>InverseMatrix</code> class
-                    // that represents the inverse of 
-                    // the matrix A. As described in
-                    // the introduction, the inverse
-                    // is generated with the help
-                    // of an inner preconditioner of
-                    // type <code>InnerPreconditioner<dim></code>.
+                                // After the discussion in the introduction
+                                // and the definition of the respective
+                                // classes above, the implementation of the
+                                // <code>solve</code> function is rather
+                                // straigt-forward and done in a similar way
+                                // as in step-20. To start with, we need an
+                                // object of the <code>InverseMatrix</code>
+                                // class that represents the inverse of the
+                                // matrix A. As described in the
+                                // introduction, the inverse is generated
+                                // with the help of an inner preconditioner
+                                // of type
+                                // <code>InnerPreconditioner<dim></code>.
 template <int dim>
 void StokesProblem<dim>::solve () 
 {
@@ -899,16 +884,14 @@ void StokesProblem<dim>::solve ()
     A_inverse (system_matrix.block(0,0), *A_preconditioner);
   Vector<double> tmp (solution.block(0).size());
   
-                    // This is as in step-20. We generate
-                    // the right hand side 
-                    // B A^{-1} F  Ð  G for the
-                    // Schur complement and an object
-                    // that represents the respective
-                    // linear operation B A^{-1} B^T,
-                    // now with a template parameter
-                    // indicating the preconditioner -
-                    // in accordance with the definition
-                    // of the class.
+                                  // This is as in step-20. We generate the
+                                  // right hand side $B A^{-1} F - G$ for the
+                                  // Schur complement and an object that
+                                  // represents the respective linear
+                                  // operation $B A^{-1} B^T$, now with a
+                                  // template parameter indicating the
+                                  // preconditioner - in accordance with the
+                                  // definition of the class.
   {
     Vector<double> schur_rhs (solution.block(1).size());
     A_inverse.vmult (tmp, system_rhs.block(0));
@@ -918,90 +901,79 @@ void StokesProblem<dim>::solve ()
     SchurComplement<typename InnerPreconditioner<dim>::type>
       schur_complement (system_matrix, A_inverse);
     
-                    // The usual control structures for
-                    // the solver call are created...
+                                    // The usual control structures for
+                                    // the solver call are created...
     SolverControl solver_control (system_matrix.block(0,0).m(),
-                                  1e-6*schur_rhs.l2_norm());
+                                 1e-6*schur_rhs.l2_norm());
     SolverCG<>    cg (solver_control);
     
-                    // Now to the preconditioner to the
-                    // Schur complement. As explained in the
-                    // introduction, the preconditioning
-                    // is done by a mass matrix in the
-                    // pressure variable. 
-                    // It is stored in the (1,1) block
-                    // of the system matrix (that is not
-                    // used elsewhere in this function).
-                    
-                    // Actually, the solver needs to have
-                    // the preconditioner in the form
-                    // P^{-1}, so we need to create 
-                    // an inverse operation. Once again,
-                    // we use an object of the class
-                    // <code>InverseMatrix</code>, which
-                    // implements the <code>vmult</code>
-                    // operation that is needed by the solver.
-                    // In this case, we have to invert
-                    // the pressure mass matrix. As it 
-                    // already turned out in earlier tutorial
-                    // programs, the inversion of a mass
-                    // matrix is a rather cheap and
-                    // straight-forward operation (compared
-                    // to, e.g., a Laplace matrix). The CG
-                    // method with SSOR preconditioning
-                    // converges in 10-20 steps,
-                    // independently on the mesh size.
-                    // This is precisely what we do here:
-                    // We choose an SSOR preconditioner
-                    // with parameter 1.2 and take it along
-                    // to the InverseMatrix object via
-                    // the corresponding template parameter.
-                    // A CG solver is then called within
-                    // the vmult operation of the inverse
-                    // matrix.
+                                    // Now to the preconditioner to the Schur
+                                    // complement. As explained in the
+                                    // introduction, the preconditioning is
+                                    // done by a mass matrix in the pressure
+                                    // variable.  It is stored in the $(1,1)$
+                                    // block of the system matrix (that is
+                                    // not used elsewhere in this function).
+                                    //
+                                    // Actually, the solver needs to have the
+                                    // preconditioner in the form $P^{-1}$, so
+                                    // we need to create an inverse
+                                    // operation. Once again, we use an
+                                    // object of the class
+                                    // <code>InverseMatrix</code>, which
+                                    // implements the <code>vmult</code>
+                                    // operation that is needed by the
+                                    // solver.  In this case, we have to
+                                    // invert the pressure mass matrix. As it
+                                    // already turned out in earlier tutorial
+                                    // programs, the inversion of a mass
+                                    // matrix is a rather cheap and
+                                    // straight-forward operation (compared
+                                    // to, e.g., a Laplace matrix). The CG
+                                    // method with SSOR preconditioning
+                                    // converges in 10-20 steps,
+                                    // independently on the mesh size.  This
+                                    // is precisely what we do here: We
+                                    // choose an SSOR preconditioner with
+                                    // parameter 1.2 and take it along to the
+                                    // InverseMatrix object via the
+                                    // corresponding template parameter.  A
+                                    // CG solver is then called within the
+                                    // vmult operation of the inverse matrix.
     PreconditionSSOR<> preconditioner;
     preconditioner.initialize (system_matrix.block(1,1), 1.2);
   
     InverseMatrix<SparseMatrix<double>,PreconditionSSOR<> >
       m_inverse (system_matrix.block(1,1), preconditioner);
     
-                    // With the Schur complement and an
-                    // efficient preconditioner at hand,
-                    // we can solve the respective
-                    // equation in the usual way.
-    try
-      {
-         cg.solve (schur_complement, solution.block(1), schur_rhs,
-                   m_inverse);
-      }
-    catch (...)
-      {
-         abort ();
-      }
+                                    // With the Schur complement and an
+                                    // efficient preconditioner at hand,
+                                    // we can solve the respective
+                                    // equation in the usual way.
+    cg.solve (schur_complement, solution.block(1), schur_rhs,
+             m_inverse);
   
-                     // After this first solution step,
-                     // the hanging node constraints have
-                     // to be distributed to the solution -
-                     // in order to achieve a consistent 
-                     // pressure field.
+                                    // After this first solution step,
+                                    // the hanging node constraints have
+                                    // to be distributed to the solution -
+                                    // in order to achieve a consistent 
+                                    // pressure field.
     hanging_node_constraints.distribute (solution);
   
     std::cout << "  "
-              << solver_control.last_step()
-              << " outer CG Schur complement iterations for pressure"
-              << std::flush
-              << std::endl;    
+             << solver_control.last_step()
+             << " outer CG Schur complement iterations for pressure"
+             << std::flush
+             << std::endl;    
   }
     
-                    // As in step-20, we finally need to
-                    // solve for the velocity equation
-                    // where we plug in the the solution 
-                    // to the pressure equation. This involves 
-                    // only objects we already know - so
-                    // we simply
-                    // multiply p by B^T, subtract the 
-                    // right hand side and multiply
-                    // by the inverse of A.
+                                  // As in step-20, we finally need to solve
+                                  // for the velocity equation where we plug
+                                  // in the the solution to the pressure
+                                  // equation. This involves only objects we
+                                  // already know - so we simply multiply p
+                                  // by $B^T$, subtract the right hand side and
+                                  // multiply by the inverse of A.
   {
     system_matrix.block(0,1).vmult (tmp, solution.block(1));
     tmp *= -1;
@@ -1009,34 +981,30 @@ void StokesProblem<dim>::solve ()
   
     A_inverse.vmult (solution.block(0), tmp);
   
-                     // Again, we need to distribute
-                     // the constraints from hanging nodes
-                     // in order to obtain a constistent
-                     // flow field.
+                                    // Again, we need to distribute the
+                                    // constraints from hanging nodes in
+                                    // order to obtain a constistent flow
+                                    // field.
     hanging_node_constraints.distribute (solution);
   }
 }
 
 
-                        // @sect4{StokesProblem::output_results}
+                                // @sect4{StokesProblem::output_results}
                         
-                    // The next function generates graphical
-                    // output. In this example, we are going
-                    // to use the VTK file format.
-                    // We attach names to the individual
-                    // variables in the problem - 
-                    // <code>velocity</code> to the dim
-                    // components of velocity and
-                    // <code>p</code> to the pressure.
-                    // In order to tell the VTK file 
-                    // which components are vectors
-                    // and which scalars, we need to
-                    // add that information as well -
-                    // achieved by the 
-                    // <code>DataComponentInterpretation</code>
-                    // class.
-                    // The rest of the function is 
-                    // then the same as in step-20.
+                                // The next function generates graphical
+                                // output. In this example, we are going to
+                                // use the VTK file format.  We attach names
+                                // to the individual variables in the problem
+                                // - <code>velocity</code> to the dim
+                                // components of velocity and <code>p</code>
+                                // to the pressure.  In order to tell the VTK
+                                // file which components are vectors and
+                                // which scalars, we need to add that
+                                // information as well - achieved by the
+                                // <code>DataComponentInterpretation</code>
+                                // class.  The rest of the function is then
+                                // the same as in step-20.
 template <int dim>
 void
 StokesProblem<dim>::output_results (const unsigned int refinement_cycle)  const
@@ -1056,8 +1024,8 @@ StokesProblem<dim>::output_results (const unsigned int refinement_cycle)  const
       = DataComponentInterpretation::component_is_part_of_vector;
   
   data_out.add_data_vector (solution, solution_names,
-                DataOut<dim>::type_dof_data,
-                data_component_interpretation);
+                           DataOut<dim>::type_dof_data,
+                           data_component_interpretation);
   
   data_out.build_patches ();
   
@@ -1071,21 +1039,21 @@ StokesProblem<dim>::output_results (const unsigned int refinement_cycle)  const
 }
 
 
-                        // @sect4{StokesProblem::refine_mesh}
+                                // @sect4{StokesProblem::refine_mesh}
                         
-                    // This is the last interesting function
-                    // of the <code>StokesProblem</code> class.
-                    // As indicated by its name, it takes the
-                    // solution to the problem and
-                    // refines the mesh where this is
-                    // needed. The procedure is the same
-                    // as in the respective step in
-                    // step-6, with the exception that
-                    // we base the refinement only on the
-                    // change in pressure, i.e., we call
-                    // the Kelly error estimator with a
-                    // mask object. Additionally, we do
-                    // not coarsen the grid again.
+                                // This is the last interesting function
+                                // of the <code>StokesProblem</code> class.
+                                // As indicated by its name, it takes the
+                                // solution to the problem and
+                                // refines the mesh where this is
+                                // needed. The procedure is the same
+                                // as in the respective step in
+                                // step-6, with the exception that
+                                // we base the refinement only on the
+                                // change in pressure, i.e., we call
+                                // the Kelly error estimator with a
+                                // mask object. Additionally, we do
+                                // not coarsen the grid again.
 template <int dim>
 void
 StokesProblem<dim>::refine_mesh () 
@@ -1108,23 +1076,25 @@ StokesProblem<dim>::refine_mesh ()
 }
 
 
-                        // @sect4{StokesProblem::run}
+                                // @sect4{StokesProblem::run}
                         
-                    // The last step in the Stokes class
-                    // is, as usual, the program that generates
-                    // the initial grid and calls the other
-                    // functions in the respective order.
+                                // The last step in the Stokes class
+                                // is, as usual, the program that generates
+                                // the initial grid and calls the other
+                                // functions in the respective order.
 template <int dim>
 void StokesProblem<dim>::run () 
 {
-                    // We start off with a rectangle of
-                    // size 4 x 1 (x 1), placed in R^2/R^3
-                    // as (-2,2)x(-1,0) or (-2,2)x(0,1)x(-1,1),
-                    // respectively. It is natural to start 
-                    // with equal mesh size in each direction,
-                    // so we subdivide the initial rectangle
-                    // four times in the first coordinate
-                    // direction.
+                                  // We start off with a rectangle of size $4
+                                  // \times 1$ (in 2d) or $4 \times 1 times
+                                  // 1$ (in 3d), placed in $R^2/R^3$ as
+                                  // $(-2,2)times(-1,0)$ or
+                                  // $(-2,2)\times(0,1)\times(-1,1)$,
+                                  // respectively. It is natural to start
+                                  // with equal mesh size in each direction,
+                                  // so we subdivide the initial rectangle
+                                  // four times in the first coordinate
+                                  // direction.
   std::vector<unsigned int> subdivisions (dim, 1);
   subdivisions[0] = 4;
     
@@ -1137,38 +1107,37 @@ void StokesProblem<dim>::run ()
                                               Point<dim>(2,0) :
                                               Point<dim>(2,1,0)));
   
-                    // A boundary indicator is set to all 
-                    // boundaries that are subject to 
-                    // Dirichlet boundary conditions, i.e.
-                    // to faces that are located at 0 in
-                    // the last coordinate direction. See
-                    // the example description above for
-                    // details.
+                                  // A boundary indicator is set to all
+                                  // boundaries that are subject to Dirichlet
+                                  // boundary conditions, i.e.  to faces that
+                                  // are located at 0 in the last coordinate
+                                  // direction. See the example description
+                                  // above for details.
   for (typename Triangulation<dim>::active_cell_iterator
-       cell = triangulation.begin_active();
+        cell = triangulation.begin_active();
        cell != triangulation.end(); ++cell)
     for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
       if (cell->face(f)->center()[dim-1] == 0)
-      {
-        cell->face(f)->set_boundary_indicator(1);
-        /*
-        for (unsigned int e=0; e<GeometryInfo<dim>::lines_per_face; ++e)
-          cell->face(f)->line(e)->set_boundary_indicator (1);
-        */
-      }
+       {
+         cell->face(f)->set_boundary_indicator(1);
+                                          /*
+                                            for (unsigned int e=0; e<GeometryInfo<dim>::lines_per_face; ++e)
+                                            cell->face(f)->line(e)->set_boundary_indicator (1);
+                                          */
+       }
   
   
-                    // We employ an initial refinement before
-                    // solving for the first time. In 3D,
-                    // there are going to be more dofs, so
-                    // we refine less there.
+                                  // We employ an initial refinement before
+                                  // solving for the first time. In 3D, there
+                                  // are going to be more dofs, so we refine
+                                  // less there.
   triangulation.refine_global (4-dim);
 
-                    // As first seen in step-6, we cycle
-                    // over the different refinement levels
-                    // and refine (if not the first step),
-                    // setup the dofs and matrices, assemble,
-                    // solve and create an output.
+                                  // As first seen in step-6, we cycle over
+                                  // the different refinement levels and
+                                  // refine (if not the first step), setup
+                                  // the dofs and matrices, assemble, solve
+                                  // and create an output.
   for (unsigned int refinement_cycle = 0; refinement_cycle<7;
        ++refinement_cycle)
     {
@@ -1192,13 +1161,12 @@ void StokesProblem<dim>::run ()
 }
 
 
-                        // @sect3{The <code>main</code> function}
+                                // @sect3{The <code>main</code> function}
 
-                    // The main function is the same as 
-                    // in step-20. We pass the element
-                    // degree as a parameter and
-                    // choose the space dimension at the
-                    // well-known template slot.
+                                // The main function is the same as in
+                                // step-20. We pass the element degree as a
+                                // parameter and choose the space dimension
+                                // at the well-known template slot.
 int main () 
 {
   try

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.