]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
First version of new tutorial example demonstrating hybridizable discontinuous Galerk...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 6 Aug 2013 13:30:49 +0000 (13:30 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 6 Aug 2013 13:30:49 +0000 (13:30 +0000)
git-svn-id: https://svn.dealii.org/trunk@30236 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-51/CMakeLists.txt [new file with mode: 0644]
deal.II/examples/step-51/doc/builds-on [new file with mode: 0644]
deal.II/examples/step-51/doc/intro.dox [new file with mode: 0644]
deal.II/examples/step-51/doc/kind [new file with mode: 0644]
deal.II/examples/step-51/doc/results.dox [new file with mode: 0644]
deal.II/examples/step-51/doc/tooltip [new file with mode: 0644]
deal.II/examples/step-51/step-51.cc [new file with mode: 0644]

diff --git a/deal.II/examples/step-51/CMakeLists.txt b/deal.II/examples/step-51/CMakeLists.txt
new file mode 100644 (file)
index 0000000..6d7622f
--- /dev/null
@@ -0,0 +1,31 @@
+##
+#  CMake script for the step-7 tutorial program:
+##
+
+# Set the name of the project and target:
+SET(TARGET "step-51")
+
+# Declare all source files the target consists of:
+SET(TARGET_SRC
+  ${TARGET}.cc
+  # You can specify additional files here!
+  )
+
+# Usually, you will not need to modify anything beyond this point...
+
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
+
+FIND_PACKAGE(deal.II 8.0 QUIET
+  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+  )
+IF(NOT ${deal.II_FOUND})
+  MESSAGE(FATAL_ERROR "\n"
+    "*** Could not locate deal.II. ***\n\n"
+    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+    )
+ENDIF()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()
diff --git a/deal.II/examples/step-51/doc/builds-on b/deal.II/examples/step-51/doc/builds-on
new file mode 100644 (file)
index 0000000..943473d
--- /dev/null
@@ -0,0 +1 @@
+step-7
diff --git a/deal.II/examples/step-51/doc/intro.dox b/deal.II/examples/step-51/doc/intro.dox
new file mode 100644 (file)
index 0000000..a5458f7
--- /dev/null
@@ -0,0 +1,5 @@
+<a name="Intro"></a>
+<h1>Introduction</h1>
+
+This tutorial program presents the implementation of a hybridizable
+discontinuous Galkerin method for the convection-diffusion equation.
\ No newline at end of file
diff --git a/deal.II/examples/step-51/doc/kind b/deal.II/examples/step-51/doc/kind
new file mode 100644 (file)
index 0000000..c1d9154
--- /dev/null
@@ -0,0 +1 @@
+techniques
diff --git a/deal.II/examples/step-51/doc/results.dox b/deal.II/examples/step-51/doc/results.dox
new file mode 100644 (file)
index 0000000..ae4bf1c
--- /dev/null
@@ -0,0 +1,139 @@
+<h1>Results</h1>
+
+
+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
+@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.
+
+
+
diff --git a/deal.II/examples/step-51/doc/tooltip b/deal.II/examples/step-51/doc/tooltip
new file mode 100644 (file)
index 0000000..784e257
--- /dev/null
@@ -0,0 +1 @@
+Convection-diffusion equation. Hybridizable discontinuous Galerkin methods. Face elements.
diff --git a/deal.II/examples/step-51/step-51.cc b/deal.II/examples/step-51/step-51.cc
new file mode 100644 (file)
index 0000000..256e48d
--- /dev/null
@@ -0,0 +1,893 @@
+/* ---------------------------------------------------------------------
+ * $Id$
+ *
+ * Copyright (C) 2013 - 2013 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE at
+ * the top level of the deal.II distribution.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ * Author: Martin Kronbichler, TU Muenchen,
+ *         Scott T. Miller, xxx, 2013
+ */
+
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/tensor_function.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/convergence_table.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_face.h>
+
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <deal.II/numerics/data_out_faces.h>
+
+namespace Step51
+{
+  using namespace dealii;
+
+  // @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;
+  };
+
+
+  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<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>
+  class ConvectionVelocity : public TensorFunction<1,dim>
+  {
+  public:
+    ConvectionVelocity() : TensorFunction<1,dim>() {}
+
+    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>
+  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 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));
+      }
+
+    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;
+
+    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 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) {}
+
+    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>() {}
+
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
+
+  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)));
+      }
+
+    return return_value / Utilities::fixed_power<dim>(std::sqrt(2 * numbers::PI)
+                                                      * this->width);
+  }
+
+
+
+  template <int dim>
+  class Step51
+  {
+  public:
+    enum RefinementMode
+    {
+      global_refinement, adaptive_refinement
+    };
+
+    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);
+
+    Triangulation<dim>   triangulation;
+
+    const MappingQ<dim>  mapping;
+
+    FESystem<dim>        fe_local;
+    DoFHandler<dim>      dof_handler_local;
+
+    FE_FaceQ<dim>        fe;
+    DoFHandler<dim>      dof_handler;
+
+    FE_DGQ<dim>          fe_u_post;
+    DoFHandler<dim>      dof_handler_u_post;
+
+    ConstraintMatrix     constraints;
+    ChunkSparsityPattern sparsity_pattern;
+    ChunkSparseMatrix<double> system_matrix;
+
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+
+    Vector<double>       solution_local;
+    Vector<double>       solution_u_post;
+
+    const RefinementMode refinement_mode;
+
+    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>
+  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 ();
+
+    {
+      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);
+  }
+
+
+
+  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)
+                    {
+                      double sum_lf = 0., sum_fl = 0.;
+                      for (unsigned int q=0; q<n_face_q_points; ++q)
+                        {
+                          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];
+                        }
+                      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))
+                  {
+                    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;
+                      }
+                  }
+              }
+            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)
+                  {
+                    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;
+                  }
+              }
+          }
+
+        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);
+          }
+      }
+  }
+
+
+
+  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;
+
+    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;
+          }
+        {
+          double sum = 0;
+          for (unsigned int q=0; q<quadrature.size(); ++q)
+            sum += u_values[q] * fe_values.JxW(q);
+          cell_rhs(0) = sum;
+        }
+
+        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>::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);
+      }
+
+
+
+    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);
+  }
+}
+
+
+int main (int argc, char** argv)
+{
+  const unsigned int dim = 2;
+
+  try
+    {
+      using namespace dealii;
+      using namespace Step51;
+
+      deallog.depth_console (0);
+
+      // Now for the three calls to the main class in complete analogy to
+      // step-7.
+      {
+        std::cout << "Solving with Q1 elements, adaptive refinement" << std::endl
+                  << "=============================================" << std::endl
+                  << std::endl;
+
+        Step51<dim> hdg_problem (1, Step51<dim>::adaptive_refinement);
+        hdg_problem.run ();
+
+        std::cout << std::endl;
+      }
+
+      {
+        std::cout << "Solving with Q1 elements, global refinement" << std::endl
+                  << "===========================================" << std::endl
+                  << std::endl;
+
+        Step51<dim> hdg_problem (1, Step51<dim>::global_refinement);
+        hdg_problem.run ();
+
+        std::cout << std::endl;
+      }
+
+      {
+        std::cout << "Solving with Q3 elements, global refinement" << std::endl
+                  << "===========================================" << std::endl
+                  << std::endl;
+
+        Step51<dim> hdg_problem (3, Step51<dim>::global_refinement);
+        hdg_problem.run ();
+
+        std::cout << std::endl;
+      }
+
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+
+  return 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.