]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Use homogenous boundary values, to avoid an assertion in the
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 6 Jan 2010 10:34:24 +0000 (10:34 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 6 Jan 2010 10:34:24 +0000 (10:34 +0000)
constraint matrix.

git-svn-id: https://svn.dealii.org/trunk@20309 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 794db3544891fc4485419a3780faa7375f8a38a6..3d97a248d69882d3e0a0b4181b0d53da76123e00 100644 (file)
@@ -12,8 +12,8 @@
 
 
                                 // @sect3{Include files}
-                        
-                                // As usual, we start by including 
+
+                                // As usual, we start by including
                                 // some well-known files:
 #include <base/quadrature_lib.h>
 #include <base/logstream.h>
@@ -51,7 +51,7 @@
 #include <numerics/data_out.h>
 #include <numerics/error_estimator.h>
 
-                               
+
 #include <lac/sparse_direct.h>
 
 #include <lac/sparse_ilu.h>
 #include <multigrid/mg_coarse.h>
 #include <multigrid/mg_smoother.h>
 #include <multigrid/mg_matrix.h>
-                               
+
 #include <fstream>
 #include <sstream>
 
-                               
+
 using namespace dealii;
 
-                       
+
 template <int dim>
 struct InnerPreconditioner;
 
-                       
+
 template <>
-struct InnerPreconditioner<2> 
+struct InnerPreconditioner<2>
 {
     typedef SparseDirectUMFPACK type;
 };
 
-                               
+
 template <>
-struct InnerPreconditioner<3> 
+struct InnerPreconditioner<3>
 {
     typedef SparseILU<double> type;
 };
 
 
 template <int dim>
-class StokesProblem 
+class StokesProblem
 {
   public:
     StokesProblem (const unsigned int degree);
     void run ();
-    
+
   private:
     void setup_dofs ();
     void assemble_system ();
@@ -105,19 +105,19 @@ class StokesProblem
     void solve_block ();
     void solve_mg ();
 
-    void find_dofs_on_lower_level (std::vector<std::vector<bool> > &lower_dofs, 
+    void find_dofs_on_lower_level (std::vector<std::vector<bool> > &lower_dofs,
         std::vector<std::vector<bool> > &boundary_dofs);
     void output_results (const unsigned int refinement_cycle) const;
     void refine_mesh ();
-    
+
     const unsigned int   degree;
-    
+
     Triangulation<dim>   triangulation;
     FESystem<dim>        fe;
     MGDoFHandler<dim>    dof_handler;
 
     ConstraintMatrix     constraints;
-    
+
     BlockSparsityPattern      sparsity_pattern;
     BlockSparseMatrix<double> system_matrix;
 
@@ -130,22 +130,22 @@ class StokesProblem
 
     MGLevelObject<BlockSparseMatrix<double> > mg_interface_matrices;
     std::vector<std::vector<unsigned int> >   mg_dofs_per_component;
-                                 
+
     std_cxx1x::shared_ptr<typename InnerPreconditioner<dim>::type> A_preconditioner;
     std::vector<std_cxx1x::shared_ptr<typename InnerPreconditioner<dim>::type> > mg_A_preconditioner;
 };
 
-                       
+
 template <int dim>
-class BoundaryValues : public Function<dim> 
+class BoundaryValues : public Function<dim>
 {
   public:
     BoundaryValues () : Function<dim>(dim+1) {}
-    
+
     virtual double value (const Point<dim>   &p,
                           const unsigned int  component = 0) const;
 
-    virtual void vector_value (const Point<dim> &p, 
+    virtual void vector_value (const Point<dim> &p,
                                Vector<double>   &value) const;
 };
 
@@ -153,11 +153,11 @@ 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
 {
   Assert (component < this->n_components,
          ExcIndexRange (component, 0, this->n_components));
-  
+
    if (component == 0 && p[0] == 0)
     return (dim == 2 ? - p[1]*(p[1]-1.) : p[1]*(p[1]-1.) * p[2]*(p[2]-1.));
   return 0;
@@ -167,7 +167,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);
@@ -175,43 +175,43 @@ BoundaryValues<dim>::vector_value (const Point<dim> &p,
 
 
 
-                       
+
 template <int dim>
-class RightHandSide : public Function<dim> 
+class RightHandSide : public Function<dim>
 {
   public:
     RightHandSide () : Function<dim>(dim+1) {}
-    
+
     virtual double value (const Point<dim>   &p,
                           const unsigned int  component = 0) const;
 
-    virtual void vector_value (const Point<dim> &p, 
+    virtual void vector_value (const Point<dim> &p,
                                Vector<double>   &value) const;
-    
+
 };
 
 
 template <int dim>
 double
 RightHandSide<dim>::value (const Point<dim>  &/*p*/,
-                           const unsigned int /*component*/) const 
+                           const unsigned int component) const
 {
-  return 0;
+  return (component == 0 ? 1 : 0);
 }
 
 
 template <int dim>
 void
 RightHandSide<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) = RightHandSide<dim>::value (p, c);
 }
 
 
-                       
-                               
+
+
 template <class Matrix, class Preconditioner>
 class InverseMatrix : public Subscriptor
 {
@@ -264,10 +264,10 @@ class BlockSchurPreconditioner : public Subscriptor
 
   private:
     const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
-    const SmartPointer<const InverseMatrix<SparseMatrix<double>, 
+    const SmartPointer<const InverseMatrix<SparseMatrix<double>,
                        PreconditionerMp > > m_inverse;
     const PreconditionerA &a_preconditioner;
-    
+
     mutable Vector<double> tmp;
 
 };
@@ -286,7 +286,7 @@ BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::BlockSchurPrecondit
 {}
 
         // Now the interesting function, the multiplication of
-        // the preconditioner with a BlockVector. 
+        // the preconditioner with a BlockVector.
 template <class PreconditionerA, class PreconditionerMp>
 void BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::vmult (
                                      BlockVector<double>       &dst,
@@ -294,13 +294,13 @@ void BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::vmult (
 {
         // Form u_new = A^{-1} u
   a_preconditioner.vmult (dst.block(0), src.block(0));
-        // Form tmp = - B u_new + p 
+        // Form tmp = - B u_new + p
         // (<code>SparseMatrix::residual</code>
         // does precisely this)
   system_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1));
         // Change sign in tmp
   tmp *= -1;
-        // Multiply by approximate Schur complement 
+        // Multiply by approximate Schur complement
         // (i.e. a pressure mass matrix)
   m_inverse->vmult (dst.block(1), tmp);
 }
@@ -319,7 +319,7 @@ class SchurComplement : public Subscriptor
   private:
     const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
     const SmartPointer<const InverseMatrix<SparseMatrix<double>, Preconditioner> > A_inverse;
-    
+
     mutable Vector<double> tmp1, tmp2;
 };
 
@@ -347,7 +347,7 @@ void SchurComplement<Preconditioner>::vmult (Vector<double>       &dst,
 }
 
 
-                       
+
 template <int dim>
 StokesProblem<dim>::StokesProblem (const unsigned int degree)
                 :
@@ -359,31 +359,31 @@ StokesProblem<dim>::StokesProblem (const unsigned int degree)
 {}
 
 
-                               
-                       
-                               
+
+
+
 template <int dim>
 void StokesProblem<dim>::setup_dofs ()
 {
   A_preconditioner.reset ();
   mg_A_preconditioner.resize (0);
   system_matrix.clear ();
-  
-  dof_handler.distribute_dofs (fe);  
+
+  dof_handler.distribute_dofs (fe);
 //  DoFRenumbering::Cuthill_McKee (dof_handler);
 
   std::vector<unsigned int> block_component (dim+1,0);
   block_component[dim] = 1;
   DoFRenumbering::component_wise (dof_handler, block_component);
 
-                               
+
   {
     constraints.clear ();
     std::vector<bool> component_mask (dim+1, true);
     component_mask[dim] = false;
     VectorTools::interpolate_boundary_values (dof_handler,
                                              0,
-                                             BoundaryValues<dim>(),
+                                             ZeroFunction<dim>(dim+1), //TODO: go back to BoundaryValues
                                              constraints,
                                              component_mask);
     DoFTools::make_hanging_node_constraints (dof_handler,
@@ -392,10 +392,10 @@ void StokesProblem<dim>::setup_dofs ()
 
   constraints.close ();
 
-                               
+
   std::vector<unsigned int> dofs_per_block (2);
   DoFTools::count_dofs_per_block (dof_handler, dofs_per_block,
-                                 block_component);  
+                                 block_component);
   const unsigned int n_u = dofs_per_block[0],
                      n_p = dofs_per_block[1];
 
@@ -406,10 +406,10 @@ void StokesProblem<dim>::setup_dofs ()
             << dof_handler.n_dofs()
             << " (" << n_u << '+' << n_p << ')'
             << std::endl;
-      
-                               
-                                
-                                
+
+
+
+
   {
     BlockCompressedSimpleSparsityPattern csp (2,2);
 
@@ -417,23 +417,23 @@ void StokesProblem<dim>::setup_dofs ()
     csp.block(1,0).reinit (n_p, n_u);
     csp.block(0,1).reinit (n_u, n_p);
     csp.block(1,1).reinit (n_p, n_p);
-  
-    csp.collect_sizes();    
-  
+
+    csp.collect_sizes();
+
     DoFTools::make_sparsity_pattern (
-        static_cast<const DoFHandler<dim>&>(dof_handler), 
+        static_cast<const DoFHandler<dim>&>(dof_handler),
         csp, constraints, false);
     sparsity_pattern.copy_from (csp);
   }
-  
-                                
+
+
   system_matrix.reinit (sparsity_pattern);
-                                   
+
   solution.reinit (2);
   solution.block(0).reinit (n_u);
   solution.block(1).reinit (n_p);
   solution.collect_sizes ();
-  
+
   system_rhs.reinit (2);
   system_rhs.block(0).reinit (n_u);
   system_rhs.block(1).reinit (n_p);
@@ -460,7 +460,7 @@ void StokesProblem<dim>::setup_dofs ()
 
   for (unsigned int level=0; level<nlevels; ++level)
   {
-    BlockCompressedSparsityPattern bcsp (mg_dofs_per_component[level], 
+    BlockCompressedSparsityPattern bcsp (mg_dofs_per_component[level],
         mg_dofs_per_component[level]);
     MGTools::make_sparsity_pattern(dof_handler, bcsp, level);
     mg_sparsity[level].copy_from (bcsp);
@@ -470,13 +470,13 @@ void StokesProblem<dim>::setup_dofs ()
 }
 
 
-                       
+
 template <int dim>
-void StokesProblem<dim>::assemble_system () 
+void StokesProblem<dim>::assemble_system ()
 {
   system_matrix=0;
   system_rhs=0;
-  
+
   QGauss<dim>   quadrature_formula(degree+2);
 
   FEValues<dim> fe_values (fe, quadrature_formula,
@@ -484,42 +484,42 @@ void StokesProblem<dim>::assemble_system ()
                            update_quadrature_points  |
                            update_JxW_values |
                            update_gradients);
-  
+
   const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
-  
+
   const unsigned int   n_q_points      = quadrature_formula.size();
 
   FullMatrix<double>   local_matrix (dofs_per_cell, dofs_per_cell);
   Vector<double>       local_rhs (dofs_per_cell);
 
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-  
+
   const RightHandSide<dim>          right_hand_side;
   std::vector<Vector<double> >      rhs_values (n_q_points,
                                                 Vector<double>(dim+1));
 
-                                
+
   const FEValuesExtractors::Vector velocities (0);
   const FEValuesExtractors::Scalar pressure (dim);
 
-                                
-                                 
+
+
   std::vector<Tensor<2,dim> >          phi_grads_u (dofs_per_cell);
   std::vector<double>                  div_phi_u   (dofs_per_cell);
   std::vector<double>                  phi_p       (dofs_per_cell);
-                                  
+
   typename MGDoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
   for (; cell!=endc; ++cell)
-    { 
+    {
       fe_values.reinit (cell);
       local_matrix = 0;
       local_rhs = 0;
-      
+
       right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
                                         rhs_values);
-      
+
       for (unsigned int q=0; q<n_q_points; ++q)
        {
          for (unsigned int k=0; k<dofs_per_cell; ++k)
@@ -537,30 +537,29 @@ void StokesProblem<dim>::assemble_system ()
                                        - div_phi_u[i] * phi_p[j]
                                        - phi_p[i] * div_phi_u[j]
                                        + phi_p[i] * phi_p[j])
-                                      * fe_values.JxW(q);     
-
+                                      * 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) * 
+             local_rhs(i) += fe_values.shape_value(i,q) *
                              rhs_values[q](component_i) *
                              fe_values.JxW(q);
            }
        }
 
-                                     
+
 
 
       cell->get_dof_indices (local_dof_indices);
       constraints.distribute_local_to_global (local_matrix, local_rhs,
-                                             local_dof_indices, 
+                                             local_dof_indices,
                                              system_matrix, system_rhs);
     }
-  
-                               
+
+
   std::cout << "   Computing preconditioner..." << std::endl << std::flush;
-      
+
   A_preconditioner
     = std_cxx1x::shared_ptr<typename InnerPreconditioner<dim>::type>(new typename InnerPreconditioner<dim>::type());
   A_preconditioner->initialize (system_matrix.block(0,0),
@@ -570,7 +569,7 @@ void StokesProblem<dim>::assemble_system ()
 
 
 template <int dim>
-void StokesProblem<dim>::assemble_multigrid () 
+void StokesProblem<dim>::assemble_multigrid ()
 {
   QGauss<dim>   quadrature_formula(degree+2);
   FEValues<dim> fe_values (fe, quadrature_formula,
@@ -578,22 +577,22 @@ void StokesProblem<dim>::assemble_multigrid ()
                            update_quadrature_points  |
                            update_JxW_values |
                            update_gradients);
-  
+
   const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
   const unsigned int   n_q_points      = quadrature_formula.size();
 
   FullMatrix<double>   local_matrix (dofs_per_cell, dofs_per_cell);
 
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-                                
+
   const FEValuesExtractors::Vector velocities (0);
   const FEValuesExtractors::Scalar pressure (dim);
-                                
-                                 
+
+
   std::vector<Tensor<2,dim> >          phi_grads_u (dofs_per_cell);
   std::vector<double>                  div_phi_u   (dofs_per_cell);
   std::vector<double>                  phi_p       (dofs_per_cell);
-                                  
+
   std::vector<std::vector<bool> > interface_dofs;
   std::vector<std::vector<bool> > boundary_interface_dofs;
   for (unsigned int level = 0; level<triangulation.n_levels(); ++level)
@@ -632,7 +631,7 @@ void StokesProblem<dim>::assemble_multigrid ()
     cell = dof_handler.begin(),
     endc = dof_handler.end();
   for (; cell!=endc; ++cell)
-    { 
+    {
                                       // Remember the level of the
                                       // current cell.
       const unsigned int level = cell->level();
@@ -640,7 +639,7 @@ void StokesProblem<dim>::assemble_multigrid ()
                                       // by update flags above.
       fe_values.reinit (cell);
       local_matrix = 0;
-      
+
       for (unsigned int q=0; q<n_q_points; ++q)
       {
          for (unsigned int k=0; k<dofs_per_cell; ++k)
@@ -657,17 +656,17 @@ void StokesProblem<dim>::assemble_multigrid ()
                                        - div_phi_u[i] * phi_p[j]
                                        - phi_p[i] * div_phi_u[j]
                                        + phi_p[i] * phi_p[j])
-                                      * fe_values.JxW(q);     
+                                      * fe_values.JxW(q);
 
       }
-                                     
+
           cell->get_mg_dof_indices (local_dof_indices);
 //       for (unsigned int i=0; i<dofs_per_cell; ++i)
 //         for (unsigned int j=0; j<dofs_per_cell; ++j)
 //           mg_matrices[level].add (local_dof_indices[i],
 //                                   local_dof_indices[j],
 //                                   local_matrix(i,j));
-         
+
          boundary_constraints[level]
             .distribute_local_to_global (local_matrix,
                 local_dof_indices,
@@ -675,7 +674,7 @@ void StokesProblem<dim>::assemble_multigrid ()
 
     for (unsigned int i=0; i<dofs_per_cell; ++i)
       for (unsigned int j=0; j<dofs_per_cell; ++j)
-       if( !(interface_dofs[level][local_dof_indices[i]]==true && 
+       if( !(interface_dofs[level][local_dof_indices[i]]==true &&
              interface_dofs[level][local_dof_indices[j]]==false))
          local_matrix(i,j) = 0;
 
@@ -696,9 +695,9 @@ void StokesProblem<dim>::assemble_multigrid ()
     }
 }
 
-                               
+
 template <int dim>
-void StokesProblem<dim>::solve_block () 
+void StokesProblem<dim>::solve_block ()
 {
   SparseMatrix<double> pressure_mass_matrix;
   pressure_mass_matrix.reinit(sparsity_pattern.block(1,1));
@@ -706,14 +705,14 @@ void StokesProblem<dim>::solve_block ()
   system_matrix.block(1,1) = 0;
 
   SparseILU<double> pmass_preconditioner;
-  pmass_preconditioner.initialize (pressure_mass_matrix, 
+  pmass_preconditioner.initialize (pressure_mass_matrix,
       SparseILU<double>::AdditionalData());
 
   InverseMatrix<SparseMatrix<double>,SparseILU<double> >
     m_inverse (pressure_mass_matrix, pmass_preconditioner);
 
   BlockSchurPreconditioner<typename InnerPreconditioner<dim>::type,
-    SparseILU<double> > 
+    SparseILU<double> >
       preconditioner (system_matrix, m_inverse, *A_preconditioner);
 
   SolverControl solver_control (system_matrix.m(),
@@ -740,14 +739,14 @@ template <typename InnerPreconditioner>
 class SchurComplementSmoother
 {
   public:
-    struct AdditionalData 
+    struct AdditionalData
     {
        const InnerPreconditioner *A_preconditioner;
     };
 
     void initialize (const BlockSparseMatrix<double> &system_matrix,
                     const AdditionalData            &data);
-    
+
     void vmult (BlockVector<double> &dst,
                const BlockVector<double> &src) const;
 
@@ -755,7 +754,7 @@ class SchurComplementSmoother
                 const BlockVector<double> &src) const;
 
     void clear ();
-    
+
   private:
     SmartPointer<const BlockSparseMatrix<double> > system_matrix;
     SmartPointer<const InnerPreconditioner>        A_preconditioner;
@@ -784,51 +783,51 @@ vmult (BlockVector<double> &dst,
   const InverseMatrix<SparseMatrix<double>,InnerPreconditioner>
     A_inverse (system_matrix->block(0,0), *A_preconditioner);
   Vector<double> tmp (dst.block(0).size());
-  
-                                 
+
+
   {
     Vector<double> schur_rhs (dst.block(1).size());
     A_inverse.vmult (tmp, src.block(0));
     system_matrix->block(1,0).vmult (schur_rhs, tmp);
     schur_rhs -= src.block(1);
-  
+
     SchurComplement<InnerPreconditioner>
       schur_complement (*system_matrix, A_inverse);
-    
+
                                     // The usual control structures for
                                     // the solver call are created...
     SolverControl solver_control (dst.block(1).size(),
                                  1e-6*schur_rhs.l2_norm());
     SolverCG<>    cg (solver_control);
-    
-                                   
-                                  
+
+
+
     SparseILU<double> preconditioner;
-    preconditioner.initialize (system_matrix->block(1,1), 
+    preconditioner.initialize (system_matrix->block(1,1),
                               SparseILU<double>::AdditionalData());
-  
+
     InverseMatrix<SparseMatrix<double>,SparseILU<double> >
       m_inverse (system_matrix->block(1,1), preconditioner);
-    
-                                   
+
+
     cg.solve (schur_complement, dst.block(1), schur_rhs,
              m_inverse);
-  
+
 // no constraints to be taken care of here
-    
+
     std::cout << "  "
              << solver_control.last_step()
              << " CG Schur complement iterations in smoother"
              << std::flush
-             << std::endl;    
+             << std::endl;
   }
-    
-                                
+
+
   {
     system_matrix->block(0,1).vmult (tmp, dst.block(1));
     tmp *= -1;
     tmp += src.block(0);
-  
+
     A_inverse.vmult (dst.block(0), tmp);
 
 // no constraints here either
@@ -856,7 +855,7 @@ Tvmult (BlockVector<double> &,
 
 
 template <int dim>
-void StokesProblem<dim>::solve_mg () 
+void StokesProblem<dim>::solve_mg ()
 {
   assemble_multigrid ();
   typedef PreconditionMG<dim, BlockVector<double>, MGTransferPrebuilt<BlockVector<double> > >
@@ -881,7 +880,7 @@ void StokesProblem<dim>::solve_mg ()
   typedef
     SchurComplementSmoother<typename InnerPreconditioner<dim>::type>
     Smoother;
-    
+
   MGSmootherPrecondition<BlockSparseMatrix<double>,
     Smoother,
     BlockVector<double> >
@@ -930,58 +929,58 @@ void StokesProblem<dim>::solve_mg ()
 
 
 template <int dim>
-void StokesProblem<dim>::solve_schur () 
+void StokesProblem<dim>::solve_schur ()
 {
   const InverseMatrix<SparseMatrix<double>,
                       typename InnerPreconditioner<dim>::type>
     A_inverse (system_matrix.block(0,0), *A_preconditioner);
   Vector<double> tmp (solution.block(0).size());
-  
-                                 
+
+
   {
     Vector<double> schur_rhs (solution.block(1).size());
     A_inverse.vmult (tmp, system_rhs.block(0));
     system_matrix.block(1,0).vmult (schur_rhs, tmp);
     schur_rhs -= system_rhs.block(1);
-  
+
     SchurComplement<typename InnerPreconditioner<dim>::type>
       schur_complement (system_matrix, A_inverse);
-    
+
                                     // The usual control structures for
                                     // the solver call are created...
     SolverControl solver_control (solution.block(1).size(),
                                  1e-6*schur_rhs.l2_norm());
     SolverCG<>    cg (solver_control);
-    
-                                   
-                                  
+
+
+
     SparseILU<double> preconditioner;
-    preconditioner.initialize (system_matrix.block(1,1), 
+    preconditioner.initialize (system_matrix.block(1,1),
       SparseILU<double>::AdditionalData());
-  
+
     InverseMatrix<SparseMatrix<double>,SparseILU<double> >
       m_inverse (system_matrix.block(1,1), preconditioner);
-    
-                                   
+
+
     cg.solve (schur_complement, solution.block(1), schur_rhs,
              m_inverse);
-  
-                                  
+
+
     constraints.distribute (solution);
-  
+
     std::cout << "  "
              << solver_control.last_step()
              << " outer CG Schur complement iterations for pressure"
              << std::flush
-             << std::endl;    
+             << std::endl;
   }
-    
-                                
+
+
   {
     system_matrix.block(0,1).vmult (tmp, solution.block(1));
     tmp *= -1;
     tmp += system_rhs.block(0);
-  
+
     A_inverse.vmult (solution.block(0), tmp);
 
     constraints.distribute (solution);
@@ -1010,14 +1009,14 @@ void StokesProblem<dim>::find_dofs_on_lower_level (std::vector<std::vector<bool>
     std::vector<bool> boundary_cell_dofs(dofs_per_cell);
     const unsigned int level = cell->level();
     cell->get_mg_dof_indices (local_dof_indices);
-    for (unsigned int face_nr=0; 
+    for (unsigned int face_nr=0;
         face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
     {
       typename DoFHandler<dim>::face_iterator face = cell->face(face_nr);
       if(!cell->at_boundary(face_nr))
       {
         //interior face
-        typename MGDoFHandler<dim>::cell_iterator neighbor 
+        typename MGDoFHandler<dim>::cell_iterator neighbor
           = cell->neighbor(face_nr);
         // Do refinement face
         // from the coarse side
@@ -1038,7 +1037,7 @@ void StokesProblem<dim>::find_dofs_on_lower_level (std::vector<std::vector<bool>
           if(!cell->at_boundary(neighbor_face_nr))
           {
             //other face is interior
-            typename MGDoFHandler<dim>::cell_iterator neighbor 
+            typename MGDoFHandler<dim>::cell_iterator neighbor
               = cell->neighbor(neighbor_face_nr);
             // Do refinement face
             // from the coarse side
@@ -1063,27 +1062,27 @@ void StokesProblem<dim>::find_dofs_on_lower_level (std::vector<std::vector<bool>
     }
   }
 }
-                       
+
 template <int dim>
 void
 StokesProblem<dim>::output_results (const unsigned int refinement_cycle)  const
 {
   std::vector<std::string> solution_names (dim, "velocity");
   solution_names.push_back ("pressure");
-  
+
   std::vector<DataComponentInterpretation::DataComponentInterpretation>
     data_component_interpretation
     (dim, DataComponentInterpretation::component_is_part_of_vector);
   data_component_interpretation
     .push_back (DataComponentInterpretation::component_is_scalar);
-      
+
   DataOut<dim> data_out;
-  data_out.attach_dof_handler (dof_handler);  
+  data_out.attach_dof_handler (dof_handler);
   data_out.add_data_vector (solution, solution_names,
                            DataOut<dim>::type_dof_data,
                            data_component_interpretation);
   data_out.build_patches ();
-  
+
   std::ostringstream filename;
   filename << "solution-"
            << Utilities::int_to_string (refinement_cycle, 2)
@@ -1094,10 +1093,10 @@ StokesProblem<dim>::output_results (const unsigned int refinement_cycle)  const
 }
 
 
-                       
+
 template <int dim>
 void
-StokesProblem<dim>::refine_mesh () 
+StokesProblem<dim>::refine_mesh ()
 {
   Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
 
@@ -1117,9 +1116,9 @@ StokesProblem<dim>::refine_mesh ()
 }
 
 
-                               
+
 template <int dim>
-void StokesProblem<dim>::run () 
+void StokesProblem<dim>::run ()
 {
   {
     std::vector<unsigned int> subdivisions (dim, 1);
@@ -1131,42 +1130,42 @@ void StokesProblem<dim>::run ()
     const Point<dim> top_right   = (dim == 2 ?
                                    Point<dim>(1,1) :
                                    Point<dim>(1,1,1));
-    
+
     GridGenerator::subdivided_hyper_rectangle (triangulation,
                                               subdivisions,
                                               bottom_left,
                                               top_right);
   }
-  
-                                
+
+
   for (typename Triangulation<dim>::active_cell_iterator
         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()[0] == 1)
        cell->face(f)->set_all_boundary_indicators(1);
-  
-  
-                               
+
+
+
   triangulation.refine_global (4-dim);
 
-                                 
+
   for (unsigned int refinement_cycle = 0; refinement_cycle<6;
        ++refinement_cycle)
     {
       std::cout << "Refinement cycle " << refinement_cycle << std::endl;
-      
+
       if (refinement_cycle > 0)
         refine_mesh ();
-      
+
       setup_dofs ();
 
       std::cout << "   Assembling..." << std::endl << std::flush;
-      assemble_system ();      
+      assemble_system ();
 
       std::cout << "   Solving..." << std::flush;
       solve_mg ();
-      
+
       output_results (refinement_cycle);
 
       std::cout << std::endl;
@@ -1174,12 +1173,12 @@ void StokesProblem<dim>::run ()
 }
 
 
-                       
-int main () 
+
+int main ()
 {
   try
     {
-//      deallog.depth_console (0);
+      deallog.depth_console (0);
 
       StokesProblem<2> flow_problem(1);
       flow_problem.run ();
@@ -1194,10 +1193,10 @@ int main ()
                 << "Aborting!" << std::endl
                 << "----------------------------------------------------"
                 << std::endl;
-      
+
       return 1;
     }
-  catch (...) 
+  catch (...)
     {
       std::cerr << std::endl << std::endl
                 << "----------------------------------------------------"

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.