From: Guido Kanschat Date: Fri, 12 Mar 2010 20:54:11 +0000 (+0000) Subject: step 12 removed X-Git-Tag: v8.0.0~6332 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=76fe88d206bbffe97cce881cd13ff1a91a13c065;p=dealii.git step 12 removed git-svn-id: https://svn.dealii.org/trunk@20816 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-12/Makefile b/deal.II/examples/step-12/Makefile deleted file mode 100644 index db1de21066..0000000000 --- a/deal.II/examples/step-12/Makefile +++ /dev/null @@ -1,154 +0,0 @@ -# $Id$ - - -# For the small projects Makefile, you basically need to fill in only -# four fields. -# -# The first is the name of the application. It is assumed that the -# application name is the same as the base file name of the single C++ -# file from which the application is generated. -target = $(basename $(shell echo step-*.cc)) - -# The second field determines whether you want to run your program in -# debug or optimized mode. The latter is significantly faster, but no -# run-time checking of parameters and internal states is performed, so -# you should set this value to `on' while you develop your program, -# and to `off' when running production computations. -debug-mode = on - - -# As third field, we need to give the path to the top-level deal.II -# directory. You need to adjust this to your needs. Since this path is -# probably the most often needed one in the Makefile internals, it is -# designated by a single-character variable, since that can be -# reference using $D only, i.e. without the parentheses that are -# required for most other parameters, as e.g. in $(target). -D = ../../ - - -# The last field specifies the names of data and other files that -# shall be deleted when calling `make clean'. Object and backup files, -# executables and the like are removed anyway. Here, we give a list of -# files in the various output formats that deal.II supports. -clean-up-files = *gmv *gnuplot *gpl *eps *pov - - - - -# -# -# Usually, you will not need to change anything beyond this point. -# -# -# The next statement tell the `make' program where to find the -# deal.II top level directory and to include the file with the global -# settings -include $D/common/Make.global_options - - -# Since the whole project consists of only one file, we need not -# consider difficult dependencies. We only have to declare the -# libraries which we want to link to the object file, and there need -# to be two sets of libraries: one for the debug mode version of the -# application and one for the optimized mode. Here we have selected -# the versions for 2d. Note that the order in which the libraries are -# given here is important and that your applications won't link -# properly if they are given in another order. -# -# You may need to augment the lists of libraries when compiling your -# program for other dimensions, or when using third party libraries -libs.g = $(lib-deal2-2d.g) \ - $(lib-lac.g) \ - $(lib-base.g) -libs.o = $(lib-deal2-2d.o) \ - $(lib-lac.o) \ - $(lib-base.o) - - -# We now use the variable defined above which switch between debug and -# optimized mode to select the set of libraries to link with. Included -# in the list of libraries is the name of the object file which we -# will produce from the single C++ file. Note that by default we use -# the extension .g.o for object files compiled in debug mode and .o for -# object files in optimized mode (or whatever the local default on your -# system is instead of .o). -ifeq ($(debug-mode),on) - libraries = $(target).g.$(OBJEXT) $(libs.g) -else - libraries = $(target).$(OBJEXT) $(libs.o) -endif - - -# Now comes the first production rule: how to link the single object -# file produced from the single C++ file into the executable. Since -# this is the first rule in the Makefile, it is the one `make' selects -# if you call it without arguments. -$(target) : $(libraries) - @echo ============================ Linking $@ - @$(CXX) -o $@$(EXEEXT) $^ $(LIBS) $(LDFLAGS) - - -# To make running the application somewhat independent of the actual -# program name, we usually declare a rule `run' which simply runs the -# program. You can then run it by typing `make run'. This is also -# useful if you want to call the executable with arguments which do -# not change frequently. You may then want to add them to the -# following rule: -run: $(target) - @echo ============================ Running $< - @./$(target)$(EXEEXT) - - -# As a last rule to the `make' program, we define what to do when -# cleaning up a directory. This usually involves deleting object files -# and other automatically created files such as the executable itself, -# backup files, and data files. Since the latter are not usually quite -# diverse, you needed to declare them at the top of this file. -clean: - -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files) - - -# Since we have not yet stated how to make an object file from a C++ -# file, we should do so now. Since the many flags passed to the -# compiler are usually not of much interest, we suppress the actual -# command line using the `at' sign in the first column of the rules -# and write the string indicating what we do instead. -./%.g.$(OBJEXT) : - @echo ==============debug========= $( $@ \ - || (rm -f $@ ; false) - @if test -s $@ ; then : else rm $@ ; fi - - -# To make the dependencies known to `make', we finally have to include -# them: -include Makefile.dep - - diff --git a/deal.II/examples/step-12/doc/intro.dox b/deal.II/examples/step-12/doc/intro.dox deleted file mode 100644 index d96538c885..0000000000 --- a/deal.II/examples/step-12/doc/intro.dox +++ /dev/null @@ -1,266 +0,0 @@ - -

Introduction

- - -

Overview

- -This example is devoted to the discontinuous Galerkin method, or -in short: DG method. It includes the following topics. -
    -
  1. Discretization of the linear transport equation with the DG method -
  2. Two different assembling routines for the system matrix based on - face terms given as a sum of integrals that -
      -
    1. loops over all cell and all their faces, or that -
    2. loops over all faces, whereas each face is treated only once. -
    -
  3. Time comparison of the two assembling routines. -
- - -

Problem

- -The DG method was first introduced to discretize simple transport -equations. Over the past years DG methods have been applied to a -variety of problems and many different schemes were introduced -employing a big zoo of different convective and diffusive fluxes. As -this example's purpose is to illustrate some implementational issues -of the DG discretization only, here we simply consider the linear -transport equation -@f[ - \nabla\cdot \left({\mathbf \beta} u\right)=f \qquad\mbox{in }\Omega, -\qquad\qquad\qquad\mathrm{[transport-equation]}@f] -subject to the boundary conditions -@f[ -u=g\quad\mbox{on }\Gamma_-, -@f] -on the inflow part $\Gamma_-$ of the boundary $\Gamma=\partial\Omega$ -of the domain. Here, ${\mathbf \beta}={\mathbf \beta}({\bf x})$ denotes a -vector field, $f$ a source function, $u$ the (scalar) solution -function, $g$ a boundary value function, -@f[ -\Gamma_-:=\{{\bf x}\in\Gamma, {\mathbf \beta}({\bf x})\cdot{\bf n}({\bf x})<0\} -@f] -the inflow part of the boundary of the domain and ${\bf n}$ denotes -the unit outward normal to the boundary $\Gamma$. Equation -[transport-equation] is the conservative version of the -transport equation already considered in step 9 of this tutorial. - -In particular, we consider problem [transport-equation] on -$\Omega=[0,1]^2$ with ${\mathbf \beta}=\frac{1}{|x|}(-x_2, x_1)$ -representing a circular counterclockwise flow field, $f=0$ and $g=1$ -on ${\bf x}\in\Gamma_-^1:=[0,0.5]\times\{0\}$ and $g=0$ on ${\bf x}\in -\Gamma_-\setminus \Gamma_-^1$. - - -

Discretization

- -Following the general paradigm of deriving DG discretizations for -purely hyperbolic equations, we first consider the general hyperbolic -problem -@f[ - \nabla\cdot {\mathcal F}(u)=f \qquad\mbox{in }\Omega, -@f] -subject to appropriate boundary conditions. Here ${\mathcal F}$ -denotes the flux function of the equation under consideration that in -our case, see equation [transport-equation], is represented by -${\mathcal F}(u)={\mathbf \beta} u$. For deriving the DG -discretization we start with a variational, mesh-dependent -formulation of the problem, -@f[ - \sum_\kappa\left\{-({\mathcal F}(u),\nabla v)_\kappa+({\mathcal - F}(u)\cdot{\bf n}, v)_{\partial\kappa}\right\}=(f,v)_\Omega, -@f] -that originates from [transport-equation] by multiplication with -a test function $v$ and integration by parts on each cell $\kappa$ of -the triangulation. Here $(\cdot, \cdot)_\kappa$ and $(\cdot, -\cdot)_{\partial\kappa}$ simply denote the integrals over the cell -$\kappa$ and the boundary $\partial\kappa$ of the cell, -respectively. To discretize the problem, the functions $u$ and $v$ are -replaced by discrete functions $u_h$ and $v_h$ that in the case of -discontinuous Galerkin methods belong to the space $V_h$ of -discontinuous piecewise polynomial functions of some degree $p$. Due -to the discontinuity of the discrete function $u_h$ on interelement -faces, the flux ${\mathcal F}(u)\cdot{\bf n}$ must be replaced by a -numerical flux function ${\mathcal H}(u_h^+, u_h^-, {\bf n})$, -where $u_h^+|_{\partial\kappa}$ denotes the inner trace (w.r.t. the -cell $\kappa$) of $u_h$ and $u_h^-|_{\partial\kappa}$ the outer trace, -i.e. the value of $u_h$ on the neighboring cell. Furthermore the -numerical flux function ${\mathcal H}$, among other things, must be -consistent, i.e. -@f[ -{\mathcal H}(u,u,{\bf n})={\mathcal F}(u)\cdot{\bf n}, -@f] -and conservative, i.e. -@f[ -{\mathcal H}(v,w,{\bf n})=-{\mathcal H}(w,v,-{\bf n}). -\qquad\qquad\qquad\mathrm{[conservativity]}@f] -This yields the following discontinuous Galerkin - discretization: find $u_h\in V_h$ such that -@f[ - \sum_\kappa\left\{-({\mathcal F}(u_h),\nabla v_h)_\kappa+({\mathcal H}(u_h^+,u_h^-,{\bf n}), v_h)_{\partial\kappa}\right\}=(f,v_h)_\Omega, \quad\forall v_h\in V_h. -\qquad\qquad\qquad\mathrm{[dg-general1]}@f] -%Boundary conditions are realized by replacing $u_h^-$ on the inflow boundary $\Gamma_-$ by the boundary function $g$. -In the special case of the -[transport-equation] the numerical flux in its simplest form -is given by -@f[ - {\mathcal H}(u_h^+,u_h^-,{\bf n})({\bf x})=\left\{\begin{array}{ll} - ({\mathbf \beta}\cdot{\bf n}\, u_h^-)({\bf x}),&\mbox{for } {\mathbf \beta}({\bf x})\cdot{\bf n}({\bf x})<0,\\ - ({\mathbf \beta}\cdot{\bf n}\, u_h^+)({\bf x}),&\mbox{for } {\mathbf \beta}({\bf x})\cdot{\bf n}({\bf x})\geq 0, -\end{array} -\right. -\qquad\qquad\qquad\mathrm{[upwind-flux]}@f] -where on the inflow part of the cell the value is taken from the -neighboring cell, $u_h^-$, and on the outflow part the value is -taken from the current cell, $u_h^+$. Hence, the discontinuous Galerkin -scheme for the [transport-equation] is given -by: find $u_h\in V_h$ such that for all $v_h\in V_h$ following -equation holds: -@f[ - \sum_\kappa\left\{-(u_h,{\mathbf \beta}\cdot\nabla v_h)_\kappa - +({\mathbf \beta}\cdot{\bf n}\, u_h, v_h)_{\partial\kappa_+} - +({\mathbf \beta}\cdot{\bf n}\, u_h^-, v_h)_{\partial\kappa_-\setminus\Gamma}\right\} - =(f,v_h)_\Omega-({\mathbf \beta}\cdot{\bf n}\, g, v_h)_{\Gamma_-}, -\qquad\qquad\qquad\mathrm{[dg-transport1]}@f] -where $\partial\kappa_-:=\{{\bf x}\in\partial\kappa, -{\mathbf \beta}({\bf x})\cdot{\bf n}({\bf x})<0\}$ denotes the inflow boundary -and $\partial\kappa_+=\partial\kappa\setminus \partial \kappa_-$ the -outflow part of cell $\kappa$. Below, this equation will be referred -to as first version of the DG method. We note that after a -second integration by parts, we obtain: find $u_h\in V_h$ such that -@f[ - \sum_\kappa\left\{(\nabla\cdot\{{\mathbf \beta} u_h\},v_h)_\kappa - -({\mathbf \beta}\cdot{\bf n} [u_h], v_h)_{\partial\kappa_-}\right\} - =(f,v_h)_\Omega, \quad\forall v_h\in V_h, -@f] -where $[u_h]=u_h^+-u_h^-$ denotes the jump of the discrete function -between two neighboring cells and is defined to be $[u_h]=u_h^+-g$ on -the boundary of the domain. This is the discontinuous Galerkin scheme -for the transport equation given in its original notation. -Nevertheless, we will base the implementation of the scheme on the -form given by [dg-general1] and [upwind-flux], -or [dg-transport1], respectively. - -Finally, we rewrite [dg-general1] in terms of a summation over all -faces where each face $e=\partial \kappa\cap\partial \kappa'$ -between two neighboring cells $\kappa$ and $\kappa'$ occurs twice: -Find $u_h\in V_h$ such that -@f[ - -\sum_\kappa({\mathcal F}(u_h),\nabla v_h)_\kappa+\sum_e\left\{({\mathcal H}(u_h^+,u_h^-,{\bf n}), v_h)_e+({\mathcal H}(u_h^-, u_h^+,-{\bf n}), v_h^-)_{e\setminus\Gamma}\right\}=(f,v_h)_\Omega \quad\forall v_h\in V_h, -\qquad\qquad\qquad\mathrm{[dg-general2]}@f] -By employing [conservativity] of the numerical flux -this equation simplifies to: find $u_h\in V_h$ such that -@f[ - -\sum_\kappa({\mathcal F}(u_h),\nabla v_h)_\kappa+\sum_e({\mathcal H}(u_h^+,u_h^-,{\bf n}), [v_h])_{e\setminus\Gamma}+({\mathcal H}(u_h,g,{\bf n}), v_h)_{\Gamma}=(f,v_h)_\Omega \quad\forall v_h\in V_h. -\qquad\qquad\qquad\mathrm{[dg-general3]}@f] -Whereas the outer unit normal ${\bf n}|_{\partial\kappa}$ is uniquely -defined this is not so for ${\bf n}_e$ as the latter might be the -normal from either side of the face. Hence, we need to fix the normal -${\bf n}$ on the face to be one of the two normals and denote the -other normal by $-{\bf n}$. This way we get $-{\bf n}$ in the second -face term in [dg-general2] that finally produces the -minus sign in the jump $[v_h]$ in equation [dg-general3]. - -For the linear [transport-equation] -equation [dg-general3] simplifies to -@f[ - -\sum_\kappa(u_h,{\mathbf \beta}\cdot\nabla v_h)_\kappa+\sum_e\left\{({\mathbf \beta}\cdot{\bf n}\, u_h, [v_h])_{e_+\setminus\Gamma}+({\mathbf \beta}\cdot{\bf n}\, u_h^-, [v_h])_{e_-\setminus\Gamma}\right\}=(f,v_h)_\Omega-({\mathbf \beta}\cdot{\bf n}\, g, v_h)_{\Gamma_-}, -\qquad\qquad\qquad\mathrm{[dg-transport2]}@f] -which will be refered to as second version of the DG method. - - -

Implementation

- - -As already mentioned at the beginning of this example we will -implement assembling the system matrix in two different ways. -The first one will be based on the first version [dg-transport1] -of the DG method that includes a sum of integrals over all cell -boundaries $\partial\kappa$. This is realized by a loop over all cells and -a nested loop over all faces of each cell. Thereby each inner face -$e=\partial\kappa\cap\partial \kappa'$ is treated twice, the first -time when the outer loop treats cell $\kappa$ and the second time when it -treats cell $\kappa'$. This way some values like the shape function -values at quadrature points on faces need to be computed twice. - -To overcome this overhead and for comparison, we implement -assembling of matrix also in a second and different way. This will -be based on the second version [dg-transport2] that -includes a sum of integrals over all faces $e$. Here, several -difficulties occurs. -
    -
  1. As degrees of freedom are associated with cells (and not to faces) - and as a normal is only defined w.r.t. a cell adjacent to the face we - cannot simply run over all faces of the triangulation but need to - perform the nested loop over all cells and all faces of each cell - like in the first implementation. This, because in deal.II - faces are accessible from cells but not visa versa. -
  2. Due to the nested loop we arrive twice at each face. In order to - assemble face terms only once we either need to track which - faces we have treated before, or we introduce a simple rule that decides - which of the two adjacent cells the face should be accessed and - treated from. Here, we employ the second approach and define the - following rule: -
      -
    1. If the two cells adjacent to a face are of the same refinement level we access and treat the face from the cell with lower index on this level. -
    2. If the two cells are of different refinement levels we access - and treat the face from the coarser cell. -
    -
-Before we start with the description of the code we first introduce -its main ingredients. The main class is called -DGMethod. It comprises all basic objects like the -triangulation, the dofhandler, the system matrix and solution vectors. -Furthermore it has got some member functions, the most prominent of -which are the assemble_system1 and assemble_system2 -functions that implement the two different ways mentioned above for -assembling the system matrix. Within these assembling routines several -different cases must be distinguished while performing the nested -loops over all cells and all faces of each cell and assembling the -respective face terms. While sitting on the current cell and looking -at a specific face there are the cases -
    -
  1. face is at boundary, -
  2. neighboring cell is finer, -
  3. neighboring cell is of the same refinement level, and -
  4. neighboring cell is coarser -
-where the `neighboring cell' and the current cell have the mentioned -faces in common. In last three cases the assembling of the face terms -are almost the same. Hence, we can implement the assembling of the -face terms either by `copy and paste' (the lazy way, whose -disadvantages come up when the scheme or the equation might want to be -changed afterwards) or by calling a separate function that covers all -three cases. To be kind of educational within this tutorial we perform -the latter approach, of course. We go even further and encapsulate -this function and everything that is needed for assembling the -specific equation under consideration within a class called -DGTransportEquation. This class includes objects of all -equation--specific functions, the RHS and the -BoundaryValues class, both derived from the Function -class, and the Beta class representing the vector field. -Furthermore, the DGTransportEquation class comprises member -functions assemble_face_terms1 and -assemble_face_terms2 that are invoked by the -assemble_system1 and assemble_system2 functions of the -DGMethod, respectively, and the functions -assemble_cell_term and assemble_boundary_term that -are the same for both assembling routines. Due to the encapsulation of -all equation- and scheme-specific functions, the -DGTransportEquation class can easily be replaced by a similar -class that implements a different equation and a different DG method. -Indeed, the implementation of the assemble_system1 and -assemble_system2 functions of the DGMethod class will -be general enough to serve for different DG methods, different -equations, even for systems of equations (!) and, under small -modifications, for nonlinear problems. Finally, we note that the -program is dimension independent, i.e. after replacing -DGMethod<2> by DGMethod<3> the code runs in 3d. - - - - - - - diff --git a/deal.II/examples/step-12/doc/results.dox b/deal.II/examples/step-12/doc/results.dox deleted file mode 100644 index ecf8b2c67f..0000000000 --- a/deal.II/examples/step-12/doc/results.dox +++ /dev/null @@ -1,91 +0,0 @@ -

Results

- - -The output of this program consist of the console output, the eps -files including the grids, and the solutions given in gnuplot format. -@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 ... -Writing solution to ... - -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 ... -Writing solution to ... - -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 ... -Writing solution to ... - -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 ... -Writing solution to ... - -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 ... -Writing solution to ... - -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 ... -Writing solution to ... -@endcode - -We see that, as expected, on each refinement step the two solutions -coincide. The difference measured in time of treating each face only -once (second version of the DG method) in comparison with treating -each face twice within a nested loop over all cells and all faces of -each cell (first version), is much less than one might have -expected. The gain is less than 20% on the last few refinement steps. - - - First we show the solutions on the initial mesh, the mesh after two -and after five adaptive refinement steps. - -@image html step-12.sol-0.png -@image html step-12.sol-2.png -@image html step-12.sol-5.png - - -Then we show the final grid (after 5 refinement steps) and the solution again, -this time with a nicer 3d rendering (obtained using the DataOutBase::write_vtk -function and the VTK-based VisIt visualization program) that better shows the -sharpness of the jump on the refined mesh and the over- and undershoots of the -solution along the interface: - -@image html step-12.grid-5.png -@image html step-12.3d-solution.png - - -And finally we show a plot of a 3d computation. - -@image html step-12.sol-5-3d.png - diff --git a/deal.II/examples/step-12/step-12.cc b/deal.II/examples/step-12/step-12.cc deleted file mode 100644 index 4d8f246278..0000000000 --- a/deal.II/examples/step-12/step-12.cc +++ /dev/null @@ -1,1682 +0,0 @@ -/* $Id$ */ -/* Author: Ralf Hartmann, University of Heidelberg, 2001 */ - -/* $Id$ */ -/* */ -/* Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors */ -/* */ -/* This file is subject to QPL and may not be distributed */ -/* without copyright and license information. Please refer */ -/* to the file deal.II/doc/license.html for the text and */ -/* further information on this license. */ - - // The first few files have already - // been covered in previous examples - // and will thus not be further - // commented on. -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - // This is the first new file. It - // declares the MappingQ1 class - // that gives the standard bilinear - // mapping. For bilinear mappings use - // an object of this class rather - // than an object of the - // MappingQ(1) class, as the - // MappingQ1 class is optimized - // due to the pre-knowledge of the - // actual polynomial degree 1. -#include - // Here the discontinuous finite - // elements are defined. They are - // used in the same way as all other - // finite elements, though -- as you - // have seen in previous tutorial - // programs -- there isn't much user - // interaction with finite element - // classes at all: the are passed to - // DoFHandler and FEValues - // objects, and that is about it. -#include - // We are going to use the simplest - // possible solver, called Richardson - // iteration, that represents a - // simple defect correction. This, in - // combination with a block SSOR - // preconditioner (defined in - // precondition_block.h), that uses - // the special block matrix structure - // of system matrices arising from DG - // discretizations. -#include -#include - // We are going to use gradients as - // refinement indicator. -#include - // Finally we do some time comparison - // using the Timer class. -#include - - // And this again is C++: -#include -#include - - // 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, RHS - // and BoundaryValues, are - // derived from the Function - // class. Only the value_list - // function are implemented because - // only lists of function values are - // computed rather than single - // values. -template -class RHS: public Function -{ - public: - virtual void value_list (const std::vector > &points, - std::vector &values, - const unsigned int component=0) const; -}; - - -template -class BoundaryValues: public Function -{ - public: - virtual void value_list (const std::vector > &points, - std::vector &values, - const unsigned int component=0) const; -}; - - - // The class Beta represents the - // vector valued flow field of the - // linear transport equation and is - // not derived from the Function - // class as we prefer to get function - // values of type Point rather - // than of type - // Vector@. This, because - // there exist scalar products - // between Point and Point as - // well as between Point and - // Tensor, simplifying terms like - // $\beta\cdot n$ and - // $\beta\cdot\nabla v$. - // - // An unnecessary empty constructor - // is added to the class to work - // around a bug in Compaq's cxx - // compiler which otherwise reports - // an error about an omitted - // initializer for an object of - // this class further down. -template -class Beta -{ - public: - Beta () {} - void value_list (const std::vector > &points, - std::vector > &values) const; -}; - - - // The implementation of the - // value_list functions of these - // classes are rather simple. For - // simplicity the right hand side is - // set to be zero but will be - // assembled anyway. -template -void RHS::value_list(const std::vector > &points, - std::vector &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 -void Beta::value_list(const std::vector > &points, - std::vector > &values) const -{ - Assert(values.size()==points.size(), - ExcDimensionMismatch(values.size(),points.size())); - - for (unsigned int i=0; i -void BoundaryValues::value_list(const std::vector > &points, - std::vector &values, - const unsigned int) const -{ - Assert(values.size()==points.size(), - ExcDimensionMismatch(values.size(),points.size())); - - for (unsigned int i=0; iDGTransportEquation. Its - // member functions were already - // mentioned in the Introduction and - // will be explained - // below. Furthermore it includes - // objects of the previously defined - // Beta, RHS and - // BoundaryValues function - // classes. -template -class DGTransportEquation -{ - public: - DGTransportEquation(); - - void assemble_cell_term(const FEValues& fe_v, - FullMatrix &ui_vi_matrix, - Vector &cell_vector) const; - - void assemble_boundary_term(const FEFaceValues& fe_v, - FullMatrix &ui_vi_matrix, - Vector &cell_vector) const; - - void assemble_face_term1(const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, - FullMatrix &ui_vi_matrix, - FullMatrix &ue_vi_matrix) const; - - void assemble_face_term2(const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, - FullMatrix &ui_vi_matrix, - FullMatrix &ue_vi_matrix, - FullMatrix &ui_ve_matrix, - FullMatrix &ue_ve_matrix) const; - private: - const Beta beta_function; - const RHS rhs_function; - const BoundaryValues boundary_function; -}; - - -template -DGTransportEquation::DGTransportEquation () - : - beta_function (), - rhs_function (), - boundary_function () -{} - - - // @sect4{Function: assemble_cell_term} - // - // The assemble_cell_term - // function assembles the cell terms - // of the discretization. - // ui_vi_matrix is a cell matrix, - // i.e. for a DG method of degree 1, - // it is of size 4 times 4, and - // cell_vector is of size 4. - // When this function is invoked, - // fe_v is already reinit'ed with the - // current cell before and includes - // all shape values needed. -template -void DGTransportEquation::assemble_cell_term( - const FEValues &fe_v, - FullMatrix &ui_vi_matrix, - Vector &cell_vector) const -{ - // First we ask fe_v for the - // quadrature weights, - const std::vector &JxW = fe_v.get_JxW_values (); - - // Then the flow field beta and the - // rhs_function are evaluated at - // the quadrature points, - std::vector > beta (fe_v.n_quadrature_points); - std::vector 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; pointassemble_boundary_term - // function assembles the face terms - // at boundary faces. When this - // function is invoked, fe_v is - // already reinit'ed with the current - // cell and current face. Hence it - // provides the shape values on that - // boundary face. -template -void DGTransportEquation::assemble_boundary_term( - const FEFaceValues& fe_v, - FullMatrix &ui_vi_matrix, - Vector &cell_vector) const -{ - // Again, as in the previous - // function, we ask the - // FEValues object for the - // quadrature weights - const std::vector &JxW = fe_v.get_JxW_values (); - // but here also for the normals. - const std::vector > &normals = fe_v.get_normal_vectors (); - - // We evaluate the flow field - // and the boundary values at the - // quadrature points. - std::vector > beta (fe_v.n_quadrature_points); - std::vector 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; point0) - for (unsigned int i=0; iassemble_face_term1 - // function assembles the face terms - // corresponding to the first version - // of the DG method, cf. above. For - // that case, the face terms are - // given as a sum of integrals over - // all cell boundaries. - // - // When this function is invoked, - // fe_v and fe_v_neighbor are - // already reinit'ed with the current - // cell and the neighoring cell, - // respectively, as well as with the - // current face. Hence they provide - // the inner and outer shape values - // on the face. - // - // In addition to the cell matrix - // ui_vi_matrix this function - // gets a new argument - // ue_vi_matrix, that stores - // contributions to the system matrix - // that are based on exterior values - // of $u$ and interior values of - // $v$. Here we note that ue is - // the short notation for u - // exterior and represents $u_h^-$, - // see the introduction. -template -void DGTransportEquation::assemble_face_term1( - const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, - FullMatrix &ui_vi_matrix, - FullMatrix &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 &JxW = fe_v.get_JxW_values (); - const std::vector > &normals = fe_v.get_normal_vectors (); - - // and we evaluate the flow field - // at the quadrature points. - std::vector > 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; point0) - for (unsigned int i=0; iassemble_face_term2 function - // that assembles the face terms - // corresponding to the second - // version of the DG method, - // cf. above. For that case the face - // terms are given as a sum of - // integrals over all faces. Here we - // need two additional cell matrices - // ui_ve_matrix and - // ue_ve_matrix that will store - // contributions due to terms - // involving ui and ve as well as ue - // and ve. -template -void DGTransportEquation::assemble_face_term2( - const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, - FullMatrix &ui_vi_matrix, - FullMatrix &ue_vi_matrix, - FullMatrix &ui_ve_matrix, - FullMatrix &ue_ve_matrix) const -{ - // the first few lines are the same - const std::vector &JxW = fe_v.get_JxW_values (); - const std::vector > &normals = fe_v.get_normal_vectors (); - - std::vector > beta (fe_v.n_quadrature_points); - - beta_function.value_list (fe_v.get_quadrature_points(), beta); - - for (unsigned int point=0; point0) - { - // This term we've already - // seen. - for (unsigned int i=0; iDGMethod 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 -class DGMethod -{ - public: - DGMethod (); - ~DGMethod (); - - void run (); - - private: - void setup_system (); - void assemble_system1 (); - void assemble_system2 (); - void solve (Vector &solution); - void refine_grid (); - void output_results (const unsigned int cycle) const; - - Triangulation triangulation; - const MappingQ1 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 fe; - DoFHandler dof_handler; - - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; - - // We define the quadrature - // formulae for the cell and the - // face terms of the - // discretization. - const QGauss quadrature; - const QGauss face_quadrature; - - // And there are two solution - // vectors, that store the - // solutions to the problems - // corresponding to the two - // different assembling routines - // assemble_system1 and - // assemble_system2; - Vector solution1; - Vector solution2; - Vector right_hand_side; - - // Finally this class includes an - // object of the - // DGTransportEquations class - // described above. - const DGTransportEquation dg; -}; - - -template -DGMethod::DGMethod () - : - mapping (), - // Change here for DG - // methods of - // different degrees. - fe (1), - dof_handler (triangulation), - quadrature (4), - face_quadrature (4), - dg () -{} - - -template -DGMethod::~DGMethod () -{ - dof_handler.clear (); -} - - -template -void DGMethod::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::faces_per_cell - *GeometryInfo::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 - // assemble_system1 function that - // implements the DG discretization - // in its first version. This - // function repeatedly calls the - // assemble_cell_term, - // assemble_boundary_term and - // assemble_face_term1 functions - // of the DGTransportEquation - // object. The - // assemble_boundary_term covers - // the first case mentioned in the - // introduction. - // - // 1. face is at boundary - // - // This function takes a - // FEFaceValues object as - // argument. In contrast to that - // assemble_face_term1 - // takes two FEFaceValuesBase - // objects; one for the shape - // functions on the current cell and - // the other for shape functions on - // the neighboring cell under - // consideration. Both objects are - // either of class FEFaceValues - // or of class FESubfaceValues - // (both derived from - // FEFaceValuesBase) according to - // the remaining cases mentioned - // in the introduction: - // - // 2. neighboring cell is finer - // (current cell: FESubfaceValues, - // neighboring cell: FEFaceValues); - // - // 3. neighboring cell is of the same - // refinement level (both, current - // and neighboring cell: - // FEFaceValues); - // - // 4. neighboring cell is coarser - // (current cell: FEFaceValues, - // neighboring cell: - // FESubfaceValues). - // - // If we considered globally refined - // meshes then only case 3 would - // occur. But as we consider also - // locally refined meshes we need to - // distinguish all four cases making - // the following assembling function - // a bit longish. -template -void DGMethod::assemble_system1 () -{ - const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; - std::vector dofs (dofs_per_cell); - std::vector dofs_neighbor (dofs_per_cell); - - // First we create the - // update_flags for the - // FEValues and the - // FEFaceValues objects. - const UpdateFlags update_flags = update_values - | update_gradients - | update_quadrature_points - | update_JxW_values; - - // Note, that on faces we do not - // need gradients but we need - // normal vectors. - const UpdateFlags face_update_flags = update_values - | update_quadrature_points - | update_JxW_values - | update_normal_vectors; - - // On the neighboring cell we only - // need the shape values. Given a - // specific face, the quadrature - // points and `JxW values' are the - // same as for the current cells, - // the normal vectors are known to - // be the negative of the normal - // vectors of the current cell. - const UpdateFlags neighbor_face_update_flags = update_values; - - // Then we create the FEValues - // object. Note, that since version - // 3.2.0 of deal.II the constructor - // of this class takes a - // Mapping object as first - // argument. Although the - // constructor without Mapping - // argument is still supported it - // is recommended to use the new - // constructor. This reduces the - // effect of `hidden magic' (the - // old constructor implicitely - // assumes a MappingQ1 mapping) - // and makes it easier to change - // the mapping object later. - FEValues fe_v ( - mapping, fe, quadrature, update_flags); - - // Similarly we create the - // FEFaceValues and - // FESubfaceValues objects for - // both, the current and the - // neighboring cell. Within the - // following nested loop over all - // cells and all faces of the cell - // they will be reinited to the - // current cell and the face (and - // subface) number. - FEFaceValues fe_v_face ( - mapping, fe, face_quadrature, face_update_flags); - FESubfaceValues fe_v_subface ( - mapping, fe, face_quadrature, face_update_flags); - FEFaceValues fe_v_face_neighbor ( - mapping, fe, face_quadrature, neighbor_face_update_flags); - FESubfaceValues fe_v_subface_neighbor ( - mapping, fe, face_quadrature, neighbor_face_update_flags); - - // Now we create the cell matrices - // and vectors. Here we need two - // cell matrices, both for face - // terms that include test - // functions vi (internal shape - // functions, i.e. shape functions - // of the current cell). To be more - // precise, the first matrix will - // include the `ui and vi terms' - // and the second will include the - // `ue and vi terms'. Here we - // recall the convention that `ui' - // is the shortcut for $u_h^+$ and - // `ue' represents $u_h^-$, see the - // introduction. - FullMatrix ui_vi_matrix (dofs_per_cell, dofs_per_cell); - FullMatrix ue_vi_matrix (dofs_per_cell, dofs_per_cell); - - Vector cell_vector (dofs_per_cell); - - // Furthermore we need some cell - // iterators. - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - // Now we start the loop over all - // active cells. - for (;cell!=endc; ++cell) - { - // In the - // assemble_face_term1 - // function contributions to - // the cell matrices and the - // cell vector are only - // ADDED. Therefore on each - // cell we need to reset the - // ui_vi_matrix and - // cell_vector to zero, - // before assembling the cell terms. - ui_vi_matrix = 0; - cell_vector = 0; - - // Now we reinit the FEValues - // object for the current cell - fe_v.reinit (cell); - - // and call the function - // that assembles the cell - // terms. The first argument is - // the FEValues that was - // previously reinit'ed on the - // current cell. - dg.assemble_cell_term(fe_v, - ui_vi_matrix, - cell_vector); - - // As in previous examples the - // vector `dofs' includes the - // dof_indices. - cell->get_dof_indices (dofs); - - // This is the start of the - // nested loop over all faces. - for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) - { - // First we set the face - // iterator - typename DoFHandler::face_iterator face=cell->face(face_no); - - // and clear the - // ue_vi_matrix on each - // face. - ue_vi_matrix = 0; - - // Now we distinguish the - // four different cases in - // the ordering mentioned - // above. We start with - // faces belonging to the - // boundary of the domain. - if (face->at_boundary()) - { - // We reinit the - // FEFaceValues - // object to the - // current face - fe_v_face.reinit (cell, face_no); - - // and assemble the - // corresponding face - // terms. - dg.assemble_boundary_term(fe_v_face, - ui_vi_matrix, - cell_vector); - } - else - { - // Now we are not on - // the boundary of the - // domain, therefore - // there must exist a - // neighboring cell. - typename DoFHandler::cell_iterator neighbor= - cell->neighbor(face_no);; - - // We proceed with the - // second and most - // complicated case: - // the neighboring cell - // is more refined than - // the current cell. As - // in deal.II - // neighboring cells - // are restricted to - // have a level - // difference of not - // more than one, the - // neighboring cell is - // known to be at most - // ONCE more refined - // than the current - // cell. Furthermore - // also the face is - // more refined, - // i.e. it has - // children. Here we - // note that the - // following part of - // code will not work - // for dim==1. - if (face->has_children()) - { - // First we store - // which number the - // current cell has - // in the list of - // neighbors of the - // neighboring - // cell. Hence, - // neighbor-@>neighbor(neighbor2) - // equals the - // current cell - // cell. - const unsigned int neighbor2= - cell->neighbor_of_neighbor(face_no); - - - // We loop over - // subfaces - for (unsigned int subface_no=0; - subface_non_children(); ++subface_no) - { - // and set the - // cell - // iterator - // neighbor_child - // to the cell - // placed - // `behind' the - // current - // subface. - typename DoFHandler::active_cell_iterator - neighbor_child - = cell->neighbor_child_on_subface (face_no, subface_no); - - Assert (!neighbor_child->has_children(), ExcInternalError()); - - // We need to - // reset the - // ue_vi_matrix - // on each - // subface - // because on - // each subface - // the un - // belong to - // different - // neighboring - // cells. - ue_vi_matrix = 0; - - // As already - // mentioned - // above for - // the current - // case (case - // 2) we employ - // the - // FESubfaceValues - // of the - // current - // cell (here - // reinited for - // the current - // cell, face - // and subface) - // and we - // employ the - // FEFaceValues - // of the - // neighboring - // child cell. - fe_v_subface.reinit (cell, face_no, subface_no); - fe_v_face_neighbor.reinit (neighbor_child, neighbor2); - - dg.assemble_face_term1(fe_v_subface, - fe_v_face_neighbor, - ui_vi_matrix, - ue_vi_matrix); - - // Then we get - // the dof - // indices of - // the - // neighbor_child - // cell - neighbor_child->get_dof_indices (dofs_neighbor); - - // and - // distribute - // ue_vi_matrix - // to the - // system_matrix - for (unsigned int i=0; iif - // (face-@>has_children()) - } - else - { - // We proceed with - // case 3, - // i.e. neighboring - // cell is of the - // same refinement - // level as the - // current cell. - if (neighbor->level() == cell->level()) - { - // Like before - // we store - // which number - // the current - // cell has in - // the list of - // neighbors of - // the - // neighboring - // cell. - const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no); - - // We reinit - // the - // FEFaceValues - // of the - // current and - // neighboring - // cell to the - // current face - // and assemble - // the - // corresponding - // face terms. - fe_v_face.reinit (cell, face_no); - fe_v_face_neighbor.reinit (neighbor, neighbor2); - - dg.assemble_face_term1(fe_v_face, - fe_v_face_neighbor, - ui_vi_matrix, - ue_vi_matrix); - // End of if - // (neighbor-@>level() - // == - // cell-@>level()) - } - else - { - // Finally we - // consider - // case 4. When - // the - // neighboring - // cell is not - // finer and - // not of the - // same - // refinement - // level as the - // current cell - // it must be - // coarser. - Assert(neighbor->level() < cell->level(), ExcInternalError()); - - // Find out the - // how many'th - // face_no and - // subface_no - // the current - // face is - // w.r.t. the - // neighboring - // cell. - const std::pair faceno_subfaceno= - cell->neighbor_of_coarser_neighbor(face_no); - const unsigned int neighbor_face_no=faceno_subfaceno.first, - neighbor_subface_no=faceno_subfaceno.second; - - Assert (neighbor->neighbor_child_on_subface (neighbor_face_no, - neighbor_subface_no) - == cell, - ExcInternalError()); - - // Reinit the - // appropriate - // FEFaceValues - // and assemble - // the face - // terms. - fe_v_face.reinit (cell, face_no); - fe_v_subface_neighbor.reinit (neighbor, neighbor_face_no, - neighbor_subface_no); - - dg.assemble_face_term1(fe_v_face, - fe_v_subface_neighbor, - ui_vi_matrix, - ue_vi_matrix); - } - - // Now we get the - // dof indices of - // the - // neighbor_child - // cell, - neighbor->get_dof_indices (dofs_neighbor); - - // and distribute the - // ue_vi_matrix. - for (unsigned int i=0; iface not at boundary: - } - // End of loop over all faces: - } - - // Finally we distribute the - // ui_vi_matrix - for (unsigned int i=0; iassemble_system2 function that - // implements the DG discretization - // in its second version. This - // function is very similar to the - // assemble_system1 - // function. Therefore, here we only - // discuss the differences between - // the two functions. This function - // repeatedly calls the - // assemble_face_term2 function - // of the DGTransportEquation object, - // that assembles the face terms - // written as a sum of integrals over - // all faces. Therefore, we need to - // make sure that each face is - // treated only once. This is achieved - // by introducing the rule: - // - // a) If the current and the - // neighboring cells are of the same - // refinement level we access and - // treat the face from the cell with - // lower index. - // - // b) If the two cells are of - // different refinement levels we - // access and treat the face from the - // coarser cell. - // - // Due to rule b) we do not need to - // consider case 4 (neighboring cell - // is coarser) any more. - -template -void DGMethod::assemble_system2 () -{ - const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; - std::vector dofs (dofs_per_cell); - std::vector dofs_neighbor (dofs_per_cell); - - const UpdateFlags update_flags = update_values - | update_gradients - | update_quadrature_points - | update_JxW_values; - - const UpdateFlags face_update_flags = update_values - | update_quadrature_points - | update_JxW_values - | update_normal_vectors; - - const UpdateFlags neighbor_face_update_flags = update_values; - - // Here we do not need - // fe_v_face_neighbor as case 4 - // does not occur. - FEValues fe_v ( - mapping, fe, quadrature, update_flags); - FEFaceValues fe_v_face ( - mapping, fe, face_quadrature, face_update_flags); - FESubfaceValues fe_v_subface ( - mapping, fe, face_quadrature, face_update_flags); - FEFaceValues fe_v_face_neighbor ( - mapping, fe, face_quadrature, neighbor_face_update_flags); - - - FullMatrix ui_vi_matrix (dofs_per_cell, dofs_per_cell); - FullMatrix ue_vi_matrix (dofs_per_cell, dofs_per_cell); - - // Additionally we need the - // following two cell matrices, - // both for face term that include - // test function ve (external - // shape functions, i.e. shape - // functions of the neighboring - // cell). To be more precise, the - // first matrix will include the `u - // and vn terms' and the second - // that will include the `un and vn - // terms'. - FullMatrix ui_ve_matrix (dofs_per_cell, dofs_per_cell); - FullMatrix ue_ve_matrix (dofs_per_cell, dofs_per_cell); - - Vector cell_vector (dofs_per_cell); - - // The following lines are roughly - // the same as in the previous - // function. - typename DoFHandler::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::faces_per_cell; ++face_no) - { - typename DoFHandler::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::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_non_children(); ++subface_no) - { - typename DoFHandler::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; ilevel() == 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 -void DGMethod::solve (Vector &solution) -{ - SolverControl solver_control (1000, 1e-12, false, false); - SolverRichardson<> solver (solver_control); - - // Here we create the - // preconditioner, - PreconditionBlockSSOR > preconditioner; - - // we assigned the matrix to it and - // set the right block size. - preconditioner.initialize(system_matrix, fe.dofs_per_cell); - - // After these preparations we are - // ready to start the linear solver. - solver.solve (system_matrix, solution, right_hand_side, - preconditioner); -} - - - // We refine the grid according to a - // very simple refinement criterion, - // namely an approximation to the - // gradient of the solution. As here - // we consider the DG(1) method - // (i.e. we use piecewise bilinear - // shape functions) we could simply - // compute the gradients on each - // cell. But we do not want to base - // our refinement indicator on the - // gradients on each cell only, but - // want to base them also on jumps of - // the discontinuous solution - // function over faces between - // neighboring cells. The simpliest - // way of doing that is to compute - // approximative gradients by - // difference quotients including the - // cell under consideration and its - // neighbors. This is done by the - // DerivativeApproximation class - // that computes the approximate - // gradients in a way similar to the - // GradientEstimation described - // in step-9 of this tutorial. In - // fact, the - // DerivativeApproximation class - // was developed following the - // GradientEstimation class of - // step-9. Relating to the - // discussion in step-9, here we - // consider $h^{1+d/2}|\nabla_h - // u_h|$. Futhermore we note that we - // do not consider approximate second - // derivatives because solutions to - // the linear advection equation are - // in general not in $H^2$ but in $H^1$ - // (to be more precise, in $H^1_\beta$) - // only. -template -void DGMethod::refine_grid () -{ - // The DerivativeApproximation - // class computes the gradients to - // float precision. This is - // sufficient as they are - // approximate and serve as - // refinement indicators only. - Vector 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::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 -void DGMethod::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 data_out; - data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (solution2, "u"); - - data_out.build_patches (); - - data_out.write_gnuplot(gnuplot_output); -} - - - // The following run function is - // similar to previous examples. The - // only difference is that the - // problem is assembled and solved - // twice on each refinement step; - // first by assemble_system1 that - // implements the first version and - // then by assemble_system2 that - // implements the second version of - // writing the DG - // discretization. Furthermore the - // time needed by each of the two - // assembling routines is measured. -template -void DGMethod::run () -{ - for (unsigned int cycle=0; cycle<6; ++cycle) - { - std::cout << "Cycle " << cycle << ':' << std::endl; - - if (cycle == 0) - { - GridGenerator::hyper_cube (triangulation); - - triangulation.refine_global (3); - } - else - refine_grid (); - - - std::cout << " Number of active cells: " - << triangulation.n_active_cells() - << std::endl; - - setup_system (); - - std::cout << " Number of degrees of freedom: " - << dof_handler.n_dofs() - << std::endl; - - // The constructor of the Timer - // class automatically starts - // the time measurement. - Timer assemble_timer; - // First assembling routine. - assemble_system1 (); - // The operator () accesses the - // current time without - // disturbing the time - // measurement. - std::cout << "Time of assemble_system1: " - << assemble_timer() - << std::endl; - solve (solution1); - - // As preparation for the - // second assembling routine we - // reinit the system matrix, the - // right hand side vector and - // the Timer object. - system_matrix = 0; - right_hand_side = 0; - assemble_timer.reset(); - - // We start the Timer, - assemble_timer.start(); - // call the second assembling routine - assemble_system2 (); - // and access the current time. - std::cout << "Time of assemble_system2: " - << assemble_timer() - << std::endl; - solve (solution2); - - // To make sure that both - // versions of the DG method - // yield the same - // discretization and hence the - // same solution we check the - // two solutions for equality. - solution1-=solution2; - const double difference=solution1.linfty_norm(); - if (difference>1e-13) - std::cout << "solution1 and solution2 differ!!" << std::endl; - else - std::cout << "solution1 and solution2 coincide." << std::endl; - - // Finally we perform the - // output. - output_results (cycle); - } -} - - // The following main function is - // similar to previous examples and - // need not to be commented on. -int main () -{ - try - { - DGMethod<2> dgmethod; - dgmethod.run (); - } - catch (std::exception &exc) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - catch (...) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - }; - - return 0; -} - -