]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
add example 38, based on 12
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 18 Sep 2009 18:38:23 +0000 (18:38 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 18 Sep 2009 18:38:23 +0000 (18:38 +0000)
git-svn-id: https://svn.dealii.org/trunk@19475 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-38/Makefile [new file with mode: 0644]
deal.II/examples/step-38/doc/intro.dox [new file with mode: 0644]
deal.II/examples/step-38/doc/results.dox [new file with mode: 0644]
deal.II/examples/step-38/doc/step-12.3d-solution.png [new file with mode: 0644]
deal.II/examples/step-38/doc/step-12.grid-5.png [new file with mode: 0644]
deal.II/examples/step-38/doc/step-12.sol-0.png [new file with mode: 0644]
deal.II/examples/step-38/doc/step-12.sol-2.png [new file with mode: 0644]
deal.II/examples/step-38/doc/step-12.sol-5-3d.png [new file with mode: 0644]
deal.II/examples/step-38/doc/step-12.sol-5.png [new file with mode: 0644]
deal.II/examples/step-38/step-38.cc [new file with mode: 0644]

diff --git a/deal.II/examples/step-38/Makefile b/deal.II/examples/step-38/Makefile
new file mode 100644 (file)
index 0000000..db1de21
--- /dev/null
@@ -0,0 +1,154 @@
+# $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-38/doc/intro.dox b/deal.II/examples/step-38/doc/intro.dox
new file mode 100644 (file)
index 0000000..2f9248c
--- /dev/null
@@ -0,0 +1,13 @@
+<a name="Intro"></a>
+<h1>This will be an example of the new MeshWorker framework and will be available soon</h1>
+
+<h3>Overview</h3>
+
+This example is devoted to the <em>discontinuous Galerkin method</em>,
+or in short: DG method. It solves the same problem as step-12 (see
+there for a description of the problem and discretization), but here
+we use the MeshWorker framework in order to save reprogramming the
+cell/face loops.
+
+<h3>Implementation</h3>
+
diff --git a/deal.II/examples/step-38/doc/results.dox b/deal.II/examples/step-38/doc/results.dox
new file mode 100644 (file)
index 0000000..ecf8b2c
--- /dev/null
@@ -0,0 +1,91 @@
+<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-38/doc/step-12.3d-solution.png b/deal.II/examples/step-38/doc/step-12.3d-solution.png
new file mode 100644 (file)
index 0000000..5520d74
Binary files /dev/null and b/deal.II/examples/step-38/doc/step-12.3d-solution.png differ
diff --git a/deal.II/examples/step-38/doc/step-12.grid-5.png b/deal.II/examples/step-38/doc/step-12.grid-5.png
new file mode 100644 (file)
index 0000000..a2383b9
Binary files /dev/null and b/deal.II/examples/step-38/doc/step-12.grid-5.png differ
diff --git a/deal.II/examples/step-38/doc/step-12.sol-0.png b/deal.II/examples/step-38/doc/step-12.sol-0.png
new file mode 100644 (file)
index 0000000..16ae3fb
Binary files /dev/null and b/deal.II/examples/step-38/doc/step-12.sol-0.png differ
diff --git a/deal.II/examples/step-38/doc/step-12.sol-2.png b/deal.II/examples/step-38/doc/step-12.sol-2.png
new file mode 100644 (file)
index 0000000..081f3af
Binary files /dev/null and b/deal.II/examples/step-38/doc/step-12.sol-2.png differ
diff --git a/deal.II/examples/step-38/doc/step-12.sol-5-3d.png b/deal.II/examples/step-38/doc/step-12.sol-5-3d.png
new file mode 100644 (file)
index 0000000..6f80120
Binary files /dev/null and b/deal.II/examples/step-38/doc/step-12.sol-5-3d.png differ
diff --git a/deal.II/examples/step-38/doc/step-12.sol-5.png b/deal.II/examples/step-38/doc/step-12.sol-5.png
new file mode 100644 (file)
index 0000000..d9cfd87
Binary files /dev/null and b/deal.II/examples/step-38/doc/step-12.sol-5.png differ
diff --git a/deal.II/examples/step-38/step-38.cc b/deal.II/examples/step-38/step-38.cc
new file mode 100644 (file)
index 0000000..2e79616
--- /dev/null
@@ -0,0 +1,1650 @@
+/* $Id$ */
+/* Author: Guido Kanschat, Texas A&M University, 2009 */
+
+/*    $Id$       */
+/*                                                                */
+/*    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 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 example 12
+                                // 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>
+#include <fe/mapping_q1.h>
+#include <fe/fe_dgq.h>
+#include <lac/solver_richardson.h>
+#include <lac/precondition_block.h>
+#include <numerics/derivative_approximation.h>
+#include <base/timer.h>
+
+                                // Here come the new include files
+                                // for using the MeshWorker framework
+
+
+#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.