+++ /dev/null
-<a name="Intro"></a>
-<h1>Introduction</h1>
-
-
-<h3>Overview</h3>
-
-This example is devoted to the <em>discontinuous 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> Two different assembling routines for the system matrix based on
- face terms given as a sum of integrals that
- <ol>
- <li> loops over all cell and all their faces, or that
- <li> loops over all faces, whereas each face is treated only once.
- </ol>
- <li> Time comparison of the two assembling routines.
-</ol>
-
-
-<h3>Problem</h3>
-
-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
-<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>
-
-Following the general paradigm of deriving DG discretizations for
-purely hyperbolic equations, we first consider the general hyperbolic
-problem
-@f[
- \nabla\cdot {\mathcal F}(u)=f \qquad\mbox{in }\Omega,
-@f]
-subject to appropriate boundary conditions. Here ${\mathcal F}$
-denotes the flux function of the equation under consideration that in
-our case, see equation <a href="#step-12.transport-equation">[transport-equation]</a>, 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,
-@f[
- \sum_\kappa\left\{-({\mathcal F}(u),\nabla v)_\kappa+({\mathcal
- F}(u)\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}$ 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
-<em>numerical flux</em> 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.
-@f[
-{\mathcal H}(u,u,{\bf n})={\mathcal F}(u)\cdot{\bf n},
-@f]
-and conservative, i.e.
-<a name="step-12.conservativity">@f[
-{\mathcal H}(v,w,{\bf n})=-{\mathcal H}(w,v,-{\bf n}).
-\qquad\qquad\qquad\mathrm{[conservativity]}@f]</a>
-This yields the following <em>discontinuous Galerkin
- discretization</em>: find $u_h\in V_h$ such that
-<a name="step-12.dg-general1">@f[
- \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]}@f]</a>
-%Boundary conditions are realized by replacing $u_h^-$ on the inflow boundary $\Gamma_-$ by the boundary function $g$.
-In the special case of the
-<a href="#step-12.transport-equation">[transport-equation]</a> the numerical flux in its simplest form
-is given by
-<a name="step-12.upwind-flux">@f[
- {\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]}@f]</a>
-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 <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:
-Find $u_h\in V_h$ such that
-<a name="step-12.dg-general2">@f[
- -\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]}@f]</a>
-By employing <a href="#step-12.conservativity">[conservativity]</a> of the numerical flux
-this equation simplifies to: find $u_h\in V_h$ such that
-<a name="step-12.dg-general3">@f[
- -\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]}@f]</a>
-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 <a href="#step-12.dg-general2">[dg-general2]</a> that finally produces the
-minus sign in the jump $[v_h]$ in equation <a href="#step-12.dg-general3">[dg-general3]</a>.
-
-For the linear <a href="#step-12.transport-equation">[transport-equation]</a>
-equation <a href="#step-12.dg-general3">[dg-general3]</a> simplifies to
-<a name="step-12.dg-transport2">@f[
- -\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]}@f]</a>
-which will be refered to as <em>second version</em> of the DG method.
-
-
-<h3>Implementation</h3>
-
-
-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 <a href="#step-12.dg-transport1">[dg-transport1]</a>
-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 <a href="#step-12.dg-transport2">[dg-transport2]</a> that
-includes a sum of integrals over all faces $e$. Here, several
-difficulties occurs.
-<ol>
-<li> 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 <code>deal.II</code>
- faces are accessible from cells but not visa versa.
-<li> 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:
- <ol>
- <li> 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.
- <li> If the two cells are of different refinement levels we access
- and treat the face from the coarser cell.
- </ol>
-</ol>
-Before we start with the description of the code we first introduce
-its main ingredients. The main class is called
-<code>DGMethod</code>. 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 <code>assemble_system1</code> and <code>assemble_system2</code>
-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
-<ol>
-<li> face is at boundary,
-<li> neighboring cell is finer,
-<li> neighboring cell is of the same refinement level, and
-<li> neighboring cell is coarser
-</ol>
-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
-<code>DGTransportEquation</code>. This class includes objects of all
-equation--specific functions, the <code>RHS</code> and the
-<code>BoundaryValues</code> class, both derived from the <code>Function</code>
-class, and the <code>Beta</code> class representing the vector field.
-Furthermore, the <code>DGTransportEquation</code> class comprises member
-functions <code>assemble_face_terms1</code> and
-<code>assemble_face_terms2</code> that are invoked by the
-<code>assemble_system1</code> and <code>assemble_system2</code> functions of the
-<code>DGMethod</code>, respectively, and the functions
-<code>assemble_cell_term</code> and <code>assemble_boundary_term</code> that
-are the same for both assembling routines. Due to the encapsulation of
-all equation- and scheme-specific functions, the
-<code>DGTransportEquation</code> class can easily be replaced by a similar
-class that implements a different equation and a different DG method.
-Indeed, the implementation of the <code>assemble_system1</code> and
-<code>assemble_system2</code> functions of the <code>DGMethod</code> 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
-<code>DGMethod<2></code> by <code>DGMethod<3></code> the code runs in 3d.
-
-
-
-
-
-
-
+++ /dev/null
-/* $Id$ */
-/* Author: Ralf Hartmann, University of Heidelberg, 2001 */
-
-/* $Id$ */
-/* */
-/* 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. */
-
- // 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 <code>MappingQ1</code> class
- // that gives the standard bilinear
- // mapping. For bilinear mappings use
- // an object of this class rather
- // than an object of the
- // <code>MappingQ(1)</code> class, as the
- // <code>MappingQ1</code> 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
- // <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>
- // Finally we do some time comparison
- // using the <code>Timer</code> 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;
-
- // @sect3{Equation data}
- //
- // First we define the classes
- // representing the equation-specific
- // functions. Both classes, <code>RHS</code>
- // and <code>BoundaryValues</code>, are
- // derived from the <code>Function</code>
- // class. Only the <code>value_list</code>
- // 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 <code>Beta</code> represents the
- // vector valued flow field of the
- // linear transport equation and is
- // not derived from the <code>Function</code>
- // class as we prefer to get function
- // values of type <code>Point</code> rather
- // than of type
- // <code>Vector@<double@></code>. This, because
- // there exist scalar products
- // between <code>Point</code> and <code>Point</code> as
- // well as between <code>Point</code> and
- // <code>Tensor</code>, 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
- // <code>value_list</code> 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.;
- }
-}
-
-
- // @sect3{Class: DGTransportEquation}
- //
- // Next we define the
- // equation-dependent and
- // DG-method-dependent class
- // <code>DGTransportEquation</code>. Its
- // member functions were already
- // mentioned in the Introduction and
- // will be explained
- // below. Furthermore it includes
- // objects of the previously defined
- // <code>Beta</code>, <code>RHS</code> and
- // <code>BoundaryValues</code> 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 ()
-{}
-
-
- // @sect4{Function: assemble_cell_term}
- //
- // The <code>assemble_cell_term</code>
- // function assembles the cell terms
- // of the discretization.
- // <code>ui_vi_matrix</code> is a cell matrix,
- // i.e. for a DG method of degree 1,
- // it is of size 4 times 4, and
- // <code>cell_vector</code> is of size 4.
- // When this function is invoked,
- // <code>fe_v</code> 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 <code>fe_v</code> for the
- // quadrature weights,
- const std::vector<double> &JxW = fe_v.get_JxW_values ();
-
- // Then the flow field beta and the
- // <code>rhs_function</code> 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];
- }
-}
-
-
- // @sect4{Function: assemble_boundary_term}
- //
- // The <code>assemble_boundary_term</code>
- // function assembles the face terms
- // at boundary faces. When this
- // function is invoked, <code>fe_v</code> 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
- // <code>FEValues</code> 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];
- }
-}
-
-
- // @sect4{Function: assemble_face_term1}
- //
- // The <code>assemble_face_term1</code>
- // 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,
- // <code>fe_v</code> and <code>fe_v_neighbor</code> 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
- // <code>ui_vi_matrix</code> this function
- // gets a new argument
- // <code>ue_vi_matrix</code>, that stores
- // contributions to the system matrix
- // that are based on exterior values
- // of $u$ and interior values of
- // $v$. Here we note that <code>ue</code> is
- // the short notation for <code>u
- // exterior</code> 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];
- }
-}
-
-
- // @sect4{Function: assemble_face_term2}
- //
- // Now we look at the
- // <code>assemble_face_term2</code> 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
- // <code>ui_ve_matrix</code> and
- // <code>ue_ve_matrix</code> 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];
- }
- }
-}
-
-
- // @sect3{Class: DGMethod}
- //
- // 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,
- // 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
- // <code>assemble_system1</code> and
- // <code>assemble_system2</code>;
- 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());
-}
-
-
- // @sect4{Function: assemble_system1}
- //
- // We proceed with the
- // <code>assemble_system1</code> function that
- // implements the DG discretization
- // in its first version. This
- // function repeatedly calls the
- // <code>assemble_cell_term</code>,
- // <code>assemble_boundary_term</code> and
- // <code>assemble_face_term1</code> functions
- // of the <code>DGTransportEquation</code>
- // object. The
- // <code>assemble_boundary_term</code> covers
- // the first case mentioned in the
- // introduction.
- //
- // 1. face is at boundary
- //
- // This function takes a
- // <code>FEFaceValues</code> object as
- // argument. In contrast to that
- // <code>assemble_face_term1</code>
- // takes two <code>FEFaceValuesBase</code>
- // 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 <code>FEFaceValues</code>
- // or of class <code>FESubfaceValues</code>
- // (both derived from
- // <code>FEFaceValuesBase</code>) according to
- // the remaining cases mentioned
- // in the introduction:
- //
- // 2. neighboring cell is finer
- // (current cell: <code>FESubfaceValues</code>,
- // neighboring cell: <code>FEFaceValues</code>);
- //
- // 3. neighboring cell is of the same
- // refinement level (both, current
- // and neighboring cell:
- // <code>FEFaceValues</code>);
- //
- // 4. neighboring cell is coarser
- // (current cell: <code>FEFaceValues</code>,
- // neighboring cell:
- // <code>FESubfaceValues</code>).
- //
- // 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
- // <code>update_flags</code> for the
- // <code>FEValues</code> and the
- // <code>FEFaceValues</code> 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 <code>FEValues</code>
- // object. Note, that since version
- // 3.2.0 of deal.II the constructor
- // of this class takes a
- // <code>Mapping</code> object as first
- // argument. Although the
- // constructor without <code>Mapping</code>
- // 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 <code>MappingQ1</code> mapping)
- // and makes it easier to change
- // the mapping object later.
- FEValues<dim> fe_v (
- mapping, fe, quadrature, update_flags);
-
- // Similarly we create the
- // <code>FEFaceValues</code> and
- // <code>FESubfaceValues</code> 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 <code>vi</code> (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
- // <code>assemble_face_term1</code>
- // function contributions to
- // the cell matrices and the
- // cell vector are only
- // ADDED. Therefore on each
- // cell we need to reset the
- // <code>ui_vi_matrix</code> and
- // <code>cell_vector</code> to zero,
- // before assembling the cell terms.
- ui_vi_matrix = 0;
- cell_vector = 0;
-
- // Now we reinit the <code>FEValues</code>
- // object for the current cell
- fe_v.reinit (cell);
-
- // and call the function
- // that assembles the cell
- // terms. The first argument is
- // the <code>FEValues</code> 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
- // <code>ue_vi_matrix</code> 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
- // <code>FEFaceValues</code>
- // 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 <code>dim==1</code>.
- 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
- // <code>cell</code>.
- 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
- // <code>neighbor_child</code>
- // 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
- // <code>ue_vi_matrix</code>
- // on each
- // subface
- // because on
- // each subface
- // the <code>un</code>
- // belong to
- // different
- // neighboring
- // cells.
- ue_vi_matrix = 0;
-
- // As already
- // mentioned
- // above for
- // the current
- // case (case
- // 2) we employ
- // the
- // <code>FESubfaceValues</code>
- // 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
- // <code>ue_vi_matrix</code>
- // 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 <code>if
- // (face-@>has_children())</code>
- }
- 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
- // <code>FEFaceValues</code>
- // 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 <code>if
- // (neighbor-@>level()
- // ==
- // cell-@>level())</code>
- }
- 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
- // <code>FEFaceValues</code>
- // 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
- // <code>neighbor_child</code>
- // cell,
- neighbor->get_dof_indices (dofs_neighbor);
-
- // and distribute the
- // <code>ue_vi_matrix</code>.
- 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 <code>face not at boundary</code>:
- }
- // End of loop over all faces:
- }
-
- // Finally we distribute the
- // <code>ui_vi_matrix</code>
- 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);
- }
-}
-
-
- // @sect4{Function: assemble_system2}
- //
- // We proceed with the
- // <code>assemble_system2</code> function that
- // implements the DG discretization
- // in its second version. This
- // function is very similar to the
- // <code>assemble_system1</code>
- // function. Therefore, here we only
- // discuss the differences between
- // the two functions. This function
- // repeatedly calls the
- // <code>assemble_face_term2</code> 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
- // <code>fe_v_face_neighbor</code> 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 <code>ve</code> (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);
- }
-}
-
-
- // @sect3{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
- // <code>DerivativeApproximation</code> class
- // that computes the approximate
- // gradients in a way similar to the
- // <code>GradientEstimation</code> described
- // in step-9 of this tutorial. In
- // fact, the
- // <code>DerivativeApproximation</code> class
- // was developed following the
- // <code>GradientEstimation</code> 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 <code>DerivativeApproximation</code>
- // 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 <code>run</code> function is
- // similar to previous examples. The
- // only difference is that the
- // problem is assembled and solved
- // twice on each refinement step;
- // first by <code>assemble_system1</code> that
- // implements the first version and
- // then by <code>assemble_system2</code> 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 <code>main</code> 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;
-}
-
-
projects / dealii.git / commitdiff