The step-12 tutorial program

Table of contents
  1. Introduction
  2. The commented program
  1. Results
  2. The plain program

Introduction

Overview

This example is devoted to the discontinuous Galerkin method, or in short: DG method. It includes the following topics.

  1. Discretization of the linear transport equation with the DG method
  2. Two different assembling routines for the system matrix based on face terms given as a sum of integrals that
    1. loops over all cell and all their faces, or that
    2. loops over all faces, whereas each face is treated only once.
  3. Time comparison of the two assembling routines.

Problem

The DG method was first introduced to discretize simple transport equations. Over the past years DG methods have been applied to a variety of problems and many different schemes were introduced employing a big zoo of different convective and diffusive fluxes. As this example's purpose is to illustrate some implementational issues of the DG discretization only, here we simply consider the linear transport equation

\[ \nabla\cdot \left({\mathbf \beta} u\right)=f \qquad\mbox{in }\Omega, \qquad\qquad\qquad\mathrm{[transport-equation]}\]

subject to the boundary conditions

\[ u=g\quad\mbox{on }\Gamma_-, \]

on the inflow part $\Gamma_-$ of the boundary $\Gamma=\partial\Omega$ of the domain. Here, ${\mathbf \beta}={\mathbf \beta}({\bf x})$ denotes a vector field, $f$ a source function, $u$ the (scalar) solution function, $g$ a boundary value function,

\[ \Gamma_-:=\{{\bf x}\in\Gamma, {\mathbf \beta}({\bf x})\cdot{\bf n}({\bf x})<0\} \]

the inflow part of the boundary of the domain and ${\bf n}$ denotes the unit outward normal to the boundary $\Gamma$. Equation [transport-equation] is the conservative version of the transport equation already considered in step 9 of this tutorial.

In particular, we consider problem [transport-equation] on $\Omega=[0,1]^2$ with ${\mathbf \beta}=\frac{1}{|x|}(-x_2, x_1)$ representing a circular counterclockwise flow field, $f=0$ and $g=1$ on ${\bf x}\in\Gamma_-^1:=[0,0.5]\times\{0\}$ and $g=0$ on ${\bf x}\in \Gamma_-\setminus \Gamma_-^1$.

Discretization

Following the general paradigm of deriving DG discretizations for purely hyperbolic equations, we first consider the general hyperbolic problem

\[ \nabla\cdot {\mathcal F}(u)=f \qquad\mbox{in }\Omega, \]

subject to appropriate boundary conditions. Here ${\mathcal F}$ denotes the flux function of the equation under consideration that in our case, see equation [transport-equation], is represented by ${\mathcal F}(u)={\mathbf \beta} u$. For deriving the DG discretization we start with a variational, mesh-dependent formulation of the problem,

\[ \sum_\kappa\left\{-({\mathcal F}(u),\nabla v)_\kappa+({\mathcal F}(u)\cdot{\bf n}, v)_{\partial\kappa}\right\}=(f,v)_\Omega, \]

that originates from [transport-equation] by multiplication with a test function $v$ and integration by parts on each cell $\kappa$ of the triangulation. Here $(\cdot, \cdot)_\kappa$ and $(\cdot, \cdot)_{\partial\kappa}$ simply denote the integrals over the cell $\kappa$ and the boundary $\partial\kappa$ of the cell, respectively. To discretize the problem, the functions $u$ and $v$ are replaced by discrete functions $u_h$ and $v_h$ that in the case of discontinuous Galerkin methods belong to the space $V_h$ of discontinuous piecewise polynomial functions of some degree $p$. Due to the discontinuity of the discrete function $u_h$ on interelement faces, the flux ${\mathcal F}(u)\cdot{\bf n}$ must be replaced by a numerical flux function ${\mathcal H}(u_h^+, u_h^-, {\bf n})$, where $u_h^+|_{\partial\kappa}$ denotes the inner trace (w.r.t. the cell $\kappa$) of $u_h$ and $u_h^-|_{\partial\kappa}$ the outer trace, i.e. the value of $u_h$ on the neighboring cell. Furthermore the numerical flux function ${\mathcal H}$, among other things, must be consistent, i.e.

\[ {\mathcal H}(u,u,{\bf n})={\mathcal F}(u)\cdot{\bf n}, \]

and conservative, i.e.

\[ {\mathcal H}(v,w,{\bf n})=-{\mathcal H}(w,v,-{\bf n}). \qquad\qquad\qquad\mathrm{[conservativity]}\]

This yields the following discontinuous Galerkin discretization: find $u_h\in V_h$ such that

\[ \sum_\kappa\left\{-({\mathcal F}(u_h),\nabla v_h)_\kappa+({\mathcal H}(u_h^+,u_h^-,{\bf n}), v_h)_{\partial\kappa}\right\}=(f,v_h)_\Omega, \quad\forall v_h\in V_h. \qquad\qquad\qquad\mathrm{[dg-general1]}\]

Boundary conditions are realized by replacing $u_h^-$ on the inflow boundary $\Gamma_-$ by the boundary function $g$. In the special case of the [transport-equation] the numerical flux in its simplest form is given by

\[ {\mathcal H}(u_h^+,u_h^-,{\bf n})({\bf x})=\left\{\begin{array}{ll} ({\mathbf \beta}\cdot{\bf n}\, u_h^-)({\bf x}),&\mbox{for } {\mathbf \beta}({\bf x})\cdot{\bf n}({\bf x})<0,\\ ({\mathbf \beta}\cdot{\bf n}\, u_h^+)({\bf x}),&\mbox{for } {\mathbf \beta}({\bf x})\cdot{\bf n}({\bf x})\geq 0, \end{array} \right. \qquad\qquad\qquad\mathrm{[upwind-flux]}\]

where on the inflow part of the cell the value is taken from the neighboring cell, $u_h^-$, and on the outflow part the value is taken from the current cell, $u_h^+$. Hence, the discontinuous Galerkin scheme for the [transport-equation] is given by: find $u_h\in V_h$ such that for all $v_h\in V_h$ following equation holds:

\[ \sum_\kappa\left\{-(u_h,{\mathbf \beta}\cdot\nabla v_h)_\kappa +({\mathbf \beta}\cdot{\bf n}\, u_h, v_h)_{\partial\kappa_+} +({\mathbf \beta}\cdot{\bf n}\, u_h^-, v_h)_{\partial\kappa_-\setminus\Gamma}\right\} =(f,v_h)_\Omega-({\mathbf \beta}\cdot{\bf n}\, g, v_h)_{\Gamma_-}, \qquad\qquad\qquad\mathrm{[dg-transport1]}\]

where $\partial\kappa_-:=\{{\bf x}\in\partial\kappa, {\mathbf \beta}({\bf x})\cdot{\bf n}({\bf x})<0\}$ denotes the inflow boundary and $\partial\kappa_+=\partial\kappa\setminus \partial \kappa_-$ the outflow part of cell $\kappa$. Below, this equation will be referred to as first version of the DG method. We note that after a second integration by parts, we obtain: find $u_h\in V_h$ such that

\[ \sum_\kappa\left\{(\nabla\cdot\{{\mathbf \beta} u_h\},v_h)_\kappa -({\mathbf \beta}\cdot{\bf n} [u_h], v_h)_{\partial\kappa_-}\right\} =(f,v_h)_\Omega, \quad\forall v_h\in V_h, \]

where $[u_h]=u_h^+-u_h^-$ denotes the jump of the discrete function between two neighboring cells and is defined to be $[u_h]=u_h^+-g$ on the boundary of the domain. This is the discontinuous Galerkin scheme for the transport equation given in its original notation. Nevertheless, we will base the implementation of the scheme on the form given by [dg-general1] and [upwind-flux], or [dg-transport1], respectively.

Finally, we rewrite [dg-general1] in terms of a summation over all faces where each face $e=\partial \kappa\cap\partial \kappa'$ between two neighboring cells $\kappa$ and $\kappa'$ occurs twice: Find $u_h\in V_h$ such that

\[ -\sum_\kappa({\mathcal F}(u_h),\nabla v_h)_\kappa+\sum_e\left\{({\mathcal H}(u_h^+,u_h^-,{\bf n}), v_h)_e+({\mathcal H}(u_h^-, u_h^+,-{\bf n}), v_h^-)_{e\setminus\Gamma}\right\}=(f,v_h)_\Omega \quad\forall v_h\in V_h, \qquad\qquad\qquad\mathrm{[dg-general2]}\]

By employing [conservativity] of the numerical flux this equation simplifies to: find $u_h\in V_h$ such that

\[ -\sum_\kappa({\mathcal F}(u_h),\nabla v_h)_\kappa+\sum_e({\mathcal H}(u_h^+,u_h^-,{\bf n}), [v_h])_{e\setminus\Gamma}+({\mathcal H}(u_h,g,{\bf n}), v_h)_{\Gamma}=(f,v_h)_\Omega \quad\forall v_h\in V_h. \qquad\qquad\qquad\mathrm{[dg-general3]}\]

Whereas the outer unit normal ${\bf n}|_{\partial\kappa}$ is uniquely defined this is not so for ${\bf n}_e$ as the latter might be the normal from either side of the face. Hence, we need to fix the normal ${\bf n}$ on the face to be one of the two normals and denote the other normal by $-{\bf n}$. This way we get $-{\bf n}$ in the second face term in [dg-general2] that finally produces the minus sign in the jump $[v_h]$ in equation [dg-general3].

For the linear [transport-equation] equation [dg-general3] simplifies to

\[ -\sum_\kappa(u_h,{\mathbf \beta}\cdot\nabla v_h)_\kappa+\sum_e\left\{({\mathbf \beta}\cdot{\bf n}\, u_h, [v_h])_{e_+\setminus\Gamma}+({\mathbf \beta}\cdot{\bf n}\, u_h^-, [v_h])_{e_-\setminus\Gamma}\right\}=(f,v_h)_\Omega-({\mathbf \beta}\cdot{\bf n}\, g, v_h)_{\Gamma_-}, \qquad\qquad\qquad\mathrm{[dg-transport2]}\]

which will be refered to as second version of the DG method.

Implementation

As already mentioned at the beginning of this example we will implement assembling the system matrix in two different ways. The first one will be based on the first version [dg-transport1] of the DG method that includes a sum of integrals over all cell boundaries $\partial\kappa$. This is realized by a loop over all cells and a nested loop over all faces of each cell. Thereby each inner face $e=\partial\kappa\cap\partial \kappa'$ is treated twice, the first time when the outer loop treats cell $\kappa$ and the second time when it treats cell $\kappa'$. This way some values like the shape function values at quadrature points on faces need to be computed twice.

To overcome this overhead and for comparison, we implement assembling of matrix also in a second and different way. This will be based on the second version [dg-transport2] that includes a sum of integrals over all faces $e$. Here, several difficulties occurs.

  1. As degrees of freedom are associated with cells (and not to faces) and as a normal is only defined w.r.t. a cell adjacent to the face we cannot simply run over all faces of the triangulation but need to perform the nested loop over all cells and all faces of each cell like in the first implementation. This, because in deal.II faces are accessible from cells but not visa versa.
  2. Due to the nested loop we arrive twice at each face. In order to assemble face terms only once we either need to track which faces we have treated before, or we introduce a simple rule that decides which of the two adjacent cells the face should be accessed and treated from. Here, we employ the second approach and define the following rule:
    1. If the two cells adjacent to a face are of the same refinement level we access and treat the face from the cell with lower index on this level.
    2. If the two cells are of different refinement levels we access and treat the face from the coarser cell.

Before we start with the description of the code we first introduce its main ingredients. The main class is called DGMethod. It comprises all basic objects like the triangulation, the dofhandler, the system matrix and solution vectors. Furthermore it has got some member functions, the most prominent of which are the assemble_system1 and assemble_system2 functions that implement the two different ways mentioned above for assembling the system matrix. Within these assembling routines several different cases must be distinguished while performing the nested loops over all cells and all faces of each cell and assembling the respective face terms. While sitting on the current cell and looking at a specific face there are the cases

  1. face is at boundary,
  2. neighboring cell is finer,
  3. neighboring cell is of the same refinement level, and
  4. neighboring cell is coarser

where the `neighboring cell' and the current cell have the mentioned faces in common. In last three cases the assembling of the face terms are almost the same. Hence, we can implement the assembling of the face terms either by `copy and paste' (the lazy way, whose disadvantages come up when the scheme or the equation might want to be changed afterwards) or by calling a separate function that covers all three cases. To be kind of educational within this tutorial we perform the latter approach, of course. We go even further and encapsulate this function and everything that is needed for assembling the specific equation under consideration within a class called DGTransportEquation. This class includes objects of all equation--specific functions, the RHS and the BoundaryValues class, both derived from the Function class, and the Beta class representing the vector field. Furthermore, the DGTransportEquation class comprises member functions assemble_face_terms1 and assemble_face_terms2 that are invoked by the assemble_system1 and assemble_system2 functions of the DGMethod, respectively, and the functions assemble_cell_term and assemble_boundary_term that are the same for both assembling routines. Due to the encapsulation of all equation- and scheme-specific functions, the DGTransportEquation class can easily be replaced by a similar class that implements a different equation and a different DG method. Indeed, the implementation of the assemble_system1 and assemble_system2 functions of the DGMethod class will be general enough to serve for different DG methods, different equations, even for systems of equations (!) and, under small modifications, for nonlinear problems. Finally, we note that the program is dimension independent, i.e. after replacing DGMethod<2> by DGMethod<3> the code runs in 3d.

The commented program

The first few files have already been covered in previous examples and will thus not be further commented on.

 #include <base/quadrature_lib.h>
 #include <base/function.h>
 #include <lac/vector.h>
 #include <lac/sparse_matrix.h>
 #include <grid/tria.h>
 #include <grid/grid_generator.h>
 #include <grid/grid_out.h>
 #include <grid/grid_refinement.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
 #include <fe/fe_values.h>
 #include <dofs/dof_handler.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
 #include <numerics/data_out.h>

This is the first new file. It declares the MappingQ1 class that gives the standard bilinear mapping. For bilinear mappings use an object of this class rather than an object of the MappingQ(1) class, as the MappingQ1 class is optimized due to the pre-knowledge of the actual polynomial degree 1.

 #include <fe/mapping_q1.h>

Here the discontinuous finite elements are defined. They are used in the same way as all other finite elements, though -- as you have seen in previous tutorial programs -- there isn't much user interaction with finite element classes at all: the are passed to DoFHandler and FEValues objects, and that is about it.

 #include <fe/fe_dgq.h>

We are going to use the simplest possible solver, called Richardson iteration, that represents a simple defect correction. This, in combination with a block SSOR preconditioner (defined in precondition_block.h), that uses the special block matrix structure of system matrices arising from DG discretizations.

 #include <lac/solver_richardson.h>
 #include <lac/precondition_block.h>

We are going to use gradients as refinement indicator.

 #include <numerics/derivative_approximation.h>

Finally we do some time comparison using the Timer class.

 #include <base/timer.h>

And this again is C++:

 #include <iostream>
 #include <fstream>

The last step is as in all previous programs:

 using namespace dealii;

Equation data

First we define the classes representing the equation-specific functions. Both classes, RHS and BoundaryValues, are derived from the Function class. Only the value_list function are implemented because only lists of function values are computed rather than single values.

 template <int dim>
 class RHS:  public Function<dim>
 {
   public:
     virtual void value_list (const std::vector<Point<dim> > &points,
                              std::vector<double> &values,
                              const unsigned int component=0) const;
 };
 
 
 template <int dim>
 class BoundaryValues:  public Function<dim>
 {
   public:
     virtual void value_list (const std::vector<Point<dim> > &points,
                              std::vector<double> &values,
                              const unsigned int component=0) const;
 };

The class Beta represents the vector valued flow field of the linear transport equation and is not derived from the Function class as we prefer to get function values of type Point rather than of type Vector<double>. This, because there exist scalar products between Point and Point as well as between Point and Tensor, simplifying terms like $\beta\cdot n$ and $\beta\cdot\nabla v$.

An unnecessary empty constructor is added to the class to work around a bug in Compaq's cxx compiler which otherwise reports an error about an omitted initializer for an object of this class further down.

 template <int dim>
 class Beta
 {
   public:
     Beta () {}
     void value_list (const std::vector<Point<dim> > &points,
                      std::vector<Point<dim> > &values) const;
 };

The implementation of the value_list functions of these classes are rather simple. For simplicity the right hand side is set to be zero but will be assembled anyway.

 template <int dim>
 void RHS<dim>::value_list(const std::vector<Point<dim> > &points,
                           std::vector<double> &values,
                           const unsigned int) const
 {

Usually we check whether input parameter have the right sizes.

   Assert(values.size()==points.size(),
          ExcDimensionMismatch(values.size(),points.size()));
 
   for (unsigned int i=0; i<values.size(); ++i)
     values[i]=0;
 }

The flow field is chosen to be circular, counterclockwise, and with the origin as midpoint.

 template <int dim>
 void Beta<dim>::value_list(const std::vector<Point<dim> > &points,
                            std::vector<Point<dim> > &values) const
 {
   Assert(values.size()==points.size(),
          ExcDimensionMismatch(values.size(),points.size()));
 
   for (unsigned int i=0; i<points.size(); ++i)
     {
       values[i](0) = -points[i](1);
       values[i](1) = points[i](0);
       values[i] /= std::sqrt(values[i].square());
     }
 }

Hence the inflow boundary of the unit square [0,1]^2 are the right and the lower boundaries. We prescribe discontinuous boundary values 1 and 0 on the x-axis and value 0 on the right boundary. The values of this function on the outflow boundaries will not be used within the DG scheme.

 template <int dim>
 void BoundaryValues<dim>::value_list(const std::vector<Point<dim> > &points,
                                        std::vector<double> &values,
                                        const unsigned int) const
 {
   Assert(values.size()==points.size(),
          ExcDimensionMismatch(values.size(),points.size()));
 
   for (unsigned int i=0; i<values.size(); ++i)
     {
       if (points[i](0)<0.5)
         values[i]=1.;
       else
         values[i]=0.;
     }
 }

Class: DGTransportEquation

Next we define the equation-dependent and DG-method-dependent class DGTransportEquation. Its member functions were already mentioned in the Introduction and will be explained below. Furthermore it includes objects of the previously defined Beta, RHS and BoundaryValues function classes.

 template <int dim>
 class DGTransportEquation
 {
   public:
     DGTransportEquation();
 
     void assemble_cell_term(const FEValues<dim>& fe_v,
                             FullMatrix<double> &ui_vi_matrix,
                             Vector<double> &cell_vector) const;
     
     void assemble_boundary_term(const FEFaceValues<dim>& fe_v,
                                 FullMatrix<double> &ui_vi_matrix,
                                 Vector<double> &cell_vector) const;
     
     void assemble_face_term1(const FEFaceValuesBase<dim>& fe_v,
                              const FEFaceValuesBase<dim>& fe_v_neighbor,
                              FullMatrix<double> &ui_vi_matrix,
                              FullMatrix<double> &ue_vi_matrix) const;
 
     void assemble_face_term2(const FEFaceValuesBase<dim>& fe_v,
                              const FEFaceValuesBase<dim>& fe_v_neighbor,
                              FullMatrix<double> &ui_vi_matrix,
                              FullMatrix<double> &ue_vi_matrix,
                              FullMatrix<double> &ui_ve_matrix,
                              FullMatrix<double> &ue_ve_matrix) const;
   private:
     const Beta<dim> beta_function;
     const RHS<dim> rhs_function;
     const BoundaryValues<dim> boundary_function;
 };
 
 
 template <int dim>
 DGTransportEquation<dim>::DGTransportEquation ()
                 :
                 beta_function (),
                 rhs_function (),
                 boundary_function ()
 {}

Function: assemble_cell_term

The assemble_cell_term function assembles the cell terms of the discretization. ui_vi_matrix is a cell matrix, i.e. for a DG method of degree 1, it is of size 4 times 4, and cell_vector is of size 4. When this function is invoked, fe_v is already reinit'ed with the current cell before and includes all shape values needed.

 template <int dim>
 void DGTransportEquation<dim>::assemble_cell_term(
   const FEValues<dim> &fe_v,
   FullMatrix<double> &ui_vi_matrix,
   Vector<double> &cell_vector) const
 {

First we ask fe_v for the quadrature weights,

   const std::vector<double> &JxW = fe_v.get_JxW_values ();

Then the flow field beta and the rhs_function are evaluated at the quadrature points,

   std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
   std::vector<double> rhs (fe_v.n_quadrature_points);
   
   beta_function.value_list (fe_v.get_quadrature_points(), beta);
   rhs_function.value_list (fe_v.get_quadrature_points(), rhs);

and the cell matrix and cell vector are assembled due to the terms $-(u,\beta\cdot\nabla v)_K$ and $(f,v)_K$.

   for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
     for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) 
       {
         for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
           ui_vi_matrix(i,j) -= beta[point]*fe_v.shape_grad(i,point)*
                               fe_v.shape_value(j,point) *
                               JxW[point];
         
         cell_vector(i) += rhs[point] * fe_v.shape_value(i,point) * JxW[point];
       }
 }

Function: assemble_boundary_term

The assemble_boundary_term function assembles the face terms at boundary faces. When this function is invoked, fe_v is already reinit'ed with the current cell and current face. Hence it provides the shape values on that boundary face.

 template <int dim>
 void DGTransportEquation<dim>::assemble_boundary_term(
   const FEFaceValues<dim>& fe_v,    
   FullMatrix<double> &ui_vi_matrix,
   Vector<double> &cell_vector) const
 {

Again, as in the previous function, we ask the FEValues object for the quadrature weights

   const std::vector<double> &JxW = fe_v.get_JxW_values ();

but here also for the normals.

   const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();

We evaluate the flow field and the boundary values at the quadrature points.

   std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
   std::vector<double> g(fe_v.n_quadrature_points);
   
   beta_function.value_list (fe_v.get_quadrature_points(), beta);
   boundary_function.value_list (fe_v.get_quadrature_points(), g);

Then we assemble cell vector and cell matrix according to the DG method given in the introduction.

   for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
     {
       const double beta_n=beta[point] * normals[point];      

We assemble the term $(\beta\cdot n u,v)_{\partial\kappa_+}$,

       if (beta_n>0)
         for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
           for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
             ui_vi_matrix(i,j) += beta_n *
                                fe_v.shape_value(j,point) *
                                fe_v.shape_value(i,point) *
                                JxW[point];
       else

and the term $(\beta\cdot n g,v)_{\partial \kappa_-\cap\partial\Omega}$,

         for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
           cell_vector(i) -= beta_n *
                             g[point] *
                             fe_v.shape_value(i,point) *
                             JxW[point];
     }
 }

Function: assemble_face_term1

The assemble_face_term1 function assembles the face terms corresponding to the first version of the DG method, cf. above. For that case, the face terms are given as a sum of integrals over all cell boundaries.

When this function is invoked, fe_v and fe_v_neighbor are already reinit'ed with the current cell and the neighoring cell, respectively, as well as with the current face. Hence they provide the inner and outer shape values on the face.

In addition to the cell matrix ui_vi_matrix this function gets a new argument ue_vi_matrix, that stores contributions to the system matrix that are based on exterior values of $u$ and interior values of $v$. Here we note that ue is the short notation for u exterior and represents $u_h^-$, see the introduction.

 template <int dim>
 void DGTransportEquation<dim>::assemble_face_term1(
   const FEFaceValuesBase<dim>& fe_v,
   const FEFaceValuesBase<dim>& fe_v_neighbor,      
   FullMatrix<double> &ui_vi_matrix,
   FullMatrix<double> &ue_vi_matrix) const
 {

Again, as in the previous function, we ask the FEValues objects for the quadrature weights and the normals

   const std::vector<double> &JxW = fe_v.get_JxW_values ();
   const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();

and we evaluate the flow field at the quadrature points.

   std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
   beta_function.value_list (fe_v.get_quadrature_points(), beta);

Then we assemble the cell matrices according to the DG method given in the introduction.

   for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
     {
       const double beta_n=beta[point] * normals[point];

We assemble the term $(\beta\cdot n u,v)_{\partial\kappa_+}$,

       if (beta_n>0)
         for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
           for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
             ui_vi_matrix(i,j) += beta_n *
                                fe_v.shape_value(j,point) *
                                fe_v.shape_value(i,point) *
                                JxW[point];
       else

and the term $(\beta\cdot n \hat u,v)_{\partial \kappa_-}$.

         for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
           for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
             ue_vi_matrix(i,k) += beta_n *
                                 fe_v_neighbor.shape_value(k,point) *
                                 fe_v.shape_value(i,point) *
                                 JxW[point];
     }
 }

Function: assemble_face_term2

Now we look at the assemble_face_term2 function that assembles the face terms corresponding to the second version of the DG method, cf. above. For that case the face terms are given as a sum of integrals over all faces. Here we need two additional cell matrices ui_ve_matrix and ue_ve_matrix that will store contributions due to terms involving ui and ve as well as ue and ve.

 template <int dim>
 void DGTransportEquation<dim>::assemble_face_term2(
   const FEFaceValuesBase<dim>& fe_v,
   const FEFaceValuesBase<dim>& fe_v_neighbor,      
   FullMatrix<double> &ui_vi_matrix,
   FullMatrix<double> &ue_vi_matrix,
   FullMatrix<double> &ui_ve_matrix,
   FullMatrix<double> &ue_ve_matrix) const
 {

the first few lines are the same

   const std::vector<double> &JxW = fe_v.get_JxW_values ();
   const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
 
   std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
   
   beta_function.value_list (fe_v.get_quadrature_points(), beta);
 
   for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
     {
       const double beta_n=beta[point] * normals[point];
       if (beta_n>0)
         {

This term we've already seen.

           for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
             for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
               ui_vi_matrix(i,j) += beta_n *
                                  fe_v.shape_value(j,point) *
                                  fe_v.shape_value(i,point) *
                                  JxW[point];

We additionally assemble the term $(\beta\cdot n u,\hat v)_{\partial \kappa_+}$,

           for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
             for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
               ui_ve_matrix(k,j) -= beta_n *
                                   fe_v.shape_value(j,point) *
                                   fe_v_neighbor.shape_value(k,point) *
                                   JxW[point];
         }
       else
         {

This one we've already seen, too.

           for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
             for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l)
               ue_vi_matrix(i,l) += beta_n *
                                   fe_v_neighbor.shape_value(l,point) *
                                   fe_v.shape_value(i,point) *
                                   JxW[point];

And this is another new one: $(\beta\cdot n \hat u,\hat v)_{\partial \kappa_-}$.

           for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
             for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l)
               ue_ve_matrix(k,l) -= beta_n *
                                    fe_v_neighbor.shape_value(l,point) *
                                    fe_v_neighbor.shape_value(k,point) *
                                    JxW[point];
         }
     }
 }

Class: DGMethod

After these preparations, we proceed with the main part of this program. The main class, here called DGMethod is basically the main class of step-6. One of the differences is that there's no ConstraintMatrix object. This is, because there are no hanging node constraints in DG discretizations.

 template <int dim>
 class DGMethod
 {
   public:
     DGMethod ();
     ~DGMethod ();
 
     void run ();
     
   private:
     void setup_system ();
     void assemble_system1 ();
     void assemble_system2 ();
     void solve (Vector<double> &solution);
     void refine_grid ();
     void output_results (const unsigned int cycle) const;
     
     Triangulation<dim>   triangulation;
     const MappingQ1<dim> mapping;

Furthermore we want to use DG elements of degree 1 (but this is only specified in the constructor). If you want to use a DG method of a different degree the whole program stays the same, only replace 1 in the constructor by the wanted degree.

     FE_DGQ<dim>          fe;
     DoFHandler<dim>      dof_handler;
 
     SparsityPattern      sparsity_pattern;
     SparseMatrix<double> system_matrix;

We define the quadrature formulae for the cell and the face terms of the discretization.

     const QGauss<dim>   quadrature;
     const QGauss<dim-1> face_quadrature;

And there are two solution vectors, that store the solutions to the problems corresponding to the two different assembling routines assemble_system1 and assemble_system2;

     Vector<double>       solution1;
     Vector<double>       solution2;
     Vector<double>       right_hand_side;

Finally this class includes an object of the DGTransportEquations class described above.

     const DGTransportEquation<dim> dg;
 };
 
 
 template <int dim>
 DGMethod<dim>::DGMethod ()
                 :
                 mapping (),

Change here for DG methods of different degrees.

                 fe (1),
                 dof_handler (triangulation),
                 quadrature (4),
                 face_quadrature (4),
                 dg ()
 {}
 
 
 template <int dim>
 DGMethod<dim>::~DGMethod () 
 {
   dof_handler.clear ();
 }
 
 
 template <int dim>
 void DGMethod<dim>::setup_system ()
 {

First we need to distribute the DoFs.

   dof_handler.distribute_dofs (fe);

The DoFs of a cell are coupled with all DoFs of all neighboring cells. Therefore the maximum number of matrix entries per row is needed when all neighbors of a cell are once more refined than the cell under consideration.

   sparsity_pattern.reinit (dof_handler.n_dofs(),
                            dof_handler.n_dofs(),
                            (GeometryInfo<dim>::faces_per_cell
                             *GeometryInfo<dim>::max_children_per_face+1)*fe.dofs_per_cell);

For DG discretizations we call the function analogue to DoFTools::make_sparsity_pattern.

   DoFTools::make_flux_sparsity_pattern (dof_handler, sparsity_pattern);

All following function calls are already known.

   sparsity_pattern.compress();
   
   system_matrix.reinit (sparsity_pattern);
 
   solution1.reinit (dof_handler.n_dofs());
   solution2.reinit (dof_handler.n_dofs());
   right_hand_side.reinit (dof_handler.n_dofs());
 }

Function: assemble_system1

We proceed with the assemble_system1 function that implements the DG discretization in its first version. This function repeatedly calls the assemble_cell_term, assemble_boundary_term and assemble_face_term1 functions of the DGTransportEquation object. The assemble_boundary_term covers the first case mentioned in the introduction.

1. face is at boundary

This function takes a FEFaceValues object as argument. In contrast to that assemble_face_term1 takes two FEFaceValuesBase objects; one for the shape functions on the current cell and the other for shape functions on the neighboring cell under consideration. Both objects are either of class FEFaceValues or of class FESubfaceValues (both derived from FEFaceValuesBase) according to the remaining cases mentioned in the introduction:

2. neighboring cell is finer (current cell: FESubfaceValues, neighboring cell: FEFaceValues);

3. neighboring cell is of the same refinement level (both, current and neighboring cell: FEFaceValues);

4. neighboring cell is coarser (current cell: FEFaceValues, neighboring cell: FESubfaceValues).

If we considered globally refined meshes then only case 3 would occur. But as we consider also locally refined meshes we need to distinguish all four cases making the following assembling function a bit longish.

 template <int dim>
 void DGMethod<dim>::assemble_system1 () 
 {
   const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
   std::vector<unsigned int> dofs (dofs_per_cell);
   std::vector<unsigned int> dofs_neighbor (dofs_per_cell);

First we create the update_flags for the FEValues and the FEFaceValues objects.

   const UpdateFlags update_flags = update_values
                                    | update_gradients
                                    | update_quadrature_points
                                    | update_JxW_values;

Note, that on faces we do not need gradients but we need normal vectors.

   const UpdateFlags face_update_flags = update_values
                                         | update_quadrature_points
                                         | update_JxW_values
                                         | update_normal_vectors;

On the neighboring cell we only need the shape values. Given a specific face, the quadrature points and `JxW values' are the same as for the current cells, the normal vectors are known to be the negative of the normal vectors of the current cell.

   const UpdateFlags neighbor_face_update_flags = update_values;

Then we create the FEValues object. Note, that since version 3.2.0 of deal.II the constructor of this class takes a Mapping object as first argument. Although the constructor without Mapping argument is still supported it is recommended to use the new constructor. This reduces the effect of `hidden magic' (the old constructor implicitely assumes a MappingQ1 mapping) and makes it easier to change the mapping object later.

   FEValues<dim> fe_v (
     mapping, fe, quadrature, update_flags);

Similarly we create the FEFaceValues and FESubfaceValues objects for both, the current and the neighboring cell. Within the following nested loop over all cells and all faces of the cell they will be reinited to the current cell and the face (and subface) number.

   FEFaceValues<dim> fe_v_face (
     mapping, fe, face_quadrature, face_update_flags);
   FESubfaceValues<dim> fe_v_subface (
     mapping, fe, face_quadrature, face_update_flags);
   FEFaceValues<dim> fe_v_face_neighbor (
     mapping, fe, face_quadrature, neighbor_face_update_flags);
   FESubfaceValues<dim> fe_v_subface_neighbor (
     mapping, fe, face_quadrature, neighbor_face_update_flags);

Now we create the cell matrices and vectors. Here we need two cell matrices, both for face terms that include test functions vi (internal shape functions, i.e. shape functions of the current cell). To be more precise, the first matrix will include the `ui and vi terms' and the second will include the `ue and vi terms'. Here we recall the convention that `ui' is the shortcut for $u_h^+$ and `ue' represents $u_h^-$, see the introduction.

   FullMatrix<double> ui_vi_matrix (dofs_per_cell, dofs_per_cell);
   FullMatrix<double> ue_vi_matrix (dofs_per_cell, dofs_per_cell);
 
   Vector<double>  cell_vector (dofs_per_cell);

Furthermore we need some cell iterators.

   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();

Now we start the loop over all active cells.

   for (;cell!=endc; ++cell) 
     {

In the assemble_face_term1 function contributions to the cell matrices and the cell vector are only ADDED. Therefore on each cell we need to reset the ui_vi_matrix and cell_vector to zero, before assembling the cell terms.

       ui_vi_matrix = 0;
       cell_vector = 0;

Now we reinit the FEValues object for the current cell

       fe_v.reinit (cell);

and call the function that assembles the cell terms. The first argument is the FEValues that was previously reinit'ed on the current cell.

       dg.assemble_cell_term(fe_v,
                             ui_vi_matrix,
                             cell_vector);

As in previous examples the vector `dofs' includes the dof_indices.

       cell->get_dof_indices (dofs);

This is the start of the nested loop over all faces.

       for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
         {

First we set the face iterator

           typename DoFHandler<dim>::face_iterator face=cell->face(face_no);

and clear the ue_vi_matrix on each face.

           ue_vi_matrix = 0;

Now we distinguish the four different cases in the ordering mentioned above. We start with faces belonging to the boundary of the domain.

           if (face->at_boundary())
             {

We reinit the FEFaceValues object to the current face

               fe_v_face.reinit (cell, face_no);

and assemble the corresponding face terms.

               dg.assemble_boundary_term(fe_v_face,
                                         ui_vi_matrix,
                                         cell_vector);
             }
           else
             {

Now we are not on the boundary of the domain, therefore there must exist a neighboring cell.

               typename DoFHandler<dim>::cell_iterator neighbor=
                 cell->neighbor(face_no);;

We proceed with the second and most complicated case: the neighboring cell is more refined than the current cell. As in deal.II neighboring cells are restricted to have a level difference of not more than one, the neighboring cell is known to be at most ONCE more refined than the current cell. Furthermore also the face is more refined, i.e. it has children. Here we note that the following part of code will not work for dim==1.

               if (face->has_children())
                 {

First we store which number the current cell has in the list of neighbors of the neighboring cell. Hence, neighbor->neighbor(neighbor2) equals the current cell cell.

                   const unsigned int neighbor2=
                     cell->neighbor_of_neighbor(face_no);

We loop over subfaces

                   for (unsigned int subface_no=0;
                        subface_no<face->n_children(); ++subface_no)
                     {

and set the cell iterator neighbor_child to the cell placed `behind' the current subface.

                       typename DoFHandler<dim>::active_cell_iterator
                         neighbor_child
                         = cell->neighbor_child_on_subface (face_no, subface_no);
                       
                       Assert (!neighbor_child->has_children(), ExcInternalError());

We need to reset the ue_vi_matrix on each subface because on each subface the un belong to different neighboring cells.

                       ue_vi_matrix = 0;

As already mentioned above for the current case (case 2) we employ the FESubfaceValues of the current cell (here reinited for the current cell, face and subface) and we employ the FEFaceValues of the neighboring child cell.

                       fe_v_subface.reinit (cell, face_no, subface_no);
                       fe_v_face_neighbor.reinit (neighbor_child, neighbor2);
 
                       dg.assemble_face_term1(fe_v_subface,
                                              fe_v_face_neighbor,
                                              ui_vi_matrix,
                                              ue_vi_matrix);

Then we get the dof indices of the neighbor_child cell

                       neighbor_child->get_dof_indices (dofs_neighbor);

and distribute ue_vi_matrix to the system_matrix

                       for (unsigned int i=0; i<dofs_per_cell; ++i)
                         for (unsigned int k=0; k<dofs_per_cell; ++k)
                           system_matrix.add(dofs[i], dofs_neighbor[k],
                                             ue_vi_matrix(i,k));
                     }

End of if (face->has_children())

                 }
               else
                 {

We proceed with case 3, i.e. neighboring cell is of the same refinement level as the current cell.

                   if (neighbor->level() == cell->level()) 
                     {

Like before we store which number the current cell has in the list of neighbors of the neighboring cell.

                       const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);

We reinit the FEFaceValues of the current and neighboring cell to the current face and assemble the corresponding face terms.

                       fe_v_face.reinit (cell, face_no);
                       fe_v_face_neighbor.reinit (neighbor, neighbor2);
                       
                       dg.assemble_face_term1(fe_v_face,
                                              fe_v_face_neighbor,
                                              ui_vi_matrix,
                                              ue_vi_matrix);

End of if (neighbor->level() == cell->level())

                     }
                   else
                     {

Finally we consider case 4. When the neighboring cell is not finer and not of the same refinement level as the current cell it must be coarser.

                       Assert(neighbor->level() < cell->level(), ExcInternalError());

Find out the how many'th face_no and subface_no the current face is w.r.t. the neighboring cell.

                       const std::pair<unsigned int, unsigned int> faceno_subfaceno=
                         cell->neighbor_of_coarser_neighbor(face_no);
                       const unsigned int neighbor_face_no=faceno_subfaceno.first,
                                       neighbor_subface_no=faceno_subfaceno.second;
 
                       Assert (neighbor->neighbor_child_on_subface (neighbor_face_no,
                                                                    neighbor_subface_no)
                               == cell,
                               ExcInternalError());

Reinit the appropriate FEFaceValues and assemble the face terms.

                       fe_v_face.reinit (cell, face_no);
                       fe_v_subface_neighbor.reinit (neighbor, neighbor_face_no,
                                                     neighbor_subface_no);
                       
                       dg.assemble_face_term1(fe_v_face,
                                              fe_v_subface_neighbor,
                                              ui_vi_matrix,
                                              ue_vi_matrix);
                     }

Now we get the dof indices of the neighbor_child cell,

                   neighbor->get_dof_indices (dofs_neighbor);

and distribute the ue_vi_matrix.

                   for (unsigned int i=0; i<dofs_per_cell; ++i)
                     for (unsigned int k=0; k<dofs_per_cell; ++k)
                       system_matrix.add(dofs[i], dofs_neighbor[k],
                                         ue_vi_matrix(i,k));
                 }

End of face not at boundary:

             }

End of loop over all faces:

         }

Finally we distribute the ui_vi_matrix

       for (unsigned int i=0; i<dofs_per_cell; ++i)
         for (unsigned int j=0; j<dofs_per_cell; ++j)
           system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i,j));

and the cell vector.

       for (unsigned int i=0; i<dofs_per_cell; ++i)
         right_hand_side(dofs[i]) += cell_vector(i);
     }
 }

Function: assemble_system2

We proceed with the assemble_system2 function that implements the DG discretization in its second version. This function is very similar to the assemble_system1 function. Therefore, here we only discuss the differences between the two functions. This function repeatedly calls the assemble_face_term2 function of the DGTransportEquation object, that assembles the face terms written as a sum of integrals over all faces. Therefore, we need to make sure that each face is treated only once. This is achieved by introducing the rule:

a) If the current and the neighboring cells are of the same refinement level we access and treat the face from the cell with lower index.

b) If the two cells are of different refinement levels we access and treat the face from the coarser cell.

Due to rule b) we do not need to consider case 4 (neighboring cell is coarser) any more.

 template <int dim>
 void DGMethod<dim>::assemble_system2 () 
 {
   const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
   std::vector<unsigned int> dofs (dofs_per_cell);
   std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
 
   const UpdateFlags update_flags = update_values
                                    | update_gradients
                                    | update_quadrature_points
                                    | update_JxW_values;
   
   const UpdateFlags face_update_flags = update_values
                                         | update_quadrature_points
                                         | update_JxW_values
                                         | update_normal_vectors;
   
   const UpdateFlags neighbor_face_update_flags = update_values;

Here we do not need fe_v_face_neighbor as case 4 does not occur.

   FEValues<dim> fe_v (
     mapping, fe, quadrature, update_flags);
   FEFaceValues<dim> fe_v_face (
     mapping, fe, face_quadrature, face_update_flags);
   FESubfaceValues<dim> fe_v_subface (
     mapping, fe, face_quadrature, face_update_flags);
   FEFaceValues<dim> fe_v_face_neighbor (
     mapping, fe, face_quadrature, neighbor_face_update_flags);
 
 
   FullMatrix<double> ui_vi_matrix (dofs_per_cell, dofs_per_cell);
   FullMatrix<double> ue_vi_matrix (dofs_per_cell, dofs_per_cell);

Additionally we need the following two cell matrices, both for face term that include test function ve (external shape functions, i.e. shape functions of the neighboring cell). To be more precise, the first matrix will include the `u and vn terms' and the second that will include the `un and vn terms'.

   FullMatrix<double> ui_ve_matrix (dofs_per_cell, dofs_per_cell);
   FullMatrix<double> ue_ve_matrix (dofs_per_cell, dofs_per_cell);
   
   Vector<double>  cell_vector (dofs_per_cell);

The following lines are roughly the same as in the previous function.

   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
   for (;cell!=endc; ++cell) 
     {
       ui_vi_matrix = 0;
       cell_vector = 0;
 
       fe_v.reinit (cell);
 
       dg.assemble_cell_term(fe_v,
                             ui_vi_matrix,
                             cell_vector);
       
       cell->get_dof_indices (dofs);
 
       for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
         {
           typename DoFHandler<dim>::face_iterator face=
             cell->face(face_no);

Case 1:

           if (face->at_boundary())
             {
               fe_v_face.reinit (cell, face_no);
 
               dg.assemble_boundary_term(fe_v_face,
                                         ui_vi_matrix,
                                         cell_vector);
             }
           else
             {
               Assert (cell->neighbor(face_no).state() == IteratorState::valid,
                       ExcInternalError());
               typename DoFHandler<dim>::cell_iterator neighbor=
                 cell->neighbor(face_no);

Case 2:

               if (face->has_children())
                 {
                   const unsigned int neighbor2=
                     cell->neighbor_of_neighbor(face_no);
                   
                   for (unsigned int subface_no=0;
                        subface_no<face->n_children(); ++subface_no)
                     {
                       typename DoFHandler<dim>::cell_iterator neighbor_child
                         = cell->neighbor_child_on_subface (face_no, subface_no);
                       Assert (!neighbor_child->has_children(), ExcInternalError());
                       
                       ue_vi_matrix = 0;
                       ui_ve_matrix = 0;
                       ue_ve_matrix = 0;
                       
                       fe_v_subface.reinit (cell, face_no, subface_no);
                       fe_v_face_neighbor.reinit (neighbor_child, neighbor2);
 
                       dg.assemble_face_term2(fe_v_subface,
                                              fe_v_face_neighbor,
                                              ui_vi_matrix,
                                              ue_vi_matrix,
                                              ui_ve_matrix,
                                              ue_ve_matrix);
                   
                       neighbor_child->get_dof_indices (dofs_neighbor);
                                                                       
                       for (unsigned int i=0; i<dofs_per_cell; ++i)
                         for (unsigned int j=0; j<dofs_per_cell; ++j)
                           {
                             system_matrix.add(dofs[i], dofs_neighbor[j],
                                               ue_vi_matrix(i,j));
                             system_matrix.add(dofs_neighbor[i], dofs[j],
                                               ui_ve_matrix(i,j));
                             system_matrix.add(dofs_neighbor[i], dofs_neighbor[j],
                                               ue_ve_matrix(i,j));
                           }
                     }
                 }
               else
                 {

Case 3, with the additional rule a)

                   if (neighbor->level() == cell->level() &&
                       neighbor->index() > cell->index()) 
                     {
                       const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
                       
                       ue_vi_matrix = 0;
                       ui_ve_matrix = 0;
                       ue_ve_matrix = 0;
                       
                       fe_v_face.reinit (cell, face_no);
                       fe_v_face_neighbor.reinit (neighbor, neighbor2);
                       
                       dg.assemble_face_term2(fe_v_face,
                                              fe_v_face_neighbor,
                                              ui_vi_matrix,
                                              ue_vi_matrix,
                                              ui_ve_matrix,
                                              ue_ve_matrix);
 
                       neighbor->get_dof_indices (dofs_neighbor);
 
                       for (unsigned int i=0; i<dofs_per_cell; ++i)
                         for (unsigned int j=0; j<dofs_per_cell; ++j)
                           {
                             system_matrix.add(dofs[i], dofs_neighbor[j],
                                               ue_vi_matrix(i,j));
                             system_matrix.add(dofs_neighbor[i], dofs[j],
                                               ui_ve_matrix(i,j));
                             system_matrix.add(dofs_neighbor[i], dofs_neighbor[j],
                                               ue_ve_matrix(i,j));
                           }
                     }

Due to rule b) we do not need to consider case 4.

                 }
             }
         }
       
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         for (unsigned int j=0; j<dofs_per_cell; ++j)
           system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i,j));
       
       for (unsigned int i=0; i<dofs_per_cell; ++i)
         right_hand_side(dofs[i]) += cell_vector(i);
     }
 }

All the rest

For this simple problem we use the simplest possible solver, called Richardson iteration, that represents a simple defect correction. This, in combination with a block SSOR preconditioner, that uses the special block matrix structure of system matrices arising from DG discretizations. The size of these blocks are the number of DoFs per cell. Here, we use a SSOR preconditioning as we have not renumbered the DoFs according to the flow field. If the DoFs are renumbered downstream the flow, then a block Gauss-Seidel preconditioner (see the PreconditionBlockSOR class with relaxation=1) makes a much better job.

 template <int dim>
 void DGMethod<dim>::solve (Vector<double> &solution) 
 {
   SolverControl           solver_control (1000, 1e-12, false, false);
   SolverRichardson<>      solver (solver_control);

Here we create the preconditioner,

   PreconditionBlockSSOR<SparseMatrix<double> > preconditioner;

we assigned the matrix to it and set the right block size.

   preconditioner.initialize(system_matrix, fe.dofs_per_cell);

After these preparations we are ready to start the linear solver.

   solver.solve (system_matrix, solution, right_hand_side,
                 preconditioner);
 }

We refine the grid according to a very simple refinement criterion, namely an approximation to the gradient of the solution. As here we consider the DG(1) method (i.e. we use piecewise bilinear shape functions) we could simply compute the gradients on each cell. But we do not want to base our refinement indicator on the gradients on each cell only, but want to base them also on jumps of the discontinuous solution function over faces between neighboring cells. The simpliest way of doing that is to compute approximative gradients by difference quotients including the cell under consideration and its neighbors. This is done by the DerivativeApproximation class that computes the approximate gradients in a way similar to the GradientEstimation described in step-9 of this tutorial. In fact, the DerivativeApproximation class was developed following the GradientEstimation class of step-9. Relating to the discussion in step-9, here we consider $h^{1+d/2}|\nabla_h u_h|$. Futhermore we note that we do not consider approximate second derivatives because solutions to the linear advection equation are in general not in $H^2$ but in $H^1$ (to be more precise, in $H^1_\beta$) only.

 template <int dim>
 void DGMethod<dim>::refine_grid ()
 {

The DerivativeApproximation class computes the gradients to float precision. This is sufficient as they are approximate and serve as refinement indicators only.

   Vector<float> gradient_indicator (triangulation.n_active_cells());

Now the approximate gradients are computed

   DerivativeApproximation::approximate_gradient (mapping,
                                                  dof_handler,
                                                  solution2,
                                                  gradient_indicator);

and they are cell-wise scaled by the factor $h^{1+d/2}$

   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
   for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no)
     gradient_indicator(cell_no)*=std::pow(cell->diameter(), 1+1.0*dim/2);

Finally they serve as refinement indicator.

   GridRefinement::refine_and_coarsen_fixed_number (triangulation,
                                                    gradient_indicator,
                                                    0.3, 0.1);
 
   triangulation.execute_coarsening_and_refinement ();
 }

The output of this program consists of eps-files of the adaptively refined grids and the numerical solutions given in gnuplot format. This was covered in previous examples and will not be further commented on.

 template <int dim>
 void DGMethod<dim>::output_results (const unsigned int cycle) const
 {

Write the grid in eps format.

   std::string filename = "grid-";
   filename += ('0' + cycle);
   Assert (cycle < 10, ExcInternalError());
   
   filename += ".eps";
   std::cout << "Writing grid to <" << filename << ">..." << std::endl;
   std::ofstream eps_output (filename.c_str());
 
   GridOut grid_out;
   grid_out.write_eps (triangulation, eps_output);

Output of the solution in gnuplot format.

   filename = "sol-";
   filename += ('0' + cycle);
   Assert (cycle < 10, ExcInternalError());
   
   filename += ".gnuplot";
   std::cout << "Writing solution to <" << filename << ">..."
             << std::endl << std::endl;
   std::ofstream gnuplot_output (filename.c_str());
   
   DataOut<dim> data_out;
   data_out.attach_dof_handler (dof_handler);
   data_out.add_data_vector (solution2, "u");
 
   data_out.build_patches ();
   
   data_out.write_gnuplot(gnuplot_output);
 }

The following run function is similar to previous examples. The only difference is that the problem is assembled and solved twice on each refinement step; first by assemble_system1 that implements the first version and then by assemble_system2 that implements the second version of writing the DG discretization. Furthermore the time needed by each of the two assembling routines is measured.

 template <int dim>
 void DGMethod<dim>::run () 
 {
   for (unsigned int cycle=0; cycle<6; ++cycle)
     {
       std::cout << "Cycle " << cycle << ':' << std::endl;
 
       if (cycle == 0)
         {
           GridGenerator::hyper_cube (triangulation);
 
           triangulation.refine_global (3);
         }
       else
         refine_grid ();
       
 
       std::cout << "   Number of active cells:       "
                 << triangulation.n_active_cells()
                 << std::endl;
 
       setup_system ();
 
       std::cout << "   Number of degrees of freedom: "
                 << dof_handler.n_dofs()
                 << std::endl;

The constructor of the Timer class automatically starts the time measurement.

       Timer assemble_timer;

First assembling routine.

       assemble_system1 ();

The operator () accesses the current time without disturbing the time measurement.

       std::cout << "Time of assemble_system1: "
                 << assemble_timer()
                 << std::endl;
       solve (solution1);

As preparation for the second assembling routine we reinit the system matrix, the right hand side vector and the Timer object.

       system_matrix = 0;
       right_hand_side = 0;
       assemble_timer.reset();

We start the Timer,

       assemble_timer.start();

call the second assembling routine

       assemble_system2 ();

and access the current time.

       std::cout << "Time of assemble_system2: "
                 << assemble_timer()
                 << std::endl;
       solve (solution2);

To make sure that both versions of the DG method yield the same discretization and hence the same solution we check the two solutions for equality.

       solution1-=solution2;
       const double difference=solution1.linfty_norm();
       if (difference>1e-13)
         std::cout << "solution1 and solution2 differ!!" << std::endl;
       else
         std::cout << "solution1 and solution2 coincide." << std::endl;

Finally we perform the output.

       output_results (cycle);
     }
 }

The following main function is similar to previous examples and need not to be commented on.

 int main () 
 {
   try
     {
       DGMethod<2> dgmethod;
       dgmethod.run ();
     }
   catch (std::exception &exc)
     {
       std::cerr << std::endl << std::endl
                 << "----------------------------------------------------"
                 << std::endl;
       std::cerr << "Exception on processing: " << std::endl
                 << exc.what() << std::endl
                 << "Aborting!" << std::endl
                 << "----------------------------------------------------"
                 << std::endl;
       return 1;
     }
   catch (...) 
     {
       std::cerr << std::endl << std::endl
                 << "----------------------------------------------------"
                 << std::endl;
       std::cerr << "Unknown exception!" << std::endl
                 << "Aborting!" << std::endl
                 << "----------------------------------------------------"
                 << std::endl;
       return 1;
     };
   
   return 0;
 }

Results

The output of this program consist of the console output, the eps files including the grids, and the solutions given in gnuplot format.

Cycle 0:
   Number of active cells:       64
   Number of degrees of freedom: 256
Time of assemble_system1: 0.05
Time of assemble_system2: 0.04
solution1 and solution2 coincide.
Writing grid to <grid-0.eps>...
Writing solution to <sol-0.gnuplot>...

Cycle 1:
   Number of active cells:       112
   Number of degrees of freedom: 448
Time of assemble_system1: 0.09
Time of assemble_system2: 0.07
solution1 and solution2 coincide.
Writing grid to <grid-1.eps>...
Writing solution to <sol-1.gnuplot>...

Cycle 2:
   Number of active cells:       214
   Number of degrees of freedom: 856
Time of assemble_system1: 0.17
Time of assemble_system2: 0.14
solution1 and solution2 coincide.
Writing grid to <grid-2.eps>...
Writing solution to <sol-2.gnuplot>...

Cycle 3:
   Number of active cells:       415
   Number of degrees of freedom: 1660
Time of assemble_system1: 0.32
Time of assemble_system2: 0.28
solution1 and solution2 coincide.
Writing grid to <grid-3.eps>...
Writing solution to <sol-3.gnuplot>...

Cycle 4:
   Number of active cells:       796
   Number of degrees of freedom: 3184
Time of assemble_system1: 0.62
Time of assemble_system2: 0.52
solution1 and solution2 coincide.
Writing grid to <grid-4.eps>...
Writing solution to <sol-4.gnuplot>...

Cycle 5:
   Number of active cells:       1561
   Number of degrees of freedom: 6244
Time of assemble_system1: 1.23
Time of assemble_system2: 1.03
solution1 and solution2 coincide.
Writing grid to <grid-5.eps>...
Writing solution to <sol-5.gnuplot>...

We see that, as expected, on each refinement step the two solutions coincide. The difference measured in time of treating each face only once (second version of the DG method) in comparison with treating each face twice within a nested loop over all cells and all faces of each cell (first version), is much less than one might have expected. The gain is less than 20% on the last few refinement steps.

First we show the solutions on the initial mesh, the mesh after two and after five adaptive refinement steps.

step-12.sol-0.png
step-12.sol-2.png
step-12.sol-5.png

Then we show the final grid (after 5 refinement steps) and the solution again, this time with a nicer 3d rendering (obtained using the DataOutBase::write_vtk function and the VTK-based VisIt visualization program) that better shows the sharpness of the jump on the refined mesh and the over- and undershoots of the solution along the interface:

step-12.grid-5.png
step-12.3d-solution.png

And finally we show a plot of a 3d computation.

step-12.sol-5-3d.png

The plain program

(If you are looking at a locally installed deal.II version, then the program can be found at /build /buildd /deal.ii-6.2.1 /examples /step-12 /step-12.cc . Otherwise, this is only the path on some remote server.)

 / * Subversion Id:  step-12.cc 16430 2008-07-08 15:25:01Z hartmann  * /
 / * Author: Ralf Hartmann, University of Heidelberg, 2001 * /
 
 / *    Subversion Id:  step-12.cc 16430 2008-07-08 15:25:01Z hartmann        * /
 / *                                                                * /
 / *    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors * /
 / *                                                                * /
 / *    This file is subject to QPL and may not be  distributed     * /
 / *    without copyright and license information. Please refer     * /
 / *    to the file deal.II/doc/license.html for the  text  and     * /
 / *    further information on this license.                        * /
 
 #include <base/quadrature_lib.h>
 #include <base/function.h>
 #include <lac/vector.h>
 #include <lac/sparse_matrix.h>
 #include <grid/tria.h>
 #include <grid/grid_generator.h>
 #include <grid/grid_out.h>
 #include <grid/grid_refinement.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
 #include <fe/fe_values.h>
 #include <dofs/dof_handler.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
 #include <numerics/data_out.h>
 
 #include <fe/mapping_q1.h>
 #include <fe/fe_dgq.h>
 #include <lac/solver_richardson.h>
 #include <lac/precondition_block.h>
 #include <numerics/derivative_approximation.h>
 #include <base/timer.h>
 
 #include <iostream>
 #include <fstream>
 
 using namespace dealii;

 template <int dim>
 class RHS:  public Function<dim>
 {
   public:
     virtual void value_list (const std::vector<Point<dim> > &points,
                             std::vector<double> &values,
                             const unsigned int component=0) const;
 };
 
 
 template <int dim>
 class BoundaryValues:  public Function<dim>
 {
   public:
     virtual void value_list (const std::vector<Point<dim> > &points,
                             std::vector<double> &values,
                             const unsigned int component=0) const;
 };
 
 
 template <int dim>
 class Beta
 {
   public:
     Beta () {}
     void value_list (const std::vector<Point<dim> > &points,
                     std::vector<Point<dim> > &values) const;
 };
 
 
 template <int dim>
 void RHS<dim>::value_list(const std::vector<Point<dim> > &points,
                          std::vector<double> &values,
                          const unsigned int) const
 {
   Assert(values.size()==points.size(),
         ExcDimensionMismatch(values.size(),points.size()));
 
   for (unsigned int i=0; i<values.size(); ++i)
     values[i]=0;
 }
 
 
 template <int dim>
 void Beta<dim>::value_list(const std::vector<Point<dim> > &points,
                           std::vector<Point<dim> > &values) const
 {
   Assert(values.size()==points.size(),
         ExcDimensionMismatch(values.size(),points.size()));
 
   for (unsigned int i=0; i<points.size(); ++i)
     {
       values[i](0) = -points[i](1);
       values[i](1) = points[i](0);
       values[i] /= std::sqrt(values[i].square());
     }
 }
 
 
 template <int dim>
 void BoundaryValues<dim>::value_list(const std::vector<Point<dim> > &points,
                                       std::vector<double> &values,
                                       const unsigned int) const
 {
   Assert(values.size()==points.size(),
         ExcDimensionMismatch(values.size(),points.size()));
 
   for (unsigned int i=0; i<values.size(); ++i)
     {
       if (points[i](0)<0.5)
        values[i]=1.;
       else
        values[i]=0.;
     }
 }

 template <int dim>
 class DGTransportEquation
 {
   public:
     DGTransportEquation();
 
     void assemble_cell_term(const FEValues<dim>& fe_v,
                            FullMatrix<double> &ui_vi_matrix,
                            Vector<double> &cell_vector) const;
     
     void assemble_boundary_term(const FEFaceValues<dim>& fe_v,
                                FullMatrix<double> &ui_vi_matrix,
                                Vector<double> &cell_vector) const;
     
     void assemble_face_term1(const FEFaceValuesBase<dim>& fe_v,
                             const FEFaceValuesBase<dim>& fe_v_neighbor,
                             FullMatrix<double> &ui_vi_matrix,
                             FullMatrix<double> &ue_vi_matrix) const;
 
     void assemble_face_term2(const FEFaceValuesBase<dim>& fe_v,
                             const FEFaceValuesBase<dim>& fe_v_neighbor,
                             FullMatrix<double> &ui_vi_matrix,
                             FullMatrix<double> &ue_vi_matrix,
                             FullMatrix<double> &ui_ve_matrix,
                             FullMatrix<double> &ue_ve_matrix) const;
   private:
     const Beta<dim> beta_function;
     const RHS<dim> rhs_function;
     const BoundaryValues<dim> boundary_function;
 };
 
 
 template <int dim>
 DGTransportEquation<dim>::DGTransportEquation ()
                :
                beta_function (),
                rhs_function (),
                boundary_function ()
 {}

 template <int dim>
 void DGTransportEquation<dim>::assemble_cell_term(
   const FEValues<dim> &fe_v,
   FullMatrix<double> &ui_vi_matrix,
   Vector<double> &cell_vector) const
 {
   const std::vector<double> &JxW = fe_v.get_JxW_values ();
 
   std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
   std::vector<double> rhs (fe_v.n_quadrature_points);
   
   beta_function.value_list (fe_v.get_quadrature_points(), beta);
   rhs_function.value_list (fe_v.get_quadrature_points(), rhs);
   
   for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
     for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) 
       {
        for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
          ui_vi_matrix(i,j) -= beta[point]*fe_v.shape_grad(i,point)*
                              fe_v.shape_value(j,point) *
                              JxW[point];
        
        cell_vector(i) += rhs[point] * fe_v.shape_value(i,point) * JxW[point];
       }
 }

 template <int dim>
 void DGTransportEquation<dim>::assemble_boundary_term(
   const FEFaceValues<dim>& fe_v,    
   FullMatrix<double> &ui_vi_matrix,
   Vector<double> &cell_vector) const
 {
   const std::vector<double> &JxW = fe_v.get_JxW_values ();
   const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
 
   std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
   std::vector<double> g(fe_v.n_quadrature_points);
   
   beta_function.value_list (fe_v.get_quadrature_points(), beta);
   boundary_function.value_list (fe_v.get_quadrature_points(), g);
 
   for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
     {
       const double beta_n=beta[point] * normals[point];      
       if (beta_n>0)
        for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
          for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
            ui_vi_matrix(i,j) += beta_n *
                               fe_v.shape_value(j,point) *
                               fe_v.shape_value(i,point) *
                               JxW[point];
       else
        for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
          cell_vector(i) -= beta_n *
                            g[point] *
                            fe_v.shape_value(i,point) *
                            JxW[point];
     }
 }

 template <int dim>
 void DGTransportEquation<dim>::assemble_face_term1(
   const FEFaceValuesBase<dim>& fe_v,
   const FEFaceValuesBase<dim>& fe_v_neighbor,      
   FullMatrix<double> &ui_vi_matrix,
   FullMatrix<double> &ue_vi_matrix) const
 {
   const std::vector<double> &JxW = fe_v.get_JxW_values ();
   const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
 
   std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
   beta_function.value_list (fe_v.get_quadrature_points(), beta);
 
   for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
     {
       const double beta_n=beta[point] * normals[point];
       if (beta_n>0)
        for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
          for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
            ui_vi_matrix(i,j) += beta_n *
                               fe_v.shape_value(j,point) *
                               fe_v.shape_value(i,point) *
                               JxW[point];
       else
        for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
          for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
            ue_vi_matrix(i,k) += beta_n *
                                fe_v_neighbor.shape_value(k,point) *
                                fe_v.shape_value(i,point) *
                                JxW[point];
     }
 }

 template <int dim>
 void DGTransportEquation<dim>::assemble_face_term2(
   const FEFaceValuesBase<dim>& fe_v,
   const FEFaceValuesBase<dim>& fe_v_neighbor,      
   FullMatrix<double> &ui_vi_matrix,
   FullMatrix<double> &ue_vi_matrix,
   FullMatrix<double> &ui_ve_matrix,
   FullMatrix<double> &ue_ve_matrix) const
 {
   const std::vector<double> &JxW = fe_v.get_JxW_values ();
   const std::vector<Point<dim> > &normals = fe_v.get_normal_vectors ();
 
   std::vector<Point<dim> > beta (fe_v.n_quadrature_points);
   
   beta_function.value_list (fe_v.get_quadrature_points(), beta);
 
   for (unsigned int point=0; point<fe_v.n_quadrature_points; ++point)
     {
       const double beta_n=beta[point] * normals[point];
       if (beta_n>0)
        {
          for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
            for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
              ui_vi_matrix(i,j) += beta_n *
                                 fe_v.shape_value(j,point) *
                                 fe_v.shape_value(i,point) *
                                 JxW[point];
 
          for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
            for (unsigned int j=0; j<fe_v.dofs_per_cell; ++j)
              ui_ve_matrix(k,j) -= beta_n *
                                  fe_v.shape_value(j,point) *
                                  fe_v_neighbor.shape_value(k,point) *
                                  JxW[point];
        }
       else
        {
          for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
            for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l)
              ue_vi_matrix(i,l) += beta_n *
                                  fe_v_neighbor.shape_value(l,point) *
                                  fe_v.shape_value(i,point) *
                                  JxW[point];
 
          for (unsigned int k=0; k<fe_v_neighbor.dofs_per_cell; ++k)
            for (unsigned int l=0; l<fe_v_neighbor.dofs_per_cell; ++l)
              ue_ve_matrix(k,l) -= beta_n *
                                   fe_v_neighbor.shape_value(l,point) *
                                   fe_v_neighbor.shape_value(k,point) *
                                   JxW[point];
        }
     }
 }

 template <int dim>
 class DGMethod
 {
   public:
     DGMethod ();
     ~DGMethod ();
 
     void run ();
     
   private:
     void setup_system ();
     void assemble_system1 ();
     void assemble_system2 ();
     void solve (Vector<double> &solution);
     void refine_grid ();
     void output_results (const unsigned int cycle) const;
     
     Triangulation<dim>   triangulation;
     const MappingQ1<dim> mapping;
     
     FE_DGQ<dim>          fe;
     DoFHandler<dim>      dof_handler;
 
     SparsityPattern      sparsity_pattern;
     SparseMatrix<double> system_matrix;
 
     const QGauss<dim>   quadrature;
     const QGauss<dim-1> face_quadrature;
     
     Vector<double>       solution1;
     Vector<double>       solution2;
     Vector<double>       right_hand_side;
     
     const DGTransportEquation<dim> dg;
 };
 
 
 template <int dim>
 DGMethod<dim>::DGMethod ()
                :
                mapping (),
                 fe (1),
                dof_handler (triangulation),
                quadrature (4),
                face_quadrature (4),
                dg ()
 {}
 
 
 template <int dim>
 DGMethod<dim>::~DGMethod () 
 {
   dof_handler.clear ();
 }
 
 
 template <int dim>
 void DGMethod<dim>::setup_system ()
 {
   dof_handler.distribute_dofs (fe);
 
   sparsity_pattern.reinit (dof_handler.n_dofs(),
                           dof_handler.n_dofs(),
                           (GeometryInfo<dim>::faces_per_cell
                            *GeometryInfo<dim>::max_children_per_face+1)*fe.dofs_per_cell);
   
   DoFTools::make_flux_sparsity_pattern (dof_handler, sparsity_pattern);
   
   sparsity_pattern.compress();
   
   system_matrix.reinit (sparsity_pattern);
 
   solution1.reinit (dof_handler.n_dofs());
   solution2.reinit (dof_handler.n_dofs());
   right_hand_side.reinit (dof_handler.n_dofs());
 }

 template <int dim>
 void DGMethod<dim>::assemble_system1 () 
 {
   const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
   std::vector<unsigned int> dofs (dofs_per_cell);
   std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
 
   const UpdateFlags update_flags = update_values
                                    | update_gradients
                                    | update_quadrature_points
                                    | update_JxW_values;
 
   const UpdateFlags face_update_flags = update_values
                                         | update_quadrature_points
                                         | update_JxW_values
                                         | update_normal_vectors;
   
   const UpdateFlags neighbor_face_update_flags = update_values;
    
   FEValues<dim> fe_v (
     mapping, fe, quadrature, update_flags);
   
   FEFaceValues<dim> fe_v_face (
     mapping, fe, face_quadrature, face_update_flags);
   FESubfaceValues<dim> fe_v_subface (
     mapping, fe, face_quadrature, face_update_flags);
   FEFaceValues<dim> fe_v_face_neighbor (
     mapping, fe, face_quadrature, neighbor_face_update_flags);
   FESubfaceValues<dim> fe_v_subface_neighbor (
     mapping, fe, face_quadrature, neighbor_face_update_flags);
 
   FullMatrix<double> ui_vi_matrix (dofs_per_cell, dofs_per_cell);
   FullMatrix<double> ue_vi_matrix (dofs_per_cell, dofs_per_cell);
 
   Vector<double>  cell_vector (dofs_per_cell);
 
   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
 
   for (;cell!=endc; ++cell) 
     {
       ui_vi_matrix = 0;
       cell_vector = 0;
       
       fe_v.reinit (cell);
 
       dg.assemble_cell_term(fe_v,
                            ui_vi_matrix,
                            cell_vector);
 
       cell->get_dof_indices (dofs);
 
       for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
        {
          typename DoFHandler<dim>::face_iterator face=cell->face(face_no);
          
          ue_vi_matrix = 0;
                  
          if (face->at_boundary())
            {
              fe_v_face.reinit (cell, face_no);
 
              dg.assemble_boundary_term(fe_v_face,
                                        ui_vi_matrix,
                                        cell_vector);
            }
          else
            {
              typename DoFHandler<dim>::cell_iterator neighbor=
                cell->neighbor(face_no);;
 
              if (face->has_children())
                {
                  const unsigned int neighbor2=
                    cell->neighbor_of_neighbor(face_no);
                  
                  
                  for (unsigned int subface_no=0;
                       subface_no<face->n_children(); ++subface_no)
                    {
                      typename DoFHandler<dim>::active_cell_iterator
                         neighbor_child
                         = cell->neighbor_child_on_subface (face_no, subface_no);
                      
                      Assert (!neighbor_child->has_children(), ExcInternalError());
 
                      ue_vi_matrix = 0;
                      
                      fe_v_subface.reinit (cell, face_no, subface_no);
                      fe_v_face_neighbor.reinit (neighbor_child, neighbor2);
 
                      dg.assemble_face_term1(fe_v_subface,
                                             fe_v_face_neighbor,
                                             ui_vi_matrix,
                                             ue_vi_matrix);
                      
                      neighbor_child->get_dof_indices (dofs_neighbor);
                                                                
                      for (unsigned int i=0; i<dofs_per_cell; ++i)
                        for (unsigned int k=0; k<dofs_per_cell; ++k)
                          system_matrix.add(dofs[i], dofs_neighbor[k],
                                            ue_vi_matrix(i,k));
                    }
                }
              else
                {
                  if (neighbor->level() == cell->level()) 
                    {
                      const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
 
                      fe_v_face.reinit (cell, face_no);
                      fe_v_face_neighbor.reinit (neighbor, neighbor2);
                      
                      dg.assemble_face_term1(fe_v_face,
                                             fe_v_face_neighbor,
                                             ui_vi_matrix,
                                             ue_vi_matrix);
                    }
                  else
                    {
                      Assert(neighbor->level() < cell->level(), ExcInternalError());
 
                      const std::pair<unsigned int, unsigned int> faceno_subfaceno=
                        cell->neighbor_of_coarser_neighbor(face_no);
                      const unsigned int neighbor_face_no=faceno_subfaceno.first,
                                      neighbor_subface_no=faceno_subfaceno.second;
 
                      Assert (neighbor->neighbor_child_on_subface (neighbor_face_no,
                                                                    neighbor_subface_no)
                               == cell,
                               ExcInternalError());
 
                      fe_v_face.reinit (cell, face_no);
                      fe_v_subface_neighbor.reinit (neighbor, neighbor_face_no,
                                                    neighbor_subface_no);
                      
                      dg.assemble_face_term1(fe_v_face,
                                             fe_v_subface_neighbor,
                                             ui_vi_matrix,
                                             ue_vi_matrix);
                    }
 
                  neighbor->get_dof_indices (dofs_neighbor);
                                                                
                  for (unsigned int i=0; i<dofs_per_cell; ++i)
                    for (unsigned int k=0; k<dofs_per_cell; ++k)
                      system_matrix.add(dofs[i], dofs_neighbor[k],
                                        ue_vi_matrix(i,k));
                }
            }
        }
       
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        for (unsigned int j=0; j<dofs_per_cell; ++j)
          system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i,j));
       
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        right_hand_side(dofs[i]) += cell_vector(i);
     }
 }

 template <int dim>
 void DGMethod<dim>::assemble_system2 () 
 {
   const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
   std::vector<unsigned int> dofs (dofs_per_cell);
   std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
 
   const UpdateFlags update_flags = update_values
                                    | update_gradients
                                    | update_quadrature_points
                                    | update_JxW_values;
   
   const UpdateFlags face_update_flags = update_values
                                         | update_quadrature_points
                                         | update_JxW_values
                                         | update_normal_vectors;
   
   const UpdateFlags neighbor_face_update_flags = update_values;
 
   FEValues<dim> fe_v (
     mapping, fe, quadrature, update_flags);
   FEFaceValues<dim> fe_v_face (
     mapping, fe, face_quadrature, face_update_flags);
   FESubfaceValues<dim> fe_v_subface (
     mapping, fe, face_quadrature, face_update_flags);
   FEFaceValues<dim> fe_v_face_neighbor (
     mapping, fe, face_quadrature, neighbor_face_update_flags);
 
 
   FullMatrix<double> ui_vi_matrix (dofs_per_cell, dofs_per_cell);
   FullMatrix<double> ue_vi_matrix (dofs_per_cell, dofs_per_cell);
   
   FullMatrix<double> ui_ve_matrix (dofs_per_cell, dofs_per_cell);
   FullMatrix<double> ue_ve_matrix (dofs_per_cell, dofs_per_cell);
   
   Vector<double>  cell_vector (dofs_per_cell);
 
   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
   for (;cell!=endc; ++cell) 
     {
       ui_vi_matrix = 0;
       cell_vector = 0;
 
       fe_v.reinit (cell);
 
       dg.assemble_cell_term(fe_v,
                            ui_vi_matrix,
                            cell_vector);
       
       cell->get_dof_indices (dofs);
 
       for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
        {
          typename DoFHandler<dim>::face_iterator face=
            cell->face(face_no);
 
          if (face->at_boundary())
            {
              fe_v_face.reinit (cell, face_no);
 
              dg.assemble_boundary_term(fe_v_face,
                                        ui_vi_matrix,
                                        cell_vector);
            }
          else
            {
              Assert (cell->neighbor(face_no).state() == IteratorState::valid,
                      ExcInternalError());
              typename DoFHandler<dim>::cell_iterator neighbor=
                cell->neighbor(face_no);
              if (face->has_children())
                {
                  const unsigned int neighbor2=
                    cell->neighbor_of_neighbor(face_no);
                  
                  for (unsigned int subface_no=0;
                       subface_no<face->n_children(); ++subface_no)
                    {
                      typename DoFHandler<dim>::cell_iterator neighbor_child
                         = cell->neighbor_child_on_subface (face_no, subface_no);
                      Assert (!neighbor_child->has_children(), ExcInternalError());
                      
                      ue_vi_matrix = 0;
                      ui_ve_matrix = 0;
                      ue_ve_matrix = 0;
                      
                      fe_v_subface.reinit (cell, face_no, subface_no);
                      fe_v_face_neighbor.reinit (neighbor_child, neighbor2);
 
                      dg.assemble_face_term2(fe_v_subface,
                                             fe_v_face_neighbor,
                                             ui_vi_matrix,
                                             ue_vi_matrix,
                                             ui_ve_matrix,
                                             ue_ve_matrix);
                  
                      neighbor_child->get_dof_indices (dofs_neighbor);
                                                                
                      for (unsigned int i=0; i<dofs_per_cell; ++i)
                        for (unsigned int j=0; j<dofs_per_cell; ++j)
                          {
                            system_matrix.add(dofs[i], dofs_neighbor[j],
                                              ue_vi_matrix(i,j));
                            system_matrix.add(dofs_neighbor[i], dofs[j],
                                              ui_ve_matrix(i,j));
                            system_matrix.add(dofs_neighbor[i], dofs_neighbor[j],
                                              ue_ve_matrix(i,j));
                          }
                    }
                }
              else
                {
                  if (neighbor->level() == cell->level() &&
                      neighbor->index() > cell->index()) 
                    {
                      const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
                      
                      ue_vi_matrix = 0;
                      ui_ve_matrix = 0;
                      ue_ve_matrix = 0;
                      
                      fe_v_face.reinit (cell, face_no);
                      fe_v_face_neighbor.reinit (neighbor, neighbor2);
                      
                      dg.assemble_face_term2(fe_v_face,
                                             fe_v_face_neighbor,
                                             ui_vi_matrix,
                                             ue_vi_matrix,
                                             ui_ve_matrix,
                                             ue_ve_matrix);
 
                      neighbor->get_dof_indices (dofs_neighbor);
 
                      for (unsigned int i=0; i<dofs_per_cell; ++i)
                        for (unsigned int j=0; j<dofs_per_cell; ++j)
                          {
                            system_matrix.add(dofs[i], dofs_neighbor[j],
                                              ue_vi_matrix(i,j));
                            system_matrix.add(dofs_neighbor[i], dofs[j],
                                              ui_ve_matrix(i,j));
                            system_matrix.add(dofs_neighbor[i], dofs_neighbor[j],
                                              ue_ve_matrix(i,j));
                          }
                    }
 
                }
            }
        }
       
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        for (unsigned int j=0; j<dofs_per_cell; ++j)
          system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i,j));
       
       for (unsigned int i=0; i<dofs_per_cell; ++i)
        right_hand_side(dofs[i]) += cell_vector(i);
     }
 }

 template <int dim>
 void DGMethod<dim>::solve (Vector<double> &solution) 
 {
   SolverControl           solver_control (1000, 1e-12, false, false);
   SolverRichardson<>      solver (solver_control);
 
   PreconditionBlockSSOR<SparseMatrix<double> > preconditioner;
 
   preconditioner.initialize(system_matrix, fe.dofs_per_cell);
 
   solver.solve (system_matrix, solution, right_hand_side,
                preconditioner);
 }
 
 
 template <int dim>
 void DGMethod<dim>::refine_grid ()
 {
   Vector<float> gradient_indicator (triangulation.n_active_cells());
 
   DerivativeApproximation::approximate_gradient (mapping,
                                                 dof_handler,
                                                 solution2,
                                                 gradient_indicator);
 
   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
   for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no)
     gradient_indicator(cell_no)*=std::pow(cell->diameter(), 1+1.0*dim/2);
 
   GridRefinement::refine_and_coarsen_fixed_number (triangulation,
                                                   gradient_indicator,
                                                   0.3, 0.1);
 
   triangulation.execute_coarsening_and_refinement ();
 }
 
 
 template <int dim>
 void DGMethod<dim>::output_results (const unsigned int cycle) const
 {
   std::string filename = "grid-";
   filename += ('0' + cycle);
   Assert (cycle < 10, ExcInternalError());
   
   filename += ".eps";
   std::cout << "Writing grid to <" << filename << ">..." << std::endl;
   std::ofstream eps_output (filename.c_str());
 
   GridOut grid_out;
   grid_out.write_eps (triangulation, eps_output);
   
   filename = "sol-";
   filename += ('0' + cycle);
   Assert (cycle < 10, ExcInternalError());
   
   filename += ".gnuplot";
   std::cout << "Writing solution to <" << filename << ">..."
            << std::endl << std::endl;
   std::ofstream gnuplot_output (filename.c_str());
   
   DataOut<dim> data_out;
   data_out.attach_dof_handler (dof_handler);
   data_out.add_data_vector (solution2, "u");
 
   data_out.build_patches ();
   
   data_out.write_gnuplot(gnuplot_output);
 }
 
 
 template <int dim>
 void DGMethod<dim>::run () 
 {
   for (unsigned int cycle=0; cycle<6; ++cycle)
     {
       std::cout << "Cycle " << cycle << ':' << std::endl;
 
       if (cycle == 0)
        {
          GridGenerator::hyper_cube (triangulation);
 
          triangulation.refine_global (3);
        }
       else
        refine_grid ();
       
 
       std::cout << "   Number of active cells:       "
                << triangulation.n_active_cells()
                << std::endl;
 
       setup_system ();
 
       std::cout << "   Number of degrees of freedom: "
                << dof_handler.n_dofs()
                << std::endl;
 
       Timer assemble_timer;
       assemble_system1 ();
       std::cout << "Time of assemble_system1: "
                << assemble_timer()
                << std::endl;
       solve (solution1);
 
       system_matrix = 0;
       right_hand_side = 0;
       assemble_timer.reset();
 
       assemble_timer.start();
       assemble_system2 ();
       std::cout << "Time of assemble_system2: "
                << assemble_timer()
                << std::endl;
       solve (solution2);
 
       solution1-=solution2;
       const double difference=solution1.linfty_norm();
       if (difference>1e-13)
        std::cout << "solution1 and solution2 differ!!" << std::endl;
       else
        std::cout << "solution1 and solution2 coincide." << std::endl;
        
       output_results (cycle);
     }
 }
 
 int main () 
 {
   try
     {
       DGMethod<2> dgmethod;
       dgmethod.run ();
     }
   catch (std::exception &exc)
     {
       std::cerr << std::endl << std::endl
                << "----------------------------------------------------"
                << std::endl;
       std::cerr << "Exception on processing: " << std::endl
                << exc.what() << std::endl
                << "Aborting!" << std::endl
                << "----------------------------------------------------"
                << std::endl;
       return 1;
     }
   catch (...) 
     {
       std::cerr << std::endl << std::endl
                << "----------------------------------------------------"
                << std::endl;
       std::cerr << "Unknown exception!" << std::endl
                << "Aborting!" << std::endl
                << "----------------------------------------------------"
                << std::endl;
       return 1;
     };
   
   return 0;
 }

deal.II documentation generated on Mon Nov 23 22:56:38 2009 by doxygen 1.6.1