<h3>Overview</h3>
This example is devoted to the MeshWorker framework and the <em>discontinuous
-Galerkin method</em>, or in short: DG method. It solves the same problem as
-@ref step_12 "step-12" (see there for a description of the problem and
-discretization), but here we use the MeshWorker framework in order to save
-programming the cell/face loops that are often rather. The aim of the
-MeshWorker framework is to simplify this process, by putting the majority of
-the boring setup into a framework class and leaving to user code only things
-that are specific to the application. We have tried to strip this example of
-peripheral information such that the structure becomes more clear.
+Galerkin method</em>, or in short: DG method. It includes the following topics.
+<ol>
+ <li> Discretization of the linear transport equation with the DG method.
+ <li> Assembling of the system matrix using the MeshWorker::loop().
+</ol>
The particular concern of this program are the loops of DG methods. These turn
out to be especially complex, primarily because for the face terms, we have to
classes from the MeshWorker namespace and combine them to achieve your
goal.
+<h3>Problem</h3>
+
+The model problem solved in this example is the linear advection equation
+<a name="step-12.transport-equation">@f[
+ \nabla\cdot \left({\mathbf \beta} u\right)=f \qquad\mbox{in }\Omega,
+\qquad\qquad\qquad\mathrm{[transport-equation]}@f]</a>
+subject to the boundary conditions
+@f[
+u=g\quad\mbox{on }\Gamma_-,
+@f]
+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,
+@f[
+\Gamma_-:=\{{\bf x}\in\Gamma, {\mathbf \beta}({\bf x})\cdot{\bf n}({\bf x})<0\}
+@f]
+the inflow part of the boundary of the domain and ${\bf n}$ denotes
+the unit outward normal to the boundary $\Gamma$. Equation
+<a href="#step-12.transport-equation">[transport-equation]</a> is the conservative version of the
+transport equation already considered in step 9 of this tutorial.
+
+In particular, we consider problem <a href="#step-12.transport-equation">[transport-equation]</a> 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$.
+
+
+<h3>Discretization</h3>
+
+For deriving the DG
+discretization we start with a variational, mesh-dependent
+formulation of the problem,
+@f[
+ \sum_\kappa\left\{-\beta u,\nabla v)_\kappa+(u^+ \beta\cdot{\bf n}, v)_{\partial\kappa}\right\}=(f,v)_\Omega,
+@f]
+
+that originates from <a
+href="#step-12.transport-equation">[transport-equation]</a> 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}$ denote the
+<i>L<sup>2</sup></i>-inner products on the cell $\kappa$ and the
+boundary $\partial\kappa$ of the cell, respectively. $u^+$ is the
+value of <i>u</i> taken from the upwind cell with respect to $\beta$
+of the face, that is, the cell $\beta$ points away from. 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$.
+
+Hence, the discontinuous Galerkin
+scheme for the <a href="#step-12.transport-equation">[transport-equation]</a> is given
+by: find $u_h\in V_h$ such that for all $v_h\in V_h$ following
+equation holds:
+<a name="step-12.dg-transport1">@f[
+ \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]}@f]</a>
+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 <em>first version</em> of the DG method. We note that after a
+second integration by parts, we obtain: find $u_h\in V_h$ such that
+@f[
+ \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,
+@f]
+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 <a href="#step-12.dg-general1">[dg-general1]</a> and <a href="#step-12.upwind-flux">[upwind-flux]</a>,
+or <a href="#step-12.dg-transport1">[dg-transport1]</a>, respectively.
+
+Finally, we rewrite <a href="#step-12.dg-general1">[dg-general1]</a> 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, obtaining
+
+@f[
+ -\sum_\kappa(u_h,{\mathbf \beta}\cdot\nabla v_h)_\kappa
+ +\sum_{E\in\mathbb E_h^i} (u_h^-, \beta\cdot[v_h\mathbf n])_{E}
+ =(f,v_h)_\Omega-(g, \beta\cdot\mathbf n v_h)_{\Gamma_-}.
+@f]
+
+In this form, we need to implement a
/* further information on this license. */
// The first few files have already
- // been covered in step-12
+ // 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/compressed_sparsity_pattern.h>
#include <lac/sparse_matrix.h>
#include <grid/tria.h>
#include <grid/grid_generator.h>
#include <dofs/dof_tools.h>
#include <numerics/data_out.h>
#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
+ // <code>DoFHandler</code> and <code>FEValues</code>
+ // 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>
-#include <base/timer.h>
// Here come the new include files
// for using the MeshWorker framework:
#include <numerics/mesh_worker_info.h>
#include <numerics/mesh_worker_loop.h>
+ // Like in all programs, we finish
+ // this section by including the
+ // needed C++ headers and declaring
+ // we want to use objects in the
+ // dealii namespace without prefix.
#include <iostream>
#include <fstream>
- // The last step is as in all
- // previous programs:
using namespace dealii;
// @sect3{Equation data}
//
- // First, we need to describe the
- // coefficients in the equation. Here, this
- // concerns the boundary values, which we
- // choose in the same way as for step-12:
+ // First, we define a class
+ // describing the inhomogeneous
+ // boundary data. Since only its
+ // values are used, we implement
+ // value_list(), but leave all other
+ // functions of Function undefined.
template <int dim>
class BoundaryValues: public Function<dim>
{
values[i]=0.;
}
}
-
-
- // @sect3{Integrating cell and face matrices}
- // @sect3{Class: DGMethod}
+ // @sect3{Class: Step12}
//
- // After these preparations, we
- // proceed with the main part of this
- // program. The main class, here
- // called <code>DGMethod</code> is basically
- // the main class of step-6. One of
- // the differences is that there's no
- // ConstraintMatrix object. This is,
+ // After this preparations, we
+ // proceed with the main class of
+ // this program,
+ // called Step12. It is basically
+ // the main class of step-6. We do
+ // not have a ConstraintMatrix,
// because there are no hanging node
// constraints in DG discretizations.
+
+ // Major differences will only come
+ // up in the implementation of the
+ // assemble functions, since here, we
+ // not only need to cover the flux
+ // integrals over faces, we also use
+ // the MeshWorker interface to
+ // simplify the loops involved.
template <int dim>
-class DGMethod
+class Step12
{
public:
- DGMethod ();
- ~DGMethod ();
-
+ Step12 ();
void run ();
private:
FE_DGQ<dim> fe;
DoFHandler<dim> dof_handler;
+ // The next four members
+ // represent the linear system to
+ // be solved. #system_matrix and
+ // #right_hand_side are generated
+ // by assemble_system(), the
+ // #solution is computed in
+ // solve(). The #sparsity_pattern
+ // is used to determine the
+ // location of nonzero elements
+ // in #system_matrix.
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
-
- // In step-12 we had two solution vectors
- // that stored the solutions to the
- // problems corresponding to the two
- // different assembling routines
- // <code>assemble_system1</code> and
- // <code>assemble_system2</code>. In this
- // program, the goal is only to show the
- // MeshWorker framework, so we only
- // assemble the system in one of the two
- // ways, and consequently we have only
- // one solution vector along with the
- // single <code>assemble_system</code>
- // function declared above:
+
Vector<double> solution;
Vector<double> right_hand_side;
// then work on intermediate
// objects for which first, we
// here define typedefs to the
- // two info objects handed to the
+ // info objects handed to the
// local integration functions in
// order to make our life easier
// below.
+ typedef MeshWorker::DoFInfo<dim> DoFInfo;
typedef typename MeshWorker::IntegrationWorker<dim>::CellInfo CellInfo;
typedef typename MeshWorker::IntegrationWorker<dim>::FaceInfo FaceInfo;
// types of arguments, but have
// in fact other arguments
// already bound.
- static void integrate_cell_term (MeshWorker::DoFInfo<dim>& dinfo, CellInfo& info);
- static void integrate_boundary_term (MeshWorker::DoFInfo<dim>& dinfo, FaceInfo& info);
- static void integrate_face_term (MeshWorker::DoFInfo<dim>& dinfo1,
- MeshWorker::DoFInfo<dim>& dinfo2, FaceInfo& info1,
- FaceInfo& info2);
+ static void integrate_cell_term (DoFInfo& dinfo, CellInfo& info);
+ static void integrate_boundary_term (DoFInfo& dinfo, FaceInfo& info);
+ static void integrate_face_term (DoFInfo& dinfo1, DoFInfo& dinfo2,
+ FaceInfo& info1, FaceInfo& info2);
};
- // We start with the
- // constructor. This is the
- // place to change the
- // polynomial degree of the
- // finite element shape
- // functions.
+ // We start with the constructor. The
+ // 1 in the constructor call of #fe
+ // is the polynomial degree.
template <int dim>
-DGMethod<dim>::DGMethod ()
+Step12<dim>::Step12 ()
:
fe (1),
dof_handler (triangulation)
template <int dim>
-DGMethod<dim>::~DGMethod ()
+void Step12<dim>::setup_system ()
{
- dof_handler.clear ();
-}
-
-
// In the function that sets up the usual
// finite element data structures, we first
// need to distribute the DoFs.
-template <int dim>
-void DGMethod<dim>::setup_system ()
-{
dof_handler.distribute_dofs (fe);
- // The DoFs of a cell are coupled with all
- // DoFs of all neighboring cells, along
- // with all of its siblings on the current
- // cell. 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);
-
+ // We start by generating the
+ // sparsity pattern. To this end,
+ // we first fill an intermediate
+ // object of type
+ // CompressedSparsityPattern with
+ // the couplings appearing in the
+ // system. After building the
+ // pattern, this object is copied
+ // to #sparsity_pattern and can be
+ // discarded.
+
// To build the sparsity pattern for DG
// discretizations, we can call the
// function analogue to
// DoFTools::make_sparsity_pattern, which
// is called
// DoFTools::make_flux_sparsity_pattern:
- DoFTools::make_flux_sparsity_pattern (dof_handler, sparsity_pattern);
-
- // All following function calls are
- // already known.
- sparsity_pattern.compress();
+ CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
+ DoFTools::make_flux_sparsity_pattern (dof_handler, c_sparsity);
+ sparsity_pattern.copy_from(c_sparsity);
+ // Finally, we set up the structure
+ // of all components of the linear system.
system_matrix.reinit (sparsity_pattern);
-
solution.reinit (dof_handler.n_dofs());
right_hand_side.reinit (dof_handler.n_dofs());
}
// loops over cells and faces, we leave all
// this to the MeshWorker framework. In order
// to do so, we just have to define local
- // integration objects and use one of the
+ // integration functions and use one of the
// classes in namespace MeshWorker::Assembler
// to build the global system.
template <int dim>
-void DGMethod<dim>::assemble_system ()
+void Step12<dim>::assemble_system ()
{
// This is the magic object, which
// knows everything about the data
// object distributes these into
// the global sparse matrix and the
// right hand side vector.
- //
- // MeshWorker::AssemblingIntegrator
- // is not all that clever by
- // itself, but its capabilities are
- // provided the arguments provided
- // to the constructor and by its
- // second template argument. By
- // exchanging
- // MeshWorker::Assembler::SystemSimple,
- // we could for instance assemble a
- // BlockMatrix or just a Vector
- // instead.
- //
- // As noted in the discussion when
- // declaring the local integration
- // functions in the class
- // declaration, the arguments
- // expected by the assembling
- // integrator class are not
- // actually function
- // pointers. Rather, they are
- // objects that can be called like
- // functions with a certain number
- // of arguments. Consequently, we
- // could also pass objects with
- // appropriate operator()
- // implementations here, or the
- // result of std::bind if the local
- // integrators were, for example,
- // non-static member functions.
MeshWorker::IntegrationWorker<dim> integration_worker;
// First, we initialize the
// (determined by the first
// argument, which is an active
// iterator).
+ //
+ // As noted in the discussion when
+ // declaring the local integration
+ // functions in the class
+ // declaration, the arguments
+ // expected by the assembling
+ // integrator class are not
+ // actually function
+ // pointers. Rather, they are
+ // objects that can be called like
+ // functions with a certain number
+ // of arguments. Consequently, we
+ // could also pass objects with
+ // appropriate operator()
+ // implementations here, or the
+ // result of std::bind if the local
+ // integrators were, for example,
+ // non-static member functions.
MeshWorker::integration_loop<CellInfo, FaceInfo, dim>
(dof_handler.begin_active(), dof_handler.end(),
info_box,
- &DGMethod<dim>::integrate_cell_term,
- &DGMethod<dim>::integrate_boundary_term,
- &DGMethod<dim>::integrate_face_term,
+ &Step12<dim>::integrate_cell_term,
+ &Step12<dim>::integrate_boundary_term,
+ &Step12<dim>::integrate_face_term,
assembler, true);
}
// added soon).
template <int dim>
-void DGMethod<dim>::integrate_cell_term (MeshWorker::DoFInfo<dim>& dinfo, CellInfo& info)
+void Step12<dim>::integrate_cell_term (DoFInfo& dinfo, CellInfo& info)
{
// First, let us retrieve some of
// the objects used here from
// FESubfaceValues, in order to get access to
// normal vectors.
template <int dim>
-void DGMethod<dim>::integrate_boundary_term (MeshWorker::DoFInfo<dim>& dinfo, FaceInfo& info)
+void Step12<dim>::integrate_boundary_term (DoFInfo& dinfo, FaceInfo& info)
{
const FEFaceValuesBase<dim>& fe_v = info.fe_values();
FullMatrix<double>& local_matrix = dinfo.matrix(0).matrix;
// for each cell and two for coupling
// back and forth.
template <int dim>
-void DGMethod<dim>::integrate_face_term (MeshWorker::DoFInfo<dim>& dinfo1,
- MeshWorker::DoFInfo<dim>& dinfo2,
- FaceInfo& info1,
- FaceInfo& info2)
+void Step12<dim>::integrate_face_term (DoFInfo& dinfo1, DoFInfo& dinfo2,
+ FaceInfo& info1, FaceInfo& info2)
{
// For quadrature points, weights,
// etc., we use the
// PreconditionBlockSOR class with
// relaxation=1) does a much better job.
template <int dim>
-void DGMethod<dim>::solve (Vector<double> &solution)
+void Step12<dim>::solve (Vector<double> &solution)
{
SolverControl solver_control (1000, 1e-12, false, false);
SolverRichardson<> solver (solver_control);
// (to be more precise, in $H^1_\beta$)
// only.
template <int dim>
-void DGMethod<dim>::refine_grid ()
+void Step12<dim>::refine_grid ()
{
// The <code>DerivativeApproximation</code>
// class computes the gradients to
// in previous examples and will not
// be further commented on.
template <int dim>
-void DGMethod<dim>::output_results (const unsigned int cycle) const
+void Step12<dim>::output_results (const unsigned int cycle) const
{
// Write the grid in eps format.
std::string filename = "grid-";
// The following <code>run</code> function is
// similar to previous examples.
template <int dim>
-void DGMethod<dim>::run ()
+void Step12<dim>::run ()
{
for (unsigned int cycle=0; cycle<6; ++cycle)
{
{
try
{
- DGMethod<2> dgmethod;
+ Step12<2> dgmethod;
dgmethod.run ();
}
catch (std::exception &exc)