]> https://gitweb.dealii.org/ - dealii.git/commitdiff
New tutorial step 42 created
authorBaerbel Jannsen <baerbel.janssen@gmail.com>
Mon, 4 Jan 2010 13:08:22 +0000 (13:08 +0000)
committerBaerbel Jannsen <baerbel.janssen@gmail.com>
Mon, 4 Jan 2010 13:08:22 +0000 (13:08 +0000)
git-svn-id: https://svn.dealii.org/trunk@20276 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-42/Makefile [new file with mode: 0644]
deal.II/examples/step-42/doc/intro.dox [new file with mode: 0644]
deal.II/examples/step-42/doc/results.dox [new file with mode: 0644]
deal.II/examples/step-42/step-22.cc [new file with mode: 0644]

diff --git a/deal.II/examples/step-42/Makefile b/deal.II/examples/step-42/Makefile
new file mode 100644 (file)
index 0000000..a443603
--- /dev/null
@@ -0,0 +1,156 @@
+# $Id: Makefile,v 1.4 2006/02/10 17:53:05 wolf Exp $
+
+
+# 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 *vtk
+
+
+
+
+#
+#
+# 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-deal2-3d.g) \
+          $(lib-lac.g)      \
+           $(lib-base.g)
+libs.o   = $(lib-deal2-2d.o) \
+          $(lib-deal2-3d.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-42/doc/intro.dox b/deal.II/examples/step-42/doc/intro.dox
new file mode 100644 (file)
index 0000000..24c3d91
--- /dev/null
@@ -0,0 +1,783 @@
+<br>
+
+<i>This program was contributed by Martin Kronbichler and Wolfgang
+Bangerth. 
+<br>
+This material is based upon work partly supported by the National
+Science Foundation under Award No. EAR-0426271 and The California Institute of 
+Technology. Any opinions, findings, and conclusions or recommendations
+expressed in this publication are those of the author and do not
+necessarily reflect the views of the National Science Foundation or of The
+California Institute of Technology.
+</i> 
+
+
+
+<a name="Intro"></a>
+<h1>Introduction</h1>
+
+This program deals with the Stokes system of equations which reads as
+follows in non-dimensionalized form:
+@f{eqnarray*}
+  -\textrm{div}\; \varepsilon(\textbf{u}) + \nabla p &=& \textbf{f},
+  \\
+  -\textrm{div}\; \textbf{u} &=& 0,
+@f}
+where $\textbf u$ denotes the velocity of a fluid, $p$ is its
+pressure, $\textbf f$ are external forces, and
+$\varepsilon(\textbf{u})= \nabla^s{\textbf{u}}= \frac 12 \left[
+(\nabla \textbf{u}) + (\nabla \textbf{u})^T\right]$  is the
+rank-2 tensor of symmetrized gradients; a component-wise definition
+of it is $\varepsilon(\textbf{u})_{ij}=\frac
+12\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right)$.
+
+The Stokes equations describe the steady-state motion of a
+slow-moving, viscous fluid such as honey, rocks in the earth mantle,
+or other cases where inertia does not play a significant role. If a
+fluid is moving fast enough that inertia forces are significant
+compared to viscous friction, the Stokes equations are no longer
+valid; taking into account inertia effects then leads to the
+nonlinear Navier-Stokes equations. However, in this tutorial program,
+we will focus on the simpler Stokes system.
+
+To be well-posed, we will have to add boundary conditions to the
+equations. What boundary conditions are readily possible here will
+become clear once we discuss the weak form of the equations.
+
+The equations covered here fall into the class of vector-valued problems. A
+toplevel overview of this topic can be found in the @ref vector_valued module.
+
+
+<h3>Weak form</h3>
+
+The weak form of the equations is obtained by writing it in vector
+form as
+@f{eqnarray*}
+  \left(
+  {-\textrm{div}\; \varepsilon(\textbf{u}) + \nabla p}
+  \atop
+  {-\textrm{div}\; \textbf{u}}
+  \right)
+  =
+  \left(
+  {\textbf{f}}
+  \atop
+  0
+  \right),
+@f}
+forming the dot product from the left with a vector-valued test
+function $\phi = \left({\textbf v \atop q}\right)$ and integrating
+over the domain $\Omega$, yielding the following set of equations:
+@f{eqnarray*}
+  (\mathrm v,
+   -\textrm{div}\; \varepsilon(\textbf{u}) + \nabla p)_{\Omega}
+  -
+  (q,\textrm{div}\; \textbf{u})_{\Omega}
+  =
+  (\textbf{v}, \textbf{f})_\Omega,
+@f}
+which has to hold for all test functions $\phi = \left({\textbf v
+\atop q}\right)$.
+
+In practice, one wants to impose as little regularity on the pressure
+variable as possible; consequently, we integrate by parts the second term:
+@f{eqnarray*}
+  (\mathrm v, -\textrm{div}\; \varepsilon(\textbf{u}))_{\Omega}
+  - (\textrm{div}\; \textbf{v}, p)_{\Omega}
+  + (\textbf{n}\cdot\textbf{v}, p)_{\partial\Omega}
+  -
+  (q,\textrm{div}\; \textbf{u})_{\Omega}
+  =
+  (\textbf{v}, \textbf{f})_\Omega.
+@f}
+Likewise, we integrate by parts the first term to obtain
+@f{eqnarray*}
+  (\nabla \mathrm v,\varepsilon(\textbf{u}))_{\Omega}
+  -
+  (\textbf{n} \otimes \mathrm v,\varepsilon(\textbf{u}))_{\partial\Omega}
+  - (\textrm{div}\; \textbf{v}, p)_{\Omega}
+  + (\textbf{n}\cdot\textbf{v}, p)_{\partial\Omega}
+  -
+  (q,\textrm{div}\; \textbf{u})_{\Omega}
+  =
+  (\textbf{v}, \textbf{f})_\Omega,
+@f}
+where the scalar product between two tensor-valued quantities is here
+defined as
+@f{eqnarray*}
+  (\nabla \mathrm v,\varepsilon(\textbf{u}))_{\Omega}
+  =
+  \int_\Omega \sum_{i,j=1}^d \frac{\partial v_j}{\partial x_i}
+  \varepsilon(\textbf{u})_{ij} \ dx.
+@f}
+Because the scalar product between a general tensor like
+$\nabla\mathrm v$ and a symmetric tensor like
+$\varepsilon(\textbf{u})$ equals the scalar product between the
+symmetrized forms of the two, we can also write the bilinear form
+above as follows:
+@f{eqnarray*}
+  (\varepsilon(\mathrm v),\varepsilon(\textbf{u}))_{\Omega}
+  -
+  (\textbf{n} \otimes \mathrm v,\varepsilon(\textbf{u}))_{\partial\Omega}
+  - (\textrm{div}\; \textbf{v}, p)_{\Omega}
+  + (\textbf{n}\cdot\textbf{v}, p)_{\partial\Omega}
+  -
+  (q,\textrm{div}\; \textbf{u})_{\Omega}
+  =
+  (\textbf{v}, \textbf{f})_\Omega,
+@f}
+We will deal with the boundary terms in the next section, but it is already
+clear from the domain terms
+@f{eqnarray*}
+  (\varepsilon(\mathrm v),\varepsilon(\textbf{u}))_{\Omega}
+  - (\textrm{div}\; \textbf{v}, p)_{\Omega}
+  -
+  (q,\textrm{div}\; \textbf{u})_{\Omega}
+@f}
+of the bilinear form that the Stokes equations yield a symmetric bilinear
+form, and consequently a symmetric (if indefinite) system matrix.
+
+
+<h3>%Boundary conditions</h3>
+
+The weak form just derived immediately presents us with different
+possibilities for imposing boundary conditions:
+<ol>
+<li>Dirichlet velocity boundary conditions: On a part
+    $\Gamma_D\subset\partial\Omega$ we may impose Dirichlet conditions
+    on the velocity $\textbf u$:
+
+    @f{eqnarray*}
+        \textbf u = \textbf g_D \qquad\qquad \textrm{on}\ \Gamma_D.
+    @f}
+    Because test functions $\textbf v$ come from the tangent space of
+    the solution variable, we have that $\textbf v=0$ on $\Gamma_D$
+    and consequently that
+    @f{eqnarray*}
+      -(\textbf{n} \otimes \mathrm
+        v,\varepsilon(\textbf{u}))_{\Gamma_D}
+      +
+      (\textbf{n}\cdot\textbf{v}, p)_{\Gamma_D}
+      = 0.
+    @f}
+    In other words, as usual, strongly imposed boundary values do not
+    appear in the weak form.
+    It is noteworthy that if we impose Dirichlet boundary values on the entire
+    boundary, then the pressure is only determined up to a constant. An
+    algorithmic realization of that would use similar tools as have been seen in
+    step-11.
+
+<li>Neumann-type or natural boundary conditions: On the rest of the boundary
+    $\Gamma_N=\partial\Omega\backslash\Gamma_D$, let us re-write the
+    boundary terms as follows:
+    @f{eqnarray*}
+      -(\textbf{n} \otimes \mathrm
+        v,\varepsilon(\textbf{u}))_{\Gamma_N}
+      +
+      (\textbf{n}\cdot\textbf{v}, p)_{\Gamma_N}
+      &=&
+      \sum_{i,j=1}^d
+      -(n_i v_j,\varepsilon(\textbf{u})_{ij})_{\Gamma_N}
+      +
+      \sum_{i=1}^d 
+      (n_i v_i, p)_{\Gamma_N}  
+      \\
+      &=&
+      \sum_{i,j=1}^d
+      -(n_i v_j,\varepsilon(\textbf{u})_{ij})_{\Gamma_N}
+      +
+      \sum_{i,j=1}^d 
+      (n_i v_j, p \delta_{ij})_{\Gamma_N}      
+      \\
+      &=&
+      \sum_{i,j=1}^d
+      (n_i v_j,p \delta_{ij} - \varepsilon(\textbf{u})_{ij})_{\Gamma_N}
+      \\
+      &=&
+      (\textbf{n} \otimes \mathrm v,
+      p \textbf{1} - \varepsilon(\textbf{u}))_{\Gamma_N}.
+      \\
+      &=&
+      (\mathrm v,
+       \textbf{n}\cdot [p \textbf{1} - \varepsilon(\textbf{u})])_{\Gamma_N}.
+    @f}
+    In other words, on the Neumann part of the boundary we can
+    prescribe values for the total stress:
+    @f{eqnarray*}
+      \textbf{n}\cdot [p \textbf{1} - \varepsilon(\textbf{u})]
+      = 
+      \textbf g_N \qquad\qquad \textrm{on}\ \Gamma_N.
+    @f}
+    If the boundary is subdivided into Dirichlet and Neumann parts
+    $\Gamma_D,\Gamma_N$, this then leads to the following weak form:
+    @f{eqnarray*}
+      (\varepsilon(\mathrm v),\varepsilon(\textbf{u}))_{\Omega}
+      - (\textrm{div}\; \textbf{v}, p)_{\Omega}
+      -
+      (q,\textrm{div}\; \textbf{u})_{\Omega}
+      =
+      (\textbf{v}, \textbf{f})_\Omega
+      -
+      (\textbf{v}, \textbf g_N)_{\Gamma_N}.
+    @f}
+    
+    
+<li>Robin-type boundary conditions: Robin boundary conditions are a mixture of
+    Dirichlet and Neumann boundary conditions. They would read 
+    @f{eqnarray*}
+      \textbf{n}\cdot [p \textbf{1} - \varepsilon(\textbf{u})]
+      = 
+      \textbf S \textbf u \qquad\qquad \textrm{on}\ \Gamma_R,
+    @f}
+    with a rank-2 tensor (matrix) $\textbf S$. The associated weak form is
+    @f{eqnarray*}
+      (\varepsilon(\mathrm v),\varepsilon(\textbf{u}))_{\Omega}
+      - (\textrm{div}\; \textbf{v}, p)_{\Omega}
+      -
+      (q,\textrm{div}\; \textbf{u})_{\Omega}
+      +
+      (\textbf S \textbf u, \textbf{v})_{\Gamma_R}
+      =
+      (\textbf{v}, \textbf{f})_\Omega.
+    @f}
+
+<li>Partial boundary conditions: It is possible to combine Dirichlet and
+    Neumann boundary conditions by only enforcing each of them for certain
+    components of the velocity. For example, one way to impose artificial
+    boundary conditions is to require that the flow is perpendicular to the
+    boundary, i.e. the tangential component $\textbf u_{\textbf t}=(\textbf
+    1-\textbf n\otimes\textbf n)\textbf u$ be zero, thereby constraining
+    <code>dim</code>-1 components of the velocity. The remaining component can
+    be constrained by requiring that the normal component of the normal 
+    stress be zero, yielding the following set of boundary conditions:
+    @f{eqnarray*}
+      \textbf u_{\textbf t} &=& 0,
+      \\
+      \textbf n \cdot \left(\textbf{n}\cdot [p \textbf{1} -
+      \varepsilon(\textbf{u})] \right)
+      &=& 
+      0.
+    @f}
+
+    An alternative to this is when one wants the flow to be <i>parallel</i>
+    rather than perpendicular to the boundary (in deal.II, the
+    VectorTools::compute_no_normal_flux_constraints function can do this for
+    you). This is frequently the case for problems with a free boundary
+    (e.g. at the surface of a river or lake if vertical forces of the flow are
+    not large enough to actually deform the surface), or if no significant
+    friction is exerted by the boundary on the fluid (e.g. at the interface
+    between earth mantle and earth core where two fluids meet that are
+    stratified by different densities but that both have small enough
+    viscosities to not introduce much tangential stress on each other). 
+    In formulas, this means that 
+    @f{eqnarray*}
+      \textbf{n}\cdot\textbf u &=& 0,
+      \\
+      (\textbf 1-\textbf n\otimes\textbf n)
+      \left(\textbf{n}\cdot [p \textbf{1} -
+      \varepsilon(\textbf{u})] \right)
+      &=& 
+      0,
+    @f}
+    the first condition (which needs to be imposed strongly) fixing a single
+    component of the velocity, with the second (which would be enforced in the
+    weak form) fixing the remaining two components.
+</ol>
+
+Despite this wealth of possibilities, we will only use Dirichlet and
+(homogenous) Neumann boundary conditions in this tutorial program.
+
+
+<h3>Discretization</h3>
+
+As developed above, the weak form of the equations with Dirichlet and Neumann
+boundary conditions on $\Gamma_D$ and $\Gamma_N$ reads like this: find
+$\textbf u\in \textbf V_g = \{\varphi \in H^1(\Omega)^d: \varphi_{\Gamma_D}=\textbf
+g_D\}, p\in Q=L^2(\Omega)$ so that
+@f{eqnarray*}
+  (\varepsilon(\mathrm v),\varepsilon(\textbf{u}))_{\Omega}
+  - (\textrm{div}\; \textbf{v}, p)_{\Omega}
+  -
+  (q,\textrm{div}\; \textbf{u})_{\Omega}
+  =
+  (\textbf{v}, \textbf{f})_\Omega
+  -
+  (\textbf{v}, \textbf g_N)_{\Gamma_N}
+@f}
+for all test functions 
+$\textbf v\in \textbf V_0 = \{\varphi \in H^1(\Omega)^d: \varphi_{\Gamma_D}=0\},q\in
+Q$.
+
+These equations represent a symmetric saddle point problem. It is well known
+that then a solution only exists if the function spaces in which we search for
+a solution have to satisfy certain conditions, typically referred to as the
+Babuska-Brezzi or Ladyzhenskaya-Babuska-Brezzi (LBB) conditions. The continuous
+function spaces above satisfy them. However, when we discretize the equations by
+replacing the continuous variables and test functions by finite element
+functions in finite dimensional spaces $\textbf V_{g,h}\subset \textbf V_g,
+Q_h\subset Q$, we have to make sure that $\textbf V_h,Q_h$ also satisfy the LBB
+conditions. This is similar to what we had to do in @ref step_20 "step-20".
+
+For the Stokes equations, there are a number of possible choices to ensure
+that the finite element spaces are compatible with the LBB condition. A simple
+and accurate choice that we will use here is $\textbf u_h\in Q_{p+1}^d,
+p_h\in Q_p$, i.e. use elements one order higher for the velocities than for the
+pressures. 
+
+This then leads to the following discrete problem: find $\textbf u_h,p_h$ so
+that
+@f{eqnarray*}
+  (\varepsilon(\mathrm v_h),\varepsilon(\textbf u_h))_{\Omega}
+  - (\textrm{div}\; \textbf{v}_h, p_h)_{\Omega}
+  -
+  (q_h,\textrm{div}\; \textbf{u}_h)_{\Omega}
+  =
+  (\textbf{v}_h, \textbf{f})_\Omega
+  -
+  (\textbf{v}_h, \textbf g_N)_{\Gamma_N}
+@f}
+for all test functions $\textbf v_h, q_h$. Assembling the linear system
+associated with this problem follows the same lines used in @ref step_20
+"step-20", @ref step_21 "step-21", and explained in detail in the @ref
+vector_valued module.
+
+
+
+<h3>Linear solver and preconditioning issues</h3>
+
+The weak form of the discrete equations naturally leads to the following
+linear system for the nodal values of the velocity and pressure fields:
+@f{eqnarray*}
+  \left(\begin{array}{cc}
+    A & B^T \\ B & 0
+  \end{array}\right)
+  \left(\begin{array}{c}
+    U \\ P
+  \end{array}\right)
+  =
+  \left(\begin{array}{c}
+    F \\ G
+  \end{array}\right),
+@f}
+Like in @ref step_20 "step-20" and @ref step_21 "step-21", we will solve this
+system of equations by forming the Schur complement, i.e. we will first find
+the solution $P$ of 
+@f{eqnarray*}
+  BA^{-1}B^T P &=& BA^{-1} F - G, \\
+@f}
+and then
+@f{eqnarray*}
+  AU &=& F - B^TP.
+@f}
+The way we do this is pretty much exactly like we did in these previous
+tutorial programs, i.e. we use the same classes <code>SchurComplement</code>
+and <code>InverseMatrix</code> again. There are two significant differences,
+however:
+
+<ol>
+<li>
+First, in the mixed Laplace equation we had to deal with the question of how
+to precondition the Schur complement $B^TM^{-1}B$, which was spectrally
+equivalent to the Laplace operator on the pressure space (because $B$
+represents the gradient operator, $B^T$ its adjoint $-\textrm{div}$, and $M$
+the identity (up to the material parameter $K^{-1}$), so $B^TM^{-1}B$ is
+something like $-\textrm{div} \mathbf 1 \nabla = -\Delta$). Consequently, the
+matrix is badly conditioned for small mesh sizes and we had to come up with an
+elaborate preconditioning scheme for the Schur complement.
+
+<li>
+Second, every time we multiplied with $B^TM^{-1}B$ we had to solve with the
+mass matrix $M$. This wasn't particularly difficult, however, since the mass
+matrix is always well conditioned and so simple to invert using CG and a
+little bit of preconditioning.
+</ol>
+In other words, preconditioning the inner solver for $M$ was simple whereas
+preconditioning the outer solver for $B^TM^{-1}B$ was complicated.
+
+Here, the situation is pretty much exactly the opposite. The difference stems
+from the fact that the matrix at the heart of the Schur complement does not
+stem from the identity operator but from a variant of the Laplace operator,
+$-\textrm{div} \nabla^s$ (where $\nabla^s$ is the symmetric gradient)
+acting on a vector field. In the investigation of this issue
+we largely follow the paper D. Silvester and A. Wathen:
+"Fast iterative solution of stabilised Stokes systems part II. Using
+general block preconditioners." (SIAM J. Numer. Anal., 31 (1994),
+pp. 1352-1367), which is available online <a
+href="http://siamdl.aip.org/getabs/servlet/GetabsServlet?prog=normal&id=SJNAAM000031000005001352000001&idtype=cvips&gifs=Yes" target="_top">here</a>.
+Principally, the difference in the matrix at the heart of the Schur
+complement has two consequences:
+
+<ol>
+<li>
+First, it makes the outer preconditioner simple: the Schur complement
+corresponds to the operator $-\textrm{div} (-\textrm{div} \nabla^s)^{-1}
+\nabla$ on the pressure space; forgetting about the fact that we deal with
+symmetric gradients instead of the regular one, the Schur complement is
+something like $-\textrm{div} (-\textrm{div} \nabla)^{-1} \nabla =
+-\textrm{div} (-\Delta)^{-1} \nabla$, which, even if not mathematically
+entirely concise, is spectrally equivalent to the identity operator (a
+heuristic argument would be to commute the operators into
+$-\textrm{div}(-\Delta)^{-1} \nabla = -\textrm{div}\nabla(-\Delta)^{-1} =
+-\Delta(-\Delta)^{-1} = \mathbf 1$). It turns out that it isn't easy to solve
+this Schur complement in a straightforward way with the CG method:
+using no preconditioner, the condition number of the Schur complement matrix
+depends on the size ratios of the largest to the smallest cells, and one still
+needs on the order of 50-100 CG iterations. However, there is a simple cure:
+precondition with the mass matrix on the pressure space and we get down to a
+number between 5-15 CG iterations, pretty much independently of the structure
+of the mesh (take a look at the <a href="#Results">results section</a> of this
+program to see that indeed the number of CG iterations does not change as we
+refine the mesh).
+
+So all we need in addition to what we already have is the mass matrix on the
+pressure variables. We could do that by building this matrix on the
+side in a separate data structure. However, it is worth remembering
+that although we build the system matrix
+@f{eqnarray*}
+  \left(\begin{array}{cc}
+    A & B^T \\ B & 0
+  \end{array}\right)
+@f}
+as one object (of type BlockSparseMatrix), we never actually do
+matrix-vector products with this matrix, or any other operations that
+consider the entire matrix. Rather, we only build it in this form for
+convenience (because it reflects the structure of the FESystem finite
+element and associated DoFHandler object) but later only operate on
+the $(0,0),(0,1)$, and $(1,0)$ blocks of this matrix. In other words,
+our algorithm so far entirely ignores the $(1,1)$ (pressure-pressure)
+block as it is empty anyway. 
+
+Now, as mentioned, we need a pressure mass matrix to precondition the
+Schur complement and that conveniently the pressure-pressure block of
+the matrix we build anyway is currently empty and ignored. So what we
+will do is to assemble the needed mass matrix in this space; this does
+change the global system matrix but since our algorithm never operates
+on the global matrix and instead only considers individual blocks,
+this fact does not affect what we actually compute. Later, when
+solving, we then precondition the Schur complement with $M_p^{-1}$ by
+doing a few CG iterations on the well-conditioned pressure mass matrix
+$M_p$ stored in the $(1,1)$ block.
+
+
+
+<li>
+While the outer preconditioner has become simpler compared to the
+mixed Laplace case discussed in @ref step_20 "step-20", the issue of
+the inner solver has become more complicated. In the mixed Laplace
+discretization, the Schur complement has the form $B^TM^{-1}B$. Thus,
+every time we multiplied with the Schur complement, we had to solve a
+linear system $M_uz=y$; this isn't too complicated there, however,
+since the mass matrix $M_u$ on the pressure space is well-conditioned.
+
+
+On the other hand, for the Stokes equation we consider here, the Schur
+complement is $BA^{-1}B^T$ where the matrix $A$ is related to the
+Laplace operator (it is, in fact, the matrix corresponding to the
+bilinear form $(\nabla^s \varphi_i, \nabla^s\varphi_j)$). Thus,
+solving with $A$ is a lot more complicated: the matrix is badly
+conditioned and we know that we need many iterations unless we have a
+very good preconditioner. What is worse, we have to solve with $A$
+every time we multiply with the Schur complement, which is 5-15 times
+using the preconditioner described above.
+
+Because we have to solve with $A$ several times, it pays off to spend
+a bit more time once to create a good preconditioner for this
+matrix. So here's what we're going to do: if in 2d, we use the
+ultimate preconditioner, namely a direct sparse LU decomposition of
+the matrix. This is implemented using the SparseDirectUMFPACK class
+that uses the UMFPACK direct solver to compute the decomposition. To
+use it, you will have to specify the <code>--enable-umfpack</code>
+switch when configuring the deal.II library, see the <a
+href="../../readme.html#optional-software">ReadMe file</a> for
+instructions. With this, the inner solver converges in one iteration.
+
+In 2d, we can do this sort of thing because even reasonably large problems
+rarely have more than a few 100,000 unknowns with relatively few nonzero
+entries per row. Furthermore, the bandwidth of matrices in 2d is ${\cal
+O}(\sqrt{N})$ and therefore moderate. For such matrices, sparse factors can be
+computed in a matter of a few seconds. (As a point of reference, computing the
+sparse factors of a matrix of size $N$ and bandwidth $B$ takes ${\cal
+O}(NB^2)$ operations. In 2d, this is ${\cal O}(N^2)$; though this is a higher
+complexity than, for example, assembling the linear system which takes ${\cal
+O}(N)$, the constant for computing the decomposition is so small that it
+doesn't become the dominating factor in the entire program until we get to
+very large %numbers of unknowns in the high 100,000s or more.)
+
+The situation changes in 3d, because there we quickly have many more
+unknowns and the bandwidth of matrices (which determines the number of
+nonzero entries in sparse LU factors) is ${\cal O}(N^{2/3})$, and there
+are many more entries per row as well. This makes using a sparse
+direct solver such as UMFPACK inefficient: only for problem sizes of a
+few 10,000 to maybe 100,000 unknowns can a sparse decomposition be
+computed using reasonable time and memory resources.
+
+What we do in that case is to use an incomplete LU decomposition (ILU) as a
+preconditioner, rather than actually computing complete LU factors. As it so
+happens, deal.II has a class that does this: SparseILU. Computing the ILU
+takes a time that only depends on the number of nonzero entries in the sparse
+matrix (or that we are willing to fill in the LU factors, if these should be
+more than the ones in the matrix), but is independent of the bandwidth of the
+matrix. It is therefore an operation that can efficiently also be computed in
+3d. On the other hand, an incomplete LU decomposition, by definition, does not
+represent an exact inverse of the matrix $A$. Consequently, preconditioning
+with the ILU will still require more than one iteration, unlike
+preconditioning with the sparse direct solver. The inner solver will therefore
+take more time when multiplying with the Schur complement, a tradeoff
+unavoidable.
+</ol>
+
+In the program below, we will make use of the fact that the SparseILU and
+SparseDirectUMFPACK classes have a very similar interface and can be used
+interchangeably. All that we need is a switch class that, depending on the
+dimension, provides a type that is either of the two classes mentioned
+above. This is how we do that:
+@code
+template <int dim>
+struct InnerPreconditioner;
+
+template <>
+struct InnerPreconditioner<2> 
+{
+    typedef SparseDirectUMFPACK type;
+};
+
+template <>
+struct InnerPreconditioner<3> 
+{
+    typedef SparseILU<double> type;
+};
+@endcode
+
+From hereon, we can refer to the type <code>typename
+InnerPreconditioner@<dim@>::type</code> and automatically get the correct
+preconditioner class. Because of the similarity of the interfaces of the two
+classes, we will be able to use them interchangeably using the same syntax in
+all places.
+
+
+<h3>The testcase</h3>
+
+The domain, right hand side and boundary conditions we implement below relate
+to a problem in geophysics: there, one wants to compute the flow field of
+magma in the earth's interior under a mid-ocean rift. Rifts are places where
+two continental plates are very slowly drifting apart (a few centimeters per
+year at most), leaving a crack in the earth crust that is filled with magma
+from below. Without trying to be entirely realistic, we model this situation
+by solving the following set of equations and boundary conditions on the
+domain $\Omega=[-2,2]\times[0,1]\times[-1,0]$:
+@f{eqnarray*}
+  -\textrm{div}\; \varepsilon(\textbf{u}) + \nabla p &=& 0,
+  \\
+  -\textrm{div}\; \textbf{u} &=& 0,
+  \\
+  \mathbf u &=&   \left(\begin{array}{c}
+    -1 \\ 0 \\0
+  \end{array}\right)
+  \qquad\qquad \textrm{at}\ z=0, x<0,
+  \\
+  \mathbf u &=&   \left(\begin{array}{c}
+    +1 \\ 0 \\0
+  \end{array}\right)
+  \qquad\qquad \textrm{at}\ z=0, x>0,
+  \\
+  \mathbf u &=&   \left(\begin{array}{c}
+    0 \\ 0 \\0
+  \end{array}\right)
+  \qquad\qquad \textrm{at}\ z=0, x=0,
+@f}
+and using natural boundary conditions $\textbf{n}\cdot [p \textbf{1} -
+\varepsilon(\textbf{u})] = 0$ everywhere else. In other words, at the
+left part of the top surface we prescribe that the fluid moves with the
+continental plate to the left at speed $-1$, that it moves to the right on the
+right part of the top surface, and impose natural flow conditions everywhere
+else. If we are in 2d, the description is essentially the same, with the
+exception that we omit the second component of all vectors stated above.
+
+As will become apparent in the <a href="#Results">results section</a>, the
+flow field will pull material from below and move it to the left and right
+ends of the domain, as expected. The discontinuity of velocity boundary
+conditions will produce a singularity in the pressure at the center of the top
+surface that sucks material all the way to the top surface to fill the gap
+left by the outward motion of material at this location.
+
+
+<h3>Implementation</h3>
+
+<h4>Using imhomogeneous constraints for implementing Dirichlet boundary conditions</h4>
+
+In all the previous tutorial programs, we used the ConstraintMatrix merely
+for handling hanging node constraints (with exception of step-11). However,
+the class can also be used to implement Dirichlet boundary conditions, as we
+will show in this program, by fixing some node values $x_i = b_i$. Note that
+these are inhomogeneous constraints, and we have to pay some special
+attention to that. The way we are going to implement this is to first read
+in the boundary values into the ConstraintMatrix object by using the call
+
+@code
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                            1,
+                                            BoundaryValues<dim>(),
+                                            constraints);
+@endcode 
+
+very similar to how we were making the list of boundary nodes
+before (note that we set Dirichlet conditions only on boundaries with
+boundary flag 1). The actual application of the boundary values is then
+handled by the ConstraintMatrix object directly, without any additional
+interference.
+
+We could then proceed as before, namely by filling the matrix, and then
+calling a condense function on the constraints object of the form
+@code
+  constraints.condense (system_matrix, system_rhs);
+@endcode
+
+Note that we call this on the system matrix and system right hand side
+simultaneously, since resolving inhomogeneous constraints requires knowledge
+about both the matrix entries and the right hand side. For efficiency
+reasons, though, we choose another strategy: all the constraints collected
+in the ConstraintMatrix can be resolved on the fly while writing local data
+into the global matrix, by using the call
+@code
+  constraints.distribute_local_to_global (local_matrix, local_rhs,
+                                          local_dof_indices,
+                                          system_matrix, system_rhs);
+@endcode
+
+This technique is further discussed in the @ref step_27 "step-27" tutorial
+program. All we need to know here is that this functions does three things
+at once: it writes the local data into the global matrix and right hand
+side, it distributes the hanging node constraints and additionally
+implements (inhomogeneous) Dirichlet boundary conditions. That's nice, isn't
+it?
+
+We can conclude that the ConstraintMatrix provides an alternative to using
+MatrixTools::apply_boundary_values for implementing Dirichlet boundary
+conditions.
+
+
+<a name="constraint-matrix">
+<h4>Using ConstraintMatrix for increasing performance</h4>
+</a>
+
+Frequently, a sparse matrix contains a substantial amount of elements that
+acutally are zero when we are about to start a linear solve. Such elements are
+introduced when we eliminate constraints or implement Dirichlet conditions,
+where we usually delete all entries in constrained rows and columns, i.e., we
+set them to zero. The fraction of elements that are present in the sparsity
+pattern, but do not really contain any information, can be up to one fourth
+of the total number of elements in the matrix for the 3D application
+considered in this tutorial program. Remember that matrix-vector products or
+preconditioners operate on all the elements of a sparse matrix (even those 
+that are zero), which is an inefficiency we will avoid here.
+
+An advantage of directly resolving constrained degrees of freedom is that we
+can avoid having most of the entries that are going to be zero in our sparse
+matrix &mdash; we do not need constrained entries during matrix construction 
+(as opposed to the traditional algorithms, which first fill the matrix, and 
+only resolve constraints afterwards). This will save both memory and time 
+when forming matrix-vector products. The way we are going to do that is to 
+pass the information about constraints to the function that generates the 
+sparsity pattern, and then set a <tt>false</tt> argument specifying that we 
+do not intend to use constrained entries:
+@code
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern,
+                                   constraints, false);
+@endcode
+This functions obviates, by the way, also the call to the
+<tt>condense()</tt> function on the sparsity pattern.
+
+
+<h4>Performance optimizations</h4>
+
+The program developed below has seen a lot of TLC. We have run it over and
+over under profiling tools (mainly <a
+href="http://www.valgrind.org/">valgrind</a>'s cachegrind and callgrind
+tools, as well as the KDE <a
+href="http://kcachegrind.sourceforge.net/">KCachegrind</a> program for
+visualization) to see where the bottlenecks are. This has paid off: through
+this effort, the program has become about four times as fast when
+considering the runtime of the refinement cycles zero through three,
+reducing the overall number of CPU instructions executed from
+869,574,060,348 to 199,853,005,625. For higher refinement levels, the gain
+is probably even larger since some algorithms that are not ${\cal O}(N)$
+have been eliminated.
+
+Essentially, there are currently two algorithms in the program that do not
+scale linearly with the number of degrees of freedom: renumbering of degrees
+of freedom (which is ${\cal O}(N \log N)$, and the linear solver (which is
+${\cal O}(N^{4/3})$). As for the first, while reordering degrees of freedom
+may not scale linearly, it is an indispensible part of the overall algorithm
+as it greatly improves the quality of the sparse ILU, easily making up for
+the time spent on computing the renumbering; graphs and timings to
+demonstrate this are shown in the documentation of the DoFRenumbering
+namespace, also underlining the choice of the Cuthill-McKee reordering
+algorithm chosen below.
+
+As for the linear solver: as mentioned above, our implementation here uses a
+Schur complement formulation. This is not necessarily the very best choice
+but demonstrates various important techniques available in deal.II. The
+question of which solver is best is again discussed in the <a
+href="#improved-solver">section on improved solvers in the results part</a>
+of this program, along with code showing alternative solvers and a
+comparison of their results.
+
+Apart from this, many other algorithms have been tested and improved during
+the creation of this program. For example, in building the sparsity pattern,
+we originally used a BlockCompressedSparsityPattern object that added one
+element at a time; however, its data structures are poorly adapted for the
+large numbers of nonzero entries per row created by our discretization in
+3d, leading to a quadratic behavior. Replacing the internal algorithms in
+deal.II to set many elements at a time, and using a
+BlockCompressedSimpleSparsityPattern as a better adapted data structure,
+removed this bottleneck at the price of a slightly higher memory
+consumption. Likewise, the implementation of the decomposition step in the
+SparseILU class was very inefficient and has been replaced by one that is
+about 10 times faster. Even the vmult function of the SparseILU has been
+improved to save about twenty precent of time. Small improvements were
+applied here and there. Moreover, the ConstraintMatrix object has been used
+to eliminate a lot of entries in the sparse matrix that are eventually going
+to be zero, see <a href="#constraint-matrix">the section on using advanced
+features of the ConstraintMatrix class</a>.
+
+A profile of how many CPU instructions are spent at the various
+different places in the program during refinement cycles
+zero through three in 3d is shown here:
+
+@image html step-22.profile-3.png
+
+As can be seen, at this refinement level approximately three quarters of the
+instruction count is spent on the actual solver (the SparseILU::vmult calls
+on the left, the SparseMatrix::vmult call in the middle for the Schur
+complement solve, and another box representing the multiplications with
+SparseILU and SparseMatrix in the solve for <i>U</i>). About one fifth of
+the instruction count is spent on matrix assembly and sparse ILU computation
+(box in the lower right corner) and the rest on other things. Since floating
+point operations such as in the SparseILU::vmult calls typically take much
+longer than many of the logical operations and table lookups in matrix
+assembly, the fraction of the run time taken up by matrix assembly is
+actually significantly less than the fraction of instructions, as will
+become apparent in the comparison we make in the results section.
+
+For higher refinement levels, the boxes representing the solver as well as
+the blue box at the top right stemming from reordering algorithm are going
+to grow at the expense of the other parts of the program, since they don't
+scale linearly. The fact that at this moderate refinement level (3168 cells
+and 93176 degrees of freedom) the linear solver already makes up about three
+quarters of the instructions is a good sign that most of the algorithms used
+in this program are well-tuned and that major improvements in speeding up
+the program are most likely not to come from hand-optimizing individual
+aspects but by changing solver algorithms. We will address this point in the
+discussion of results below as well.
+
+As a final point, and as a point of reference, the following picture also
+shows how the profile looked at an early stage of optimizing this program:
+
+@image html step-22.profile-3.original.png
+
+As mentioned above, the runtime of this version was about four times as long
+as for the first profile, with the SparseILU decomposition taking up about
+30% of the instruction count, and operations on the ill-suited
+CompressedSparsityPattern about 10%. Both these bottlenecks have since been
+completely removed.
diff --git a/deal.II/examples/step-42/doc/results.dox b/deal.II/examples/step-42/doc/results.dox
new file mode 100644 (file)
index 0000000..2afe3ca
--- /dev/null
@@ -0,0 +1,809 @@
+<a name="Results"></a>
+<h1>Results</h1>
+
+<h3>Output of the program and graphical visualization</h3>
+
+<h4>2D calculations</h4>
+
+Running the program with the space dimension set to 2 in <code>main()</code>
+yields the following output (when the flag is set to optimized in the
+Makefile):
+@code
+examples/step-22> make run
+============================ Remaking Makefile.dep
+==============optimized===== step-22.cc
+============================ Linking step-22
+============================ Running step-22
+Refinement cycle 0
+   Number of active cells: 64
+   Number of degrees of freedom: 679 (594+85)
+   Assembling...
+   Computing preconditioner...
+   Solving...  11 outer CG Schur complement iterations for pressure
+
+Refinement cycle 1
+   Number of active cells: 160
+   Number of degrees of freedom: 1683 (1482+201)
+   Assembling...
+   Computing preconditioner...
+   Solving...  11 outer CG Schur complement iterations for pressure
+
+Refinement cycle 2
+   Number of active cells: 376
+   Number of degrees of freedom: 3813 (3370+443)
+   Assembling...
+   Computing preconditioner...
+   Solving...  11 outer CG Schur complement iterations for pressure
+
+Refinement cycle 3
+   Number of active cells: 880
+   Number of degrees of freedom: 8723 (7722+1001)
+   Assembling...
+   Computing preconditioner...
+   Solving...  11 outer CG Schur complement iterations for pressure
+
+Refinement cycle 4
+   Number of active cells: 2008
+   Number of degrees of freedom: 19383 (17186+2197)
+   Assembling...
+   Computing preconditioner...
+   Solving...  11 outer CG Schur complement iterations for pressure
+
+Refinement cycle 5
+   Number of active cells: 4288
+   Number of degrees of freedom: 40855 (36250+4605)
+   Assembling...
+   Computing preconditioner...
+   Solving...  11 outer CG Schur complement iterations for pressure
+@endcode
+
+The entire computation above takes about 20 seconds on a reasonably
+quick (for 2007 standards) machine.
+
+What we see immediately from this is that the number of (outer)
+iterations does not increase as we refine the mesh. This confirms the
+statement in the introduction that preconditioning the Schur
+complement with the mass matrix indeed yields a matrix spectrally
+equivalent to the identity matrix (i.e. with eigenvalues bounded above
+and below independently of the mesh size or the relative sizes of
+cells). In other words, the mass matrix and the Schur complement are
+spectrally equivalent.
+
+In the images below, we show the grids for the first six refinement
+steps in the program.  Observe how the grid is refined in regions
+where the solution rapidly changes: On the upper boundary, we have
+Dirichlet boundary conditions that are -1 in the left half of the line
+and 1 in the right one, so there is an aprupt change at $x=0$. Likewise,
+there are changes from Dirichlet to Neumann data in the two upper
+corners, so there is need for refinement there as well:
+
+<TABLE WIDTH="60%" ALIGN="center">
+  <tr>
+    <td ALIGN="center">
+      @image html step-22.2d.mesh-0.png
+    </td>
+
+    <td ALIGN="center">
+      @image html step-22.2d.mesh-1.png
+    </td>
+  </tr>
+
+  <tr>
+    <td ALIGN="center">
+      @image html step-22.2d.mesh-2.png
+    </td>
+
+    <td ALIGN="center">
+      @image html step-22.2d.mesh-3.png
+    </td>
+  </tr>
+
+  <tr>
+    <td ALIGN="center">
+      @image html step-22.2d.mesh-4.png
+    </td>
+
+    <td ALIGN="center">
+      @image html step-22.2d.mesh-5.png
+    </td>
+  </tr>
+</table>
+
+Finally, following is a plot of the flow field. It shows fluid
+transported along with the moving upper boundary and being replaced by
+material coming from below:
+
+@image html step-22.2d.solution.png
+
+This plot uses the capability of VTK-based visualization programs (in
+this case of VisIt) to show vector data; this is the result of us
+declaring the velocity components of the finite element in use to be a
+set of vector components, rather than independent scalar components in
+the <code>StokesProblem@<dim@>::output_results</code> function of this
+tutorial program.
+
+
+
+<h4>3D calculations</h4>
+
+In 3d, the screen output of the program looks like this:
+
+@code
+Refinement cycle 0
+   Number of active cells: 32
+   Number of degrees of freedom: 1356 (1275+81)
+   Assembling...
+   Computing preconditioner...
+   Solving...  13 outer CG Schur complement iterations for pressure.
+
+Refinement cycle 1
+   Number of active cells: 144
+   Number of degrees of freedom: 5088 (4827+261)
+   Assembling...
+   Computing preconditioner...
+   Solving...  14 outer CG Schur complement iterations for pressure.
+
+Refinement cycle 2
+   Number of active cells: 704
+   Number of degrees of freedom: 22406 (21351+1055)
+   Assembling...
+   Computing preconditioner...
+   Solving...  14 outer CG Schur complement iterations for pressure.
+
+Refinement cycle 3
+   Number of active cells: 3168
+   Number of degrees of freedom: 93176 (89043+4133)
+   Assembling...
+   Computing preconditioner...
+   Solving...  15 outer CG Schur complement iterations for pressure.
+
+Refinement cycle 4
+   Number of active cells: 11456
+   Number of degrees of freedom: 327808 (313659+14149)
+   Assembling...
+   Computing preconditioner...
+   Solving...  15 outer CG Schur complement iterations for pressure.
+
+Refinement cycle 5
+   Number of active cells: 45056
+   Number of degrees of freedom: 1254464 (1201371+53093)
+   Assembling...
+   Computing preconditioner...
+   Solving...  14 outer CG Schur complement iterations for pressure.
+@endcode
+
+Again, we see that the number of outer iterations does not increase as
+we refine the mesh. Nevertheless, the compute time increases
+significantly: for each of the iterations above separately, it takes a
+few seconds, a few seconds, 30sec, 4min, 15min, and 1h18min. This overall
+superlinear (in the number of unknowns) increase in runtime is due to the fact
+that our inner solver is not ${\cal O}(N)$: a simple experiment shows
+that as we keep refining the mesh, the average number of
+ILU-preconditioned CG iterations to invert the velocity-velocity block
+$A$ increases.
+
+We will address the question of how possibly to improve our solver <a
+href="#improved-solver">below</a>.
+
+As for the graphical output, the grids generated during the solution
+look as follow: 
+
+<TABLE WIDTH="60%" ALIGN="center">
+  <tr>
+    <td ALIGN="center">
+      @image html step-22.3d.mesh-0.png
+    </td>
+
+    <td ALIGN="center">
+      @image html step-22.3d.mesh-1.png
+    </td>
+  </tr>
+
+  <tr>
+    <td ALIGN="center">
+      @image html step-22.3d.mesh-2.png
+    </td>
+
+    <td ALIGN="center">
+      @image html step-22.3d.mesh-3.png
+    </td>
+  </tr>
+
+  <tr>
+    <td ALIGN="center">
+      @image html step-22.3d.mesh-4.png
+    </td>
+
+    <td ALIGN="center">
+      @image html step-22.3d.mesh-5.png
+    </td>
+  </tr>
+</table>
+
+Again, they show essentially the location of singularities introduced
+by boundary conditions. The vector field computed makes for an
+interesting graph:
+
+@image html step-22.3d.solution.png
+
+The isocountours shown here as well are those of the pressure
+variable, showing the singularity at the point of discontinuous
+velocity boundary conditions.
+
+
+
+<h3>Sparsity pattern</h3>
+
+As explained during the generation of the sparsity pattern, it is
+important to have the numbering of degrees of freedom in mind when
+using preconditioners like incomplete LU decompositions. This is most
+conveniently visualized using the distribution of nonzero elements in
+the stiffness matrix.
+
+If we don't do anything special to renumber degrees of freedom (i.e.,
+without using DoFRenumbering::Cuthill_McKee, but with using
+DoFRenumbering::component_wise to ensure that degrees of freedom are
+appropriately sorted into their corresponding blocks of the matrix and
+vector), then we get the following image after the first adaptive
+refinement in two dimensions:
+
+@image html step-22.2d.sparsity-nor.png
+
+In order to generate such a graph, you have to insert a piece of
+code like the following to the end of the setup step. 
+@code
+  {
+    std::ofstream out ("sparsity_pattern.gpl");
+    sparsity_pattern.print_gnuplot(out);
+  }
+@endcode
+
+It is clearly visible that the nonzero entries are spread over almost the
+whole matrix.  This makes preconditioning by ILU inefficient: ILU generates a
+Gaussian elimination (LU decomposition) without fill-in elements, which means
+that more tentative fill-ins left out will result in a worse approximation of
+the complete decomposition.
+
+In this program, we have thus chosen a more advanced renumbering of
+components.  The renumbering with DoFRenumbering::Cuthill_McKee and grouping
+the components into velocity and pressure yields the following output:
+
+@image html step-22.2d.sparsity-ren.png
+
+It is apparent that the situation has improved a lot. Most of the elements are
+now concentrated around the diagonal in the (0,0) block in the matrix. Similar
+effects are also visible for the other blocks. In this case, the ILU
+decomposition will be much closer to the full LU decomposition, which improves
+the quality of the preconditioner. (It may be interesting to note that the
+sparse direct solver UMFPACK does some %internal renumbering of the equations
+before actually generating a sparse LU decomposition; that procedure leads to
+a very similar pattern to the one we got from the Cuthill-McKee algorithm.)
+
+Finally, we want to have a closer
+look at a sparsity pattern in 3D. We show only the (0,0) block of the
+matrix, again after one adaptive refinement. Apart from the fact that the matrix
+size has increased, it is also visible that there are many more entries
+in the matrix. Moreover, even for the optimized renumbering, there will be a
+considerable amount of tentative fill-in elements. This illustrates why UMFPACK 
+is not a good choice in 3D - a full decomposition needs many new entries that
+ eventually won't fit into the physical memory (RAM):
+
+@image html step-22.3d.sparsity_uu-ren.png
+
+
+
+<h3>Possible Extensions</h3>
+
+<a name="improved-solver">
+<h4>Improved linear solver in 3D</h4>
+</a>
+
+We have seen in the section of computational results that the number of outer
+iterations does not depend on the mesh size, which is optimal in a sense of
+scalability. This does, however, not apply to the solver as a whole, as
+mentioned above:
+We did not look at the number of inner iterations when generating the inverse of
+the matrix $A$ and the mass matrix $M_p$. Of course, this is unproblematic in
+the 2D case where we precondition $A$ with a direct solver and the
+<code>vmult</code> operation of the inverse matrix structure will converge in
+one single CG step, but this changes in 3D where we only use an ILU
+preconditioner.  There, the number of required preconditioned CG steps to
+invert $A$ increases as the mesh is refined, and each <code>vmult</code> 
+operation involves on average approximately 14, 23, 36, 59, 75 and 101 inner
+CG iterations in the refinement steps shown above. (On the other hand,
+the number of iterations for applying the inverse pressure mass matrix is
+always around five, both in two and three dimensions.)  To summarize, most work 
+is spent on solving linear systems with the same matrix $A$ over and over again.
+What makes this look even worse is the fact that we 
+actually invert a matrix that is about 95 precent the size of the total system
+matrix and stands for 85 precent of the non-zero entries in the sparsity
+pattern. Hence, the natural question is whether it is reasonable to solve a
+linear system with matrix $A$ for about 15 times when calculating the solution 
+to the block system.
+  
+The answer is, of course, that we can do that in a few other (most of the time
+better) ways.
+Nevertheless, it has to be remarked that an indefinite system as the one
+at hand puts indeed much higher
+demands on the linear algebra than standard elliptic problems as we have seen
+in the early tutorial programs. The improvements are still rather 
+unsatisfactory, if one compares with an elliptic problem of similar
+size. Either way, we will introduce below a number of improvements to the
+linear solver, a discussion that we will re-consider again with additional
+options in the @ref step_31 "step-31" program.
+
+<a name="improved-ilu">
+<h5>Better ILU decomposition by smart reordering</h5>
+</a>
+A first attempt to improve the speed of the linear solution process is to choose
+a dof reordering that makes the ILU being closer to a full LU decomposition, as
+already mentioned in the in-code comments. The DoFRenumbering namespace compares
+several choices for the renumbering of dofs for the Stokes equations. The best
+result regarding the computing time was found for the King ordering, which is
+accessed through the call DoFRenumbering::boost::king_ordering. With that
+program, the inner solver needs considerably less operations, e.g. about 62 
+inner CG iterations for the inversion of $A$ at cycle 4 compared to about 75 
+iterations with the standard Cuthill-McKee-algorithm. Also, the computing time 
+at cycle 4 decreased from about 17 to 11 minutes for the <code>solve()</code> 
+call. However, the King ordering (and the orderings provided by the
+DoFRenumbering::boost namespace in general) has a serious drawback - it uses
+much more memory than the in-build deal versions, since it acts on abstract
+graphs rather than the geometry provided by the triangulation. In the present
+case, the renumbering takes about 5 times as much memory, which yields an
+infeasible algorithm for the last cycle in 3D with 1.2 million
+unknowns.
+
+<h5>Better preconditioner for the inner CG solver</h5>
+Another idea to improve the situation even more would be to choose a 
+preconditioner that makes CG for the (0,0) matrix $A$ converge in a 
+mesh-independent number of iterations, say 10 to 30. We have seen such a
+canditate in @ref step_16 "step-16": multigrid.
+
+<h5>Block Schur complement preconditioner</h5>
+Even with a good preconditioner for $A$, we still 
+need to solve of the same linear system repeatedly (with different 
+right hand sides, though) in order to make the Schur complement solve 
+converge. The approach we are going to discuss here is how inner iteration 
+and outer iteration can be combined. If we persist in calculating the Schur 
+complement, there is no other possibility.
+
+The alternative is to attack the block system at once and use an approximate
+Schur complement as efficient preconditioner. The idea is as
+follows: If we find a block preconditioner $P$ such that the matrix
+@f{eqnarray*}
+  P^{-1}\left(\begin{array}{cc}
+    A & B^T \\ B & 0
+  \end{array}\right)
+@f}
+is simple, then an iterative solver with that preconditioner will converge in a
+few iterations. Using the Schur complement $S = B A^{-1} B^T$, one finds that
+@f{eqnarray*}
+  P^{-1}
+  = 
+  \left(\begin{array}{cc}
+    A^{-1} & 0 \\ S^{-1} B A^{-1} & -S^{-1}
+  \end{array}\right)
+@f}
+would appear to be a good choice since
+@f{eqnarray*}
+  P^{-1}\left(\begin{array}{cc}
+    A & B^T \\ B & 0
+  \end{array}\right) 
+  = 
+  \left(\begin{array}{cc}
+    A^{-1} & 0 \\ S^{-1} B A^{-1} & -S^{-1}
+  \end{array}\right)\cdot \left(\begin{array}{cc}
+    A & B^T \\ B & 0
+  \end{array}\right) 
+  = 
+  \left(\begin{array}{cc}
+    I & A^{-1} B^T \\ 0 & I
+  \end{array}\right).
+@f}
+This is the approach taken by the paper by Silvester and Wathen referenced
+to in the introduction (with the exception that Silvester and Wathen use
+right preconditioning). In this case, a Krylov-based iterative method would
+converge in one step only if exact inverses of $A$ and $S$ were applied,
+since all the eigenvalues are one (and the number of iterations in such a
+method is bounded by the number of distinct eigenvalues). Below, we will
+discuss the choice of an adequate solver for this problem. First, we are
+going to have a closer look at the implementation of the preconditioner.
+
+Since $P$ is aimed to be a preconditioner only, we shall use approximations to 
+the inverse of the Schur complement $S$ and the matrix $A$. Hence, the Schur
+complement will be approximated by the pressure mass matrix $M_p$, and we use
+a preconditioner to $A$ (without an InverseMatrix class around it) for
+approximating $A^{-1}$.
+
+Here comes the class that implements the block Schur
+complement preconditioner. The <code>vmult</code> operation for block vectors
+according to the derivation above can be specified by three successive
+operations:
+@code
+template <class PreconditionerA, class PreconditionerMp>
+class BlockSchurPreconditioner : public Subscriptor
+{
+  public:
+    BlockSchurPreconditioner (const BlockSparseMatrix<double>         &S,
+          const InverseMatrix<SparseMatrix<double>,PreconditionerMp>  &Mpinv,
+          const PreconditionerA &Apreconditioner);
+
+  void vmult (BlockVector<double>       &dst,
+              const BlockVector<double> &src) const;
+
+  private:
+    const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
+    const SmartPointer<const InverseMatrix<SparseMatrix<double>, 
+                       PreconditionerMp > > m_inverse;
+    const PreconditionerA &a_preconditioner;
+    
+    mutable Vector<double> tmp;
+
+};
+
+template <class PreconditionerA, class PreconditionerMp>
+BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::BlockSchurPreconditioner(
+          const BlockSparseMatrix<double>                            &S,
+          const InverseMatrix<SparseMatrix<double>,PreconditionerMp> &Mpinv,
+          const PreconditionerA &Apreconditioner
+          )
+                :
+                system_matrix           (&S),
+                m_inverse               (&Mpinv),
+                a_preconditioner        (Apreconditioner),
+                tmp                     (S.block(1,1).m())
+{}
+
+        // Now the interesting function, the multiplication of
+        // the preconditioner with a BlockVector. 
+template <class PreconditionerA, class PreconditionerMp>
+void BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::vmult (
+                                     BlockVector<double>       &dst,
+                                     const BlockVector<double> &src) const
+{
+        // Form u_new = A^{-1} u
+  a_preconditioner.vmult (dst.block(0), src.block(0));
+        // Form tmp = - B u_new + p 
+        // (<code>SparseMatrix::residual</code>
+        // does precisely this)
+  system_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1));
+        // Change sign in tmp
+  tmp *= -1;
+        // Multiply by approximate Schur complement 
+        // (i.e. a pressure mass matrix)
+  m_inverse->vmult (dst.block(1), tmp);
+}
+@endcode
+
+Since we act on the whole block system now, we have to live with one
+disadvantage: we need to perform the solver iterations on
+the full block system instead of the smaller pressure space.
+
+Now we turn to the question which solver we should use for the block
+system. The first observation is that the resulting preconditioned matrix cannot
+be solved with CG since it is neither positive definite nor symmetric.
+
+The deal.II libraries implement several solvers that are appropriate for the
+problem at hand. One choice is the solver @ref SolverBicgstab "BiCGStab", which
+was used for the solution of the unsymmetric advection problem in step-9. The
+second option, the one we are going to choose, is @ref SolverGMRES "GMRES"
+(generalized minimum residual). Both methods have their pros and cons - there
+are problems where one of the two candidates clearly outperforms the other, and 
+vice versa.
+<a href="http://en.wikipedia.org/wiki/GMRES#Comparison_with_other_solvers">Wikipedia</a>'s 
+article on the GMRES method gives a comparative presentation.
+A more comprehensive and well-founded comparsion can be read e.g. in the book by
+J.W. Demmel (Applied Numerical Linear Algebra, SIAM, 1997, section 6.6.6).
+
+For our specific problem with the ILU preconditioner for $A$, we certainly need
+to perform hundreds of iterations on the block system for large problem sizes
+(we won't beat CG!). Actually, this disfavors GMRES: During the GMRES
+iterations, a basis of Krylov vectors is successively built up and some
+operations are performed on these vectors. The more vectors are in this basis, 
+the more operations and memory will be needed. The number of operations scales
+as ${\cal O}(n + k^2)$ and memory as ${\cal O}(kn)$, where $k$ is the number of 
+vectors in the Krylov basis and $n$ the size of the (block) matrix.
+To not let these demands grow excessively, deal.II limits the size $k$ of the
+basis to 30 vectors by default.
+Then, the basis is rebuilt. This implementation of the GMRES method is called
+GMRES(k), with default $k=30$. What we have gained by this restriction, 
+namely a bound on operations and memory requirements, will be compensated by
+the fact that we use an incomplete basis - this will increase the number of
+required iterations. 
+
+BiCGStab, on the other hand, won't get slower when many iterations are needed
+(one iteration uses only results from one preceeding step and
+not all the steps as GMRES). Besides the fact the BiCGStab is more expensive per
+step since two matrix-vector products are needed (compared to one for
+CG or GMRES), there is one main reason which makes BiCGStab not appropriate for
+this problem: The preconditioner applies the inverse of the pressure
+mass matrix by using the InverseMatrix class. Since the application of the
+inverse matrix to a vector is done only in approximative way (an exact inverse
+is too expensive), this will also affect the solver. In the case of BiCGStab, 
+the Krylov vectors will not be orthogonal due to that perturbation. While
+this is uncritical for a small number of steps (up to about 50), it ruins the
+performance of the solver when these perturbations have grown to a significant
+magnitude in the coarse of iterations.
+
+We did some experiments with BiCGStab and found it to
+be faster than GMRES up to refinement cycle 3 (in 3D), but it became very slow 
+for cycles 4 and 5 (even slower than the original Schur complement), so the
+solver is useless in this situation. Choosing a sharper tolerance for the
+inverse matrix class (<code>1e-10*src.l2_norm()</code> instead of 
+<code>1e-6*src.l2_norm()</code>) made BiCGStab perform well also for cycle 4, 
+but did not change the failure on the very large problems.
+
+GMRES is of course also effected by the approximate inverses, but it is not as
+sensitive to orthogonality and retains a relatively good performance also for
+large sizes, see the results below.
+
+With this said, we turn to the realization of the solver call with GMRES with 
+$k=100$ temporary vectors:
+
+@code
+      SparseMatrix<double> pressure_mass_matrix;
+      pressure_mass_matrix.reinit(sparsity_pattern.block(1,1));
+      pressure_mass_matrix.copy_from(system_matrix.block(1,1));
+      system_matrix.block(1,1) = 0;
+
+      SparseILU<double> pmass_preconditioner;
+      pmass_preconditioner.initialize (pressure_mass_matrix, 
+        SparseILU<double>::AdditionalData());
+      
+      InverseMatrix<SparseMatrix<double>,SparseILU<double> >
+        m_inverse (pressure_mass_matrix, pmass_preconditioner);
+
+      BlockSchurPreconditioner<typename InnerPreconditioner<dim>::type,
+                               SparseILU<double> > 
+        preconditioner (system_matrix, m_inverse, *A_preconditioner);
+      
+      SolverControl solver_control (system_matrix.m(),
+                                    1e-6*system_rhs.l2_norm());
+      GrowingVectorMemory<BlockVector<double> > vector_memory;
+      SolverGMRES<BlockVector<double> >::AdditionalData gmres_data;
+      gmres_data.max_n_tmp_vectors = 100;
+      
+      SolverGMRES<BlockVector<double> > gmres(solver_control, vector_memory,
+                                              gmres_data);
+      
+      gmres.solve(system_matrix, solution, system_rhs,
+                  preconditioner);
+      
+      constraints.distribute (solution);
+      
+      std::cout << " "
+                << solver_control.last_step()
+                << " block GMRES iterations";
+@endcode
+
+Obviously, one needs to add the include file @ref SolverGMRES 
+"<lac/solver_gmres.h>" in order to make this run.
+We call the solver with a BlockVector template in order to enable
+GMRES to operate on block vectors and matrices.
+Note also that we need to set the (1,1) block in the system
+matrix to zero (we saved the pressure mass matrix there which is not part of the
+problem) after we copied the information to another matrix.
+
+Using the Timer class, we collect some statistics that compare the runtime
+of the block solver with the one from the problem implementation above.
+Besides the solution with the two options we also check if the solutions
+of the two variants are close to each other (i.e. this solver gives indeed the
+same solution as we had before) and calculate the infinity
+norm of the vector difference.
+
+Let's first see the results in 2D:
+@code
+Refinement cycle 0
+   Number of active cells: 64
+   Number of degrees of freedom: 679 (594+85) [0.005999 s] 
+   Assembling... [0.002 s]
+   Computing preconditioner... [0.003 s]
+   Solving...
+      Schur complement:  11 outer CG iterations for p  [0.007999 s]
+      Block Schur preconditioner:  12 GMRES iterations [0.008998 s]
+   difference l_infty between solution vectors: 8.18909e-07
+
+Refinement cycle 1
+   Number of active cells: 160
+   Number of degrees of freedom: 1683 (1482+201) [0.013998 s] 
+   Assembling... [0.005999 s]
+   Computing preconditioner... [0.012998 s]
+   Solving...
+      Schur complement:  11 outer CG iterations for p  [0.029995 s]
+      Block Schur preconditioner:  12 GMRES iterations [0.030995 s]
+   difference l_infty between solution vectors: 9.32504e-06
+
+Refinement cycle 2
+   Number of active cells: 376
+   Number of degrees of freedom: 3813 (3370+443) [0.031995 s] 
+   Assembling... [0.014998 s]
+   Computing preconditioner... [0.044994 s]
+   Solving...
+      Schur complement:  11 outer CG iterations for p  [0.079987 s]
+      Block Schur preconditioner:  13 GMRES iterations [0.092986 s]
+   difference l_infty between solution vectors: 5.40689e-06
+
+Refinement cycle 3
+   Number of active cells: 880
+   Number of degrees of freedom: 8723 (7722+1001) [0.074988 s] 
+   Assembling... [0.035995 s]
+   Computing preconditioner... [0.110983 s]
+   Solving...
+      Schur complement:  11 outer CG iterations for p  [0.19697 s]
+      Block Schur preconditioner:  13 GMRES iterations [0.242963 s]
+   difference l_infty between solution vectors: 1.14676e-05
+
+Refinement cycle 4
+   Number of active cells: 2008
+   Number of degrees of freedom: 19383 (17186+2197) [0.180973 s] 
+   Assembling... [0.081987 s]
+   Computing preconditioner... [0.315952 s]
+   Solving...
+      Schur complement:  11 outer CG iterations for p  [0.673898 s]
+      Block Schur preconditioner:  13 GMRES iterations [0.778882 s]
+   difference l_infty between solution vectors: 3.13142e-05
+
+Refinement cycle 5
+   Number of active cells: 4288
+   Number of degrees of freedom: 40855 (36250+4605) [0.386941 s] 
+   Assembling... [0.171974 s]
+   Computing preconditioner... [0.766883 s]
+   Solving...
+      Schur complement:  11 outer CG iterations for p  [1.65275 s]
+      Block Schur preconditioner:  13 GMRES iterations [1.81372 s]
+   difference l_infty between solution vectors: 8.59668e-05
+@endcode
+
+We see that there is no huge difference in the solution time between the
+block Schur complement preconditioner solver and the Schur complement
+itself. The reason is simple: we used a direct solve as preconditioner for
+$A$ - so we cannot expect any gain by avoiding the inner iterations. We see
+that the number of iterations has slightly increased for GMRES, but all in
+all the two choices are fairly similar.
+
+The picture of course changes in 3D:
+
+@code
+Refinement cycle 0
+   Number of active cells: 32
+   Number of degrees of freedom: 1356 (1275+81) [0.025996 s] 
+   Assembling... [0.056992 s]
+   Computing preconditioner... [0.027995 s]
+   Solving...
+      Schur complement:  13 outer CG iterations for p  [0.275958 s]
+      Block Schur preconditioner:  23 GMRES iterations [0.042994 s]
+   difference l_infty between solution vectors: 1.11307e-05
+
+Refinement cycle 1
+   Number of active cells: 144
+   Number of degrees of freedom: 5088 (4827+261) [0.102984 s] 
+   Assembling... [0.254961 s]
+   Computing preconditioner... [0.161976 s]
+   Solving...
+      Schur complement:  14 outer CG iterations for p  [2.43963 s]
+      Block Schur preconditioner:  42 GMRES iterations [0.352946 s]
+   difference l_infty between solution vectors: 9.07409e-06
+
+Refinement cycle 2
+   Number of active cells: 704
+   Number of degrees of freedom: 22406 (21351+1055) [0.52592 s] 
+   Assembling... [1.24481 s]
+   Computing preconditioner... [0.948856 s]
+   Solving...
+      Schur complement:  14 outer CG iterations for p  [22.2056 s]
+      Block Schur preconditioner:  78 GMRES iterations [4.75928 s]
+   difference l_infty between solution vectors: 2.48042e-05
+
+Refinement cycle 3
+   Number of active cells: 3168
+   Number of degrees of freedom: 93176 (89043+4133) [2.66759 s] 
+   Assembling... [5.66014 s]
+   Computing preconditioner... [4.69529 s]
+   Solving...
+      Schur complement:  15 outer CG iterations for p  [235.74 s]
+      Block Schur preconditioner:  162 GMRES iterations [63.7883 s]
+   difference l_infty between solution vectors: 5.62978e-05
+
+Refinement cycle 4
+   Number of active cells: 11456
+   Number of degrees of freedom: 327808 (313659+14149) [12.0242 s] 
+   Assembling... [20.2669 s]
+   Computing preconditioner... [17.3384 s]
+   Solving...
+      Schur complement:  15 outer CG iterations for p  [817.287 s]
+      Block Schur preconditioner:  294 GMRES iterations [347.307 s]
+   difference l_infty between solution vectors: 0.000107536
+
+Refinement cycle 5
+   Number of active cells: 45056
+   Number of degrees of freedom: 1254464 (1201371+53093) [89.8533 s] 
+   Assembling... [80.3588 s]
+   Computing preconditioner... [73.0849 s]
+   Solving...
+      Schur complement:  14 outer CG iterations for p  [4401.66 s]
+      Block Schur preconditioner:  587 GMRES iterations [3083.21 s]
+   difference l_infty between solution vectors: 0.00025531
+@endcode
+
+Here, the block preconditioned solver is clearly superior to the Schur
+complement, but the advantage gets less for more mesh points. This is
+because GMRES(k) scales worse with the problem size than CG, as we discussed
+above.  Nonetheless, the improvement by a factor of 3-5 for moderate problem
+sizes is quite impressive.
+
+<h5>Combining block preconditioner and multigrid</h5> 
+An ultimate linear solver for this problem could be imagined as a 
+combination of an optimal
+preconditioner for $A$ (e.g. multigrid) and the block preconditioner
+described above, which is the approach taken in the @ref step_31 "step-31"
+tutorial program.
+
+<h5>No block matrices and vectors</h5>
+Another possibility that can be taken into account is to not set up a block
+system, but rather solve the system of velocity and pressure all at once. The
+options are direct solve with UMFPACK (2D) or GMRES with ILU
+preconditioning (3D). It should be straightforward to try that.
+
+
+
+<h4>More interesting testcases</h4>
+
+The program can of course also serve as a basis to compute the flow in more
+interesting cases. The original motivation to write this program was for it to
+be a starting point for some geophysical flow problems, such as the
+movement of magma under places where continental plates drift apart (for
+example mid-ocean ridges). Of course, in such places, the geometry is more
+complicated than the examples shown above, but it is not hard to accomodate
+for that.
+
+For example, by using the folllowing modification of the boundary values
+function
+@code
+template <int dim>
+double
+BoundaryValues<dim>::value (const Point<dim>  &p,
+                            const unsigned int component) const
+{
+  Assert (component < this->n_components,
+          ExcIndexRange (component, 0, this->n_components));
+
+  const double x_offset = std::atan(p[1]*4)/3;
+                                                                   
+  if (component == 0)
+    return (p[0] < x_offset ? -1 : (p[0] > x_offset ? 1 : 0));
+  return 0;
+}
+@endcode
+and the following way to generate the mesh as the domain
+$[-2,2]\times[-2,2]\times[-1,0]$ 
+@code
+    std::vector<unsigned int> subdivisions (dim, 1);
+    subdivisions[0] = 4;
+    if (dim>2)
+      subdivisions[1] = 4;
+
+    const Point<dim> bottom_left = (dim == 2 ?
+                                    Point<dim>(-2,-1) :
+                                    Point<dim>(-2,-2,-1));
+    const Point<dim> top_right   = (dim == 2 ?
+                                    Point<dim>(2,0) :
+                                    Point<dim>(2,2,0));
+    
+    GridGenerator::subdivided_hyper_rectangle (triangulation,
+                                               subdivisions,
+                                               bottom_left,
+                                               top_right);
+@endcode
+then we get images where the the fault line is curved:
+<TABLE WIDTH="60%" ALIGN="center">
+  <tr>
+    <td ALIGN="center">
+      @image html step-22.3d-extension.png
+    </td>
+
+    <td ALIGN="center">
+      @image html step-22.3d-grid-extension.png
+    </td>
+  </tr>
+</table>
+
diff --git a/deal.II/examples/step-42/step-22.cc b/deal.II/examples/step-42/step-22.cc
new file mode 100644 (file)
index 0000000..5eeadb5
--- /dev/null
@@ -0,0 +1,776 @@
+/* $Id$ */
+/* Author: Wolfgang Bangerth, Texas A&M University, 2008 */
+
+/*    $Id$       */
+/*                                                                */
+/*    Copyright (C) 2008, 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.                        */
+
+
+                                // @sect3{Include files}
+                        
+                                // As usual, we start by including 
+                                // some well-known files:
+#include <base/quadrature_lib.h>
+#include <base/logstream.h>
+#include <base/function.h>
+#include <base/utilities.h>
+
+#include <lac/block_vector.h>
+#include <lac/full_matrix.h>
+#include <lac/block_sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/solver_gmres.h>
+#include <lac/precondition.h>
+#include <lac/constraint_matrix.h>
+
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <grid/grid_tools.h>
+#include <grid/grid_refinement.h>
+
+#include <dofs/dof_handler.h>
+#include <dofs/dof_renumbering.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+
+#include <fe/fe_q.h>
+#include <fe/fe_system.h>
+#include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <numerics/error_estimator.h>
+
+                               
+#include <lac/sparse_direct.h>
+
+#include <lac/sparse_ilu.h>
+
+                               
+#include <fstream>
+#include <sstream>
+
+                               
+using namespace dealii;
+
+                       
+template <int dim>
+struct InnerPreconditioner;
+
+                       
+template <>
+struct InnerPreconditioner<2> 
+{
+    typedef SparseDirectUMFPACK type;
+};
+
+                               
+template <>
+struct InnerPreconditioner<3> 
+{
+    typedef SparseILU<double> type;
+};
+
+
+template <int dim>
+class StokesProblem 
+{
+  public:
+    StokesProblem (const unsigned int degree);
+    void run ();
+    
+  private:
+    void setup_dofs ();
+    void assemble_system ();
+    void solve_schur ();
+    void solve_block ();
+
+    void output_results (const unsigned int refinement_cycle) const;
+    void refine_mesh ();
+    
+    const unsigned int   degree;
+    
+    Triangulation<dim>   triangulation;
+    FESystem<dim>        fe;
+    DoFHandler<dim>      dof_handler;
+
+    ConstraintMatrix     constraints;
+    
+    BlockSparsityPattern      sparsity_pattern;
+    BlockSparseMatrix<double> system_matrix;
+
+    BlockVector<double> solution;
+    BlockVector<double> system_rhs;
+
+                                 
+    std_cxx1x::shared_ptr<typename InnerPreconditioner<dim>::type> A_preconditioner;
+};
+
+                       
+template <int dim>
+class BoundaryValues : public Function<dim> 
+{
+  public:
+    BoundaryValues () : Function<dim>(dim+1) {}
+    
+    virtual double value (const Point<dim>   &p,
+                          const unsigned int  component = 0) const;
+
+    virtual void vector_value (const Point<dim> &p, 
+                               Vector<double>   &value) const;
+};
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value (const Point<dim>  &p,
+                           const unsigned int component) const 
+{
+  Assert (component < this->n_components,
+         ExcIndexRange (component, 0, this->n_components));
+  
+   if (component == 0 && p[0] == 0)
+    return (dim == 2 ? - p[1]*(p[1]-1.) : p[1]*(p[1]-1.) * p[2]*(p[2]-1.));
+  return 0;
+}
+
+
+template <int dim>
+void
+BoundaryValues<dim>::vector_value (const Point<dim> &p,
+                                  Vector<double>   &values) const 
+{
+  for (unsigned int c=0; c<this->n_components; ++c)
+    values(c) = BoundaryValues<dim>::value (p, c);
+}
+
+
+
+                       
+template <int dim>
+class RightHandSide : public Function<dim> 
+{
+  public:
+    RightHandSide () : Function<dim>(dim+1) {}
+    
+    virtual double value (const Point<dim>   &p,
+                          const unsigned int  component = 0) const;
+
+    virtual void vector_value (const Point<dim> &p, 
+                               Vector<double>   &value) const;
+    
+};
+
+
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim>  &/*p*/,
+                           const unsigned int /*component*/) const 
+{
+  return 0;
+}
+
+
+template <int dim>
+void
+RightHandSide<dim>::vector_value (const Point<dim> &p,
+                                  Vector<double>   &values) const 
+{
+  for (unsigned int c=0; c<this->n_components; ++c)
+    values(c) = RightHandSide<dim>::value (p, c);
+}
+
+
+                       
+                               
+template <class Matrix, class Preconditioner>
+class InverseMatrix : public Subscriptor
+{
+  public:
+    InverseMatrix (const Matrix         &m,
+                   const Preconditioner &preconditioner);
+
+    void vmult (Vector<double>       &dst,
+                const Vector<double> &src) const;
+
+  private:
+    const SmartPointer<const Matrix> matrix;
+    const SmartPointer<const Preconditioner> preconditioner;
+};
+
+
+template <class Matrix, class Preconditioner>
+InverseMatrix<Matrix,Preconditioner>::InverseMatrix (const Matrix &m,
+                                                    const Preconditioner &preconditioner)
+               :
+               matrix (&m),
+               preconditioner (&preconditioner)
+{}
+
+
+template <class Matrix, class Preconditioner>
+void InverseMatrix<Matrix,Preconditioner>::vmult (Vector<double>       &dst,
+                                                 const Vector<double> &src) const
+{
+  SolverControl solver_control (src.size(), 1.0e-12);
+  SolverCG<>    cg (solver_control);
+
+  dst = 0;
+
+  cg.solve (*matrix, dst, src, *preconditioner);
+
+  std::cout << "CG steps: " << solver_control.last_step() << std::endl;
+}
+
+template <class PreconditionerA, class PreconditionerMp>
+class BlockSchurPreconditioner : public Subscriptor
+{
+  public:
+    BlockSchurPreconditioner (const BlockSparseMatrix<double>         &S,
+          const InverseMatrix<SparseMatrix<double>,PreconditionerMp>  &Mpinv,
+          const PreconditionerA &Apreconditioner);
+
+  void vmult (BlockVector<double>       &dst,
+              const BlockVector<double> &src) const;
+
+  private:
+    const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
+    const SmartPointer<const InverseMatrix<SparseMatrix<double>, 
+                       PreconditionerMp > > m_inverse;
+    const PreconditionerA &a_preconditioner;
+    
+    mutable Vector<double> tmp;
+
+};
+
+template <class PreconditionerA, class PreconditionerMp>
+BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::BlockSchurPreconditioner(
+          const BlockSparseMatrix<double>                            &S,
+          const InverseMatrix<SparseMatrix<double>,PreconditionerMp> &Mpinv,
+          const PreconditionerA &Apreconditioner
+          )
+                :
+                system_matrix           (&S),
+                m_inverse               (&Mpinv),
+                a_preconditioner        (Apreconditioner),
+                tmp                     (S.block(1,1).m())
+{}
+
+        // Now the interesting function, the multiplication of
+        // the preconditioner with a BlockVector. 
+template <class PreconditionerA, class PreconditionerMp>
+void BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::vmult (
+                                     BlockVector<double>       &dst,
+                                     const BlockVector<double> &src) const
+{
+        // Form u_new = A^{-1} u
+  a_preconditioner.vmult (dst.block(0), src.block(0));
+        // Form tmp = - B u_new + p 
+        // (<code>SparseMatrix::residual</code>
+        // does precisely this)
+  system_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1));
+        // Change sign in tmp
+  tmp *= -1;
+        // Multiply by approximate Schur complement 
+        // (i.e. a pressure mass matrix)
+  m_inverse->vmult (dst.block(1), tmp);
+}
+
+
+template <class Preconditioner>
+class SchurComplement : public Subscriptor
+{
+  public:
+    SchurComplement (const BlockSparseMatrix<double> &system_matrix,
+                    const InverseMatrix<SparseMatrix<double>, Preconditioner> &A_inverse);
+
+    void vmult (Vector<double>       &dst,
+               const Vector<double> &src) const;
+
+  private:
+    const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
+    const SmartPointer<const InverseMatrix<SparseMatrix<double>, Preconditioner> > A_inverse;
+    
+    mutable Vector<double> tmp1, tmp2;
+};
+
+
+
+template <class Preconditioner>
+SchurComplement<Preconditioner>::
+SchurComplement (const BlockSparseMatrix<double> &system_matrix,
+                const InverseMatrix<SparseMatrix<double>,Preconditioner> &A_inverse)
+               :
+               system_matrix (&system_matrix),
+               A_inverse (&A_inverse),
+               tmp1 (system_matrix.block(0,0).m()),
+               tmp2 (system_matrix.block(0,0).m())
+{}
+
+
+template <class Preconditioner>
+void SchurComplement<Preconditioner>::vmult (Vector<double>       &dst,
+                                            const Vector<double> &src) const
+{
+  system_matrix->block(0,1).vmult (tmp1, src);
+  A_inverse->vmult (tmp2, tmp1);
+  system_matrix->block(1,0).vmult (dst, tmp2);
+}
+
+
+                       
+template <int dim>
+StokesProblem<dim>::StokesProblem (const unsigned int degree)
+                :
+                degree (degree),
+                triangulation (Triangulation<dim>::maximum_smoothing),
+                fe (FE_Q<dim>(degree+1), dim,
+                    FE_Q<dim>(degree), 1),
+                dof_handler (triangulation)
+{}
+
+
+                               
+                       
+                               
+template <int dim>
+void StokesProblem<dim>::setup_dofs ()
+{
+  A_preconditioner.reset ();
+  system_matrix.clear ();
+  
+  dof_handler.distribute_dofs (fe);  
+  DoFRenumbering::Cuthill_McKee (dof_handler);
+
+  std::vector<unsigned int> block_component (dim+1,0);
+  block_component[dim] = 1;
+  DoFRenumbering::component_wise (dof_handler, block_component);
+
+                               
+  {
+    constraints.clear ();
+    std::vector<bool> component_mask (dim+1, true);
+    component_mask[dim] = false;
+    VectorTools::interpolate_boundary_values (dof_handler,
+                                             0,
+                                             BoundaryValues<dim>(),
+                                             constraints,
+                                             component_mask);
+    DoFTools::make_hanging_node_constraints (dof_handler,
+                                            constraints);
+  }
+
+  constraints.close ();
+
+                               
+  std::vector<unsigned int> dofs_per_block (2);
+  DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);  
+  const unsigned int n_u = dofs_per_block[0],
+                     n_p = dofs_per_block[1];
+
+  std::cout << "   Number of active cells: "
+            << triangulation.n_active_cells()
+            << std::endl
+            << "   Number of degrees of freedom: "
+            << dof_handler.n_dofs()
+            << " (" << n_u << '+' << n_p << ')'
+            << std::endl;
+      
+                               
+                                
+                                
+  {
+    BlockCompressedSimpleSparsityPattern csp (2,2);
+
+    csp.block(0,0).reinit (n_u, n_u);
+    csp.block(1,0).reinit (n_p, n_u);
+    csp.block(0,1).reinit (n_u, n_p);
+    csp.block(1,1).reinit (n_p, n_p);
+  
+    csp.collect_sizes();    
+  
+    DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false);
+    sparsity_pattern.copy_from (csp);
+  }
+  
+                                
+  system_matrix.reinit (sparsity_pattern);
+                                   
+  solution.reinit (2);
+  solution.block(0).reinit (n_u);
+  solution.block(1).reinit (n_p);
+  solution.collect_sizes ();
+  
+  system_rhs.reinit (2);
+  system_rhs.block(0).reinit (n_u);
+  system_rhs.block(1).reinit (n_p);
+  system_rhs.collect_sizes ();
+}
+
+
+                       
+template <int dim>
+void StokesProblem<dim>::assemble_system () 
+{
+  system_matrix=0;
+  system_rhs=0;
+  
+  QGauss<dim>   quadrature_formula(degree+2);
+
+  FEValues<dim> fe_values (fe, quadrature_formula,
+                           update_values    |
+                           update_quadrature_points  |
+                           update_JxW_values |
+                           update_gradients);
+  
+  const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+  
+  const unsigned int   n_q_points      = quadrature_formula.size();
+
+  FullMatrix<double>   local_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       local_rhs (dofs_per_cell);
+
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+  
+  const RightHandSide<dim>          right_hand_side;
+  std::vector<Vector<double> >      rhs_values (n_q_points,
+                                                Vector<double>(dim+1));
+
+                                
+  const FEValuesExtractors::Vector velocities (0);
+  const FEValuesExtractors::Scalar pressure (dim);
+
+                                
+                                 
+  std::vector<Tensor<2,dim> >          phi_grads_u (dofs_per_cell);
+  std::vector<double>                  div_phi_u   (dofs_per_cell);
+  std::vector<double>                  phi_p       (dofs_per_cell);
+                                  
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    { 
+      fe_values.reinit (cell);
+      local_matrix = 0;
+      local_rhs = 0;
+      
+      right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
+                                        rhs_values);
+      
+      for (unsigned int q=0; q<n_q_points; ++q)
+       {
+         for (unsigned int k=0; k<dofs_per_cell; ++k)
+           {
+             phi_grads_u[k] = fe_values[velocities].gradient (k, q);
+             div_phi_u[k]   = fe_values[velocities].divergence (k, q);
+             phi_p[k]       = fe_values[pressure].value (k, q);
+           }
+
+         for (unsigned int i=0; i<dofs_per_cell; ++i)
+           {
+             for (unsigned int j=0; j<dofs_per_cell; ++j)
+               {
+                 local_matrix(i,j) += (scalar_product(phi_grads_u[i], phi_grads_u[j])
+                                       - div_phi_u[i] * phi_p[j]
+                                       - phi_p[i] * div_phi_u[j]
+                                       + phi_p[i] * phi_p[j])
+                                      * fe_values.JxW(q);     
+
+               }
+
+             const unsigned int component_i =
+               fe.system_to_component_index(i).first;
+             local_rhs(i) += fe_values.shape_value(i,q) * 
+                             rhs_values[q](component_i) *
+                             fe_values.JxW(q);
+           }
+       }
+
+                                     
+
+
+      cell->get_dof_indices (local_dof_indices);
+      constraints.distribute_local_to_global (local_matrix, local_rhs,
+                                             local_dof_indices, 
+                                             system_matrix, system_rhs);
+    }
+  
+                               
+  std::cout << "   Computing preconditioner..." << std::endl << std::flush;
+      
+  A_preconditioner
+    = std_cxx1x::shared_ptr<typename InnerPreconditioner<dim>::type>(new typename InnerPreconditioner<dim>::type());
+  A_preconditioner->initialize (system_matrix.block(0,0),
+                               typename InnerPreconditioner<dim>::type::AdditionalData());
+
+}
+
+
+
+                               
+template <int dim>
+void StokesProblem<dim>::solve_block () 
+{
+ SparseMatrix<double> pressure_mass_matrix;
+      pressure_mass_matrix.reinit(sparsity_pattern.block(1,1));
+      pressure_mass_matrix.copy_from(system_matrix.block(1,1));
+      system_matrix.block(1,1) = 0;
+
+      SparseILU<double> pmass_preconditioner;
+      pmass_preconditioner.initialize (pressure_mass_matrix, 
+        SparseILU<double>::AdditionalData());
+      
+      InverseMatrix<SparseMatrix<double>,SparseILU<double> >
+        m_inverse (pressure_mass_matrix, pmass_preconditioner);
+
+      BlockSchurPreconditioner<typename InnerPreconditioner<dim>::type,
+                               SparseILU<double> > 
+        preconditioner (system_matrix, m_inverse, *A_preconditioner);
+      
+      SolverControl solver_control (system_matrix.m(),
+                                    1e-6*system_rhs.l2_norm());
+      GrowingVectorMemory<BlockVector<double> > vector_memory;
+      SolverGMRES<BlockVector<double> >::AdditionalData gmres_data;
+      gmres_data.max_n_tmp_vectors = 100;
+      
+      SolverGMRES<BlockVector<double> > gmres(solver_control, vector_memory,
+                                              gmres_data);
+      
+      gmres.solve(system_matrix, solution, system_rhs,
+                  preconditioner);
+      
+      constraints.distribute (solution);
+      
+      std::cout << " "
+                << solver_control.last_step()
+                << " block GMRES iterations";
+
+
+
+
+}
+
+template <int dim>
+void StokesProblem<dim>::solve_schur () 
+{
+
+
+
+  const InverseMatrix<SparseMatrix<double>,
+                      typename InnerPreconditioner<dim>::type>
+    A_inverse (system_matrix.block(0,0), *A_preconditioner);
+  Vector<double> tmp (solution.block(0).size());
+  
+                                 
+  {
+    Vector<double> schur_rhs (solution.block(1).size());
+    A_inverse.vmult (tmp, system_rhs.block(0));
+    system_matrix.block(1,0).vmult (schur_rhs, tmp);
+    schur_rhs -= system_rhs.block(1);
+  
+    SchurComplement<typename InnerPreconditioner<dim>::type>
+      schur_complement (system_matrix, A_inverse);
+    
+                                    // The usual control structures for
+                                    // the solver call are created...
+    SolverControl solver_control (solution.block(1).size(),
+                                 1e-6*schur_rhs.l2_norm());
+    SolverCG<>    cg (solver_control);
+    
+                                   
+                                  
+    SparseILU<double> preconditioner;
+    preconditioner.initialize (system_matrix.block(1,1), 
+      SparseILU<double>::AdditionalData());
+  
+    InverseMatrix<SparseMatrix<double>,SparseILU<double> >
+      m_inverse (system_matrix.block(1,1), preconditioner);
+    
+                                   
+    cg.solve (schur_complement, solution.block(1), schur_rhs,
+             m_inverse);
+  
+                                  
+    constraints.distribute (solution);
+  
+    std::cout << "  "
+             << solver_control.last_step()
+             << " outer CG Schur complement iterations for pressure"
+             << std::flush
+             << std::endl;    
+  }
+    
+                                
+  {
+    system_matrix.block(0,1).vmult (tmp, solution.block(1));
+    tmp *= -1;
+    tmp += system_rhs.block(0);
+  
+    A_inverse.vmult (solution.block(0), tmp);
+
+    constraints.distribute (solution);
+  }
+
+
+}
+
+
+                       
+template <int dim>
+void
+StokesProblem<dim>::output_results (const unsigned int refinement_cycle)  const
+{
+  std::vector<std::string> solution_names (dim, "velocity");
+  solution_names.push_back ("pressure");
+  
+  std::vector<DataComponentInterpretation::DataComponentInterpretation>
+    data_component_interpretation
+    (dim, DataComponentInterpretation::component_is_part_of_vector);
+  data_component_interpretation
+    .push_back (DataComponentInterpretation::component_is_scalar);
+      
+  DataOut<dim> data_out;
+  data_out.attach_dof_handler (dof_handler);  
+  data_out.add_data_vector (solution, solution_names,
+                           DataOut<dim>::type_dof_data,
+                           data_component_interpretation);
+  data_out.build_patches ();
+  
+  std::ostringstream filename;
+  filename << "solution-"
+           << Utilities::int_to_string (refinement_cycle, 2)
+           << ".vtk";
+
+  std::ofstream output (filename.str().c_str());
+  data_out.write_vtk (output);
+}
+
+
+                       
+template <int dim>
+void
+StokesProblem<dim>::refine_mesh () 
+{
+  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+  std::vector<bool> component_mask (dim+1, false);
+  component_mask[dim] = true;
+  KellyErrorEstimator<dim>::estimate (dof_handler,
+                                      QGauss<dim-1>(degree+1),
+                                      typename FunctionMap<dim>::type(),
+                                      solution,
+                                      estimated_error_per_cell,
+                                      component_mask);
+
+  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                   estimated_error_per_cell,
+                                                   0.3, 0.0);
+  triangulation.execute_coarsening_and_refinement ();
+}
+
+
+                               
+template <int dim>
+void StokesProblem<dim>::run () 
+{
+  {
+    std::vector<unsigned int> subdivisions (dim, 1);
+    subdivisions[0] = 4;
+
+    const Point<dim> bottom_left = (dim == 2 ?
+                                   Point<dim>(0,0) :
+                                   Point<dim>(0,0,0));
+    const Point<dim> top_right   = (dim == 2 ?
+                                   Point<dim>(1,1) :
+                                   Point<dim>(1,1,1));
+    
+    GridGenerator::subdivided_hyper_rectangle (triangulation,
+                                              subdivisions,
+                                              bottom_left,
+                                              top_right);
+  }
+  
+                                
+  for (typename Triangulation<dim>::active_cell_iterator
+        cell = triangulation.begin_active();
+       cell != triangulation.end(); ++cell)
+    for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+      if (cell->face(f)->center()[0] == 1)
+       cell->face(f)->set_all_boundary_indicators(1);
+  
+  
+                               
+  triangulation.refine_global (4-dim);
+
+                                 
+  for (unsigned int refinement_cycle = 0; refinement_cycle<6;
+       ++refinement_cycle)
+    {
+      std::cout << "Refinement cycle " << refinement_cycle << std::endl;
+      
+      if (refinement_cycle > 0)
+        refine_mesh ();
+      
+      setup_dofs ();
+
+      std::cout << "   Assembling..." << std::endl << std::flush;
+      assemble_system ();      
+
+      std::cout << "   Solving..." << std::flush;
+      solve_block ();
+      
+      output_results (refinement_cycle);
+
+      std::cout << std::endl;
+    }
+}
+
+
+                       
+int main () 
+{
+  try
+    {
+      deallog.depth_console (0);
+
+      StokesProblem<2> flow_problem(1);
+      flow_problem.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.