]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Remove namespace step-51 as the class name already contains step-51. Speed up assembl...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 6 Aug 2013 17:03:45 +0000 (17:03 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 6 Aug 2013 17:03:45 +0000 (17:03 +0000)
git-svn-id: https://svn.dealii.org/trunk@30237 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-51/doc/results.dox
deal.II/examples/step-51/step-51.cc

index ae4bf1cc53e70d7484d62e2f6c95c65f95f009bb..1b7193f9dc1681a09ae3c4ac0edfacedcf409253 100644 (file)
@@ -5,135 +5,38 @@ The program writes convergence tables to the screen while running:
 
 
 @code
-examples/\step-51> make run
-============================ Running \step-51
-Solving with Q1 elements, adaptive refinement
-=============================================
-
-Cycle 0:
-   Number of active cells:       4
-   Number of degrees of freedom: 9
-Cycle 1:
-   Number of active cells:       13
-   Number of degrees of freedom: 22
-Cycle 2:
-   Number of active cells:       31
-   Number of degrees of freedom: 46
-Cycle 3:
-   Number of active cells:       64
-   Number of degrees of freedom: 87
-Cycle 4:
-   Number of active cells:       127
-   Number of degrees of freedom: 160
-Cycle 5:
-   Number of active cells:       244
-   Number of degrees of freedom: 297
-Cycle 6:
-   Number of active cells:       466
-   Number of degrees of freedom: 543
-
-cycle cells dofs    L2        H1      Linfty
-    0     4    9 1.198e+00 2.732e+00 1.383e+00
-    1    13   22 8.795e-02 1.193e+00 1.816e-01
-    2    31   46 8.147e-02 1.167e+00 1.654e-01
-    3    64   87 7.702e-02 1.077e+00 1.310e-01
-    4   127  160 4.643e-02 7.988e-01 6.745e-02
-    5   244  297 2.470e-02 5.568e-01 3.668e-02
-    6   466  543 1.622e-02 4.107e-01 2.966e-02
-
-Solving with Q1 elements, global refinement
-===========================================
-
-Cycle 0:
-   Number of active cells:       4
-   Number of degrees of freedom: 9
-Cycle 1:
-   Number of active cells:       16
-   Number of degrees of freedom: 25
-Cycle 2:
-   Number of active cells:       64
-   Number of degrees of freedom: 81
-Cycle 3:
-   Number of active cells:       256
-   Number of degrees of freedom: 289
-Cycle 4:
-   Number of active cells:       1024
-   Number of degrees of freedom: 1089
-Cycle 5:
-   Number of active cells:       4096
-   Number of degrees of freedom: 4225
-Cycle 6:
-   Number of active cells:       16384
-   Number of degrees of freedom: 16641
-
-cycle cells dofs     L2        H1      Linfty
-    0     4     9 1.198e+00 2.732e+00 1.383e+00
-    1    16    25 8.281e-02 1.190e+00 1.808e-01
-    2    64    81 8.142e-02 1.129e+00 1.294e-01
-    3   256   289 2.113e-02 5.828e-01 4.917e-02
-    4  1024  1089 5.319e-03 2.934e-01 1.359e-02
-    5  4096  4225 1.332e-03 1.469e-01 3.482e-03
-    6 16384 16641 3.332e-04 7.350e-02 8.758e-04
-
-n cells  H1      L2
-      0     4 2.732e+00    - 1.198e+00     -    -
-      1    16 1.190e+00 1.20 8.281e-02 14.47 3.86
-      2    64 1.129e+00 0.08 8.142e-02  1.02 0.02
-      3   256 5.828e-01 0.95 2.113e-02  3.85 1.95
-      4  1024 2.934e-01 0.99 5.319e-03  3.97 1.99
-      5  4096 1.469e-01 1.00 1.332e-03  3.99 2.00
-      6 16384 7.350e-02 1.00 3.332e-04  4.00 2.00
-
-Solving with Q2 elements, global refinement
-===========================================
-
-Cycle 0:
-   Number of active cells:       4
-   Number of degrees of freedom: 25
-Cycle 1:
-   Number of active cells:       16
-   Number of degrees of freedom: 81
-Cycle 2:
-   Number of active cells:       64
-   Number of degrees of freedom: 289
-Cycle 3:
-   Number of active cells:       256
-   Number of degrees of freedom: 1089
-Cycle 4:
-   Number of active cells:       1024
-   Number of degrees of freedom: 4225
-Cycle 5:
-   Number of active cells:       4096
-   Number of degrees of freedom: 16641
-Cycle 6:
-   Number of active cells:       16384
-   Number of degrees of freedom: 66049
-
-cycle cells dofs     L2        H1      Linfty
-    0     4    25 1.433e+00 2.445e+00 1.286e+00
-    1    16    81 7.912e-02 1.168e+00 1.728e-01
-    2    64   289 7.755e-03 2.511e-01 1.991e-02
-    3   256  1089 9.969e-04 6.235e-02 2.764e-03
-    4  1024  4225 1.265e-04 1.571e-02 3.527e-04
-    5  4096 16641 1.587e-05 3.937e-03 4.343e-05
-    6 16384 66049 1.986e-06 9.847e-04 5.402e-06
-
-n cells  H1      L2
-      0     4 2.445e+00    - 1.433e+00     -    -
-      1    16 1.168e+00 1.07 7.912e-02 18.11 4.18
-      2    64 2.511e-01 2.22 7.755e-03 10.20 3.35
-      3   256 6.235e-02 2.01 9.969e-04  7.78 2.96
-      4  1024 1.571e-02 1.99 1.265e-04  7.88 2.98
-      5  4096 3.937e-03 2.00 1.587e-05  7.97 2.99
-      6 16384 9.847e-04 2.00 1.986e-06  7.99 3.00
+Q1 elements:
+   16    80 4.570e+00    - 1.221e+01    - 4.333e+00    -
+   36   168 1.869e+00 2.20 7.299e+00 1.27 1.734e+00 2.26
+   64   288 7.177e-01 3.33 4.218e+00 1.91 2.538e-01 6.68
+  144   624 2.729e-01 2.38 1.867e+00 2.01 6.110e-02 3.51
+  256  1088 1.493e-01 2.10 1.046e+00 2.01 2.878e-02 2.62
+  576  2400 6.964e-02 1.88 4.847e-01 1.90 9.202e-03 2.81
+ 1024  4224 4.018e-02 1.91 2.785e-01 1.93 4.027e-03 2.87
+ 2304  9408 1.831e-02 1.94 1.264e-01 1.95 1.236e-03 2.91
+ 4096 16640 1.043e-02 1.96 7.185e-02 1.96 5.306e-04 2.94
+ 9216 37248 4.690e-03 1.97 3.228e-02 1.97 1.600e-04 2.96
+
+Q3 elements:
+   16   160 2.398e-01    - 1.873e+00    - 1.354e-01    -
+   36   336 5.843e-02 3.48 5.075e-01 3.22 1.882e-02 4.87
+   64   576 3.466e-02 1.82 2.534e-01 2.41 4.326e-03 5.11
+  144  1248 8.297e-03 3.53 5.925e-02 3.58 6.330e-04 4.74
+  256  2176 2.254e-03 4.53 1.636e-02 4.47 1.403e-04 5.24
+  576  4800 4.558e-04 3.94 3.278e-03 3.96 1.844e-05 5.01
+ 1024  8448 1.471e-04 3.93 1.052e-03 3.95 4.378e-06 5.00
+ 2304 18816 2.956e-05 3.96 2.104e-04 3.97 5.751e-07 5.01
+ 4096 33280 9.428e-06 3.97 6.697e-05 3.98 1.362e-07 5.01
+ 9216 74496 1.876e-06 3.98 1.330e-05 3.99 1.817e-08 4.97
 @endcode
 
 
-One can see the error reduction upon grid refinement, and for the
-cases where global refinement was performed, also the convergence
-rates can be seen. The linear and quadratic convergence rates of Q1
-and Q2 elements in the $H^1$ norm can clearly be seen, as
-are the quadratic and cubic rates in the $L_2$ norm.
-
-
+One can see the error reduction upon grid refinement, and for the cases where
+global refinement was performed, also the convergence rates can be seen. The
+quadratic convergence rates of Q1 elements in the $L_2$ norm for both the
+scalar variable and the gradient variable can clearly be seen, as is the cubic
+rate for the postprocessed scalar variable in the $L_2$ norm. Likewise, the
+scalar variable and gradient for Q3 elements converge at fourth order and the
+postprocessed scalar variable at fifth order.
 
+The same convergence rates are observed in 3d.
\ No newline at end of file
index 256e48dbba1dd29cd6bcd4d0697273741beaf760..41d3975d36f20a60dc93aa52ad1693297ffd673e 100644 (file)
 #include <deal.II/lac/chunk_sparse_matrix.h>
 #include <deal.II/numerics/data_out_faces.h>
 
-namespace Step51
-{
-  using namespace dealii;
+using namespace dealii;
 
-  // @sect3{Equation data}
+// @sect3{Equation data}
 
-  // The structure of the analytic solution is the same as in step-7. There
-  // are two exceptions. Firstly, we also create a solution for the 3d case,
-  // and secondly, we take into account the convection velocity in the right
-  // hand side that is variable in this case.
-  template <int dim>
-  class SolutionBase
-  {
-  protected:
-    static const unsigned int  n_source_centers = 3;
-    static const Point<dim>    source_centers[n_source_centers];
-    static const double        width;
-  };
+// The structure of the analytic solution is the same as in step-7. There
+// are two exceptions. Firstly, we also create a solution for the 3d case,
+// and secondly, we take into account the convection velocity in the right
+// hand side that is variable in this case.
+template <int dim>
+class SolutionBase
+{
+protected:
+  static const unsigned int  n_source_centers = 3;
+  static const Point<dim>    source_centers[n_source_centers];
+  static const double        width;
+};
 
 
-  template <>
-  const Point<1>
-  SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers]
-    = { Point<1>(-1.0 / 3.0),
-        Point<1>(0.0),
-        Point<1>(+1.0 / 3.0)
-      };
+template <>
+const Point<1>
+SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers]
+= { Point<1>(-1.0 / 3.0),
+    Point<1>(0.0),
+    Point<1>(+1.0 / 3.0)
+};
 
 
-  template <>
-  const Point<2>
-  SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers]
-    = { Point<2>(-0.5, +0.5),
-        Point<2>(-0.5, -0.5),
-        Point<2>(+0.5, -0.5)
-      };
+template <>
+const Point<2>
+SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers]
+= { Point<2>(-0.5, +0.5),
+    Point<2>(-0.5, -0.5),
+    Point<2>(+0.5, -0.5)
+};
 
-  template <>
-  const Point<3>
-  SolutionBase<3>::source_centers[SolutionBase<3>::n_source_centers]
-  = { Point<3>(-0.5, +0.5, 0.25),
-      Point<3>(-0.6, -0.5, -0.125),
-      Point<3>(+0.5, -0.5, 0.5)   };
+template <>
+const Point<3>
+SolutionBase<3>::source_centers[SolutionBase<3>::n_source_centers]
+= { Point<3>(-0.5, +0.5, 0.25),
+    Point<3>(-0.6, -0.5, -0.125),
+    Point<3>(+0.5, -0.5, 0.5)   };
 
-  template <int dim>
-  const double SolutionBase<dim>::width = 1./5.;
+template <int dim>
+const double SolutionBase<dim>::width = 1./5.;
 
 
 
-  template <int dim>
-  class ConvectionVelocity : public TensorFunction<1,dim>
-  {
-  public:
-    ConvectionVelocity() : TensorFunction<1,dim>() {}
+template <int dim>
+class ConvectionVelocity : public TensorFunction<1,dim>
+{
+public:
+  ConvectionVelocity() : TensorFunction<1,dim>() {}
 
-    virtual Tensor<1,dim> value (const Point<dim> &p) const;
-  };
+  virtual Tensor<1,dim> value (const Point<dim> &p) const;
+};
 
 
 
-  template <int dim>
-  Tensor<1,dim>
-  ConvectionVelocity<dim>::value(const Point<dim> &p) const
-  {
-    Tensor<1,dim> convection;
-    switch (dim)
-      {
-      case 1:
-        convection[0] = 1;
-        break;
-      case 2:
-        convection[0] = p[1];
-        convection[1] = -p[0];
-        break;
-      case 3:
-        convection[0] = p[1];
-        convection[1] = -p[0];
-        convection[2] = 1;
-        break;
-      default:
-        Assert(false, ExcNotImplemented());
-      }
-    return convection;
-  }
+template <int dim>
+Tensor<1,dim>
+ConvectionVelocity<dim>::value(const Point<dim> &p) const
+{
+  Tensor<1,dim> convection;
+  switch (dim)
+    {
+    case 1:
+      convection[0] = 1;
+      break;
+    case 2:
+      convection[0] = p[1];
+      convection[1] = -p[0];
+      break;
+    case 3:
+      convection[0] = p[1];
+      convection[1] = -p[0];
+      convection[2] = 1;
+      break;
+    default:
+      Assert(false, ExcNotImplemented());
+    }
+  return convection;
+}
 
 
-  template <int dim>
-  class Solution : public Function<dim>,
-                   protected SolutionBase<dim>
-  {
-  public:
-    Solution () : Function<dim>() {}
+template <int dim>
+class Solution : public Function<dim>,
+                 protected SolutionBase<dim>
+{
+public:
+  Solution () : Function<dim>() {}
 
-    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 Tensor<1,dim> gradient (const Point<dim>   &p,
-                                   const unsigned int  component = 0) const;
-  };
+  virtual Tensor<1,dim> gradient (const Point<dim>   &p,
+                                  const unsigned int  component = 0) const;
+};
 
 
 
-  template <int dim>
-  double Solution<dim>::value (const Point<dim>   &p,
-                               const unsigned int) const
-  {
-    double return_value = 0;
-    for (unsigned int i=0; i<this->n_source_centers; ++i)
-      {
-        const Point<dim> x_minus_xi = p - this->source_centers[i];
-        return_value += std::exp(-x_minus_xi.square() /
-                                 (this->width * this->width));
-      }
+template <int dim>
+double Solution<dim>::value (const Point<dim>   &p,
+                             const unsigned int) const
+{
+  double return_value = 0;
+  for (unsigned int i=0; i<this->n_source_centers; ++i)
+    {
+      const Point<dim> x_minus_xi = p - this->source_centers[i];
+      return_value += std::exp(-x_minus_xi.square() /
+                               (this->width * this->width));
+    }
 
-    return return_value /
-      Utilities::fixed_power<dim>(std::sqrt(2. * numbers::PI) * this->width);
-  }
+  return return_value /
+    Utilities::fixed_power<dim>(std::sqrt(2. * numbers::PI) * this->width);
+}
 
 
 
-  template <int dim>
-  Tensor<1,dim> Solution<dim>::gradient (const Point<dim>   &p,
-                                         const unsigned int) const
-  {
-    Tensor<1,dim> return_value;
+template <int dim>
+Tensor<1,dim> Solution<dim>::gradient (const Point<dim>   &p,
+                                       const unsigned int) const
+{
+  Tensor<1,dim> return_value;
 
-    for (unsigned int i=0; i<this->n_source_centers; ++i)
-      {
-        const Point<dim> x_minus_xi = p - this->source_centers[i];
+  for (unsigned int i=0; i<this->n_source_centers; ++i)
+    {
+      const Point<dim> x_minus_xi = p - this->source_centers[i];
 
-        return_value += (-2 / (this->width * this->width) *
-                         std::exp(-x_minus_xi.square() /
-                                  (this->width * this->width)) *
-                         x_minus_xi);
-      }
+      return_value += (-2 / (this->width * this->width) *
+                       std::exp(-x_minus_xi.square() /
+                                (this->width * this->width)) *
+                       x_minus_xi);
+    }
 
-    return return_value / Utilities::fixed_power<dim>(std::sqrt(2 * numbers::PI) *
-                                                      this->width);
-  }
+  return return_value / Utilities::fixed_power<dim>(std::sqrt(2 * numbers::PI) *
+                                                    this->width);
+}
 
 
 
-  template <int dim>
-  class SolutionAndGradient : public Function<dim>,
-                              protected SolutionBase<dim>
-  {
-  public:
-    SolutionAndGradient () : Function<dim>(dim) {}
+template <int dim>
+class SolutionAndGradient : public Function<dim>,
+                            protected SolutionBase<dim>
+{
+public:
+  SolutionAndGradient () : Function<dim>(dim) {}
 
-    virtual void vector_value (const Point<dim>   &p,
-                               Vector<double>     &v) const
-    {
-      AssertDimension(v.size(), dim+1);
-      Solution<dim> solution;
-      Tensor<1,dim> grad = solution.gradient(p);
-      for (unsigned int d=0; d<dim; ++d)
-        v[d] = -grad[d];
-      v[dim] = solution.value(p);
-    }
-  };
+  virtual void vector_value (const Point<dim>   &p,
+                             Vector<double>     &v) const
+  {
+    AssertDimension(v.size(), dim+1);
+    Solution<dim> solution;
+    Tensor<1,dim> grad = solution.gradient(p);
+    for (unsigned int d=0; d<dim; ++d)
+      v[d] = -grad[d];
+    v[dim] = solution.value(p);
+  }
+};
 
 
 
-  template <int dim>
-  class RightHandSide : public Function<dim>,
-                        protected SolutionBase<dim>
-  {
-  public:
-    RightHandSide () : Function<dim>() {}
+template <int dim>
+class RightHandSide : public Function<dim>,
+                      protected SolutionBase<dim>
+{
+public:
+  RightHandSide () : Function<dim>() {}
 
-    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;
 
-  private:
-    const ConvectionVelocity<dim> convection_velocity;
-  };
+private:
+  const ConvectionVelocity<dim> convection_velocity;
+};
 
 
-  template <int dim>
-  double RightHandSide<dim>::value (const Point<dim>   &p,
-                                    const unsigned int) const
-  {
-    Tensor<1,dim> convection = convection_velocity.value(p);
-    double return_value = 0;
-    for (unsigned int i=0; i<this->n_source_centers; ++i)
-      {
-        const Point<dim> x_minus_xi = p - this->source_centers[i];
-
-        return_value +=
-          ((2*dim - 2*convection*x_minus_xi - 4*x_minus_xi.square()/
-            (this->width * this->width)) /
-           (this->width * this->width) *
-           std::exp(-x_minus_xi.square() /
-                    (this->width * this->width)));
-      }
+template <int dim>
+double RightHandSide<dim>::value (const Point<dim>   &p,
+                                  const unsigned int) const
+{
+  Tensor<1,dim> convection = convection_velocity.value(p);
+  double return_value = 0;
+  for (unsigned int i=0; i<this->n_source_centers; ++i)
+    {
+      const Point<dim> x_minus_xi = p - this->source_centers[i];
+
+      return_value +=
+        ((2*dim - 2*convection*x_minus_xi - 4*x_minus_xi.square()/
+          (this->width * this->width)) /
+         (this->width * this->width) *
+         std::exp(-x_minus_xi.square() /
+                  (this->width * this->width)));
+    }
 
-    return return_value / Utilities::fixed_power<dim>(std::sqrt(2 * numbers::PI)
-                                                      * this->width);
-  }
+  return return_value / Utilities::fixed_power<dim>(std::sqrt(2 * numbers::PI)
+                                                    * this->width);
+}
 
 
 
-  template <int dim>
-  class Step51
-  {
-  public:
-    enum RefinementMode
+template <int dim>
+class Step51
+{
+public:
+  enum RefinementMode
     {
       global_refinement, adaptive_refinement
     };
 
-    Step51 (const unsigned int degree,
-            const RefinementMode refinement_mode);
-    void run ();
+  Step51 (const unsigned int degree,
+          const RefinementMode refinement_mode);
+  void run ();
+
+private:
+  void setup_system ();
+  void assemble_system (const bool reconstruct_trace = false);
+  void solve ();
+  void postprocess ();
+  void refine_mesh ();
+  void output_results (const unsigned int cycle);
 
-  private:
-    void setup_system ();
-    void assemble_system (const bool reconstruct_trace = false);
-    void solve ();
-    void postprocess ();
-    void refine_mesh ();
-    void output_results (const unsigned int cycle);
+  Triangulation<dim>   triangulation;
 
-    Triangulation<dim>   triangulation;
+  const MappingQ<dim>  mapping;
 
-    const MappingQ<dim>  mapping;
+  FESystem<dim>        fe_local;
+  DoFHandler<dim>      dof_handler_local;
 
-    FESystem<dim>        fe_local;
-    DoFHandler<dim>      dof_handler_local;
+  FE_FaceQ<dim>        fe;
+  DoFHandler<dim>      dof_handler;
 
-    FE_FaceQ<dim>        fe;
-    DoFHandler<dim>      dof_handler;
+  FE_DGQ<dim>          fe_u_post;
+  DoFHandler<dim>      dof_handler_u_post;
 
-    FE_DGQ<dim>          fe_u_post;
-    DoFHandler<dim>      dof_handler_u_post;
+  ConstraintMatrix     constraints;
+  ChunkSparsityPattern sparsity_pattern;
+  ChunkSparseMatrix<double> system_matrix;
 
-    ConstraintMatrix     constraints;
-    ChunkSparsityPattern sparsity_pattern;
-    ChunkSparseMatrix<double> system_matrix;
+  Vector<double>       solution;
+  Vector<double>       system_rhs;
 
-    Vector<double>       solution;
-    Vector<double>       system_rhs;
+  Vector<double>       solution_local;
+  Vector<double>       solution_u_post;
 
-    Vector<double>       solution_local;
-    Vector<double>       solution_u_post;
+  const RefinementMode refinement_mode;
 
-    const RefinementMode refinement_mode;
+  ConvergenceTable     convergence_table;
+};
 
-    ConvergenceTable     convergence_table;
-  };
 
 
+template <int dim>
+Step51<dim>::Step51 (const unsigned int degree,
+                     const RefinementMode refinement_mode) :
+  mapping  (3),
+  fe_local (FE_DGQ<dim>(degree), dim,
+            FE_DGQ<dim>(degree), 1),
+  dof_handler_local (triangulation),
+  fe (degree),
+  dof_handler (triangulation),
+  fe_u_post (degree+1),
+  dof_handler_u_post (triangulation),
+  refinement_mode (refinement_mode)
+{}
 
-  template <int dim>
-  Step51<dim>::Step51 (const unsigned int degree,
-                       const RefinementMode refinement_mode) :
-    mapping  (3),
-    fe_local (FE_DGQ<dim>(degree), dim,
-              FE_DGQ<dim>(degree), 1),
-    dof_handler_local (triangulation),
-    fe (degree),
-    dof_handler (triangulation),
-    fe_u_post (degree+1),
-    dof_handler_u_post (triangulation),
-    refinement_mode (refinement_mode)
-  {}
 
 
+template <int dim>
+void
+Step51<dim>::setup_system ()
+{
+  dof_handler_local.distribute_dofs(fe_local);
+  dof_handler.distribute_dofs(fe);
+  dof_handler_u_post.distribute_dofs(fe_u_post);
+
+  std::cout << "   Number of degrees of freedom: "
+            << dof_handler.n_dofs()
+            << std::endl;
+
+  solution.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+
+  solution_local.reinit (dof_handler_local.n_dofs());
+  solution_u_post.reinit (dof_handler_u_post.n_dofs());
+
+  constraints.clear ();
+  DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+  std::map<unsigned int,double> boundary_values;
+  typename FunctionMap<dim>::type boundary_functions;
+  Solution<dim> solution;
+  boundary_functions[0] = &solution;
+  VectorTools::project_boundary_values (mapping, dof_handler,
+                                        boundary_functions,
+                                        QGauss<dim-1>(fe.degree+1),
+                                        boundary_values);
+  for (std::map<unsigned int,double>::iterator it = boundary_values.begin();
+       it != boundary_values.end(); ++it)
+    if (constraints.is_constrained(it->first) == false)
+      {
+        constraints.add_line(it->first);
+        constraints.set_inhomogeneity(it->first, it->second);
+      }
+  constraints.close ();
 
-  template <int dim>
-  void
-  Step51<dim>::setup_system ()
   {
-    dof_handler_local.distribute_dofs(fe_local);
-    dof_handler.distribute_dofs(fe);
-    dof_handler_u_post.distribute_dofs(fe_u_post);
+    CompressedSimpleSparsityPattern csp (dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern (dof_handler, csp,
+                                     constraints, false);
+    sparsity_pattern.copy_from(csp, fe.dofs_per_face);
+  }
+  system_matrix.reinit (sparsity_pattern);
+}
 
-    std::cout << "   Number of degrees of freedom: "
-              << dof_handler.n_dofs()
-              << std::endl;
 
-    solution.reinit (dof_handler.n_dofs());
-    system_rhs.reinit (dof_handler.n_dofs());
 
-    solution_local.reinit (dof_handler_local.n_dofs());
-    solution_u_post.reinit (dof_handler_u_post.n_dofs());
+template <int dim>
+void
+Step51<dim>::assemble_system (const bool trace_reconstruct)
+{
+  QGauss<dim>   quadrature_formula(fe.degree+1);
+  QGauss<dim-1> face_quadrature_formula(fe.degree+1);
+
+  FEValues<dim> fe_values_local (mapping, fe_local, quadrature_formula,
+                                 update_values | update_gradients |
+                                 update_JxW_values | update_quadrature_points);
+  FEFaceValues<dim> fe_face_values (mapping, fe, face_quadrature_formula,
+                                    update_values | update_normal_vectors |
+                                    update_quadrature_points |
+                                    update_JxW_values);
+  FEFaceValues<dim> fe_face_values_local (mapping, fe_local,
+                                          face_quadrature_formula,
+                                          update_values);
+
+  const unsigned int n_q_points    = quadrature_formula.size();
+  const unsigned int n_face_q_points = face_quadrature_formula.size();
+
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int loc_dofs_per_cell = fe_local.dofs_per_cell;
+
+  FullMatrix<double> ll_matrix (loc_dofs_per_cell, loc_dofs_per_cell);
+  FullMatrix<double> lf_matrix (loc_dofs_per_cell, dofs_per_cell);
+  FullMatrix<double> fl_matrix (dofs_per_cell, loc_dofs_per_cell);
+  FullMatrix<double> tmp_matrix (dofs_per_cell, loc_dofs_per_cell);
+  FullMatrix<double> ff_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>     l_rhs (loc_dofs_per_cell);
+  Vector<double>     f_rhs (dofs_per_cell);
+  Vector<double>     tmp_rhs (loc_dofs_per_cell);
+
+  std::vector<types::global_dof_index> dof_indices (dofs_per_cell);
+  std::vector<types::global_dof_index> loc_dof_indices (loc_dofs_per_cell);
+
+  std::vector<Tensor<1,dim> > q_phi (loc_dofs_per_cell);
+  std::vector<double>         q_phi_div (loc_dofs_per_cell);
+  std::vector<double>         u_phi (loc_dofs_per_cell);
+  std::vector<Tensor<1,dim> > u_phi_grad (loc_dofs_per_cell);
+  std::vector<double>         tr_phi (dofs_per_cell);
+
+  std::vector<double> trace_values(n_face_q_points);
+
+  // Choose stabilization parameter to be 5 * diffusion = 5
+  const double tau_stab_diffusion = 5.;
+
+  ConvectionVelocity<dim> convection_velocity;
+  RightHandSide<dim> right_hand_side;
+  const Solution<dim> exact_solution;
+
+  const FEValuesExtractors::Vector fluxes (0);
+  const FEValuesExtractors::Scalar scalar (dim);
+
+  std::vector<std::vector<unsigned int> >
+    fe_local_support_on_face(GeometryInfo<dim>::faces_per_cell);
+  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+    for (unsigned int i=0; i<loc_dofs_per_cell; ++i)
+      if (fe_local.has_support_on_face(i,face))
+        fe_local_support_on_face[face].push_back(i);
+  std::vector<std::vector<unsigned int> >
+    fe_support_on_face(GeometryInfo<dim>::faces_per_cell);
+  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+    for (unsigned int i=0; i<dofs_per_cell; ++i)
+      if (fe.has_support_on_face(i,face))
+        fe_support_on_face[face].push_back(i);
+
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    loc_cell = dof_handler_local.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell, ++loc_cell)
+    {
+      ll_matrix = 0;
+      l_rhs = 0;
+      if (!trace_reconstruct)
+        {
+          lf_matrix = 0;
+          fl_matrix = 0;
+          ff_matrix = 0;
+          f_rhs = 0;
+        }
+      fe_values_local.reinit (loc_cell);
 
-    constraints.clear ();
-    DoFTools::make_hanging_node_constraints (dof_handler, constraints);
-    std::map<unsigned int,double> boundary_values;
-    typename FunctionMap<dim>::type boundary_functions;
-    Solution<dim> solution;
-    boundary_functions[0] = &solution;
-    VectorTools::project_boundary_values (mapping, dof_handler,
-                                          boundary_functions,
-                                          QGauss<dim-1>(fe.degree+1),
-                                          boundary_values);
-    for (std::map<unsigned int,double>::iterator it = boundary_values.begin();
-         it != boundary_values.end(); ++it)
-      if (constraints.is_constrained(it->first) == false)
+      for (unsigned int q=0; q<n_q_points; ++q)
         {
-          constraints.add_line(it->first);
-          constraints.set_inhomogeneity(it->first, it->second);
+          const double rhs_value
+            = right_hand_side.value(fe_values_local.quadrature_point(q));
+          const Tensor<1,dim> convection
+            = convection_velocity.value(fe_values_local.quadrature_point(q));
+          const double JxW = fe_values_local.JxW(q);
+          for (unsigned int k=0; k<loc_dofs_per_cell; ++k)
+            {
+              q_phi[k] = fe_values_local[fluxes].value(k,q);
+              q_phi_div[k] = fe_values_local[fluxes].divergence(k,q);
+              u_phi[k] = fe_values_local[scalar].value(k,q);
+              u_phi_grad[k] = fe_values_local[scalar].gradient(k,q);
+            }
+          for (unsigned int i=0; i<loc_dofs_per_cell; ++i)
+            {
+              for (unsigned int j=0; j<loc_dofs_per_cell; ++j)
+                ll_matrix(i,j) += (
+                                   q_phi[i] * q_phi[j]
+                                   -
+                                   q_phi_div[i] * u_phi[j]
+                                   +
+                                   u_phi[i] * q_phi_div[j]
+                                   -
+                                   (u_phi_grad[i] * convection) * u_phi[j]
+                                   ) * JxW;
+              l_rhs(i) += u_phi[i] * rhs_value * JxW;
+            }
         }
-    constraints.close ();
 
-    {
-      CompressedSimpleSparsityPattern csp (dof_handler.n_dofs());
-      DoFTools::make_sparsity_pattern (dof_handler, csp,
-                                       constraints, false);
-      sparsity_pattern.copy_from(csp, fe.dofs_per_face);
-    }
-    system_matrix.reinit (sparsity_pattern);
-  }
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+        {
+          fe_face_values_local.reinit(loc_cell, face);
+          fe_face_values.reinit(cell, face);
+          if (trace_reconstruct)
+            fe_face_values.get_function_values (solution, trace_values);
 
+          for (unsigned int q=0; q<n_face_q_points; ++q)
+            {
+              const double JxW = fe_face_values.JxW(q);
+              const Point<dim> normal = fe_face_values.normal_vector(q);
+              const Tensor<1,dim> convection
+                = convection_velocity.value(fe_face_values.quadrature_point(q));
+              const double tau_stab = (tau_stab_diffusion +
+                                       std::abs(convection * normal));
+
+              for (unsigned int k=0; k<fe_local_support_on_face[face].size(); ++k)
+                {
+                  const unsigned int kk=fe_local_support_on_face[face][k];
+                  q_phi[k] = fe_face_values_local[fluxes].value(kk,q);
+                  u_phi[k] = fe_face_values_local[scalar].value(kk,q);
+                }
+              if (!trace_reconstruct)
+                {
+                  for (unsigned int k=0; k<fe_support_on_face[face].size(); ++k)
+                    tr_phi[k] =
+                      fe_face_values.shape_value(fe_support_on_face[face][k],q);
+                  for (unsigned int i=0; i<fe_local_support_on_face[face].size(); ++i)
+                    for (unsigned int j=0; j<fe_support_on_face[face].size(); ++j)
+                      {
+                        const unsigned int ii=fe_local_support_on_face[face][i];
+                        const unsigned int jj=fe_support_on_face[face][j];
+                        lf_matrix(ii,jj) += (
+                                             (q_phi[i] * normal
+                                              +
+                                              (convection * normal -
+                                               tau_stab) * u_phi[i])
+                                             * tr_phi[j]
+                                             ) * JxW;
+                        fl_matrix(jj,ii) -= (
+                                             (q_phi[i] * normal
+                                              +
+                                              tau_stab * u_phi[i])
+                                             * tr_phi[j]
+                                             ) * JxW;
+                      }
 
+                  for (unsigned int i=0; i<fe_support_on_face[face].size(); ++i)
+                    for (unsigned int j=0; j<fe_support_on_face[face].size(); ++j)
+                      {
+                        const unsigned int ii=fe_support_on_face[face][i];
+                        const unsigned int jj=fe_support_on_face[face][j];
+                        ff_matrix(ii,jj) += (
+                                             (convection * normal - tau_stab) *
+                                             tr_phi[i] * tr_phi[j]
+                                             ) * JxW;
+                      }
 
-  template <int dim>
-  void
-  Step51<dim>::assemble_system (const bool trace_reconstruct)
-  {
-    QGauss<dim>   quadrature_formula(fe.degree+1);
-    QGauss<dim-1> face_quadrature_formula(fe.degree+1);
-
-    FEValues<dim> fe_values_local (mapping, fe_local, quadrature_formula,
-                                   update_values | update_gradients |
-                                   update_JxW_values | update_quadrature_points);
-    FEFaceValues<dim> fe_face_values (mapping, fe, face_quadrature_formula,
-                                      update_values | update_normal_vectors |
-                                      update_quadrature_points |
-                                      update_JxW_values);
-    FEFaceValues<dim> fe_face_values_local (mapping, fe_local,
-                                            face_quadrature_formula,
-                                            update_values);
-
-    const unsigned int n_q_points    = quadrature_formula.size();
-    const unsigned int n_face_q_points = face_quadrature_formula.size();
-
-    const unsigned int dofs_per_cell = fe.dofs_per_cell;
-    const unsigned int loc_dofs_per_cell = fe_local.dofs_per_cell;
-
-    FullMatrix<double> ll_matrix (loc_dofs_per_cell, loc_dofs_per_cell);
-    FullMatrix<double> lf_matrix (loc_dofs_per_cell, dofs_per_cell);
-    FullMatrix<double> fl_matrix (dofs_per_cell, loc_dofs_per_cell);
-    FullMatrix<double> tmp_matrix (dofs_per_cell, loc_dofs_per_cell);
-    FullMatrix<double> ff_matrix (dofs_per_cell, dofs_per_cell);
-    Vector<double>     l_rhs (loc_dofs_per_cell);
-    Vector<double>     f_rhs (dofs_per_cell);
-    Vector<double>     tmp_rhs (loc_dofs_per_cell);
-
-    std::vector<types::global_dof_index> dof_indices (dofs_per_cell);
-    std::vector<types::global_dof_index> loc_dof_indices (loc_dofs_per_cell);
-
-    ConvectionVelocity<dim> convection;
-    std::vector<Tensor<1,dim> > convection_values (n_q_points);
-    std::vector<Tensor<1,dim> > convection_values_face (n_face_q_points);
-
-    std::vector<double> trace_values(n_face_q_points);
-
-    // Choose stabilization parameter to be 5 * diffusion = 5
-    const double tau_stab_diffusion = 5.;
-    std::vector<double> tau_stab (n_q_points);
-
-    RightHandSide<dim> right_hand_side;
-    std::vector<double> rhs_values (n_q_points);
-
-    const Solution<dim> exact_solution;
-    std::vector<double> neumann_values (n_face_q_points);
-
-    const FEValuesExtractors::Vector gradients (0);
-    const FEValuesExtractors::Scalar values (dim);
-
-    typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    loc_cell = dof_handler_local.begin_active(),
-    endc = dof_handler.end();
-    for (; cell!=endc; ++cell, ++loc_cell)
-      {
-        if (!trace_reconstruct)
-          {
-            lf_matrix = 0;
-            fl_matrix = 0;
-            ff_matrix = 0;
-            f_rhs = 0;
-          }
-        fe_values_local.reinit (loc_cell);
-        right_hand_side.value_list (fe_values_local.get_quadrature_points(),
-                                    rhs_values);
-        convection.value_list(fe_values_local.get_quadrature_points(),
-                              convection_values);
-
-        for (unsigned int i=0; i<loc_dofs_per_cell; ++i)
-          for (unsigned int j=0; j<loc_dofs_per_cell; ++j)
-            {
-              double sum = 0;
-              for (unsigned int q=0; q<n_q_points; ++q)
-                sum += (fe_values_local[gradients].value(i,q) *
-                        fe_values_local[gradients].value(j,q)
-                        -
-                        fe_values_local[gradients].divergence(i,q) *
-                        fe_values_local[values].value(j,q)
-                        +
-                        fe_values_local[values].value(i,q) *
-                        fe_values_local[gradients].divergence(j,q)
-                        -
-                        fe_values_local[values].value(j,q) *
-                        (fe_values_local[values].gradient(i,q) *
-                         convection_values[q])
-                        ) * fe_values_local.JxW(q);
-              ll_matrix(i,j) = sum;
-            }
-        for (unsigned int i=0; i<loc_dofs_per_cell; ++i)
-          {
-            double sum = 0.;
-            for (unsigned int q=0; q<n_q_points; ++q)
-              sum += rhs_values[q] * fe_values_local.JxW(q) *
-                fe_values_local[values].value(i,q);
-            l_rhs(i) = sum;
-          }
-
-        for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-          {
-            fe_face_values_local.reinit(loc_cell, face);
-            fe_face_values.reinit(cell, face);
-            const std::vector<double> &JxW = fe_face_values.get_JxW_values();
-            const std::vector<Point<dim> > &normals =
-              fe_face_values.get_normal_vectors();
-            convection.value_list(fe_face_values.get_quadrature_points(),
-                                  convection_values_face);
-            for (unsigned int q=0; q<n_face_q_points; ++q)
-              tau_stab[q] = (tau_stab_diffusion +
-                             std::abs(convection_values_face[q] * normals[q]));
-            if (!trace_reconstruct)
-              {
-                for (unsigned int i=0; i<dofs_per_cell; ++i)
-                  for (unsigned int j=0; j<loc_dofs_per_cell; ++j)
+                  if (cell->face(face)->at_boundary()
+                      &&
+                      (cell->face(face)->boundary_indicator() == 1))
                     {
-                      double sum_lf = 0., sum_fl = 0.;
-                      for (unsigned int q=0; q<n_face_q_points; ++q)
+                      const double neumann_value =
+                        exact_solution.value(fe_face_values.quadrature_point(q));
+                      for (unsigned int i=0; i<fe_support_on_face[face].size(); ++i)
                         {
-                          sum_lf += (fe_face_values.shape_value(i,q) *
-                                     (fe_face_values_local[gradients].value(j,q) *
-                                      normals[q]
-                                      +
-                                      (convection_values_face[q] *
-                                       normals[q]
-                                       -
-                                       tau_stab[q]) *
-                                      fe_face_values_local[values].value(j,q))
-                                     ) * JxW[q];
-                          sum_fl += (fe_face_values.shape_value(i,q) *
-                                     (fe_face_values_local[gradients].value(j,q) *
-                                      normals[q]
-                                      +
-                                      tau_stab[q] *
-                                      fe_face_values_local[values].value(j,q))
-                                     ) * JxW[q];
+                          const unsigned int ii=fe_support_on_face[face][i];
+                          f_rhs(ii) -= tr_phi[i] * neumann_value * JxW;
                         }
-                      lf_matrix(j,i) += sum_lf;
-                      fl_matrix(i,j) -= sum_fl;
                     }
-                for (unsigned int i=0; i<dofs_per_cell; ++i)
-                  for (unsigned int j=0; j<dofs_per_cell; ++j)
-                    {
-                      double sum = 0;
-                      for (unsigned int q=0; q<n_face_q_points; ++q)
-                        sum += ((convection_values_face[q] * normals[q]
-                                 -
-                                 tau_stab[q]
-                                 ) *
-                                fe_face_values.shape_value(i,q) *
-                                fe_face_values.shape_value(j,q)
-                                ) * JxW[q];
-                      ff_matrix(i,j) += sum;
-                    }
-                if (cell->face(face)->at_boundary()
-                    &&
-                    (cell->face(face)->boundary_indicator() == 1))
+                }
+
+              for (unsigned int i=0; i<fe_local_support_on_face[face].size(); ++i)
+                for (unsigned int j=0; j<fe_local_support_on_face[face].size(); ++j)
                   {
-                    exact_solution.value_list(fe_face_values.get_quadrature_points(),
-                                              neumann_values);
-                    for (unsigned int i=0; i<dofs_per_cell; ++i)
-                      {
-                        double sum = 0;
-                        for (unsigned int q=0; q<n_face_q_points; ++q)
-                          sum -= (fe_face_values.shape_value(i,q) *
-                                  neumann_values[q]) * JxW[q];
-                        f_rhs(i) += sum;
-                      }
+                    const unsigned int ii=fe_local_support_on_face[face][i];
+                    const unsigned int jj=fe_local_support_on_face[face][j];
+                    ll_matrix(ii,jj) += tau_stab * u_phi[i] * u_phi[j] * JxW;
                   }
-              }
-            for (unsigned int i=0; i<loc_dofs_per_cell; ++i)
-              for (unsigned int j=0; j<loc_dofs_per_cell; ++j)
-                {
-                  double sum = 0;
-                  for (unsigned int q=0; q<n_face_q_points; ++q)
-                    sum += (tau_stab[q] *
-                            fe_face_values_local[values].value(i,q) *
-                            fe_face_values_local[values].value(j,q)) * JxW[q];
-                  ll_matrix(i,j) += sum;
-                }
 
-            // compute the local right hand side contributions from trace
-            if (trace_reconstruct)
-              {
-                fe_face_values.get_function_values (solution, trace_values);
-                for (unsigned int i=0; i<loc_dofs_per_cell; ++i)
+              // compute the local right hand side contributions from trace
+              if (trace_reconstruct)
+                for (unsigned int i=0; i<fe_local_support_on_face[face].size(); ++i)
                   {
-                    double sum = 0;
-                    for (unsigned int q=0; q<n_face_q_points; ++q)
-                      sum += ((fe_face_values_local[gradients].value(i,q) *
-                               normals[q]) *
-                              trace_values[q]
-                              +
-                              fe_face_values_local[values].value(i,q) *
-                              (convection_values_face[q] * normals[q]
-                               -
-                               tau_stab[q]) * trace_values[q]) * JxW[q];
-                    l_rhs(i) -= sum;
+                    const unsigned int ii=fe_local_support_on_face[face][i];
+                    l_rhs(ii) -= (q_phi[i] * normal
+                                  +
+                                  u_phi[i] * (convection * normal - tau_stab)
+                                  ) * trace_values[q] * JxW;
                   }
-              }
-          }
-
-        ll_matrix.gauss_jordan();
-        if (!trace_reconstruct)
-          {
-            fl_matrix.mmult(tmp_matrix, ll_matrix);
-            tmp_matrix.vmult_add(f_rhs, l_rhs);
-            tmp_matrix.mmult(ff_matrix, lf_matrix, true);
-            cell->get_dof_indices(dof_indices);
-            constraints.distribute_local_to_global (ff_matrix, f_rhs,
-                                                    dof_indices,
-                                                    system_matrix, system_rhs);
-          }
-        else
-          {
-            ll_matrix.vmult(tmp_rhs, l_rhs);
-            loc_cell->set_dof_values(tmp_rhs, solution_local);
-          }
-      }
-  }
+            }
+        }
 
+      ll_matrix.gauss_jordan();
+      if (trace_reconstruct == false)
+        {
+          fl_matrix.mmult(tmp_matrix, ll_matrix);
+          tmp_matrix.vmult_add(f_rhs, l_rhs);
+          tmp_matrix.mmult(ff_matrix, lf_matrix, true);
+          cell->get_dof_indices(dof_indices);
+          constraints.distribute_local_to_global (ff_matrix, f_rhs,
+                                                  dof_indices,
+                                                  system_matrix, system_rhs);
+        }
+      else
+        {
+          ll_matrix.vmult(tmp_rhs, l_rhs);
+          loc_cell->set_dof_values(tmp_rhs, solution_local);
+        }
+    }
+}
 
 
-  template <int dim>
-  void Step51<dim>::solve ()
-  {
-    SolverControl solver_control (system_matrix.m()*10,
-                                  1e-10*system_rhs.l2_norm());
-    SolverGMRES<> solver (solver_control, 50);
-    solver.solve (system_matrix, solution, system_rhs,
-                  PreconditionIdentity());
 
-    std::cout << "   Number of GMRES iterations: " << solver_control.last_step()
-              << std::endl;
+template <int dim>
+void Step51<dim>::solve ()
+{
+  SolverControl solver_control (system_matrix.m()*10,
+                                1e-10*system_rhs.l2_norm());
+  SolverGMRES<> solver (solver_control, 50);
+  solver.solve (system_matrix, solution, system_rhs,
+                PreconditionIdentity());
 
-    system_matrix.clear();
-    sparsity_pattern.reinit(0,0,0,1);
-    constraints.distribute(solution);
+  std::cout << "   Number of GMRES iterations: " << solver_control.last_step()
+            << std::endl;
 
-    // update local values
-    assemble_system(true);
-  }
+  system_matrix.clear();
+  sparsity_pattern.reinit(0,0,0,1);
+  constraints.distribute(solution);
 
+  // update local values
+  assemble_system(true);
+}
 
 
-  template <int dim>
-  void
-  Step51<dim>::postprocess()
-  {
-    const unsigned int n_active_cells=triangulation.n_active_cells();
-    Vector<float> difference_per_cell (triangulation.n_active_cells());
-
-    ComponentSelectFunction<dim> value_select (dim, dim+1);
-    VectorTools::integrate_difference (mapping, dof_handler_local,
-                                       solution_local,
-                                       SolutionAndGradient<dim>(),
-                                       difference_per_cell,
-                                       QGauss<dim>(fe.degree+2),
-                                       VectorTools::L2_norm,
-                                       &value_select);
-    const double L2_error = difference_per_cell.l2_norm();
-
-    ComponentSelectFunction<dim> gradient_select (std::pair<unsigned int,unsigned int>(0, dim),
-                                                  dim+1);
-    VectorTools::integrate_difference (mapping, dof_handler_local,
-                                       solution_local,
-                                       SolutionAndGradient<dim>(),
-                                       difference_per_cell,
-                                       QGauss<dim>(fe.degree+2),
-                                       VectorTools::L2_norm,
-                                       &gradient_select);
-    const double grad_error = difference_per_cell.l2_norm();
-
-    convergence_table.add_value("cells", n_active_cells);
-    convergence_table.add_value("dofs", dof_handler.n_dofs());
-    convergence_table.add_value("val L2", L2_error);
-    convergence_table.add_value("grad L2", grad_error);
-
-    // construct post-processed solution with (hopefully) higher order of
-    // accuracy
-    QGauss<dim> quadrature(fe_u_post.degree+1);
-    FEValues<dim> fe_values(mapping, fe_u_post, quadrature,
-                            update_values | update_JxW_values |
-                            update_gradients);
-
-    const unsigned int n_q_points = quadrature.size();
-    std::vector<double> u_values(n_q_points);
-    std::vector<Tensor<1,dim> > u_gradients(n_q_points);
-    FEValuesExtractors::Vector gradients(0);
-    FEValuesExtractors::Scalar values(dim);
-    FEValues<dim> fe_values_local(mapping, fe_local, quadrature, update_values);
-    FullMatrix<double> cell_matrix(fe_u_post.dofs_per_cell,
-                                   fe_u_post.dofs_per_cell);
-    Vector<double> cell_rhs(fe_u_post.dofs_per_cell);
-    Vector<double> cell_sol(fe_u_post.dofs_per_cell);
-
-    typename DoFHandler<dim>::active_cell_iterator
-      cell_loc = dof_handler_local.begin_active(),
-      cell = dof_handler_u_post.begin_active(),
-      endc = dof_handler_u_post.end();
-    for ( ; cell != endc; ++cell, ++cell_loc)
-      {
-        fe_values.reinit(cell);
-        fe_values_local.reinit(cell_loc);
-
-        fe_values_local[values].get_function_values(solution_local, u_values);
-        fe_values_local[gradients].get_function_values(solution_local, u_gradients);
-        for (unsigned int i=1; i<fe_u_post.dofs_per_cell; ++i)
-          {
-            for (unsigned int j=0; j<fe_u_post.dofs_per_cell; ++j)
-              {
-                double sum = 0;
-                for (unsigned int q=0; q<quadrature.size(); ++q)
-                  sum += (fe_values.shape_grad(i,q) *
-                          fe_values.shape_grad(j,q)
-                          ) * fe_values.JxW(q);
-                cell_matrix(i,j) = sum;
-              }
-            double sum = 0;
-            for (unsigned int q=0; q<quadrature.size(); ++q)
-              sum -= (fe_values.shape_grad(i,q) * u_gradients[q]
-                      ) * fe_values.JxW(q);
-            cell_rhs(i) = sum;
-          }
-        for (unsigned int j=0; j<fe_u_post.dofs_per_cell; ++j)
-          {
-            double sum = 0;
-            for (unsigned int q=0; q<quadrature.size(); ++q)
-              sum += fe_values.shape_value(j,q) * fe_values.JxW(q);
-            cell_matrix(0,j) = sum;
-          }
+
+template <int dim>
+void
+Step51<dim>::postprocess()
+{
+  const unsigned int n_active_cells=triangulation.n_active_cells();
+  Vector<float> difference_per_cell (triangulation.n_active_cells());
+
+  ComponentSelectFunction<dim> value_select (dim, dim+1);
+  VectorTools::integrate_difference (mapping, dof_handler_local,
+                                     solution_local,
+                                     SolutionAndGradient<dim>(),
+                                     difference_per_cell,
+                                     QGauss<dim>(fe.degree+2),
+                                     VectorTools::L2_norm,
+                                     &value_select);
+  const double L2_error = difference_per_cell.l2_norm();
+
+  ComponentSelectFunction<dim> gradient_select (std::pair<unsigned int,unsigned int>(0, dim),
+                                                dim+1);
+  VectorTools::integrate_difference (mapping, dof_handler_local,
+                                     solution_local,
+                                     SolutionAndGradient<dim>(),
+                                     difference_per_cell,
+                                     QGauss<dim>(fe.degree+2),
+                                     VectorTools::L2_norm,
+                                     &gradient_select);
+  const double grad_error = difference_per_cell.l2_norm();
+
+  convergence_table.add_value("cells", n_active_cells);
+  convergence_table.add_value("dofs", dof_handler.n_dofs());
+  convergence_table.add_value("val L2", L2_error);
+  convergence_table.add_value("grad L2", grad_error);
+
+  // construct post-processed solution with (hopefully) higher order of
+  // accuracy
+  QGauss<dim> quadrature(fe_u_post.degree+1);
+  FEValues<dim> fe_values(mapping, fe_u_post, quadrature,
+                          update_values | update_JxW_values |
+                          update_gradients);
+
+  const unsigned int n_q_points = quadrature.size();
+  std::vector<double> u_values(n_q_points);
+  std::vector<Tensor<1,dim> > u_gradients(n_q_points);
+  FEValuesExtractors::Vector gradients(0);
+  FEValuesExtractors::Scalar values(dim);
+  FEValues<dim> fe_values_local(mapping, fe_local, quadrature, update_values);
+  FullMatrix<double> cell_matrix(fe_u_post.dofs_per_cell,
+                                 fe_u_post.dofs_per_cell);
+  Vector<double> cell_rhs(fe_u_post.dofs_per_cell);
+  Vector<double> cell_sol(fe_u_post.dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator
+    cell_loc = dof_handler_local.begin_active(),
+    cell = dof_handler_u_post.begin_active(),
+    endc = dof_handler_u_post.end();
+  for ( ; cell != endc; ++cell, ++cell_loc)
+    {
+      fe_values.reinit(cell);
+      fe_values_local.reinit(cell_loc);
+
+      fe_values_local[values].get_function_values(solution_local, u_values);
+      fe_values_local[gradients].get_function_values(solution_local, u_gradients);
+      for (unsigned int i=1; i<fe_u_post.dofs_per_cell; ++i)
+        {
+          for (unsigned int j=0; j<fe_u_post.dofs_per_cell; ++j)
+            {
+              double sum = 0;
+              for (unsigned int q=0; q<quadrature.size(); ++q)
+                sum += (fe_values.shape_grad(i,q) *
+                        fe_values.shape_grad(j,q)
+                        ) * fe_values.JxW(q);
+              cell_matrix(i,j) = sum;
+            }
+          double sum = 0;
+          for (unsigned int q=0; q<quadrature.size(); ++q)
+            sum -= (fe_values.shape_grad(i,q) * u_gradients[q]
+                    ) * fe_values.JxW(q);
+          cell_rhs(i) = sum;
+        }
+      for (unsigned int j=0; j<fe_u_post.dofs_per_cell; ++j)
         {
           double sum = 0;
           for (unsigned int q=0; q<quadrature.size(); ++q)
-            sum += u_values[q] * fe_values.JxW(q);
-          cell_rhs(0) = sum;
+            sum += fe_values.shape_value(j,q) * fe_values.JxW(q);
+          cell_matrix(0,j) = sum;
         }
-
-        cell_matrix.gauss_jordan();
-        cell_matrix.vmult(cell_sol, cell_rhs);
-        cell->distribute_local_to_global(cell_sol, solution_u_post);
+      {
+        double sum = 0;
+        for (unsigned int q=0; q<quadrature.size(); ++q)
+          sum += u_values[q] * fe_values.JxW(q);
+        cell_rhs(0) = sum;
       }
 
-    VectorTools::integrate_difference (mapping, dof_handler_u_post,
-                                       solution_u_post,
-                                       Solution<dim>(),
-                                       difference_per_cell,
-                                       QGauss<dim>(fe.degree+3),
-                                       VectorTools::L2_norm);
-    double post_error = difference_per_cell.l2_norm();
-    convergence_table.add_value("val L2-post", post_error);
-  }
+      cell_matrix.gauss_jordan();
+      cell_matrix.vmult(cell_sol, cell_rhs);
+      cell->distribute_local_to_global(cell_sol, solution_u_post);
+    }
 
+  VectorTools::integrate_difference (mapping, dof_handler_u_post,
+                                     solution_u_post,
+                                     Solution<dim>(),
+                                     difference_per_cell,
+                                     QGauss<dim>(fe.degree+3),
+                                     VectorTools::L2_norm);
+  double post_error = difference_per_cell.l2_norm();
+  convergence_table.add_value("val L2-post", post_error);
+}
 
 
-  template <int dim>
-  void Step51<dim>::output_results (const unsigned int cycle)
-  {
-    std::string filename;
-    switch (refinement_mode)
-      {
-      case global_refinement:
-        filename = "solution-global";
-        break;
-      case adaptive_refinement:
-        filename = "solution-adaptive";
-        break;
-      default:
-        Assert (false, ExcNotImplemented());
-      }
 
-    filename += "-q" + Utilities::int_to_string(fe.degree,1);
-    filename += "-" + Utilities::int_to_string(cycle,2);
-    filename += ".vtk";
-    std::ofstream output (filename.c_str());
-
-    DataOut<dim> data_out;
-    std::vector<std::string> names (dim, "gradient");
-    names.push_back ("solution");
-    std::vector<DataComponentInterpretation::DataComponentInterpretation>
-      component_interpretation
-      (dim+1, DataComponentInterpretation::component_is_part_of_vector);
-    component_interpretation[dim]
-      = DataComponentInterpretation::component_is_scalar;
-    data_out.add_data_vector (dof_handler_local, solution_local,
-                              names, component_interpretation);
-
-    data_out.build_patches (fe.degree);
-    data_out.write_vtk (output);
-  }
+template <int dim>
+void Step51<dim>::output_results (const unsigned int cycle)
+{
+  std::string filename;
+  switch (refinement_mode)
+    {
+    case global_refinement:
+      filename = "solution-global";
+      break;
+    case adaptive_refinement:
+      filename = "solution-adaptive";
+      break;
+    default:
+      Assert (false, ExcNotImplemented());
+    }
 
+  filename += "-q" + Utilities::int_to_string(fe.degree,1);
+  filename += "-" + Utilities::int_to_string(cycle,2);
+  filename += ".vtk";
+  std::ofstream output (filename.c_str());
+
+  DataOut<dim> data_out;
+  std::vector<std::string> names (dim, "gradient");
+  names.push_back ("solution");
+  std::vector<DataComponentInterpretation::DataComponentInterpretation>
+    component_interpretation
+    (dim+1, DataComponentInterpretation::component_is_part_of_vector);
+  component_interpretation[dim]
+    = DataComponentInterpretation::component_is_scalar;
+  data_out.add_data_vector (dof_handler_local, solution_local,
+                            names, component_interpretation);
+
+  data_out.build_patches (fe.degree);
+  data_out.write_vtk (output);
+}
 
 
 
-  template <int dim>
-  void Step51<dim>::run ()
-  {
-    const bool do_cube = true;
-    if (!do_cube)
-      {
-        GridGenerator::hyper_ball (triangulation);
-        static const HyperBallBoundary<dim> boundary;
-        triangulation.set_boundary(0, boundary);
-        triangulation.refine_global(6-2*dim);
-      }
 
-    for (unsigned int cycle=0; cycle<10; ++cycle)
-      {
-        std::cout << "Cycle " << cycle << ':' << std::endl;
-
-        if (do_cube)
-          {
-            triangulation.clear();
-            GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1);
-            triangulation.refine_global(3-dim+cycle/2);
-          }
-        else triangulation.refine_global(1);
-
-        setup_system ();
-        assemble_system (false);
-        solve ();
-        postprocess();
-        output_results (cycle);
-      }
+template <int dim>
+void Step51<dim>::run ()
+{
+  const bool do_cube = true;
+  if (!do_cube)
+    {
+      GridGenerator::hyper_ball (triangulation);
+      static const HyperBallBoundary<dim> boundary;
+      triangulation.set_boundary(0, boundary);
+      triangulation.refine_global(6-2*dim);
+    }
 
+  for (unsigned int cycle=0; cycle<10; ++cycle)
+    {
+      std::cout << "Cycle " << cycle << ':' << std::endl;
 
+      if (do_cube)
+        {
+          triangulation.clear();
+          GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1);
+          triangulation.refine_global(3-dim+cycle/2);
+        }
+      else triangulation.refine_global(1);
 
-    convergence_table.set_precision("val L2", 3);
-    convergence_table.set_scientific("val L2", true);
-    convergence_table.set_precision("grad L2", 3);
-    convergence_table.set_scientific("grad L2", true);
-    convergence_table.set_precision("val L2-post", 3);
-    convergence_table.set_scientific("val L2-post", true);
+      setup_system ();
+      assemble_system (false);
+      solve ();
+      postprocess();
+      output_results (cycle);
+    }
 
-    convergence_table
-      .evaluate_convergence_rates("val L2", "cells", ConvergenceTable::reduction_rate_log2, dim);
-    convergence_table
-      .evaluate_convergence_rates("grad L2", "cells", ConvergenceTable::reduction_rate_log2, dim);
-    convergence_table
-      .evaluate_convergence_rates("val L2-post", "cells", ConvergenceTable::reduction_rate_log2, dim);
-    convergence_table.write_text(std::cout);
-  }
+
+
+  convergence_table.set_precision("val L2", 3);
+  convergence_table.set_scientific("val L2", true);
+  convergence_table.set_precision("grad L2", 3);
+  convergence_table.set_scientific("grad L2", true);
+  convergence_table.set_precision("val L2-post", 3);
+  convergence_table.set_scientific("val L2-post", true);
+
+  convergence_table
+    .evaluate_convergence_rates("val L2", "cells", ConvergenceTable::reduction_rate_log2, dim);
+  convergence_table
+    .evaluate_convergence_rates("grad L2", "cells", ConvergenceTable::reduction_rate_log2, dim);
+  convergence_table
+    .evaluate_convergence_rates("val L2-post", "cells", ConvergenceTable::reduction_rate_log2, dim);
+  convergence_table.write_text(std::cout);
 }
 
 
@@ -825,7 +823,6 @@ int main (int argc, char** argv)
   try
     {
       using namespace dealii;
-      using namespace Step51;
 
       deallog.depth_console (0);
 

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.