]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Wrap everything into a namespace Step32 as was done for all other
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 23 Sep 2011 00:04:16 +0000 (00:04 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 23 Sep 2011 00:04:16 +0000 (00:04 +0000)
tutorial programs.

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

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

index 5981134d3cd2dc48b81c59e8af604f420ca19413..8f7c09449043fd028c50e725108af649ab178369 100644 (file)
@@ -94,7 +94,9 @@
 
 
 
-using namespace dealii;
+namespace Step32
+{
+  using namespace dealii;
 
 // @sect3{Equation data}
 
@@ -103,150 +105,150 @@ using namespace dealii;
 // these are exhaustively discussed in the
 // description of the testcase in the
 // introduction:
-namespace EquationData
-{
-  const double eta                   = 1e21;    /* Pa s       */
-  const double kappa                 = 1e-6;
-  const double reference_density     = 3300;    /* kg / m^3   */
-  const double reference_temperature = 293;     /* K          */
-  const double expansion_coefficient = 2e-5;    /* 1/K        */
-  const double specific_heat         = 1250;    /* J / K / kg */  //??
-  const double radiogenic_heating    = 7.4e-12; /* W / kg     */  //??
+  namespace EquationData
+  {
+    const double eta                   = 1e21;    /* Pa s       */
+    const double kappa                 = 1e-6;
+    const double reference_density     = 3300;    /* kg / m^3   */
+    const double reference_temperature = 293;     /* K          */
+    const double expansion_coefficient = 2e-5;    /* 1/K        */
+    const double specific_heat         = 1250;    /* J / K / kg */  //??
+    const double radiogenic_heating    = 7.4e-12; /* W / kg     */  //??
 
 
-  const double R0      = 6371000.-2890000.;     /* m          */
-  const double R1      = 6371000.-  35000.;     /* m          */
+    const double R0      = 6371000.-2890000.;     /* m          */
+    const double R1      = 6371000.-  35000.;     /* m          */
 
-  const double T0      = 4000+273;              /* K          */
-  const double T1      =  700+273;              /* K          */
+    const double T0      = 4000+273;              /* K          */
+    const double T1      =  700+273;              /* K          */
 
-  const double year_in_seconds  = 60*60*24*365.2425;
+    const double year_in_seconds  = 60*60*24*365.2425;
 
 //TODO: document in intro.dox
-  // scale not by R1-R0, but by a
-  // typical length scale, say 10km,
-  // of variation ("plume
-  // diameter"). this choice also
-  // roughly equilibrates the sizes
-  // of the velocity and pressure
-  // components of the solution
-  // vectors
-  const double pressure_scaling = eta / 10000;
-
-
-  double density (const double temperature)
-  {
-    return (reference_density *
-            (1 - expansion_coefficient * (temperature -
-                                          reference_temperature)));
-  }
+                                    // scale not by R1-R0, but by a
+                                    // typical length scale, say 10km,
+                                    // of variation ("plume
+                                    // diameter"). this choice also
+                                    // roughly equilibrates the sizes
+                                    // of the velocity and pressure
+                                    // components of the solution
+                                    // vectors
+    const double pressure_scaling = eta / 10000;
+
+
+    double density (const double temperature)
+    {
+      return (reference_density *
+             (1 - expansion_coefficient * (temperature -
+                                           reference_temperature)));
+    }
 
 
-  template <int dim>
-  Tensor<1,dim> gravity_vector (const Point<dim> &p)
-  {
+    template <int dim>
+    Tensor<1,dim> gravity_vector (const Point<dim> &p)
+    {
 // interpolate the following values with a physically realistic model:
 //    const double g0      = 10.7;                  /* m / s^2    */
 //    const double g1      = 9.81;                  /* m / s^2    */
 
-    const double r = p.norm();
-    return -(1.245e-6 * r + 7.714e13/r/r) * p / p.norm();
-  }
+      const double r = p.norm();
+      return -(1.245e-6 * r + 7.714e13/r/r) * p / p.norm();
+    }
 
 
 
-  template <int dim>
-  double adiabatic_pressure (const Point<dim> &p)
-  {
-    // the static, adiabatic pressure
-    // satisfies
-    //    dP/dr = -g rho
-
-    // assuming a constant density,
-    // we can integrate the pressure
-    // equation in depth to get that
-    // the adiabatic pressure equals
-    //   $P(r) = rho_0 \int_r^{R_1} g(r) dr$
-    //
-    // using the model for the
-    // gravity vector above, this
-    // yields the following formula:
-    const double r = p.norm();
-    return reference_density * (1./2 * 1.245e-6 * (R1*R1 - r*r) - 7.714e13 * (1./R1 - 1./r));
-  }
+    template <int dim>
+    double adiabatic_pressure (const Point<dim> &p)
+    {
+                                      // the static, adiabatic pressure
+                                      // satisfies
+                                      //    dP/dr = -g rho
+
+                                      // assuming a constant density,
+                                      // we can integrate the pressure
+                                      // equation in depth to get that
+                                      // the adiabatic pressure equals
+                                      //   $P(r) = rho_0 \int_r^{R_1} g(r) dr$
+                                      //
+                                      // using the model for the
+                                      // gravity vector above, this
+                                      // yields the following formula:
+      const double r = p.norm();
+      return reference_density * (1./2 * 1.245e-6 * (R1*R1 - r*r) - 7.714e13 * (1./R1 - 1./r));
+    }
 
 
-  template <int dim>
-  double adiabatic_temperature (const Point<dim> &p)
-  {
-    // the static, adiabatic
-    // temperature satisfies
-    //    $dT/dr = -T alpha/c_P g$
+    template <int dim>
+    double adiabatic_temperature (const Point<dim> &p)
+    {
+                                      // the static, adiabatic
+                                      // temperature satisfies
+                                      //    $dT/dr = -T alpha/c_P g$
 
-    // let's assume constant gravity,
-    // then we get by integration
-    const double r = p.norm();
+                                      // let's assume constant gravity,
+                                      // then we get by integration
+      const double r = p.norm();
 
-    return T1 * std::exp(-expansion_coefficient * 9.81 / specific_heat * (r-R1));
-  }
+      return T1 * std::exp(-expansion_coefficient * 9.81 / specific_heat * (r-R1));
+    }
 
 
 
-  template <int dim>
-  class TemperatureInitialValues : public Function<dim>
-  {
-    public:
-      TemperatureInitialValues () : Function<dim>(1) {}
+    template <int dim>
+    class TemperatureInitialValues : public Function<dim>
+    {
+      public:
+       TemperatureInitialValues () : Function<dim>(1) {}
 
-      virtual double value (const Point<dim>   &p,
-                            const unsigned int  component = 0) const;
+       virtual double value (const Point<dim>   &p,
+                             const unsigned int  component = 0) const;
 
-      virtual void vector_value (const Point<dim> &p,
-                                 Vector<double>   &value) const;
-  };
+       virtual void vector_value (const Point<dim> &p,
+                                  Vector<double>   &value) const;
+    };
 
 
 
-  template <int dim>
-  double
-  TemperatureInitialValues<dim>::value (const Point<dim>  &p,
-                                        const unsigned int) const
-  {
-    const double r = p.norm();
-    const double h = R1-R0;
-
-    // s = fraction of the way from
-    // the inner to the outer
-    // boundary; 0<=s<=1
-    const double s = (r-R0)/h;
-
-    /* now compute an angular variation of the linear temperature field by
-       stretching the variable s appropriately. note that the following
-       formula leaves the end points s=0 and s=1 fixed, but stretches the
-       region in between depending on the angle phi=atan2(x,y).
-
-       For a plot, see
-       http://www.wolframalpha.com/input/?i=plot+%28%282*sqrt%28x^2%2By^2%29-1%29%2B0.2*%282*sqrt%28x^2%2By^2%29-1%29*%281-%282*sqrt%28x^2%2By^2%29-1%29%29*sin%286*atan2%28x%2Cy%29%29%29%2C+x%3D-1+to+1%2C+y%3D-1+to+1
-    */
-    const double scale = (dim==3)?std::max(0.0,cos(3.14159*abs(p(2)/R1))):1.0;
-    const double phi   = std::atan2(p(0),p(1));
-    const double s_mod = s
-                         +
-                         0.2 * s * (1-s) * std::sin(6*phi) * scale;
-
-    return T0*(1.0-s_mod) + T1*s_mod;
-  }
+    template <int dim>
+    double
+    TemperatureInitialValues<dim>::value (const Point<dim>  &p,
+                                         const unsigned int) const
+    {
+      const double r = p.norm();
+      const double h = R1-R0;
+
+                                      // s = fraction of the way from
+                                      // the inner to the outer
+                                      // boundary; 0<=s<=1
+      const double s = (r-R0)/h;
+
+                                      /* now compute an angular variation of the linear temperature field by
+                                         stretching the variable s appropriately. note that the following
+                                         formula leaves the end points s=0 and s=1 fixed, but stretches the
+                                         region in between depending on the angle phi=atan2(x,y).
+
+                                         For a plot, see
+                                         http://www.wolframalpha.com/input/?i=plot+%28%282*sqrt%28x^2%2By^2%29-1%29%2B0.2*%282*sqrt%28x^2%2By^2%29-1%29*%281-%282*sqrt%28x^2%2By^2%29-1%29%29*sin%286*atan2%28x%2Cy%29%29%29%2C+x%3D-1+to+1%2C+y%3D-1+to+1
+                                      */
+      const double scale = (dim==3)?std::max(0.0,cos(3.14159*abs(p(2)/R1))):1.0;
+      const double phi   = std::atan2(p(0),p(1));
+      const double s_mod = s
+                          +
+                          0.2 * s * (1-s) * std::sin(6*phi) * scale;
+
+      return T0*(1.0-s_mod) + T1*s_mod;
+    }
 
 
-  template <int dim>
-  void
-  TemperatureInitialValues<dim>::vector_value (const Point<dim> &p,
-                                               Vector<double>   &values) const
-  {
-    for (unsigned int c=0; c<this->n_components; ++c)
-      values(c) = TemperatureInitialValues<dim>::value (p, c);
+    template <int dim>
+    void
+    TemperatureInitialValues<dim>::vector_value (const Point<dim> &p,
+                                                Vector<double>   &values) const
+    {
+      for (unsigned int c=0; c<this->n_components; ++c)
+       values(c) = TemperatureInitialValues<dim>::value (p, c);
+    }
   }
-}
 
 
 
@@ -275,73 +277,73 @@ namespace EquationData
 // respective variable in the
 // BlockSchurPreconditioner class another
 // name.
-namespace LinearSolvers
-{
-  template <class PreconditionerA, class PreconditionerMp>
-  class RightPrecond : public Subscriptor
+  namespace LinearSolvers
   {
-    public:
-      RightPrecond (
-        const TrilinosWrappers::BlockSparseMatrix  &S,
-        const TrilinosWrappers::BlockSparseMatrix  &Spre,
-        const PreconditionerMp                     &Mppreconditioner,
-        const PreconditionerA                      &Apreconditioner,
-        const bool                                  do_solve_A_in = true)
-        :
-        stokes_matrix     (&S),
-        stokes_preconditioner_matrix     (&Spre),
-        mp_preconditioner (Mppreconditioner),
-        a_preconditioner  (Apreconditioner),
-        do_solve_A        (do_solve_A_in)
-      {}
-
-      void solve_S(TrilinosWrappers::MPI::Vector &dst,
-                   const TrilinosWrappers::MPI::Vector &src) const
-      {
-        SolverControl cn(5000, 1e-5);
-
-        TrilinosWrappers::SolverCG solver(cn);
-
-        solver.solve(stokes_preconditioner_matrix->block(1,1),
-                     dst, src,
-                     mp_preconditioner);
-
-        dst*=-1.0;
-      }
-
-      void solve_A(TrilinosWrappers::MPI::Vector &dst,
-                   const TrilinosWrappers::MPI::Vector &src) const
-      {
-        SolverControl cn(5000, src.l2_norm()*1e-2);
-        TrilinosWrappers::SolverCG solver(cn);
-        solver.solve(stokes_matrix->block(0,0), dst, src, a_preconditioner);
-      }
-
-      void vmult (TrilinosWrappers::MPI::BlockVector       &dst,
-                  const TrilinosWrappers::MPI::BlockVector &src) const
-      {
-        TrilinosWrappers::MPI::Vector utmp(src.block(0));
-
-        solve_S(dst.block(1), src.block(1));
-
-        stokes_matrix->block(0,1).vmult(utmp, dst.block(1)); //B^T
-        utmp*=-1.0;
-        utmp.add(src.block(0));
-
-        if (do_solve_A == true)
-          solve_A(dst.block(0), utmp);
-        else
-          a_preconditioner.vmult (dst.block(0), utmp);
-      }
-
-    private:
-      const SmartPointer<const TrilinosWrappers::BlockSparseMatrix> stokes_matrix;
-      const SmartPointer<const TrilinosWrappers::BlockSparseMatrix> stokes_preconditioner_matrix;
-      const PreconditionerMp &mp_preconditioner;
-      const PreconditionerA  &a_preconditioner;
-      const bool do_solve_A;
-  };
-}
+    template <class PreconditionerA, class PreconditionerMp>
+    class RightPrecond : public Subscriptor
+    {
+      public:
+       RightPrecond (
+         const TrilinosWrappers::BlockSparseMatrix  &S,
+         const TrilinosWrappers::BlockSparseMatrix  &Spre,
+         const PreconditionerMp                     &Mppreconditioner,
+         const PreconditionerA                      &Apreconditioner,
+         const bool                                  do_solve_A_in = true)
+                       :
+                       stokes_matrix     (&S),
+                       stokes_preconditioner_matrix     (&Spre),
+                       mp_preconditioner (Mppreconditioner),
+                       a_preconditioner  (Apreconditioner),
+                       do_solve_A        (do_solve_A_in)
+         {}
+
+       void solve_S(TrilinosWrappers::MPI::Vector &dst,
+                    const TrilinosWrappers::MPI::Vector &src) const
+         {
+           SolverControl cn(5000, 1e-5);
+
+           TrilinosWrappers::SolverCG solver(cn);
+
+           solver.solve(stokes_preconditioner_matrix->block(1,1),
+                        dst, src,
+                        mp_preconditioner);
+
+           dst*=-1.0;
+         }
+
+       void solve_A(TrilinosWrappers::MPI::Vector &dst,
+                    const TrilinosWrappers::MPI::Vector &src) const
+         {
+           SolverControl cn(5000, src.l2_norm()*1e-2);
+           TrilinosWrappers::SolverCG solver(cn);
+           solver.solve(stokes_matrix->block(0,0), dst, src, a_preconditioner);
+         }
+
+       void vmult (TrilinosWrappers::MPI::BlockVector       &dst,
+                   const TrilinosWrappers::MPI::BlockVector &src) const
+         {
+           TrilinosWrappers::MPI::Vector utmp(src.block(0));
+
+           solve_S(dst.block(1), src.block(1));
+
+           stokes_matrix->block(0,1).vmult(utmp, dst.block(1)); //B^T
+           utmp*=-1.0;
+           utmp.add(src.block(0));
+
+           if (do_solve_A == true)
+             solve_A(dst.block(0), utmp);
+           else
+             a_preconditioner.vmult (dst.block(0), utmp);
+         }
+
+      private:
+       const SmartPointer<const TrilinosWrappers::BlockSparseMatrix> stokes_matrix;
+       const SmartPointer<const TrilinosWrappers::BlockSparseMatrix> stokes_preconditioner_matrix;
+       const PreconditionerMp &mp_preconditioner;
+       const PreconditionerA  &a_preconditioner;
+       const bool do_solve_A;
+    };
+  }
 
 
 
@@ -407,398 +409,398 @@ namespace LinearSolvers
 // itself), and provide some additional
 // vector fields that are used to improve
 // performance of assembly.
-namespace Assembly
-{
-  namespace Scratch
+  namespace Assembly
   {
-    template <int dim>
-    struct StokesPreconditioner
+    namespace Scratch
     {
-      StokesPreconditioner (const FiniteElement<dim> &stokes_fe,
-                            const Quadrature<dim>    &stokes_quadrature,
-                            const Mapping<dim>       &mapping,
-                            const UpdateFlags         update_flags);
-      StokesPreconditioner (const StokesPreconditioner &data);
-
-      FEValues<dim>               stokes_fe_values;
-
-      std::vector<SymmetricTensor<2,dim> > grads_phi_u;
-      std::vector<double>                  phi_p;
-    };
-
-    template <int dim>
-    StokesPreconditioner<dim>::
-    StokesPreconditioner (const FiniteElement<dim> &stokes_fe,
-                          const Quadrature<dim>    &stokes_quadrature,
-                          const Mapping<dim>       &mapping,
-                          const UpdateFlags         update_flags)
-      :
-      stokes_fe_values (mapping, stokes_fe, stokes_quadrature,
-                        update_flags),
-      grads_phi_u (stokes_fe.dofs_per_cell),
-      phi_p (stokes_fe.dofs_per_cell)
-    {}
+      template <int dim>
+      struct StokesPreconditioner
+      {
+         StokesPreconditioner (const FiniteElement<dim> &stokes_fe,
+                               const Quadrature<dim>    &stokes_quadrature,
+                               const Mapping<dim>       &mapping,
+                               const UpdateFlags         update_flags);
+         StokesPreconditioner (const StokesPreconditioner &data);
 
+         FEValues<dim>               stokes_fe_values;
 
+         std::vector<SymmetricTensor<2,dim> > grads_phi_u;
+         std::vector<double>                  phi_p;
+      };
 
-    template <int dim>
-    StokesPreconditioner<dim>::
-    StokesPreconditioner (const StokesPreconditioner &scratch)
-      :
-      stokes_fe_values (scratch.stokes_fe_values.get_mapping(),
-                        scratch.stokes_fe_values.get_fe(),
-                        scratch.stokes_fe_values.get_quadrature(),
-                        scratch.stokes_fe_values.get_update_flags()),
-      grads_phi_u (scratch.grads_phi_u),
-      phi_p (scratch.phi_p)
-    {}
-
-
-
-    // Observe that we derive the
-    // StokesSystem scratch array from the
-    // StokesPreconditioner array. We do this
-    // because all the objects that are
-    // necessary for the assembly of the
-    // preconditioner are also needed for the
-    // actual matrix system and right hand
-    // side, plus some extra data. This makes
-    // the program more compact. Note also
-    // that the assembly of the Stokes system
-    // and the temperature right hand side
-    // further down requires data from
-    // temperature and velocity,
-    // respectively, so we actually need two
-    // FEValues objects for those two cases.
-    template <int dim>
-    struct StokesSystem : public StokesPreconditioner<dim>
-    {
-      StokesSystem (const FiniteElement<dim> &stokes_fe,
-                    const Mapping<dim>       &mapping,
-                    const Quadrature<dim>    &stokes_quadrature,
-                    const UpdateFlags         stokes_update_flags,
-                    const FiniteElement<dim> &temperature_fe,
-                    const UpdateFlags         temperature_update_flags);
+      template <int dim>
+      StokesPreconditioner<dim>::
+      StokesPreconditioner (const FiniteElement<dim> &stokes_fe,
+                           const Quadrature<dim>    &stokes_quadrature,
+                           const Mapping<dim>       &mapping,
+                           const UpdateFlags         update_flags)
+                     :
+                     stokes_fe_values (mapping, stokes_fe, stokes_quadrature,
+                                       update_flags),
+                     grads_phi_u (stokes_fe.dofs_per_cell),
+                     phi_p (stokes_fe.dofs_per_cell)
+      {}
 
-      StokesSystem (const StokesSystem<dim> &data);
 
-      FEValues<dim>  temperature_fe_values;
 
-      std::vector<Tensor<1,dim> >          phi_u;
-      std::vector<SymmetricTensor<2,dim> > grads_phi_u;
-      std::vector<double>                  div_phi_u;
+      template <int dim>
+      StokesPreconditioner<dim>::
+      StokesPreconditioner (const StokesPreconditioner &scratch)
+                     :
+                     stokes_fe_values (scratch.stokes_fe_values.get_mapping(),
+                                       scratch.stokes_fe_values.get_fe(),
+                                       scratch.stokes_fe_values.get_quadrature(),
+                                       scratch.stokes_fe_values.get_update_flags()),
+                     grads_phi_u (scratch.grads_phi_u),
+                     phi_p (scratch.phi_p)
+      {}
 
-      std::vector<double>                  old_temperature_values;
-    };
 
 
-    template <int dim>
-    StokesSystem<dim>::
-    StokesSystem (const FiniteElement<dim> &stokes_fe,
-                  const Mapping<dim>       &mapping,
-                  const Quadrature<dim>    &stokes_quadrature,
-                  const UpdateFlags         stokes_update_flags,
-                  const FiniteElement<dim> &temperature_fe,
-                  const UpdateFlags         temperature_update_flags)
-      :
-      StokesPreconditioner<dim> (stokes_fe, stokes_quadrature,
-                                 mapping,
-                                 stokes_update_flags),
-      temperature_fe_values (mapping, temperature_fe, stokes_quadrature,
-                             temperature_update_flags),
-      phi_u (stokes_fe.dofs_per_cell),
-      grads_phi_u (stokes_fe.dofs_per_cell),
-      div_phi_u (stokes_fe.dofs_per_cell),
-      old_temperature_values (stokes_quadrature.size())
-    {}
+                                      // Observe that we derive the
+                                      // StokesSystem scratch array from the
+                                      // StokesPreconditioner array. We do this
+                                      // because all the objects that are
+                                      // necessary for the assembly of the
+                                      // preconditioner are also needed for the
+                                      // actual matrix system and right hand
+                                      // side, plus some extra data. This makes
+                                      // the program more compact. Note also
+                                      // that the assembly of the Stokes system
+                                      // and the temperature right hand side
+                                      // further down requires data from
+                                      // temperature and velocity,
+                                      // respectively, so we actually need two
+                                      // FEValues objects for those two cases.
+      template <int dim>
+      struct StokesSystem : public StokesPreconditioner<dim>
+      {
+         StokesSystem (const FiniteElement<dim> &stokes_fe,
+                       const Mapping<dim>       &mapping,
+                       const Quadrature<dim>    &stokes_quadrature,
+                       const UpdateFlags         stokes_update_flags,
+                       const FiniteElement<dim> &temperature_fe,
+                       const UpdateFlags         temperature_update_flags);
 
+         StokesSystem (const StokesSystem<dim> &data);
 
-    template <int dim>
-    StokesSystem<dim>::
-    StokesSystem (const StokesSystem<dim> &scratch)
-      :
-      StokesPreconditioner<dim> (scratch),
-      temperature_fe_values (scratch.temperature_fe_values.get_mapping(),
-                             scratch.temperature_fe_values.get_fe(),
-                             scratch.temperature_fe_values.get_quadrature(),
-                             scratch.temperature_fe_values.get_update_flags()),
-      phi_u (scratch.phi_u),
-      grads_phi_u (scratch.grads_phi_u),
-      div_phi_u (scratch.div_phi_u),
-      old_temperature_values (scratch.old_temperature_values)
-    {}
+         FEValues<dim>  temperature_fe_values;
 
+         std::vector<Tensor<1,dim> >          phi_u;
+         std::vector<SymmetricTensor<2,dim> > grads_phi_u;
+         std::vector<double>                  div_phi_u;
 
+         std::vector<double>                  old_temperature_values;
+      };
 
-    template <int dim>
-    struct TemperatureMatrix
-    {
-      TemperatureMatrix (const FiniteElement<dim> &temperature_fe,
-                         const Mapping<dim>       &mapping,
-                         const Quadrature<dim>    &temperature_quadrature);
-      TemperatureMatrix (const TemperatureMatrix &data);
 
-      FEValues<dim>               temperature_fe_values;
+      template <int dim>
+      StokesSystem<dim>::
+      StokesSystem (const FiniteElement<dim> &stokes_fe,
+                   const Mapping<dim>       &mapping,
+                   const Quadrature<dim>    &stokes_quadrature,
+                   const UpdateFlags         stokes_update_flags,
+                   const FiniteElement<dim> &temperature_fe,
+                   const UpdateFlags         temperature_update_flags)
+                     :
+                     StokesPreconditioner<dim> (stokes_fe, stokes_quadrature,
+                                                mapping,
+                                                stokes_update_flags),
+                     temperature_fe_values (mapping, temperature_fe, stokes_quadrature,
+                                            temperature_update_flags),
+                     phi_u (stokes_fe.dofs_per_cell),
+                     grads_phi_u (stokes_fe.dofs_per_cell),
+                     div_phi_u (stokes_fe.dofs_per_cell),
+                     old_temperature_values (stokes_quadrature.size())
+      {}
 
-      std::vector<double>         phi_T;
-      std::vector<Tensor<1,dim> > grad_phi_T;
-    };
 
-    template <int dim>
-    TemperatureMatrix<dim>::
-    TemperatureMatrix (const FiniteElement<dim> &temperature_fe,
-                       const Mapping<dim>       &mapping,
-                       const Quadrature<dim>    &temperature_quadrature)
-      :
-      temperature_fe_values (mapping,
-                             temperature_fe, temperature_quadrature,
-                             update_values    | update_gradients |
-                             update_JxW_values),
-      phi_T (temperature_fe.dofs_per_cell),
-      grad_phi_T (temperature_fe.dofs_per_cell)
-    {}
+      template <int dim>
+      StokesSystem<dim>::
+      StokesSystem (const StokesSystem<dim> &scratch)
+                     :
+                     StokesPreconditioner<dim> (scratch),
+                     temperature_fe_values (scratch.temperature_fe_values.get_mapping(),
+                                            scratch.temperature_fe_values.get_fe(),
+                                            scratch.temperature_fe_values.get_quadrature(),
+                                            scratch.temperature_fe_values.get_update_flags()),
+                     phi_u (scratch.phi_u),
+                     grads_phi_u (scratch.grads_phi_u),
+                     div_phi_u (scratch.div_phi_u),
+                     old_temperature_values (scratch.old_temperature_values)
+      {}
 
 
-    template <int dim>
-    TemperatureMatrix<dim>::
-    TemperatureMatrix (const TemperatureMatrix &scratch)
-      :
-      temperature_fe_values (scratch.temperature_fe_values.get_mapping(),
-                             scratch.temperature_fe_values.get_fe(),
-                             scratch.temperature_fe_values.get_quadrature(),
-                             scratch.temperature_fe_values.get_update_flags()),
-      phi_T (scratch.phi_T),
-      grad_phi_T (scratch.grad_phi_T)
-    {}
 
+      template <int dim>
+      struct TemperatureMatrix
+      {
+         TemperatureMatrix (const FiniteElement<dim> &temperature_fe,
+                            const Mapping<dim>       &mapping,
+                            const Quadrature<dim>    &temperature_quadrature);
+         TemperatureMatrix (const TemperatureMatrix &data);
 
-    template <int dim>
-    struct TemperatureRHS
-    {
-      TemperatureRHS (const FiniteElement<dim> &temperature_fe,
-                      const FiniteElement<dim> &stokes_fe,
-                      const Mapping<dim>       &mapping,
-                      const Quadrature<dim>    &quadrature);
-      TemperatureRHS (const TemperatureRHS &data);
+         FEValues<dim>               temperature_fe_values;
 
-      FEValues<dim>               temperature_fe_values;
-      FEValues<dim>               stokes_fe_values;
+         std::vector<double>         phi_T;
+         std::vector<Tensor<1,dim> > grad_phi_T;
+      };
 
-      std::vector<double>         phi_T;
-      std::vector<Tensor<1,dim> > grad_phi_T;
+      template <int dim>
+      TemperatureMatrix<dim>::
+      TemperatureMatrix (const FiniteElement<dim> &temperature_fe,
+                        const Mapping<dim>       &mapping,
+                        const Quadrature<dim>    &temperature_quadrature)
+                     :
+                     temperature_fe_values (mapping,
+                                            temperature_fe, temperature_quadrature,
+                                            update_values    | update_gradients |
+                                            update_JxW_values),
+                     phi_T (temperature_fe.dofs_per_cell),
+                     grad_phi_T (temperature_fe.dofs_per_cell)
+      {}
 
-      std::vector<Tensor<1,dim> > old_velocity_values;
-      std::vector<Tensor<1,dim> > old_old_velocity_values;
 
-      std::vector<SymmetricTensor<2,dim> > old_strain_rates;
-      std::vector<SymmetricTensor<2,dim> > old_old_strain_rates;
+      template <int dim>
+      TemperatureMatrix<dim>::
+      TemperatureMatrix (const TemperatureMatrix &scratch)
+                     :
+                     temperature_fe_values (scratch.temperature_fe_values.get_mapping(),
+                                            scratch.temperature_fe_values.get_fe(),
+                                            scratch.temperature_fe_values.get_quadrature(),
+                                            scratch.temperature_fe_values.get_update_flags()),
+                     phi_T (scratch.phi_T),
+                     grad_phi_T (scratch.grad_phi_T)
+      {}
 
-      std::vector<double>         old_temperature_values;
-      std::vector<double>         old_old_temperature_values;
-      std::vector<Tensor<1,dim> > old_temperature_grads;
-      std::vector<Tensor<1,dim> > old_old_temperature_grads;
-      std::vector<double>         old_temperature_laplacians;
-      std::vector<double>         old_old_temperature_laplacians;
-    };
 
-    template <int dim>
-    TemperatureRHS<dim>::
-    TemperatureRHS (const FiniteElement<dim> &temperature_fe,
-                    const FiniteElement<dim> &stokes_fe,
-                    const Mapping<dim>       &mapping,
-                    const Quadrature<dim>    &quadrature)
-      :
-      temperature_fe_values (mapping,
-                             temperature_fe, quadrature,
-                             update_values    |
-                             update_gradients |
-                             update_hessians  |
-                             update_quadrature_points |
-                             update_JxW_values),
-      stokes_fe_values (mapping,
-                        stokes_fe, quadrature,
-                        update_values | update_gradients),
-      phi_T (temperature_fe.dofs_per_cell),
-      grad_phi_T (temperature_fe.dofs_per_cell),
-
-      old_velocity_values (quadrature.size()),
-      old_old_velocity_values (quadrature.size()),
-      old_strain_rates (quadrature.size()),
-      old_old_strain_rates (quadrature.size()),
-
-      old_temperature_values (quadrature.size()),
-      old_old_temperature_values(quadrature.size()),
-      old_temperature_grads(quadrature.size()),
-      old_old_temperature_grads(quadrature.size()),
-      old_temperature_laplacians(quadrature.size()),
-      old_old_temperature_laplacians(quadrature.size())
-    {}
+      template <int dim>
+      struct TemperatureRHS
+      {
+         TemperatureRHS (const FiniteElement<dim> &temperature_fe,
+                         const FiniteElement<dim> &stokes_fe,
+                         const Mapping<dim>       &mapping,
+                         const Quadrature<dim>    &quadrature);
+         TemperatureRHS (const TemperatureRHS &data);
+
+         FEValues<dim>               temperature_fe_values;
+         FEValues<dim>               stokes_fe_values;
+
+         std::vector<double>         phi_T;
+         std::vector<Tensor<1,dim> > grad_phi_T;
+
+         std::vector<Tensor<1,dim> > old_velocity_values;
+         std::vector<Tensor<1,dim> > old_old_velocity_values;
+
+         std::vector<SymmetricTensor<2,dim> > old_strain_rates;
+         std::vector<SymmetricTensor<2,dim> > old_old_strain_rates;
+
+         std::vector<double>         old_temperature_values;
+         std::vector<double>         old_old_temperature_values;
+         std::vector<Tensor<1,dim> > old_temperature_grads;
+         std::vector<Tensor<1,dim> > old_old_temperature_grads;
+         std::vector<double>         old_temperature_laplacians;
+         std::vector<double>         old_old_temperature_laplacians;
+      };
+
+      template <int dim>
+      TemperatureRHS<dim>::
+      TemperatureRHS (const FiniteElement<dim> &temperature_fe,
+                     const FiniteElement<dim> &stokes_fe,
+                     const Mapping<dim>       &mapping,
+                     const Quadrature<dim>    &quadrature)
+                     :
+                     temperature_fe_values (mapping,
+                                            temperature_fe, quadrature,
+                                            update_values    |
+                                            update_gradients |
+                                            update_hessians  |
+                                            update_quadrature_points |
+                                            update_JxW_values),
+                     stokes_fe_values (mapping,
+                                       stokes_fe, quadrature,
+                                       update_values | update_gradients),
+                     phi_T (temperature_fe.dofs_per_cell),
+                     grad_phi_T (temperature_fe.dofs_per_cell),
+
+                     old_velocity_values (quadrature.size()),
+                     old_old_velocity_values (quadrature.size()),
+                     old_strain_rates (quadrature.size()),
+                     old_old_strain_rates (quadrature.size()),
+
+                     old_temperature_values (quadrature.size()),
+                     old_old_temperature_values(quadrature.size()),
+                     old_temperature_grads(quadrature.size()),
+                     old_old_temperature_grads(quadrature.size()),
+                     old_temperature_laplacians(quadrature.size()),
+                     old_old_temperature_laplacians(quadrature.size())
+      {}
 
 
-    template <int dim>
-    TemperatureRHS<dim>::
-    TemperatureRHS (const TemperatureRHS &scratch)
-      :
-      temperature_fe_values (scratch.temperature_fe_values.get_mapping(),
-                             scratch.temperature_fe_values.get_fe(),
-                             scratch.temperature_fe_values.get_quadrature(),
-                             scratch.temperature_fe_values.get_update_flags()),
-      stokes_fe_values (scratch.stokes_fe_values.get_mapping(),
-                        scratch.stokes_fe_values.get_fe(),
-                        scratch.stokes_fe_values.get_quadrature(),
-                        scratch.stokes_fe_values.get_update_flags()),
-      phi_T (scratch.phi_T),
-      grad_phi_T (scratch.grad_phi_T),
-
-      old_velocity_values (scratch.old_velocity_values),
-      old_old_velocity_values (scratch.old_old_velocity_values),
-      old_strain_rates (scratch.old_strain_rates),
-      old_old_strain_rates (scratch.old_old_strain_rates),
-
-      old_temperature_values (scratch.old_temperature_values),
-      old_old_temperature_values (scratch.old_old_temperature_values),
-      old_temperature_grads (scratch.old_temperature_grads),
-      old_old_temperature_grads (scratch.old_old_temperature_grads),
-      old_temperature_laplacians (scratch.old_temperature_laplacians),
-      old_old_temperature_laplacians (scratch.old_old_temperature_laplacians)
-    {}
-  }
+      template <int dim>
+      TemperatureRHS<dim>::
+      TemperatureRHS (const TemperatureRHS &scratch)
+                     :
+                     temperature_fe_values (scratch.temperature_fe_values.get_mapping(),
+                                            scratch.temperature_fe_values.get_fe(),
+                                            scratch.temperature_fe_values.get_quadrature(),
+                                            scratch.temperature_fe_values.get_update_flags()),
+                     stokes_fe_values (scratch.stokes_fe_values.get_mapping(),
+                                       scratch.stokes_fe_values.get_fe(),
+                                       scratch.stokes_fe_values.get_quadrature(),
+                                       scratch.stokes_fe_values.get_update_flags()),
+                     phi_T (scratch.phi_T),
+                     grad_phi_T (scratch.grad_phi_T),
+
+                     old_velocity_values (scratch.old_velocity_values),
+                     old_old_velocity_values (scratch.old_old_velocity_values),
+                     old_strain_rates (scratch.old_strain_rates),
+                     old_old_strain_rates (scratch.old_old_strain_rates),
+
+                     old_temperature_values (scratch.old_temperature_values),
+                     old_old_temperature_values (scratch.old_old_temperature_values),
+                     old_temperature_grads (scratch.old_temperature_grads),
+                     old_old_temperature_grads (scratch.old_old_temperature_grads),
+                     old_temperature_laplacians (scratch.old_temperature_laplacians),
+                     old_old_temperature_laplacians (scratch.old_old_temperature_laplacians)
+      {}
+    }
 
 
-  // The CopyData arrays are similar to the
-  // Scratch arrays. They provide a
-  // constructor, a copy operation, and
-  // some arrays for local matrix, local
-  // vectors and the relation between local
-  // and global degrees of freedom (a.k.a.
-  // <code>local_dof_indices</code>).
-  namespace CopyData
-  {
-    template <int dim>
-    struct StokesPreconditioner
+                                    // The CopyData arrays are similar to the
+                                    // Scratch arrays. They provide a
+                                    // constructor, a copy operation, and
+                                    // some arrays for local matrix, local
+                                    // vectors and the relation between local
+                                    // and global degrees of freedom (a.k.a.
+                                    // <code>local_dof_indices</code>).
+    namespace CopyData
     {
-      StokesPreconditioner (const FiniteElement<dim> &stokes_fe);
-      StokesPreconditioner (const StokesPreconditioner &data);
-
-      FullMatrix<double>          local_matrix;
-      std::vector<unsigned int>   local_dof_indices;
-    };
-
-    template <int dim>
-    StokesPreconditioner<dim>::
-    StokesPreconditioner (const FiniteElement<dim> &stokes_fe)
-      :
-      local_matrix (stokes_fe.dofs_per_cell,
-                    stokes_fe.dofs_per_cell),
-      local_dof_indices (stokes_fe.dofs_per_cell)
-    {}
-
-
-
-    template <int dim>
-    StokesPreconditioner<dim>::
-    StokesPreconditioner (const StokesPreconditioner &data)
-      :
-      local_matrix (data.local_matrix),
-      local_dof_indices (data.local_dof_indices)
-    {}
-
+      template <int dim>
+      struct StokesPreconditioner
+      {
+         StokesPreconditioner (const FiniteElement<dim> &stokes_fe);
+         StokesPreconditioner (const StokesPreconditioner &data);
+
+         FullMatrix<double>          local_matrix;
+         std::vector<unsigned int>   local_dof_indices;
+      };
+
+      template <int dim>
+      StokesPreconditioner<dim>::
+      StokesPreconditioner (const FiniteElement<dim> &stokes_fe)
+                     :
+                     local_matrix (stokes_fe.dofs_per_cell,
+                                   stokes_fe.dofs_per_cell),
+                     local_dof_indices (stokes_fe.dofs_per_cell)
+      {}
 
 
-    template <int dim>
-    struct StokesSystem : public StokesPreconditioner<dim>
-    {
-      StokesSystem (const FiniteElement<dim> &stokes_fe);
-      StokesSystem (const StokesSystem<dim> &data);
 
-      Vector<double> local_rhs;
-    };
+      template <int dim>
+      StokesPreconditioner<dim>::
+      StokesPreconditioner (const StokesPreconditioner &data)
+                     :
+                     local_matrix (data.local_matrix),
+                     local_dof_indices (data.local_dof_indices)
+      {}
 
 
-    template <int dim>
-    StokesSystem<dim>::
-    StokesSystem (const FiniteElement<dim> &stokes_fe)
-      :
-      StokesPreconditioner<dim> (stokes_fe),
-      local_rhs (stokes_fe.dofs_per_cell)
-    {}
 
+      template <int dim>
+      struct StokesSystem : public StokesPreconditioner<dim>
+      {
+         StokesSystem (const FiniteElement<dim> &stokes_fe);
+         StokesSystem (const StokesSystem<dim> &data);
 
-    template <int dim>
-    StokesSystem<dim>::
-    StokesSystem (const StokesSystem<dim> &data)
-      :
-      StokesPreconditioner<dim> (data),
-      local_rhs (data.local_rhs)
-    {}
+         Vector<double> local_rhs;
+      };
 
 
+      template <int dim>
+      StokesSystem<dim>::
+      StokesSystem (const FiniteElement<dim> &stokes_fe)
+                     :
+                     StokesPreconditioner<dim> (stokes_fe),
+                     local_rhs (stokes_fe.dofs_per_cell)
+      {}
 
-    template <int dim>
-    struct TemperatureMatrix
-    {
-      TemperatureMatrix (const FiniteElement<dim> &temperature_fe);
-      TemperatureMatrix (const TemperatureMatrix &data);
 
-      FullMatrix<double>          local_mass_matrix;
-      FullMatrix<double>          local_stiffness_matrix;
-      std::vector<unsigned int>   local_dof_indices;
-    };
+      template <int dim>
+      StokesSystem<dim>::
+      StokesSystem (const StokesSystem<dim> &data)
+                     :
+                     StokesPreconditioner<dim> (data),
+                     local_rhs (data.local_rhs)
+      {}
 
-    template <int dim>
-    TemperatureMatrix<dim>::
-    TemperatureMatrix (const FiniteElement<dim> &temperature_fe)
-      :
-      local_mass_matrix (temperature_fe.dofs_per_cell,
-                         temperature_fe.dofs_per_cell),
-      local_stiffness_matrix (temperature_fe.dofs_per_cell,
-                              temperature_fe.dofs_per_cell),
-      local_dof_indices (temperature_fe.dofs_per_cell)
-    {}
 
 
-    template <int dim>
-    TemperatureMatrix<dim>::
-    TemperatureMatrix (const TemperatureMatrix &data)
-      :
-      local_mass_matrix (data.local_mass_matrix),
-      local_stiffness_matrix (data.local_stiffness_matrix),
-      local_dof_indices (data.local_dof_indices)
-    {}
+      template <int dim>
+      struct TemperatureMatrix
+      {
+         TemperatureMatrix (const FiniteElement<dim> &temperature_fe);
+         TemperatureMatrix (const TemperatureMatrix &data);
+
+         FullMatrix<double>          local_mass_matrix;
+         FullMatrix<double>          local_stiffness_matrix;
+         std::vector<unsigned int>   local_dof_indices;
+      };
+
+      template <int dim>
+      TemperatureMatrix<dim>::
+      TemperatureMatrix (const FiniteElement<dim> &temperature_fe)
+                     :
+                     local_mass_matrix (temperature_fe.dofs_per_cell,
+                                        temperature_fe.dofs_per_cell),
+                     local_stiffness_matrix (temperature_fe.dofs_per_cell,
+                                             temperature_fe.dofs_per_cell),
+                     local_dof_indices (temperature_fe.dofs_per_cell)
+      {}
 
 
-    template <int dim>
-    struct TemperatureRHS
-    {
-      TemperatureRHS (const FiniteElement<dim> &temperature_fe);
-      TemperatureRHS (const TemperatureRHS &data);
+      template <int dim>
+      TemperatureMatrix<dim>::
+      TemperatureMatrix (const TemperatureMatrix &data)
+                     :
+                     local_mass_matrix (data.local_mass_matrix),
+                     local_stiffness_matrix (data.local_stiffness_matrix),
+                     local_dof_indices (data.local_dof_indices)
+      {}
 
-      Vector<double>              local_rhs;
-      std::vector<unsigned int>   local_dof_indices;
-      FullMatrix<double>          matrix_for_bc;
-    };
 
-    template <int dim>
-    TemperatureRHS<dim>::
-    TemperatureRHS (const FiniteElement<dim> &temperature_fe)
-      :
-      local_rhs (temperature_fe.dofs_per_cell),
-      local_dof_indices (temperature_fe.dofs_per_cell),
-      matrix_for_bc (temperature_fe.dofs_per_cell,
-                     temperature_fe.dofs_per_cell)
-    {}
+      template <int dim>
+      struct TemperatureRHS
+      {
+         TemperatureRHS (const FiniteElement<dim> &temperature_fe);
+         TemperatureRHS (const TemperatureRHS &data);
+
+         Vector<double>              local_rhs;
+         std::vector<unsigned int>   local_dof_indices;
+         FullMatrix<double>          matrix_for_bc;
+      };
+
+      template <int dim>
+      TemperatureRHS<dim>::
+      TemperatureRHS (const FiniteElement<dim> &temperature_fe)
+                     :
+                     local_rhs (temperature_fe.dofs_per_cell),
+                     local_dof_indices (temperature_fe.dofs_per_cell),
+                     matrix_for_bc (temperature_fe.dofs_per_cell,
+                                    temperature_fe.dofs_per_cell)
+      {}
 
 
-    template <int dim>
-    TemperatureRHS<dim>::
-    TemperatureRHS (const TemperatureRHS &data)
-      :
-      local_rhs (data.local_rhs),
-      local_dof_indices (data.local_dof_indices),
-      matrix_for_bc (data.matrix_for_bc)
-    {}
+      template <int dim>
+      TemperatureRHS<dim>::
+      TemperatureRHS (const TemperatureRHS &data)
+                     :
+                     local_rhs (data.local_rhs),
+                     local_dof_indices (data.local_dof_indices),
+                     matrix_for_bc (data.matrix_for_bc)
+      {}
+    }
   }
-}
 
 
 
@@ -886,321 +888,321 @@ namespace Assembly
 // program (this choice is made in the
 // from this variable in the results
 // section of this tutorial program.
-template <int dim>
-class BoussinesqFlowProblem
-{
-  public:
-    struct Parameters;
-    BoussinesqFlowProblem (Parameters &parameters);
-    void run ();
-
-  private:
-    void setup_dofs ();
-    void assemble_stokes_preconditioner ();
-    void build_stokes_preconditioner ();
-    void assemble_stokes_system ();
-    void assemble_temperature_matrix ();
-    void assemble_temperature_system (const double maximal_velocity);
-    void project_temperature_field ();
-    double get_maximal_velocity () const;
-    double get_cfl_number () const;
-    double get_entropy_variation (const double average_temperature) const;
-    std::pair<double,double> get_extrapolated_temperature_range () const;
-    void solve ();
-    void output_results ();
-    void refine_mesh (const unsigned int max_grid_level);
+  template <int dim>
+  class BoussinesqFlowProblem
+  {
+    public:
+      struct Parameters;
+      BoussinesqFlowProblem (Parameters &parameters);
+      void run ();
 
-    double
-    compute_viscosity(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<double>          &old_temperature_laplacians,
-                      const std::vector<double>          &old_old_temperature_laplacians,
-                      const std::vector<Tensor<1,dim> >  &old_velocity_values,
-                      const std::vector<Tensor<1,dim> >  &old_old_velocity_values,
-                      const std::vector<SymmetricTensor<2,dim> >  &old_strain_rates,
-                      const std::vector<SymmetricTensor<2,dim> >  &old_old_strain_rates,
-                      const double                        global_u_infty,
-                      const double                        global_T_variation,
-                      const double                        average_temperature,
-                      const double                        global_entropy_variation,
-                      const double                        cell_diameter) const;
-
-  public:
-    struct Parameters
-    {
-      Parameters (const std::string &parameter_filename);
+    private:
+      void setup_dofs ();
+      void assemble_stokes_preconditioner ();
+      void build_stokes_preconditioner ();
+      void assemble_stokes_system ();
+      void assemble_temperature_matrix ();
+      void assemble_temperature_system (const double maximal_velocity);
+      void project_temperature_field ();
+      double get_maximal_velocity () const;
+      double get_cfl_number () const;
+      double get_entropy_variation (const double average_temperature) const;
+      std::pair<double,double> get_extrapolated_temperature_range () const;
+      void solve ();
+      void output_results ();
+      void refine_mesh (const unsigned int max_grid_level);
+
+      double
+      compute_viscosity(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<double>          &old_temperature_laplacians,
+                       const std::vector<double>          &old_old_temperature_laplacians,
+                       const std::vector<Tensor<1,dim> >  &old_velocity_values,
+                       const std::vector<Tensor<1,dim> >  &old_old_velocity_values,
+                       const std::vector<SymmetricTensor<2,dim> >  &old_strain_rates,
+                       const std::vector<SymmetricTensor<2,dim> >  &old_old_strain_rates,
+                       const double                        global_u_infty,
+                       const double                        global_T_variation,
+                       const double                        average_temperature,
+                       const double                        global_entropy_variation,
+                       const double                        cell_diameter) const;
 
-      static void declare_parameters (ParameterHandler &prm);
-      void parse_parameters (ParameterHandler &prm);
+    public:
+      struct Parameters
+      {
+         Parameters (const std::string &parameter_filename);
 
-      double end_time;
+         static void declare_parameters (ParameterHandler &prm);
+         void parse_parameters (ParameterHandler &prm);
 
-      unsigned int initial_global_refinement;
-      unsigned int initial_adaptive_refinement;
+         double end_time;
 
-      bool         generate_graphical_output;
-      unsigned int graphical_output_interval;
+         unsigned int initial_global_refinement;
+         unsigned int initial_adaptive_refinement;
 
-      unsigned int adaptive_refinement_interval;
+         bool         generate_graphical_output;
+         unsigned int graphical_output_interval;
 
-      double       stabilization_alpha;
-      double       stabilization_c_R;
-      double       stabilization_beta;
+         unsigned int adaptive_refinement_interval;
 
-      unsigned int stokes_velocity_degree;
-      bool         use_locally_conservative_discretization;
+         double       stabilization_alpha;
+         double       stabilization_c_R;
+         double       stabilization_beta;
 
-      unsigned int temperature_degree;
-    };
+         unsigned int stokes_velocity_degree;
+         bool         use_locally_conservative_discretization;
 
-  private:
-    Parameters                           &parameters;
-    ConditionalOStream                  pcout;
+         unsigned int temperature_degree;
+      };
 
-    parallel::distributed::Triangulation<dim> triangulation;
-    double                              global_Omega_diameter;
+    private:
+      Parameters                           &parameters;
+      ConditionalOStream                  pcout;
 
-    const MappingQ<dim>                 mapping;
+      parallel::distributed::Triangulation<dim> triangulation;
+      double                              global_Omega_diameter;
 
-    const FESystem<dim>                 stokes_fe;
+      const MappingQ<dim>                 mapping;
 
-    DoFHandler<dim>                     stokes_dof_handler;
-    ConstraintMatrix                    stokes_constraints;
+      const FESystem<dim>                 stokes_fe;
 
-    TrilinosWrappers::BlockSparseMatrix stokes_matrix;
-    TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix;
+      DoFHandler<dim>                     stokes_dof_handler;
+      ConstraintMatrix                    stokes_constraints;
 
-    TrilinosWrappers::MPI::BlockVector  stokes_solution;
-    TrilinosWrappers::MPI::BlockVector  old_stokes_solution;
-    TrilinosWrappers::MPI::BlockVector  stokes_rhs;
+      TrilinosWrappers::BlockSparseMatrix stokes_matrix;
+      TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix;
 
+      TrilinosWrappers::MPI::BlockVector  stokes_solution;
+      TrilinosWrappers::MPI::BlockVector  old_stokes_solution;
+      TrilinosWrappers::MPI::BlockVector  stokes_rhs;
 
-    FE_Q<dim>                           temperature_fe;
-    DoFHandler<dim>                     temperature_dof_handler;
-    ConstraintMatrix                    temperature_constraints;
 
-    TrilinosWrappers::SparseMatrix      temperature_mass_matrix;
-    TrilinosWrappers::SparseMatrix      temperature_stiffness_matrix;
-    TrilinosWrappers::SparseMatrix      temperature_matrix;
+      FE_Q<dim>                           temperature_fe;
+      DoFHandler<dim>                     temperature_dof_handler;
+      ConstraintMatrix                    temperature_constraints;
 
-    TrilinosWrappers::MPI::Vector       temperature_solution;
-    TrilinosWrappers::MPI::Vector       old_temperature_solution;
-    TrilinosWrappers::MPI::Vector       old_old_temperature_solution;
-    TrilinosWrappers::MPI::Vector       temperature_rhs;
+      TrilinosWrappers::SparseMatrix      temperature_mass_matrix;
+      TrilinosWrappers::SparseMatrix      temperature_stiffness_matrix;
+      TrilinosWrappers::SparseMatrix      temperature_matrix;
 
+      TrilinosWrappers::MPI::Vector       temperature_solution;
+      TrilinosWrappers::MPI::Vector       old_temperature_solution;
+      TrilinosWrappers::MPI::Vector       old_old_temperature_solution;
+      TrilinosWrappers::MPI::Vector       temperature_rhs;
 
-    double time_step;
-    double old_time_step;
-    unsigned int timestep_number;
 
-    std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionAMG> Amg_preconditioner;
-    std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionILU> Mp_preconditioner;
-    std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionIC>  T_preconditioner;
+      double time_step;
+      double old_time_step;
+      unsigned int timestep_number;
 
-    bool rebuild_stokes_matrix;
-    bool rebuild_stokes_preconditioner;
-    bool rebuild_temperature_matrices;
-    bool rebuild_temperature_preconditioner;
+      std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionAMG> Amg_preconditioner;
+      std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionILU> Mp_preconditioner;
+      std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionIC>  T_preconditioner;
 
-    TimerOutput computing_timer;
+      bool rebuild_stokes_matrix;
+      bool rebuild_stokes_preconditioner;
+      bool rebuild_temperature_matrices;
+      bool rebuild_temperature_preconditioner;
 
-    void setup_stokes_matrix (const std::vector<IndexSet> &stokes_partitioning);
-    void setup_stokes_preconditioner (const std::vector<IndexSet> &stokes_partitioning);
-    void setup_temperature_matrices (const IndexSet &temperature_partitioning);
+      TimerOutput computing_timer;
 
-    void
-    local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
-                                          Assembly::Scratch::StokesPreconditioner<dim> &scratch,
-                                          Assembly::CopyData::StokesPreconditioner<dim> &data);
+      void setup_stokes_matrix (const std::vector<IndexSet> &stokes_partitioning);
+      void setup_stokes_preconditioner (const std::vector<IndexSet> &stokes_partitioning);
+      void setup_temperature_matrices (const IndexSet &temperature_partitioning);
 
-    void
-    copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPreconditioner<dim> &data);
+      void
+      local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
+                                           Assembly::Scratch::StokesPreconditioner<dim> &scratch,
+                                           Assembly::CopyData::StokesPreconditioner<dim> &data);
 
+      void
+      copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPreconditioner<dim> &data);
 
-    void
-    local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
-                                  Assembly::Scratch::StokesSystem<dim>  &scratch,
-                                  Assembly::CopyData::StokesSystem<dim> &data);
 
-    void
-    copy_local_to_global_stokes_system (const Assembly::CopyData::StokesSystem<dim> &data);
+      void
+      local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
+                                   Assembly::Scratch::StokesSystem<dim>  &scratch,
+                                   Assembly::CopyData::StokesSystem<dim> &data);
 
+      void
+      copy_local_to_global_stokes_system (const Assembly::CopyData::StokesSystem<dim> &data);
 
-    void
-    local_assemble_temperature_matrix (const typename DoFHandler<dim>::active_cell_iterator &cell,
-                                       Assembly::Scratch::TemperatureMatrix<dim>  &scratch,
-                                       Assembly::CopyData::TemperatureMatrix<dim> &data);
 
-    void
-    copy_local_to_global_temperature_matrix (const Assembly::CopyData::TemperatureMatrix<dim> &data);
+      void
+      local_assemble_temperature_matrix (const typename DoFHandler<dim>::active_cell_iterator &cell,
+                                        Assembly::Scratch::TemperatureMatrix<dim>  &scratch,
+                                        Assembly::CopyData::TemperatureMatrix<dim> &data);
 
+      void
+      copy_local_to_global_temperature_matrix (const Assembly::CopyData::TemperatureMatrix<dim> &data);
 
 
-    void
-    local_assemble_temperature_rhs (const std::pair<double,double> global_T_range,
-                                    const double                   global_max_velocity,
-                                    const double                   global_entropy_variation,
-                                    const typename DoFHandler<dim>::active_cell_iterator &cell,
-                                    Assembly::Scratch::TemperatureRHS<dim> &scratch,
-                                    Assembly::CopyData::TemperatureRHS<dim> &data);
 
-    void
-    copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS<dim> &data);
+      void
+      local_assemble_temperature_rhs (const std::pair<double,double> global_T_range,
+                                     const double                   global_max_velocity,
+                                     const double                   global_entropy_variation,
+                                     const typename DoFHandler<dim>::active_cell_iterator &cell,
+                                     Assembly::Scratch::TemperatureRHS<dim> &scratch,
+                                     Assembly::CopyData::TemperatureRHS<dim> &data);
 
-    class Postprocessor;
-};
+      void
+      copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS<dim> &data);
+
+      class Postprocessor;
+  };
 
 
 // @sect3{BoussinesqFlowProblem class implementation}
 
 // @sect4{BoussinesqFlowProblem::Parameters}
-template <int dim>
-BoussinesqFlowProblem<dim>::Parameters::Parameters (const std::string &parameter_filename)
-  :
-  end_time (1e8),
-  initial_global_refinement (2),
-  initial_adaptive_refinement (2),
-  adaptive_refinement_interval (10),
-  stabilization_alpha (2),
-  stabilization_c_R (0.11),
-  stabilization_beta (0.078),
-  stokes_velocity_degree (2),
-  use_locally_conservative_discretization (true),
-  temperature_degree (2)
-{
-  ParameterHandler prm;
-  BoussinesqFlowProblem<dim>::Parameters::declare_parameters (prm);
+  template <int dim>
+  BoussinesqFlowProblem<dim>::Parameters::Parameters (const std::string &parameter_filename)
+                 :
+                 end_time (1e8),
+                 initial_global_refinement (2),
+                 initial_adaptive_refinement (2),
+                 adaptive_refinement_interval (10),
+                 stabilization_alpha (2),
+                 stabilization_c_R (0.11),
+                 stabilization_beta (0.078),
+                 stokes_velocity_degree (2),
+                 use_locally_conservative_discretization (true),
+                 temperature_degree (2)
+  {
+    ParameterHandler prm;
+    BoussinesqFlowProblem<dim>::Parameters::declare_parameters (prm);
 
-  std::ifstream parameter_file (parameter_filename.c_str());
+    std::ifstream parameter_file (parameter_filename.c_str());
 
-  if (!parameter_file)
-    {
-      parameter_file.close ();
+    if (!parameter_file)
+      {
+       parameter_file.close ();
 
-      std::ostringstream message;
-      message << "Input parameter file <"
-              << parameter_filename << "> not found. Creating a"
-              << std::endl
-              << "template file of the same name."
-              << std::endl;
+       std::ostringstream message;
+       message << "Input parameter file <"
+               << parameter_filename << "> not found. Creating a"
+               << std::endl
+               << "template file of the same name."
+               << std::endl;
 
-      std::ofstream parameter_out (parameter_filename.c_str());
-      prm.print_parameters (parameter_out,
-                            ParameterHandler::Text);
+       std::ofstream parameter_out (parameter_filename.c_str());
+       prm.print_parameters (parameter_out,
+                             ParameterHandler::Text);
 
-      AssertThrow (false, ExcMessage (message.str().c_str()));
-    }
+       AssertThrow (false, ExcMessage (message.str().c_str()));
+      }
 
-  const bool success = prm.read_input (parameter_file);
-  AssertThrow (success, ExcMessage ("Invalid input parameter file."));
+    const bool success = prm.read_input (parameter_file);
+    AssertThrow (success, ExcMessage ("Invalid input parameter file."));
 
-  parse_parameters (prm);
-}
+    parse_parameters (prm);
+  }
 
 
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::Parameters::
-declare_parameters (ParameterHandler &prm)
-{
-  prm.declare_entry ("End time", "1e8",
-                     Patterns::Double (0),
-                     "The end time of the simulation in years.");
-  prm.declare_entry ("Initial global refinement", "2",
-                     Patterns::Integer (0),
-                     "The number of global refinement steps performed on "
-                     "the initial coarse mesh, before the problem is first "
-                     "solved there.");
-  prm.declare_entry ("Initial adaptive refinement", "2",
-                     Patterns::Integer (0),
-                     "The number of adaptive refinement steps performed after "
-                     "initial global refinement.");
-  prm.declare_entry ("Time steps between mesh refinement", "10",
-                     Patterns::Integer (1),
-                     "The number of time steps after which the mesh is to be "
-                     "adapted based on computed error indicators.");
-  prm.declare_entry ("Generate graphical output", "false",
-                     Patterns::Bool (),
-                     "Whether graphical output is to be generated or not. "
-                     "You may not want to get graphical output if the number "
-                     "of processors is large.");
-  prm.declare_entry ("Time steps between graphical output", "50",
-                     Patterns::Integer (1),
-                     "The number of time steps between each generation of "
-                     "graphical output files.");
-
-  prm.enter_subsection ("Stabilization parameters");
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::Parameters::
+  declare_parameters (ParameterHandler &prm)
   {
-    prm.declare_entry ("alpha", "2",
-                       Patterns::Double (1, 2),
-                       "The exponent in the entropy viscosity stabilization.");
-    prm.declare_entry ("c_R", "0.11",
-                       Patterns::Double (0),
-                       "The c_R factor in the entropy viscosity "
-                       "stabilization.");
-    prm.declare_entry ("beta", "0.078",
-                       Patterns::Double (0),
-                       "The beta factor in the artificial viscosity "
-                       "stabilization. An appropriate value for 2d is 0.052 "
-                       "and 0.078 for 3d.");
-  }
-  prm.leave_subsection ();
+    prm.declare_entry ("End time", "1e8",
+                      Patterns::Double (0),
+                      "The end time of the simulation in years.");
+    prm.declare_entry ("Initial global refinement", "2",
+                      Patterns::Integer (0),
+                      "The number of global refinement steps performed on "
+                      "the initial coarse mesh, before the problem is first "
+                      "solved there.");
+    prm.declare_entry ("Initial adaptive refinement", "2",
+                      Patterns::Integer (0),
+                      "The number of adaptive refinement steps performed after "
+                      "initial global refinement.");
+    prm.declare_entry ("Time steps between mesh refinement", "10",
+                      Patterns::Integer (1),
+                      "The number of time steps after which the mesh is to be "
+                      "adapted based on computed error indicators.");
+    prm.declare_entry ("Generate graphical output", "false",
+                      Patterns::Bool (),
+                      "Whether graphical output is to be generated or not. "
+                      "You may not want to get graphical output if the number "
+                      "of processors is large.");
+    prm.declare_entry ("Time steps between graphical output", "50",
+                      Patterns::Integer (1),
+                      "The number of time steps between each generation of "
+                      "graphical output files.");
+
+    prm.enter_subsection ("Stabilization parameters");
+    {
+      prm.declare_entry ("alpha", "2",
+                        Patterns::Double (1, 2),
+                        "The exponent in the entropy viscosity stabilization.");
+      prm.declare_entry ("c_R", "0.11",
+                        Patterns::Double (0),
+                        "The c_R factor in the entropy viscosity "
+                        "stabilization.");
+      prm.declare_entry ("beta", "0.078",
+                        Patterns::Double (0),
+                        "The beta factor in the artificial viscosity "
+                        "stabilization. An appropriate value for 2d is 0.052 "
+                        "and 0.078 for 3d.");
+    }
+    prm.leave_subsection ();
 
-  prm.enter_subsection ("Discretization");
-  {
-    prm.declare_entry ("Stokes velocity polynomial degree", "2",
-                       Patterns::Integer (1),
-                       "The polynomial degree to use for the velocity variables "
-                       "in the Stokes system.");
-    prm.declare_entry ("Temperature polynomial degree", "2",
-                       Patterns::Integer (1),
-                       "The polynomial degree to use for the temperature variable.");
-    prm.declare_entry ("Use locally conservative discretization", "true",
-                       Patterns::Bool (),
-                       "Whether to use a Stokes discretization that is locally "
-                       "conservative at the expense of a larger number of degrees "
-                       "of freedom, or to go with a cheaper discretization "
-                       "that does not locally conserve mass (although it is "
-                       "globally conservative.");
+    prm.enter_subsection ("Discretization");
+    {
+      prm.declare_entry ("Stokes velocity polynomial degree", "2",
+                        Patterns::Integer (1),
+                        "The polynomial degree to use for the velocity variables "
+                        "in the Stokes system.");
+      prm.declare_entry ("Temperature polynomial degree", "2",
+                        Patterns::Integer (1),
+                        "The polynomial degree to use for the temperature variable.");
+      prm.declare_entry ("Use locally conservative discretization", "true",
+                        Patterns::Bool (),
+                        "Whether to use a Stokes discretization that is locally "
+                        "conservative at the expense of a larger number of degrees "
+                        "of freedom, or to go with a cheaper discretization "
+                        "that does not locally conserve mass (although it is "
+                        "globally conservative.");
+    }
+    prm.leave_subsection ();
   }
-  prm.leave_subsection ();
-}
 
 
 
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::Parameters::
-parse_parameters (ParameterHandler &prm)
-{
-  end_time                    = prm.get_double ("End time");
-  initial_global_refinement   = prm.get_integer ("Initial global refinement");
-  initial_adaptive_refinement = prm.get_integer ("Initial adaptive refinement");
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::Parameters::
+  parse_parameters (ParameterHandler &prm)
+  {
+    end_time                    = prm.get_double ("End time");
+    initial_global_refinement   = prm.get_integer ("Initial global refinement");
+    initial_adaptive_refinement = prm.get_integer ("Initial adaptive refinement");
 
-  adaptive_refinement_interval= prm.get_integer ("Time steps between mesh refinement");
+    adaptive_refinement_interval= prm.get_integer ("Time steps between mesh refinement");
 
-  generate_graphical_output   = prm.get_bool ("Generate graphical output");
-  graphical_output_interval   = prm.get_integer ("Time steps between graphical output");
+    generate_graphical_output   = prm.get_bool ("Generate graphical output");
+    graphical_output_interval   = prm.get_integer ("Time steps between graphical output");
 
-  prm.enter_subsection ("Stabilization parameters");
-  {
-    stabilization_alpha = prm.get_double ("alpha");
-    stabilization_c_R   = prm.get_double ("c_R");
-    stabilization_beta  = prm.get_double ("beta");
-  }
-  prm.leave_subsection ();
+    prm.enter_subsection ("Stabilization parameters");
+    {
+      stabilization_alpha = prm.get_double ("alpha");
+      stabilization_c_R   = prm.get_double ("c_R");
+      stabilization_beta  = prm.get_double ("beta");
+    }
+    prm.leave_subsection ();
 
-  prm.enter_subsection ("Discretization");
-  {
-    stokes_velocity_degree = prm.get_integer ("Stokes velocity polynomial degree");
-    temperature_degree     = prm.get_integer ("Temperature polynomial degree");
-    use_locally_conservative_discretization
-      = prm.get_bool ("Use locally conservative discretization");
+    prm.enter_subsection ("Discretization");
+    {
+      stokes_velocity_degree = prm.get_integer ("Stokes velocity polynomial degree");
+      temperature_degree     = prm.get_integer ("Temperature polynomial degree");
+      use_locally_conservative_discretization
+       = prm.get_bool ("Use locally conservative discretization");
+    }
+    prm.leave_subsection ();
   }
-  prm.leave_subsection ();
-}
 
 
 
@@ -1248,49 +1250,49 @@ parse_parameters (ParameterHandler &prm)
 // get a summary table in the end of the
 // program which shows us wallclock times
 // (as opposed to CPU times).
-template <int dim>
-BoussinesqFlowProblem<dim>::BoussinesqFlowProblem (Parameters &parameters_)
-  :
-  parameters (parameters_),
-  pcout (std::cout,
-         (Utilities::System::
-          get_this_mpi_process(MPI_COMM_WORLD)
-          == 0)),
-
-  triangulation (MPI_COMM_WORLD,
-                 typename Triangulation<dim>::MeshSmoothing
-                 (Triangulation<dim>::smoothing_on_refinement |
-                  Triangulation<dim>::smoothing_on_coarsening)),
-
-  mapping (4),
-
-  stokes_fe (FE_Q<dim>(parameters.stokes_velocity_degree),
-             dim,
-             (parameters.use_locally_conservative_discretization
-              ?
-              static_cast<const FiniteElement<dim> &>
-              (FE_DGP<dim>(parameters.stokes_velocity_degree-1))
-              :
-              static_cast<const FiniteElement<dim> &>
-              (FE_Q<dim>(parameters.stokes_velocity_degree-1))),
-             1),
-
-  stokes_dof_handler (triangulation),
-
-  temperature_fe (parameters.temperature_degree),
-  temperature_dof_handler (triangulation),
-
-  time_step (0),
-  old_time_step (0),
-  timestep_number (0),
-  rebuild_stokes_matrix (true),
-  rebuild_stokes_preconditioner (true),
-  rebuild_temperature_matrices (true),
-  rebuild_temperature_preconditioner (true),
-
-  computing_timer (pcout, TimerOutput::summary,
-                   TimerOutput::wall_times)
-{}
+  template <int dim>
+  BoussinesqFlowProblem<dim>::BoussinesqFlowProblem (Parameters &parameters_)
+                 :
+                 parameters (parameters_),
+                 pcout (std::cout,
+                        (Utilities::System::
+                         get_this_mpi_process(MPI_COMM_WORLD)
+                         == 0)),
+
+                 triangulation (MPI_COMM_WORLD,
+                                typename Triangulation<dim>::MeshSmoothing
+                                (Triangulation<dim>::smoothing_on_refinement |
+                                 Triangulation<dim>::smoothing_on_coarsening)),
+
+                 mapping (4),
+
+                 stokes_fe (FE_Q<dim>(parameters.stokes_velocity_degree),
+                            dim,
+                            (parameters.use_locally_conservative_discretization
+                             ?
+                             static_cast<const FiniteElement<dim> &>
+                             (FE_DGP<dim>(parameters.stokes_velocity_degree-1))
+                             :
+                             static_cast<const FiniteElement<dim> &>
+                             (FE_Q<dim>(parameters.stokes_velocity_degree-1))),
+                            1),
+
+                 stokes_dof_handler (triangulation),
+
+                 temperature_fe (parameters.temperature_degree),
+                 temperature_dof_handler (triangulation),
+
+                 time_step (0),
+                 old_time_step (0),
+                 timestep_number (0),
+                 rebuild_stokes_matrix (true),
+                 rebuild_stokes_preconditioner (true),
+                 rebuild_temperature_matrices (true),
+                 rebuild_temperature_preconditioner (true),
+
+                 computing_timer (pcout, TimerOutput::summary,
+                                  TimerOutput::wall_times)
+  {}
 
 
 
@@ -1342,46 +1344,46 @@ BoussinesqFlowProblem<dim>::BoussinesqFlowProblem (Parameters &parameters_)
 // global maximum (output), and an
 // integer value one that says that
 // we only work on one double.
-template <int dim>
-double BoussinesqFlowProblem<dim>::get_maximal_velocity () const
-{
-  const QIterated<dim> quadrature_formula (QTrapez<1>(),
-                                           parameters.stokes_velocity_degree);
-  const unsigned int n_q_points = quadrature_formula.size();
-
-  FEValues<dim> fe_values (mapping, stokes_fe, quadrature_formula, update_values);
-  std::vector<Tensor<1,dim> > velocity_values(n_q_points);
+  template <int dim>
+  double BoussinesqFlowProblem<dim>::get_maximal_velocity () const
+  {
+    const QIterated<dim> quadrature_formula (QTrapez<1>(),
+                                            parameters.stokes_velocity_degree);
+    const unsigned int n_q_points = quadrature_formula.size();
 
-  const FEValuesExtractors::Vector velocities (0);
+    FEValues<dim> fe_values (mapping, stokes_fe, quadrature_formula, update_values);
+    std::vector<Tensor<1,dim> > velocity_values(n_q_points);
 
-  double max_local_velocity = 0;
+    const FEValuesExtractors::Vector velocities (0);
 
-  typename DoFHandler<dim>::active_cell_iterator
-  cell = stokes_dof_handler.begin_active(),
-  endc = stokes_dof_handler.end();
-  for (; cell!=endc; ++cell)
-    if (cell->subdomain_id() ==
-        Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
-      {
-        fe_values.reinit (cell);
-        fe_values[velocities].get_function_values (stokes_solution,
-                                                   velocity_values);
-
-        for (unsigned int q=0; q<n_q_points; ++q)
-          max_local_velocity = std::max (max_local_velocity,
-                                         velocity_values[q].norm());
-      }
+    double max_local_velocity = 0;
 
-  double max_velocity = 0.;
+    typename DoFHandler<dim>::active_cell_iterator
+      cell = stokes_dof_handler.begin_active(),
+      endc = stokes_dof_handler.end();
+    for (; cell!=endc; ++cell)
+      if (cell->subdomain_id() ==
+         Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
+       {
+         fe_values.reinit (cell);
+         fe_values[velocities].get_function_values (stokes_solution,
+                                                    velocity_values);
+
+         for (unsigned int q=0; q<n_q_points; ++q)
+           max_local_velocity = std::max (max_local_velocity,
+                                          velocity_values[q].norm());
+       }
+
+    double max_velocity = 0.;
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
-  MPI_Allreduce (&max_local_velocity, &max_velocity, 1, MPI_DOUBLE,
-                 MPI_MAX, MPI_COMM_WORLD);
+    MPI_Allreduce (&max_local_velocity, &max_velocity, 1, MPI_DOUBLE,
+                  MPI_MAX, MPI_COMM_WORLD);
 #else
-  max_velocity = max_local_velocity;
+    max_velocity = max_local_velocity;
 #endif
 
-  return max_velocity;
-}
+    return max_velocity;
+  }
 
 
 
@@ -1389,130 +1391,130 @@ double BoussinesqFlowProblem<dim>::get_maximal_velocity () const
 // compute the cfl number, i.e., maximal
 // velocity on a cell divided by the cell
 // diameter
-template <int dim>
-double BoussinesqFlowProblem<dim>::get_cfl_number () const
-{
-  const QIterated<dim> quadrature_formula (QTrapez<1>(),
-                                           parameters.stokes_velocity_degree);
-  const unsigned int n_q_points = quadrature_formula.size();
-
-  FEValues<dim> fe_values (mapping, stokes_fe, quadrature_formula, update_values);
-  std::vector<Tensor<1,dim> > velocity_values(n_q_points);
+  template <int dim>
+  double BoussinesqFlowProblem<dim>::get_cfl_number () const
+  {
+    const QIterated<dim> quadrature_formula (QTrapez<1>(),
+                                            parameters.stokes_velocity_degree);
+    const unsigned int n_q_points = quadrature_formula.size();
 
-  const FEValuesExtractors::Vector velocities (0);
+    FEValues<dim> fe_values (mapping, stokes_fe, quadrature_formula, update_values);
+    std::vector<Tensor<1,dim> > velocity_values(n_q_points);
 
-  double max_local_cfl = 0;
+    const FEValuesExtractors::Vector velocities (0);
 
-  typename DoFHandler<dim>::active_cell_iterator
-  cell = stokes_dof_handler.begin_active(),
-  endc = stokes_dof_handler.end();
-  for (; cell!=endc; ++cell)
-    if (cell->subdomain_id() ==
-        Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
-      {
-        fe_values.reinit (cell);
-        fe_values[velocities].get_function_values (stokes_solution,
-                                                   velocity_values);
-
-        double max_local_velocity = 1e-10;
-        for (unsigned int q=0; q<n_q_points; ++q)
-          max_local_velocity = std::max (max_local_velocity,
-                                         velocity_values[q].norm());
-        max_local_cfl = std::max(max_local_cfl,
-                                 max_local_velocity / cell->diameter());
-      }
+    double max_local_cfl = 0;
 
-  double max_cfl_number = 0.;
+    typename DoFHandler<dim>::active_cell_iterator
+      cell = stokes_dof_handler.begin_active(),
+      endc = stokes_dof_handler.end();
+    for (; cell!=endc; ++cell)
+      if (cell->subdomain_id() ==
+         Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
+       {
+         fe_values.reinit (cell);
+         fe_values[velocities].get_function_values (stokes_solution,
+                                                    velocity_values);
+
+         double max_local_velocity = 1e-10;
+         for (unsigned int q=0; q<n_q_points; ++q)
+           max_local_velocity = std::max (max_local_velocity,
+                                          velocity_values[q].norm());
+         max_local_cfl = std::max(max_local_cfl,
+                                  max_local_velocity / cell->diameter());
+       }
+
+    double max_cfl_number = 0.;
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
-  MPI_Allreduce (&max_local_cfl, &max_cfl_number, 1, MPI_DOUBLE,
-                 MPI_MAX, MPI_COMM_WORLD);
+    MPI_Allreduce (&max_local_cfl, &max_cfl_number, 1, MPI_DOUBLE,
+                  MPI_MAX, MPI_COMM_WORLD);
 #else
-  max_cfl_number = max_local_cfl;
+    max_cfl_number = max_local_cfl;
 #endif
 
-  return max_cfl_number;
-}
+    return max_cfl_number;
+  }
 
 
 
-template <int dim>
-double
-BoussinesqFlowProblem<dim>::get_entropy_variation (const double average_temperature) const
-{
-  // only do this if we really need entropy
-  // variation
-  if (parameters.stabilization_alpha != 2)
-    return 1.;
-
-  // record maximal entropy on Gauss quadrature
-  // points
-  const QGauss<dim> quadrature_formula (parameters.temperature_degree+1);
-  const unsigned int n_q_points = quadrature_formula.size();
-
-  FEValues<dim> fe_values (temperature_fe, quadrature_formula,
-                           update_values | update_JxW_values);
-  std::vector<double> old_temperature_values(n_q_points);
-  std::vector<double> old_old_temperature_values(n_q_points);
-
-  double min_entropy = std::numeric_limits<double>::max(),
-         max_entropy = -std::numeric_limits<double>::max(),
-         area = 0,
-         entropy_integrated = 0;
-
-  typename DoFHandler<dim>::active_cell_iterator
-  cell = temperature_dof_handler.begin_active(),
-  endc = temperature_dof_handler.end();
-  for (; cell!=endc; ++cell)
-    if (cell->subdomain_id() ==
-        Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
-      {
-        fe_values.reinit (cell);
-        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)
-          {
-            const double T = (old_temperature_values[q] +
-                              old_old_temperature_values[q]) / 2;
-            const double entropy = ((T-average_temperature) *
-                                    (T-average_temperature));
-
-            min_entropy = std::min (min_entropy, entropy);
-            max_entropy = std::max (max_entropy, entropy);
-            area += fe_values.JxW(q);
-            entropy_integrated += fe_values.JxW(q) * entropy;
-          }
-      }
+  template <int dim>
+  double
+  BoussinesqFlowProblem<dim>::get_entropy_variation (const double average_temperature) const
+  {
+                                    // only do this if we really need entropy
+                                    // variation
+    if (parameters.stabilization_alpha != 2)
+      return 1.;
+
+                                    // record maximal entropy on Gauss quadrature
+                                    // points
+    const QGauss<dim> quadrature_formula (parameters.temperature_degree+1);
+    const unsigned int n_q_points = quadrature_formula.size();
+
+    FEValues<dim> fe_values (temperature_fe, quadrature_formula,
+                            update_values | update_JxW_values);
+    std::vector<double> old_temperature_values(n_q_points);
+    std::vector<double> old_old_temperature_values(n_q_points);
+
+    double min_entropy = std::numeric_limits<double>::max(),
+          max_entropy = -std::numeric_limits<double>::max(),
+                 area = 0,
+    entropy_integrated = 0;
 
-  // do MPI data exchange: we need to sum over
-  // the two integrals (area,
-  // entropy_integrated), and get the extrema
-  // for maximum and minimum. combine
-  // MPI_Allreduce for two values since that is
-  // an expensive operation
-  double local_for_sum[2], global_for_sum[2];
-  double local_for_max[2], global_for_max[2];
-  local_for_sum[0] = entropy_integrated;
-  local_for_sum[1] = area;
-  local_for_max[0] = -min_entropy;
-  local_for_max[1] = max_entropy;
+    typename DoFHandler<dim>::active_cell_iterator
+      cell = temperature_dof_handler.begin_active(),
+      endc = temperature_dof_handler.end();
+    for (; cell!=endc; ++cell)
+      if (cell->subdomain_id() ==
+         Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
+       {
+         fe_values.reinit (cell);
+         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)
+           {
+             const double T = (old_temperature_values[q] +
+                               old_old_temperature_values[q]) / 2;
+             const double entropy = ((T-average_temperature) *
+                                     (T-average_temperature));
+
+             min_entropy = std::min (min_entropy, entropy);
+             max_entropy = std::max (max_entropy, entropy);
+             area += fe_values.JxW(q);
+             entropy_integrated += fe_values.JxW(q) * entropy;
+           }
+       }
+
+                                    // do MPI data exchange: we need to sum over
+                                    // the two integrals (area,
+                                    // entropy_integrated), and get the extrema
+                                    // for maximum and minimum. combine
+                                    // MPI_Allreduce for two values since that is
+                                    // an expensive operation
+    double local_for_sum[2], global_for_sum[2];
+    double local_for_max[2], global_for_max[2];
+    local_for_sum[0] = entropy_integrated;
+    local_for_sum[1] = area;
+    local_for_max[0] = -min_entropy;
+    local_for_max[1] = max_entropy;
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
-  MPI_Allreduce (&local_for_sum[0], &global_for_sum[0], 2, MPI_DOUBLE,
-                 MPI_SUM, MPI_COMM_WORLD);
-  MPI_Allreduce (&local_for_max[0], &global_for_max[0], 2, MPI_DOUBLE,
-                 MPI_MAX, MPI_COMM_WORLD);
+    MPI_Allreduce (&local_for_sum[0], &global_for_sum[0], 2, MPI_DOUBLE,
+                  MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce (&local_for_max[0], &global_for_max[0], 2, MPI_DOUBLE,
+                  MPI_MAX, MPI_COMM_WORLD);
 #else
-  global_for_sum[0] = local_for_sum[0];
-  global_for_sum[1] = local_for_sum[1];
-  global_for_max[0] = local_for_max[0];
-  global_for_max[1] = local_for_max[1];
+    global_for_sum[0] = local_for_sum[0];
+    global_for_sum[1] = local_for_sum[1];
+    global_for_max[0] = local_for_max[0];
+    global_for_max[1] = local_for_max[1];
 #endif
-  const double average_entropy = global_for_sum[0] / global_for_sum[1];
-  const double entropy_diff = std::max(global_for_max[1] - average_entropy,
-                                       average_entropy - (-global_for_max[0]));
-  return entropy_diff;
-}
+    const double average_entropy = global_for_sum[0] / global_for_sum[1];
+    const double entropy_diff = std::max(global_for_max[1] - average_entropy,
+                                        average_entropy - (-global_for_max[0]));
+    return entropy_diff;
+  }
 
 
 
@@ -1525,182 +1527,182 @@ BoussinesqFlowProblem<dim>::get_entropy_variation (const double average_temperat
 // partition. Two global communication
 // steps synchronize the data among the
 // processors.
-template <int dim>
-std::pair<double,double>
-BoussinesqFlowProblem<dim>::get_extrapolated_temperature_range () const
-{
-  const QIterated<dim> quadrature_formula (QTrapez<1>(),
-                                           parameters.temperature_degree);
-  const unsigned int n_q_points = quadrature_formula.size();
-
-  FEValues<dim> fe_values (mapping, temperature_fe, quadrature_formula,
-                           update_values);
-  std::vector<double> old_temperature_values(n_q_points);
-  std::vector<double> old_old_temperature_values(n_q_points);
-
-  // This presets the minimum with a bigger
-  // and the maximum with a smaller number
-  // than one that is going to appear. Will
-  // be overwritten in the cell loop or in
-  // the communication step at the
-  // latest.
-  double min_local_temperature = std::numeric_limits<double>::max(),
-         max_local_temperature = -std::numeric_limits<double>::max();
-
-  if (timestep_number != 0)
-    {
-      typename DoFHandler<dim>::active_cell_iterator
-      cell = temperature_dof_handler.begin_active(),
-      endc = temperature_dof_handler.end();
-      for (; cell!=endc; ++cell)
-        if (cell->subdomain_id() ==
-            Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
-          {
-            fe_values.reinit (cell);
-            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)
-              {
-                const double temperature =
-                  (1. + time_step/old_time_step) * old_temperature_values[q]-
-                  time_step/old_time_step * old_old_temperature_values[q];
-
-                min_local_temperature = std::min (min_local_temperature,
-                                                  temperature);
-                max_local_temperature = std::max (max_local_temperature,
-                                                  temperature);
-              }
-          }
-    }
-  else
-    {
-      typename DoFHandler<dim>::active_cell_iterator
-      cell = temperature_dof_handler.begin_active(),
-      endc = temperature_dof_handler.end();
-      for (; cell!=endc; ++cell)
-        if (cell->subdomain_id() ==
-            Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
-          {
-            fe_values.reinit (cell);
-            fe_values.get_function_values (old_temperature_solution,
-                                           old_temperature_values);
-
-            for (unsigned int q=0; q<n_q_points; ++q)
-              {
-                const double temperature = old_temperature_values[q];
-
-                min_local_temperature = std::min (min_local_temperature,
-                                                  temperature);
-                max_local_temperature = std::max (max_local_temperature,
-                                                  temperature);
-              }
-          }
-    }
+  template <int dim>
+  std::pair<double,double>
+  BoussinesqFlowProblem<dim>::get_extrapolated_temperature_range () const
+  {
+    const QIterated<dim> quadrature_formula (QTrapez<1>(),
+                                            parameters.temperature_degree);
+    const unsigned int n_q_points = quadrature_formula.size();
+
+    FEValues<dim> fe_values (mapping, temperature_fe, quadrature_formula,
+                            update_values);
+    std::vector<double> old_temperature_values(n_q_points);
+    std::vector<double> old_old_temperature_values(n_q_points);
+
+                                    // This presets the minimum with a bigger
+                                    // and the maximum with a smaller number
+                                    // than one that is going to appear. Will
+                                    // be overwritten in the cell loop or in
+                                    // the communication step at the
+                                    // latest.
+    double min_local_temperature = std::numeric_limits<double>::max(),
+          max_local_temperature = -std::numeric_limits<double>::max();
+
+    if (timestep_number != 0)
+      {
+       typename DoFHandler<dim>::active_cell_iterator
+         cell = temperature_dof_handler.begin_active(),
+         endc = temperature_dof_handler.end();
+       for (; cell!=endc; ++cell)
+         if (cell->subdomain_id() ==
+             Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
+           {
+             fe_values.reinit (cell);
+             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)
+               {
+                 const double temperature =
+                   (1. + time_step/old_time_step) * old_temperature_values[q]-
+                   time_step/old_time_step * old_old_temperature_values[q];
+
+                 min_local_temperature = std::min (min_local_temperature,
+                                                   temperature);
+                 max_local_temperature = std::max (max_local_temperature,
+                                                   temperature);
+               }
+           }
+      }
+    else
+      {
+       typename DoFHandler<dim>::active_cell_iterator
+         cell = temperature_dof_handler.begin_active(),
+         endc = temperature_dof_handler.end();
+       for (; cell!=endc; ++cell)
+         if (cell->subdomain_id() ==
+             Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
+           {
+             fe_values.reinit (cell);
+             fe_values.get_function_values (old_temperature_solution,
+                                            old_temperature_values);
+
+             for (unsigned int q=0; q<n_q_points; ++q)
+               {
+                 const double temperature = old_temperature_values[q];
+
+                 min_local_temperature = std::min (min_local_temperature,
+                                                   temperature);
+                 max_local_temperature = std::max (max_local_temperature,
+                                                   temperature);
+               }
+           }
+      }
 
-  double min_temperature, max_temperature;
+    double min_temperature, max_temperature;
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
-  MPI_Allreduce (&max_local_temperature, &max_temperature, 1, MPI_DOUBLE,
-                 MPI_MAX, MPI_COMM_WORLD);
-  MPI_Allreduce (&min_local_temperature, &min_temperature, 1, MPI_DOUBLE,
-                 MPI_MIN, MPI_COMM_WORLD);
+    MPI_Allreduce (&max_local_temperature, &max_temperature, 1, MPI_DOUBLE,
+                  MPI_MAX, MPI_COMM_WORLD);
+    MPI_Allreduce (&min_local_temperature, &min_temperature, 1, MPI_DOUBLE,
+                  MPI_MIN, MPI_COMM_WORLD);
 #else
-  min_temperature = min_local_temperature;
-  max_temperature = max_local_temperature;
+    min_temperature = min_local_temperature;
+    max_temperature = max_local_temperature;
 #endif
 
-  return std::make_pair(min_temperature, max_temperature);
-}
+    return std::make_pair(min_temperature, max_temperature);
+  }
 
 
 
 // The function that calculates the
 // viscosity is purely local, so this is
 // the same code as in step-31.
-template <int dim>
-double
-BoussinesqFlowProblem<dim>::
-compute_viscosity (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<double>          &old_temperature_laplacians,
-                   const std::vector<double>          &old_old_temperature_laplacians,
-                   const std::vector<Tensor<1,dim> >  &old_velocity_values,
-                   const std::vector<Tensor<1,dim> >  &old_old_velocity_values,
-                   const std::vector<SymmetricTensor<2,dim> >  &old_strain_rates,
-                   const std::vector<SymmetricTensor<2,dim> >  &old_old_strain_rates,
-                   const double                        global_u_infty,
-                   const double                        global_T_variation,
-                   const double                        average_temperature,
-                   const double                        global_entropy_variation,
-                   const double                        cell_diameter) const
-{
-  if (global_u_infty == 0)
-    return 5e-3 * cell_diameter;
+  template <int dim>
+  double
+  BoussinesqFlowProblem<dim>::
+  compute_viscosity (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<double>          &old_temperature_laplacians,
+                    const std::vector<double>          &old_old_temperature_laplacians,
+                    const std::vector<Tensor<1,dim> >  &old_velocity_values,
+                    const std::vector<Tensor<1,dim> >  &old_old_velocity_values,
+                    const std::vector<SymmetricTensor<2,dim> >  &old_strain_rates,
+                    const std::vector<SymmetricTensor<2,dim> >  &old_old_strain_rates,
+                    const double                        global_u_infty,
+                    const double                        global_T_variation,
+                    const double                        average_temperature,
+                    const double                        global_entropy_variation,
+                    const double                        cell_diameter) const
+  {
+    if (global_u_infty == 0)
+      return 5e-3 * cell_diameter;
 
-  const unsigned int n_q_points = old_temperature.size();
+    const unsigned int n_q_points = old_temperature.size();
 
-  double max_residual = 0;
-  double max_velocity = 0;
+    double max_residual = 0;
+    double max_velocity = 0;
 
-  for (unsigned int q=0; q < n_q_points; ++q)
-    {
-      const Tensor<1,dim> u = (old_velocity_values[q] +
-                               old_old_velocity_values[q]) / 2;
-
-      const SymmetricTensor<2,dim> strain_rate = (old_strain_rates[q] +
-                                                  old_old_strain_rates[q]) / 2;
-
-      const double T = (old_temperature[q] + old_old_temperature[q]) / 2;
-      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 = EquationData::kappa
-                                   * (old_temperature_laplacians[q] +
-                                      old_old_temperature_laplacians[q]) / 2;
-      const double gamma
-        = ((EquationData::radiogenic_heating * EquationData::density(T)
-            +
-            2 * EquationData::eta * strain_rate * strain_rate) /
-           (EquationData::density(T) * EquationData::specific_heat));
-
-      double residual
-        = std::abs(dT_dt + u_grad_T - kappa_Delta_T - gamma);
-      if (parameters.stabilization_alpha == 2)
-        residual *= std::abs(T - average_temperature);
-
-      max_residual = std::max (residual,        max_residual);
-      max_velocity = std::max (std::sqrt (u*u), max_velocity);
-    }
-
-  const double max_viscosity = (parameters.stabilization_beta *
-                                max_velocity * cell_diameter);
-  if (timestep_number == 0)
-    return max_viscosity;
-  else
-    {
-      Assert (old_time_step > 0, ExcInternalError());
-
-      double entropy_viscosity;
-      if (parameters.stabilization_alpha == 2)
-        entropy_viscosity = (parameters.stabilization_c_R *
-                             cell_diameter * cell_diameter *
-                             max_residual /
-                             global_entropy_variation);
-      else
-        entropy_viscosity = (parameters.stabilization_c_R *
-                             cell_diameter * global_Omega_diameter *
-                             max_velocity * max_residual /
-                             (global_u_infty * global_T_variation));
+    for (unsigned int q=0; q < n_q_points; ++q)
+      {
+       const Tensor<1,dim> u = (old_velocity_values[q] +
+                                old_old_velocity_values[q]) / 2;
+
+       const SymmetricTensor<2,dim> strain_rate = (old_strain_rates[q] +
+                                                   old_old_strain_rates[q]) / 2;
+
+       const double T = (old_temperature[q] + old_old_temperature[q]) / 2;
+       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 = EquationData::kappa
+                                    * (old_temperature_laplacians[q] +
+                                       old_old_temperature_laplacians[q]) / 2;
+       const double gamma
+         = ((EquationData::radiogenic_heating * EquationData::density(T)
+             +
+             2 * EquationData::eta * strain_rate * strain_rate) /
+            (EquationData::density(T) * EquationData::specific_heat));
+
+       double residual
+         = std::abs(dT_dt + u_grad_T - kappa_Delta_T - gamma);
+       if (parameters.stabilization_alpha == 2)
+         residual *= std::abs(T - average_temperature);
+
+       max_residual = std::max (residual,        max_residual);
+       max_velocity = std::max (std::sqrt (u*u), max_velocity);
+      }
 
-      return std::min (max_viscosity, entropy_viscosity);
-    }
-}
+    const double max_viscosity = (parameters.stabilization_beta *
+                                 max_velocity * cell_diameter);
+    if (timestep_number == 0)
+      return max_viscosity;
+    else
+      {
+       Assert (old_time_step > 0, ExcInternalError());
+
+       double entropy_viscosity;
+       if (parameters.stabilization_alpha == 2)
+         entropy_viscosity = (parameters.stabilization_c_R *
+                              cell_diameter * cell_diameter *
+                              max_residual /
+                              global_entropy_variation);
+       else
+         entropy_viscosity = (parameters.stabilization_c_R *
+                              cell_diameter * global_Omega_diameter *
+                              max_velocity * max_residual /
+                              (global_u_infty * global_T_variation));
+
+       return std::min (max_viscosity, entropy_viscosity);
+      }
+  }
 
 
 
@@ -1788,116 +1790,116 @@ compute_viscosity (const std::vector<double>          &old_temperature,
 // exactly the entry that we would fill the
 // matrix with &mdash; otherwise, this will
 // not work.
-template <int dim>
-void BoussinesqFlowProblem<dim>::project_temperature_field ()
-{
-  assemble_temperature_matrix ();
+  template <int dim>
+  void BoussinesqFlowProblem<dim>::project_temperature_field ()
+  {
+    assemble_temperature_matrix ();
 
-  QGauss<dim> quadrature(parameters.temperature_degree+2);
-  UpdateFlags update_flags = UpdateFlags(update_values   |
-                                         update_quadrature_points |
-                                         update_JxW_values);
-  FEValues<dim> fe_values (mapping, temperature_fe, quadrature, update_flags);
+    QGauss<dim> quadrature(parameters.temperature_degree+2);
+    UpdateFlags update_flags = UpdateFlags(update_values   |
+                                          update_quadrature_points |
+                                          update_JxW_values);
+    FEValues<dim> fe_values (mapping, temperature_fe, quadrature, update_flags);
 
-  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
-                     n_q_points    = fe_values.n_quadrature_points;
+    const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+                      n_q_points    = fe_values.n_quadrature_points;
 
-  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-  Vector<double> cell_vector (dofs_per_cell);
-  FullMatrix<double> matrix_for_bc (dofs_per_cell, dofs_per_cell);
+    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+    Vector<double> cell_vector (dofs_per_cell);
+    FullMatrix<double> matrix_for_bc (dofs_per_cell, dofs_per_cell);
 
-  std::vector<double> rhs_values(n_q_points);
+    std::vector<double> rhs_values(n_q_points);
 
-  TrilinosWrappers::MPI::Vector
-  rhs (temperature_mass_matrix.row_partitioner()),
+    TrilinosWrappers::MPI::Vector
+      rhs (temperature_mass_matrix.row_partitioner()),
       solution (temperature_mass_matrix.row_partitioner());
 
-  typename DoFHandler<dim>::active_cell_iterator
-  cell = temperature_dof_handler.begin_active(),
-  endc = temperature_dof_handler.end();
-
-  for (; cell!=endc; ++cell)
-    if (cell->subdomain_id() ==
-        Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
-      {
-        cell->get_dof_indices (local_dof_indices);
-        fe_values.reinit(cell);
-
-        EquationData::TemperatureInitialValues<dim>().value_list
-        (fe_values.get_quadrature_points(), rhs_values);
-
-        cell_vector = 0;
-        matrix_for_bc = 0;
-        for (unsigned int point=0; point<n_q_points; ++point)
-          for (unsigned int i=0; i<dofs_per_cell; ++i)
-            {
-              cell_vector(i) += rhs_values[point] *
-                                fe_values.shape_value(i,point) *
-                                fe_values.JxW(point);
-              if (temperature_constraints.is_inhomogeneously_constrained(local_dof_indices[i]))
-                {
-                  for (unsigned int j=0; j<dofs_per_cell; ++j)
-                    matrix_for_bc(j,i) += fe_values.shape_value(i,point) *
-                                          fe_values.shape_value(j,point) *
-                                          fe_values.JxW(point);
-                }
-            }
-
-        temperature_constraints.distribute_local_to_global (cell_vector,
-                                                            local_dof_indices,
-                                                            rhs,
-                                                            matrix_for_bc);
-      }
+    typename DoFHandler<dim>::active_cell_iterator
+      cell = temperature_dof_handler.begin_active(),
+      endc = temperature_dof_handler.end();
 
-  rhs.compress (Add);
-
-  SolverControl solver_control(5*rhs.size(), 1e-12*rhs.l2_norm());
-  SolverCG<TrilinosWrappers::MPI::Vector> cg(solver_control);
-
-  TrilinosWrappers::PreconditionJacobi preconditioner_mass;
-  preconditioner_mass.initialize(temperature_mass_matrix, 1.3);
-
-  cg.solve (temperature_mass_matrix, solution, rhs, preconditioner_mass);
-
-  temperature_constraints.distribute (solution);
-                                  // Having so computed the current
-                                  // temperature field, let us set
-                                  // the member variable that holds
-                                  // the temperature nodes. Strictly
-                                  // speaking, we really only need to
-                                  // set
-                                  // <code>old_temperature_solution</code>
-                                  // since the first thing we will do
-                                  // is to compute the Stokes
-                                  // solution that only requires the
-                                  // previous time step's temperature
-                                  // field. That said, nothing good
-                                  // can come from not initializing
-                                  // the other vectors as well
-                                  // (especially since it's a
-                                  // relatively cheap operation and
-                                  // we only have to do it once at
-                                  // the beginning of the program) if
-                                  // we ever want to extend our
-                                  // numerical method or physical
-                                  // model, and so we initialize
-                                  // <code>temperature_solution</code>
-                                  // and
-                                  // <code>old_old_temperature_solution</code>
-                                  // as well. The second and third
-                                  // argument to the reinit function
-                                  // indicates that the elements of
-                                  // the first argument shall be
-                                  // copied and that while the
-                                  // <code>solution</code> vector is
-                                  // strictly distributed, the
-                                  // vectors here initialized can
-                                  // (and indeed do) have ghost
-                                  // elements:
-  temperature_solution.reinit(solution, false, true);
-  old_temperature_solution.reinit(solution, false, true);
-  old_old_temperature_solution.reinit(solution, false, true);
-}
+    for (; cell!=endc; ++cell)
+      if (cell->subdomain_id() ==
+         Utilities::System::get_this_mpi_process(MPI_COMM_WORLD))
+       {
+         cell->get_dof_indices (local_dof_indices);
+         fe_values.reinit(cell);
+
+         EquationData::TemperatureInitialValues<dim>().value_list
+           (fe_values.get_quadrature_points(), rhs_values);
+
+         cell_vector = 0;
+         matrix_for_bc = 0;
+         for (unsigned int point=0; point<n_q_points; ++point)
+           for (unsigned int i=0; i<dofs_per_cell; ++i)
+             {
+               cell_vector(i) += rhs_values[point] *
+                                 fe_values.shape_value(i,point) *
+                                 fe_values.JxW(point);
+               if (temperature_constraints.is_inhomogeneously_constrained(local_dof_indices[i]))
+                 {
+                   for (unsigned int j=0; j<dofs_per_cell; ++j)
+                     matrix_for_bc(j,i) += fe_values.shape_value(i,point) *
+                                           fe_values.shape_value(j,point) *
+                                           fe_values.JxW(point);
+                 }
+             }
+
+         temperature_constraints.distribute_local_to_global (cell_vector,
+                                                             local_dof_indices,
+                                                             rhs,
+                                                             matrix_for_bc);
+       }
+
+    rhs.compress (Add);
+
+    SolverControl solver_control(5*rhs.size(), 1e-12*rhs.l2_norm());
+    SolverCG<TrilinosWrappers::MPI::Vector> cg(solver_control);
+
+    TrilinosWrappers::PreconditionJacobi preconditioner_mass;
+    preconditioner_mass.initialize(temperature_mass_matrix, 1.3);
+
+    cg.solve (temperature_mass_matrix, solution, rhs, preconditioner_mass);
+
+    temperature_constraints.distribute (solution);
+                                    // Having so computed the current
+                                    // temperature field, let us set
+                                    // the member variable that holds
+                                    // the temperature nodes. Strictly
+                                    // speaking, we really only need to
+                                    // set
+                                    // <code>old_temperature_solution</code>
+                                    // since the first thing we will do
+                                    // is to compute the Stokes
+                                    // solution that only requires the
+                                    // previous time step's temperature
+                                    // field. That said, nothing good
+                                    // can come from not initializing
+                                    // the other vectors as well
+                                    // (especially since it's a
+                                    // relatively cheap operation and
+                                    // we only have to do it once at
+                                    // the beginning of the program) if
+                                    // we ever want to extend our
+                                    // numerical method or physical
+                                    // model, and so we initialize
+                                    // <code>temperature_solution</code>
+                                    // and
+                                    // <code>old_old_temperature_solution</code>
+                                    // as well. The second and third
+                                    // argument to the reinit function
+                                    // indicates that the elements of
+                                    // the first argument shall be
+                                    // copied and that while the
+                                    // <code>solution</code> vector is
+                                    // strictly distributed, the
+                                    // vectors here initialized can
+                                    // (and indeed do) have ghost
+                                    // elements:
+    temperature_solution.reinit(solution, false, true);
+    old_temperature_solution.reinit(solution, false, true);
+    old_old_temperature_solution.reinit(solution, false, true);
+  }
 
 
 
@@ -1979,88 +1981,88 @@ void BoussinesqFlowProblem<dim>::project_temperature_field ()
 // patterns, so we can safely release the
 // <code>sp</code> variable once the matrix
 // has been given the sparsity structure.
-template <int dim>
-void BoussinesqFlowProblem<dim>::
-setup_stokes_matrix (const std::vector<IndexSet> &stokes_partitioning)
-{
-  stokes_matrix.clear ();
-
-  TrilinosWrappers::BlockSparsityPattern sp (stokes_partitioning,
-                                             MPI_COMM_WORLD);
-
-  Table<2,DoFTools::Coupling> coupling (dim+1, dim+1);
-
-  for (unsigned int c=0; c<dim+1; ++c)
-    for (unsigned int d=0; d<dim+1; ++d)
-      if (! ((c==dim) && (d==dim)))
-        coupling[c][d] = DoFTools::always;
-      else
-        coupling[c][d] = DoFTools::none;
+  template <int dim>
+  void BoussinesqFlowProblem<dim>::
+  setup_stokes_matrix (const std::vector<IndexSet> &stokes_partitioning)
+  {
+    stokes_matrix.clear ();
 
-  DoFTools::make_sparsity_pattern (stokes_dof_handler,
-                                   coupling, sp,
-                                   stokes_constraints, false,
-                                   Utilities::System::
-                                   get_this_mpi_process(MPI_COMM_WORLD));
-  sp.compress();
+    TrilinosWrappers::BlockSparsityPattern sp (stokes_partitioning,
+                                              MPI_COMM_WORLD);
 
-  stokes_matrix.reinit (sp);
-}
+    Table<2,DoFTools::Coupling> coupling (dim+1, dim+1);
 
+    for (unsigned int c=0; c<dim+1; ++c)
+      for (unsigned int d=0; d<dim+1; ++d)
+       if (! ((c==dim) && (d==dim)))
+         coupling[c][d] = DoFTools::always;
+       else
+         coupling[c][d] = DoFTools::none;
 
+    DoFTools::make_sparsity_pattern (stokes_dof_handler,
+                                    coupling, sp,
+                                    stokes_constraints, false,
+                                    Utilities::System::
+                                    get_this_mpi_process(MPI_COMM_WORLD));
+    sp.compress();
 
-template <int dim>
-void BoussinesqFlowProblem<dim>::
-setup_stokes_preconditioner (const std::vector<IndexSet> &stokes_partitioning)
-{
-  Amg_preconditioner.reset ();
-  Mp_preconditioner.reset ();
-
-  stokes_preconditioner_matrix.clear ();
-
-  TrilinosWrappers::BlockSparsityPattern sp (stokes_partitioning,
-                                             MPI_COMM_WORLD);
+    stokes_matrix.reinit (sp);
+  }
 
-  Table<2,DoFTools::Coupling> coupling (dim+1, dim+1);
-  for (unsigned int c=0; c<dim+1; ++c)
-    for (unsigned int d=0; d<dim+1; ++d)
-      if (c == d)
-        coupling[c][d] = DoFTools::always;
-      else
-        coupling[c][d] = DoFTools::none;
 
-  DoFTools::make_sparsity_pattern (stokes_dof_handler,
-                                   coupling, sp,
-                                   stokes_constraints, false,
-                                   Utilities::System::
-                                   get_this_mpi_process(MPI_COMM_WORLD));
-  sp.compress();
 
-  stokes_preconditioner_matrix.reinit (sp);
-}
+  template <int dim>
+  void BoussinesqFlowProblem<dim>::
+  setup_stokes_preconditioner (const std::vector<IndexSet> &stokes_partitioning)
+  {
+    Amg_preconditioner.reset ();
+    Mp_preconditioner.reset ();
+
+    stokes_preconditioner_matrix.clear ();
+
+    TrilinosWrappers::BlockSparsityPattern sp (stokes_partitioning,
+                                              MPI_COMM_WORLD);
+
+    Table<2,DoFTools::Coupling> coupling (dim+1, dim+1);
+    for (unsigned int c=0; c<dim+1; ++c)
+      for (unsigned int d=0; d<dim+1; ++d)
+       if (c == d)
+         coupling[c][d] = DoFTools::always;
+       else
+         coupling[c][d] = DoFTools::none;
+
+    DoFTools::make_sparsity_pattern (stokes_dof_handler,
+                                    coupling, sp,
+                                    stokes_constraints, false,
+                                    Utilities::System::
+                                    get_this_mpi_process(MPI_COMM_WORLD));
+    sp.compress();
+
+    stokes_preconditioner_matrix.reinit (sp);
+  }
 
 
-template <int dim>
-void BoussinesqFlowProblem<dim>::
-setup_temperature_matrices (const IndexSet &temperature_partitioner)
-{
-  T_preconditioner.reset ();
-  temperature_mass_matrix.clear ();
-  temperature_stiffness_matrix.clear ();
-  temperature_matrix.clear ();
-
-  TrilinosWrappers::SparsityPattern sp (temperature_partitioner,
-                                        MPI_COMM_WORLD);
-  DoFTools::make_sparsity_pattern (temperature_dof_handler, sp,
-                                   temperature_constraints, false,
-                                   Utilities::System::
-                                   get_this_mpi_process(MPI_COMM_WORLD));
-  sp.compress();
-
-  temperature_matrix.reinit (sp);
-  temperature_mass_matrix.reinit (sp);
-  temperature_stiffness_matrix.reinit (sp);
-}
+  template <int dim>
+  void BoussinesqFlowProblem<dim>::
+  setup_temperature_matrices (const IndexSet &temperature_partitioner)
+  {
+    T_preconditioner.reset ();
+    temperature_mass_matrix.clear ();
+    temperature_stiffness_matrix.clear ();
+    temperature_matrix.clear ();
+
+    TrilinosWrappers::SparsityPattern sp (temperature_partitioner,
+                                         MPI_COMM_WORLD);
+    DoFTools::make_sparsity_pattern (temperature_dof_handler, sp,
+                                    temperature_constraints, false,
+                                    Utilities::System::
+                                    get_this_mpi_process(MPI_COMM_WORLD));
+    sp.compress();
+
+    temperature_matrix.reinit (sp);
+    temperature_mass_matrix.reinit (sp);
+    temperature_stiffness_matrix.reinit (sp);
+  }
 
 
 
@@ -2142,130 +2144,130 @@ setup_temperature_matrices (const IndexSet &temperature_partitioner)
 // get_this_mpi_process
 // (MPI_COMM_WORLD) == 0)</code>,
 // hardly a pretty solution.
-template <int dim>
-void BoussinesqFlowProblem<dim>::setup_dofs ()
-{
-  computing_timer.enter_section("Setup dof systems");
-
-  std::vector<unsigned int> stokes_sub_blocks (dim+1,0);
-  stokes_sub_blocks[dim] = 1;
-  stokes_dof_handler.distribute_dofs (stokes_fe);
-  DoFRenumbering::component_wise (stokes_dof_handler, stokes_sub_blocks);
-
-  temperature_dof_handler.distribute_dofs (temperature_fe);
-
-  std::vector<unsigned int> stokes_dofs_per_block (2);
-  DoFTools::count_dofs_per_block (stokes_dof_handler, stokes_dofs_per_block,
-                                  stokes_sub_blocks);
-
-  const unsigned int n_u = stokes_dofs_per_block[0],
-                     n_p = stokes_dofs_per_block[1],
-                     n_T = temperature_dof_handler.n_dofs();
-
-  // print dof numbers with 1000s
-  // separator since they are frequently
-  // large
-  std::locale s = pcout.get_stream().getloc();
-  pcout.get_stream().imbue(std::locale(""));
-  pcout << "Number of active cells: "
-        << triangulation.n_global_active_cells()
-        << " (on "
-        << triangulation.n_levels()
-        << " levels)"
-        << std::endl
-        << "Number of degrees of freedom: "
-        << n_u + n_p + n_T
-        << " (" << n_u << '+' << n_p << '+'<< n_T <<')'
-        << std::endl
-        << std::endl;
-  pcout.get_stream().imbue(s);
-
-
-
-  std::vector<IndexSet> stokes_partitioning, stokes_relevant_partitioning;
-  IndexSet temperature_partitioning (n_T), temperature_relevant_partitioning (n_T);
-  IndexSet stokes_relevant_set;
+  template <int dim>
+  void BoussinesqFlowProblem<dim>::setup_dofs ()
   {
-    const unsigned int my_id =
-      Utilities::System::get_this_mpi_process(MPI_COMM_WORLD);
-    IndexSet stokes_index_set = stokes_dof_handler.locally_owned_dofs();
-    stokes_partitioning.push_back(stokes_index_set.get_view(0,n_u));
-    stokes_partitioning.push_back(stokes_index_set.get_view(n_u,n_u+n_p));
-
-    DoFTools::extract_locally_relevant_dofs (stokes_dof_handler,
-                                             stokes_relevant_set);
-    stokes_relevant_partitioning.push_back(stokes_relevant_set.get_view(0,n_u));
-    stokes_relevant_partitioning.push_back(stokes_relevant_set.get_view(n_u,n_u+n_p));
-
-    temperature_partitioning = temperature_dof_handler.locally_owned_dofs();
-    DoFTools::extract_locally_relevant_dofs (temperature_dof_handler,
-                                             temperature_relevant_partitioning);
-  }
+    computing_timer.enter_section("Setup dof systems");
+
+    std::vector<unsigned int> stokes_sub_blocks (dim+1,0);
+    stokes_sub_blocks[dim] = 1;
+    stokes_dof_handler.distribute_dofs (stokes_fe);
+    DoFRenumbering::component_wise (stokes_dof_handler, stokes_sub_blocks);
+
+    temperature_dof_handler.distribute_dofs (temperature_fe);
+
+    std::vector<unsigned int> stokes_dofs_per_block (2);
+    DoFTools::count_dofs_per_block (stokes_dof_handler, stokes_dofs_per_block,
+                                   stokes_sub_blocks);
+
+    const unsigned int n_u = stokes_dofs_per_block[0],
+                      n_p = stokes_dofs_per_block[1],
+                      n_T = temperature_dof_handler.n_dofs();
+
+                                    // print dof numbers with 1000s
+                                    // separator since they are frequently
+                                    // large
+    std::locale s = pcout.get_stream().getloc();
+    pcout.get_stream().imbue(std::locale(""));
+    pcout << "Number of active cells: "
+         << triangulation.n_global_active_cells()
+         << " (on "
+         << triangulation.n_levels()
+         << " levels)"
+         << std::endl
+         << "Number of degrees of freedom: "
+         << n_u + n_p + n_T
+         << " (" << n_u << '+' << n_p << '+'<< n_T <<')'
+         << std::endl
+         << std::endl;
+    pcout.get_stream().imbue(s);
+
+
+
+    std::vector<IndexSet> stokes_partitioning, stokes_relevant_partitioning;
+    IndexSet temperature_partitioning (n_T), temperature_relevant_partitioning (n_T);
+    IndexSet stokes_relevant_set;
+    {
+      const unsigned int my_id =
+       Utilities::System::get_this_mpi_process(MPI_COMM_WORLD);
+      IndexSet stokes_index_set = stokes_dof_handler.locally_owned_dofs();
+      stokes_partitioning.push_back(stokes_index_set.get_view(0,n_u));
+      stokes_partitioning.push_back(stokes_index_set.get_view(n_u,n_u+n_p));
+
+      DoFTools::extract_locally_relevant_dofs (stokes_dof_handler,
+                                              stokes_relevant_set);
+      stokes_relevant_partitioning.push_back(stokes_relevant_set.get_view(0,n_u));
+      stokes_relevant_partitioning.push_back(stokes_relevant_set.get_view(n_u,n_u+n_p));
+
+      temperature_partitioning = temperature_dof_handler.locally_owned_dofs();
+      DoFTools::extract_locally_relevant_dofs (temperature_dof_handler,
+                                              temperature_relevant_partitioning);
+    }
 
-  {
+    {
 
-    stokes_constraints.clear ();
+      stokes_constraints.clear ();
 //    IndexSet stokes_la;
 //    DoFTools::extract_locally_active_dofs (stokes_dof_handler,
 //             stokes_la);
-    stokes_constraints.reinit(stokes_relevant_set);
-
-    DoFTools::make_hanging_node_constraints (stokes_dof_handler,
-                                             stokes_constraints);
-
-    std::vector<bool> velocity_mask (dim+1, true);
-    velocity_mask[dim] = false;
-    VectorTools::interpolate_boundary_values (stokes_dof_handler,
-                                              0,
-                                              ZeroFunction<dim>(dim+1),
-                                              stokes_constraints,
-                                              velocity_mask);
-
-    std::set<unsigned char> no_normal_flux_boundaries;
-    no_normal_flux_boundaries.insert (1);
-    VectorTools::compute_no_normal_flux_constraints (stokes_dof_handler, 0,
-                                                     no_normal_flux_boundaries,
-                                                     stokes_constraints,
-                                                     mapping);
-    stokes_constraints.close ();
-  }
-  {
-    temperature_constraints.clear ();
-    temperature_constraints.reinit(temperature_relevant_partitioning);//temp_locally_active);
-
-    DoFTools::make_hanging_node_constraints (temperature_dof_handler,
-                                             temperature_constraints);
-    VectorTools::interpolate_boundary_values (temperature_dof_handler,
-                                              0,
-                                              EquationData::TemperatureInitialValues<dim>(),
-                                              temperature_constraints);
-    VectorTools::interpolate_boundary_values (temperature_dof_handler,
-                                              1,
-                                              EquationData::TemperatureInitialValues<dim>(),
-                                              temperature_constraints);
-    temperature_constraints.close ();
-  }
+      stokes_constraints.reinit(stokes_relevant_set);
+
+      DoFTools::make_hanging_node_constraints (stokes_dof_handler,
+                                              stokes_constraints);
+
+      std::vector<bool> velocity_mask (dim+1, true);
+      velocity_mask[dim] = false;
+      VectorTools::interpolate_boundary_values (stokes_dof_handler,
+                                               0,
+                                               ZeroFunction<dim>(dim+1),
+                                               stokes_constraints,
+                                               velocity_mask);
+
+      std::set<unsigned char> no_normal_flux_boundaries;
+      no_normal_flux_boundaries.insert (1);
+      VectorTools::compute_no_normal_flux_constraints (stokes_dof_handler, 0,
+                                                      no_normal_flux_boundaries,
+                                                      stokes_constraints,
+                                                      mapping);
+      stokes_constraints.close ();
+    }
+    {
+      temperature_constraints.clear ();
+      temperature_constraints.reinit(temperature_relevant_partitioning);//temp_locally_active);
+
+      DoFTools::make_hanging_node_constraints (temperature_dof_handler,
+                                              temperature_constraints);
+      VectorTools::interpolate_boundary_values (temperature_dof_handler,
+                                               0,
+                                               EquationData::TemperatureInitialValues<dim>(),
+                                               temperature_constraints);
+      VectorTools::interpolate_boundary_values (temperature_dof_handler,
+                                               1,
+                                               EquationData::TemperatureInitialValues<dim>(),
+                                               temperature_constraints);
+      temperature_constraints.close ();
+    }
 
-  setup_stokes_matrix (stokes_partitioning);
-  setup_stokes_preconditioner (stokes_partitioning);
-  setup_temperature_matrices (temperature_partitioning);
+    setup_stokes_matrix (stokes_partitioning);
+    setup_stokes_preconditioner (stokes_partitioning);
+    setup_temperature_matrices (temperature_partitioning);
 
-  stokes_rhs.reinit (stokes_partitioning, MPI_COMM_WORLD);
-  stokes_solution.reinit (stokes_relevant_partitioning, MPI_COMM_WORLD);
-  old_stokes_solution.reinit (stokes_solution);
+    stokes_rhs.reinit (stokes_partitioning, MPI_COMM_WORLD);
+    stokes_solution.reinit (stokes_relevant_partitioning, MPI_COMM_WORLD);
+    old_stokes_solution.reinit (stokes_solution);
 
-  temperature_rhs.reinit (temperature_partitioning, MPI_COMM_WORLD);
-  temperature_solution.reinit (temperature_relevant_partitioning, MPI_COMM_WORLD);
-  old_temperature_solution.reinit (temperature_solution);
-  old_old_temperature_solution.reinit (temperature_solution);
+    temperature_rhs.reinit (temperature_partitioning, MPI_COMM_WORLD);
+    temperature_solution.reinit (temperature_relevant_partitioning, MPI_COMM_WORLD);
+    old_temperature_solution.reinit (temperature_solution);
+    old_old_temperature_solution.reinit (temperature_solution);
 
-  rebuild_stokes_matrix              = true;
-  rebuild_stokes_preconditioner      = true;
-  rebuild_temperature_matrices       = true;
-  rebuild_temperature_preconditioner = true;
+    rebuild_stokes_matrix              = true;
+    rebuild_stokes_preconditioner      = true;
+    rebuild_temperature_matrices       = true;
+    rebuild_temperature_preconditioner = true;
 
-  computing_timer.exit_section();
-}
+    computing_timer.exit_section();
+  }
 
 
 
@@ -2320,60 +2322,60 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
 // scratch data object, which is allocated
 // only once before we start the loop over
 // all cells:
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::
-local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
-                                      Assembly::Scratch::StokesPreconditioner<dim> &scratch,
-                                      Assembly::CopyData::StokesPreconditioner<dim> &data)
-{
-  const unsigned int   dofs_per_cell   = stokes_fe.dofs_per_cell;
-  const unsigned int   n_q_points      = scratch.stokes_fe_values.n_quadrature_points;
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::
+  local_assemble_stokes_preconditioner (const typename DoFHandler<dim>::active_cell_iterator &cell,
+                                       Assembly::Scratch::StokesPreconditioner<dim> &scratch,
+                                       Assembly::CopyData::StokesPreconditioner<dim> &data)
+  {
+    const unsigned int   dofs_per_cell   = stokes_fe.dofs_per_cell;
+    const unsigned int   n_q_points      = scratch.stokes_fe_values.n_quadrature_points;
 
-  const FEValuesExtractors::Vector velocities (0);
-  const FEValuesExtractors::Scalar pressure (dim);
+    const FEValuesExtractors::Vector velocities (0);
+    const FEValuesExtractors::Scalar pressure (dim);
 
-  scratch.stokes_fe_values.reinit (cell);
-  cell->get_dof_indices (data.local_dof_indices);
+    scratch.stokes_fe_values.reinit (cell);
+    cell->get_dof_indices (data.local_dof_indices);
 
-  data.local_matrix = 0;
+    data.local_matrix = 0;
 
-  for (unsigned int q=0; q<n_q_points; ++q)
-    {
-      for (unsigned int k=0; k<dofs_per_cell; ++k)
-        {
-          scratch.grads_phi_u[k] = scratch.stokes_fe_values[velocities].symmetric_gradient(k,q);
-          scratch.phi_p[k]       = scratch.stokes_fe_values[pressure].value (k, q);
-        }
-
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-        for (unsigned int j=0; j<dofs_per_cell; ++j)
-          if (stokes_fe.system_to_component_index(i).first
-              ==
-              stokes_fe.system_to_component_index(j).first)
-            data.local_matrix(i,j) += (EquationData::eta *
-                                       (scratch.grads_phi_u[i] *
-                                        scratch.grads_phi_u[j])
-                                       +
-                                       (1./EquationData::eta) *
-                                       EquationData::pressure_scaling *
-                                       EquationData::pressure_scaling *
-                                       (scratch.phi_p[i] * scratch.phi_p[j]))
-                                      * scratch.stokes_fe_values.JxW(q);
-    }
-}
+    for (unsigned int q=0; q<n_q_points; ++q)
+      {
+       for (unsigned int k=0; k<dofs_per_cell; ++k)
+         {
+           scratch.grads_phi_u[k] = scratch.stokes_fe_values[velocities].symmetric_gradient(k,q);
+           scratch.phi_p[k]       = scratch.stokes_fe_values[pressure].value (k, q);
+         }
+
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         for (unsigned int j=0; j<dofs_per_cell; ++j)
+           if (stokes_fe.system_to_component_index(i).first
+               ==
+               stokes_fe.system_to_component_index(j).first)
+             data.local_matrix(i,j) += (EquationData::eta *
+                                        (scratch.grads_phi_u[i] *
+                                         scratch.grads_phi_u[j])
+                                        +
+                                        (1./EquationData::eta) *
+                                        EquationData::pressure_scaling *
+                                        EquationData::pressure_scaling *
+                                        (scratch.phi_p[i] * scratch.phi_p[j]))
+                                       * scratch.stokes_fe_values.JxW(q);
+      }
+  }
 
 
 
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::
-copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPreconditioner<dim> &data)
-{
-  stokes_constraints.distribute_local_to_global (data.local_matrix,
-                                                 data.local_dof_indices,
-                                                 stokes_preconditioner_matrix);
-}
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::
+  copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPreconditioner<dim> &data)
+  {
+    stokes_constraints.distribute_local_to_global (data.local_matrix,
+                                                  data.local_dof_indices,
+                                                  stokes_preconditioner_matrix);
+  }
 
 
 
@@ -2429,46 +2431,46 @@ copy_local_to_global_stokes_preconditioner (const Assembly::CopyData::StokesPrec
 // send it to the owner at the end of
 // assembly, by calling the
 // <code>compress()</code> command.
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::assemble_stokes_preconditioner ()
-{
-  stokes_preconditioner_matrix = 0;
-
-  const QGauss<dim> quadrature_formula(parameters.stokes_velocity_degree+1);
-
-  typedef
-  FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
-  SubdomainFilter;
-
-  WorkStream::
-  run (SubdomainFilter (IteratorFilters::SubdomainEqualTo
-                        (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
-                        stokes_dof_handler.begin_active()),
-       SubdomainFilter (IteratorFilters::SubdomainEqualTo
-                        (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
-                        stokes_dof_handler.end()),
-       std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
-                        local_assemble_stokes_preconditioner,
-                        this,
-                        std_cxx1x::_1,
-                        std_cxx1x::_2,
-                        std_cxx1x::_3),
-       std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
-                        copy_local_to_global_stokes_preconditioner,
-                        this,
-                        std_cxx1x::_1),
-       Assembly::Scratch::
-       StokesPreconditioner<dim> (stokes_fe, quadrature_formula,
-                                  mapping,
-                                  update_JxW_values |
-                                  update_values |
-                                  update_gradients),
-       Assembly::CopyData::
-       StokesPreconditioner<dim> (stokes_fe));
-
-  stokes_preconditioner_matrix.compress();
-}
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::assemble_stokes_preconditioner ()
+  {
+    stokes_preconditioner_matrix = 0;
+
+    const QGauss<dim> quadrature_formula(parameters.stokes_velocity_degree+1);
+
+    typedef
+      FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
+      SubdomainFilter;
+
+    WorkStream::
+      run (SubdomainFilter (IteratorFilters::SubdomainEqualTo
+                           (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
+                           stokes_dof_handler.begin_active()),
+          SubdomainFilter (IteratorFilters::SubdomainEqualTo
+                           (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
+                           stokes_dof_handler.end()),
+          std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
+                           local_assemble_stokes_preconditioner,
+                           this,
+                           std_cxx1x::_1,
+                           std_cxx1x::_2,
+                           std_cxx1x::_3),
+          std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
+                           copy_local_to_global_stokes_preconditioner,
+                           this,
+                           std_cxx1x::_1),
+          Assembly::Scratch::
+          StokesPreconditioner<dim> (stokes_fe, quadrature_formula,
+                                     mapping,
+                                     update_JxW_values |
+                                     update_values |
+                                     update_gradients),
+          Assembly::CopyData::
+          StokesPreconditioner<dim> (stokes_fe));
+
+    stokes_preconditioner_matrix.compress();
+  }
 
 
 
@@ -2481,45 +2483,45 @@ BoussinesqFlowProblem<dim>::assemble_stokes_preconditioner ()
 // preconditioner for the pressure mass
 // matrix instead of IC, as discussed in the
 // introduction.
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
-{
-  if (rebuild_stokes_preconditioner == false)
-    return;
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
+  {
+    if (rebuild_stokes_preconditioner == false)
+      return;
 
-  computing_timer.enter_section ("   Build Stokes preconditioner");
-  pcout << "   Rebuilding Stokes preconditioner..." << std::flush;
+    computing_timer.enter_section ("   Build Stokes preconditioner");
+    pcout << "   Rebuilding Stokes preconditioner..." << std::flush;
 
-  assemble_stokes_preconditioner ();
+    assemble_stokes_preconditioner ();
 
-  std::vector<std::vector<bool> > constant_modes;
-  std::vector<bool>  velocity_components (dim+1,true);
-  velocity_components[dim] = false;
-  DoFTools::extract_constant_modes (stokes_dof_handler, velocity_components,
-                                    constant_modes);
+    std::vector<std::vector<bool> > constant_modes;
+    std::vector<bool>  velocity_components (dim+1,true);
+    velocity_components[dim] = false;
+    DoFTools::extract_constant_modes (stokes_dof_handler, velocity_components,
+                                     constant_modes);
 
-  Mp_preconditioner  = std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionILU>
-                       (new TrilinosWrappers::PreconditionILU());
-  Amg_preconditioner = std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionAMG>
-                       (new TrilinosWrappers::PreconditionAMG());
+    Mp_preconditioner  = std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionILU>
+                        (new TrilinosWrappers::PreconditionILU());
+    Amg_preconditioner = std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionAMG>
+                        (new TrilinosWrappers::PreconditionAMG());
 
-  TrilinosWrappers::PreconditionAMG::AdditionalData Amg_data;
-  Amg_data.constant_modes = constant_modes;
-  Amg_data.elliptic = true;
-  Amg_data.higher_order_elements = true;
-  Amg_data.smoother_sweeps = 2;
-  Amg_data.aggregation_threshold = 0.02;
+    TrilinosWrappers::PreconditionAMG::AdditionalData Amg_data;
+    Amg_data.constant_modes = constant_modes;
+    Amg_data.elliptic = true;
+    Amg_data.higher_order_elements = true;
+    Amg_data.smoother_sweeps = 2;
+    Amg_data.aggregation_threshold = 0.02;
 
-  Mp_preconditioner->initialize (stokes_preconditioner_matrix.block(1,1));
-  Amg_preconditioner->initialize (stokes_preconditioner_matrix.block(0,0),
-                                  Amg_data);
+    Mp_preconditioner->initialize (stokes_preconditioner_matrix.block(1,1));
+    Amg_preconditioner->initialize (stokes_preconditioner_matrix.block(0,0),
+                                   Amg_data);
 
-  rebuild_stokes_preconditioner = false;
+    rebuild_stokes_preconditioner = false;
 
-  pcout << std::endl;
-  computing_timer.exit_section();
-}
+    pcout << std::endl;
+    computing_timer.exit_section();
+  }
 
 // @sect5{Stokes system assembly}
 
@@ -2542,152 +2544,152 @@ BoussinesqFlowProblem<dim>::build_stokes_preconditioner ()
 // functions in the end of the assembly in
 // order to send non-local data to the
 // owner process.
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::
-local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
-                              Assembly::Scratch::StokesSystem<dim> &scratch,
-                              Assembly::CopyData::StokesSystem<dim> &data)
-{
-  const unsigned int dofs_per_cell = scratch.stokes_fe_values.get_fe().dofs_per_cell;
-  const unsigned int n_q_points    = scratch.stokes_fe_values.n_quadrature_points;
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::
+  local_assemble_stokes_system (const typename DoFHandler<dim>::active_cell_iterator &cell,
+                               Assembly::Scratch::StokesSystem<dim> &scratch,
+                               Assembly::CopyData::StokesSystem<dim> &data)
+  {
+    const unsigned int dofs_per_cell = scratch.stokes_fe_values.get_fe().dofs_per_cell;
+    const unsigned int n_q_points    = scratch.stokes_fe_values.n_quadrature_points;
 
-  const FEValuesExtractors::Vector velocities (0);
-  const FEValuesExtractors::Scalar pressure (dim);
+    const FEValuesExtractors::Vector velocities (0);
+    const FEValuesExtractors::Scalar pressure (dim);
 
-  scratch.stokes_fe_values.reinit (cell);
+    scratch.stokes_fe_values.reinit (cell);
 
-  typename DoFHandler<dim>::active_cell_iterator
-  temperature_cell (&triangulation,
-                    cell->level(),
-                    cell->index(),
-                    &temperature_dof_handler);
-  scratch.temperature_fe_values.reinit (temperature_cell);
+    typename DoFHandler<dim>::active_cell_iterator
+      temperature_cell (&triangulation,
+                       cell->level(),
+                       cell->index(),
+                       &temperature_dof_handler);
+    scratch.temperature_fe_values.reinit (temperature_cell);
 
-  if (rebuild_stokes_matrix)
-    data.local_matrix = 0;
-  data.local_rhs = 0;
+    if (rebuild_stokes_matrix)
+      data.local_matrix = 0;
+    data.local_rhs = 0;
 
-  scratch.temperature_fe_values.get_function_values (old_temperature_solution,
-                                                     scratch.old_temperature_values);
+    scratch.temperature_fe_values.get_function_values (old_temperature_solution,
+                                                      scratch.old_temperature_values);
 
-  for (unsigned int q=0; q<n_q_points; ++q)
-    {
-      const double old_temperature = scratch.old_temperature_values[q];
-
-      for (unsigned int k=0; k<dofs_per_cell; ++k)
-        {
-          scratch.phi_u[k] = scratch.stokes_fe_values[velocities].value (k,q);
-          if (rebuild_stokes_matrix)
-            {
-              scratch.grads_phi_u[k] = scratch.stokes_fe_values[velocities].symmetric_gradient(k,q);
-              scratch.div_phi_u[k]   = scratch.stokes_fe_values[velocities].divergence (k, q);
-              scratch.phi_p[k]       = scratch.stokes_fe_values[pressure].value (k, q);
-            }
-        }
-
-      if (rebuild_stokes_matrix)
-        for (unsigned int i=0; i<dofs_per_cell; ++i)
-          for (unsigned int j=0; j<dofs_per_cell; ++j)
-            data.local_matrix(i,j) += (EquationData::eta * 2 *
-                                       (scratch.grads_phi_u[i] * scratch.grads_phi_u[j])
-                                       - (EquationData::pressure_scaling *
-                                          scratch.div_phi_u[i] * scratch.phi_p[j])
-                                       - (EquationData::pressure_scaling *
-                                          scratch.phi_p[i] * scratch.div_phi_u[j]))
-                                      * scratch.stokes_fe_values.JxW(q);
-
-      const Tensor<1,dim>
-      gravity = EquationData::gravity_vector (scratch.stokes_fe_values
-                                              .quadrature_point(q));
-
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-        data.local_rhs(i) += (EquationData::density(old_temperature) *
-                              gravity  *
-                              scratch.phi_u[i]) *
-                             scratch.stokes_fe_values.JxW(q);
-    }
+    for (unsigned int q=0; q<n_q_points; ++q)
+      {
+       const double old_temperature = scratch.old_temperature_values[q];
+
+       for (unsigned int k=0; k<dofs_per_cell; ++k)
+         {
+           scratch.phi_u[k] = scratch.stokes_fe_values[velocities].value (k,q);
+           if (rebuild_stokes_matrix)
+             {
+               scratch.grads_phi_u[k] = scratch.stokes_fe_values[velocities].symmetric_gradient(k,q);
+               scratch.div_phi_u[k]   = scratch.stokes_fe_values[velocities].divergence (k, q);
+               scratch.phi_p[k]       = scratch.stokes_fe_values[pressure].value (k, q);
+             }
+         }
+
+       if (rebuild_stokes_matrix)
+         for (unsigned int i=0; i<dofs_per_cell; ++i)
+           for (unsigned int j=0; j<dofs_per_cell; ++j)
+             data.local_matrix(i,j) += (EquationData::eta * 2 *
+                                        (scratch.grads_phi_u[i] * scratch.grads_phi_u[j])
+                                        - (EquationData::pressure_scaling *
+                                           scratch.div_phi_u[i] * scratch.phi_p[j])
+                                        - (EquationData::pressure_scaling *
+                                           scratch.phi_p[i] * scratch.div_phi_u[j]))
+                                       * scratch.stokes_fe_values.JxW(q);
+
+       const Tensor<1,dim>
+         gravity = EquationData::gravity_vector (scratch.stokes_fe_values
+                                                 .quadrature_point(q));
+
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         data.local_rhs(i) += (EquationData::density(old_temperature) *
+                               gravity  *
+                               scratch.phi_u[i]) *
+                              scratch.stokes_fe_values.JxW(q);
+      }
 
-  cell->get_dof_indices (data.local_dof_indices);
-}
+    cell->get_dof_indices (data.local_dof_indices);
+  }
 
 
 
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::
-copy_local_to_global_stokes_system (const Assembly::CopyData::StokesSystem<dim> &data)
-{
-  if (rebuild_stokes_matrix == true)
-    stokes_constraints.distribute_local_to_global (data.local_matrix,
-                                                   data.local_rhs,
-                                                   data.local_dof_indices,
-                                                   stokes_matrix,
-                                                   stokes_rhs);
-  else
-    stokes_constraints.distribute_local_to_global (data.local_rhs,
-                                                   data.local_dof_indices,
-                                                   stokes_rhs);
-}
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::
+  copy_local_to_global_stokes_system (const Assembly::CopyData::StokesSystem<dim> &data)
+  {
+    if (rebuild_stokes_matrix == true)
+      stokes_constraints.distribute_local_to_global (data.local_matrix,
+                                                    data.local_rhs,
+                                                    data.local_dof_indices,
+                                                    stokes_matrix,
+                                                    stokes_rhs);
+    else
+      stokes_constraints.distribute_local_to_global (data.local_rhs,
+                                                    data.local_dof_indices,
+                                                    stokes_rhs);
+  }
 
 
 
-template <int dim>
-void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
-{
-  computing_timer.enter_section ("   Assemble Stokes system");
-
-  if (rebuild_stokes_matrix == true)
-    stokes_matrix=0;
-
-  stokes_rhs=0;
-
-  const QGauss<dim> quadrature_formula(parameters.stokes_velocity_degree+1);
-
-  typedef
-  FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
-  SubdomainFilter;
-
-  WorkStream::
-  run (SubdomainFilter (IteratorFilters::SubdomainEqualTo
-                        (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
-                        stokes_dof_handler.begin_active()),
-       SubdomainFilter (IteratorFilters::SubdomainEqualTo
-                        (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
-                        stokes_dof_handler.end()),
-       std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
-                        local_assemble_stokes_system,
-                        this,
-                        std_cxx1x::_1,
-                        std_cxx1x::_2,
-                        std_cxx1x::_3),
-       std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
-                        copy_local_to_global_stokes_system,
-                        this,
-                        std_cxx1x::_1),
-       Assembly::Scratch::
-       StokesSystem<dim> (stokes_fe, mapping, quadrature_formula,
-                          (update_values    |
-                           update_quadrature_points  |
-                           update_JxW_values |
-                           (rebuild_stokes_matrix == true
-                            ?
-                            update_gradients
-                            :
-                            UpdateFlags(0))),
-                          temperature_fe,
-                          update_values),
-       Assembly::CopyData::
-       StokesSystem<dim> (stokes_fe));
-
-  stokes_matrix.compress();
-  stokes_rhs.compress(Add);
-
-  rebuild_stokes_matrix = false;
-
-  pcout << std::endl;
-  computing_timer.exit_section();
-}
+  template <int dim>
+  void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
+  {
+    computing_timer.enter_section ("   Assemble Stokes system");
+
+    if (rebuild_stokes_matrix == true)
+      stokes_matrix=0;
+
+    stokes_rhs=0;
+
+    const QGauss<dim> quadrature_formula(parameters.stokes_velocity_degree+1);
+
+    typedef
+      FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
+      SubdomainFilter;
+
+    WorkStream::
+      run (SubdomainFilter (IteratorFilters::SubdomainEqualTo
+                           (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
+                           stokes_dof_handler.begin_active()),
+          SubdomainFilter (IteratorFilters::SubdomainEqualTo
+                           (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
+                           stokes_dof_handler.end()),
+          std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
+                           local_assemble_stokes_system,
+                           this,
+                           std_cxx1x::_1,
+                           std_cxx1x::_2,
+                           std_cxx1x::_3),
+          std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
+                           copy_local_to_global_stokes_system,
+                           this,
+                           std_cxx1x::_1),
+          Assembly::Scratch::
+          StokesSystem<dim> (stokes_fe, mapping, quadrature_formula,
+                             (update_values    |
+                              update_quadrature_points  |
+                              update_JxW_values |
+                              (rebuild_stokes_matrix == true
+                               ?
+                               update_gradients
+                               :
+                               UpdateFlags(0))),
+                             temperature_fe,
+                             update_values),
+          Assembly::CopyData::
+          StokesSystem<dim> (stokes_fe));
+
+    stokes_matrix.compress();
+    stokes_rhs.compress(Add);
+
+    rebuild_stokes_matrix = false;
+
+    pcout << std::endl;
+    computing_timer.exit_section();
+  }
 
 
 // @sect5{Temperature matrix assembly}
@@ -2705,106 +2707,106 @@ void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
 //
 // The two following functions perform
 // similar services as the ones above.
-template <int dim>
-void BoussinesqFlowProblem<dim>::
-local_assemble_temperature_matrix (const typename DoFHandler<dim>::active_cell_iterator &cell,
-                                   Assembly::Scratch::TemperatureMatrix<dim> &scratch,
-                                   Assembly::CopyData::TemperatureMatrix<dim> &data)
-{
-  const unsigned int dofs_per_cell = scratch.temperature_fe_values.get_fe().dofs_per_cell;
-  const unsigned int n_q_points    = scratch.temperature_fe_values.n_quadrature_points;
+  template <int dim>
+  void BoussinesqFlowProblem<dim>::
+  local_assemble_temperature_matrix (const typename DoFHandler<dim>::active_cell_iterator &cell,
+                                    Assembly::Scratch::TemperatureMatrix<dim> &scratch,
+                                    Assembly::CopyData::TemperatureMatrix<dim> &data)
+  {
+    const unsigned int dofs_per_cell = scratch.temperature_fe_values.get_fe().dofs_per_cell;
+    const unsigned int n_q_points    = scratch.temperature_fe_values.n_quadrature_points;
 
-  scratch.temperature_fe_values.reinit (cell);
-  cell->get_dof_indices (data.local_dof_indices);
+    scratch.temperature_fe_values.reinit (cell);
+    cell->get_dof_indices (data.local_dof_indices);
 
-  data.local_mass_matrix = 0;
-  data.local_stiffness_matrix = 0;
+    data.local_mass_matrix = 0;
+    data.local_stiffness_matrix = 0;
 
-  for (unsigned int q=0; q<n_q_points; ++q)
-    {
-      for (unsigned int k=0; k<dofs_per_cell; ++k)
-        {
-          scratch.grad_phi_T[k] = scratch.temperature_fe_values.shape_grad (k,q);
-          scratch.phi_T[k]      = scratch.temperature_fe_values.shape_value (k, q);
-        }
-
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-        for (unsigned int j=0; j<dofs_per_cell; ++j)
-          {
-            data.local_mass_matrix(i,j)
-            += (scratch.phi_T[i] * scratch.phi_T[j]
-                *
-                scratch.temperature_fe_values.JxW(q));
-            data.local_stiffness_matrix(i,j)
-            += (EquationData::kappa * scratch.grad_phi_T[i] * scratch.grad_phi_T[j]
-                *
-                scratch.temperature_fe_values.JxW(q));
-          }
-    }
-}
+    for (unsigned int q=0; q<n_q_points; ++q)
+      {
+       for (unsigned int k=0; k<dofs_per_cell; ++k)
+         {
+           scratch.grad_phi_T[k] = scratch.temperature_fe_values.shape_grad (k,q);
+           scratch.phi_T[k]      = scratch.temperature_fe_values.shape_value (k, q);
+         }
+
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         for (unsigned int j=0; j<dofs_per_cell; ++j)
+           {
+             data.local_mass_matrix(i,j)
+               += (scratch.phi_T[i] * scratch.phi_T[j]
+                   *
+                   scratch.temperature_fe_values.JxW(q));
+             data.local_stiffness_matrix(i,j)
+               += (EquationData::kappa * scratch.grad_phi_T[i] * scratch.grad_phi_T[j]
+                   *
+                   scratch.temperature_fe_values.JxW(q));
+           }
+      }
+  }
 
 
 
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::
-copy_local_to_global_temperature_matrix (const Assembly::CopyData::TemperatureMatrix<dim> &data)
-{
-  temperature_constraints.distribute_local_to_global (data.local_mass_matrix,
-                                                      data.local_dof_indices,
-                                                      temperature_mass_matrix);
-  temperature_constraints.distribute_local_to_global (data.local_stiffness_matrix,
-                                                      data.local_dof_indices,
-                                                      temperature_stiffness_matrix);
-}
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::
+  copy_local_to_global_temperature_matrix (const Assembly::CopyData::TemperatureMatrix<dim> &data)
+  {
+    temperature_constraints.distribute_local_to_global (data.local_mass_matrix,
+                                                       data.local_dof_indices,
+                                                       temperature_mass_matrix);
+    temperature_constraints.distribute_local_to_global (data.local_stiffness_matrix,
+                                                       data.local_dof_indices,
+                                                       temperature_stiffness_matrix);
+  }
 
 
-template <int dim>
-void BoussinesqFlowProblem<dim>::assemble_temperature_matrix ()
-{
-  if (rebuild_temperature_matrices == false)
-    return;
-
-  computing_timer.enter_section ("   Assemble temperature matrices");
-  temperature_mass_matrix = 0;
-  temperature_stiffness_matrix = 0;
-
-  const QGauss<dim> quadrature_formula(parameters.temperature_degree+2);
-
-  typedef
-  FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
-  SubdomainFilter;
-
-  WorkStream::
-  run (SubdomainFilter (IteratorFilters::SubdomainEqualTo
-                        (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
-                        temperature_dof_handler.begin_active()),
-       SubdomainFilter (IteratorFilters::SubdomainEqualTo
-                        (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
-                        temperature_dof_handler.end()),
-       std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
-                        local_assemble_temperature_matrix,
-                        this,
-                        std_cxx1x::_1,
-                        std_cxx1x::_2,
-                        std_cxx1x::_3),
-       std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
-                        copy_local_to_global_temperature_matrix,
-                        this,
-                        std_cxx1x::_1),
-       Assembly::Scratch::
-       TemperatureMatrix<dim> (temperature_fe, mapping, quadrature_formula),
-       Assembly::CopyData::
-       TemperatureMatrix<dim> (temperature_fe));
-
-  temperature_mass_matrix.compress();
-  temperature_stiffness_matrix.compress();
-
-  rebuild_temperature_matrices = false;
-  rebuild_temperature_preconditioner = true;
-
-  computing_timer.exit_section();
-}
+  template <int dim>
+  void BoussinesqFlowProblem<dim>::assemble_temperature_matrix ()
+  {
+    if (rebuild_temperature_matrices == false)
+      return;
+
+    computing_timer.enter_section ("   Assemble temperature matrices");
+    temperature_mass_matrix = 0;
+    temperature_stiffness_matrix = 0;
+
+    const QGauss<dim> quadrature_formula(parameters.temperature_degree+2);
+
+    typedef
+      FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
+      SubdomainFilter;
+
+    WorkStream::
+      run (SubdomainFilter (IteratorFilters::SubdomainEqualTo
+                           (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
+                           temperature_dof_handler.begin_active()),
+          SubdomainFilter (IteratorFilters::SubdomainEqualTo
+                           (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
+                           temperature_dof_handler.end()),
+          std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
+                           local_assemble_temperature_matrix,
+                           this,
+                           std_cxx1x::_1,
+                           std_cxx1x::_2,
+                           std_cxx1x::_3),
+          std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
+                           copy_local_to_global_temperature_matrix,
+                           this,
+                           std_cxx1x::_1),
+          Assembly::Scratch::
+          TemperatureMatrix<dim> (temperature_fe, mapping, quadrature_formula),
+          Assembly::CopyData::
+          TemperatureMatrix<dim> (temperature_fe));
+
+    temperature_mass_matrix.compress();
+    temperature_stiffness_matrix.compress();
+
+    rebuild_temperature_matrices = false;
+    rebuild_temperature_preconditioner = true;
+
+    computing_timer.exit_section();
+  }
 
 
 // @sect5{Temperature right hand side assembly}
@@ -2831,170 +2833,170 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_matrix ()
 // dofs. That will account for the correct
 // balance of the right hand side vector
 // with the matrix system of temperature.
-template <int dim>
-void BoussinesqFlowProblem<dim>::
-local_assemble_temperature_rhs (const std::pair<double,double> global_T_range,
-                                const double                   global_max_velocity,
-                                const double                   global_entropy_variation,
-                                const typename DoFHandler<dim>::active_cell_iterator &cell,
-                                Assembly::Scratch::TemperatureRHS<dim> &scratch,
-                                Assembly::CopyData::TemperatureRHS<dim> &data)
-{
-  const bool use_bdf2_scheme = (timestep_number != 0);
-
-  const unsigned int dofs_per_cell = scratch.temperature_fe_values.get_fe().dofs_per_cell;
-  const unsigned int n_q_points    = scratch.temperature_fe_values.n_quadrature_points;
-
-  const FEValuesExtractors::Vector velocities (0);
-
-  data.local_rhs = 0;
-  data.matrix_for_bc = 0;
-  cell->get_dof_indices (data.local_dof_indices);
-
-  scratch.temperature_fe_values.reinit (cell);
-
-  typename DoFHandler<dim>::active_cell_iterator
-  stokes_cell (&triangulation,
-               cell->level(),
-               cell->index(),
-               &stokes_dof_handler);
-  scratch.stokes_fe_values.reinit (stokes_cell);
-
-  scratch.temperature_fe_values.get_function_values (old_temperature_solution,
-                                                     scratch.old_temperature_values);
-  scratch.temperature_fe_values.get_function_values (old_old_temperature_solution,
-                                                     scratch.old_old_temperature_values);
-
-  scratch.temperature_fe_values.get_function_gradients (old_temperature_solution,
-                                                        scratch.old_temperature_grads);
-  scratch.temperature_fe_values.get_function_gradients (old_old_temperature_solution,
-                                                        scratch.old_old_temperature_grads);
-
-  scratch.temperature_fe_values.get_function_laplacians (old_temperature_solution,
-                                                         scratch.old_temperature_laplacians);
-  scratch.temperature_fe_values.get_function_laplacians (old_old_temperature_solution,
-                                                         scratch.old_old_temperature_laplacians);
-
-  scratch.stokes_fe_values[velocities].get_function_values (stokes_solution,
-                                                            scratch.old_velocity_values);
-  scratch.stokes_fe_values[velocities].get_function_values (old_stokes_solution,
-                                                            scratch.old_old_velocity_values);
-  scratch.stokes_fe_values[velocities].get_function_symmetric_gradients (stokes_solution,
-                                                                         scratch.old_strain_rates);
-  scratch.stokes_fe_values[velocities].get_function_symmetric_gradients (old_stokes_solution,
-                                                                         scratch.old_old_strain_rates);
-
-  const double nu
-    = compute_viscosity (scratch.old_temperature_values,
-                         scratch.old_old_temperature_values,
-                         scratch.old_temperature_grads,
-                         scratch.old_old_temperature_grads,
-                         scratch.old_temperature_laplacians,
-                         scratch.old_old_temperature_laplacians,
-                         scratch.old_velocity_values,
-                         scratch.old_old_velocity_values,
-                         scratch.old_strain_rates,
-                         scratch.old_old_strain_rates,
-                         global_max_velocity,
-                         global_T_range.second - global_T_range.first,
-                         0.5 * (global_T_range.second + global_T_range.first),
-                         global_entropy_variation,
-                         cell->diameter());
-
-  for (unsigned int q=0; q<n_q_points; ++q)
-    {
-      for (unsigned int k=0; k<dofs_per_cell; ++k)
-        {
-          scratch.phi_T[k]      = scratch.temperature_fe_values.shape_value (k, q);
-          scratch.grad_phi_T[k] = scratch.temperature_fe_values.shape_grad (k,q);
-        }
-
-
-      const double old_Ts
-        = (use_bdf2_scheme ?
-           (scratch.old_temperature_values[q] *
-            (time_step + old_time_step) / old_time_step
-            -
-            scratch.old_old_temperature_values[q] *
-            (time_step * time_step) /
-            (old_time_step * (time_step + old_time_step)))
-           :
-           scratch.old_temperature_values[q]);
-
-      const Tensor<1,dim> ext_grad_T
-        = (use_bdf2_scheme ?
-           (scratch.old_temperature_grads[q] *
-            (1+time_step/old_time_step)
-            -
-            scratch.old_old_temperature_grads[q] *
-            time_step / old_time_step)
-           :
-           scratch.old_temperature_grads[q]);
-
-      const Tensor<1,dim> extrapolated_u
-        = (use_bdf2_scheme ?
-           (scratch.old_velocity_values[q] * (1+time_step/old_time_step) -
-            scratch.old_old_velocity_values[q] * time_step/old_time_step)
-           :
-           scratch.old_velocity_values[q]);
-      const SymmetricTensor<2,dim> extrapolated_strain_rate
-        = (use_bdf2_scheme ?
-           (scratch.old_strain_rates[q] * (1+time_step/old_time_step) -
-            scratch.old_old_strain_rates[q] * time_step/old_time_step)
-           :
-           scratch.old_strain_rates[q]);
-
-      const double gamma
-        = ((EquationData::radiogenic_heating * EquationData::density(old_Ts) //?????? why old_Ts?
-            +
-            2 * EquationData::eta * extrapolated_strain_rate * extrapolated_strain_rate) /
-           (EquationData::density(old_Ts) * EquationData::specific_heat));
-
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-        {
-          data.local_rhs(i) += (old_Ts * scratch.phi_T[i]
-                                -
-                                time_step *
-                                extrapolated_u * ext_grad_T * scratch.phi_T[i]
-                                -
-                                time_step *
-                                nu * ext_grad_T * scratch.grad_phi_T[i]
-                                +
-                                time_step *
-                                gamma * scratch.phi_T[i])
-                               *
-                               scratch.temperature_fe_values.JxW(q);
-
-          if (temperature_constraints.is_inhomogeneously_constrained(data.local_dof_indices[i]))
-            {
-              for (unsigned int j=0; j<dofs_per_cell; ++j)
-                data.matrix_for_bc(j,i) += (scratch.phi_T[i] * scratch.phi_T[j] *
-                                            (use_bdf2_scheme ?
-                                             ((2*time_step + old_time_step) /
-                                              (time_step + old_time_step)) : 1.)
-                                            +
-                                            scratch.grad_phi_T[i] *
-                                            scratch.grad_phi_T[j] *
-                                            EquationData::kappa *
-                                            time_step)
-                                           *
-                                           scratch.temperature_fe_values.JxW(q);
-            }
-        }
-    }
-}
+  template <int dim>
+  void BoussinesqFlowProblem<dim>::
+  local_assemble_temperature_rhs (const std::pair<double,double> global_T_range,
+                                 const double                   global_max_velocity,
+                                 const double                   global_entropy_variation,
+                                 const typename DoFHandler<dim>::active_cell_iterator &cell,
+                                 Assembly::Scratch::TemperatureRHS<dim> &scratch,
+                                 Assembly::CopyData::TemperatureRHS<dim> &data)
+  {
+    const bool use_bdf2_scheme = (timestep_number != 0);
 
+    const unsigned int dofs_per_cell = scratch.temperature_fe_values.get_fe().dofs_per_cell;
+    const unsigned int n_q_points    = scratch.temperature_fe_values.n_quadrature_points;
 
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::
-copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS<dim> &data)
-{
-  temperature_constraints.distribute_local_to_global (data.local_rhs,
-                                                      data.local_dof_indices,
-                                                      temperature_rhs,
-                                                      data.matrix_for_bc);
-}
+    const FEValuesExtractors::Vector velocities (0);
+
+    data.local_rhs = 0;
+    data.matrix_for_bc = 0;
+    cell->get_dof_indices (data.local_dof_indices);
+
+    scratch.temperature_fe_values.reinit (cell);
+
+    typename DoFHandler<dim>::active_cell_iterator
+      stokes_cell (&triangulation,
+                  cell->level(),
+                  cell->index(),
+                  &stokes_dof_handler);
+    scratch.stokes_fe_values.reinit (stokes_cell);
+
+    scratch.temperature_fe_values.get_function_values (old_temperature_solution,
+                                                      scratch.old_temperature_values);
+    scratch.temperature_fe_values.get_function_values (old_old_temperature_solution,
+                                                      scratch.old_old_temperature_values);
+
+    scratch.temperature_fe_values.get_function_gradients (old_temperature_solution,
+                                                         scratch.old_temperature_grads);
+    scratch.temperature_fe_values.get_function_gradients (old_old_temperature_solution,
+                                                         scratch.old_old_temperature_grads);
+
+    scratch.temperature_fe_values.get_function_laplacians (old_temperature_solution,
+                                                          scratch.old_temperature_laplacians);
+    scratch.temperature_fe_values.get_function_laplacians (old_old_temperature_solution,
+                                                          scratch.old_old_temperature_laplacians);
+
+    scratch.stokes_fe_values[velocities].get_function_values (stokes_solution,
+                                                             scratch.old_velocity_values);
+    scratch.stokes_fe_values[velocities].get_function_values (old_stokes_solution,
+                                                             scratch.old_old_velocity_values);
+    scratch.stokes_fe_values[velocities].get_function_symmetric_gradients (stokes_solution,
+                                                                          scratch.old_strain_rates);
+    scratch.stokes_fe_values[velocities].get_function_symmetric_gradients (old_stokes_solution,
+                                                                          scratch.old_old_strain_rates);
+
+    const double nu
+      = compute_viscosity (scratch.old_temperature_values,
+                          scratch.old_old_temperature_values,
+                          scratch.old_temperature_grads,
+                          scratch.old_old_temperature_grads,
+                          scratch.old_temperature_laplacians,
+                          scratch.old_old_temperature_laplacians,
+                          scratch.old_velocity_values,
+                          scratch.old_old_velocity_values,
+                          scratch.old_strain_rates,
+                          scratch.old_old_strain_rates,
+                          global_max_velocity,
+                          global_T_range.second - global_T_range.first,
+                          0.5 * (global_T_range.second + global_T_range.first),
+                          global_entropy_variation,
+                          cell->diameter());
+
+    for (unsigned int q=0; q<n_q_points; ++q)
+      {
+       for (unsigned int k=0; k<dofs_per_cell; ++k)
+         {
+           scratch.phi_T[k]      = scratch.temperature_fe_values.shape_value (k, q);
+           scratch.grad_phi_T[k] = scratch.temperature_fe_values.shape_grad (k,q);
+         }
+
+
+       const double old_Ts
+         = (use_bdf2_scheme ?
+            (scratch.old_temperature_values[q] *
+             (time_step + old_time_step) / old_time_step
+             -
+             scratch.old_old_temperature_values[q] *
+             (time_step * time_step) /
+             (old_time_step * (time_step + old_time_step)))
+            :
+            scratch.old_temperature_values[q]);
+
+       const Tensor<1,dim> ext_grad_T
+         = (use_bdf2_scheme ?
+            (scratch.old_temperature_grads[q] *
+             (1+time_step/old_time_step)
+             -
+             scratch.old_old_temperature_grads[q] *
+             time_step / old_time_step)
+            :
+            scratch.old_temperature_grads[q]);
+
+       const Tensor<1,dim> extrapolated_u
+         = (use_bdf2_scheme ?
+            (scratch.old_velocity_values[q] * (1+time_step/old_time_step) -
+             scratch.old_old_velocity_values[q] * time_step/old_time_step)
+            :
+            scratch.old_velocity_values[q]);
+       const SymmetricTensor<2,dim> extrapolated_strain_rate
+         = (use_bdf2_scheme ?
+            (scratch.old_strain_rates[q] * (1+time_step/old_time_step) -
+             scratch.old_old_strain_rates[q] * time_step/old_time_step)
+            :
+            scratch.old_strain_rates[q]);
+
+       const double gamma
+         = ((EquationData::radiogenic_heating * EquationData::density(old_Ts) //?????? why old_Ts?
+             +
+             2 * EquationData::eta * extrapolated_strain_rate * extrapolated_strain_rate) /
+            (EquationData::density(old_Ts) * EquationData::specific_heat));
+
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         {
+           data.local_rhs(i) += (old_Ts * scratch.phi_T[i]
+                                 -
+                                 time_step *
+                                 extrapolated_u * ext_grad_T * scratch.phi_T[i]
+                                 -
+                                 time_step *
+                                 nu * ext_grad_T * scratch.grad_phi_T[i]
+                                 +
+                                 time_step *
+                                 gamma * scratch.phi_T[i])
+                                *
+                                scratch.temperature_fe_values.JxW(q);
+
+           if (temperature_constraints.is_inhomogeneously_constrained(data.local_dof_indices[i]))
+             {
+               for (unsigned int j=0; j<dofs_per_cell; ++j)
+                 data.matrix_for_bc(j,i) += (scratch.phi_T[i] * scratch.phi_T[j] *
+                                             (use_bdf2_scheme ?
+                                              ((2*time_step + old_time_step) /
+                                               (time_step + old_time_step)) : 1.)
+                                             +
+                                             scratch.grad_phi_T[i] *
+                                             scratch.grad_phi_T[j] *
+                                             EquationData::kappa *
+                                             time_step)
+                                            *
+                                            scratch.temperature_fe_values.JxW(q);
+             }
+         }
+      }
+  }
+
+
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::
+  copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS<dim> &data)
+  {
+    temperature_constraints.distribute_local_to_global (data.local_rhs,
+                                                       data.local_dof_indices,
+                                                       temperature_rhs,
+                                                       data.matrix_for_bc);
+  }
 
 
 
@@ -3024,82 +3026,82 @@ copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS<d
 // too big a problem because we remesh every
 // fifth time step (and regenerate the
 // preconditioner then).
-template <int dim>
-void BoussinesqFlowProblem<dim>::assemble_temperature_system (const double maximal_velocity)
-{
-  const bool use_bdf2_scheme = (timestep_number != 0);
+  template <int dim>
+  void BoussinesqFlowProblem<dim>::assemble_temperature_system (const double maximal_velocity)
+  {
+    const bool use_bdf2_scheme = (timestep_number != 0);
 
-  if (use_bdf2_scheme == true)
-    {
-      temperature_matrix.copy_from (temperature_mass_matrix);
-      temperature_matrix *= (2*time_step + old_time_step) /
-                            (time_step + old_time_step);
-      temperature_matrix.add (time_step, temperature_stiffness_matrix);
-    }
-  else
-    {
-      temperature_matrix.copy_from (temperature_mass_matrix);
-      temperature_matrix.add (time_step, temperature_stiffness_matrix);
-    }
-  temperature_matrix.compress();
+    if (use_bdf2_scheme == true)
+      {
+       temperature_matrix.copy_from (temperature_mass_matrix);
+       temperature_matrix *= (2*time_step + old_time_step) /
+                             (time_step + old_time_step);
+       temperature_matrix.add (time_step, temperature_stiffness_matrix);
+      }
+    else
+      {
+       temperature_matrix.copy_from (temperature_mass_matrix);
+       temperature_matrix.add (time_step, temperature_stiffness_matrix);
+      }
+    temperature_matrix.compress();
 
-  if (rebuild_temperature_preconditioner == true)
-    {
-      T_preconditioner =  std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionIC>
-                          (new TrilinosWrappers::PreconditionIC());
-      T_preconditioner->initialize (temperature_matrix);
-      rebuild_temperature_preconditioner = false;
-    }
+    if (rebuild_temperature_preconditioner == true)
+      {
+       T_preconditioner =  std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionIC>
+                           (new TrilinosWrappers::PreconditionIC());
+       T_preconditioner->initialize (temperature_matrix);
+       rebuild_temperature_preconditioner = false;
+      }
 
-  temperature_rhs = 0;
-
-  const QGauss<dim> quadrature_formula(parameters.temperature_degree+2);
-  const std::pair<double,double>
-  global_T_range = get_extrapolated_temperature_range();
-
-  // use midpoint between maximum and minimum
-  // temperature for definition of average
-  // temperature in entropy viscosity. Could
-  // also use the integral average, but the
-  // results are not very sensitive to this
-  // choice.
-  const double average_temperature = 0.5 * (global_T_range.first +
-                                            global_T_range.second);
-  const double global_entropy_variation =
-    get_entropy_variation (average_temperature);
-
-  typedef
-  FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
-  SubdomainFilter;
-
-  WorkStream::
-  run (SubdomainFilter (IteratorFilters::SubdomainEqualTo
-                        (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
-                        temperature_dof_handler.begin_active()),
-       SubdomainFilter (IteratorFilters::SubdomainEqualTo
-                        (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
-                        temperature_dof_handler.end()),
-       std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
-                        local_assemble_temperature_rhs,
-                        this,
-                        global_T_range,
-                        maximal_velocity,
-                        global_entropy_variation,
-                        std_cxx1x::_1,
-                        std_cxx1x::_2,
-                        std_cxx1x::_3),
-       std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
-                        copy_local_to_global_temperature_rhs,
-                        this,
-                        std_cxx1x::_1),
-       Assembly::Scratch::
-       TemperatureRHS<dim> (temperature_fe, stokes_fe, mapping,
-                            quadrature_formula),
-       Assembly::CopyData::
-       TemperatureRHS<dim> (temperature_fe));
-
-  temperature_rhs.compress(Add);
-}
+    temperature_rhs = 0;
+
+    const QGauss<dim> quadrature_formula(parameters.temperature_degree+2);
+    const std::pair<double,double>
+      global_T_range = get_extrapolated_temperature_range();
+
+                                    // use midpoint between maximum and minimum
+                                    // temperature for definition of average
+                                    // temperature in entropy viscosity. Could
+                                    // also use the integral average, but the
+                                    // results are not very sensitive to this
+                                    // choice.
+    const double average_temperature = 0.5 * (global_T_range.first +
+                                             global_T_range.second);
+    const double global_entropy_variation =
+      get_entropy_variation (average_temperature);
+
+    typedef
+      FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
+      SubdomainFilter;
+
+    WorkStream::
+      run (SubdomainFilter (IteratorFilters::SubdomainEqualTo
+                           (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
+                           temperature_dof_handler.begin_active()),
+          SubdomainFilter (IteratorFilters::SubdomainEqualTo
+                           (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)),
+                           temperature_dof_handler.end()),
+          std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
+                           local_assemble_temperature_rhs,
+                           this,
+                           global_T_range,
+                           maximal_velocity,
+                           global_entropy_variation,
+                           std_cxx1x::_1,
+                           std_cxx1x::_2,
+                           std_cxx1x::_3),
+          std_cxx1x::bind (&BoussinesqFlowProblem<dim>::
+                           copy_local_to_global_temperature_rhs,
+                           this,
+                           std_cxx1x::_1),
+          Assembly::Scratch::
+          TemperatureRHS<dim> (temperature_fe, stokes_fe, mapping,
+                               quadrature_formula),
+          Assembly::CopyData::
+          TemperatureRHS<dim> (temperature_fe));
+
+    temperature_rhs.compress(Add);
+  }
 
 
 
@@ -3150,333 +3152,333 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system (const double maxim
 // Apart from these two changes, everything
 // is the same as in step-31, so we don't
 // need to further comment on it.
-template <int dim>
-void BoussinesqFlowProblem<dim>::solve ()
-{
-  computing_timer.enter_section ("   Solve Stokes system");
-
+  template <int dim>
+  void BoussinesqFlowProblem<dim>::solve ()
   {
-    pcout << "   Solving Stokes system... " << std::flush;
-
-    TrilinosWrappers::MPI::BlockVector
-    distributed_stokes_solution (stokes_rhs);
-//    distributed_stokes_solution = stokes_solution;
-    distributed_stokes_solution.block(0).reinit(stokes_solution.block(0),false,true);
-    distributed_stokes_solution.block(1).reinit(stokes_solution.block(1),false,true);
-
-
-    const unsigned int
-    start = (distributed_stokes_solution.block(0).size() +
-             distributed_stokes_solution.block(1).local_range().first),
-            end   = (distributed_stokes_solution.block(0).size() +
-                     distributed_stokes_solution.block(1).local_range().second);
-    for (unsigned int i=start; i<end; ++i)
-      if (stokes_constraints.is_constrained (i))
-        distributed_stokes_solution(i) = 0;
-
-
-    PrimitiveVectorMemory< TrilinosWrappers::MPI::BlockVector > mem;
-
-    // step 1: try if the simple and fast solver
-    // succeeds in 30 steps or less.
-    unsigned int n_iterations = 0;
-    double reduction = 0;
-    const double solver_tolerance = 1e-7 * stokes_rhs.l2_norm();
-    SolverControl solver_control (30, solver_tolerance);
-
-    try
-      {
-        const LinearSolvers::RightPrecond<TrilinosWrappers::PreconditionAMG,
-              TrilinosWrappers::PreconditionILU>
-              preconditioner (stokes_matrix, stokes_preconditioner_matrix,
-                              *Mp_preconditioner, *Amg_preconditioner,
-                              false);
-
-        SolverFGMRES<TrilinosWrappers::MPI::BlockVector>
-        solver(solver_control, mem,
-               SolverFGMRES<TrilinosWrappers::MPI::BlockVector>::
-               AdditionalData(30, true));
-        solver.solve(stokes_matrix, distributed_stokes_solution, stokes_rhs,
-                     preconditioner);
-
-        n_iterations = solver_control.last_step();
-        reduction = solver_control.last_value()/solver_control.initial_value();
-      }
-
-    // step 2: take the stronger solver in case
-    // the simple solver failed
-    catch (SolverControl::NoConvergence)
-      {
-        const LinearSolvers::RightPrecond<TrilinosWrappers::PreconditionAMG,
-              TrilinosWrappers::PreconditionILU>
-              preconditioner (stokes_matrix, stokes_preconditioner_matrix,
-                              *Mp_preconditioner, *Amg_preconditioner,
-                              true);
-
-        SolverControl solver_control_refined (stokes_matrix.m(), solver_tolerance);
-        SolverFGMRES<TrilinosWrappers::MPI::BlockVector>
-        solver(solver_control_refined, mem,
-               SolverFGMRES<TrilinosWrappers::MPI::BlockVector>::
-               AdditionalData(50, true));
-        solver.solve(stokes_matrix, distributed_stokes_solution, stokes_rhs,
-                     preconditioner);
-
-        n_iterations = (solver_control.last_step() +
-                        solver_control_refined.last_step());
-        reduction = (solver_control_refined.last_value()/
-                     std::max(solver_control.initial_value(),
-                              solver_control_refined.initial_value()));
-      }
+    computing_timer.enter_section ("   Solve Stokes system");
 
+    {
+      pcout << "   Solving Stokes system... " << std::flush;
 
-    stokes_constraints.distribute (distributed_stokes_solution);
-    stokes_solution.block(0).reinit(distributed_stokes_solution.block(0),
-                                    false, true);
-    stokes_solution.block(1).reinit(distributed_stokes_solution.block(1),
-                                    false, true);
-
-    pcout << n_iterations  << " iterations."
-          << " Reduced residual by " << reduction
-          << std::endl;
-
-    TrilinosWrappers::MPI::Vector tmp;
-    tmp.reinit (stokes_rhs.block(1));
-    pcout << "   Relative divergence residual: "
-          << stokes_matrix.block(1,0).residual (tmp,
-                                                distributed_stokes_solution.block(0),
-                                                stokes_rhs.block(1))
-          /
-          distributed_stokes_solution.block(0).l2_norm() / EquationData::pressure_scaling
-          << std::endl;
-
-    pcout << "   Relative vector sizes: "
-          << distributed_stokes_solution.block(0).linfty_norm() << ' '
-          << distributed_stokes_solution.block(1).linfty_norm() << std::endl;
-  }
-  computing_timer.exit_section();
+      TrilinosWrappers::MPI::BlockVector
+       distributed_stokes_solution (stokes_rhs);
+//    distributed_stokes_solution = stokes_solution;
+      distributed_stokes_solution.block(0).reinit(stokes_solution.block(0),false,true);
+      distributed_stokes_solution.block(1).reinit(stokes_solution.block(1),false,true);
+
+
+      const unsigned int
+       start = (distributed_stokes_solution.block(0).size() +
+                distributed_stokes_solution.block(1).local_range().first),
+       end   = (distributed_stokes_solution.block(0).size() +
+                distributed_stokes_solution.block(1).local_range().second);
+      for (unsigned int i=start; i<end; ++i)
+       if (stokes_constraints.is_constrained (i))
+         distributed_stokes_solution(i) = 0;
+
+
+      PrimitiveVectorMemory< TrilinosWrappers::MPI::BlockVector > mem;
+
+                                      // step 1: try if the simple and fast solver
+                                      // succeeds in 30 steps or less.
+      unsigned int n_iterations = 0;
+      double reduction = 0;
+      const double solver_tolerance = 1e-7 * stokes_rhs.l2_norm();
+      SolverControl solver_control (30, solver_tolerance);
+
+      try
+       {
+         const LinearSolvers::RightPrecond<TrilinosWrappers::PreconditionAMG,
+           TrilinosWrappers::PreconditionILU>
+           preconditioner (stokes_matrix, stokes_preconditioner_matrix,
+                           *Mp_preconditioner, *Amg_preconditioner,
+                           false);
+
+         SolverFGMRES<TrilinosWrappers::MPI::BlockVector>
+           solver(solver_control, mem,
+                  SolverFGMRES<TrilinosWrappers::MPI::BlockVector>::
+                  AdditionalData(30, true));
+         solver.solve(stokes_matrix, distributed_stokes_solution, stokes_rhs,
+                      preconditioner);
+
+         n_iterations = solver_control.last_step();
+         reduction = solver_control.last_value()/solver_control.initial_value();
+       }
+
+                                      // step 2: take the stronger solver in case
+                                      // the simple solver failed
+      catch (SolverControl::NoConvergence)
+       {
+         const LinearSolvers::RightPrecond<TrilinosWrappers::PreconditionAMG,
+           TrilinosWrappers::PreconditionILU>
+           preconditioner (stokes_matrix, stokes_preconditioner_matrix,
+                           *Mp_preconditioner, *Amg_preconditioner,
+                           true);
+
+         SolverControl solver_control_refined (stokes_matrix.m(), solver_tolerance);
+         SolverFGMRES<TrilinosWrappers::MPI::BlockVector>
+           solver(solver_control_refined, mem,
+                  SolverFGMRES<TrilinosWrappers::MPI::BlockVector>::
+                  AdditionalData(50, true));
+         solver.solve(stokes_matrix, distributed_stokes_solution, stokes_rhs,
+                      preconditioner);
+
+         n_iterations = (solver_control.last_step() +
+                         solver_control_refined.last_step());
+         reduction = (solver_control_refined.last_value()/
+                      std::max(solver_control.initial_value(),
+                               solver_control_refined.initial_value()));
+       }
+
+
+      stokes_constraints.distribute (distributed_stokes_solution);
+      stokes_solution.block(0).reinit(distributed_stokes_solution.block(0),
+                                     false, true);
+      stokes_solution.block(1).reinit(distributed_stokes_solution.block(1),
+                                     false, true);
+
+      pcout << n_iterations  << " iterations."
+           << " Reduced residual by " << reduction
+           << std::endl;
+
+      TrilinosWrappers::MPI::Vector tmp;
+      tmp.reinit (stokes_rhs.block(1));
+      pcout << "   Relative divergence residual: "
+           << stokes_matrix.block(1,0).residual (tmp,
+                                                 distributed_stokes_solution.block(0),
+                                                 stokes_rhs.block(1))
+       /
+       distributed_stokes_solution.block(0).l2_norm() / EquationData::pressure_scaling
+           << std::endl;
+
+      pcout << "   Relative vector sizes: "
+           << distributed_stokes_solution.block(0).linfty_norm() << ' '
+           << distributed_stokes_solution.block(1).linfty_norm() << std::endl;
+    }
+    computing_timer.exit_section();
 
 
-  computing_timer.enter_section ("   Assemble temperature rhs");
-  {
-    old_time_step = time_step;
-    const double cfl_number = get_cfl_number();
-
-    // we found out that we need
-    // approximately a quarter the time step
-    // size in 3d
-    double scaling = (dim==3)?0.25:1.0;
-    time_step = (scaling/(2.1*dim*std::sqrt(1.*dim)) /
-                 (parameters.temperature_degree *
-                  cfl_number));
-
-    const double maximal_velocity = get_maximal_velocity();
-    pcout << "   Maximal velocity: "
-          << maximal_velocity *EquationData::year_in_seconds * 100
-          << " cm/year"
-          << std::endl;
-    pcout << "   " << "Time step: "
-          << time_step/EquationData::year_in_seconds
-          << " years"
-          << std::endl;
-
-    temperature_solution = old_temperature_solution;
-    assemble_temperature_system (maximal_velocity);
-  }
-  computing_timer.exit_section ();
+    computing_timer.enter_section ("   Assemble temperature rhs");
+    {
+      old_time_step = time_step;
+      const double cfl_number = get_cfl_number();
+
+                                      // we found out that we need
+                                      // approximately a quarter the time step
+                                      // size in 3d
+      double scaling = (dim==3)?0.25:1.0;
+      time_step = (scaling/(2.1*dim*std::sqrt(1.*dim)) /
+                  (parameters.temperature_degree *
+                   cfl_number));
+
+      const double maximal_velocity = get_maximal_velocity();
+      pcout << "   Maximal velocity: "
+           << maximal_velocity *EquationData::year_in_seconds * 100
+           << " cm/year"
+           << std::endl;
+      pcout << "   " << "Time step: "
+           << time_step/EquationData::year_in_seconds
+           << " years"
+           << std::endl;
+
+      temperature_solution = old_temperature_solution;
+      assemble_temperature_system (maximal_velocity);
+    }
+    computing_timer.exit_section ();
 
-  computing_timer.enter_section ("   Solve temperature system");
-  {
-    SolverControl solver_control (temperature_matrix.m(),
-                                  1e-12*temperature_rhs.l2_norm());
-    SolverCG<TrilinosWrappers::MPI::Vector>   cg (solver_control);
+    computing_timer.enter_section ("   Solve temperature system");
+    {
+      SolverControl solver_control (temperature_matrix.m(),
+                                   1e-12*temperature_rhs.l2_norm());
+      SolverCG<TrilinosWrappers::MPI::Vector>   cg (solver_control);
 
-    TrilinosWrappers::MPI::Vector
-    distributed_temperature_solution (temperature_rhs);
+      TrilinosWrappers::MPI::Vector
+       distributed_temperature_solution (temperature_rhs);
 //    distributed_temperature_solution = temperature_solution;
-    distributed_temperature_solution.reinit(temperature_solution, false, true);
+      distributed_temperature_solution.reinit(temperature_solution, false, true);
 
-    cg.solve (temperature_matrix, distributed_temperature_solution,
-              temperature_rhs, *T_preconditioner);
+      cg.solve (temperature_matrix, distributed_temperature_solution,
+               temperature_rhs, *T_preconditioner);
 
-    temperature_constraints.distribute (distributed_temperature_solution);
+      temperature_constraints.distribute (distributed_temperature_solution);
 //    temperature_solution = distributed_temperature_solution;
-    temperature_solution.reinit(distributed_temperature_solution, false, true);
-
-    pcout << "   "
-          << solver_control.last_step()
-          << " CG iterations for temperature" << std::endl;
-    computing_timer.exit_section();
-
-    // extract temperature range
-    std::vector<double> temperature (2), global_temperature (2);
-    temperature[0] = std::numeric_limits<double>::max(),
-                     temperature[1] = -std::numeric_limits<double>::max();
-    for (unsigned int i=0; i<distributed_temperature_solution.local_size(); ++i)
-      {
-        temperature[0] = std::min<double> (temperature[0],
-                                           distributed_temperature_solution.trilinos_vector()[0][i]);
-        temperature[1] = std::max<double> (temperature[1],
-                                           distributed_temperature_solution.trilinos_vector()[0][i]);
-      }
+      temperature_solution.reinit(distributed_temperature_solution, false, true);
+
+      pcout << "   "
+           << solver_control.last_step()
+           << " CG iterations for temperature" << std::endl;
+      computing_timer.exit_section();
+
+                                      // extract temperature range
+      std::vector<double> temperature (2), global_temperature (2);
+      temperature[0] = std::numeric_limits<double>::max(),
+      temperature[1] = -std::numeric_limits<double>::max();
+      for (unsigned int i=0; i<distributed_temperature_solution.local_size(); ++i)
+       {
+         temperature[0] = std::min<double> (temperature[0],
+                                            distributed_temperature_solution.trilinos_vector()[0][i]);
+         temperature[1] = std::max<double> (temperature[1],
+                                            distributed_temperature_solution.trilinos_vector()[0][i]);
+       }
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
-    temperature[0] *= -1.0;
-    MPI_Allreduce (&temperature[0], &global_temperature[0],
-                   2, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
-    global_temperature[0] *= -1.0;
+      temperature[0] *= -1.0;
+      MPI_Allreduce (&temperature[0], &global_temperature[0],
+                    2, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+      global_temperature[0] *= -1.0;
 #else
-    global_temperature = local_temperature;
+      global_temperature = local_temperature;
 #endif
 
-    pcout << "   Temperature range: "
-          << global_temperature[0] << ' ' << global_temperature[1]
-          << std::endl;
+      pcout << "   Temperature range: "
+           << global_temperature[0] << ' ' << global_temperature[1]
+           << std::endl;
+    }
   }
-}
 
 
 // @sect4{BoussinesqFlowProblem::output_results}
 
-template <int dim>
-class BoussinesqFlowProblem<dim>::Postprocessor : public DataPostprocessor<dim>
-{
-  public:
-    Postprocessor (const unsigned int partition,
-                   const double       minimal_pressure);
+  template <int dim>
+  class BoussinesqFlowProblem<dim>::Postprocessor : public DataPostprocessor<dim>
+  {
+    public:
+      Postprocessor (const unsigned int partition,
+                    const double       minimal_pressure);
 
-    virtual
-    void
-    compute_derived_quantities_vector (const std::vector<Vector<double> >              &uh,
-                                       const std::vector<std::vector<Tensor<1,dim> > > &duh,
-                                       const std::vector<std::vector<Tensor<2,dim> > > &dduh,
-                                       const std::vector<Point<dim> >                  &normals,
-                                       const std::vector<Point<dim> >                  &evaluation_points,
-                                       std::vector<Vector<double> >                    &computed_quantities) const;
+      virtual
+      void
+      compute_derived_quantities_vector (const std::vector<Vector<double> >              &uh,
+                                        const std::vector<std::vector<Tensor<1,dim> > > &duh,
+                                        const std::vector<std::vector<Tensor<2,dim> > > &dduh,
+                                        const std::vector<Point<dim> >                  &normals,
+                                        const std::vector<Point<dim> >                  &evaluation_points,
+                                        std::vector<Vector<double> >                    &computed_quantities) const;
 
-    virtual std::vector<std::string> get_names () const;
+      virtual std::vector<std::string> get_names () const;
 
-    virtual unsigned int n_output_variables() const;
+      virtual unsigned int n_output_variables() const;
 
-    virtual
-    std::vector<DataComponentInterpretation::DataComponentInterpretation>
-    get_data_component_interpretation () const;
+      virtual
+      std::vector<DataComponentInterpretation::DataComponentInterpretation>
+      get_data_component_interpretation () const;
 
-    virtual UpdateFlags get_needed_update_flags () const;
+      virtual UpdateFlags get_needed_update_flags () const;
 
-  private:
-    const unsigned int partition;
-    const double       minimal_pressure;
-};
+    private:
+      const unsigned int partition;
+      const double       minimal_pressure;
+  };
 
 
-template <int dim>
-BoussinesqFlowProblem<dim>::Postprocessor::
-Postprocessor (const unsigned int partition,
-               const double       minimal_pressure)
-  :
-  partition (partition),
-  minimal_pressure (minimal_pressure)
-{}
+  template <int dim>
+  BoussinesqFlowProblem<dim>::Postprocessor::
+  Postprocessor (const unsigned int partition,
+                const double       minimal_pressure)
+                 :
+                 partition (partition),
+                 minimal_pressure (minimal_pressure)
+  {}
 
 
-template <int dim>
-std::vector<std::string>
-BoussinesqFlowProblem<dim>::Postprocessor::get_names() const
-{
-  std::vector<std::string> solution_names (dim, "velocity");
-  solution_names.push_back ("p");
-  solution_names.push_back ("T");
-  solution_names.push_back ("friction_heating");
-  solution_names.push_back ("partition");
-  solution_names.push_back ("non_adiabatic_pressure");
-  solution_names.push_back ("non_adiabatic_temperature");
-  return solution_names;
-}
+  template <int dim>
+  std::vector<std::string>
+  BoussinesqFlowProblem<dim>::Postprocessor::get_names() const
+  {
+    std::vector<std::string> solution_names (dim, "velocity");
+    solution_names.push_back ("p");
+    solution_names.push_back ("T");
+    solution_names.push_back ("friction_heating");
+    solution_names.push_back ("partition");
+    solution_names.push_back ("non_adiabatic_pressure");
+    solution_names.push_back ("non_adiabatic_temperature");
+    return solution_names;
+  }
 
 
-template <int dim>
-unsigned int
-BoussinesqFlowProblem<dim>::Postprocessor::n_output_variables() const
-{
-  return get_names().size();
-}
+  template <int dim>
+  unsigned int
+  BoussinesqFlowProblem<dim>::Postprocessor::n_output_variables() const
+  {
+    return get_names().size();
+  }
 
 
-template <int dim>
-std::vector<DataComponentInterpretation::DataComponentInterpretation>
-BoussinesqFlowProblem<dim>::Postprocessor::
-get_data_component_interpretation () const
-{
+  template <int dim>
   std::vector<DataComponentInterpretation::DataComponentInterpretation>
-  interpretation (dim,
-                  DataComponentInterpretation::component_is_part_of_vector);
+  BoussinesqFlowProblem<dim>::Postprocessor::
+  get_data_component_interpretation () const
+  {
+    std::vector<DataComponentInterpretation::DataComponentInterpretation>
+      interpretation (dim,
+                     DataComponentInterpretation::component_is_part_of_vector);
 
-  interpretation.push_back (DataComponentInterpretation::component_is_scalar);
-  interpretation.push_back (DataComponentInterpretation::component_is_scalar);
-  interpretation.push_back (DataComponentInterpretation::component_is_scalar);
-  interpretation.push_back (DataComponentInterpretation::component_is_scalar);
-  interpretation.push_back (DataComponentInterpretation::component_is_scalar);
-  interpretation.push_back (DataComponentInterpretation::component_is_scalar);
+    interpretation.push_back (DataComponentInterpretation::component_is_scalar);
+    interpretation.push_back (DataComponentInterpretation::component_is_scalar);
+    interpretation.push_back (DataComponentInterpretation::component_is_scalar);
+    interpretation.push_back (DataComponentInterpretation::component_is_scalar);
+    interpretation.push_back (DataComponentInterpretation::component_is_scalar);
+    interpretation.push_back (DataComponentInterpretation::component_is_scalar);
 
-  return interpretation;
-}
+    return interpretation;
+  }
 
 
-template <int dim>
-UpdateFlags
-BoussinesqFlowProblem<dim>::Postprocessor::get_needed_update_flags() const
-{
-  return update_values | update_gradients | update_q_points;
-}
+  template <int dim>
+  UpdateFlags
+  BoussinesqFlowProblem<dim>::Postprocessor::get_needed_update_flags() const
+  {
+    return update_values | update_gradients | update_q_points;
+  }
 
 
-template <int dim>
-void
-BoussinesqFlowProblem<dim>::Postprocessor::
-compute_derived_quantities_vector (const std::vector<Vector<double> >              &uh,
-                                   const std::vector<std::vector<Tensor<1,dim> > > &duh,
-                                   const std::vector<std::vector<Tensor<2,dim> > > &/*dduh*/,
-                                   const std::vector<Point<dim> >                  &/*normals*/,
-                                   const std::vector<Point<dim> >                  &evaluation_points,
-                                   std::vector<Vector<double> >                    &computed_quantities) const
-{
-  const unsigned int n_quadrature_points = uh.size();
-  Assert (duh.size() == n_quadrature_points,                  ExcInternalError());
-  Assert (computed_quantities.size() == n_quadrature_points,  ExcInternalError());
-  Assert (uh[0].size() == dim+2,                              ExcInternalError());
-  Assert (computed_quantities[0].size()==n_output_variables(),ExcInternalError());
+  template <int dim>
+  void
+  BoussinesqFlowProblem<dim>::Postprocessor::
+  compute_derived_quantities_vector (const std::vector<Vector<double> >              &uh,
+                                    const std::vector<std::vector<Tensor<1,dim> > > &duh,
+                                    const std::vector<std::vector<Tensor<2,dim> > > &/*dduh*/,
+                                    const std::vector<Point<dim> >                  &/*normals*/,
+                                    const std::vector<Point<dim> >                  &evaluation_points,
+                                    std::vector<Vector<double> >                    &computed_quantities) const
+  {
+    const unsigned int n_quadrature_points = uh.size();
+    Assert (duh.size() == n_quadrature_points,                  ExcInternalError());
+    Assert (computed_quantities.size() == n_quadrature_points,  ExcInternalError());
+    Assert (uh[0].size() == dim+2,                              ExcInternalError());
+    Assert (computed_quantities[0].size()==n_output_variables(),ExcInternalError());
 
-  for (unsigned int q=0; q<n_quadrature_points; ++q)
-    {
-      // velocity; rescale in cm/year
-      for (unsigned int d=0; d<dim; ++d)
-        computed_quantities[q](d)
-          = (uh[q](d) *  EquationData::year_in_seconds * 100);
-
-      // pressure
-      const double pressure = (uh[q](dim)-minimal_pressure) * EquationData::pressure_scaling;
-      computed_quantities[q](dim) = pressure;
-
-      // temperature
-      const double temperature = uh[q](dim+1);
-      computed_quantities[q](dim+1) = temperature;
-
-      // friction heating
-      Tensor<2,dim> grad_u;
-      for (unsigned int d=0; d<dim; ++d)
-        grad_u[d] = duh[q][d];
-      const SymmetricTensor<2,dim> strain_rate = symmetrize (grad_u);
-      computed_quantities[q](dim+2) = 2 * EquationData::eta *
-                                      strain_rate * strain_rate;
-
-      computed_quantities[q](dim+3) = partition;
-
-      computed_quantities[q](dim+4) = pressure -
-                                      EquationData::adiabatic_pressure (evaluation_points[q]);
-
-      computed_quantities[q](dim+5) = temperature -
-                                      EquationData::adiabatic_temperature (evaluation_points[q]);
-    }
-}
+    for (unsigned int q=0; q<n_quadrature_points; ++q)
+      {
+                                        // velocity; rescale in cm/year
+       for (unsigned int d=0; d<dim; ++d)
+         computed_quantities[q](d)
+           = (uh[q](d) *  EquationData::year_in_seconds * 100);
+
+                                        // pressure
+       const double pressure = (uh[q](dim)-minimal_pressure) * EquationData::pressure_scaling;
+       computed_quantities[q](dim) = pressure;
+
+                                        // temperature
+       const double temperature = uh[q](dim+1);
+       computed_quantities[q](dim+1) = temperature;
+
+                                        // friction heating
+       Tensor<2,dim> grad_u;
+       for (unsigned int d=0; d<dim; ++d)
+         grad_u[d] = duh[q][d];
+       const SymmetricTensor<2,dim> strain_rate = symmetrize (grad_u);
+       computed_quantities[q](dim+2) = 2 * EquationData::eta *
+                                       strain_rate * strain_rate;
+
+       computed_quantities[q](dim+3) = partition;
+
+       computed_quantities[q](dim+4) = pressure -
+                                       EquationData::adiabatic_pressure (evaluation_points[q]);
+
+       computed_quantities[q](dim+5) = temperature -
+                                       EquationData::adiabatic_temperature (evaluation_points[q]);
+      }
+  }
 
 
 // This function does mostly what the
@@ -3545,121 +3547,121 @@ compute_derived_quantities_vector (const std::vector<Vector<double> >
 // therefore rescale things into
 // centimeters per year, the unit
 // commonly used in geophysics.
-template <int dim>
-void BoussinesqFlowProblem<dim>::output_results ()
-{
-  computing_timer.enter_section ("Postprocessing");
+  template <int dim>
+  void BoussinesqFlowProblem<dim>::output_results ()
+  {
+    computing_timer.enter_section ("Postprocessing");
 
-  const FESystem<dim> joint_fe (stokes_fe, 1,
-                                temperature_fe, 1);
+    const FESystem<dim> joint_fe (stokes_fe, 1,
+                                 temperature_fe, 1);
 
-  DoFHandler<dim> joint_dof_handler (triangulation);
-  joint_dof_handler.distribute_dofs (joint_fe);
-  Assert (joint_dof_handler.n_dofs() ==
-          stokes_dof_handler.n_dofs() + temperature_dof_handler.n_dofs(),
-          ExcInternalError());
+    DoFHandler<dim> joint_dof_handler (triangulation);
+    joint_dof_handler.distribute_dofs (joint_fe);
+    Assert (joint_dof_handler.n_dofs() ==
+           stokes_dof_handler.n_dofs() + temperature_dof_handler.n_dofs(),
+           ExcInternalError());
 
-  TrilinosWrappers::MPI::Vector joint_solution;
-  joint_solution.reinit (joint_dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
+    TrilinosWrappers::MPI::Vector joint_solution;
+    joint_solution.reinit (joint_dof_handler.locally_owned_dofs(), MPI_COMM_WORLD);
 
-  {
-    std::vector<unsigned int> local_joint_dof_indices (joint_fe.dofs_per_cell);
-    std::vector<unsigned int> local_stokes_dof_indices (stokes_fe.dofs_per_cell);
-    std::vector<unsigned int> local_temperature_dof_indices (temperature_fe.dofs_per_cell);
+    {
+      std::vector<unsigned int> local_joint_dof_indices (joint_fe.dofs_per_cell);
+      std::vector<unsigned int> local_stokes_dof_indices (stokes_fe.dofs_per_cell);
+      std::vector<unsigned int> local_temperature_dof_indices (temperature_fe.dofs_per_cell);
 
-    typename DoFHandler<dim>::active_cell_iterator
-    joint_cell       = joint_dof_handler.begin_active(),
-    joint_endc       = joint_dof_handler.end(),
-    stokes_cell      = stokes_dof_handler.begin_active(),
-    temperature_cell = temperature_dof_handler.begin_active();
-    for (; joint_cell!=joint_endc;
-         ++joint_cell, ++stokes_cell, ++temperature_cell)
-      if (joint_cell->is_locally_owned())
-        {
-          joint_cell->get_dof_indices (local_joint_dof_indices);
-          stokes_cell->get_dof_indices (local_stokes_dof_indices);
-          temperature_cell->get_dof_indices (local_temperature_dof_indices);
-
-          for (unsigned int i=0; i<joint_fe.dofs_per_cell; ++i)
-            if (joint_fe.system_to_base_index(i).first.first == 0)
-              {
-                Assert (joint_fe.system_to_base_index(i).second
-                        <
-                        local_stokes_dof_indices.size(),
-                        ExcInternalError());
-
-                joint_solution(local_joint_dof_indices[i])
-                  = stokes_solution(local_stokes_dof_indices
-                                    [joint_fe.system_to_base_index(i).second]);
-              }
-            else
-              {
-                Assert (joint_fe.system_to_base_index(i).first.first == 1,
-                        ExcInternalError());
-                Assert (joint_fe.system_to_base_index(i).second
-                        <
-                        local_temperature_dof_indices.size(),
-                        ExcInternalError());
-                joint_solution(local_joint_dof_indices[i])
-                  = temperature_solution(local_temperature_dof_indices
-                                         [joint_fe.system_to_base_index(i).second]);
-              }
-        }
-  }
+      typename DoFHandler<dim>::active_cell_iterator
+       joint_cell       = joint_dof_handler.begin_active(),
+       joint_endc       = joint_dof_handler.end(),
+       stokes_cell      = stokes_dof_handler.begin_active(),
+       temperature_cell = temperature_dof_handler.begin_active();
+      for (; joint_cell!=joint_endc;
+          ++joint_cell, ++stokes_cell, ++temperature_cell)
+       if (joint_cell->is_locally_owned())
+         {
+           joint_cell->get_dof_indices (local_joint_dof_indices);
+           stokes_cell->get_dof_indices (local_stokes_dof_indices);
+           temperature_cell->get_dof_indices (local_temperature_dof_indices);
+
+           for (unsigned int i=0; i<joint_fe.dofs_per_cell; ++i)
+             if (joint_fe.system_to_base_index(i).first.first == 0)
+               {
+                 Assert (joint_fe.system_to_base_index(i).second
+                         <
+                         local_stokes_dof_indices.size(),
+                         ExcInternalError());
+
+                 joint_solution(local_joint_dof_indices[i])
+                   = stokes_solution(local_stokes_dof_indices
+                                     [joint_fe.system_to_base_index(i).second]);
+               }
+             else
+               {
+                 Assert (joint_fe.system_to_base_index(i).first.first == 1,
+                         ExcInternalError());
+                 Assert (joint_fe.system_to_base_index(i).second
+                         <
+                         local_temperature_dof_indices.size(),
+                         ExcInternalError());
+                 joint_solution(local_joint_dof_indices[i])
+                   = temperature_solution(local_temperature_dof_indices
+                                          [joint_fe.system_to_base_index(i).second]);
+               }
+         }
+    }
 
 
-  IndexSet locally_relevant_joint_dofs(joint_dof_handler.n_dofs());
-  DoFTools::extract_locally_relevant_dofs (joint_dof_handler, locally_relevant_joint_dofs);
-  TrilinosWrappers::MPI::Vector locally_relevant_joint_solution;
-  locally_relevant_joint_solution.reinit (locally_relevant_joint_dofs, MPI_COMM_WORLD);
-  locally_relevant_joint_solution = joint_solution;
-
-  Postprocessor postprocessor (Utilities::System::
-                               get_this_mpi_process(MPI_COMM_WORLD),
-                               stokes_solution.block(1).minimal_value());
-
-  DataOut<dim> data_out;
-  data_out.attach_dof_handler (joint_dof_handler);
-  data_out.add_data_vector (locally_relevant_joint_solution, postprocessor);
-  data_out.build_patches ();
-
-  static int out_index=0;
-  const std::string filename = ("solution-" +
-                                Utilities::int_to_string (out_index, 5) +
-                                "." +
-                                Utilities::int_to_string
-                                (triangulation.locally_owned_subdomain(), 4) +
-                                ".vtu");
-  std::ofstream output (filename.c_str());
-  data_out.write_vtu (output);
-
-  if (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD) == 0)
-    {
-      std::vector<std::string> filenames;
-      for (unsigned int i=0; i<Utilities::System::get_n_mpi_processes(MPI_COMM_WORLD); ++i)
-        filenames.push_back (std::string("solution-") +
-                             Utilities::int_to_string (out_index, 5) +
-                             "." +
-                             Utilities::int_to_string(i, 4) +
-                             ".vtu");
-      const std::string
-      pvtu_master_filename = ("solution-" +
-                              Utilities::int_to_string (out_index, 5) +
-                              ".pvtu");
-      std::ofstream pvtu_master (pvtu_master_filename.c_str());
-      data_out.write_pvtu_record (pvtu_master, filenames);
-
-      const std::string
-      visit_master_filename = ("solution-" +
-                               Utilities::int_to_string (out_index, 5) +
-                               ".visit");
-      std::ofstream visit_master (visit_master_filename.c_str());
-      data_out.write_visit_record (visit_master, filenames);
-    }
+    IndexSet locally_relevant_joint_dofs(joint_dof_handler.n_dofs());
+    DoFTools::extract_locally_relevant_dofs (joint_dof_handler, locally_relevant_joint_dofs);
+    TrilinosWrappers::MPI::Vector locally_relevant_joint_solution;
+    locally_relevant_joint_solution.reinit (locally_relevant_joint_dofs, MPI_COMM_WORLD);
+    locally_relevant_joint_solution = joint_solution;
+
+    Postprocessor postprocessor (Utilities::System::
+                                get_this_mpi_process(MPI_COMM_WORLD),
+                                stokes_solution.block(1).minimal_value());
+
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler (joint_dof_handler);
+    data_out.add_data_vector (locally_relevant_joint_solution, postprocessor);
+    data_out.build_patches ();
+
+    static int out_index=0;
+    const std::string filename = ("solution-" +
+                                 Utilities::int_to_string (out_index, 5) +
+                                 "." +
+                                 Utilities::int_to_string
+                                 (triangulation.locally_owned_subdomain(), 4) +
+                                 ".vtu");
+    std::ofstream output (filename.c_str());
+    data_out.write_vtu (output);
+
+    if (Utilities::System::get_this_mpi_process(MPI_COMM_WORLD) == 0)
+      {
+       std::vector<std::string> filenames;
+       for (unsigned int i=0; i<Utilities::System::get_n_mpi_processes(MPI_COMM_WORLD); ++i)
+         filenames.push_back (std::string("solution-") +
+                              Utilities::int_to_string (out_index, 5) +
+                              "." +
+                              Utilities::int_to_string(i, 4) +
+                              ".vtu");
+       const std::string
+         pvtu_master_filename = ("solution-" +
+                                 Utilities::int_to_string (out_index, 5) +
+                                 ".pvtu");
+       std::ofstream pvtu_master (pvtu_master_filename.c_str());
+       data_out.write_pvtu_record (pvtu_master, filenames);
+
+       const std::string
+         visit_master_filename = ("solution-" +
+                                  Utilities::int_to_string (out_index, 5) +
+                                  ".visit");
+       std::ofstream visit_master (visit_master_filename.c_str());
+       data_out.write_visit_record (visit_master, filenames);
+      }
 
-  computing_timer.exit_section ();
-  out_index++;
-}
+    computing_timer.exit_section ();
+    out_index++;
+  }
 
 
 
@@ -3753,93 +3755,93 @@ void BoussinesqFlowProblem<dim>::output_results ()
 // complete set of refinement
 // indicators, and the rest of the
 // function proceeds as before.
-template <int dim>
-void BoussinesqFlowProblem<dim>::refine_mesh (const unsigned int max_grid_level)
-{
-  computing_timer.enter_section ("Refine mesh structure, part 1");
-  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
-
-  KellyErrorEstimator<dim>::estimate (temperature_dof_handler,
-                                      QGauss<dim-1>(parameters.temperature_degree+1),
-                                      typename FunctionMap<dim>::type(),
-                                      temperature_solution,
-                                      estimated_error_per_cell,
-                                      std::vector<bool>(),
-                                      0,
-                                      0,
-                                      triangulation.locally_owned_subdomain());
-
-  parallel::distributed::GridRefinement::
-  refine_and_coarsen_fixed_fraction (triangulation,
-                                     estimated_error_per_cell,
-                                     0.3, 0.1);
-
-  // limit maximum refinement level
-  if (triangulation.n_levels() > max_grid_level)
-    for (typename Triangulation<dim>::active_cell_iterator
-         cell = triangulation.begin_active(max_grid_level);
-         cell != triangulation.end(); ++cell)
-      cell->clear_refine_flag ();
-
-  std::vector<const TrilinosWrappers::MPI::Vector *> x_temperature (2);
-  x_temperature[0] = &temperature_solution;
-  x_temperature[1] = &old_temperature_solution;
-  std::vector<const TrilinosWrappers::MPI::BlockVector *> x_stokes (2);
-  x_stokes[0] = &stokes_solution;
-  x_stokes[1] = &old_stokes_solution;
-
-  parallel::distributed::SolutionTransfer<dim,TrilinosWrappers::MPI::Vector>
-  temperature_trans(temperature_dof_handler);
-  parallel::distributed::SolutionTransfer<dim,TrilinosWrappers::MPI::BlockVector>
-  stokes_trans(stokes_dof_handler);
-
-  triangulation.prepare_coarsening_and_refinement();
-  temperature_trans.prepare_for_coarsening_and_refinement(x_temperature);
-  stokes_trans.prepare_for_coarsening_and_refinement(x_stokes);
-
-  triangulation.execute_coarsening_and_refinement ();
-  computing_timer.exit_section();
-
-  setup_dofs ();
-
-  computing_timer.enter_section ("Refine mesh structure, part 2");
-
+  template <int dim>
+  void BoussinesqFlowProblem<dim>::refine_mesh (const unsigned int max_grid_level)
   {
-    TrilinosWrappers::MPI::Vector
-    distributed_temp1 (temperature_rhs);
-    TrilinosWrappers::MPI::Vector
-    distributed_temp2 (temperature_rhs);
+    computing_timer.enter_section ("Refine mesh structure, part 1");
+    Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+    KellyErrorEstimator<dim>::estimate (temperature_dof_handler,
+                                       QGauss<dim-1>(parameters.temperature_degree+1),
+                                       typename FunctionMap<dim>::type(),
+                                       temperature_solution,
+                                       estimated_error_per_cell,
+                                       std::vector<bool>(),
+                                       0,
+                                       0,
+                                       triangulation.locally_owned_subdomain());
+
+    parallel::distributed::GridRefinement::
+      refine_and_coarsen_fixed_fraction (triangulation,
+                                        estimated_error_per_cell,
+                                        0.3, 0.1);
+
+                                    // limit maximum refinement level
+    if (triangulation.n_levels() > max_grid_level)
+      for (typename Triangulation<dim>::active_cell_iterator
+            cell = triangulation.begin_active(max_grid_level);
+          cell != triangulation.end(); ++cell)
+       cell->clear_refine_flag ();
+
+    std::vector<const TrilinosWrappers::MPI::Vector *> x_temperature (2);
+    x_temperature[0] = &temperature_solution;
+    x_temperature[1] = &old_temperature_solution;
+    std::vector<const TrilinosWrappers::MPI::BlockVector *> x_stokes (2);
+    x_stokes[0] = &stokes_solution;
+    x_stokes[1] = &old_stokes_solution;
+
+    parallel::distributed::SolutionTransfer<dim,TrilinosWrappers::MPI::Vector>
+      temperature_trans(temperature_dof_handler);
+    parallel::distributed::SolutionTransfer<dim,TrilinosWrappers::MPI::BlockVector>
+      stokes_trans(stokes_dof_handler);
+
+    triangulation.prepare_coarsening_and_refinement();
+    temperature_trans.prepare_for_coarsening_and_refinement(x_temperature);
+    stokes_trans.prepare_for_coarsening_and_refinement(x_stokes);
+
+    triangulation.execute_coarsening_and_refinement ();
+    computing_timer.exit_section();
 
-    std::vector<TrilinosWrappers::MPI::Vector *> tmp (2);
-    tmp[0] = &(distributed_temp1);
-    tmp[1] = &(distributed_temp2);
-    temperature_trans.interpolate(tmp);
+    setup_dofs ();
 
-    //  temperature_solution = distributed_temp1;
-    temperature_solution.reinit(distributed_temp1, false, true);
-    //  old_temperature_solution = distributed_temp2;
-    old_temperature_solution.reinit(distributed_temp2, false, true);
-  }
+    computing_timer.enter_section ("Refine mesh structure, part 2");
 
-  {
-    TrilinosWrappers::MPI::BlockVector
-    distributed_stokes (stokes_rhs);
-    TrilinosWrappers::MPI::BlockVector
-    old_distributed_stokes (stokes_rhs);
-    std::vector<TrilinosWrappers::MPI::BlockVector *> stokes_tmp (2);
-    stokes_tmp[0] = &(distributed_stokes);
-    stokes_tmp[1] = &(old_distributed_stokes);
-
-    stokes_trans.interpolate (stokes_tmp);
-    //  stokes_solution = distributed_stokes;
-    stokes_solution.block(0).reinit(distributed_stokes.block(0), false, true);
-    stokes_solution.block(1).reinit(distributed_stokes.block(1), false, true);
-    old_stokes_solution.block(0).reinit(old_distributed_stokes.block(0), false, true);
-    old_stokes_solution.block(1).reinit(old_distributed_stokes.block(1), false, true);
-  }
+    {
+      TrilinosWrappers::MPI::Vector
+       distributed_temp1 (temperature_rhs);
+      TrilinosWrappers::MPI::Vector
+       distributed_temp2 (temperature_rhs);
+
+      std::vector<TrilinosWrappers::MPI::Vector *> tmp (2);
+      tmp[0] = &(distributed_temp1);
+      tmp[1] = &(distributed_temp2);
+      temperature_trans.interpolate(tmp);
+
+                                      //  temperature_solution = distributed_temp1;
+      temperature_solution.reinit(distributed_temp1, false, true);
+                                      //  old_temperature_solution = distributed_temp2;
+      old_temperature_solution.reinit(distributed_temp2, false, true);
+    }
 
-  computing_timer.exit_section();
-}
+    {
+      TrilinosWrappers::MPI::BlockVector
+       distributed_stokes (stokes_rhs);
+      TrilinosWrappers::MPI::BlockVector
+       old_distributed_stokes (stokes_rhs);
+      std::vector<TrilinosWrappers::MPI::BlockVector *> stokes_tmp (2);
+      stokes_tmp[0] = &(distributed_stokes);
+      stokes_tmp[1] = &(old_distributed_stokes);
+
+      stokes_trans.interpolate (stokes_tmp);
+                                      //  stokes_solution = distributed_stokes;
+      stokes_solution.block(0).reinit(distributed_stokes.block(0), false, true);
+      stokes_solution.block(1).reinit(distributed_stokes.block(1), false, true);
+      old_stokes_solution.block(0).reinit(old_distributed_stokes.block(0), false, true);
+      old_stokes_solution.block(1).reinit(old_distributed_stokes.block(1), false, true);
+    }
+
+    computing_timer.exit_section();
+  }
 
 
 
@@ -3854,111 +3856,112 @@ void BoussinesqFlowProblem<dim>::refine_mesh (const unsigned int max_grid_level)
 // function instead of the library function
 // <code>VectorTools::project</code>, the
 // rest is as before.
-template <int dim>
-void BoussinesqFlowProblem<dim>::run ()
-{
-  GridGenerator::hyper_shell (triangulation,
-                              Point<dim>(),
-                              EquationData::R0,
-                              EquationData::R1,
-                              (dim==3) ? 96 : 12,
-                              true);
-  static HyperShellBoundary<dim> boundary;
-  triangulation.set_boundary (0, boundary);
-  triangulation.set_boundary (1, boundary);
+  template <int dim>
+  void BoussinesqFlowProblem<dim>::run ()
+  {
+    GridGenerator::hyper_shell (triangulation,
+                               Point<dim>(),
+                               EquationData::R0,
+                               EquationData::R1,
+                               (dim==3) ? 96 : 12,
+                               true);
+    static HyperShellBoundary<dim> boundary;
+    triangulation.set_boundary (0, boundary);
+    triangulation.set_boundary (1, boundary);
 
-  global_Omega_diameter = GridTools::diameter (triangulation);
+    global_Omega_diameter = GridTools::diameter (triangulation);
 
-  triangulation.refine_global (parameters.initial_global_refinement);
+    triangulation.refine_global (parameters.initial_global_refinement);
 
-  setup_dofs();
+    setup_dofs();
 
-  unsigned int pre_refinement_step = 0;
+    unsigned int pre_refinement_step = 0;
 
-start_time_iteration:
+    start_time_iteration:
 
-  project_temperature_field ();
+    project_temperature_field ();
 
-  timestep_number           = 0;
-  time_step = old_time_step = 0;
+    timestep_number           = 0;
+    time_step = old_time_step = 0;
 
-  double time = 0;
+    double time = 0;
 
-  do
-    {
-      pcout << "Timestep " << timestep_number
-            << ":  t=" << time/EquationData::year_in_seconds
-            << " years"
-            << std::endl;
-
-      assemble_stokes_system ();
-      build_stokes_preconditioner ();
-      assemble_temperature_matrix ();
-
-      solve ();
-
-      pcout << std::endl;
-
-      if ((timestep_number == 0) &&
-          (pre_refinement_step < parameters.initial_adaptive_refinement))
-        {
-          refine_mesh (parameters.initial_global_refinement +
-                       parameters.initial_adaptive_refinement);
-          ++pre_refinement_step;
-          goto start_time_iteration;
-        }
-      else if ((timestep_number > 0)
-               &&
-               (timestep_number % parameters.adaptive_refinement_interval == 0))
-        refine_mesh (parameters.initial_global_refinement +
-                     parameters.initial_adaptive_refinement);
-
-      if ((parameters.generate_graphical_output == true)
-          &&
-          (timestep_number % parameters.graphical_output_interval == 0))
-        output_results ();
-
-      time += time_step;
-      ++timestep_number;
-
-      // if we are at the end of
-      // time, stop now
-      if (time > parameters.end_time * EquationData::year_in_seconds)
-        break;
-
-      // otherwise prepare for the
-      // next time step
-      TrilinosWrappers::MPI::BlockVector old_old_stokes_solution;
-      old_old_stokes_solution      = old_stokes_solution;
-      old_stokes_solution          = stokes_solution;
-      old_old_temperature_solution = old_temperature_solution;
-      old_temperature_solution     = temperature_solution;
-      if (old_time_step > 0)
-        {
-          stokes_solution.sadd (1.+time_step/old_time_step, -time_step/old_time_step,
-                                old_old_stokes_solution);
-          temperature_solution.sadd (1.+time_step/old_time_step,
-                                     -time_step/old_time_step,
-                                     old_old_temperature_solution);
-        }
-
-      // every 100 time steps output
-      // a summary of the current
-      // timing information
-      if (timestep_number % 100 == 0)
-        computing_timer.print_summary ();
-    }
-  while (true);
-
-  // if we are generating graphical
-  // output, do so also for the last
-  // time step unless we had just
-  // done so before we left the
-  // do-while loop
-  if ((parameters.generate_graphical_output == true)
-      &&
-      !((timestep_number-1) % parameters.graphical_output_interval == 0))
-    output_results ();
+    do
+      {
+       pcout << "Timestep " << timestep_number
+             << ":  t=" << time/EquationData::year_in_seconds
+             << " years"
+             << std::endl;
+
+       assemble_stokes_system ();
+       build_stokes_preconditioner ();
+       assemble_temperature_matrix ();
+
+       solve ();
+
+       pcout << std::endl;
+
+       if ((timestep_number == 0) &&
+           (pre_refinement_step < parameters.initial_adaptive_refinement))
+         {
+           refine_mesh (parameters.initial_global_refinement +
+                        parameters.initial_adaptive_refinement);
+           ++pre_refinement_step;
+           goto start_time_iteration;
+         }
+       else if ((timestep_number > 0)
+                &&
+                (timestep_number % parameters.adaptive_refinement_interval == 0))
+         refine_mesh (parameters.initial_global_refinement +
+                      parameters.initial_adaptive_refinement);
+
+       if ((parameters.generate_graphical_output == true)
+           &&
+           (timestep_number % parameters.graphical_output_interval == 0))
+         output_results ();
+
+       time += time_step;
+       ++timestep_number;
+
+                                        // if we are at the end of
+                                        // time, stop now
+       if (time > parameters.end_time * EquationData::year_in_seconds)
+         break;
+
+                                        // otherwise prepare for the
+                                        // next time step
+       TrilinosWrappers::MPI::BlockVector old_old_stokes_solution;
+       old_old_stokes_solution      = old_stokes_solution;
+       old_stokes_solution          = stokes_solution;
+       old_old_temperature_solution = old_temperature_solution;
+       old_temperature_solution     = temperature_solution;
+       if (old_time_step > 0)
+         {
+           stokes_solution.sadd (1.+time_step/old_time_step, -time_step/old_time_step,
+                                 old_old_stokes_solution);
+           temperature_solution.sadd (1.+time_step/old_time_step,
+                                      -time_step/old_time_step,
+                                      old_old_temperature_solution);
+         }
+
+                                        // every 100 time steps output
+                                        // a summary of the current
+                                        // timing information
+       if (timestep_number % 100 == 0)
+         computing_timer.print_summary ();
+      }
+    while (true);
+
+                                    // if we are generating graphical
+                                    // output, do so also for the last
+                                    // time step unless we had just
+                                    // done so before we left the
+                                    // do-while loop
+    if ((parameters.generate_graphical_output == true)
+       &&
+       !((timestep_number-1) % parameters.graphical_output_interval == 0))
+      output_results ();
+  }
 }
 
 
@@ -3968,6 +3971,9 @@ start_time_iteration:
 // This is copied verbatim from step-31:
 int main (int argc, char *argv[])
 {
+  using namespace Step32;
+  using namespace dealii;
+
   Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv);
 
   try

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.