]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Also create 1D version of FE_FaceQ/P to allow step-51 to run in 1D (in principle...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 11 Feb 2014 08:33:28 +0000 (08:33 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 11 Feb 2014 08:33:28 +0000 (08:33 +0000)
git-svn-id: https://svn.dealii.org/trunk@32454 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/fe/fe_face.h
deal.II/include/deal.II/fe/fe_poly_face.templates.h
deal.II/source/fe/fe_face.cc
deal.II/source/fe/fe_face.inst.in
tests/bits/step-51.cc [new file with mode: 0644]
tests/bits/step-51.output [new file with mode: 0644]
tests/bits/step-51p.cc [new file with mode: 0644]
tests/bits/step-51p.output [new file with mode: 0644]
tests/fe/fe_data_test.cc

index 6a1402798326ec3a90981f293ca1e30aa571bfbb..e62fcd2e49b8341267dc364e132a1fc6c212c7c5 100644 (file)
@@ -124,6 +124,16 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li>New: FE_FaceQ and FE_FaceP now also work in 1D (with a single dof
+  on each vertex).
+  <br>
+  (Martin Kronbichler, 2014/02/11)
+
+  <li>Fixed: FE_DGQ::has_support_on_face returned a wrong number for element
+  degree larger than 1 in 1D. This is now fixed.
+  <br>
+  (Martin Kronbichler, 2014/02/10)
+
   <li>Changed: DerivativeApproximation used to be a class that only had
   static members. It is now a namespace.
   <br>
index 6e806085489863f2681e12ea63f36381bdd0f0db..0831942687342191325c3137863eca94061af449 100644 (file)
@@ -124,7 +124,6 @@ public:
   FiniteElementDomination::Domination
   compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const;
 
-
   /**
    * Returns a list of constant modes of the element. For this element, it
    * simply returns one row with all entries set to true.
@@ -140,6 +139,197 @@ private:
 
 
 
+/**
+ * Specialization of FE_FaceQ for 1D. In that case, the finite element only
+ * consists of one degree of freedom in each of the two faces (= vertices) of
+ * a cell, irrespective of the degree. However, this element still accepts a
+ * degree in its constructor and also returns that degree. This way,
+ * dimension-independent programming with trace elements is also possible in
+ * 1D (even though there is no computational benefit at all from it in 1D).
+ *
+ * @ingroup fe
+ * @author Guido Kanschat, Martin Kronbichler
+ * @date 2014
+ */
+template <int spacedim>
+class FE_FaceQ<1,spacedim> : public FiniteElement<1,spacedim>
+{
+public:
+  /**
+   * Constructor.
+   */
+  FE_FaceQ (const unsigned int p);
+
+  /**
+   * Clone method.
+   */
+  virtual FiniteElement<1,spacedim> *clone() const;
+
+  /**
+   * Return a string that uniquely identifies a finite element. This class
+   * returns <tt>FE_FaceQ<dim>(degree)</tt>, with <tt>dim</tt> and
+   * <tt>degree</tt> replaced by appropriate values.
+   */
+  virtual std::string get_name () const;
+
+  /**
+   * Return the matrix interpolating from a face of of one element to the face
+   * of the neighboring element.  The size of the matrix is then
+   * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. This
+   * element only provides interpolation matrices for elements of the same
+   * type and FE_Nothing. For all other elements, an exception of type
+   * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is thrown.
+   */
+  virtual void
+  get_face_interpolation_matrix (const FiniteElement<1,spacedim> &source,
+                                 FullMatrix<double>       &matrix) const;
+
+  /**
+   * Return the matrix interpolating from a face of of one element to the face
+   * of the neighboring element.  The size of the matrix is then
+   * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. This
+   * element only provides interpolation matrices for elements of the same
+   * type and FE_Nothing. For all other elements, an exception of type
+   * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is thrown.
+   */
+  virtual void
+  get_subface_interpolation_matrix (const FiniteElement<1,spacedim> &source,
+                                    const unsigned int        subface,
+                                    FullMatrix<double>       &matrix) const;
+
+  /**
+   * Check for non-zero values on a face.
+   *
+   * This function returns @p true, if the shape function @p shape_index has
+   * non-zero values on the face @p face_index.
+   *
+   * Implementation of the interface in FiniteElement
+   */
+  virtual bool has_support_on_face (const unsigned int shape_index,
+                                    const unsigned int face_index) const;
+
+  /**
+   * Return whether this element implements its hanging node constraints in
+   * the new way, which has to be used to make elements "hp compatible".
+   */
+  virtual bool hp_constraints_are_implemented () const;
+
+  /**
+   * Return whether this element dominates the one given as argument when they
+   * meet at a common face, whether it is the other way around, whether
+   * neither dominates, or if either could dominate.
+   *
+   * For a definition of domination, see FiniteElementBase::Domination and in
+   * particular the @ref hp_paper "hp paper".
+   */
+  virtual
+  FiniteElementDomination::Domination
+  compare_for_face_domination (const FiniteElement<1,spacedim> &fe_other) const;
+
+  /**
+   * Returns a list of constant modes of the element. For this element, it
+   * simply returns one row with all entries set to true.
+   */
+  virtual Table<2,bool> get_constant_modes () const;
+
+protected:
+  virtual
+  typename Mapping<1,spacedim>::InternalDataBase *
+  get_data (const UpdateFlags,
+            const Mapping<1,spacedim> &mapping,
+            const Quadrature<1> &quadrature) const ;
+
+  typename Mapping<1,spacedim>::InternalDataBase *
+  get_face_data (const UpdateFlags,
+                 const Mapping<1,spacedim> &mapping,
+                 const Quadrature<0>& quadrature) const ;
+
+  typename Mapping<1,spacedim>::InternalDataBase *
+  get_subface_data (const UpdateFlags,
+                    const Mapping<1,spacedim> &mapping,
+                    const Quadrature<0>& quadrature) const ;
+
+  virtual void
+  fill_fe_values (const Mapping<1,spacedim>                           &mapping,
+                  const typename Triangulation<1,spacedim>::cell_iterator &cell,
+                  const Quadrature<1>                                 &quadrature,
+                  typename Mapping<1,spacedim>::InternalDataBase      &mapping_internal,
+                  typename Mapping<1,spacedim>::InternalDataBase      &fe_internal,
+                  FEValuesData<1,spacedim>                            &data,
+                  CellSimilarity::Similarity                       &cell_similarity) const;
+
+  virtual void
+  fill_fe_face_values (const Mapping<1,spacedim> &mapping,
+                       const typename Triangulation<1,spacedim>::cell_iterator &cell,
+                       const unsigned int                    face_no,
+                       const Quadrature<0>                &quadrature,
+                       typename Mapping<1,spacedim>::InternalDataBase      &mapping_internal,
+                       typename Mapping<1,spacedim>::InternalDataBase      &fe_internal,
+                       FEValuesData<1,spacedim> &data) const ;
+
+  virtual void
+  fill_fe_subface_values (const Mapping<1,spacedim> &mapping,
+                          const typename Triangulation<1,spacedim>::cell_iterator &cell,
+                          const unsigned int                    face_no,
+                          const unsigned int                    sub_no,
+                          const Quadrature<0>                &quadrature,
+                          typename Mapping<1,spacedim>::InternalDataBase      &mapping_internal,
+                          typename Mapping<1,spacedim>::InternalDataBase      &fe_internal,
+                          FEValuesData<1,spacedim> &data) const ;
+
+
+  /**
+   * Determine the values that need to be computed on the unit cell to be able
+   * to compute all values required by <tt>flags</tt>.
+   *
+   * For the purpuse of this function, refer to the documentation in
+   * FiniteElement.
+   *
+   * This class assumes that shape functions of this FiniteElement do
+   * <em>not</em> depend on the actual shape of the cells in real
+   * space. Therefore, the effect in this element is as follows: if
+   * <tt>update_values</tt> is set in <tt>flags</tt>, copy it to the
+   * result. All other flags of the result are cleared, since everything else
+   * must be computed for each cell.
+   */
+  virtual UpdateFlags update_once (const UpdateFlags flags) const;
+
+  /**
+   * Determine the values that need to be computed on every cell to be able to
+   * compute all values required by <tt>flags</tt>.
+   *
+   * For the purpuse of this function, refer to the documentation in
+   * FiniteElement.
+   *
+   * This class assumes that shape functions of this FiniteElement do
+   * <em>not</em> depend on the actual shape of the cells in real space.
+   *
+   * The effect in this element is as follows:
+   * <ul>
+   * <li> if <tt>update_gradients</tt> is set, the result will contain
+   * <tt>update_gradients</tt> and <tt>update_covariant_transformation</tt>.
+   * The latter is required to transform the gradient on the unit cell to the
+   * real cell. Remark, that the action required by
+   * <tt>update_covariant_transformation</tt> is actually performed by the
+   * Mapping object used in conjunction with this finite element.  <li> if
+   * <tt>update_hessians</tt> is set, the result will contain
+   * <tt>update_hessians</tt> and <tt>update_covariant_transformation</tt>.
+   * The rationale is the same as above and no higher derivatives of the
+   * transformation are required, since we use difference quotients for the
+   * actual computation.
+   *
+   * </ul>
+   */
+  virtual UpdateFlags update_each (const UpdateFlags flags) const;
+
+private:
+  /**
+   * Return vector with dofs per vertex, line, quad, hex.
+   */
+  static std::vector<unsigned int> get_dpo_vector (const unsigned int deg);
+};
+
+
 
 /**
  * A finite element, which is a Legendre element of complete polynomials on
@@ -257,6 +447,26 @@ private:
 };
 
 
+
+/**
+ * FE_FaceP in 1D, i.e., with degrees of freedom on the element vertices.
+ */
+template <int spacedim>
+class FE_FaceP<1,spacedim> : public FE_FaceQ<1,spacedim>
+{
+public:
+  /**
+   * Constructor.
+   */
+  FE_FaceP (const unsigned int p);
+
+  /**
+   * Returns the name of the element
+   */
+  std::string get_name() const;
+};
+
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
index 39651c017005e086f248228dced633713b7f20df..c1067d752180f38fe9ffa76d8d2406c396d99570 100644 (file)
@@ -56,11 +56,8 @@ template <class POLY, int dim, int spacedim>
 UpdateFlags
 FE_PolyFace<POLY,dim,spacedim>::update_once (const UpdateFlags) const
 {
-  // for this kind of elements, only
-  // the values can be precomputed
-  // once and for all. set this flag
-  // if the values are requested at
-  // all
+  // for this kind of elements, only the values can be precomputed once and
+  // for all. set this flag if the values are requested at all
   return update_default;
 }
 
index 728a6d24980aad63f81075629c72437e20131fab..9d62ff19d54edb7b79c31de16557adc957b70565 100644 (file)
@@ -300,6 +300,270 @@ FE_FaceQ<dim,spacedim>::get_constant_modes () const
 
 
 
+// ----------------------------- FE_FaceQ<1,spacedim> ------------------------
+
+template <int spacedim>
+FE_FaceQ<1,spacedim>::FE_FaceQ (const unsigned int degree)
+  :
+  FiniteElement<1,spacedim> (FiniteElementData<1>(get_dpo_vector(degree), 1, degree, FiniteElementData<1>::L2),
+                             std::vector<bool>(1,true),
+                             std::vector<ComponentMask> (1, ComponentMask(1,true)))
+{
+  this->unit_face_support_points.resize(1);
+
+  // initialize unit support points (this makes it possible to assign initial
+  // values to FE_FaceQ)
+  this->unit_support_points.resize(GeometryInfo<1>::faces_per_cell);
+  this->unit_support_points[1] = Point<1>(1.);
+}
+
+
+
+template <int spacedim>
+FiniteElement<1,spacedim> *
+FE_FaceQ<1,spacedim>::clone() const
+{
+  return new FE_FaceQ<1,spacedim>(this->degree);
+}
+
+
+
+template <int spacedim>
+std::string
+FE_FaceQ<1,spacedim>::get_name () const
+{
+  // note that the FETools::get_fe_from_name function depends on the
+  // particular format of the string this function returns, so they have to be
+  // kept in synch
+  std::ostringstream namebuf;
+  namebuf << "FE_FaceQ<1>(" << this->degree << ")";
+
+  return namebuf.str();
+}
+
+
+
+template <int spacedim>
+void
+FE_FaceQ<1,spacedim>::
+get_face_interpolation_matrix (const FiniteElement<1,spacedim> &source_fe,
+                               FullMatrix<double>       &interpolation_matrix) const
+{
+  get_subface_interpolation_matrix (source_fe, numbers::invalid_unsigned_int,
+                                    interpolation_matrix);
+}
+
+
+
+template <int spacedim>
+void
+FE_FaceQ<1,spacedim>::
+get_subface_interpolation_matrix (const FiniteElement<1,spacedim> &x_source_fe,
+                                  const unsigned int        subface,
+                                  FullMatrix<double>       &interpolation_matrix) const
+{
+  Assert (interpolation_matrix.n() == this->dofs_per_face,
+          ExcDimensionMismatch (interpolation_matrix.n(),
+                                this->dofs_per_face));
+  Assert (interpolation_matrix.m() == x_source_fe.dofs_per_face,
+          ExcDimensionMismatch (interpolation_matrix.m(),
+                                x_source_fe.dofs_per_face));
+  interpolation_matrix(0,0) = 1.;
+}
+
+
+
+template <int spacedim>
+bool
+FE_FaceQ<1,spacedim>::has_support_on_face (
+  const unsigned int shape_index,
+  const unsigned int face_index) const
+{
+  AssertIndexRange(shape_index, 2);
+  return (face_index == shape_index);
+}
+
+
+
+template <int spacedim>
+std::vector<unsigned int>
+FE_FaceQ<1,spacedim>::get_dpo_vector (const unsigned int)
+{
+  std::vector<unsigned int> dpo(2, 0U);
+  dpo[0] = 1;
+  return dpo;
+}
+
+
+
+template <int spacedim>
+bool
+FE_FaceQ<1,spacedim>::hp_constraints_are_implemented () const
+{
+  return true;
+}
+
+
+
+template <int spacedim>
+FiniteElementDomination::Domination
+FE_FaceQ<1,spacedim>::
+compare_for_face_domination (const FiniteElement<1,spacedim> &fe_other) const
+{
+  return FiniteElementDomination::no_requirements;
+}
+
+
+
+template <int spacedim>
+Table<2,bool>
+FE_FaceQ<1,spacedim>::get_constant_modes () const
+{
+  Table<2,bool> constant_modes(1, this->dofs_per_cell);
+  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+    constant_modes(0,i) = true;
+  return constant_modes;
+}
+
+
+
+template <int spacedim>
+UpdateFlags
+FE_FaceQ<1,spacedim>::update_once (const UpdateFlags) const
+{
+  return update_default;
+}
+
+
+
+template <int spacedim>
+UpdateFlags
+FE_FaceQ<1,spacedim>::update_each (const UpdateFlags flags) const
+{
+  UpdateFlags out = flags & update_values;
+  if (flags & update_gradients)
+    out |= update_gradients | update_covariant_transformation;
+  if (flags & update_hessians)
+    out |= update_hessians | update_covariant_transformation;
+  if (flags & update_cell_normal_vectors)
+    out |= update_cell_normal_vectors | update_JxW_values;
+
+  return out;
+}
+
+
+
+template <int spacedim>
+typename Mapping<1,spacedim>::InternalDataBase *
+FE_FaceQ<1,spacedim>::get_data (
+  const UpdateFlags,
+  const Mapping<1,spacedim> &,
+  const Quadrature<1> &) const
+{
+  return new typename Mapping<1,spacedim>::InternalDataBase;
+}
+
+
+template <int spacedim>
+typename Mapping<1,spacedim>::InternalDataBase *
+FE_FaceQ<1,spacedim>::get_face_data (
+  const UpdateFlags update_flags,
+  const Mapping<1,spacedim> &,
+  const Quadrature<0>& quadrature) const
+{
+  // generate a new data object and initialize some fields
+  typename Mapping<1,spacedim>::InternalDataBase *data =
+    new typename Mapping<1,spacedim>::InternalDataBase;
+
+  // check what needs to be initialized only once and what on every
+  // cell/face/subface we visit
+  data->update_once = update_once(update_flags);
+  data->update_each = update_each(update_flags);
+  data->update_flags = data->update_once | data->update_each;
+
+  const UpdateFlags flags(data->update_flags);
+  const unsigned int n_q_points = quadrature.size();
+  AssertDimension(n_q_points, 1);
+
+  // No derivatives of this element are implemented.
+  if (flags & update_gradients || flags & update_hessians)
+    {
+      Assert(false, ExcNotImplemented());
+    }
+
+  return data;
+}
+
+
+
+template <int spacedim>
+typename Mapping<1,spacedim>::InternalDataBase *
+FE_FaceQ<1,spacedim>::get_subface_data (
+  const UpdateFlags flags,
+  const Mapping<1,spacedim> &mapping,
+  const Quadrature<0>& quadrature) const
+{
+  return get_face_data (flags, mapping, quadrature);
+}
+
+
+
+template <int spacedim>
+void
+FE_FaceQ<1,spacedim>::fill_fe_values
+(const Mapping<1,spacedim> &,
+ const typename Triangulation<1,spacedim>::cell_iterator &,
+ const Quadrature<1> &,
+ typename Mapping<1,spacedim>::InternalDataBase &,
+ typename Mapping<1,spacedim>::InternalDataBase &,
+ FEValuesData<1,spacedim> &,
+ CellSimilarity::Similarity &) const
+{
+  // Do nothing, since we do not have values in the interior
+}
+
+
+
+template <int spacedim>
+void
+FE_FaceQ<1,spacedim>::fill_fe_face_values (
+  const Mapping<1,spacedim> &,
+  const typename Triangulation<1,spacedim>::cell_iterator &,
+  const unsigned int face,
+  const Quadrature<0>& quadrature,
+  typename Mapping<1,spacedim>::InternalDataBase &,
+  typename Mapping<1,spacedim>::InternalDataBase &fedata,
+  FEValuesData<1,spacedim> &data) const
+{
+  const UpdateFlags flags(fedata.update_once | fedata.update_each);
+
+  const unsigned int foffset = face;
+  if (flags & update_values)
+    {
+      for (unsigned int k=0; k<this->dofs_per_cell; ++k)
+        data.shape_values(k,0) = 0.;
+      data.shape_values(foffset,0) = 1;
+    }
+}
+
+
+template <int spacedim>
+void
+FE_FaceQ<1,spacedim>::fill_fe_subface_values (
+  const Mapping<1,spacedim> &,
+  const typename Triangulation<1,spacedim>::cell_iterator &,
+  const unsigned int ,
+  const unsigned int ,
+  const Quadrature<0>& ,
+  typename Mapping<1,spacedim>::InternalDataBase &,
+  typename Mapping<1,spacedim>::InternalDataBase &,
+  FEValuesData<1,spacedim> &) const
+{
+  Assert(false, ExcMessage("Should not fille subface values in 1D"));
+}
+
+
+
 // --------------------------------------- FE_FaceP --------------------------
 
 template <int dim, int spacedim>
@@ -312,6 +576,7 @@ FE_FaceP<dim,spacedim>::FE_FaceP (const unsigned int degree)
 {}
 
 
+
 template <int dim, int spacedim>
 FiniteElement<dim,spacedim> *
 FE_FaceP<dim,spacedim>::clone() const
@@ -320,6 +585,7 @@ FE_FaceP<dim,spacedim>::clone() const
 }
 
 
+
 template <int dim, int spacedim>
 std::string
 FE_FaceP<dim,spacedim>::get_name () const
@@ -525,6 +791,29 @@ FE_FaceP<dim,spacedim>::get_constant_modes () const
 
 
 
+template <int spacedim>
+FE_FaceP<1,spacedim>::FE_FaceP (const unsigned int degree)
+  :
+  FE_FaceQ<1,spacedim> (degree)
+{}
+
+
+
+template <int spacedim>
+std::string
+FE_FaceP<1,spacedim>::get_name () const
+{
+  // note that the FETools::get_fe_from_name function depends on the
+  // particular format of the string this function returns, so they have to be
+  // kept in synch
+  std::ostringstream namebuf;
+  namebuf << "FE_FaceP<1>(" << this->degree << ")";
+
+  return namebuf.str();
+}
+
+
+
 // explicit instantiations
 #include "fe_face.inst"
 
index 513eb0030bf9060d0f2f0028a595eff2da8f39b2..c1af4014419ae9f52a0f034c107cf2eef85102f1 100644 (file)
@@ -22,8 +22,8 @@ for (deal_II_dimension : DIMENSIONS)
     template class FE_PolyFace<TensorProductPolynomials<deal_II_dimension-1> >;
     template class FE_PolyFace<PolynomialSpace<deal_II_dimension-1>, deal_II_dimension>;
     //template class FE_PolyFace<PolynomialsP<deal_II_dimension>, deal_II_dimension>;
-    template class FE_FaceQ<deal_II_dimension>;
-    template class FE_FaceP<deal_II_dimension>;
 #endif
+    template class FE_FaceQ<deal_II_dimension,deal_II_dimension>;
+    template class FE_FaceP<deal_II_dimension,deal_II_dimension>;
   }
 
diff --git a/tests/bits/step-51.cc b/tests/bits/step-51.cc
new file mode 100644 (file)
index 0000000..ed203d6
--- /dev/null
@@ -0,0 +1,1005 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2005 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// tests step-51 (without output, without ChunkSparseMatrix, without
+// adaptivity) in 1d and 2d
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <fstream>
+std::ofstream logfile("output");
+
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/convergence_table.h>
+#include <deal.II/base/work_stream.h>
+#include <deal.II/base/tensor_function.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/solver_bicgstab.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_face.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iomanip>
+#include <fstream>
+#include <cmath>
+
+
+namespace Step51
+{
+
+  using namespace dealii;
+
+  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 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;
+  };
+
+  template <int dim>
+  void SolutionAndGradient<dim>::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 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 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 HDG
+  {
+  public:
+    HDG (const unsigned int degree);
+    void run ();
+
+  private:
+
+    void setup_system ();
+    void assemble_system (const bool reconstruct_trace = false);
+    void solve ();
+    void postprocess ();
+    void refine_grid (const unsigned int cylce);
+    void output_results (const unsigned int cycle);
+
+    struct PerTaskData;
+    struct ScratchData;
+
+    struct PostProcessScratchData;
+
+    void assemble_system_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+                                   ScratchData &scratch,
+                                   PerTaskData &task_data);
+
+    void copy_local_to_global(const PerTaskData &data);
+
+    void postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+                               PostProcessScratchData &scratch,
+                               unsigned int &empty_data);
+
+
+    Triangulation<dim>   triangulation;
+
+    FESystem<dim>        fe_local;
+    DoFHandler<dim>      dof_handler_local;
+    Vector<double>       solution_local;
+
+    FE_FaceQ<dim>        fe;
+    DoFHandler<dim>      dof_handler;
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+
+    FE_DGQ<dim>          fe_u_post;
+    DoFHandler<dim>      dof_handler_u_post;
+    Vector<double>       solution_u_post;
+
+    ConstraintMatrix     constraints;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+    ConvergenceTable     convergence_table;
+  };
+
+
+  template <int dim>
+  HDG<dim>::HDG (const unsigned int degree) :
+    fe_local (FE_DGQ<dim>(degree), dim+1),
+    dof_handler_local (triangulation),
+    fe (degree),
+    dof_handler (triangulation),
+    fe_u_post (degree+1),
+    dof_handler_u_post (triangulation)
+  {}
+
+
+
+  template <int dim>
+  void
+  HDG<dim>::setup_system ()
+  {
+    dof_handler_local.distribute_dofs(fe_local);
+    dof_handler.distribute_dofs(fe);
+    dof_handler_u_post.distribute_dofs(fe_u_post);
+
+    deallog << "   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);
+    typename FunctionMap<dim>::type boundary_functions;
+    Solution<dim> solution_function;
+    boundary_functions[0] = &solution_function;
+    VectorTools::project_boundary_values (dof_handler,
+                                          boundary_functions,
+                                          QGauss<dim-1>(fe.degree+1),
+                                          constraints);
+    constraints.close ();
+
+    {
+      CompressedSimpleSparsityPattern csp (dof_handler.n_dofs());
+      DoFTools::make_sparsity_pattern (dof_handler, csp,
+                                       constraints, false);
+      sparsity_pattern.copy_from(csp);
+    }
+    system_matrix.reinit (sparsity_pattern);
+  }
+
+
+
+  template <int dim>
+  struct HDG<dim>::PerTaskData
+  {
+    FullMatrix<double> cell_matrix;
+    Vector<double>     cell_vector;
+    std::vector<types::global_dof_index> dof_indices;
+
+    bool trace_reconstruct;
+
+    PerTaskData(const unsigned int n_dofs, const bool trace_reconstruct)
+      : cell_matrix(n_dofs, n_dofs),
+        cell_vector(n_dofs),
+        dof_indices(n_dofs),
+        trace_reconstruct(trace_reconstruct)
+    {}
+  };
+
+
+
+  template <int dim>
+  struct HDG<dim>::ScratchData
+  {
+    FEValues<dim>     fe_values_local;
+    FEFaceValues<dim> fe_face_values_local;
+    FEFaceValues<dim> fe_face_values;
+
+    FullMatrix<double> ll_matrix;
+    FullMatrix<double> lf_matrix;
+    FullMatrix<double> fl_matrix;
+    FullMatrix<double> tmp_matrix;
+    Vector<double>     l_rhs;
+    Vector<double>     tmp_rhs;
+
+    std::vector<Tensor<1,dim> > q_phi;
+    std::vector<double>         q_phi_div;
+    std::vector<double>         u_phi;
+    std::vector<Tensor<1,dim> > u_phi_grad;
+    std::vector<double>         tr_phi;
+    std::vector<double>         trace_values;
+
+    std::vector<std::vector<unsigned int> > fe_local_support_on_face;
+    std::vector<std::vector<unsigned int> > fe_support_on_face;
+
+    ConvectionVelocity<dim> convection_velocity;
+    RightHandSide<dim> right_hand_side;
+    const Solution<dim> exact_solution;
+
+    ScratchData(const FiniteElement<dim> &fe,
+                const FiniteElement<dim> &fe_local,
+                const QGauss<dim>   &quadrature_formula,
+                const QGauss<dim-1> &face_quadrature_formula,
+                const UpdateFlags local_flags,
+                const UpdateFlags local_face_flags,
+                const UpdateFlags flags)
+      :
+      fe_values_local (fe_local, quadrature_formula, local_flags),
+      fe_face_values_local (fe_local, face_quadrature_formula, local_face_flags),
+      fe_face_values (fe, face_quadrature_formula, flags),
+      ll_matrix (fe_local.dofs_per_cell, fe_local.dofs_per_cell),
+      lf_matrix (fe_local.dofs_per_cell, fe.dofs_per_cell),
+      fl_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell),
+      tmp_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell),
+      l_rhs (fe_local.dofs_per_cell),
+      tmp_rhs (fe_local.dofs_per_cell),
+      q_phi (fe_local.dofs_per_cell),
+      q_phi_div (fe_local.dofs_per_cell),
+      u_phi (fe_local.dofs_per_cell),
+      u_phi_grad (fe_local.dofs_per_cell),
+      tr_phi (fe.dofs_per_cell),
+      trace_values(face_quadrature_formula.size()),
+      fe_local_support_on_face(GeometryInfo<dim>::faces_per_cell),
+      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<fe_local.dofs_per_cell; ++i)
+          {
+            if (fe_local.has_support_on_face(i,face))
+              {
+                fe_local_support_on_face[face].push_back(i);
+              }
+          }
+
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+        for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+          {
+            if (fe.has_support_on_face(i,face))
+              {
+                fe_support_on_face[face].push_back(i);
+              }
+          }
+    }
+
+    ScratchData(const ScratchData &sd)
+      :
+      fe_values_local (sd.fe_values_local.get_fe(),
+                       sd.fe_values_local.get_quadrature(),
+                       sd.fe_values_local.get_update_flags()),
+      fe_face_values_local (sd.fe_face_values_local.get_fe(),
+                            sd.fe_face_values_local.get_quadrature(),
+                            sd.fe_face_values_local.get_update_flags()),
+      fe_face_values (sd.fe_face_values.get_fe(),
+                      sd.fe_face_values.get_quadrature(),
+                      sd.fe_face_values.get_update_flags()),
+      ll_matrix (sd.ll_matrix),
+      lf_matrix (sd.lf_matrix),
+      fl_matrix (sd.fl_matrix),
+      tmp_matrix (sd.tmp_matrix),
+      l_rhs (sd.l_rhs),
+      tmp_rhs (sd.tmp_rhs),
+      q_phi (sd.q_phi),
+      q_phi_div (sd.q_phi_div),
+      u_phi (sd.u_phi),
+      u_phi_grad (sd.u_phi_grad),
+      tr_phi (sd.tr_phi),
+      trace_values(sd.trace_values),
+      fe_local_support_on_face(sd.fe_local_support_on_face),
+      fe_support_on_face(sd.fe_support_on_face)
+    {}
+  };
+
+
+
+  template <int dim>
+  struct HDG<dim>::PostProcessScratchData
+  {
+    FEValues<dim> fe_values_local;
+    FEValues<dim> fe_values;
+
+    std::vector<double> u_values;
+    std::vector<Tensor<1,dim> > u_gradients;
+    FullMatrix<double> cell_matrix;
+
+    Vector<double> cell_rhs;
+    Vector<double> cell_sol;
+
+    PostProcessScratchData(const FiniteElement<dim> &fe,
+                           const FiniteElement<dim> &fe_local,
+                           const QGauss<dim>   &quadrature_formula,
+                           const UpdateFlags local_flags,
+                           const UpdateFlags flags)
+      :
+      fe_values_local (fe_local, quadrature_formula, local_flags),
+      fe_values (fe, quadrature_formula, flags),
+      u_values (quadrature_formula.size()),
+      u_gradients (quadrature_formula.size()),
+      cell_matrix (fe.dofs_per_cell, fe.dofs_per_cell),
+      cell_rhs (fe.dofs_per_cell),
+      cell_sol (fe.dofs_per_cell)
+    {}
+
+    PostProcessScratchData(const PostProcessScratchData &sd)
+      :
+      fe_values_local (sd.fe_values_local.get_fe(),
+                       sd.fe_values_local.get_quadrature(),
+                       sd.fe_values_local.get_update_flags()),
+      fe_values (sd.fe_values.get_fe(),
+                 sd.fe_values.get_quadrature(),
+                 sd.fe_values.get_update_flags()),
+      u_values (sd.u_values),
+      u_gradients (sd.u_gradients),
+      cell_matrix (sd.cell_matrix),
+      cell_rhs (sd.cell_rhs),
+      cell_sol (sd.cell_sol)
+    {}
+  };
+
+
+
+  template <int dim>
+  void
+  HDG<dim>::assemble_system (const bool trace_reconstruct)
+  {
+    const QGauss<dim>   quadrature_formula(fe.degree+1);
+    const QGauss<dim-1> face_quadrature_formula(fe.degree+1);
+
+    const UpdateFlags local_flags (update_values | update_gradients |
+                                   update_JxW_values | update_quadrature_points);
+
+    const UpdateFlags local_face_flags (update_values);
+
+    const UpdateFlags flags ( update_values | update_normal_vectors |
+                              update_quadrature_points |
+                              update_JxW_values);
+
+    PerTaskData task_data (fe.dofs_per_cell,
+                           trace_reconstruct);
+    ScratchData scratch (fe, fe_local,
+                         quadrature_formula,
+                         face_quadrature_formula,
+                         local_flags,
+                         local_face_flags,
+                         flags);
+
+    WorkStream::run(dof_handler.begin_active(),
+                    dof_handler.end(),
+                    *this,
+                    &HDG<dim>::assemble_system_one_cell,
+                    &HDG<dim>::copy_local_to_global,
+                    scratch,
+                    task_data);
+  }
+
+
+
+  template <int dim>
+  void
+  HDG<dim>::assemble_system_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+                                      ScratchData &scratch,
+                                      PerTaskData &task_data)
+  {
+    typename DoFHandler<dim>::active_cell_iterator
+    loc_cell (&triangulation,
+              cell->level(),
+              cell->index(),
+              &dof_handler_local);
+
+    const unsigned int n_q_points    = scratch.fe_values_local.get_quadrature().size();
+    const unsigned int n_face_q_points = scratch.fe_face_values_local.get_quadrature().size();
+
+    const unsigned int loc_dofs_per_cell = scratch.fe_values_local.get_fe().dofs_per_cell;
+
+    const FEValuesExtractors::Vector fluxes (0);
+    const FEValuesExtractors::Scalar scalar (dim);
+
+    scratch.ll_matrix = 0;
+    scratch.l_rhs = 0;
+    if (!task_data.trace_reconstruct)
+      {
+        scratch.lf_matrix = 0;
+        scratch.fl_matrix = 0;
+        task_data.cell_matrix = 0;
+        task_data.cell_vector = 0;
+      }
+    scratch.fe_values_local.reinit (loc_cell);
+
+    for (unsigned int q=0; q<n_q_points; ++q)
+      {
+        const double rhs_value
+        = scratch.right_hand_side.value(scratch.fe_values_local.quadrature_point(q));
+        const Tensor<1,dim> convection
+        = scratch.convection_velocity.value(scratch.fe_values_local.quadrature_point(q));
+        const double JxW = scratch.fe_values_local.JxW(q);
+        for (unsigned int k=0; k<loc_dofs_per_cell; ++k)
+          {
+            scratch.q_phi[k] = scratch.fe_values_local[fluxes].value(k,q);
+            scratch.q_phi_div[k] = scratch.fe_values_local[fluxes].divergence(k,q);
+            scratch.u_phi[k] = scratch.fe_values_local[scalar].value(k,q);
+            scratch.u_phi_grad[k] = scratch.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)
+              scratch.ll_matrix(i,j) += (
+                                          scratch.q_phi[i] * scratch.q_phi[j]
+                                          -
+                                          scratch.q_phi_div[i] * scratch.u_phi[j]
+                                          +
+                                          scratch.u_phi[i] * scratch.q_phi_div[j]
+                                          -
+                                          (scratch.u_phi_grad[i] * convection) * scratch.u_phi[j]
+                                        ) * JxW;
+            scratch.l_rhs(i) += scratch.u_phi[i] * rhs_value * JxW;
+          }
+      }
+
+    for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+      {
+        scratch.fe_face_values_local.reinit(loc_cell, face);
+        scratch.fe_face_values.reinit(cell, face);
+
+        if (task_data.trace_reconstruct)
+          scratch.fe_face_values.get_function_values (solution, scratch.trace_values);
+
+        for (unsigned int q=0; q<n_face_q_points; ++q)
+          {
+            const double JxW = scratch.fe_face_values.JxW(q);
+            const Point<dim> quadrature_point =
+              scratch.fe_face_values.quadrature_point(q);
+            const Point<dim> normal = scratch.fe_face_values.normal_vector(q);
+            const Tensor<1,dim> convection
+            = scratch.convection_velocity.value(quadrature_point);
+
+            const double tau_stab = (5. +
+                                     std::abs(convection * normal));
+
+            for (unsigned int k=0; k<scratch.fe_local_support_on_face[face].size(); ++k)
+              {
+                const unsigned int kk=scratch.fe_local_support_on_face[face][k];
+                scratch.q_phi[k] = scratch.fe_face_values_local[fluxes].value(kk,q);
+                scratch.u_phi[k] = scratch.fe_face_values_local[scalar].value(kk,q);
+              }
+
+            if (!task_data.trace_reconstruct)
+              {
+                for (unsigned int k=0; k<scratch.fe_support_on_face[face].size(); ++k)
+                  scratch.tr_phi[k] =
+                    scratch.fe_face_values.shape_value(scratch.fe_support_on_face[face][k],q);
+                for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
+                  for (unsigned int j=0; j<scratch.fe_support_on_face[face].size(); ++j)
+                    {
+                      const unsigned int ii=scratch.fe_local_support_on_face[face][i];
+                      const unsigned int jj=scratch.fe_support_on_face[face][j];
+                      scratch.lf_matrix(ii,jj) += (
+                                                    (scratch.q_phi[i] * normal
+                                                     +
+                                                     (convection * normal -
+                                                      tau_stab) * scratch.u_phi[i])
+                                                    * scratch.tr_phi[j]
+                                                  ) * JxW;
+
+                      scratch.fl_matrix(jj,ii) -= (
+                                                    (scratch.q_phi[i] * normal
+                                                     +
+                                                     tau_stab * scratch.u_phi[i])
+                                                    * scratch.tr_phi[j]
+                                                  ) * JxW;
+                    }
+
+                for (unsigned int i=0; i<scratch.fe_support_on_face[face].size(); ++i)
+                  for (unsigned int j=0; j<scratch.fe_support_on_face[face].size(); ++j)
+                    {
+                      const unsigned int ii=scratch.fe_support_on_face[face][i];
+                      const unsigned int jj=scratch.fe_support_on_face[face][j];
+                      task_data.cell_matrix(ii,jj) += (
+                                                        (convection * normal - tau_stab) *
+                                                        scratch.tr_phi[i] * scratch.tr_phi[j]
+                                                      ) * JxW;
+                    }
+
+                if (cell->face(face)->at_boundary()
+                    &&
+                    (cell->face(face)->boundary_indicator() == 1))
+                  {
+                    const double neumann_value =
+                      - scratch.exact_solution.gradient (quadrature_point) * normal
+                      + convection * normal * scratch.exact_solution.value(quadrature_point);
+                    for (unsigned int i=0; i<scratch.fe_support_on_face[face].size(); ++i)
+                      {
+                        const unsigned int ii=scratch.fe_support_on_face[face][i];
+                        task_data.cell_vector(ii) += scratch.tr_phi[i] * neumann_value * JxW;
+                      }
+                  }
+              }
+
+            for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
+              for (unsigned int j=0; j<scratch.fe_local_support_on_face[face].size(); ++j)
+                {
+                  const unsigned int ii=scratch.fe_local_support_on_face[face][i];
+                  const unsigned int jj=scratch.fe_local_support_on_face[face][j];
+                  scratch.ll_matrix(ii,jj) += tau_stab * scratch.u_phi[i] * scratch.u_phi[j] * JxW;
+                }
+  
+            if (task_data.trace_reconstruct)
+              for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
+                {
+                  const unsigned int ii=scratch.fe_local_support_on_face[face][i];
+                  scratch.l_rhs(ii) -= (scratch.q_phi[i] * normal
+                                        +
+                                        scratch.u_phi[i] * (convection * normal - tau_stab)
+                                       ) * scratch.trace_values[q] * JxW;
+                }
+          }
+      }
+
+    scratch.ll_matrix.gauss_jordan();
+
+    if (task_data.trace_reconstruct == false)
+      {
+        scratch.fl_matrix.mmult(scratch.tmp_matrix, scratch.ll_matrix);
+        scratch.tmp_matrix.vmult_add(task_data.cell_vector, scratch.l_rhs);
+        scratch.tmp_matrix.mmult(task_data.cell_matrix, scratch.lf_matrix, true);
+        cell->get_dof_indices(task_data.dof_indices);
+      }
+    else
+      {
+        scratch.ll_matrix.vmult(scratch.tmp_rhs, scratch.l_rhs);
+        loc_cell->set_dof_values(scratch.tmp_rhs, solution_local);
+      }
+  }
+
+
+
+  template <int dim>
+  void HDG<dim>::copy_local_to_global(const PerTaskData &data)
+  {
+    if (data.trace_reconstruct == false)
+      constraints.distribute_local_to_global (data.cell_matrix,
+                                              data.cell_vector,
+                                              data.dof_indices,
+                                              system_matrix, system_rhs);
+  }
+
+
+
+  template <int dim>
+  void HDG<dim>::solve ()
+  {
+    SolverControl solver_control (system_matrix.m()*10,
+                                  1e-11*system_rhs.l2_norm());
+    SolverBicgstab<> solver (solver_control, false);
+    solver.solve (system_matrix, solution, system_rhs,
+                  PreconditionIdentity());
+
+    system_matrix.clear();
+
+    constraints.distribute(solution);
+
+    assemble_system(true);
+  }
+
+
+
+
+  template <int dim>
+  void
+  HDG<dim>::postprocess()
+  {
+    {
+      const QGauss<dim>   quadrature_formula(fe_u_post.degree+1);
+      const UpdateFlags local_flags (update_values);
+      const UpdateFlags flags ( update_values | update_gradients |
+                                update_JxW_values);
+
+      PostProcessScratchData scratch (fe_u_post, fe_local,
+                                      quadrature_formula,
+                                      local_flags,
+                                      flags);
+
+      WorkStream::run(dof_handler_u_post.begin_active(),
+                      dof_handler_u_post.end(),
+                      std_cxx1x::bind (&HDG<dim>::postprocess_one_cell,
+                                       std_cxx1x::ref(*this),
+                                       std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3),
+                      std_cxx1x::function<void(const unsigned int &)>(),
+                      scratch,
+                      0U);
+    }
+
+    Vector<float> difference_per_cell (triangulation.n_active_cells());
+
+    ComponentSelectFunction<dim> value_select (dim, dim+1);
+    VectorTools::integrate_difference (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 (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();
+
+    VectorTools::integrate_difference (dof_handler_u_post,
+                                       solution_u_post,
+                                       Solution<dim>(),
+                                       difference_per_cell,
+                                       QGauss<dim>(fe.degree+3),
+                                       VectorTools::L2_norm);
+    const double post_error = difference_per_cell.l2_norm();
+
+    convergence_table.add_value("cells",     triangulation.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);
+    convergence_table.add_value("val L2-post", post_error);
+  }
+
+
+
+  template <int dim>
+  void
+  HDG<dim>::postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+                                  PostProcessScratchData &scratch,
+                                  unsigned int &)
+  {
+    typename DoFHandler<dim>::active_cell_iterator
+    loc_cell (&triangulation,
+              cell->level(),
+              cell->index(),
+              &dof_handler_local);
+
+    scratch.fe_values_local.reinit (loc_cell);
+    scratch.fe_values.reinit(cell);
+
+    FEValuesExtractors::Vector fluxes(0);
+    FEValuesExtractors::Scalar scalar(dim);
+
+    const unsigned int n_q_points = scratch.fe_values.get_quadrature().size();
+    const unsigned int dofs_per_cell = scratch.fe_values.dofs_per_cell;
+
+    scratch.fe_values_local[scalar].get_function_values(solution_local, scratch.u_values);
+    scratch.fe_values_local[fluxes].get_function_values(solution_local, scratch.u_gradients);
+
+    double sum = 0;
+    for (unsigned int i=1; i<dofs_per_cell; ++i)
+      {
+        for (unsigned int j=0; j<dofs_per_cell; ++j)
+          {
+            sum = 0;
+            for (unsigned int q=0; q<n_q_points; ++q)
+              sum += (scratch.fe_values.shape_grad(i,q) *
+                      scratch.fe_values.shape_grad(j,q)
+                     ) * scratch.fe_values.JxW(q);
+            scratch.cell_matrix(i,j) = sum;
+          }
+
+        sum = 0;
+        for (unsigned int q=0; q<n_q_points; ++q)
+          sum -= (scratch.fe_values.shape_grad(i,q) * scratch.u_gradients[q]
+                 ) * scratch.fe_values.JxW(q);
+        scratch.cell_rhs(i) = sum;
+      }
+    for (unsigned int j=0; j<dofs_per_cell; ++j)
+      {
+        sum = 0;
+        for (unsigned int q=0; q<n_q_points; ++q)
+          sum += scratch.fe_values.shape_value(j,q) * scratch.fe_values.JxW(q);
+        scratch.cell_matrix(0,j) = sum;
+      }
+    {
+      sum = 0;
+      for (unsigned int q=0; q<n_q_points; ++q)
+        sum += scratch.u_values[q] * scratch.fe_values.JxW(q);
+      scratch.cell_rhs(0) = sum;
+    }
+
+    scratch.cell_matrix.gauss_jordan();
+    scratch.cell_matrix.vmult(scratch.cell_sol, scratch.cell_rhs);
+    cell->distribute_local_to_global(scratch.cell_sol, solution_u_post);
+  }
+
+
+
+  template <int dim>
+  void HDG<dim>::output_results (const unsigned int cycle)
+  {
+    // not included in test
+  }
+
+
+
+  template <int dim>
+  void HDG<dim>::refine_grid (const unsigned int cycle)
+  {
+    // only global refinement
+    triangulation.clear();
+    GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1);
+    triangulation.refine_global(3-dim+cycle/2);
+
+    typename Triangulation<dim>::cell_iterator
+    cell = triangulation.begin (),
+    endc = triangulation.end();
+    for (; cell!=endc; ++cell)
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+        if (cell->face(face)->at_boundary())
+          for (unsigned int d=0; d<dim; ++d)
+            if ((std::fabs(cell->face(face)->center()(d) - (1)) < 1e-12))
+              cell->face(face)->set_boundary_indicator (1);
+  }
+
+
+
+  template <int dim>
+  void HDG<dim>::run ()
+  {
+    for (unsigned int cycle=0; cycle<14-4*dim; ++cycle)
+      {
+        deallog << "Cycle " << cycle << ':' << std::endl;
+
+        refine_grid (cycle);
+        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(deallog.get_file_stream());
+  }
+
+} // end of namespace Step51
+
+
+
+
+int main ()
+{
+  deallog << std::setprecision(6);
+  logfile << std::setprecision(6);
+
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  {
+    Step51::HDG<1> hdg_problem (1);
+    hdg_problem.run ();
+  }
+
+  {
+    Step51::HDG<1> hdg_problem (3);
+    hdg_problem.run ();
+  }
+
+  {
+    Step51::HDG<2> hdg_problem (1);
+    hdg_problem.run ();
+  }
+
+  {
+    Step51::HDG<2> hdg_problem (2);
+    hdg_problem.run ();
+  }
+
+  return 0;
+}
diff --git a/tests/bits/step-51.output b/tests/bits/step-51.output
new file mode 100644 (file)
index 0000000..768c5da
--- /dev/null
@@ -0,0 +1,165 @@
+
+DEAL::Cycle 0:
+DEAL::   Number of degrees of freedom: 9
+DEAL:Bicgstab::Starting value 8.19443
+DEAL:Bicgstab::Convergence step 8 value 0
+DEAL::Cycle 1:
+DEAL::   Number of degrees of freedom: 13
+DEAL:Bicgstab::Starting value 13.9539
+DEAL:Bicgstab::Convergence step 12 value 0
+DEAL::Cycle 2:
+DEAL::   Number of degrees of freedom: 17
+DEAL:Bicgstab::Starting value 13.0731
+DEAL:Bicgstab::Convergence step 16 value 0
+DEAL::Cycle 3:
+DEAL::   Number of degrees of freedom: 25
+DEAL:Bicgstab::Starting value 12.1492
+DEAL:Bicgstab::Convergence step 24 value 0
+DEAL::Cycle 4:
+DEAL::   Number of degrees of freedom: 33
+DEAL:Bicgstab::Starting value 11.0492
+DEAL:Bicgstab::Convergence step 34 value 0
+DEAL::Cycle 5:
+DEAL::   Number of degrees of freedom: 49
+DEAL:Bicgstab::Starting value 9.35310
+DEAL:Bicgstab::Convergence step 51 value 0
+DEAL::Cycle 6:
+DEAL::   Number of degrees of freedom: 65
+DEAL:Bicgstab::Starting value 8.20475
+DEAL:Bicgstab::Convergence step 74 value 0
+DEAL::Cycle 7:
+DEAL::   Number of degrees of freedom: 97
+DEAL:Bicgstab::Starting value 6.76136
+DEAL:Bicgstab::Convergence step 104 value 0
+DEAL::Cycle 8:
+DEAL::   Number of degrees of freedom: 129
+DEAL:Bicgstab::Starting value 5.87456
+DEAL:Bicgstab::Convergence step 137 value 0
+DEAL::Cycle 9:
+DEAL::   Number of degrees of freedom: 193
+DEAL:Bicgstab::Starting value 4.80772
+DEAL:Bicgstab::Convergence step 213 value 0
+cells dofs     val L2        grad L2      val L2-post   
+8     9    4.080e-01 -    1.722e+00 -    3.362e-01 -    
+12    13   1.659e-01 2.22 8.766e-01 1.67 2.151e-02 6.78 
+16    17   7.996e-02 2.54 4.531e-01 2.29 1.125e-02 2.25 
+24    25   3.705e-02 1.90 2.028e-01 1.98 3.372e-03 2.97 
+32    33   2.125e-02 1.93 1.140e-01 2.00 1.430e-03 2.98 
+48    49   9.610e-03 1.96 5.048e-02 2.01 4.242e-04 3.00 
+64    65   5.449e-03 1.97 2.831e-02 2.01 1.787e-04 3.00 
+96    97   2.440e-03 1.98 1.253e-02 2.01 5.281e-05 3.01 
+128   129  1.377e-03 1.99 7.035e-03 2.01 2.224e-05 3.01 
+192   193  6.141e-04 1.99 3.119e-03 2.01 6.573e-06 3.01 
+DEAL::Cycle 0:
+DEAL::   Number of degrees of freedom: 9
+DEAL:Bicgstab::Starting value 7.96477
+DEAL:Bicgstab::Convergence step 8 value 0
+DEAL::Cycle 1:
+DEAL::   Number of degrees of freedom: 13
+DEAL:Bicgstab::Starting value 13.4076
+DEAL:Bicgstab::Convergence step 12 value 0
+DEAL::Cycle 2:
+DEAL::   Number of degrees of freedom: 17
+DEAL:Bicgstab::Starting value 12.9549
+DEAL:Bicgstab::Convergence step 16 value 0
+DEAL::Cycle 3:
+DEAL::   Number of degrees of freedom: 25
+DEAL:Bicgstab::Starting value 12.1266
+DEAL:Bicgstab::Convergence step 24 value 0
+DEAL::Cycle 4:
+DEAL::   Number of degrees of freedom: 33
+DEAL:Bicgstab::Starting value 11.0426
+DEAL:Bicgstab::Convergence step 34 value 0
+DEAL::Cycle 5:
+DEAL::   Number of degrees of freedom: 49
+DEAL:Bicgstab::Starting value 9.35200
+DEAL:Bicgstab::Convergence step 50 value 0
+DEAL::Cycle 6:
+DEAL::   Number of degrees of freedom: 65
+DEAL:Bicgstab::Starting value 8.20444
+DEAL:Bicgstab::Convergence step 73 value 0
+DEAL::Cycle 7:
+DEAL::   Number of degrees of freedom: 97
+DEAL:Bicgstab::Starting value 6.76131
+DEAL:Bicgstab::Convergence step 105 value 0
+DEAL::Cycle 8:
+DEAL::   Number of degrees of freedom: 129
+DEAL:Bicgstab::Starting value 5.87455
+DEAL:Bicgstab::Convergence step 143 value 0
+DEAL::Cycle 9:
+DEAL::   Number of degrees of freedom: 193
+DEAL:Bicgstab::Starting value 4.80772
+DEAL:Bicgstab::Convergence step 225 value 0
+cells dofs     val L2        grad L2      val L2-post   
+8     9    2.141e-02 -    1.398e-01 -    3.487e-03 -    
+12    13   6.701e-03 2.87 3.811e-02 3.21 4.976e-04 4.80 
+16    17   1.693e-03 4.78 9.999e-03 4.65 1.085e-04 5.29 
+24    25   3.501e-04 3.89 2.009e-03 3.96 1.440e-05 4.98 
+32    33   1.126e-04 3.94 6.373e-04 3.99 3.412e-06 5.01 
+48    49   2.256e-05 3.97 1.259e-04 4.00 4.472e-07 5.01 
+64    65   7.179e-06 3.98 3.980e-05 4.00 1.057e-07 5.01 
+96    97   1.425e-06 3.99 7.848e-06 4.00 1.388e-08 5.01 
+128   129  4.519e-07 3.99 2.481e-06 4.00 3.302e-09 4.99 
+192   193  8.946e-08 3.99 4.894e-07 4.00 7.398e-10 3.69 
+DEAL::Cycle 0:
+DEAL::   Number of degrees of freedom: 80
+DEAL:Bicgstab::Starting value 12.2150
+DEAL:Bicgstab::Convergence step 37 value 0
+DEAL::Cycle 1:
+DEAL::   Number of degrees of freedom: 168
+DEAL:Bicgstab::Starting value 8.50278
+DEAL:Bicgstab::Convergence step 48 value 0
+DEAL::Cycle 2:
+DEAL::   Number of degrees of freedom: 288
+DEAL:Bicgstab::Starting value 10.5333
+DEAL:Bicgstab::Convergence step 62 value 1.01190e-10
+DEAL::Cycle 3:
+DEAL::   Number of degrees of freedom: 624
+DEAL:Bicgstab::Starting value 8.10478
+DEAL:Bicgstab::Convergence step 98 value 0
+DEAL::Cycle 4:
+DEAL::   Number of degrees of freedom: 1088
+DEAL:Bicgstab::Starting value 6.61947
+DEAL:Bicgstab::Convergence step 120 value 0
+DEAL::Cycle 5:
+DEAL::   Number of degrees of freedom: 2400
+DEAL:Bicgstab::Starting value 4.75284
+DEAL:Bicgstab::Convergence step 181 value 0
+cells dofs     val L2        grad L2      val L2-post   
+16    80   1.066e+01 -    1.571e+01 -    1.055e+01 -    
+36    168  3.727e+00 2.59 7.998e+00 1.67 3.660e+00 2.61 
+64    288  8.086e-01 5.31 4.228e+00 2.22 4.492e-01 7.29 
+144   624  2.730e-01 2.68 1.867e+00 2.02 6.110e-02 4.92 
+256   1088 1.493e-01 2.10 1.046e+00 2.01 2.878e-02 2.62 
+576   2400 6.965e-02 1.88 4.847e-01 1.90 9.202e-03 2.81 
+DEAL::Cycle 0:
+DEAL::   Number of degrees of freedom: 120
+DEAL:Bicgstab::Starting value 10.6959
+DEAL:Bicgstab::Convergence step 43 value 0
+DEAL::Cycle 1:
+DEAL::   Number of degrees of freedom: 252
+DEAL:Bicgstab::Starting value 13.2620
+DEAL:Bicgstab::Convergence step 58 value 0
+DEAL::Cycle 2:
+DEAL::   Number of degrees of freedom: 432
+DEAL:Bicgstab::Starting value 11.2567
+DEAL:Bicgstab::Convergence step 74 value 0
+DEAL::Cycle 3:
+DEAL::   Number of degrees of freedom: 936
+DEAL:Bicgstab::Starting value 9.46542
+DEAL:Bicgstab::Convergence step 105 value 0
+DEAL::Cycle 4:
+DEAL::   Number of degrees of freedom: 1632
+DEAL:Bicgstab::Starting value 7.48176
+DEAL:Bicgstab::Convergence step 138 value 0
+DEAL::Cycle 5:
+DEAL::   Number of degrees of freedom: 3600
+DEAL:Bicgstab::Starting value 5.13012
+DEAL:Bicgstab::Convergence step 189 value 0
+cells dofs     val L2        grad L2       val L2-post   
+16    120  7.645e-01 -    5.413e+00 -    2.667e-01 -     
+36    252  4.745e-01 1.18 2.375e+00 2.03 3.501e-01 -0.67 
+64    432  9.635e-02 5.54 7.772e-01 3.88 5.032e-02 6.74  
+144   936  3.487e-02 2.51 2.718e-01 2.59 5.339e-03 5.53  
+256   1632 1.920e-02 2.08 1.390e-01 2.33 1.731e-03 3.92  
+576   3600 6.021e-03 2.86 4.291e-02 2.90 3.431e-04 3.99  
diff --git a/tests/bits/step-51p.cc b/tests/bits/step-51p.cc
new file mode 100644 (file)
index 0000000..fae43af
--- /dev/null
@@ -0,0 +1,1005 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2005 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// tests step-51 (without output, without ChunkSparseMatrix, without
+// adaptivity) in 1d and 2d with FE_DGP/FE_FaceP
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <fstream>
+std::ofstream logfile("output");
+
+
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/convergence_table.h>
+#include <deal.II/base/work_stream.h>
+#include <deal.II/base/tensor_function.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/solver_bicgstab.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_face.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iomanip>
+#include <fstream>
+#include <cmath>
+
+
+namespace Step51
+{
+
+  using namespace dealii;
+
+  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 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;
+  };
+
+  template <int dim>
+  void SolutionAndGradient<dim>::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 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 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 HDG
+  {
+  public:
+    HDG (const unsigned int degree);
+    void run ();
+
+  private:
+
+    void setup_system ();
+    void assemble_system (const bool reconstruct_trace = false);
+    void solve ();
+    void postprocess ();
+    void refine_grid (const unsigned int cylce);
+    void output_results (const unsigned int cycle);
+
+    struct PerTaskData;
+    struct ScratchData;
+
+    struct PostProcessScratchData;
+
+    void assemble_system_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+                                   ScratchData &scratch,
+                                   PerTaskData &task_data);
+
+    void copy_local_to_global(const PerTaskData &data);
+
+    void postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+                               PostProcessScratchData &scratch,
+                               unsigned int &empty_data);
+
+
+    Triangulation<dim>   triangulation;
+
+    FESystem<dim>        fe_local;
+    DoFHandler<dim>      dof_handler_local;
+    Vector<double>       solution_local;
+
+    FE_FaceP<dim>        fe;
+    DoFHandler<dim>      dof_handler;
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+
+    FE_DGP<dim>          fe_u_post;
+    DoFHandler<dim>      dof_handler_u_post;
+    Vector<double>       solution_u_post;
+
+    ConstraintMatrix     constraints;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+    ConvergenceTable     convergence_table;
+  };
+
+
+  template <int dim>
+  HDG<dim>::HDG (const unsigned int degree) :
+    fe_local (FE_DGP<dim>(degree), dim+1),
+    dof_handler_local (triangulation),
+    fe (degree),
+    dof_handler (triangulation),
+    fe_u_post (degree+1),
+    dof_handler_u_post (triangulation)
+  {}
+
+
+
+  template <int dim>
+  void
+  HDG<dim>::setup_system ()
+  {
+    dof_handler_local.distribute_dofs(fe_local);
+    dof_handler.distribute_dofs(fe);
+    dof_handler_u_post.distribute_dofs(fe_u_post);
+
+    deallog << "   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);
+    typename FunctionMap<dim>::type boundary_functions;
+    Solution<dim> solution_function;
+    boundary_functions[0] = &solution_function;
+    VectorTools::project_boundary_values (dof_handler,
+                                          boundary_functions,
+                                          QGauss<dim-1>(fe.degree+1),
+                                          constraints);
+    constraints.close ();
+
+    {
+      CompressedSimpleSparsityPattern csp (dof_handler.n_dofs());
+      DoFTools::make_sparsity_pattern (dof_handler, csp,
+                                       constraints, false);
+      sparsity_pattern.copy_from(csp);
+    }
+    system_matrix.reinit (sparsity_pattern);
+  }
+
+
+
+  template <int dim>
+  struct HDG<dim>::PerTaskData
+  {
+    FullMatrix<double> cell_matrix;
+    Vector<double>     cell_vector;
+    std::vector<types::global_dof_index> dof_indices;
+
+    bool trace_reconstruct;
+
+    PerTaskData(const unsigned int n_dofs, const bool trace_reconstruct)
+      : cell_matrix(n_dofs, n_dofs),
+        cell_vector(n_dofs),
+        dof_indices(n_dofs),
+        trace_reconstruct(trace_reconstruct)
+    {}
+  };
+
+
+
+  template <int dim>
+  struct HDG<dim>::ScratchData
+  {
+    FEValues<dim>     fe_values_local;
+    FEFaceValues<dim> fe_face_values_local;
+    FEFaceValues<dim> fe_face_values;
+
+    FullMatrix<double> ll_matrix;
+    FullMatrix<double> lf_matrix;
+    FullMatrix<double> fl_matrix;
+    FullMatrix<double> tmp_matrix;
+    Vector<double>     l_rhs;
+    Vector<double>     tmp_rhs;
+
+    std::vector<Tensor<1,dim> > q_phi;
+    std::vector<double>         q_phi_div;
+    std::vector<double>         u_phi;
+    std::vector<Tensor<1,dim> > u_phi_grad;
+    std::vector<double>         tr_phi;
+    std::vector<double>         trace_values;
+
+    std::vector<std::vector<unsigned int> > fe_local_support_on_face;
+    std::vector<std::vector<unsigned int> > fe_support_on_face;
+
+    ConvectionVelocity<dim> convection_velocity;
+    RightHandSide<dim> right_hand_side;
+    const Solution<dim> exact_solution;
+
+    ScratchData(const FiniteElement<dim> &fe,
+                const FiniteElement<dim> &fe_local,
+                const QGauss<dim>   &quadrature_formula,
+                const QGauss<dim-1> &face_quadrature_formula,
+                const UpdateFlags local_flags,
+                const UpdateFlags local_face_flags,
+                const UpdateFlags flags)
+      :
+      fe_values_local (fe_local, quadrature_formula, local_flags),
+      fe_face_values_local (fe_local, face_quadrature_formula, local_face_flags),
+      fe_face_values (fe, face_quadrature_formula, flags),
+      ll_matrix (fe_local.dofs_per_cell, fe_local.dofs_per_cell),
+      lf_matrix (fe_local.dofs_per_cell, fe.dofs_per_cell),
+      fl_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell),
+      tmp_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell),
+      l_rhs (fe_local.dofs_per_cell),
+      tmp_rhs (fe_local.dofs_per_cell),
+      q_phi (fe_local.dofs_per_cell),
+      q_phi_div (fe_local.dofs_per_cell),
+      u_phi (fe_local.dofs_per_cell),
+      u_phi_grad (fe_local.dofs_per_cell),
+      tr_phi (fe.dofs_per_cell),
+      trace_values(face_quadrature_formula.size()),
+      fe_local_support_on_face(GeometryInfo<dim>::faces_per_cell),
+      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<fe_local.dofs_per_cell; ++i)
+          {
+            if (fe_local.has_support_on_face(i,face))
+              {
+                fe_local_support_on_face[face].push_back(i);
+              }
+          }
+
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+        for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+          {
+            if (fe.has_support_on_face(i,face))
+              {
+                fe_support_on_face[face].push_back(i);
+              }
+          }
+    }
+
+    ScratchData(const ScratchData &sd)
+      :
+      fe_values_local (sd.fe_values_local.get_fe(),
+                       sd.fe_values_local.get_quadrature(),
+                       sd.fe_values_local.get_update_flags()),
+      fe_face_values_local (sd.fe_face_values_local.get_fe(),
+                            sd.fe_face_values_local.get_quadrature(),
+                            sd.fe_face_values_local.get_update_flags()),
+      fe_face_values (sd.fe_face_values.get_fe(),
+                      sd.fe_face_values.get_quadrature(),
+                      sd.fe_face_values.get_update_flags()),
+      ll_matrix (sd.ll_matrix),
+      lf_matrix (sd.lf_matrix),
+      fl_matrix (sd.fl_matrix),
+      tmp_matrix (sd.tmp_matrix),
+      l_rhs (sd.l_rhs),
+      tmp_rhs (sd.tmp_rhs),
+      q_phi (sd.q_phi),
+      q_phi_div (sd.q_phi_div),
+      u_phi (sd.u_phi),
+      u_phi_grad (sd.u_phi_grad),
+      tr_phi (sd.tr_phi),
+      trace_values(sd.trace_values),
+      fe_local_support_on_face(sd.fe_local_support_on_face),
+      fe_support_on_face(sd.fe_support_on_face)
+    {}
+  };
+
+
+
+  template <int dim>
+  struct HDG<dim>::PostProcessScratchData
+  {
+    FEValues<dim> fe_values_local;
+    FEValues<dim> fe_values;
+
+    std::vector<double> u_values;
+    std::vector<Tensor<1,dim> > u_gradients;
+    FullMatrix<double> cell_matrix;
+
+    Vector<double> cell_rhs;
+    Vector<double> cell_sol;
+
+    PostProcessScratchData(const FiniteElement<dim> &fe,
+                           const FiniteElement<dim> &fe_local,
+                           const QGauss<dim>   &quadrature_formula,
+                           const UpdateFlags local_flags,
+                           const UpdateFlags flags)
+      :
+      fe_values_local (fe_local, quadrature_formula, local_flags),
+      fe_values (fe, quadrature_formula, flags),
+      u_values (quadrature_formula.size()),
+      u_gradients (quadrature_formula.size()),
+      cell_matrix (fe.dofs_per_cell, fe.dofs_per_cell),
+      cell_rhs (fe.dofs_per_cell),
+      cell_sol (fe.dofs_per_cell)
+    {}
+
+    PostProcessScratchData(const PostProcessScratchData &sd)
+      :
+      fe_values_local (sd.fe_values_local.get_fe(),
+                       sd.fe_values_local.get_quadrature(),
+                       sd.fe_values_local.get_update_flags()),
+      fe_values (sd.fe_values.get_fe(),
+                 sd.fe_values.get_quadrature(),
+                 sd.fe_values.get_update_flags()),
+      u_values (sd.u_values),
+      u_gradients (sd.u_gradients),
+      cell_matrix (sd.cell_matrix),
+      cell_rhs (sd.cell_rhs),
+      cell_sol (sd.cell_sol)
+    {}
+  };
+
+
+
+  template <int dim>
+  void
+  HDG<dim>::assemble_system (const bool trace_reconstruct)
+  {
+    const QGauss<dim>   quadrature_formula(fe.degree+1);
+    const QGauss<dim-1> face_quadrature_formula(fe.degree+1);
+
+    const UpdateFlags local_flags (update_values | update_gradients |
+                                   update_JxW_values | update_quadrature_points);
+
+    const UpdateFlags local_face_flags (update_values);
+
+    const UpdateFlags flags ( update_values | update_normal_vectors |
+                              update_quadrature_points |
+                              update_JxW_values);
+
+    PerTaskData task_data (fe.dofs_per_cell,
+                           trace_reconstruct);
+    ScratchData scratch (fe, fe_local,
+                         quadrature_formula,
+                         face_quadrature_formula,
+                         local_flags,
+                         local_face_flags,
+                         flags);
+
+    WorkStream::run(dof_handler.begin_active(),
+                    dof_handler.end(),
+                    *this,
+                    &HDG<dim>::assemble_system_one_cell,
+                    &HDG<dim>::copy_local_to_global,
+                    scratch,
+                    task_data);
+  }
+
+
+
+  template <int dim>
+  void
+  HDG<dim>::assemble_system_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+                                      ScratchData &scratch,
+                                      PerTaskData &task_data)
+  {
+    typename DoFHandler<dim>::active_cell_iterator
+    loc_cell (&triangulation,
+              cell->level(),
+              cell->index(),
+              &dof_handler_local);
+
+    const unsigned int n_q_points    = scratch.fe_values_local.get_quadrature().size();
+    const unsigned int n_face_q_points = scratch.fe_face_values_local.get_quadrature().size();
+
+    const unsigned int loc_dofs_per_cell = scratch.fe_values_local.get_fe().dofs_per_cell;
+
+    const FEValuesExtractors::Vector fluxes (0);
+    const FEValuesExtractors::Scalar scalar (dim);
+
+    scratch.ll_matrix = 0;
+    scratch.l_rhs = 0;
+    if (!task_data.trace_reconstruct)
+      {
+        scratch.lf_matrix = 0;
+        scratch.fl_matrix = 0;
+        task_data.cell_matrix = 0;
+        task_data.cell_vector = 0;
+      }
+    scratch.fe_values_local.reinit (loc_cell);
+
+    for (unsigned int q=0; q<n_q_points; ++q)
+      {
+        const double rhs_value
+        = scratch.right_hand_side.value(scratch.fe_values_local.quadrature_point(q));
+        const Tensor<1,dim> convection
+        = scratch.convection_velocity.value(scratch.fe_values_local.quadrature_point(q));
+        const double JxW = scratch.fe_values_local.JxW(q);
+        for (unsigned int k=0; k<loc_dofs_per_cell; ++k)
+          {
+            scratch.q_phi[k] = scratch.fe_values_local[fluxes].value(k,q);
+            scratch.q_phi_div[k] = scratch.fe_values_local[fluxes].divergence(k,q);
+            scratch.u_phi[k] = scratch.fe_values_local[scalar].value(k,q);
+            scratch.u_phi_grad[k] = scratch.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)
+              scratch.ll_matrix(i,j) += (
+                                          scratch.q_phi[i] * scratch.q_phi[j]
+                                          -
+                                          scratch.q_phi_div[i] * scratch.u_phi[j]
+                                          +
+                                          scratch.u_phi[i] * scratch.q_phi_div[j]
+                                          -
+                                          (scratch.u_phi_grad[i] * convection) * scratch.u_phi[j]
+                                        ) * JxW;
+            scratch.l_rhs(i) += scratch.u_phi[i] * rhs_value * JxW;
+          }
+      }
+
+    for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+      {
+        scratch.fe_face_values_local.reinit(loc_cell, face);
+        scratch.fe_face_values.reinit(cell, face);
+
+        if (task_data.trace_reconstruct)
+          scratch.fe_face_values.get_function_values (solution, scratch.trace_values);
+
+        for (unsigned int q=0; q<n_face_q_points; ++q)
+          {
+            const double JxW = scratch.fe_face_values.JxW(q);
+            const Point<dim> quadrature_point =
+              scratch.fe_face_values.quadrature_point(q);
+            const Point<dim> normal = scratch.fe_face_values.normal_vector(q);
+            const Tensor<1,dim> convection
+            = scratch.convection_velocity.value(quadrature_point);
+
+            const double tau_stab = (5. +
+                                     std::abs(convection * normal));
+
+            for (unsigned int k=0; k<scratch.fe_local_support_on_face[face].size(); ++k)
+              {
+                const unsigned int kk=scratch.fe_local_support_on_face[face][k];
+                scratch.q_phi[k] = scratch.fe_face_values_local[fluxes].value(kk,q);
+                scratch.u_phi[k] = scratch.fe_face_values_local[scalar].value(kk,q);
+              }
+
+            if (!task_data.trace_reconstruct)
+              {
+                for (unsigned int k=0; k<scratch.fe_support_on_face[face].size(); ++k)
+                  scratch.tr_phi[k] =
+                    scratch.fe_face_values.shape_value(scratch.fe_support_on_face[face][k],q);
+                for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
+                  for (unsigned int j=0; j<scratch.fe_support_on_face[face].size(); ++j)
+                    {
+                      const unsigned int ii=scratch.fe_local_support_on_face[face][i];
+                      const unsigned int jj=scratch.fe_support_on_face[face][j];
+                      scratch.lf_matrix(ii,jj) += (
+                                                    (scratch.q_phi[i] * normal
+                                                     +
+                                                     (convection * normal -
+                                                      tau_stab) * scratch.u_phi[i])
+                                                    * scratch.tr_phi[j]
+                                                  ) * JxW;
+
+                      scratch.fl_matrix(jj,ii) -= (
+                                                    (scratch.q_phi[i] * normal
+                                                     +
+                                                     tau_stab * scratch.u_phi[i])
+                                                    * scratch.tr_phi[j]
+                                                  ) * JxW;
+                    }
+
+                for (unsigned int i=0; i<scratch.fe_support_on_face[face].size(); ++i)
+                  for (unsigned int j=0; j<scratch.fe_support_on_face[face].size(); ++j)
+                    {
+                      const unsigned int ii=scratch.fe_support_on_face[face][i];
+                      const unsigned int jj=scratch.fe_support_on_face[face][j];
+                      task_data.cell_matrix(ii,jj) += (
+                                                        (convection * normal - tau_stab) *
+                                                        scratch.tr_phi[i] * scratch.tr_phi[j]
+                                                      ) * JxW;
+                    }
+
+                if (cell->face(face)->at_boundary()
+                    &&
+                    (cell->face(face)->boundary_indicator() == 1))
+                  {
+                    const double neumann_value =
+                      - scratch.exact_solution.gradient (quadrature_point) * normal
+                      + convection * normal * scratch.exact_solution.value(quadrature_point);
+                    for (unsigned int i=0; i<scratch.fe_support_on_face[face].size(); ++i)
+                      {
+                        const unsigned int ii=scratch.fe_support_on_face[face][i];
+                        task_data.cell_vector(ii) += scratch.tr_phi[i] * neumann_value * JxW;
+                      }
+                  }
+              }
+
+            for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
+              for (unsigned int j=0; j<scratch.fe_local_support_on_face[face].size(); ++j)
+                {
+                  const unsigned int ii=scratch.fe_local_support_on_face[face][i];
+                  const unsigned int jj=scratch.fe_local_support_on_face[face][j];
+                  scratch.ll_matrix(ii,jj) += tau_stab * scratch.u_phi[i] * scratch.u_phi[j] * JxW;
+                }
+  
+            if (task_data.trace_reconstruct)
+              for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
+                {
+                  const unsigned int ii=scratch.fe_local_support_on_face[face][i];
+                  scratch.l_rhs(ii) -= (scratch.q_phi[i] * normal
+                                        +
+                                        scratch.u_phi[i] * (convection * normal - tau_stab)
+                                       ) * scratch.trace_values[q] * JxW;
+                }
+          }
+      }
+
+    scratch.ll_matrix.gauss_jordan();
+
+    if (task_data.trace_reconstruct == false)
+      {
+        scratch.fl_matrix.mmult(scratch.tmp_matrix, scratch.ll_matrix);
+        scratch.tmp_matrix.vmult_add(task_data.cell_vector, scratch.l_rhs);
+        scratch.tmp_matrix.mmult(task_data.cell_matrix, scratch.lf_matrix, true);
+        cell->get_dof_indices(task_data.dof_indices);
+      }
+    else
+      {
+        scratch.ll_matrix.vmult(scratch.tmp_rhs, scratch.l_rhs);
+        loc_cell->set_dof_values(scratch.tmp_rhs, solution_local);
+      }
+  }
+
+
+
+  template <int dim>
+  void HDG<dim>::copy_local_to_global(const PerTaskData &data)
+  {
+    if (data.trace_reconstruct == false)
+      constraints.distribute_local_to_global (data.cell_matrix,
+                                              data.cell_vector,
+                                              data.dof_indices,
+                                              system_matrix, system_rhs);
+  }
+
+
+
+  template <int dim>
+  void HDG<dim>::solve ()
+  {
+    SolverControl solver_control (system_matrix.m()*10,
+                                  1e-11*system_rhs.l2_norm());
+    SolverBicgstab<> solver (solver_control, false);
+    solver.solve (system_matrix, solution, system_rhs,
+                  PreconditionIdentity());
+
+    system_matrix.clear();
+
+    constraints.distribute(solution);
+
+    assemble_system(true);
+  }
+
+
+
+
+  template <int dim>
+  void
+  HDG<dim>::postprocess()
+  {
+    {
+      const QGauss<dim>   quadrature_formula(fe_u_post.degree+1);
+      const UpdateFlags local_flags (update_values);
+      const UpdateFlags flags ( update_values | update_gradients |
+                                update_JxW_values);
+
+      PostProcessScratchData scratch (fe_u_post, fe_local,
+                                      quadrature_formula,
+                                      local_flags,
+                                      flags);
+
+      WorkStream::run(dof_handler_u_post.begin_active(),
+                      dof_handler_u_post.end(),
+                      std_cxx1x::bind (&HDG<dim>::postprocess_one_cell,
+                                       std_cxx1x::ref(*this),
+                                       std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3),
+                      std_cxx1x::function<void(const unsigned int &)>(),
+                      scratch,
+                      0U);
+    }
+
+    Vector<float> difference_per_cell (triangulation.n_active_cells());
+
+    ComponentSelectFunction<dim> value_select (dim, dim+1);
+    VectorTools::integrate_difference (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 (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();
+
+    VectorTools::integrate_difference (dof_handler_u_post,
+                                       solution_u_post,
+                                       Solution<dim>(),
+                                       difference_per_cell,
+                                       QGauss<dim>(fe.degree+3),
+                                       VectorTools::L2_norm);
+    const double post_error = difference_per_cell.l2_norm();
+
+    convergence_table.add_value("cells",     triangulation.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);
+    convergence_table.add_value("val L2-post", post_error);
+  }
+
+
+
+  template <int dim>
+  void
+  HDG<dim>::postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+                                  PostProcessScratchData &scratch,
+                                  unsigned int &)
+  {
+    typename DoFHandler<dim>::active_cell_iterator
+    loc_cell (&triangulation,
+              cell->level(),
+              cell->index(),
+              &dof_handler_local);
+
+    scratch.fe_values_local.reinit (loc_cell);
+    scratch.fe_values.reinit(cell);
+
+    FEValuesExtractors::Vector fluxes(0);
+    FEValuesExtractors::Scalar scalar(dim);
+
+    const unsigned int n_q_points = scratch.fe_values.get_quadrature().size();
+    const unsigned int dofs_per_cell = scratch.fe_values.dofs_per_cell;
+
+    scratch.fe_values_local[scalar].get_function_values(solution_local, scratch.u_values);
+    scratch.fe_values_local[fluxes].get_function_values(solution_local, scratch.u_gradients);
+
+    double sum = 0;
+    for (unsigned int i=1; i<dofs_per_cell; ++i)
+      {
+        for (unsigned int j=0; j<dofs_per_cell; ++j)
+          {
+            sum = 0;
+            for (unsigned int q=0; q<n_q_points; ++q)
+              sum += (scratch.fe_values.shape_grad(i,q) *
+                      scratch.fe_values.shape_grad(j,q)
+                     ) * scratch.fe_values.JxW(q);
+            scratch.cell_matrix(i,j) = sum;
+          }
+
+        sum = 0;
+        for (unsigned int q=0; q<n_q_points; ++q)
+          sum -= (scratch.fe_values.shape_grad(i,q) * scratch.u_gradients[q]
+                 ) * scratch.fe_values.JxW(q);
+        scratch.cell_rhs(i) = sum;
+      }
+    for (unsigned int j=0; j<dofs_per_cell; ++j)
+      {
+        sum = 0;
+        for (unsigned int q=0; q<n_q_points; ++q)
+          sum += scratch.fe_values.shape_value(j,q) * scratch.fe_values.JxW(q);
+        scratch.cell_matrix(0,j) = sum;
+      }
+    {
+      sum = 0;
+      for (unsigned int q=0; q<n_q_points; ++q)
+        sum += scratch.u_values[q] * scratch.fe_values.JxW(q);
+      scratch.cell_rhs(0) = sum;
+    }
+
+    scratch.cell_matrix.gauss_jordan();
+    scratch.cell_matrix.vmult(scratch.cell_sol, scratch.cell_rhs);
+    cell->distribute_local_to_global(scratch.cell_sol, solution_u_post);
+  }
+
+
+
+  template <int dim>
+  void HDG<dim>::output_results (const unsigned int cycle)
+  {
+    // not included in test
+  }
+
+
+
+  template <int dim>
+  void HDG<dim>::refine_grid (const unsigned int cycle)
+  {
+    // only global refinement
+    triangulation.clear();
+    GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1);
+    triangulation.refine_global(3-dim+cycle/2);
+
+    typename Triangulation<dim>::cell_iterator
+    cell = triangulation.begin (),
+    endc = triangulation.end();
+    for (; cell!=endc; ++cell)
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+        if (cell->face(face)->at_boundary())
+          for (unsigned int d=0; d<dim; ++d)
+            if ((std::fabs(cell->face(face)->center()(d) - (1)) < 1e-12))
+              cell->face(face)->set_boundary_indicator (1);
+  }
+
+
+
+  template <int dim>
+  void HDG<dim>::run ()
+  {
+    for (unsigned int cycle=0; cycle<14-4*dim; ++cycle)
+      {
+        deallog << "Cycle " << cycle << ':' << std::endl;
+
+        refine_grid (cycle);
+        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(deallog.get_file_stream());
+  }
+
+} // end of namespace Step51
+
+
+
+
+int main ()
+{
+  deallog << std::setprecision(6);
+  logfile << std::setprecision(6);
+
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  {
+    Step51::HDG<1> hdg_problem (1);
+    hdg_problem.run ();
+  }
+
+  {
+    Step51::HDG<1> hdg_problem (3);
+    hdg_problem.run ();
+  }
+
+  {
+    Step51::HDG<2> hdg_problem (1);
+    hdg_problem.run ();
+  }
+
+  {
+    Step51::HDG<2> hdg_problem (2);
+    hdg_problem.run ();
+  }
+
+  return 0;
+}
diff --git a/tests/bits/step-51p.output b/tests/bits/step-51p.output
new file mode 100644 (file)
index 0000000..a0f2912
--- /dev/null
@@ -0,0 +1,165 @@
+
+DEAL::Cycle 0:
+DEAL::   Number of degrees of freedom: 9
+DEAL:Bicgstab::Starting value 8.19443
+DEAL:Bicgstab::Convergence step 8 value 0
+DEAL::Cycle 1:
+DEAL::   Number of degrees of freedom: 13
+DEAL:Bicgstab::Starting value 13.9539
+DEAL:Bicgstab::Convergence step 12 value 0
+DEAL::Cycle 2:
+DEAL::   Number of degrees of freedom: 17
+DEAL:Bicgstab::Starting value 13.0731
+DEAL:Bicgstab::Convergence step 16 value 0
+DEAL::Cycle 3:
+DEAL::   Number of degrees of freedom: 25
+DEAL:Bicgstab::Starting value 12.1492
+DEAL:Bicgstab::Convergence step 24 value 1.07658e-10
+DEAL::Cycle 4:
+DEAL::   Number of degrees of freedom: 33
+DEAL:Bicgstab::Starting value 11.0492
+DEAL:Bicgstab::Convergence step 34 value 0
+DEAL::Cycle 5:
+DEAL::   Number of degrees of freedom: 49
+DEAL:Bicgstab::Starting value 9.35310
+DEAL:Bicgstab::Convergence step 51 value 0
+DEAL::Cycle 6:
+DEAL::   Number of degrees of freedom: 65
+DEAL:Bicgstab::Starting value 8.20475
+DEAL:Bicgstab::Convergence step 75 value 0
+DEAL::Cycle 7:
+DEAL::   Number of degrees of freedom: 97
+DEAL:Bicgstab::Starting value 6.76136
+DEAL:Bicgstab::Convergence step 104 value 0
+DEAL::Cycle 8:
+DEAL::   Number of degrees of freedom: 129
+DEAL:Bicgstab::Starting value 5.87456
+DEAL:Bicgstab::Convergence step 142 value 0
+DEAL::Cycle 9:
+DEAL::   Number of degrees of freedom: 193
+DEAL:Bicgstab::Starting value 4.80772
+DEAL:Bicgstab::Convergence step 211 value 0
+cells dofs     val L2        grad L2      val L2-post   
+8     9    4.080e-01 -    1.722e+00 -    3.362e-01 -    
+12    13   1.659e-01 2.22 8.766e-01 1.67 2.151e-02 6.78 
+16    17   7.996e-02 2.54 4.531e-01 2.29 1.125e-02 2.25 
+24    25   3.705e-02 1.90 2.028e-01 1.98 3.372e-03 2.97 
+32    33   2.125e-02 1.93 1.140e-01 2.00 1.430e-03 2.98 
+48    49   9.610e-03 1.96 5.048e-02 2.01 4.242e-04 3.00 
+64    65   5.449e-03 1.97 2.831e-02 2.01 1.787e-04 3.00 
+96    97   2.440e-03 1.98 1.253e-02 2.01 5.281e-05 3.01 
+128   129  1.377e-03 1.99 7.035e-03 2.01 2.224e-05 3.01 
+192   193  6.141e-04 1.99 3.119e-03 2.01 6.573e-06 3.01 
+DEAL::Cycle 0:
+DEAL::   Number of degrees of freedom: 9
+DEAL:Bicgstab::Starting value 7.96477
+DEAL:Bicgstab::Convergence step 8 value 0
+DEAL::Cycle 1:
+DEAL::   Number of degrees of freedom: 13
+DEAL:Bicgstab::Starting value 13.4076
+DEAL:Bicgstab::Convergence step 12 value 0
+DEAL::Cycle 2:
+DEAL::   Number of degrees of freedom: 17
+DEAL:Bicgstab::Starting value 12.9549
+DEAL:Bicgstab::Convergence step 16 value 0
+DEAL::Cycle 3:
+DEAL::   Number of degrees of freedom: 25
+DEAL:Bicgstab::Starting value 12.1266
+DEAL:Bicgstab::Convergence step 25 value 0
+DEAL::Cycle 4:
+DEAL::   Number of degrees of freedom: 33
+DEAL:Bicgstab::Starting value 11.0426
+DEAL:Bicgstab::Convergence step 33 value 0
+DEAL::Cycle 5:
+DEAL::   Number of degrees of freedom: 49
+DEAL:Bicgstab::Starting value 9.35200
+DEAL:Bicgstab::Convergence step 49 value 0
+DEAL::Cycle 6:
+DEAL::   Number of degrees of freedom: 65
+DEAL:Bicgstab::Starting value 8.20444
+DEAL:Bicgstab::Convergence step 71 value 0
+DEAL::Cycle 7:
+DEAL::   Number of degrees of freedom: 97
+DEAL:Bicgstab::Starting value 6.76131
+DEAL:Bicgstab::Convergence step 111 value 0
+DEAL::Cycle 8:
+DEAL::   Number of degrees of freedom: 129
+DEAL:Bicgstab::Starting value 5.87455
+DEAL:Bicgstab::Convergence step 144 value 0
+DEAL::Cycle 9:
+DEAL::   Number of degrees of freedom: 193
+DEAL:Bicgstab::Starting value 4.80772
+DEAL:Bicgstab::Convergence step 206 value 0
+cells dofs     val L2        grad L2      val L2-post   
+8     9    2.141e-02 -    1.398e-01 -    3.487e-03 -    
+12    13   6.701e-03 2.87 3.811e-02 3.21 4.976e-04 4.80 
+16    17   1.693e-03 4.78 9.999e-03 4.65 1.085e-04 5.29 
+24    25   3.501e-04 3.89 2.009e-03 3.96 1.440e-05 4.98 
+32    33   1.126e-04 3.94 6.373e-04 3.99 3.412e-06 5.01 
+48    49   2.256e-05 3.97 1.259e-04 4.00 4.472e-07 5.01 
+64    65   7.179e-06 3.98 3.980e-05 4.00 1.057e-07 5.01 
+96    97   1.425e-06 3.99 7.848e-06 4.00 1.388e-08 5.01 
+128   129  4.519e-07 3.99 2.481e-06 4.00 3.400e-09 4.89 
+192   193  8.946e-08 3.99 4.894e-07 4.00 1.008e-09 3.00 
+DEAL::Cycle 0:
+DEAL::   Number of degrees of freedom: 80
+DEAL:Bicgstab::Starting value 17.6593
+DEAL:Bicgstab::Convergence step 31 value 1.30726e-10
+DEAL::Cycle 1:
+DEAL::   Number of degrees of freedom: 168
+DEAL:Bicgstab::Starting value 11.4720
+DEAL:Bicgstab::Convergence step 44 value 0
+DEAL::Cycle 2:
+DEAL::   Number of degrees of freedom: 288
+DEAL:Bicgstab::Starting value 16.7917
+DEAL:Bicgstab::Convergence step 63 value 1.47629e-10
+DEAL::Cycle 3:
+DEAL::   Number of degrees of freedom: 624
+DEAL:Bicgstab::Starting value 12.5049
+DEAL:Bicgstab::Convergence step 90 value 0
+DEAL::Cycle 4:
+DEAL::   Number of degrees of freedom: 1088
+DEAL:Bicgstab::Starting value 9.93628
+DEAL:Bicgstab::Convergence step 120 value 0
+DEAL::Cycle 5:
+DEAL::   Number of degrees of freedom: 2400
+DEAL:Bicgstab::Starting value 6.94024
+DEAL:Bicgstab::Convergence step 177 value 0
+cells dofs     val L2        grad L2      val L2-post   
+16    80   1.043e+01 -    1.616e+01 -    1.034e+01 -    
+36    168  3.699e+00 2.56 8.081e+00 1.71 3.637e+00 2.58 
+64    288  8.022e-01 5.31 4.671e+00 1.91 4.622e-01 7.17 
+144   624  2.947e-01 2.47 2.255e+00 1.80 6.403e-02 4.88 
+256   1088 1.638e-01 2.04 1.343e+00 1.80 2.937e-02 2.71 
+576   2400 7.616e-02 1.89 6.522e-01 1.78 9.203e-03 2.86 
+DEAL::Cycle 0:
+DEAL::   Number of degrees of freedom: 120
+DEAL:Bicgstab::Starting value 19.9021
+DEAL:Bicgstab::Convergence step 39 value 1.15234e-10
+DEAL::Cycle 1:
+DEAL::   Number of degrees of freedom: 252
+DEAL:Bicgstab::Starting value 15.8809
+DEAL:Bicgstab::Convergence step 56 value 0
+DEAL::Cycle 2:
+DEAL::   Number of degrees of freedom: 432
+DEAL:Bicgstab::Starting value 16.3149
+DEAL:Bicgstab::Convergence step 71 value 0
+DEAL::Cycle 3:
+DEAL::   Number of degrees of freedom: 936
+DEAL:Bicgstab::Starting value 12.6721
+DEAL:Bicgstab::Convergence step 109 value 1.21847e-10
+DEAL::Cycle 4:
+DEAL::   Number of degrees of freedom: 1632
+DEAL:Bicgstab::Starting value 10.0470
+DEAL:Bicgstab::Convergence step 144 value 0
+DEAL::Cycle 5:
+DEAL::   Number of degrees of freedom: 3600
+DEAL:Bicgstab::Starting value 6.97340
+DEAL:Bicgstab::Convergence step 198 value 0
+cells dofs     val L2        grad L2      val L2-post   
+16    120  8.296e-01 -    6.543e+00 -    4.572e-01 -    
+36    252  5.045e-01 1.23 3.062e+00 1.87 3.570e-01 0.61 
+64    432  1.736e-01 3.71 1.503e+00 2.47 6.468e-02 5.94 
+144   936  5.584e-02 2.80 5.042e-01 2.69 1.006e-02 4.59 
+256   1632 2.731e-02 2.49 2.274e-01 2.77 3.120e-03 4.07 
+576   3600 8.582e-03 2.85 7.254e-02 2.82 6.213e-04 3.98 
index 7d00b8382677a7105286f6441b4b817bb861081a..a9e16506a10e6917395fad236f855145929069ba 100644 (file)
 template <int dim>
 void test_2d_3d (std::vector<FiniteElement<dim> *> &fe_datas)
 {
-                                  // Face Q elements
-  fe_datas.push_back(new FE_FaceQ<dim> (0));
-  deallog << (*fe_datas.rbegin())->get_name() << std::endl;
-  fe_datas.push_back(new FE_FaceQ<dim> (1));
-  deallog << (*fe_datas.rbegin())->get_name() << std::endl;
-  fe_datas.push_back(new FE_FaceQ<dim> (3));
-  deallog << (*fe_datas.rbegin())->get_name() << std::endl;
-                                  // Face P elements
-  fe_datas.push_back(new FE_FaceP<dim> (0));
-  deallog << (*fe_datas.rbegin())->get_name() << std::endl;
-  fe_datas.push_back(new FE_FaceP<dim> (1));
-  deallog << (*fe_datas.rbegin())->get_name() << std::endl;
-  fe_datas.push_back(new FE_FaceP<dim> (3));
-  deallog << (*fe_datas.rbegin())->get_name() << std::endl;
                                   // Vector DG elements
   fe_datas.push_back(
     new FE_DGRaviartThomas<dim>(0));
@@ -136,6 +122,21 @@ void test_fe_datas()
                                        FE_Q<dim> (2), 1));
   deallog << (*fe_datas.rbegin())->get_name() << std::endl;
 
+                                  // Face Q elements
+  fe_datas.push_back(new FE_FaceQ<dim> (0));
+  deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+  fe_datas.push_back(new FE_FaceQ<dim> (1));
+  deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+  fe_datas.push_back(new FE_FaceQ<dim> (3));
+  deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+                                  // Face P elements
+  fe_datas.push_back(new FE_FaceP<dim> (0));
+  deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+  fe_datas.push_back(new FE_FaceP<dim> (1));
+  deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+  fe_datas.push_back(new FE_FaceP<dim> (3));
+  deallog << (*fe_datas.rbegin())->get_name() << std::endl;
+
 
   // Check vector elements in 2d and higher only
   test_2d_3d (fe_datas);

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.