]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Use two DoFHandlers, matrices, sets of vectors, etc, so that we don't always have...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 19 Aug 2008 03:48:14 +0000 (03:48 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 19 Aug 2008 03:48:14 +0000 (03:48 +0000)
git-svn-id: https://svn.dealii.org/trunk@16593 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 410b88704251c4230e3b2e61182afef06bd38321..b1ccf2d3dd7753dcbfbaef866bb3a705fa6ec242 100644 (file)
@@ -319,38 +319,50 @@ class BoussinesqFlowProblem
 
   private:
     void setup_dofs ();
-    void assemble_preconditioner ();
-    void assemble_system ();
-    void assemble_rhs_T ();
+    void assemble_stokes_preconditioner ();
+    void build_stokes_preconditioner ();
+    void assemble_stokes_system ();
+    void assemble_temperature_system ();
     double get_maximal_velocity () const;
     double get_maximal_temperature () const;
     void solve ();
     void output_results () const;
     void refine_mesh (const unsigned int max_grid_level);
 
-    const unsigned int   degree;
+    Triangulation<dim>        triangulation;
 
-    Triangulation<dim>   triangulation;
-    FESystem<dim>        fe;
-    DoFHandler<dim>      dof_handler;
+    const unsigned int        degree;
+    
+    FESystem<dim>             stokes_fe;
+    DoFHandler<dim>           stokes_dof_handler;
+    ConstraintMatrix          stokes_constraints;
+    
+    BlockSparsityPattern      stokes_sparsity_pattern;
+    BlockSparseMatrix<double> stokes_matrix;
+    BlockSparsityPattern      stokes_preconditioner_sparsity_pattern;
+    BlockSparseMatrix<double> stokes_preconditioner_matrix;
+
+    BlockVector<double>       stokes_solution;
+    BlockVector<double>       stokes_rhs;
 
-    ConstraintMatrix     hanging_node_constraints;
+    
+    FE_Q<dim>                 temperature_fe;
+    DoFHandler<dim>           temperature_dof_handler;
+    ConstraintMatrix          temperature_constraints;
+    
+    SparsityPattern           temperature_sparsity_pattern;
+    SparseMatrix<double>      temperature_matrix;
 
-    BlockSparsityPattern      sparsity_pattern;
-    BlockSparseMatrix<double> system_matrix;
+    Vector<double>            temperature_solution;
+    Vector<double>            old_temperature_solution;
+    Vector<double>            old_old_temperature_solution;
+    Vector<double>            temperature_rhs;
 
-    BlockSparsityPattern      preconditioner_sparsity_pattern;
-    BlockSparseMatrix<double> preconditioner_matrix;
 
     double time_step;
     double old_time_step;
     unsigned int timestep_number;
 
-    BlockVector<double> solution;
-    BlockVector<double> old_solution;
-    BlockVector<double> old_old_solution;
-    BlockVector<double> system_rhs;
-
     boost::shared_ptr<PreconditionerTrilinosAmg>  Amg_preconditioner;
     boost::shared_ptr<SparseILU<double> > Mp_preconditioner;
 
@@ -441,10 +453,10 @@ TemperatureBoundaryValues<dim>::value (const Point<dim> &p,
 
                                 // @sect4{Initial values}
 template <int dim>
-class InitialValues : public Function<dim>
+class TemperatureInitialValues : public Function<dim>
 {
   public:
-    InitialValues () : Function<dim>(dim+2) {}
+    TemperatureInitialValues () : Function<dim>(1) {}
 
     virtual double value (const Point<dim>   &p,
                           const unsigned int  component = 0) const;
@@ -456,7 +468,7 @@ class InitialValues : public Function<dim>
 
 template <int dim>
 double
-InitialValues<dim>::value (const Point<dim>  &,
+TemperatureInitialValues<dim>::value (const Point<dim>  &,
                            const unsigned int) const
 {
   return 0;
@@ -465,11 +477,11 @@ InitialValues<dim>::value (const Point<dim>  &,
 
 template <int dim>
 void
-InitialValues<dim>::vector_value (const Point<dim> &p,
+TemperatureInitialValues<dim>::vector_value (const Point<dim> &p,
                                   Vector<double>   &values) const
 {
   for (unsigned int c=0; c<this->n_components; ++c)
-    values(c) = InitialValues<dim>::value (p, c);
+    values(c) = TemperatureInitialValues<dim>::value (p, c);
 }
 
 
@@ -488,10 +500,10 @@ InitialValues<dim>::vector_value (const Point<dim> &p,
                                 // in the problem description in the
                                 // introduction.
 template <int dim>
-class RightHandSide : public Function<dim>
+class TemperatureRightHandSide : public Function<dim>
 {
   public:
-    RightHandSide () : Function<dim>(dim+2) {}
+    TemperatureRightHandSide () : Function<dim>(1) {}
 
     virtual double value (const Point<dim>   &p,
                           const unsigned int  component = 0) const;
@@ -503,42 +515,35 @@ class RightHandSide : public Function<dim>
 
 template <int dim>
 double
-RightHandSide<dim>::value (const Point<dim>  &p,
-                           const unsigned int component) const
+TemperatureRightHandSide<dim>::value (const Point<dim>  &p,
+                                     const unsigned int /*component*/) const
 {
-                                  // the right hand side only acts on the
-                                  // temperature component
-  if (component == dim+1)
-    {
-      static const Point<dim> source_centers[3]
-       = { (dim == 2 ? Point<dim>(.3,.1) : Point<dim>(.3,.5,.1)),
-           (dim == 2 ? Point<dim>(.45,.1) : Point<dim>(.45,.5,.1)),
-           (dim == 2 ? Point<dim>(.75,.1) : Point<dim>(.75,.5,.1)) };
-      static const double source_radius
-       = (dim == 2 ? 1./32 : 1./8);
+  static const Point<dim> source_centers[3]
+    = { (dim == 2 ? Point<dim>(.3,.1) : Point<dim>(.3,.5,.1)),
+       (dim == 2 ? Point<dim>(.45,.1) : Point<dim>(.45,.5,.1)),
+       (dim == 2 ? Point<dim>(.75,.1) : Point<dim>(.75,.5,.1)) };
+  static const double source_radius
+    = (dim == 2 ? 1./32 : 1./8);
       
-      return ((source_centers[0].distance (p) < source_radius)
-             ||
-             (source_centers[1].distance (p) < source_radius)
-             ||
-             (source_centers[2].distance (p) < source_radius)
-             ?
-             1
-             :
-             0);
-    }
-  else
-    return 0;
+  return ((source_centers[0].distance (p) < source_radius)
+         ||
+         (source_centers[1].distance (p) < source_radius)
+         ||
+         (source_centers[2].distance (p) < source_radius)
+         ?
+         1
+         :
+         0);
 }
 
 
 template <int dim>
 void
-RightHandSide<dim>::vector_value (const Point<dim> &p,
-                                  Vector<double>   &values) const
+TemperatureRightHandSide<dim>::vector_value (const Point<dim> &p,
+                                            Vector<double>   &values) const
 {
   for (unsigned int c=0; c<this->n_components; ++c)
-    values(c) = RightHandSide<dim>::value (p, c);
+    values(c) = TemperatureRightHandSide<dim>::value (p, c);
 }
 
 
@@ -726,7 +731,7 @@ class BlockSchurPreconditioner : public Subscriptor
               const BlockVector<double> &src) const;
 
   private:
-    const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
+    const SmartPointer<const BlockSparseMatrix<double> > stokes_matrix;
     const SmartPointer<const InverseMatrix<SparseMatrix<double>,
                        PreconditionerMp > > m_inverse;
     const PreconditionerA &a_preconditioner;
@@ -742,7 +747,7 @@ BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::BlockSchurPrecondit
           const PreconditionerA &Apreconditioner
           )
                 :
-                system_matrix           (&S),
+                stokes_matrix           (&S),
                 m_inverse               (&Mpinv),
                 a_preconditioner        (Apreconditioner),
                 tmp                     (S.block(1,1).m())
@@ -775,7 +780,7 @@ void BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::vmult (
                                      const BlockVector<double> &src) const
 {
   a_preconditioner.vmult (dst.block(0), src.block(0));
-  system_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1));
+  stokes_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1));
   tmp *= -1;
   m_inverse->vmult (dst.block(1), tmp);
 }
@@ -801,12 +806,13 @@ void BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::vmult (
 template <int dim>
 BoussinesqFlowProblem<dim>::BoussinesqFlowProblem (const unsigned int degree)
                 :
-                degree (degree),
                triangulation (Triangulation<dim>::maximum_smoothing),
-                fe (FE_Q<dim>(degree+1), dim,
-                    FE_Q<dim>(degree), 1,
-                    FE_Q<dim>(degree), 1),
-                dof_handler (triangulation),
+                degree (degree),
+                stokes_fe (FE_Q<dim>(degree+1), dim,
+                          FE_Q<dim>(degree), 1),
+               stokes_dof_handler (triangulation),
+               temperature_fe (degree),
+                temperature_dof_handler (triangulation),
                 time_step (0),
                old_time_step (0),
                timestep_number (0),
@@ -860,32 +866,41 @@ BoussinesqFlowProblem<dim>::BoussinesqFlowProblem (const unsigned int degree)
 template <int dim>
 void BoussinesqFlowProblem<dim>::setup_dofs ()
 {
-  dof_handler.distribute_dofs (fe);
-  DoFRenumbering::Cuthill_McKee (dof_handler);
-  std::vector<unsigned int> block_component (dim+2,0);
-  block_component[dim] = 1;
-  block_component[dim+1] = 2;
-  DoFRenumbering::component_wise (dof_handler, block_component);
-
-                                // The next step is, as usual, 
-                                // to write some information
-                                // to the screen. The information
-                                // that is most interesting during
-                                // the calculations is the
-                                // number of degrees of freedom
-                                // in the individual components,
-                                // so we count them. The function 
-                                // to do this is the same as the
-                                // one used in step-22, which 
-                                // uses the grouping of all
-                                // velocity components into
-                                // one block as introduced
-                                // above.
-  std::vector<unsigned int> dofs_per_block (3);
-  DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);  
-  const unsigned int n_u = dofs_per_block[0],
-                     n_p = dofs_per_block[1],
-                    n_T = dofs_per_block[2];
+  std::vector<unsigned int> stokes_block_component (dim+1,0);
+  stokes_block_component[dim] = 1;
+  
+  {
+    stokes_dof_handler.distribute_dofs (stokes_fe);
+    DoFRenumbering::Cuthill_McKee (stokes_dof_handler);
+    DoFRenumbering::component_wise (stokes_dof_handler, stokes_block_component);
+    
+    stokes_constraints.clear ();
+    DoFTools::make_hanging_node_constraints (stokes_dof_handler,
+                                            stokes_constraints);
+    std::set<unsigned char> no_normal_flux_boundaries;
+    no_normal_flux_boundaries.insert (0);
+    VectorTools::compute_no_normal_flux_constraints (stokes_dof_handler, 0,
+                                                    no_normal_flux_boundaries,
+                                                    stokes_constraints);
+    stokes_constraints.close ();
+  }
+  {
+    temperature_dof_handler.distribute_dofs (temperature_fe);
+    DoFRenumbering::Cuthill_McKee (temperature_dof_handler);
+
+    temperature_constraints.clear ();
+    DoFTools::make_hanging_node_constraints (temperature_dof_handler,
+                                            temperature_constraints);
+    temperature_constraints.close ();
+  }
+  
+  std::vector<unsigned int> stokes_dofs_per_block (2);
+  DoFTools::count_dofs_per_block (stokes_dof_handler, stokes_dofs_per_block,
+                                 stokes_block_component);
+  
+  const unsigned int n_u = stokes_dofs_per_block[0],
+                     n_p = stokes_dofs_per_block[1],
+                    n_T = temperature_dof_handler.n_dofs();
 
   std::cout << "Number of active cells: "
             << triangulation.n_active_cells()
@@ -894,20 +909,11 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
            << " levels)"
             << std::endl
             << "Number of degrees of freedom: "
-            << dof_handler.n_dofs()
+            << n_u + n_p + n_T
             << " (" << n_u << '+' << n_p << '+'<< n_T <<')'
             << std::endl
             << std::endl;
 
-  hanging_node_constraints.clear ();
-  DoFTools::make_hanging_node_constraints (dof_handler,
-                                          hanging_node_constraints);
-  std::set<unsigned char> no_normal_flux_boundaries;
-  no_normal_flux_boundaries.insert (0);
-  VectorTools::compute_no_normal_flux_constraints (dof_handler, 0,
-                                                  no_normal_flux_boundaries,
-                                                  hanging_node_constraints);
-  hanging_node_constraints.close ();
 
   
                                 // The next step is to 
@@ -960,7 +966,7 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
                                 // Observe that we use a 
                                 // coupling argument for 
                                 // telling the function
-                                // <code>make_sparsity_pattern</code>
+                                // <code>make_stokes_sparsity_pattern</code>
                                 // which components actually
                                 // will hold data and which 
                                 // we're going to neglect.
@@ -970,23 +976,18 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
                                 // system matrix structure to
                                 // the sparsity pattern.
   {
-    system_matrix.clear ();
+    stokes_matrix.clear ();
 
-    BlockCompressedSetSparsityPattern csp (3,3);
+    BlockCompressedSetSparsityPattern csp (2,2);
  
     csp.block(0,0).reinit (n_u, n_u);
     csp.block(0,1).reinit (n_u, n_p);
-    csp.block(0,2).reinit (n_u, n_T);
     csp.block(1,0).reinit (n_p, n_u);
     csp.block(1,1).reinit (n_p, n_p);
-    csp.block(1,2).reinit (n_p, n_T);
-    csp.block(2,0).reinit (n_T, n_u);
-    csp.block(2,1).reinit (n_T, n_p);
-    csp.block(2,2).reinit (n_T, n_T);
       
     csp.collect_sizes ();
 
-    Table<2,DoFTools::Coupling> coupling (dim+2, dim+2);
+    Table<2,DoFTools::Coupling> coupling (dim+1, dim+1);
 
                                     // build the sparsity pattern. note
                                     // that all dim velocities couple with
@@ -997,15 +998,13 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
       {
        case 2:
        {
-         static const bool coupling_matrix[4][4]
-           = {{ 1, 1,   1,  0 },
-              { 1, 1,   1,  0 },
+         static const bool coupling_matrix[3][3]
+           = {{ 1, 1,   1},
+              { 1, 1,   1},
                 
-              { 1, 1,   0,  0 },
-
-              { 0, 0,   0,  1 }};
-         for (unsigned int c=0; c<dim+2; ++c)
-           for (unsigned int d=0; d<dim+2; ++d)
+              { 1, 1,   0}};
+         for (unsigned int c=0; c<dim+1; ++c)
+           for (unsigned int d=0; d<dim+1; ++d)
              if (coupling_matrix[c][d] == true)
                coupling[c][d] = DoFTools::always;
              else
@@ -1016,16 +1015,14 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
 
        case 3:
        {
-         static const bool coupling_matrix[5][5]
-           = {{ 1, 1, 1,   1,  0 },
-              { 1, 1, 1,   1,  0 },
-              { 1, 1, 1,   1,  0 },
+         static const bool coupling_matrix[4][4]
+           = {{ 1, 1, 1,   1},
+              { 1, 1, 1,   1},
+              { 1, 1, 1,   1},
                 
-              { 1, 1, 1,   0,  0 },
-
-              { 0, 0, 0,   0,  1 }};
-         for (unsigned int c=0; c<dim+2; ++c)
-           for (unsigned int d=0; d<dim+2; ++d)
+              { 1, 1, 1,   0}};
+         for (unsigned int c=0; c<dim+1; ++c)
+           for (unsigned int d=0; d<dim+1; ++d)
              if (coupling_matrix[c][d] == true)
                coupling[c][d] = DoFTools::always;
              else
@@ -1038,33 +1035,28 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
              Assert (false, ExcNotImplemented());
       } 
       
-    DoFTools::make_sparsity_pattern (dof_handler, coupling, csp);
-    hanging_node_constraints.condense (csp);
-    sparsity_pattern.copy_from (csp);
+    DoFTools::make_sparsity_pattern (stokes_dof_handler, coupling, csp);
+    stokes_constraints.condense (csp);
+    stokes_sparsity_pattern.copy_from (csp);
 
-    system_matrix.reinit (sparsity_pattern);
+    stokes_matrix.reinit (stokes_sparsity_pattern);
   }
-
+  
   {
     Amg_preconditioner.reset ();
     Mp_preconditioner.reset ();
-    preconditioner_matrix.clear ();
+    stokes_preconditioner_matrix.clear ();
 
-    BlockCompressedSetSparsityPattern csp (3,3);
+    BlockCompressedSetSparsityPattern csp (2,2);
  
     csp.block(0,0).reinit (n_u, n_u);
     csp.block(0,1).reinit (n_u, n_p);
-    csp.block(0,2).reinit (n_u, n_T);
     csp.block(1,0).reinit (n_p, n_u);
     csp.block(1,1).reinit (n_p, n_p);
-    csp.block(1,2).reinit (n_p, n_T);
-    csp.block(2,0).reinit (n_T, n_u);
-    csp.block(2,1).reinit (n_T, n_p);
-    csp.block(2,2).reinit (n_T, n_T);
       
     csp.collect_sizes ();
 
-    Table<2,DoFTools::Coupling> coupling (dim+2, dim+2);
+    Table<2,DoFTools::Coupling> coupling (dim+1, dim+1);
 
                                     // build the sparsity pattern. note
                                     // that all dim velocities couple with
@@ -1075,15 +1067,13 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
       {
        case 2:
        {
-         static const bool coupling_matrix[4][4]
-           = {{ 1, 0,   0,  0 },
-              { 0, 1,   0,  0 },
+         static const bool coupling_matrix[3][3]
+           = {{ 1, 0,   0},
+              { 0, 1,   0},
                 
-              { 0, 0,   1,  0 },
-
-              { 0, 0,   0,  0 }};
-         for (unsigned int c=0; c<dim+2; ++c)
-           for (unsigned int d=0; d<dim+2; ++d)
+              { 0, 0,   1}};
+         for (unsigned int c=0; c<dim+1; ++c)
+           for (unsigned int d=0; d<dim+1; ++d)
              if (coupling_matrix[c][d] == true)
                coupling[c][d] = DoFTools::always;
              else
@@ -1094,16 +1084,14 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
 
        case 3:
        {
-         static const bool coupling_matrix[5][5]
-           = {{ 1, 0, 0,   0,  0 },
-              { 0, 1, 0,   0,  0 },
-              { 0, 0, 1,   0,  0 },
+         static const bool coupling_matrix[4][4]
+           = {{ 1, 0, 0,   0},
+              { 0, 1, 0,   0},
+              { 0, 0, 1,   0},
                 
-              { 0, 0, 0,   1,  0 },
-
-              { 0, 0, 0,   0,  0 }};
-         for (unsigned int c=0; c<dim+2; ++c)
-           for (unsigned int d=0; d<dim+2; ++d)
+              { 0, 0, 0,   1}};
+         for (unsigned int c=0; c<dim+1; ++c)
+           for (unsigned int d=0; d<dim+1; ++d)
              if (coupling_matrix[c][d] == true)
                coupling[c][d] = DoFTools::always;
              else
@@ -1114,47 +1102,53 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
 
        default:
              Assert (false, ExcNotImplemented());
-      } 
-      
-    DoFTools::make_sparsity_pattern (dof_handler, coupling, csp);
-    hanging_node_constraints.condense (csp);
-    preconditioner_sparsity_pattern.copy_from (csp);
+      }
+
 
-    preconditioner_matrix.reinit (preconditioner_sparsity_pattern);
+    DoFTools::make_sparsity_pattern (stokes_dof_handler, coupling, csp);
+    stokes_constraints.condense (csp);
+    stokes_preconditioner_sparsity_pattern.copy_from (csp);
+
+    stokes_preconditioner_matrix.reinit (stokes_preconditioner_sparsity_pattern);
   }
 
-                                // As last action in this function,
-                                // we need to set the vectors
-                                // for the solution, the old 
-                                // solution (required for 
-                                // time stepping) and the system
-                                // right hand side to the 
-                                // three-block structure given
-                                // by velocity, pressure and
-                                // temperature.
-  solution.reinit (3);
-  solution.block(0).reinit (n_u);
-  solution.block(1).reinit (n_p);
-  solution.block(2).reinit (n_T);
-  solution.collect_sizes ();
-
-  old_solution.reinit (3);
-  old_solution.block(0).reinit (n_u);
-  old_solution.block(1).reinit (n_p);
-  old_solution.block(2).reinit (n_T);
-  old_solution.collect_sizes ();
-
-  old_old_solution.reinit (3);
-  old_old_solution.block(0).reinit (n_u);
-  old_old_solution.block(1).reinit (n_p);
-  old_old_solution.block(2).reinit (n_T);
-  old_old_solution.collect_sizes ();
+  {
+    temperature_matrix.clear ();
+
+    CompressedSetSparsityPattern csp (n_T, n_T);      
+    DoFTools::make_sparsity_pattern (temperature_dof_handler, csp);
+    temperature_constraints.condense (csp);
+    temperature_sparsity_pattern.copy_from (csp);
+
+    temperature_matrix.reinit (temperature_sparsity_pattern);
+  }
+    
+      
   
-  system_rhs.reinit (3);
-  system_rhs.block(0).reinit (n_u);
-  system_rhs.block(1).reinit (n_p);
-  system_rhs.block(2).reinit (n_T);
-  system_rhs.collect_sizes ();
+                                  // As last action in this function,
+                                  // we need to set the vectors
+                                  // for the solution, the old 
+                                  // solution (required for 
+                                  // time stepping) and the system
+                                  // right hand side to the 
+                                  // three-block structure given
+                                  // by velocity, pressure and
+                                  // temperature.
+  stokes_solution.reinit (2);
+  stokes_solution.block(0).reinit (n_u);
+  stokes_solution.block(1).reinit (n_p);
+  stokes_solution.collect_sizes ();
+
+  stokes_rhs.reinit (2);
+  stokes_rhs.block(0).reinit (n_u);
+  stokes_rhs.block(1).reinit (n_p);
+  stokes_rhs.collect_sizes ();
+  
+  temperature_solution.reinit (n_T);
+  old_temperature_solution.reinit (n_T);
+  old_old_temperature_solution.reinit (n_T);
+
+  temperature_rhs.reinit (n_T);
 }
 
 
@@ -1173,16 +1167,16 @@ double scalar_product (const Tensor<2,dim> &t1,
 
 template <int dim>
 void
-BoussinesqFlowProblem<dim>::assemble_preconditioner ()
+BoussinesqFlowProblem<dim>::assemble_stokes_preconditioner ()
 {
-  preconditioner_matrix = 0;
+  stokes_preconditioner_matrix = 0;
 
   QGauss<dim>   quadrature_formula(degree+2);
-  FEValues<dim> fe_values (fe, quadrature_formula,
-                          update_JxW_values |
-                          update_values |
-                          update_gradients);
-  const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+  FEValues<dim> stokes_fe_values (stokes_fe, quadrature_formula,
+                                 update_JxW_values |
+                                 update_values |
+                                 update_gradients);
+  const unsigned int   dofs_per_cell   = stokes_fe.dofs_per_cell;
 
   const unsigned int   n_q_points      = quadrature_formula.size();
 
@@ -1196,19 +1190,19 @@ BoussinesqFlowProblem<dim>::assemble_preconditioner ()
   const FEValuesExtractors::Scalar pressure (dim);
 
   typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
+    cell = stokes_dof_handler.begin_active(),
+    endc = stokes_dof_handler.end();
   for (; cell!=endc; ++cell)
     {
-      fe_values.reinit (cell);
+      stokes_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)
            {
-             phi_grad_u[k] = fe_values[velocities].gradient(k,q);
-             phi_p[k]      = fe_values[pressure].value (k, q);
+             phi_grad_u[k] = stokes_fe_values[velocities].gradient(k,q);
+             phi_p[k]      = stokes_fe_values[pressure].value (k, q);
            }
          
          for (unsigned int i=0; i<dofs_per_cell; ++i)
@@ -1216,24 +1210,119 @@ BoussinesqFlowProblem<dim>::assemble_preconditioner ()
              local_matrix(i,j) += (scalar_product (phi_grad_u[i], phi_grad_u[j])
                                    +
                                    phi_p[i] * phi_p[j])
-                                  * fe_values.JxW(q);
+                                  * stokes_fe_values.JxW(q);
        }
 
       cell->get_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)
-         preconditioner_matrix.add (local_dof_indices[i],
-                                    local_dof_indices[j],
-                                    local_matrix(i,j));
+         stokes_preconditioner_matrix.add (local_dof_indices[i],
+                                           local_dof_indices[j],
+                                           local_matrix(i,j));
     }
   
-  hanging_node_constraints.condense (preconditioner_matrix);
+  stokes_constraints.condense (stokes_preconditioner_matrix);
+}
+
+
+
+template <int dim>
+void
+BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
+{
+  std::cout << "   Rebuilding Stokes preconditioner..." << std::flush;
+      
+
+                                  // This last step of the assembly
+                                  // function sets up the preconditioners
+                                  // used for the solution of the
+                                  // system. We are going to use an
+                                  // ILU preconditioner for the
+                                  // velocity block (to be used
+                                  // by BlockSchurPreconditioner class)
+                                  // as well as an ILU preconditioner
+                                  // for the inversion of the 
+                                  // pressure mass matrix. Recall that
+                                  // the velocity-velocity block sits
+                                  // at position (0,0) in the 
+                                  // global system matrix, and
+                                  // the pressure mass matrix in
+                                  // (1,1). The 
+                                  // storage of these objects is
+                                  // as in step-22, that is, we
+                                  // include them using a 
+                                  // shared pointer structure from the
+                                  // boost library.
+                                  // 
+                                  // When all work is done, we 
+                                  // change the flags 
+                                  // <code>rebuild_preconditioner</code>
+                                  // and 
+                                  // <code>rebuild_matrices</code>
+                                  // to false.
+  assemble_stokes_preconditioner ();
+      
+  Amg_preconditioner = boost::shared_ptr<PreconditionerTrilinosAmg>
+                      (new PreconditionerTrilinosAmg());
+      
+  const unsigned int n_u = stokes_preconditioner_matrix.block(0,0).m();
+  std::vector<double> null_space (dim * n_u, 0.);
+      
+  std::vector<bool> precondition_dof_list (stokes_dof_handler.n_dofs(), false);
+      
+  for (unsigned int component=0; component < dim; ++component)
+    {
+      std::vector<bool> precondition_mask (dim + 1, false);
+      precondition_mask[component] = true;
+      DoFTools::extract_dofs (stokes_dof_handler, precondition_mask, 
+                             precondition_dof_list);
+
+                                      // TODO: The current implementation 
+                                      // assumes that we are working on 
+                                      // the first components of a system when
+                                      // writing into the null vector.
+                                      // Change this to the general case,
+                                      // probably use something similar as
+                                      // for block vectors.
+      unsigned int counter = 0;
+      for (unsigned int i=0; i<stokes_dof_handler.n_dofs(); ++i)
+       {
+         if (precondition_dof_list[i])
+           {
+             Assert(i < n_u,
+                    ExcMessage("Could not correctly locate "
+                               "preconditioner dofs in system!"));
+             null_space [component * n_u + i] = 1.;
+             ++counter;
+           }
+       }
+      Assert (counter == n_u / dim,
+             ExcDimensionMismatch(counter, n_u / dim));
+    }
+       
+  Amg_preconditioner->initialize(stokes_preconditioner_matrix.block(0,0),
+                                null_space, dim, true, true, false);
+
+                                  // TODO: we could throw away the (0,0)
+                                  // block here since things have been
+                                  // copied over to Trilinos. we need to
+                                  // keep the (1,1) block, though
+      
+  Mp_preconditioner
+    = boost::shared_ptr<SparseILU<double> >
+    (new SparseILU<double>);
+  Mp_preconditioner->initialize (stokes_preconditioner_matrix.block(1,1),
+                                SparseILU<double>::AdditionalData());
+      
+  std::cout << std::endl;
+
+  rebuild_preconditioner = false;
 }
 
 
 
-                                // @sect4{BoussinesqFlowProblem::assemble_system}
+                                // @sect4{BoussinesqFlowProblem::assemble_stokes_system}
                                 // 
                                 // The assembly of the Boussinesq 
                                 // system is acutally a two-step
@@ -1302,32 +1391,38 @@ BoussinesqFlowProblem<dim>::assemble_preconditioner ()
                                 // the local dofs compared to
                                 // the global system.
 template <int dim>
-void BoussinesqFlowProblem<dim>::assemble_system ()
+void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
 {
+  std::cout << "   Assembling..." << std::flush;
+
   if (rebuild_matrices == true)
-    system_matrix=0;
+    stokes_matrix=0;
 
-  system_rhs=0;
+  stokes_rhs=0;
 
   QGauss<dim>   quadrature_formula(degree+2);
   QGauss<dim-1> face_quadrature_formula(degree+2);
 
-  FEValues<dim> fe_values (fe, quadrature_formula,
-                          update_values    |
-                          update_quadrature_points  |
-                          update_JxW_values |
-                          (rebuild_matrices == true
-                           ?
-                           update_gradients
-                           :
-                           UpdateFlags(0)));
-  FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
-                                   update_values    | 
-                                   update_normal_vectors |
-                                   update_quadrature_points  | 
-                                   update_JxW_values);
-
-  const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+  FEValues<dim> stokes_fe_values (stokes_fe, quadrature_formula,
+                                 update_values    |
+                                 update_quadrature_points  |
+                                 update_JxW_values |
+                                 (rebuild_matrices == true
+                                  ?
+                                  update_gradients
+                                  :
+                                  UpdateFlags(0)));
+
+  FEValues<dim> temperature_fe_values (temperature_fe, quadrature_formula,
+                                      update_values);
+
+  FEFaceValues<dim> stokes_fe_face_values (stokes_fe, face_quadrature_formula,
+                                          update_values    | 
+                                          update_normal_vectors |
+                                          update_quadrature_points  | 
+                                          update_JxW_values);
+
+  const unsigned int   dofs_per_cell   = stokes_fe.dofs_per_cell;
 
   const unsigned int   n_q_points      = quadrature_formula.size();
   const unsigned int   n_face_q_points = face_quadrature_formula.size();
@@ -1337,43 +1432,42 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
 
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
-                                // These few declarations provide
-                                // the structures for the evaluation
-                                // of inhomogeneous Neumann boundary
-                                // conditions from the function
-                                // declaration made above.
-                                // The vector <code>old_solution_values</code>
-                                // evaluates the solution 
-                                // at the old time level, since
-                                // the temperature from the
-                                // old time level enters the 
-                                // Stokes system as a source
-                                // term in the momentum equation.
-                                // 
-                                // Then, we create a variable
-                                // to hold the Rayleigh number,
-                                // the measure of buoyancy.
-                                // 
-                                // The set of vectors we create
-                                // next hold the evaluations of
-                                // the basis functions that will
-                                // be used for creating the
-                                // matrices. This gives faster
-                                // access to that data, which
-                                // increases the performance
-                                // of the assembly. See step-22 
-                                // for details.
-                                // 
-                                // The last few declarations 
-                                // are used to extract the 
-                                // individual blocks (velocity,
-                                // pressure, temperature) from
-                                // the total FE system.
+                                  // These few declarations provide
+                                  // the structures for the evaluation
+                                  // of inhomogeneous Neumann boundary
+                                  // conditions from the function
+                                  // declaration made above.
+                                  // The vector <code>old_solution_values</code>
+                                  // evaluates the solution 
+                                  // at the old time level, since
+                                  // the temperature from the
+                                  // old time level enters the 
+                                  // Stokes system as a source
+                                  // term in the momentum equation.
+                                  // 
+                                  // Then, we create a variable
+                                  // to hold the Rayleigh number,
+                                  // the measure of buoyancy.
+                                  // 
+                                  // The set of vectors we create
+                                  // next hold the evaluations of
+                                  // the basis functions that will
+                                  // be used for creating the
+                                  // matrices. This gives faster
+                                  // access to that data, which
+                                  // increases the performance
+                                  // of the assembly. See step-22 
+                                  // for details.
+                                  // 
+                                  // The last few declarations 
+                                  // are used to extract the 
+                                  // individual blocks (velocity,
+                                  // pressure, temperature) from
+                                  // the total FE system.
   const PressureBoundaryValues<dim> pressure_boundary_values;
   std::vector<double>               boundary_values (n_face_q_points);
 
-  std::vector<Vector<double> >      old_solution_values(n_q_points,
-                                                       Vector<double>(dim+2));
+  std::vector<double>               old_temperature_values(n_q_points);
 
   const double Rayleigh_number = 10;
 
@@ -1384,178 +1478,167 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
 
   const FEValuesExtractors::Vector velocities (0);
   const FEValuesExtractors::Scalar pressure (dim);
-  const FEValuesExtractors::Scalar temperature (dim+1);
-
-                                // Now starts the loop over
-                                // all cells in the problem.
-                                // The first commands are all
-                                // very familiar, doing the
-                                // evaluations of the element
-                                // basis functions, resetting
-                                // the local arrays and 
-                                // getting the values of the
-                                // old solution at the
-                                // quadrature point. Then we
-                                // are ready to loop over
-                                // the quadrature points 
-                                // on the cell.
+
+                                  // Now start the loop over
+                                  // all cells in the problem.
+                                  // The first commands are all
+                                  // very familiar, doing the
+                                  // evaluations of the element
+                                  // basis functions, resetting
+                                  // the local arrays and 
+                                  // getting the values of the
+                                  // old solution at the
+                                  // quadrature point. Then we
+                                  // are ready to loop over
+                                  // the quadrature points 
+                                  // on the cell.
   typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-  for (; cell!=endc; ++cell)
+    cell = stokes_dof_handler.begin_active(),
+    endc = stokes_dof_handler.end();
+  typename DoFHandler<dim>::active_cell_iterator
+    temperature_cell = temperature_dof_handler.begin_active();
+  
+  for (; cell!=endc; ++cell, ++temperature_cell)
     {
-      fe_values.reinit (cell);
+      stokes_fe_values.reinit (cell);
+      temperature_fe_values.reinit (temperature_cell);
+      
       local_matrix = 0;
       local_rhs = 0;
 
-      fe_values.get_function_values (old_solution, old_solution_values);
+      temperature_fe_values.get_function_values (old_temperature_solution, old_temperature_values);
 
       for (unsigned int q=0; q<n_q_points; ++q)
        {
-         const double old_temperature = old_solution_values[q](dim+1);
-
-                               // Extract the basis relevant
-                               // terms in the inner products
-                               // once in advance as shown
-                               // in step-22 in order to 
-                               // accelerate assembly.
-                               // 
-                               // Once this is done, we 
-                               // start the loop over the
-                               // rows and columns of the
-                               // local matrix and feed
-                               // the matrix with the relevant
-                               // products. The right hand
-                               // side is filled with the 
-                               // forcing term driven by
-                               // temperature in direction
-                               // of gravity (which is 
-                               // vertical in our example).
-                               // Note that the right hand 
-                               // side term is always generated,
-                               // whereas the matrix 
-                               // contributions are only
-                               // updated when it is 
-                               // requested by the
-                               // <code>rebuild_matrices</code>
-                               // flag.
+         const double old_temperature = old_temperature_values[q];
+
+                                          // Extract the basis relevant
+                                          // terms in the inner products
+                                          // once in advance as shown
+                                          // in step-22 in order to 
+                                          // accelerate assembly.
+                                          // 
+                                          // Once this is done, we 
+                                          // start the loop over the
+                                          // rows and columns of the
+                                          // local matrix and feed
+                                          // the matrix with the relevant
+                                          // products. The right hand
+                                          // side is filled with the 
+                                          // forcing term driven by
+                                          // temperature in direction
+                                          // of gravity (which is 
+                                          // vertical in our example).
+                                          // Note that the right hand 
+                                          // side term is always generated,
+                                          // whereas the matrix 
+                                          // contributions are only
+                                          // updated when it is 
+                                          // requested by the
+                                          // <code>rebuild_matrices</code>
+                                          // flag.
          for (unsigned int k=0; k<dofs_per_cell; ++k)
            {
-             phi_u[k] = fe_values[velocities].value (k,q);
+             phi_u[k] = stokes_fe_values[velocities].value (k,q);
              if (rebuild_matrices)
                {
-                 grads_phi_u[k] = fe_values[velocities].symmetric_gradient(k,q);
-                 div_phi_u[k]   = fe_values[velocities].divergence (k, q);
-                 phi_p[k]       = fe_values[pressure].value (k, q);
+                 grads_phi_u[k] = stokes_fe_values[velocities].symmetric_gradient(k,q);
+                 div_phi_u[k]   = stokes_fe_values[velocities].divergence (k, q);
+                 phi_p[k]       = stokes_fe_values[pressure].value (k, q);
                }
            }
 
+                                          // define viscosity
+         const double eta = 1;
+         if (rebuild_matrices)
+           for (unsigned int i=0; i<dofs_per_cell; ++i)
+             for (unsigned int j=0; j<dofs_per_cell; ++j)
+               local_matrix(i,j) += (eta * grads_phi_u[i] * grads_phi_u[j]
+                                     - div_phi_u[i] * phi_p[j]
+                                     - phi_p[i] * div_phi_u[j])
+                                    * stokes_fe_values.JxW(q);
+
+         const Point<dim> gravity = ( (dim == 2) ? (Point<dim> (0,1)) : 
+                                      (Point<dim> (0,0,1)) );
          for (unsigned int i=0; i<dofs_per_cell; ++i)
-           {
-             const Tensor<1,dim> phi_i_u = fe_values[velocities].value (i, q);
-
-                                              // define viscosity
-             const double eta = 1;
-
-                                              // build Stokes part of
-                                              // the matrix. we have
-                                              // to rebuild the
-                                              // temperature part of
-                                              // it in each time step
-                                              // due to the fact that
-                                              // it depends on the
-                                              // time step size
-             if (rebuild_matrices)
-               for (unsigned int j=0; j<dofs_per_cell; ++j)
-                 local_matrix(i,j) += (eta * grads_phi_u[i] * grads_phi_u[j]
-                                       - div_phi_u[i] * phi_p[j]
-                                       - phi_p[i] * div_phi_u[j])
-                                      * fe_values.JxW(q);
-
-             const Point<dim> gravity = ( (dim == 2) ? (Point<dim> (0,1)) : 
-                                                       (Point<dim> (0,0,1)) );
-
-             local_rhs(i) += (Rayleigh_number *
-                              gravity * phi_u[i] * old_temperature)*
-                             fe_values.JxW(q);
-          }
+           local_rhs(i) += (Rayleigh_number *
+                            gravity * phi_u[i] * old_temperature)*
+                           stokes_fe_values.JxW(q);
        }
 
-                                // Next follows the assembly 
-                                // of the face terms, result
-                                // from Neumann boundary 
-                                // conditions. Since these
-                                // terms only enter the right
-                                // hand side vector and not
-                                // the matrix, there is no
-                                // substantial benefit from
-                                // extracting the data 
-                                // before using it, so 
-                                // we remain in the lines 
-                                // of step-20 at this point.
+
+                                      // Next follows the assembly 
+                                      // of the face terms, result
+                                      // from Neumann boundary 
+                                      // conditions. Since these
+                                      // terms only enter the right
+                                      // hand side vector and not
+                                      // the matrix, there is no
+                                      // substantial benefit from
+                                      // extracting the data 
+                                      // before using it, so 
+                                      // we remain in the lines 
+                                      // of step-20 at this point.
       for (unsigned int face_no=0;
-           face_no<GeometryInfo<dim>::faces_per_cell;
-           ++face_no)
-        if (cell->at_boundary(face_no))
-          {
-            fe_face_values.reinit (cell, face_no);
-
-            pressure_boundary_values
-              .value_list (fe_face_values.get_quadrature_points(),
-                           boundary_values);
-
-            for (unsigned int q=0; q<n_face_q_points; ++q)
-              for (unsigned int i=0; i<dofs_per_cell; ++i)
-                {
-                  const Tensor<1,dim>
-                    phi_i_u = fe_face_values[velocities].value (i, q);
-
-                  local_rhs(i) += -(phi_i_u *
-                                    fe_face_values.normal_vector(q) *
-                                    boundary_values[q] *
-                                    fe_face_values.JxW(q));
-                }
-          }
-
-                                // The last step in the loop 
-                                // over all cells is to
-                                // enter the local contributions
-                                // into the global matrix and 
-                                // vector structures to the
-                                // positions specified in 
-                                // <code>local_dof_indices</code>.
-                                // Again, we only add the 
-                                // matrix data when it is 
-                                // requested.
+          face_no<GeometryInfo<dim>::faces_per_cell;
+          ++face_no)
+       if (cell->at_boundary(face_no))
+         {
+           stokes_fe_face_values.reinit (cell, face_no);
+
+           pressure_boundary_values
+             .value_list (stokes_fe_face_values.get_quadrature_points(),
+                          boundary_values);
+
+           for (unsigned int q=0; q<n_face_q_points; ++q)
+             for (unsigned int i=0; i<dofs_per_cell; ++i)
+               {
+                 const Tensor<1,dim>
+                   phi_i_u = stokes_fe_face_values[velocities].value (i, q);
+
+                 local_rhs(i) += -(phi_i_u *
+                                   stokes_fe_face_values.normal_vector(q) *
+                                   boundary_values[q] *
+                                   stokes_fe_face_values.JxW(q));
+               }
+         }      
+
+                                      // The last step in the loop 
+                                      // over all cells is to
+                                      // enter the local contributions
+                                      // into the global matrix and 
+                                      // vector structures to the
+                                      // positions specified in 
+                                      // <code>local_dof_indices</code>.
+                                      // Again, we only add the 
+                                      // matrix data when it is 
+                                      // requested.
       cell->get_dof_indices (local_dof_indices);
 
       if (rebuild_matrices == true)
-       {
-         for (unsigned int i=0; i<dofs_per_cell; ++i)
-           for (unsigned int j=0; j<dofs_per_cell; ++j)
-             system_matrix.add (local_dof_indices[i],
-                                local_dof_indices[j],
-                                local_matrix(i,j));
-       }
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         for (unsigned int j=0; j<dofs_per_cell; ++j)
+           stokes_matrix.add (local_dof_indices[i],
+                              local_dof_indices[j],
+                              local_matrix(i,j));
 
       for (unsigned int i=0; i<dofs_per_cell; ++i)
-        system_rhs(local_dof_indices[i]) += local_rhs(i);
+        stokes_rhs(local_dof_indices[i]) += local_rhs(i);
     }
 
-                                // Back at the outermost
-                                // level of this function,
-                                // we continue the work
-                                // by condensing hanging
-                                // node constraints to the
-                                // right hand side and, 
-                                // possibly, to the matrix.
-  if (rebuild_matrices == true)
-    hanging_node_constraints.condense (system_matrix);
-
-  hanging_node_constraints.condense (system_rhs);
+                                  // Back at the outermost
+                                  // level of this function,
+                                  // we continue the work
+                                  // by condensing hanging
+                                  // node constraints to the
+                                  // right hand side and, 
+                                  // possibly, to the matrix.
+  stokes_constraints.condense (stokes_rhs);
 
   if (rebuild_matrices == true)
     {
+      stokes_constraints.condense (stokes_matrix);
+
 //       std::map<unsigned int,double> boundary_values;
 
 //       typename DoFHandler<dim>::active_cell_iterator
@@ -1585,107 +1668,14 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
 //                                             component_mask);
 
 //       MatrixTools::apply_boundary_values (boundary_values,
-//                                       system_matrix,
+//                                       stokes_matrix,
 //                                       solution,
 //                                       system_rhs);
-    }
-
-                                // This last step of the assembly
-                                // function sets up the preconditioners
-                                // used for the solution of the
-                                // system. We are going to use an
-                                // ILU preconditioner for the
-                                // velocity block (to be used
-                                // by BlockSchurPreconditioner class)
-                                // as well as an ILU preconditioner
-                                // for the inversion of the 
-                                // pressure mass matrix. Recall that
-                                // the velocity-velocity block sits
-                                // at position (0,0) in the 
-                                // global system matrix, and
-                                // the pressure mass matrix in
-                                // (1,1). The 
-                                // storage of these objects is
-                                // as in step-22, that is, we
-                                // include them using a 
-                                // shared pointer structure from the
-                                // boost library.
-                                // 
-                                // When all work is done, we 
-                                // change the flags 
-                                // <code>rebuild_preconditioner</code>
-                                // and 
-                                // <code>rebuild_matrices</code>
-                                // to false.
-  if (rebuild_preconditioner == true)
-    {
-      Assert (rebuild_matrices == true,
-             ExcMessage ("There is no point in rebuilding the preconditioner "
-                         "without a rebuilt matrix!"));
-
-      std::cout << "   Rebuilding preconditioner..." << std::flush;
-      
-      assemble_preconditioner ();
-      
-      Amg_preconditioner = boost::shared_ptr<PreconditionerTrilinosAmg>
-                                       (new PreconditionerTrilinosAmg());
-      
-      const unsigned int n_u = preconditioner_matrix.block(0,0).m();
-      std::vector<double> null_space (dim * n_u, 0.);
-      
-      std::vector<bool> precondition_dof_list (dof_handler.n_dofs(), false);
-      std::vector<bool> precondition_mask (dim + 2, false);
-      
-      for (unsigned int component=0; component < dim; component++)
-       {
-         precondition_mask[component] = true;
-         DoFTools::extract_dofs (dof_handler, precondition_mask, 
-                                 precondition_dof_list);
-         precondition_mask[component] = false;
-
-                               // TODO: The current implementation 
-                               // assumes that we are working on 
-                               // the first components of a system when
-                               // writing into the null vector.
-                               // Change this to the general case,
-                               // probably use something similar as
-                               // for block vectors.
-         unsigned int counter = 0;
-         for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
-           {
-             if (precondition_dof_list[i])
-               {
-                 Assert(i < n_u,
-                        ExcMessage("Could not correctly locate "
-                                   "preconditioner dofs in system!"));
-                 null_space [component * n_u + i] = 1.;
-                 ++counter;
-               }
-           }
-         Assert (counter == n_u / dim,
-                 ExcDimensionMismatch(counter, n_u / dim));
-       }
-       
-      Amg_preconditioner->initialize(preconditioner_matrix.block(0,0),
-                                    null_space, dim, true, true, false);
-
-                                      // TODO: we could throw away the (0,0)
-                                      // block here since things have been
-                                      // copied over to Trilinos. we need to
-                                      // keep the (1,1) block, though
-      
-      Mp_preconditioner
-       = boost::shared_ptr<SparseILU<double> >
-               (new SparseILU<double>);
-      Mp_preconditioner->initialize (preconditioner_matrix.block(1,1),
-                                    SparseILU<double>::AdditionalData());
-      
-      std::cout << std::endl;
-
-      rebuild_preconditioner = false;
     }
 
   rebuild_matrices = false;
+
+  std::cout << std::endl;
 }
 
 
@@ -1693,21 +1683,21 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
 
 template <int dim>
 double compute_viscosity(
-       const std::vector<Vector<double> >              present_solution,
-       const std::vector<Vector<double> >              old_solution,
-       const std::vector<Vector<double> >              old_old_solution,
-       const std::vector<std::vector<Tensor<1,dim> > > old_solution_grads,
-       const std::vector<std::vector<Tensor<1,dim> > > old_old_solution_grads,
-       const std::vector<std::vector<Tensor<2,dim> > > old_solution_hessians,
-       const std::vector<std::vector<Tensor<2,dim> > > old_old_solution_hessians,
-       const std::vector<double>                       gamma_values,
-       const double                                    kappa,
-       const double                                    global_u_infty,
-       const double                                    global_T_infty,
-       const double                                    global_Omega_diameter,
-       const double                                    cell_diameter,
-       const double                                    old_time_step
-       )
+  const std::vector<double>          &old_temperature,
+  const std::vector<double>          &old_old_temperature,
+  const std::vector<Tensor<1,dim> >  &old_temperature_grads,
+  const std::vector<Tensor<1,dim> >  &old_old_temperature_grads,
+  const std::vector<Tensor<2,dim> >  &old_temperature_hessians,
+  const std::vector<Tensor<2,dim> >  &old_old_temperature_hessians,
+  const std::vector<Vector<double> > &present_stokes_values,
+  const std::vector<double>          &gamma_values,
+  const double                        kappa,
+  const double                        global_u_infty,
+  const double                        global_T_infty,
+  const double                        global_Omega_diameter,
+  const double                        cell_diameter,
+  const double                        old_time_step
+)
 {
   const double beta = 0.1;
   const double alpha = 1;
@@ -1715,43 +1705,39 @@ double compute_viscosity(
   if (global_u_infty == 0)
     return 5e-3 * cell_diameter;
   
-  const unsigned int n_q_points = old_solution.size();
+  const unsigned int n_q_points = old_temperature.size();
   
-                                // Stage 1: calculate residual
+                                  // Stage 1: calculate residual
   double max_residual = 0;
   double max_velocity = 0;
   
   for (unsigned int q=0; q < n_q_points; ++q)
     {
-      const double dT_dt = (old_solution[q](dim+1) - old_old_solution[q](dim+1))
-                           / old_time_step;
-      double u_grad_T = 0.;
+      Tensor<1,dim> u;
       for (unsigned int d=0; d<dim; ++d)
-       u_grad_T += present_solution[q](d)*((old_solution_grads[q][dim+1][d] +
-                                            old_old_solution_grads[q][dim+1][d]) / 2);
+       u[d] = present_stokes_values[q](d);
+      
+      const double dT_dt = (old_temperature[q] - old_old_temperature[q])
+                          / old_time_step;
+      const double u_grad_T = u * (old_temperature_grads[q] +
+                                  old_old_temperature_grads[q]) / 2;
       
       const double kappa_Delta_T = kappa
-                                  * (trace(old_solution_hessians[q][dim+1]) +
-                                     trace(old_old_solution_hessians[q][dim+1])) / 2;
+                                  * (trace(old_temperature_hessians[q]) +
+                                     trace(old_old_temperature_hessians[q])) / 2;
 
       const double residual
        = std::abs((dT_dt + u_grad_T - kappa_Delta_T - gamma_values[q]) *
-                  std::pow((old_solution[q](dim+1)+old_old_solution[q](dim+1)) / 2,
+                  std::pow((old_temperature[q]+old_old_temperature[q]) / 2,
                            alpha-1.));
 
-      max_residual = std::max (residual, max_residual);
-
-      double velocity_squared = 0;
-      for (unsigned int d=0; d<dim; ++d)
-       velocity_squared += present_solution[q](d) * present_solution[q](d);
-
-      max_velocity = std::max (std::sqrt (velocity_squared),
-                              max_velocity);
+      max_residual = std::max (residual,        max_residual);
+      max_velocity = std::max (std::sqrt (u*u), max_velocity);
     }
   
   const double global_scaling = global_u_infty * global_T_infty /
-      std::pow(global_Omega_diameter, alpha - 2.);
-  
+                               std::pow(global_Omega_diameter, alpha - 2.);
+
   return (beta *
          max_velocity *
          std::min (cell_diameter,
@@ -1760,7 +1746,7 @@ double compute_viscosity(
 
 
 
-                                // @sect4{BoussinesqFlowProblem::assemble_rhs_T}
+                                // @sect4{BoussinesqFlowProblem::assemble_temperature_system}
                                 // 
                                 // This function does the second
                                 // part of the assembly work, the
@@ -1784,19 +1770,22 @@ double compute_viscosity(
                                 // flags at face level are the 
                                 // same as in step-12.
 template <int dim>
-void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
+void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
 {
   const bool use_bdf2_scheme = (timestep_number != 0);
 
-  system_matrix.block(2,2) = 0;
+  temperature_matrix = 0;
+  temperature_rhs = 0;
   
   QGauss<dim>   quadrature_formula(degree+2);
-  FEValues<dim> fe_values (fe, quadrature_formula,
-                           update_values    | update_gradients |
-                          update_hessians |
-                           update_quadrature_points  | update_JxW_values);
-
-  const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+  FEValues<dim> temperature_fe_values (temperature_fe, quadrature_formula,
+                                      update_values    | update_gradients |
+                                      update_hessians |
+                                      update_quadrature_points  | update_JxW_values);
+  FEValues<dim> stokes_fe_values (stokes_fe, quadrature_formula,
+                                 update_values);
+
+  const unsigned int   dofs_per_cell   = temperature_fe.dofs_per_cell;
   const unsigned int   n_q_points      = quadrature_formula.size();
 
   Vector<double>       local_rhs (dofs_per_cell);
@@ -1804,43 +1793,34 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
 
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
-                                // Here comes the declaration
-                                // of vectors to hold the old
-                                // and present solution values
-                                // and gradients
-                                // for both the cell as well as faces
-                                // to the cell. Next comes the
-                                // declaration of an object
-                                // to hold the temperature 
-                                // boundary values and a
-                                // well-known extractor for
-                                // accessing the temperature
-                                // part of the FE system.
-  std::vector<Vector<double> > present_solution_values (n_q_points, 
-                                                       Vector<double>(dim+2));
-  std::vector<Vector<double> > old_solution_values(n_q_points,
-                                                  Vector<double>(dim+2));
-  std::vector<Vector<double> > old_old_solution_values(n_q_points,
-                                                      Vector<double>(dim+2));
-  std::vector<std::vector<Tensor<1,dim> > >
-    old_solution_grads(n_q_points,
-                      std::vector<Tensor<1,dim> >(dim+2));
-  std::vector<std::vector<Tensor<1,dim> > >
-    old_old_solution_grads(n_q_points,
-                          std::vector<Tensor<1,dim> >(dim+2));
-  std::vector<std::vector<Tensor<2,dim> > > old_solution_hessians(
-                                 n_q_points,
-                                 std::vector<Tensor<2,dim> >(dim+2));
-  std::vector<std::vector<Tensor<2,dim> > > old_old_solution_hessians(
-                                 n_q_points,
-                                 std::vector<Tensor<2,dim> >(dim+2));
+                                  // Here comes the declaration
+                                  // of vectors to hold the old
+                                  // and present solution values
+                                  // and gradients
+                                  // for both the cell as well as faces
+                                  // to the cell. Next comes the
+                                  // declaration of an object
+                                  // to hold the temperature 
+                                  // boundary values and a
+                                  // well-known extractor for
+                                  // accessing the temperature
+                                  // part of the FE system.
+  std::vector<Vector<double> > present_stokes_values (n_q_points, 
+                                                     Vector<double>(dim+1));
+
+  
+  std::vector<double>         old_temperature_values (n_q_points);
+  std::vector<double>         old_old_temperature_values(n_q_points);
+  std::vector<Tensor<1,dim> > old_temperature_grads(n_q_points);
+  std::vector<Tensor<1,dim> > old_old_temperature_grads(n_q_points);
+  std::vector<Tensor<2,dim> > old_temperature_hessians(n_q_points);
+  std::vector<Tensor<2,dim> > old_old_temperature_hessians(n_q_points);
 
+  
   TemperatureBoundaryValues<dim> temperature_boundary_values;
-  RightHandSide<dim> right_hand_side;
+  TemperatureRightHandSide<dim>  temperature_right_hand_side;
   std::vector<double> gamma_values (n_q_points);
 
-  const FEValuesExtractors::Scalar temperature (dim+1);
-
   std::vector<double>                  phi_T       (dofs_per_cell);
   std::vector<Tensor<1,dim> >          grad_phi_T  (dofs_per_cell);
   
@@ -1848,39 +1828,43 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
   const double global_T_infty = get_maximal_temperature();
   const double global_Omega_diameter = GridTools::diameter (triangulation);
 
-                                // Now, let's start the loop
-                                // over all cells in the
-                                // triangulation. The first
-                                // actions within the loop
-                                // are, 0as usual, the evaluation
-                                // of the FE basis functions 
-                                // and the old and present
-                                // solution at the quadrature 
-                                // points.
+                                  // Now, let's start the loop
+                                  // over all cells in the
+                                  // triangulation. The first
+                                  // actions within the loop
+                                  // are, 0as usual, the evaluation
+                                  // of the FE basis functions 
+                                  // and the old and present
+                                  // solution at the quadrature 
+                                  // points.
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = temperature_dof_handler.begin_active(),
+    endc = temperature_dof_handler.end();
   typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
+    stokes_cell = stokes_dof_handler.begin_active();
 
-  for (; cell!=endc; ++cell)
+  for (; cell!=endc; ++cell, ++stokes_cell)
     {
       local_rhs = 0;
       local_matrix = 0;
 
-      fe_values.reinit (cell);
+      temperature_fe_values.reinit (cell);
+      stokes_fe_values.reinit (stokes_cell);
 
-      fe_values.get_function_values (solution, present_solution_values);
-      fe_values.get_function_values (old_solution, old_solution_values);
-      fe_values.get_function_values (old_old_solution, old_old_solution_values);
+      temperature_fe_values.get_function_values (old_temperature_solution, old_temperature_values);
+      temperature_fe_values.get_function_values (old_old_temperature_solution, old_old_temperature_values);
 
-      fe_values.get_function_gradients (old_solution, old_solution_grads);
-      fe_values.get_function_gradients (old_old_solution, old_old_solution_grads);
+      temperature_fe_values.get_function_gradients (old_temperature_solution, old_temperature_grads);
+      temperature_fe_values.get_function_gradients (old_old_temperature_solution, old_old_temperature_grads);
       
-      fe_values.get_function_hessians (old_solution, old_solution_hessians);
-      fe_values.get_function_hessians (old_old_solution, old_old_solution_hessians);
+      temperature_fe_values.get_function_hessians (old_temperature_solution, old_temperature_hessians);
+      temperature_fe_values.get_function_hessians (old_old_temperature_solution, old_old_temperature_hessians);
       
-      right_hand_side.value_list (fe_values.get_quadrature_points(),
-                                 gamma_values, dim+1);
+      temperature_right_hand_side.value_list (temperature_fe_values.get_quadrature_points(),
+                                             gamma_values);
 
+      stokes_fe_values.get_function_values (stokes_solution, present_stokes_values);
+      
                                       // build matrix contributions
 
                                       // define diffusion. take the
@@ -1890,32 +1874,37 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
                                       // impirically) to keep the
                                       // scheme stable
       const double kappa = 1e-6;
-      const double nu = 
-       compute_viscosity (present_solution_values, old_solution_values,
-                          old_old_solution_values, old_solution_grads, old_old_solution_grads,
-                          old_solution_hessians, old_old_solution_hessians, gamma_values,
-                          kappa, global_u_infty, global_T_infty,
-                          global_Omega_diameter, cell->diameter(),
-                          old_time_step);
+      const double nu
+       = compute_viscosity (old_temperature_values,
+                            old_old_temperature_values,
+                            old_temperature_grads,
+                            old_old_temperature_grads,
+                            old_temperature_hessians,
+                            old_old_temperature_hessians,
+                            present_stokes_values,
+                            gamma_values,
+                            kappa, global_u_infty, global_T_infty,
+                            global_Omega_diameter, cell->diameter(),
+                            old_time_step);
       
       for (unsigned int q=0; q<n_q_points; ++q)
        {
          for (unsigned int k=0; k<dofs_per_cell; ++k)
            {
-             grad_phi_T[k] = fe_values[temperature].gradient(k,q);
-             phi_T[k]      = fe_values[temperature].value (k, q);
+             grad_phi_T[k] = temperature_fe_values.shape_grad (k,q);
+             phi_T[k]      = temperature_fe_values.shape_value (k, q);
            }
 
-         const double        old_T      = old_solution_values[q](dim+1);
-         const double        old_old_T  = old_old_solution_values[q](dim+1);
+         const double        old_T      = old_temperature_values[q];
+         const double        old_old_T  = old_old_temperature_values[q];
 
-         const Tensor<1,dim> old_grad_T     = old_solution_grads[q][dim+1];
-         const Tensor<1,dim> old_old_grad_T = old_old_solution_grads[q][dim+1];
+         const Tensor<1,dim> old_grad_T     = old_temperature_grads[q];
+         const Tensor<1,dim> old_old_grad_T = old_old_temperature_grads[q];
 
          
          Tensor<1,dim> present_u;
          for (unsigned int d=0; d<dim; ++d)
-           present_u[d] = present_solution_values[q](d);         
+           present_u[d] = present_stokes_values[q](d);
 
          if (use_bdf2_scheme == true)
            {
@@ -1927,7 +1916,7 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
                                        +
                                        time_step *
                                        kappa * grad_phi_T[i] * grad_phi_T[j])
-                                      * fe_values.JxW(q);
+                                      * temperature_fe_values.JxW(q);
 
              for (unsigned int i=0; i<dofs_per_cell; ++i)
                local_rhs(i) += ((time_step + old_time_step) / old_time_step *
@@ -1954,7 +1943,7 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
                                 time_step *
                                 gamma_values[q] * phi_T[i])
                                *
-                               fe_values.JxW(q);
+                               temperature_fe_values.JxW(q);
            }
          else
            {
@@ -1964,7 +1953,7 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
                                        +
                                        time_step *
                                        kappa * grad_phi_T[i] * grad_phi_T[j])
-                                      * fe_values.JxW(q);
+                                      * temperature_fe_values.JxW(q);
       
              for (unsigned int i=0; i<dofs_per_cell; ++i)
                local_rhs(i) += (old_T * phi_T[i]
@@ -1979,7 +1968,7 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
                                 time_step *
                                 gamma_values[q] * phi_T[i])
                                *
-                               fe_values.JxW(q);
+                               temperature_fe_values.JxW(q);
            }
        }
       
@@ -1987,12 +1976,12 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
 
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        for (unsigned int j=0; j<dofs_per_cell; ++j)
-         system_matrix.add (local_dof_indices[i],
-                            local_dof_indices[j],
-                            local_matrix(i,j));
+         temperature_matrix.add (local_dof_indices[i],
+                                 local_dof_indices[j],
+                                 local_matrix(i,j));
 
       for (unsigned int i=0; i<dofs_per_cell; ++i)
-        system_rhs(local_dof_indices[i]) += local_rhs(i);
+        temperature_rhs(local_dof_indices[i]) += local_rhs(i);
     }
 }
 
@@ -2003,76 +1992,49 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
 template <int dim>
 void BoussinesqFlowProblem<dim>::solve ()
 {
-  solution = old_solution;
-
-                               // Use the BlockMatrixArray structure
-                               // for extracting only the upper left
-                               // 2x2 blocks from the matrix that will
-                               // be used for the solution of the
-                               // blocked system.
+  std::cout << "   Solving..." << std::endl;
+  
+                                  // Use the BlockMatrixArray structure
+                                  // for extracting only the upper left
+                                  // 2x2 blocks from the matrix that will
+                                  // be used for the solution of the
+                                  // blocked system.
   {
-    GrowingVectorMemory<Vector<double> > simple_mem;
-    BlockMatrixArray<double> stokes_submatrix(2, 2, simple_mem);
-
-    stokes_submatrix.enter(system_matrix.block(0,0),0,0);
-    stokes_submatrix.enter(system_matrix.block(0,1),0,1);
-    stokes_submatrix.enter(system_matrix.block(1,0),1,0);
-
-                               // Define some temporary vectors
-                               // for the solution process.
-                               // TODO: Can we somehow avoid copying
-                               // these vectors back and forth? I.e.
-                               // accessing the block vectors in a
-                               // similar way as the matrix with the
-                               // BlockMatrixArray class?
-    std::vector<unsigned int> block_sizes(2);
-    block_sizes[0] = solution.block(0).size();
-    block_sizes[1] = solution.block(1).size();
-
-    BlockVector<double> up_rhs(block_sizes);
-    BlockVector<double> up(block_sizes);
-
-    up_rhs.block(0) = system_rhs.block(0);
-    up_rhs.block(1) = system_rhs.block(1);
-
-                               // Set up inverse matrix for
-                               // pressure mass matrix
+                                    // Set up inverse matrix for
+                                    // pressure mass matrix
     InverseMatrix<SparseMatrix<double>,SparseILU<double> >
-      mp_inverse (preconditioner_matrix.block(1,1), *Mp_preconditioner);
+      mp_inverse (stokes_preconditioner_matrix.block(1,1), *Mp_preconditioner);
 
-                               // Set up block Schur preconditioner
-    /*BlockSchurPreconditioner<typename InnerPreconditioner<dim>::type,
-                           SparseILU<double> >
-      preconditioner (system_matrix, mp_inverse, *A_preconditioner);*/
+                                    // Set up block Schur preconditioner
+                                    /*BlockSchurPreconditioner<typename InnerPreconditioner<dim>::type,
+                                      SparseILU<double> >
+                                      preconditioner (stokes_matrix, mp_inverse, *A_preconditioner);*/
     BlockSchurPreconditioner<PreconditionerTrilinosAmg, SparseILU<double> >
-      preconditioner (system_matrix, mp_inverse, *Amg_preconditioner);
+      preconditioner (stokes_matrix, mp_inverse, *Amg_preconditioner);
 
-                               // Set up GMRES solver and
-                               // solve.
-    SolverControl solver_control (system_matrix.m(),
-                                 1e-6*system_rhs.l2_norm());
+                                    // Set up GMRES solver and
+                                    // solve.
+    SolverControl solver_control (stokes_matrix.m(),
+                                 1e-6*stokes_rhs.l2_norm());
 
     SolverGMRES<BlockVector<double> > gmres(solver_control,
-                       SolverGMRES<BlockVector<double> >::AdditionalData(100));
+                                           SolverGMRES<BlockVector<double> >::AdditionalData(100));
 
-    gmres.solve(stokes_submatrix, up, up_rhs, preconditioner);
+    gmres.solve(stokes_matrix, stokes_solution, stokes_rhs, preconditioner);
 
     std::cout << "   "
               << solver_control.last_step()
               << " GMRES iterations for Stokes subsystem."
               << std::endl;
              
-    solution.block(0) = up.block(0);
-    solution.block(1) = up.block(1);
-
-                               // Produce a constistent solution
-                               // field (we can't do this on the 'up'
-                               // vector since it does not have the
-                               // temperature component, but
-                               // hanging_node_constraints has
-                               // constraints also for the
-                               // temperature vector)
-    hanging_node_constraints.distribute (solution);
+                                    // Produce a constistent solution
+                                    // field (we can't do this on the 'up'
+                                    // vector since it does not have the
+                                    // temperature component, but
+                                    // hanging_node_constraints has
+                                    // constraints also for the
+                                    // temperature vector)
+    stokes_constraints.distribute (stokes_solution);
   }
 
                                   // TODO: determine limit of stability for
@@ -2082,29 +2044,32 @@ void BoussinesqFlowProblem<dim>::solve ()
   old_time_step = time_step;    
   time_step = GridTools::minimal_cell_diameter(triangulation) /
               std::max (get_maximal_velocity(), .05) / 4;
+  
+  temperature_solution = old_temperature_solution;
 
-  assemble_rhs_T ();
+
+  assemble_temperature_system ();
   {
 
-    SolverControl solver_control (system_matrix.block(2,2).m(),
-                                 1e-8*system_rhs.block(2).l2_norm());
+    SolverControl solver_control (temperature_matrix.m(),
+                                 1e-8*temperature_rhs.l2_norm());
     SolverCG<>   cg (solver_control);
     PreconditionSSOR<> preconditioner;
-    preconditioner.initialize (system_matrix.block(2,2), 1.2);
+    preconditioner.initialize (temperature_matrix, 1.2);
 
-    cg.solve (system_matrix.block(2,2), solution.block(2),
-             system_rhs.block(2), preconditioner);
+    cg.solve (temperature_matrix, temperature_solution,
+             temperature_rhs, preconditioner);
 
                                     // produce a consistent temperature field
-    hanging_node_constraints.distribute (solution);
+    temperature_constraints.distribute (temperature_solution);
 
     std::cout << "   "
               << solver_control.last_step()
               << " CG iterations for temperature."
               << std::endl;
     std::cout << "   Max temperature: "
-             << *std::max_element (solution.block(2).begin(),
-                                   solution.block(2).end())
+             << *std::max_element (temperature_solution.begin(),
+                                   temperature_solution.end())
              << std::endl;
   }
 }
@@ -2118,31 +2083,34 @@ void BoussinesqFlowProblem<dim>::output_results ()  const
   if (timestep_number % 10 != 0)
     return;
 
-  std::vector<std::string> solution_names (dim, "velocity");
-  solution_names.push_back ("p");
-  solution_names.push_back ("T");
+//TODO!!  
+//  Assert (false, ExcNotImplemented());
+  
+//   std::vector<std::string> solution_names (dim, "velocity");
+//   solution_names.push_back ("p");
+//   solution_names.push_back ("T");
 
-  DataOut<dim> data_out;
+//   DataOut<dim> data_out;
 
-  data_out.attach_dof_handler (dof_handler);
+//   data_out.attach_dof_handler (dof_handler);
 
-  std::vector<DataComponentInterpretation::DataComponentInterpretation>
-    data_component_interpretation
-    (dim+2, DataComponentInterpretation::component_is_scalar);
-  for (unsigned int i=0; i<dim; ++i)
-    data_component_interpretation[i]
-      = DataComponentInterpretation::component_is_part_of_vector;
+//   std::vector<DataComponentInterpretation::DataComponentInterpretation>
+//     data_component_interpretation
+//     (dim+2, DataComponentInterpretation::component_is_scalar);
+//   for (unsigned int i=0; i<dim; ++i)
+//     data_component_interpretation[i]
+//       = DataComponentInterpretation::component_is_part_of_vector;
 
-  data_out.add_data_vector (solution, solution_names,
-                           DataOut<dim>::type_dof_data,
-                           data_component_interpretation);
-  data_out.build_patches (degree);
+//   data_out.add_data_vector (solution, solution_names,
+//                         DataOut<dim>::type_dof_data,
+//                         data_component_interpretation);
+//   data_out.build_patches (degree);
 
-  std::ostringstream filename;
-  filename << "solution-" << Utilities::int_to_string(timestep_number, 4) << ".vtk";
+//   std::ostringstream filename;
+//   filename << "solution-" << Utilities::int_to_string(timestep_number, 4) << ".vtk";
 
-  std::ofstream output (filename.str().c_str());
-  data_out.write_vtk (output);
+//   std::ofstream output (filename.str().c_str());
+//   data_out.write_vtk (output);
 }
 
 
@@ -2153,14 +2121,11 @@ void BoussinesqFlowProblem<dim>::refine_mesh (const unsigned int max_grid_level)
 {
   Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
 
-  std::vector<bool> component_mask (dim+2, false);
-  component_mask[dim+1] = true;
-  KellyErrorEstimator<dim>::estimate (dof_handler,
-                                     QGauss<dim-1>(3),
+  KellyErrorEstimator<dim>::estimate (temperature_dof_handler,
+                                     QGauss<dim-1>(degree+1),
                                      typename FunctionMap<dim>::type(),
-                                     solution,
-                                     estimated_error_per_cell,
-                                     component_mask);
+                                     temperature_solution,
+                                     estimated_error_per_cell);
 
   GridRefinement::refine_and_coarsen_fixed_fraction (triangulation,
                                                     estimated_error_per_cell,
@@ -2172,15 +2137,15 @@ void BoussinesqFlowProblem<dim>::refine_mesh (const unsigned int max_grid_level)
       if (cell->has_children() == false)
        cell->clear_refine_flag ();
   
-  SolutionTransfer<dim, double> soltrans(dof_handler);
+  SolutionTransfer<dim, double> soltrans(temperature_dof_handler);
 
   triangulation.prepare_coarsening_and_refinement();
 
   std::vector<Vector<double> > x_solution (2);
-  x_solution[0].reinit (dof_handler.n_dofs());
-  x_solution[0] = solution;
-  x_solution[1].reinit (dof_handler.n_dofs());
-  x_solution[1] = old_solution;
+  x_solution[0].reinit (temperature_dof_handler.n_dofs());
+  x_solution[0] = temperature_solution;
+  x_solution[1].reinit (temperature_dof_handler.n_dofs());
+  x_solution[1] = old_temperature_solution;
 
   soltrans.prepare_for_coarsening_and_refinement(x_solution);
 
@@ -2188,15 +2153,16 @@ void BoussinesqFlowProblem<dim>::refine_mesh (const unsigned int max_grid_level)
   setup_dofs ();
 
   std::vector<Vector<double> > tmp (2);
-  tmp[0].reinit (dof_handler.n_dofs());
-  tmp[1].reinit (dof_handler.n_dofs());
+  tmp[0].reinit (temperature_dof_handler.n_dofs());
+  tmp[1].reinit (temperature_dof_handler.n_dofs());
   soltrans.interpolate(x_solution, tmp);
 
+  temperature_solution = tmp[0];
+                                  //TODO: this should be [1]!!!!!!!
+  old_temperature_solution = tmp[0];
+
   rebuild_matrices       = true;
   rebuild_preconditioner = true;
-
-  solution = tmp[0];
-  old_solution = tmp[0];
 }
 
 
@@ -2209,25 +2175,25 @@ double BoussinesqFlowProblem<dim>::get_maximal_velocity () const
   const unsigned int   n_q_points
     = quadrature_formula.size();
 
-  FEValues<dim> fe_values (fe, quadrature_formula,
+  FEValues<dim> fe_values (stokes_fe, quadrature_formula,
                            update_values);
-  std::vector<Vector<double> > solution_values(n_q_points,
-                                               Vector<double>(dim+2));
+  std::vector<Vector<double> > stokes_values(n_q_points,
+                                            Vector<double>(dim+1));
   double max_velocity = 0;
 
   typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
+    cell = stokes_dof_handler.begin_active(),
+    endc = stokes_dof_handler.end();
   for (; cell!=endc; ++cell)
     {
       fe_values.reinit (cell);
-      fe_values.get_function_values (solution, solution_values);
+      fe_values.get_function_values (stokes_solution, stokes_values);
 
       for (unsigned int q=0; q<n_q_points; ++q)
         {
           Tensor<1,dim> velocity;
           for (unsigned int i=0; i<dim; ++i)
-            velocity[i] = solution_values[q](i);
+            velocity[i] = stokes_values[q](i);
 
           max_velocity = std::max (max_velocity,
                                    velocity.norm());
@@ -2245,31 +2211,29 @@ template <int dim>
 double BoussinesqFlowProblem<dim>::get_maximal_temperature () const
 {
   QGauss<dim>   quadrature_formula(degree+2);
-  const unsigned int   n_q_points
-    = quadrature_formula.size();
+  const unsigned int   n_q_points = quadrature_formula.size();
 
-  FEValues<dim> fe_values (fe, quadrature_formula,
+  FEValues<dim> fe_values (temperature_fe, quadrature_formula,
                            update_values);
-  std::vector<Vector<double> > old_solution_values(n_q_points,
-                                                  Vector<double>(dim+2));
-  std::vector<Vector<double> > old_old_solution_values(n_q_points,
-                                                      Vector<double>(dim+2));
+  std::vector<double> old_temperature_values(n_q_points);
+  std::vector<double> old_old_temperature_values(n_q_points);
+  
   double max_temperature = 0;
 
   typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
+    cell = temperature_dof_handler.begin_active(),
+    endc = temperature_dof_handler.end();
   for (; cell!=endc; ++cell)
     {
       fe_values.reinit (cell);
-      fe_values.get_function_values (old_solution, old_solution_values);
-      fe_values.get_function_values (old_old_solution, old_old_solution_values);
+      fe_values.get_function_values (old_temperature_solution, old_temperature_values);
+      fe_values.get_function_values (old_old_temperature_solution, old_old_temperature_values);
 
       for (unsigned int q=0; q<n_q_points; ++q)
         {
           double temperature = 
-               (1. + time_step/old_time_step) * old_solution_values[q](dim+1)-
-               time_step/old_time_step * old_old_solution_values[q](dim+1);
+           (1. + time_step/old_time_step) * old_temperature_values[q]-
+           time_step/old_time_step * old_old_temperature_values[q];
 
           max_temperature = std::max (max_temperature,
                                      temperature);
@@ -2298,11 +2262,11 @@ void BoussinesqFlowProblem<dim>::run ()
   
   start_time_iteration:
 
-  VectorTools::project (dof_handler,
-                       hanging_node_constraints,
+  VectorTools::project (temperature_dof_handler,
+                       temperature_constraints,
                        QGauss<dim>(degree+2),
-                       InitialValues<dim>(),
-                       old_solution);
+                       TemperatureInitialValues<dim>(),
+                       old_temperature_solution);
   
   timestep_number = 0;
   double time = 0;
@@ -2314,10 +2278,11 @@ void BoussinesqFlowProblem<dim>::run ()
                << ", dt=" << time_step
                 << std::endl;
 
-      std::cout << "   Assembling..." << std::endl;
-      assemble_system ();
+      assemble_stokes_system ();
+
+      if (rebuild_preconditioner == true)
+       build_stokes_preconditioner ();
 
-      std::cout << "   Solving..." << std::endl;
       solve ();
 
       output_results ();
@@ -2338,8 +2303,8 @@ void BoussinesqFlowProblem<dim>::run ()
       time += time_step;
       ++timestep_number;
 
-      old_old_solution = old_solution;
-      old_solution     = solution;      
+      old_old_temperature_solution = old_temperature_solution;
+      old_temperature_solution     = temperature_solution;      
     }
   while (time <= 10);
 }

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.