]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
step 12 removed
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 12 Mar 2010 20:54:11 +0000 (20:54 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 12 Mar 2010 20:54:11 +0000 (20:54 +0000)
git-svn-id: https://svn.dealii.org/trunk@20816 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-12/Makefile [deleted file]
deal.II/examples/step-12/doc/intro.dox [deleted file]
deal.II/examples/step-12/doc/results.dox [deleted file]
deal.II/examples/step-12/doc/step-12.3d-solution.png [deleted file]
deal.II/examples/step-12/doc/step-12.grid-5.png [deleted file]
deal.II/examples/step-12/doc/step-12.sol-0.png [deleted file]
deal.II/examples/step-12/doc/step-12.sol-2.png [deleted file]
deal.II/examples/step-12/doc/step-12.sol-5-3d.png [deleted file]
deal.II/examples/step-12/doc/step-12.sol-5.png [deleted file]
deal.II/examples/step-12/step-12.cc [deleted file]

diff --git a/deal.II/examples/step-12/Makefile b/deal.II/examples/step-12/Makefile
deleted file mode 100644 (file)
index db1de21..0000000
+++ /dev/null
@@ -1,154 +0,0 @@
-# $Id$
-
-
-# For the small projects Makefile, you basically need to fill in only
-# four fields.
-#
-# The first is the name of the application. It is assumed that the
-# application name is the same as the base file name of the single C++
-# file from which the application is generated.
-target = $(basename $(shell echo step-*.cc))
-
-# The second field determines whether you want to run your program in
-# debug or optimized mode. The latter is significantly faster, but no
-# run-time checking of parameters and internal states is performed, so
-# you should set this value to `on' while you develop your program,
-# and to `off' when running production computations.
-debug-mode = on
-
-
-# As third field, we need to give the path to the top-level deal.II
-# directory. You need to adjust this to your needs. Since this path is
-# probably the most often needed one in the Makefile internals, it is
-# designated by a single-character variable, since that can be
-# reference using $D only, i.e. without the parentheses that are
-# required for most other parameters, as e.g. in $(target).
-D = ../../
-
-
-# The last field specifies the names of data and other files that
-# shall be deleted when calling `make clean'. Object and backup files,
-# executables and the like are removed anyway. Here, we give a list of
-# files in the various output formats that deal.II supports.
-clean-up-files = *gmv *gnuplot *gpl *eps *pov
-
-
-
-
-#
-#
-# Usually, you will not need to change anything beyond this point.
-#
-#
-# The next statement tell the `make' program where to find the
-# deal.II top level directory and to include the file with the global
-# settings
-include $D/common/Make.global_options
-
-
-# Since the whole project consists of only one file, we need not
-# consider difficult dependencies. We only have to declare the
-# libraries which we want to link to the object file, and there need
-# to be two sets of libraries: one for the debug mode version of the
-# application and one for the optimized mode. Here we have selected
-# the versions for 2d. Note that the order in which the libraries are
-# given here is important and that your applications won't link
-# properly if they are given in another order.
-#
-# You may need to augment the lists of libraries when compiling your
-# program for other dimensions, or when using third party libraries
-libs.g   = $(lib-deal2-2d.g) \
-          $(lib-lac.g)      \
-           $(lib-base.g)
-libs.o   = $(lib-deal2-2d.o) \
-          $(lib-lac.o)      \
-           $(lib-base.o)
-
-
-# We now use the variable defined above which switch between debug and
-# optimized mode to select the set of libraries to link with. Included
-# in the list of libraries is the name of the object file which we
-# will produce from the single C++ file. Note that by default we use
-# the extension .g.o for object files compiled in debug mode and .o for
-# object files in optimized mode (or whatever the local default on your
-# system is instead of .o).
-ifeq ($(debug-mode),on)
-  libraries = $(target).g.$(OBJEXT) $(libs.g)
-else
-  libraries = $(target).$(OBJEXT) $(libs.o)
-endif
-
-
-# Now comes the first production rule: how to link the single object
-# file produced from the single C++ file into the executable. Since
-# this is the first rule in the Makefile, it is the one `make' selects
-# if you call it without arguments.
-$(target) : $(libraries)
-       @echo ============================ Linking $@
-       @$(CXX) -o $@$(EXEEXT) $^ $(LIBS) $(LDFLAGS)
-
-
-# To make running the application somewhat independent of the actual
-# program name, we usually declare a rule `run' which simply runs the
-# program. You can then run it by typing `make run'. This is also
-# useful if you want to call the executable with arguments which do
-# not change frequently. You may then want to add them to the
-# following rule:
-run: $(target)
-       @echo ============================ Running $<
-       @./$(target)$(EXEEXT)
-
-
-# As a last rule to the `make' program, we define what to do when
-# cleaning up a directory. This usually involves deleting object files
-# and other automatically created files such as the executable itself,
-# backup files, and data files. Since the latter are not usually quite
-# diverse, you needed to declare them at the top of this file.
-clean:
-       -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files)
-
-
-# Since we have not yet stated how to make an object file from a C++
-# file, we should do so now. Since the many flags passed to the
-# compiler are usually not of much interest, we suppress the actual
-# command line using the `at' sign in the first column of the rules
-# and write the string indicating what we do instead.
-./%.g.$(OBJEXT) :
-       @echo ==============debug========= $(<F)
-       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
-./%.$(OBJEXT) :
-       @echo ==============optimized===== $(<F)
-       @$(CXX) $(CXXFLAGS.o) -c $< -o $@
-
-
-# The following statement tells make that the rules `run' and `clean'
-# are not expected to produce files of the same name as Makefile rules
-# usually do.
-.PHONY: run clean
-
-
-# Finally there is a rule which you normally need not care much about:
-# since the executable depends on some include files from the library,
-# besides the C++ application file of course, it is necessary to
-# re-generate the executable when one of the files it depends on has
-# changed. The following rule to created a dependency file
-# `Makefile.dep', which `make' uses to determine when to regenerate
-# the executable. This file is automagically remade whenever needed,
-# i.e. whenever one of the cc-/h-files changed. Make detects whether
-# to remake this file upon inclusion at the bottom of this file.
-#
-# If the creation of Makefile.dep fails, blow it away and fail
-Makefile.dep: $(target).cc Makefile \
-              $(shell echo $D/*/include/*/*.h)
-       @echo ============================ Remaking $@
-       @$D/common/scripts/make_dependencies  $(INCLUDE) -B. $(target).cc \
-               > $@ \
-         || (rm -f $@ ; false)
-       @if test -s $@ ; then : else rm $@ ; fi
-
-
-# To make the dependencies known to `make', we finally have to include
-# them:
-include Makefile.dep
-
-
diff --git a/deal.II/examples/step-12/doc/intro.dox b/deal.II/examples/step-12/doc/intro.dox
deleted file mode 100644 (file)
index d96538c..0000000
+++ /dev/null
@@ -1,266 +0,0 @@
-<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.
-
-
-
-
-
-
-
diff --git a/deal.II/examples/step-12/doc/results.dox b/deal.II/examples/step-12/doc/results.dox
deleted file mode 100644 (file)
index ecf8b2c..0000000
+++ /dev/null
@@ -1,91 +0,0 @@
-<h1>Results</h1>
-
-
-The output of this program consist of the console output, the eps
-files including the grids, and the solutions given in gnuplot format.
-@code
-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>...
-@endcode
-
-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.
-
-@image html step-12.sol-0.png
-@image html step-12.sol-2.png
-@image html 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:
-
-@image html step-12.grid-5.png
-@image html step-12.3d-solution.png
-
-
-And finally we show a plot of a 3d computation.
-
-@image html step-12.sol-5-3d.png
-
diff --git a/deal.II/examples/step-12/doc/step-12.3d-solution.png b/deal.II/examples/step-12/doc/step-12.3d-solution.png
deleted file mode 100644 (file)
index 5520d74..0000000
Binary files a/deal.II/examples/step-12/doc/step-12.3d-solution.png and /dev/null differ
diff --git a/deal.II/examples/step-12/doc/step-12.grid-5.png b/deal.II/examples/step-12/doc/step-12.grid-5.png
deleted file mode 100644 (file)
index a2383b9..0000000
Binary files a/deal.II/examples/step-12/doc/step-12.grid-5.png and /dev/null differ
diff --git a/deal.II/examples/step-12/doc/step-12.sol-0.png b/deal.II/examples/step-12/doc/step-12.sol-0.png
deleted file mode 100644 (file)
index 16ae3fb..0000000
Binary files a/deal.II/examples/step-12/doc/step-12.sol-0.png and /dev/null differ
diff --git a/deal.II/examples/step-12/doc/step-12.sol-2.png b/deal.II/examples/step-12/doc/step-12.sol-2.png
deleted file mode 100644 (file)
index 081f3af..0000000
Binary files a/deal.II/examples/step-12/doc/step-12.sol-2.png and /dev/null differ
diff --git a/deal.II/examples/step-12/doc/step-12.sol-5-3d.png b/deal.II/examples/step-12/doc/step-12.sol-5-3d.png
deleted file mode 100644 (file)
index 6f80120..0000000
Binary files a/deal.II/examples/step-12/doc/step-12.sol-5-3d.png and /dev/null differ
diff --git a/deal.II/examples/step-12/doc/step-12.sol-5.png b/deal.II/examples/step-12/doc/step-12.sol-5.png
deleted file mode 100644 (file)
index d9cfd87..0000000
Binary files a/deal.II/examples/step-12/doc/step-12.sol-5.png and /dev/null differ
diff --git a/deal.II/examples/step-12/step-12.cc b/deal.II/examples/step-12/step-12.cc
deleted file mode 100644 (file)
index 4d8f246..0000000
+++ /dev/null
@@ -1,1682 +0,0 @@
-/* $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;
-}
-
-

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.