From: heltai Date: Thu, 29 Aug 2013 10:34:57 +0000 (+0000) Subject: Merged from trunk X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4f080c6b6d159e8ce04c074dde4b32c4e1fe2148;p=dealii-svn.git Merged from trunk git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@30524 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/development/CMakeLists.txt.sample b/deal.II/doc/development/CMakeLists.txt.sample deleted file mode 100644 index 549cf0d0bb..0000000000 --- a/deal.II/doc/development/CMakeLists.txt.sample +++ /dev/null @@ -1,61 +0,0 @@ -## ## -# A simple CMakeLists.txt file for use with deal.II # -## ## - -# -# Set the name of the project and target: -# -SET(TARGET "step-1") - -# -# Declare all source files the target consists of: -# -SET(TARGET_SRC - ${TARGET}.cc - # You can specify additional files here! - ) - -# -# (Optional) -# Specify a list of files (file globs) that will be removed -# with the "make runclean" and "make distclean" targets. -# If empty, sensible default values will be used. -# -SET(CLEAN_UP_FILES - # a custom list of globs, e.g. *.log *.vtk - ) - -# -# (Optional) -# A custom command line that should be invoked by "make run". If empty, -# ./${TARGET} will be invoked. -# -SET(TARGET_RUN - # a custom command line, e.g. mpirun -np 2 ${TARGET} - ) - -# -# It is good practice to specify a version requirement: -# -CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) - -# -# Find and import the deal.II project configuration: -# -FIND_PACKAGE(deal.II 8.0 REQUIRED - HINTS - ${DEAL_II_DIRS} ../ ../../ $ENV{DEAL_II_DIR} - # - # If the deal.II library cannot be found (because it is not installed at a - # default location or your project resides at an uncommon place), you - # can specify additional hints for search paths here, e.g. - # "$ENV{HOME}/workspace/deal.II" - # - ) - -# -# And let deal.II do the rest: -# -DEAL_II_INITIALIZE_CACHED_VARIABLES() -PROJECT(${TARGET}) -DEAL_II_INVOKE_AUTOPILOT() diff --git a/deal.II/doc/development/CMakeLists.txt.sample1 b/deal.II/doc/development/CMakeLists.txt.sample1 new file mode 100644 index 0000000000..fda35efbf6 --- /dev/null +++ b/deal.II/doc/development/CMakeLists.txt.sample1 @@ -0,0 +1,20 @@ +FIND_PACKAGE(deal.II 8.0 REQUIRED + HINTS + ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + # You can specify additional hints for search paths here, e.g. + # $ENV{HOME}/workspace/deal.II + ) + +# Set the name of the project and target: +SET(TARGET "step-1") + +# Declare all source files the target consists of: +SET(TARGET_SRC + step-1.cc + # You can specify additional files here! + ) + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET} CXX) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/deal.II/doc/development/CMakeLists.txt.sample2 b/deal.II/doc/development/CMakeLists.txt.sample2 new file mode 100644 index 0000000000..2003f225e9 --- /dev/null +++ b/deal.II/doc/development/CMakeLists.txt.sample2 @@ -0,0 +1,33 @@ +FIND_PACKAGE(deal.II 8.0 REQUIRED + HINTS + ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + # You can specify additional hints for search paths here, e.g. + # $ENV{HOME}/workspace/deal.II + ) + +# Set the name of the project and target: +SET(TARGET "step-1") + +# Declare all source files the target consists of: +SET(TARGET_SRC + step-1.cc + # You can specify additional files here! + ) + +# Specify a list of files (file globs) that will be removed +# with the "make runclean" and "make distclean" targets. +# (If empty, sensible default values will be used.) +SET(CLEAN_UP_FILES + # a custom list of globs, e.g. *.log *.vtk + ) + +# A custom command line that should be invoked by "make run". +# (If empty, ./${TARGET} will be invoked.) +SET(TARGET_RUN + # a custom command line, e.g. mpirun -np 2 ${TARGET} + ) + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET} CXX) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/deal.II/doc/development/CMakeLists.txt.sample3 b/deal.II/doc/development/CMakeLists.txt.sample3 new file mode 100644 index 0000000000..7d484c182f --- /dev/null +++ b/deal.II/doc/development/CMakeLists.txt.sample3 @@ -0,0 +1,15 @@ +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) + +FIND_PACKAGE(deal.II 8.0 REQUIRED + HINTS ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(step-1 CXX) + +ADD_EXECUTABLE(step-1 step-1.cc) +DEAL_II_SETUP_TARGET(step-1) + +ADD_CUSTOM_TARGET(run COMMAND step-1 + COMMENT "Run with ${CMAKE_BUILD_TYPE} configuration" + ) diff --git a/deal.II/doc/development/cmake-internals.html b/deal.II/doc/development/cmake-internals.html index 903af06bfa..fbfdbab628 100644 --- a/deal.II/doc/development/cmake-internals.html +++ b/deal.II/doc/development/cmake-internals.html @@ -11,157 +11,155 @@ - - - -

Build system internals

- -

- This page provides details about the CMake build system. Files - processed by the top level CMakeLists.txt script are - listed in the TOC in chronological order. -

- -

- - - - -
Table of contents
- -
Configuration: - -
Target definition and installation - -
-

- - -

Coding convention

-

- Coding conventions are always a matter of choice. Nevertheless, the - following rules should be considered: -

- An example: - 
-
-    FOREACH(_build ${DEAL_II_BUILD_TYPES})
-      #
-      # Set an appropriate keyword depending on target and build type:
-      #
-      IF(NOT "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease")
-        SET(_keyword "general")
-      ELSE()
-        IF(_build MATCHES DEBUG)
-          SET(_keyword "debug")
-        ELSE()
-          SET(_keyword "optimized")
-        ENDIF()
-      ENDIF()
-    ENDFOREACH()
-
- - 
-
-    LIST(APPEND CONFIG_LIBRARIES
-      ${_keyword}
-      ${CONFIG_LIBRARIES_${_build}}
-      )
-
-    SET_TARGET_PROPERTIES(${DEAL_II_BASE_NAME}${DEAL_II_${build}_SUFFIX}
-      PROPERTIES
-      VERSION ${VERSION}
-      SOVERSION ${VERSION}
-      LINK_FLAGS "${DEAL_II_LINKER_FLAGS_${build}}"
-      COMPILE_DEFINITIONS "${DEAL_II_DEFINITIONS};${DEAL_II_DEFINITIONS_${build}}"
-      COMPILE_FLAGS "${DEAL_II_CXX_FLAGS_${build}}"
-      )
-
- CMake operates almost always with variables in global state. To guard - against accidental overwrite of variables the following naming - conventions must be followed at all times: - -

- - - -

./CMakeLists.txt and ./cmake/setup_*.cmake

- -

- The very first configuration steps after some initial setup in - ./CMakeLists.txt takes place in some - ./cmake/setup_*.cmake files: -

+

- -

./cmake/checks/check_*.cmake

+ +

./cmake/checks/check_*.cmake

-

- The next step in the configuration process is to include all - checks residing under ./cmake/checks. Currently - there are: - 

+

+  The next step in the configuration process is to include all
+  checks residing under ./cmake/checks. Currently
+  there are:
 
+
+    
     ./cmake/checks/check_01_compiler_features.cmake
       - Search for support for compiler dependent features such as stack
         trace support, demangler support, etc.
@@ -240,24 +239,27 @@
 
     ./cmake/checks/check_02_compiler_bugs.cmake
       - Check for compiler bugs
-      

-
-
  • Necessary include directories and libraries necessary for - linkage can be set in the list variables - CMAKE_REQUIRED_INCLUDES and - CMAKE_REQUIRED_LIBRARIES. It is best to append these - lists and later on reset CMAKE_REQUIRED_* (including - CMAKE_REQUIRED_FLAGS) to their default values: - 
    +  
  •  Necessary include directories and libraries necessary for
+    linkage can be set in the list variables
+    CMAKE_REQUIRED_INCLUDES and
+    CMAKE_REQUIRED_LIBRARIES. It is best to append these
+    lists and later on reset CMAKE_REQUIRED_* (including
+    CMAKE_REQUIRED_FLAGS) to their default values:
+    
+
     
     LIST(APPEND CMAKE_REQUIRED_INCLUDES <a list of includes>)
     LIST(APPEND CMAKE_REQUIRED_LIBRARIES <a list of libraries>)
     CHECK_CXX_SOURCE_COMPILES(...)
     RESET_CMAKE_REQUIRED()
-	  
    
-      
-    
    
-
-
-
-    
-    
     ./cmake/modules/Find*.cmake 
    
-
-    
    
-      These are find modules for the configure_*.cmake files
-      and the CONFIGURE_FEATURE macro as will explained later.
-      It is crucial that a find module behaves correctly. Therefore, the
-      following rules are mandatory:
-      
      
-        
    • 
-          The sole purpose of a find module is to find an external
-          library (no deal.II specific dependency checking, no
-          compatibility checking).
-        
    • 
-          It should do so by appropriate FIND_LIBRARY and
-          FIND_PATH calls. The results of this calls should be the
-          only cached variables.
-        
    • 
-          FIND_PACKAGE_HANDLE_STANDARD_ARGS should be used for
-          setting FEATURE_FOUND as this already respects the
-          QUIET and REQUIRED keywords.
-        
    • 
-          A WARNING, SEND_ERROR or
-          FATAL_ERROR must be avoided (the only exception is
-          the REQUIRED keyword).
-        
    • 
-          If successful, the following uncached variables might be set:
-          
      -
+
      
+
    
+
    
+
+
+
     ./cmake/modules/Find*.cmake 
    
+
+
    
+  These are find modules for the configure_*.cmake files
+  and the CONFIGURE_FEATURE macro as will explained later.
+  It is crucial that a find module behaves correctly. Therefore, the
+  following rules are mandatory:
+  
      
+    
    • 
+      The sole purpose of a find module is to find an external
+      library (no deal.II specific dependency checking, no
+      compatibility checking).
+    
    • 
+      It should do so by appropriate FIND_LIBRARY and
+      FIND_PATH calls. The results of this calls should be the
+      only cached variables.
+    
    • 
+      FIND_PACKAGE_HANDLE_STANDARD_ARGS should be used for
+      setting FEATURE_FOUND as this already respects the
+      QUIET and REQUIRED keywords.
+    
    • 
+      A WARNING, SEND_ERROR or
+      FATAL_ERROR must be avoided (the only exception is
+      the REQUIRED keyword).
+    
    • 
+      If successful, the following uncached variables might be set:
+
+
      +	
     FEATURE_FOUND
     FEATURE_LIBRARIES
     FEATURE_INCLUDE_DIRS
@@ -359,41 +363,43 @@
     FEATURE_VERSION_MAJOR
     FEATURE_VERSION_MINOR
     FEATURE_VERSION_SUBMINOR
-          
      
-          There are obviously valid exceptions from this rule, though.
-        
    • 
-          Only "local" variables "_<all lowercase>" or
-          "global" variables prefixed by FEATURE_ may be
-          altered. Do not set DEAL_II_* or CMAKE_*
-          variables directly!
-        
    • 
-          A hint with FEATURE_DIR can be set up for
-          convenience. It is best to start the Find module by
-          
      -
-    SET_IF_EMPTY(FEATURE_DIR "$ENV{FEATURE_DIR}")
-          
      
-          and use FEATURE_DIR as a hint. If the external
-          library could be found, hide it via
-          MARK_AS_ADVANCED, if the external library could not
-          be found, set a cache value:
-          
      -
-    SET(FEATURE_DIR "" CACHE PATH "An optional hint to a FEATURE directory")
-          
      
-      
    
+
    • +There are obviously valid exceptions from this rule, though. +
  • + Only "local" variables "_<all lowercase>" or + "global" variables prefixed by FEATURE_ may be + altered. Do not set DEAL_II_* or CMAKE_* + variables directly! +
  • + A hint with FEATURE_DIR can be set up for + convenience. It is best to start the Find module by +
     
-    
-    
     ./cmake/configure/configure_*.cmake 
    
+    SET_IF_EMPTY(FEATURE_DIR "$ENV{FEATURE_DIR}")
+
    + +and use FEATURE_DIR as a hint. If the external +library could be found, hide it via +MARK_AS_ADVANCED, if the external library could not +be found, set a cache value: +
    +  
+    SET(FEATURE_DIR "" CACHE PATH "An optional hint to a FEATURE directory")
+
    + -

    - The final step in the configuration phase is the setup of features - (which refer to external or bundled libraries - deal.II can optionally interface with.) -

    + +

    ./cmake/configure/configure_*.cmake

    + +

    + The final step in the configuration phase is the setup of features + (which refer to external or bundled libraries + deal.II can optionally interface with.) + +

    At bare minimum configure_<feature>.cmake @@ -574,7 +580,8 @@

    - + +

    Target definition and installation

    ./include/deal.II/base/config.h.in

    diff --git a/deal.II/doc/development/cmake.html b/deal.II/doc/development/cmake.html index 8139fdc110..2683cc1d7e 100644 --- a/deal.II/doc/development/cmake.html +++ b/deal.II/doc/development/cmake.html @@ -24,38 +24,34 @@ deal.II CMake build system.

    -

    - - - - -
    Table of contents
    - -
    -

    - +
    +
      +
    1. Operating cmake +
        +
      1. Manipulating the cache
        2. +
      2. ccmake and special build targets
        3. +
      3. Shortcuts
        4. +
      +
      2. +
    2. Configuration options +
        +
      1. Feature configuration
        2. +
      2. Autoconfiguration
        3. +
      3. External library locations
        4. +
      4. Manual override
        5. +
      5. Component selection
        6. +
      6. Build configuration
        7. +
      7. Selecting a compiler
        8. +
      8. Installation
        9. +
      +
      3. +
    3. Initial cache file and advanced options
      4. +
    4. Compiling only certain parts
      5. +
    +
    + -

    Operating cmake

    +

    Operating cmake

    When configuring deal.II by @@ -129,7 +125,7 @@

    -

    Manipulating the cache

    +

    Manipulating the cache

    A cached variable can be set on the command line via @@ -153,7 +149,7 @@

    -

    ccmake and special build targets

    +

    ccmake and special build targets

    A very convenient way to alter the configuration is to use the graphical @@ -177,7 +173,7 @@

    -

    Shortcuts

    +

    Shortcuts

    All variables starting with WITH_ will be automatically @@ -199,7 +195,7 @@

    -

    Configuration options

    +

    Configuration options

    The various configuration options of the @@ -212,7 +208,7 @@ -

    Feature configuration

    +

    Feature configuration

    deal.II provides (optional) interfaces to quite @@ -304,7 +300,7 @@ -

    Autoconfiguration

    +

    Autoconfiguration

    As long as DEAL_II_WITH_<FEATURE> is @@ -339,7 +335,7 @@ -

    External library locations

    +

    External library locations

    External libraries will be searched depending on hints in the following @@ -416,7 +412,7 @@

    -

    Manual override

    +

    Manual override

    Warning: Do not do this unless absolutely necessary!

    It is possible to override the CMake find mechanism for external @@ -462,7 +458,7 @@ -

    Component selection

    +

    Component selection

    The following options control which components of @@ -532,7 +528,7 @@ -

    Build configuration

    +

    Build configuration

    The cmake variable CMAKE_BUILD_TYPE @@ -595,7 +591,7 @@ -

    Selecting a compiler

    +

    Selecting a compiler

    Compilers can be switched either by command line or by setting @@ -659,7 +655,7 @@ -

    Installation

    +

    Installation

    the location, @@ -726,7 +722,7 @@ -

    Initial cache file and advanced options

    +

    Initial cache file and advanced options

    A sample configuration file for preloading the CMake cache with 
    @@ -740,7 +736,7 @@
 
 
     
-    
     Compiling only certain parts 
    
+    
     Compiling only certain parts 
    
     
    
 
     
    
diff --git a/deal.II/doc/development/cmakelists.html b/deal.II/doc/development/cmakelists.html
index 96f699ef46..44fb38ae3f 100644
--- a/deal.II/doc/development/cmakelists.html
+++ b/deal.II/doc/development/cmakelists.html
@@ -2,7 +2,7 @@
 	  "http://www.w3.org/TR/html4/loose.dtd">
 
   
-    How to interface with deal.II in your own project
+    How to use CMake to configure your projects with deal.II
     
     
     
@@ -10,469 +10,374 @@
     
     
   
+  
+
+
    How to use CMake to configure your projects with deal.II
    
 
-  
-
-    
    How to interface with deal.II in your own
-      project
    
-
-    
    
-    
-      
-      
-      
-    
    Table of contents
    
-	  
-	
    
-    
    
-
-
-    
-    
     Example CMakeLists.txt files 
    
-
-    
    
-      cmake is the configuration and build tool we use
-      in deal.II. Its advantage is not only that it makes
-      configuration of deal.II itself simpler across
-      platforms (compared to the older autoconf/make combination)
-      but also that it exports information about the deal.II
-      configuration that makes it particularly simple for projects
-      using deal.II to configure and link against it.
-    
    
+
    
+  
      
+    
    1. Simple CMakeLists.txt
+      
        
+	
      1. Adding multiple executable targets
        2. 
+	
      2. Adding libraries and common source files
        3. 
+	
      3. Adding a "run" target
        4. 
+      
      
+    
      2. 
+    
    2. Advanced CMakeLists.txt
+      
        
+	
      1.  TODO, see list
        2. 
+      
      
+    
      3. 
+    
    3. Autopilot style CMakeLists.txt
      4. 
+    
    4. deal.IIConfig.cmake
      5. 
+    
    5. Legacy Make.global_options
      6. 
+  
    
+
    
+
+
    
+  cmake is the configuration and build tool we use
+  in deal.II. Its advantage is not only that it makes
+  configuration of deal.II itself simpler across
+  platforms (compared to the older autoconf/make combination)
+  but also that it exports information about the deal.II
+  configuration that makes it particularly simple for projects
+  using deal.II to configure and link against it.
+
    
+
+
    
+  cmake is controlled by input files that by convention are
+  called CMakeLists.txt, listing both configuration commands
+  as well as dependencies between source files and targets.
+  This page presents some CMakeLists.txt examples for
+  potential use in your projects. (A detailed description of the
+  deal.II project configuration is given in the
+  deal.II CMake ReadMe.)
+
    
+
+
+
    Simple CMakeLists.txt
    
+
+
    
+  In this section, we start out with a
+  minimal CMakeLists.txt based on
+  the DEAL_II_SETUP_TARGET macro. This method gives
+  full control of what's happening and is easily extensible to
+  more complex projects, as exemplified in the subsections here an
+  later in the section on advanced
+    topics. Here is a full example
+  (plain text
+  version):
+  
+
     Fix plain text after finalizing!
    
+
+
    +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
+
+FIND_PACKAGE(deal.II 8.0 REQUIRED
+  HINTS ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+)
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+
+PROJECT(myproject CXX)
+
+ADD_EXECUTABLE(mycode mycode.cc)
+DEAL_II_SETUP_TARGET(mycode)
+
    
+
+
    
+  The first line of this code makes sure that a sufficiently high
+  version of CMake is installed. 2.8.8 is the minimal version required
+  to set up deal.II, therefore it is safe to use this
+  number here.
+
    
+
+
    
+  Next, we find our deal.II installation, in this case requiring at
+  least version 8.0, which is the first version using
+  CMake. The HINTS are a list of directories where the
+  install directory of deal.II is likely to be found. Here, we check
+  whether we are in a subdirectory (first and second level) of the
+  deal.II installation and otherwise use the
+  variable DEAL_II_DIR.
+  Matthias, what is the difference between the first
+  and the last entry?  This list can be changed according to your
+  preferences. After finding the deal.II project, we fetch its chached
+  variables. You can inspect these for instance
+  with ccmake.
+
    
+
+
    
+  Every CMAkeLists.txt must contain a project definition,
+  which we do next. Matthias, I have been happy
+  without the CXX...?
+
    
+
+
    
+  Finally, the last two lines define the executable that is to be
+  produced and its source code.
+
    
+
+
+
     Explain the macros or link to explanation below
    
+
+
+
    Adding multiple executable targets
    
+
+
    
+  In order to specify multiple executable targets, simply repeat
+  the last two lines of the simple CMakeLists.txt:
+
    
+
+
    +ADD_EXECUTABLE(mycode2 mycode2.cc)
+DEAL_II_SETUP_TARGET(mycode2)
+
+ADD_EXECUTABLE(mycode3 mycode3.cc)
+DEAL_II_SETUP_TARGET(mycode3)
+
    
+
+If the list gets longer, consider using
+a loop, possibly
+with GLOB.
+
+    
+    
    Adding libraries and common source files
    
 
     
    
-      cmake is controlled by input files that by convention are
-      called CMakeLists.txt, listing both configuration commands
-      as well as dependencies between source files and targets.
-      This page presents some CMakeLists.txt examples for
-      potential use in your projects. A detailed description of the
-      deal.II project configuration is given in the
-      deal.II CMake ReadMe.
+      Adding a library is as simple as adding an executable target. We
+      specify the library name and then have to tell cmake that the
+      executables depend on it. The code in the simple file below the
+      project definition accordingly changes for instance to:
     
    
 
-    
-    
    Simple CMakeLists.txt
    
-
-    
    
-      The easiest way to write a CMakeLists.txt file is to use
-      an "autopilot" style macro. Here is an example for the step-1
-      tutorial program (plain text version) that can be used for simple
-      projects:
-      
    -
-    #
-    # Set the name of the project and target:
-    #
-    SET(TARGET "step-1")
-
-    #
-    # Declare all source files the target consists of:
-    #
-    SET(TARGET_SRC
-      ${TARGET}.cc
-      # You can specify additional files here!
-      )
-
-    #
-    # (Optional)
-    # Specify a list of files (file globs) that will be removed
-    # with the "make runclean" and "make distclean" targets.
-    # If empty, sensible default values will be used.
-    #
-    SET(CLEAN_UP_FILES
-      # a custom list of globs, e.g. *.log *.vtk
-      )
-
-    #
-    # (Optional)
-    # A custom command line that should be invoked by "make run". If empty,
-    # ./${TARGET} will be invoked.
-    #
-    SET(TARGET_RUN
-      # a custom command line, e.g. mpirun -np 2 ${TARGET}
-      )
-
-    #
-    # It is good practice to specify a version requirement:
-    #
-    CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
-
-    #
-    # Find and import the deal.II project configuration:
-    #
-    FIND_PACKAGE(deal.II 8.0 REQUIRED
-      HINTS
-        ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
-      #
-      # If the deal.II library cannot be found (because it is not installed at a
-      # default location or your project resides at an uncommon place), you
-      # can specify additional hints for search paths here, e.g.
-      # "$ENV{HOME}/workspace/deal.II"
-      #
-      )
-
-    #
-    # And let deal.II do the rest:
-    #
-    DEAL_II_INITIALIZE_CACHED_VARIABLES()
-    PROJECT(${TARGET})
-    DEAL_II_INVOKE_AUTOPILOT()
-      
    
+
    +ADD_LIBRARY(mylib libsrc1.cc libsrc2.cc libsrc3.cc)
 
-      This CMakeLists.txt is intended for use with a small
-      project and in-source build (i.e., one does not create a separate
-      build directory as we recommend for the deal.II build
-      in the readme
-      file). Using this input file, you can run cmake in the
-      source directory as follows:
-      
    +ADD_EXECUTABLE(mycode mycode.cc)
+DEAL_II_SETUP_TARGET(mycode)
+TARGET_LINK_LIBRARIES(mycode mylib)
+
    
 
-    $ cd step-1
-    $ cmake .
-
-    [...]
-
-    ###
-    #
-    #  Successfully set up project  step-1  with  deal.II-8.0.pre  found at
-    #      /usr/local
-    #
-    #  CMAKE_BUILD_TYPE:   Debug
-    #  TARGET_SRC:         step-1.cc
-    #
-    #  You can now run
-    #      $ make            - to compile and link the program
-    #      $ make run        - to (compile, link and) run the program
-    #
-    #      $ make debug      - to switch the build type to "Debug"
-    #      $ make release    - to switch the build type to "Release"
-    #      $ make edit_cache - to change (cached) configuration variables
-    #                          and rerun the configure and generate phases of CMake
-    #
-    #      $ make clean      - to remove the generated executable as well as
-    #                          all intermediate compilation files
-    #      $ make runclean   - to remove all output generated by the program
-    #      $ make distclean  - to clean the directory from _all_ generated
-    #                          files (includes clean, runclean and the removal
-    #                          of the generated build system)
-    #
-    #  Have a nice day!
-    #
-    ###
-      
    
-    
    
+
    When you have multiple targets,
+repeat the last line of code for each of them. Accordingly,
+a loop becomes even more
+attractive.
    
 
-    
-    
    Advanced CMakeLists.txt
    
+
    If you only have a single file or few files with common source
+code, an alternative to creating a library might be the option:
    
 
-    
    
-      For large projects the simple CMakeLists.txt presented
-      above tends to be too inflexible.
-      So, if you wish to have more control about targets the
-      DEAL_II_SETUP_TARGET macro might be of interest for you.
-      Here is a full example, again using the step-1 tutorial program as a
-      template:
-      
    +
    +ADD_EXECUTABLE(mycode mycode.cc common.cc)
+DEAL_II_SETUP_TARGET(mycode)
 
-    #
-    # It is good practise to specify a version requirement:
-    #
-    CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
-
-    #
-    # Find and import the deal.II project configuration:
-    #
-    FIND_PACKAGE(deal.II 8.0 REQUIRED
-      HINTS
-        ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
-      #
-      # If the deal.II library cannot be found (because it is not installed at a
-      # default location or your project resides at an uncommon place), you
-      # can specify additional hints for search paths here, e.g.
-      # "$ENV{HOME}/workspace/deal.II"
-      #
-      )
-
-    #
-    # Initialize cached variables: This will set the compiler and
-    # compiler flags to the values deal.II was configured with, as well as,
-    # CMAKE_BUILD_TYPE to "Debug".
-    # These values can be altered by editing the cache via
-    #   $ make edit_cache
-    #
-    DEAL_II_INITIALIZE_CACHED_VARIABLES()
-
-    #
-    # PROJECT has to be called after DEAL_II_INITIALIZE_CACHED_VARIABLES
-    # otherwise CMake will have already set the compiler and compiler
-    # flags.
-    #
-    PROJECT(step-1)
-
-    #
-    # Specify a target:
-    #
-    ADD_EXECUTABLE(step-1 step-1.cc)
-
-    #
-    # DEAL_II_SETUP_TARGET will set up all necessary include paths,
-    # preprocessor definitions and the link interface:
-    #
-    DEAL_II_SETUP_TARGET(step-1)
-
-    #
-    # (Optional)
-    # If you wish to have a "run" target for make, specify one:
-    #
-    ADD_CUSTOM_TARGET(run
-      COMMAND step-1
-      COMMENT "Run ${TARGET} with ${CMAKE_BUILD_TYPE} configuration"
-      )
-      
    
-   
    
+ADD_EXECUTABLE(mycode2 mycode2.cc common.cc)
+DEAL_II_SETUP_TARGET(mycode2)
+
    
 
+
    You should be aware though that common.cc will be
+compiled for each target, not only once.
    
 
-    
-    
     Advanced setup for deal.II in a CMake project 
    
+
     Matthias, is this correct?
    
 
-    
-    
     Finding the deal.II library 
    
 
+    
+    
    Adding a "run" target
    
+    
     
    
-    Finding the deal.II library should be no more than
-      
    +      If you wish to have a "run" target for make, like in the deal.II
+      tutorial, specify one this way (obviously, a single "run" target
+      can only run a single executable):
+    
    
 
+
    +ADD_CUSTOM_TARGET(run COMMAND mycode
+  COMMENT "Run with ${CMAKE_BUILD_TYPE} configuration")
+
    
+
+
+
+
     Advanced CMakeLists.txt
    
+
+
    +  
+  A Subsection for each of the following:
+  
+  - control statements (if and foreach)
+  
+  - file glob to pick up sources
+  
+  - DEAL_II_SETUP_TARGET and DEAL_II_INITIALIZE_CACHED_VARIABLES
+  revisited
+  
+  - add include dirs and compile definitions to a directory or target
+  
+  - provide a ./run folder and target
+  
+  - custom targets to switch between release and debug
+  
+  - installation
+
    
+
+
+
     Finding the deal.II library 
    
+
+
    
+  Finding the deal.II library should be no more than
+  
    +    
     FIND_PACKAGE(deal.II REQUIRED)
-      
    
-      in your CMakeLists.txt file.  You may have to hint for the location
-      of the deal.IIConfig.cmake file.  Either by directly
-      specifying deal.II_DIR to point to the path were the
-      deal.IIConfig.cmake file is located:
-      
    -
-    cmake -Ddeal.II_DIR=/path/to/the/config/file <...>
-      
    
-      or by specifying a search path via CMAKE_PREFIX_PATH,
-      e.g.
-      
    -
+  
    
+  in your CMakeLists.txt file.  You may have to hint for the location
+  of the deal.IIConfig.cmake file.  Either by directly
+  specifying deal.II_DIR to point to the path were the
+  deal.IIConfig.cmake file is located:
+  
+
     is deal.II_DIR still accurate? We use different
+above. BTW, that is an UGLY variable name.
    
+  
+
        
+cmake -Ddeal.II_DIR=/path/to/the/config/file <...>
+
+
    
+  or by specifying a search path via CMAKE_PREFIX_PATH,
+  e.g.
+  
    +    
     cmake -DCMAKE_PREFIX_PATH=~/workspace/local
-      
    
-      In this case deal.IIConfig.cmake will be searched
-      for in
-      
    -
+  
    
+  In this case deal.IIConfig.cmake will be searched
+  for in
+  
    +    
     ~/workspace/local/
     ~/workspace/local/lib/cmake/deal.II/
-      
    
-    
    
-
-
-    
-    
      deal.IIConfig.cmake 
    
+  
    
+
    
+
+
+
    Autopilot style CMakeLists.txt
    
+
+
    
+  If you want a make interface similar to the deal.II library and
+  its tutorial, namely maker targets for debug and release
+  versions, running the code and cleaning, the easiest way to
+  write a CMakeLists.txt file may be to use
+  an "autopilot" style macro. Here is a minimalistic example for the
+  step-1 tutorial program (plain text version) that can be used for simple
+  projects:
+
    
+
+
    +FIND_PACKAGE(deal.II 8.0 REQUIRED
+  HINTS
+    ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+    # You can specify additional hints for search paths here, e.g.
+    # $ENV{HOME}/workspace/deal.II
+)
+
+# Set the name of the project and target:
+SET(TARGET "step-1")
+
+# Declare all source files the target consists of:
+SET(TARGET_SRC
+  step-1.cc
+  # You can specify additional files here!
+)
+
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET} CXX)
+DEAL_II_INVOKE_AUTOPILOT()
+
    
 
-    
    
-      Importing the deal.IIConfig.cmake file via FIND_PACKAGE
-      will set a bunch of variables and macros; all of the form
-      DEAL_II_*. There is the usual duplet:
-      
    -
-    DEAL_II_INCLUDE_DIRS
-    DEAL_II_LIBRARIES
-      
    
-      (with debug and optimized keywords. For
-      compatibility there is also DEAL_II_LIBRARIES_DEBUG
-      and DEAL_II_LIBRARIES_RELEASE only specifying the
-      respective set of libraries.)
-    
    
 
-    
    
-      Interesting additional variables might be:
-      
    -
-    DEAL_II_USER_DEFINITIONS
-    DEAL_II_USER_DEFINITIONS_DEBUG
-    DEAL_II_USER_DEFINITIONS_RELEASE
-    DEAL_II_LINKER_FLAGS
-    DEAL_II_LINKER_FLAGS_DEBUG
-    DEAL_II_LINKER_FLAGS_RELEASE
-    DEAL_II_CXX_FLAGS
-    DEAL_II_CXX_FLAGS_DEBUG
-    DEAL_II_CXX_FLAGS_RELEASE
-
-    DEAL_II_TARGET_CONFIG
-    DEAL_II_TARGET
-    DEAL_II_TARGET_DEBUG
-    DEAL_II_TARGET_RELEASE
-      
    
-    
    
-
-
-
-    
-    
     Setting up necessary configuration variables 
    
-
-    
    
-    For actually using deal.II in your CMake
-    project some configuration steps are necessary.  These can be
-    either set via macros or by hand:
-      
      
-        
    1.  Configuration with the help of
-          two convencience macros:
-        
      
-          deal.IIConfig.cmake includes some convenience macros
-          to automatically setup all necessary configuration. A fully
-          functional, minimal example for the step-1 tutorial program is:
-	  
      -
-    FIND_PACKAGE(deal.II REQUIRED)
+      This CMakeLists.txt is intended for use with a small
+      project and in-source build (i.e., one does not create a separate
+      build directory as we recommend for the deal.II build
+      in the readme
+      file). Using this input file, you can run cmake in the
+      source directory as follows:
 
-    DEAL_II_INITIALIZE_CACHED_VARIABLES()
-
-    PROJECT(step-1)
-
-    ADD_EXECUTABLE(step-1 step-1.cc)
-    DEAL_II_SETUP_TARGET(step-1)
-	  
      
-        
      
-
-
-        
    2.  Configuration by hand:
-	  
        
-	    
      •  Set up all include directories for header files:
-	      
        -
-    INCLUDE_DIRECTORIES(${DEAL_II_INCLUDE_DIRS})
-	      
        
-
-	    
      • 
-	      deal.II usually ships with an optimized Release and a Debug version
-              of the library. So it is a good idea to set up
-              CMAKE_BUILD_TYPE accordingly:
-	      
        -
-    SET(CMAKE_BUILD_TYPE "Debug" CACHE
-      "Choose the type of build, options are: Debug, Release"
-      )
-	      
        
-
-	      A cached variable ensures that we can later switch the build type
-	      by editing the cache:
-	      
        -
-    make edit_cache
-	      
        
-
-	    
      •  Often, it is a good idea to use the same compiler and linker
-	      flags as the deal.II library.
-	      
        -
-    SET(CMAKE_CXX_FLAGS ${DEAL_II_CXX_FLAGS})
-    SET(CMAKE_CXX_FLAGS_RELEASE ${DEAL_II_CXX_FLAGS_RELEASE})
-    SET(CMAKE_CXX_FLAGS_DEBUG ${DEAL_II_CXX_FLAGS_DEBUG})
-	      
        
-	      (Optionally you can set up the variables with the
-	      CACHE to be able to edit them via ccmake or
-	      make edit_cache.)
-
-	    
      • 
-	      After this, specify your project name:
-	      
        -
-    PROJECT(myProject)
-	      
        
-	      This ensures that the compiler detection and platform setup
-	      that is issued by calling PROJECT(...) will run
-	      after we have set up our cached variables. This way it is
-	      our choice of variables that will be set and not the default
-	      values determined by the PROJECT(...) call.
-
-	    
      • 
-	      After defining your targets, e.g.
-	      
        -
-    ADD_EXECUTABLE(step-1 step-1.cc)
-	      
        
-	      you have to specify to link against deal.II and all external
-	      libraries:
-	      
        -
-    TARGET_LINK_LIBRARIES(step-1 ${DEAL_II_LIBRARIES})
-	      
        
-
-	    
      •  as well as using some preprocessor definitions:
-	      
        -
-    SET_TARGET_PROPERTIES(step-1 PROPERTIES
-      COMPILE_DEFINITIONS
-      "${DEAL_II_USER_DEFINITIONS}"
-      COMPILE_DEFINITIONS_DEBUG
-      "${DEAL_II_USER_DEFINITIONS_DEBUG}"
-      COMPILE_DEFINITIONS_RELEASE
-      "${DEAL_II_USER_DEFINITIONS_RELEASE}"
-      )
-	      
        
-
-	    
      •  Optionally, you can set the link flags to the one used by the
-	      deal.II library:
-	      
        -
-    SET_TARGET_PROPERTIES(step-1 PROPERTIES
-      LINK_FLAGS
-      "${DEAL_II_LINKER_FLAGS}"
-      LINK_FLAGS_DEBUG
-      "${DEAL_II_LINKER_FLAGS_DEBUG}"
-      LINK_FLAGS_RELEASE
-      "${DEAL_II_LINKER_FLAGS_RELEASE}"
-      )
-	      
        
-	      And this it is.
-
-	  
      
-
-          
    3. 
-            Exported targets:
-
-            
      
-              For an advanced setup in a big CMake project
-              deal.IIConfig.cmake
-              provides information about the deal.II
-              installation with traditional variables, see
-              here, as well as external CMake
-              targets with link interface for direct inclusion:
-              
      -
-    INCLUDE(${DEAL_II_TARGET_CONFIG})
-
-    ADD_EXECUTABLE(step-1
-      step-1.cc
-      ${DEAL_II_TARGET}
-      )
-              
      
-            
      
-      
    
-    
    
+
    +$ cd step-1
+$ cmake .
+
+[...]
+
+###
+#
+#  Successfully set up project  step-1  with  deal.II-99.99.svn30300  found at
+#      /usr
+#
+#  CMAKE_BUILD_TYPE:         Debug
+#
+#  You can now run
+#      $ make                - to compile and link the program
+#      $ make run            - to (compile, link and) run the program
+#
+#      $ make debug          - to switch the build type to Debug
+#      $ make release        - to switch the build type to Release
+#
+#      $ make edit_cache     - to change (cached) configuration variables
+#                              and rerun the configure and generate phases of CMake
+#
+#      $ make strip_comments - strip the source files in this
+#                              directory off the documentation comments
+#      $ make clean          - to remove the generated executable as well as
+#                              all intermediate compilation files
+#      $ make runclean       - to remove all output generated by the program
+#      $ make distclean      - to clean the directory from _all_ generated
+#                              files (includes clean, runclean and the removal
+#                              of the generated build system)
+#
+#      $ make help           - to view this message again
+#
+#  Have a nice day!
+#
+###
+
    
+
+      There are two additional configuration options (in addition to
+      TARGET and TARGET_SRC) that can be set via
+      variables before DEAL_II_INVOKE_AUTOPILOT() is called
+      (plain text version):
+
+
    +# (Optional)
+# Specify a list of files (file globs) that will be removed
+# with the "make runclean" and "make distclean" targets.
+# (If empty, sensible default values will be used.)
+SET(CLEAN_UP_FILES
+  # a custom list of globs, e.g. *.log *.vtk
+)
+
+# (Optional)
+# A custom command line that should be invoked by "make run".
+# (If empty, ./${TARGET} will be invoked.)
+SET(TARGET_RUN
+  # a custom command line, e.g. mpirun -np 2 ${TARGET}
+)
+
    
+
+
+
+
      deal.IIConfig.cmake 
    
+
+
    
+  Importing the deal.IIConfig.cmake file via FIND_PACKAGE
+  will set the following variables and macros; all of the form
+  DEAL_II_*:
+
    
+
      
+A long list with detailed explanation...  
+
    
+  
 
 
 
     
-    
     Legacy Make.global_options 
    
+    
     Legacy Make.global_options 
    
 
     
    
       Before version 8.0, deal.II used the
diff --git a/deal.II/doc/development/toc.html b/deal.II/doc/development/toc.html
index a213bd675b..72d76a8232 100644
--- a/deal.II/doc/development/toc.html
+++ b/deal.II/doc/development/toc.html
@@ -41,7 +41,7 @@
          This page provides extensive information about configuration and
          installation with the CMake build system.
 
-         
  •  
    Using
+         
  •  
    Using CMake and
            deal.II in a user project:
          This page covers how to use CMake in your own project.
          An overview of how to retrieve necessary information
diff --git a/deal.II/doc/development/writing-documentation.html b/deal.II/doc/development/writing-documentation.html
index 4ab4d3e318..e9b037dffd 100644
--- a/deal.II/doc/development/writing-documentation.html
+++ b/deal.II/doc/development/writing-documentation.html
@@ -179,9 +179,9 @@
 
     
  •  Typewriter font:
          
    In order to use typewriter font for instance for function
-         arguments or variables, use the <tt> HTML
+         arguments or variables, use the <code> HTML
          tag. For a single word, you can also use the form @p
-         one_word_without_spaces.
    
+         one_word_without_spaces. The <tt> is obsolete in HTML5
    
 
          
    
          If you refer to member variables and member functions
diff --git a/deal.II/doc/documentation.html b/deal.II/doc/documentation.html
index 941beef63e..10620b4c56 100644
--- a/deal.II/doc/documentation.html
+++ b/deal.II/doc/documentation.html
@@ -2,208 +2,203 @@
           "http://www.w3.org/TR/html4/loose.dtd">
 
   
-    
     The deal.II Online Documentation
+    
+    
     
     
     
     
     
   
-  
-
-
-
-
-    
    Online documentation
    
-
-    
-      
-      
-      
-    
    Table of contents
    
-	  
-	
    
-
-
-    
    
-      deal.II comes with quite extensive online
-      documentation, which falls into several categories. They are
-      listed below.
-    
    
-
-
-    
    ReadMe
    
-
-    
    
-      First, there is the README file of
-      the library. It gives information on system requirements,
-      installation, and copyright.
-    
    
-
-    
    Tutorials
    
-
-    
    
-      The tutorial is
-      for new users of the library. It explains the basic elements of
-      finite element programs based on the library, and provides small
-      example programs. You can also reach the tutorial link from the menu bar
-      at the left.
-    
    
-
-    
    
-      Note that along with the rest of the documentation, the local HTML pages
-      of the tutorials need to be generated first. Please follow the
-      instructions in the
-      ReadMe file on how to do
-      this.
-    
    
-
-
-    
    Programming interface
    
-
-    
    
-      The programming
-	interface/manual pages are what you will need most often as a
-      reference of all classes, functions and variables in the library. They
-      are extensively documented (presently more than 5000 pages if
-      printed), and generated by Doxygen. There is also a quick link to these pages from
-      the menu bar at the left.
-    
    
-
-    
    
-      The documentation uses many of the features of Doxygen. In particular,
-      the link takes you straight to a page that lists modules, a way
-      to group classes with similar purposes. This may be what you are probably
-      interested in when you are not yet familiar with the library. If you
-      already know your way around and want to look up the signature of a
-      particular member function, for example, you may go to the class view or
-      choose any of the other ways in which Doxygen allows you to navigate
-      through the documentation.
-    
    
-
-    
    
-      Just as for the tutorials, the API docs need to be generated first,
-      if you download deal.II. Please follow the
-      instructions in the
-      ReadMe file on how to do
-      this.
-    
    
-
-
-    
    Development information
    
-
-    
    
-      If you plan to develop your own applications using
-      deal.II, it might be worth to take a look
-      at this
-	page.
-    
    
-
-
-    
    Reports
    
-
-    
    
-      There are a number of papers that describe the overall design or
-      individual aspects of deal.II. Please refer to the
-      papers listed at the top of the publications page for full references. In addition,
-      there are the following reports:
-      
      
-        
    •  
      
-            A publication on the
-	    basics of dimension independent
-              programming in deal.II (by Wolfgang
-            Bangerth). This explains some of the fundamentals of the
-            library's design, though rather briefly due to the page
-            limit imposed for the publication.
-	  
      
-
-	  
      
-	    This publication appeared in
-            the proceedings of the 16th IMACS World Congress,
-            Lausanne, Switzerland, 2000. See the
-	    publications
-	    page for the full reference as well as a way to obtain a
-            reprint of the original paper as a postscript or PDF file.
-	  
      
-
-	
    •  
      
-            A brief report about the classes involved in the
-	    iterators and accessors
-	    used to access cells in the triangulations and the data
-	    thereon, such as their geometry, the number of the
-	    degrees of freedom, data values, etc. You may want to
-            read this object or parts of it if you want to know what
-            information is available from iterators (by Wolfgang
-            Bangerth).
-	  
      
-
-        
    •  
      
-	    A brief report on mapping functions of higher polynomial
-	    degrees (by Ralf Hartmann). These allow the realisation
-	    of higher order boundary approximations. This report
-	    gives the theoretical background of the
-	    MappingQ class.
-	  
      
-
-        
    •  
      
-            A report on
-	    assembling matrices (by Wolfgang
-            Bangerth). This report looks in particular at how
-            matrices are assembled for vector-valued problems, and
-            how things change if vector-valued finite elements are
-            used for which the shape functions have more than one
-            non-zero vector component, such as Nedelec or
-            Raviart-Thomas elements. There are printable
-            versions of the report
-	    in postscript and
-	    PDF format.
-	  
      
-
-        
    •  
      
-            A very detailed report on
-	    Nedelec elements (by Anna
-            Schneebeli, University of Basel, Switzerland). It
-            explains the construction and application of Nedelec
-            edge elements for H-curl spaces, as used, for example,
-            in the numerical solution of the Maxwell equations. It
-            also gives numerical results obtained with deal.II.
-	    There is also a postscript version of the report.
-	  
      
-
-        
    •  
      
-            A report on the
-	    codimension one capabilities
-	    of the library (by Antonio DeSimone, Luca Heltai
-	    and Cataldo Manigrasso, SISSA, Trieste, Italy). It
-            explains in detail how to use the
-	    library for the solution of problems defined on codimension
-	    one manifolds, such as, for example, Boundary Element Methods.
-      
    
-
-
-    
    Publications
    
+
+
+
    Online documentation
    
+
+
    
+  
      
+    
    1.  ReadMe
      2. 
+    
    2.  Tutorials
      3. 
+    
    3.  Programming interface
      4. 
+    
    4.  Development information
      5. 
+    
    5.  Reports
      6. 
+    
    6.  Publications
      7. 
+    
    7.  Copyright and licence
      8. 
+  
    
+
    
+
+
    
+  deal.II comes with quite extensive online
+  documentation, which falls into several categories. They are
+  listed below.
+
    
+
+
    ReadMe
    
+
+
    
+  First, there is the README file of
+  the library. It gives information on system requirements,
+  installation, and copyright.
+
    
+
+
    Tutorials
    
+
+
    
+  The tutorial is
+  for new users of the library. It explains the basic elements of
+  finite element programs based on the library, and provides small
+  example programs. You can also reach the tutorial link from the menu bar
+  at the left.
+
    
+
+
    
+  Note that along with the rest of the documentation, the local HTML pages
+  of the tutorials need to be generated first. Please follow the
+  instructions in the
+  ReadMe file on how to do
+  this.
+
    
+
+
+
    Programming interface
    
+
+
    
+  The programming
+    interface/manual pages are what you will need most often as a
+  reference of all classes, functions and variables in the library. They
+  are extensively documented (presently more than 5000 pages if
+  printed), and generated by Doxygen. There is also a quick link to these pages from
+  the menu bar at the left.
+
    
+
+
    
+  The documentation uses many of the features of Doxygen. In particular,
+  the link takes you straight to a page that lists modules, a way
+  to group classes with similar purposes. This may be what you are probably
+  interested in when you are not yet familiar with the library. If you
+  already know your way around and want to look up the signature of a
+  particular member function, for example, you may go to the class view or
+  choose any of the other ways in which Doxygen allows you to navigate
+  through the documentation.
+
    
+
+
    
+  Just as for the tutorials, the API docs need to be generated first,
+  if you download deal.II. Please follow the
+  instructions in the
+  ReadMe file on how to do
+  this.
+
    
+
+
+
    Development information
    
+
+
    
+  If you plan to develop your own applications using
+  deal.II, it might be worth to take a look
+  at this
+    page.
+
    
+
+
+
    Reports
    
+
+
    
+  There are a number of papers that describe the overall design or
+  individual aspects of deal.II. Please refer to the
+  papers listed at the top of the publications page for full references. In addition,
+  there are the following reports:
+  
      
+    
    • 
+      
      
+	A publication on the
+	basics of dimension independent
+          programming in deal.II (by Wolfgang
+	Bangerth). This explains some of the fundamentals of the
+	library's design, though rather briefly due to the page
+	limit imposed for the publication.
+      
      
+      
+      
      
+	This publication appeared in
+	the proceedings of the 16th IMACS World Congress,
+	Lausanne, Switzerland, 2000. See the
+	publications
+	page for the full reference as well as a way to obtain a
+	reprint of the original paper as a postscript or PDF file.
+      
      
+    
      • 
+    
+    
    • 
+      A brief report about the classes involved in the
+      iterators and accessors
+      used to access cells in the triangulations and the data
+      thereon, such as their geometry, the number of the
+      degrees of freedom, data values, etc. You may want to
+      read this object or parts of it if you want to know what
+      information is available from iterators (by Wolfgang
+      Bangerth).
+    
      • 
+    
+    
    • 
+      A brief report on mapping functions of higher polynomial
+      degrees (by Ralf Hartmann). These allow the realisation
+      of higher order boundary approximations. This report
+      gives the theoretical background of the
+      MappingQ class.
+    
      • 
+
+    
    • 
+      A report on
+      assembling matrices (by Wolfgang
+      Bangerth). This report looks in particular at how
+      matrices are assembled for vector-valued problems, and
+      how things change if vector-valued finite elements are
+      used for which the shape functions have more than one
+      non-zero vector component, such as Nedelec or
+      Raviart-Thomas elements. There are printable
+      versions of the report
+      in postscript and
+      PDF format.
+    
      • 
+    
+    
    • 
+      A very detailed report on
+      Nedelec elements (by Anna
+      Schneebeli, University of Basel, Switzerland). It
+      explains the construction and application of Nedelec
+      edge elements for H-curl spaces, as used, for example,
+      in the numerical solution of the Maxwell equations. It
+      also gives numerical results obtained with deal.II.
+      There is also a postscript version of the report.
+    
      • 
+    
+    
    • 
+      A report on the
+      codimension one capabilities of the library
+         (by Antonio DeSimone, Luca Heltai and Cataldo Manigrasso,
+         SISSA, Trieste, Italy). It explains in detail how to use the
+         library for the solution of problems defined on codimension
+         one manifolds, such as, for example, Boundary Element
+         Methods.
+      
      • 
+  
    
+
+
+    
    Publications
    
 
     
    
       A list of
@@ -211,26 +206,28 @@
       deal.II can be found on
       this page.
     
    
-
-
-    
    Copyright and License
    
-
-    
     deal.II is copyright of the
-      deal.II authors.
-      It is a joint effort of the current maintainers, Wolfgang
-      Bangerth, Timo Heister, Guido Kanschat, and quite a number of
-      contributors.  deal.II is distributed under an
+    
+    
+    
    Copyright and License
    
+    
+    
    
+      deal.II is copyright of the
+      deal.II authors.  It is
+      a joint effort of the current maintainers, Wolfgang Bangerth,
+      Timo Heister, Guido Kanschat, and quite a number of
+      contributors.  deal.II is distributed under
+      an
       OpenSource
       license. The full text of the license can be found in
       license.html.
     
    
 
-    
    Contact
    
+    
    Contact
    
 
     
    We have a discussion group for discussion
       of issues of general interest to users and developers of
       deal.II.
    
-
+    
     
    More specific questions may be sent to the authors
       immediately at the address firstName.lastname at
        dealii.org.
    
@@ -242,5 +239,5 @@
 								 alt="Valid HTML 4.01!" height="31" width="88">
       
    
     
-  
+
 
diff --git a/deal.II/doc/readme.html b/deal.II/doc/readme.html
index b08a508c30..6bec165b65 100644
--- a/deal.II/doc/readme.html
+++ b/deal.II/doc/readme.html
@@ -1,14 +1,16 @@
-
-
+
+
+
   
     The deal.II Readme
-    
-    
-    
-    
-    
-    
+    
+    
+    
+    
+    
+    
+    
   
 
   
@@ -26,113 +28,112 @@
     
 
 
-    
    Installation instructions and further information
-      on deal.II 
    
-
-    
-      
-      
-      
-    
    Table of contents
    
-	  
-	
    
-
-
-    
-    
    System requirements
    
-
-    
-    
    Supported platforms
    
-
-    
    
-      deal.II is mostly developed on Linux using the
-      GCC compiler. However, it is not
-      platform specific and we strive to keep the source code C++
-      Standard compliant.
-    
    
-
-    
    
-      deal.II supports at least the following
-      platforms:
-      
    
+
    
+  Most other combinations of POSIX-style operating systems and C++
+  Standard compliant compilers should also work. If they don't,
+  please report it as a bug.
+
    
 
-    
    
-      Please note that the Microsoft Visual C++ compiler, that is widely
-      used on Microsoft Windows, is not supported. If you need a native
-      Windows target, use the MinGW compiler instead.
-    
    
+
    
+  Please note that the Microsoft Visual C++ compiler, that is widely
+  used on Microsoft Windows, is not supported. If you need a native
+  Windows target, use the MinGW compiler instead.
+
    
 
 
-    
-    
    Additional software requirements
    
+
+
    Additional software requirements
    
 
     
    
       In order to compile and use deal.II
       you need to have the following programs installed:
     
    
     
 
-
-
-
-    
+    
     
    Installation
    
-
-        
+    
+    
     
    Unpacking
    
 
     
    
       The whole library usually comes as a tar.gz file,
       which is a file archive compressed with gzip. After downloading it,
       unpack it using either
-      
    -
-	gunzip deal.II-X.Y.Z.tar.gz
-	tar xf deal.II-X.Y.Z.tar
-      
    
-      or, if you have GNU tar with
-      
    -
-	tar -xvf deal.II-X.Y.Z.tar.gz
-      
    
-      Note: You will want to hang on to the source files
-      of deal.II after installation as it makes developing
-      much simpler. Consequently, you should do the steps above in a permanent
-      directory, not on /tmp as one often does when installing
-      software.
     
    
+
    +  gunzip deal.II-X.Y.Z.tar.gz
+  tar xf deal.II-X.Y.Z.tar
+
    
+
    or, if you have GNU tar with
    
+
    +  tar -xvf deal.II-X.Y.Z.tar.gz
+
    
+
+
    Note: You will want to hang on to the source files
+  of deal.II after installation as it makes
+  developing much simpler. Consequently, you should do the steps above
+  in a permanent directory, not on /tmp as one often does
+  when installing software.
+
    
 
         
     
    Configuring and building the library
    
@@ -219,21 +219,23 @@
       CMake
       integrated configuration and build system. Unpacking will create a
       subdirectory deal.II/ in the current directory. Then do the
-      following steps:
-      
    -
-	mkdir build
-	cd build
-	cmake -DCMAKE_INSTALL_PREFIX=/path/to/install/dir ../deal.II
-	make install
-      
    
-      These steps compile, link and install the deal.II library. The
-      whole process should take between a few minutes and an hour, depending
-      on your machine.
-    
    
+      following steps:
    
+
+
    +  mkdir build
+  cd build
+  cmake -DCMAKE_INSTALL_PREFIX=/path/to/install/dir ../deal.II
+  make install
+
    
+
+
    
+  These steps compile, link and install the deal.II library. The
+  whole process should take between a few minutes and an hour, depending
+  on your machine.
+
    
 
     
    
-      Note:
+      Note:
    
       
      
 	
    •  /path/to/install/dir is the directory which deal.II
 	  should be installed into. This can be a directory in your home
@@ -242,6 +244,7 @@
 	  Another option is to use something like `pwd`/../installed/ (note the
 	  backticks). Make sure the installation directory is not the same
 	  as the location where you unpacked deal.II/.
+	
      • 
 
 	
    •  If your machine has multiple processors, use make
 	  -jN in the last step, where N is the
@@ -249,10 +252,12 @@
 	  to use at any given time. Allowing make to use
 	  more simultaneous build processes (assuming you have that many
 	  processor cores) will greatly lower the build time.
+	
      • 
 
 	
    •  If you do not intend to modify the deal.II
 	  sources and recompile things, then you can remove
 	  the build/ directory after the last step.
+	
      • 
 
 	
    •  In principle, after installing deal.II, you
 	  can remove the source directory as well (i.e., the directory into
@@ -262,10 +267,12 @@
 	  files around anyway, for one reason: When
 	  debugging you often end up with assertions for which you'd like to
           see the place in the library's source files that triggered it.
+	
      • 
 
 	
    •  The deal.II CMake system can accept a
 	  significant number of configuration parameters. See the
 	  discussion below.
+	
      • 
 
 	
    •  If you are changing part of the deal.II
 	code itself, you can re-compile the library using only the
@@ -276,41 +283,41 @@
 	However, this sometimes leads to surprising results and you may
 	not get exactly what you were hoping for. For more information,
 	see here.
-
+	
      • 
       
    
-    
    
 
     
    
       The commands above build and install the deal.II
       libraries in two variants:
-      
      
-	
    • 
-	  
      Debug mode: This version of the
-	    library is compiled with compiler flags so
-	    that the library contains information that can be used by debuggers.
-	  
      
-
-	  
      
-	    In addition, this library contains a great number of safety
-	    checks on most arguments of all functions you could possibly call.
-	    These assertions have proven to be an
-	    invaluable means to finding programming bugs since they will
-	    almost always abort your program if something goes wrong. In our
-	    experience, more than ninety per cent of all errors are invalid
-	    parameters (such as vectors having the wrong size, etc.) and they
-	    are usually found almost instantaneously, displaying the file
-	    name and line number of where the problem occurred.
-	  
      
-
+    
      
+    
        
+      
      • 
+	
        Debug mode: This version of the
+	  library is compiled with compiler flags so
+	  that the library contains information that can be used by debuggers.
+	
        
+	
+	
        
+	  In addition, this library contains a great number of safety
+	  checks on most arguments of all functions you could possibly call.
+	  These assertions have proven to be an
+	  invaluable means to finding programming bugs since they will
+	  almost always abort your program if something goes wrong. In our
+	  experience, more than ninety per cent of all errors are invalid
+	  parameters (such as vectors having the wrong size, etc.) and they
+	  are usually found almost instantaneously, displaying the file
+	  name and line number of where the problem occurred.
+	
        
+	
        • 
+	
 	
      •  Optimized mode: You will want to link with this version of
 	    the library once you know that your program is working as
 	    expected. It does not contain the safety checks any more
 	    and is compiled with aggressive compiler optimizations. The
 	    resulting executables are smaller and will run between 2 and 10
 	    times faster than the debug executables.
-	  
        
+	
        • 
       
      
-    
      
 
     
      
       At this point, you have generated everything necessary to write
@@ -327,17 +334,20 @@
       be found at the 
       http://www.dealii.org/ domain  can also be generated locally. To do
       so, invoke cmake in the build instructions above as follows:
-      
      -
-	cmake -DDEAL_II_COMPONENT_DOCUMENTATION=ON -DCMAKE_INSTALL_PREFIX=/path/install/dir ../deal.II
-      
      
-      For this to succeed, you will need Perl 5.x,
-      doxygen
-      and dot (which is part of
-      the GraphViz
-      package) to be installed.
     
      
+    
+
      +  cmake -DDEAL_II_COMPONENT_DOCUMENTATION=ON -DCMAKE_INSTALL_PREFIX=/path/install/dir ../deal.II
+
      
+
+
      
+  For this to succeed, you will need Perl 5.x,
+  doxygen
+  and dot (which is part of
+  the GraphViz
+  package) to be installed.
+
      
 
     
      
       Upon calling make and make install, this will
@@ -374,8 +384,8 @@
       options.
     
      
 
-    
-      
      Selecting optional behavior
      
+    
+      
      Selecting optional behavior
      
 
     
        
       
      • 
@@ -388,6 +398,7 @@
 	  threading off by passing the -DDEAL_II_WITH_THREADS=OFF
 	  argument to cmake.
 	
        
+      
        • 
 
       
      • 
 	
        
@@ -399,6 +410,7 @@
 	  specifically do not want to use it, then
 	  pass -DDEAL_II_WITH_MPI=OFF.
 	
        
+      
        • 
 
       
      • 
         
        
@@ -416,6 +428,7 @@
           -DDEAL_II_COMPONENT_PARAMETER_GUI=ON. Beware of the
           fact that parameter_gui needs development packages for Qt.
 	
        
+      
        • 
 
       
      • 
   
        
@@ -433,12 +446,12 @@
     with Trilinos, Trilinos 11.2 or later is needed (there is a bug in
     Trilinos 11.0 that prevents the use of 64bit indices).
   
        
-
+      
        • 
     
      
 
 
 
-          
+    
     
      Optional interfaces to other software packages
      
 
     
      
@@ -460,8 +473,7 @@
       can be configured.
     
      
 
-    
      
-    Notes:
+    
      Notes:
      
       
        
         
      • 
           The majority of libraries mentioned below should be readily
@@ -470,6 +482,7 @@
           e.g. ending in -dev or -devel.
           After that cmake will automatically find the
           library and use it.
+      
        • 
         
      • 
           Configuring the interface to a self compiled package,
           say foo can usually be done by specifying
@@ -480,17 +493,18 @@
           you invoke cmake in a fresh build directory.
           Any value passed on the command line wins over a value that
           may be found in an environment variable.
+      
        • 
         
      • 
           To explicitly enable or disable support for
           a library foo use the argument
           -DDEAL_II_WITH_FOO=ON resp.
           -DDEAL_II_WITH_FOO=OFF for cmake.
+      
        • 
     
      
-    
      
 
     
      
-      
-      
      ARPACK
      
+      
      
+	ARPACK
      
       
      
 	
      
 	  
       
      
 
-      
       
      
+	
 	BLAS,
 	LAPACK
       
      
@@ -531,8 +545,10 @@
 	
      
       
 
-      
-      
      FunctionParser
      
+      
      
+	
+	FunctionParser
+      
      
       
      
 	
      
 	  FunctionParser
@@ -549,8 +565,8 @@
 	
      
       
      
 
-      
-      
      HDF5
      
+      
+      
      HDF5
      
       
      
 	
      
 	  The HDF5 library
@@ -564,8 +580,8 @@
 	
      
       
      
 
-      
-      
      METIS
      
       
      
 	
      
@@ -587,8 +603,8 @@
 	
      
       
      
 
-      
-      
      MUMPS
      
+      
+      
      MUMPS
      
       
      
 	
      
 	  MUMPS
@@ -606,24 +622,23 @@
 	
      
       
      
 
-      
-      
      NetCDF
      
+      
+      
      NetCDF
      
       
      
-	
      
-	  NetCDF is a library that provides services for
-	  reading and writing mesh data (and many other
-	  things). deal.II can use it to read meshes via
-          one of the functions of the GridIn class.
-	  NetCDF should be readily packaged by most
-          Linux distributions. (Don't forget to install a development
-	  version of the library). To use a self compiled version, pass
-	  -DNETCDF_DIR=/path/to/netcdf to cmake.
+	NetCDF is a library that provides services for
+	reading and writing mesh data (and many other
+	things). deal.II can use it to read meshes via
+        one of the functions of the GridIn class.
+	NetCDF should be readily packaged by most
+        Linux distributions. (Don't forget to install a development
+	version of the library). To use a self compiled version, pass
+	-DNETCDF_DIR=/path/to/netcdf to cmake.
       
      
 
-      
-      
      p4est
      
+      
+      
      p4est
      
       
      
 	
      
 	  p4est
@@ -638,8 +653,8 @@
 	
      
       
      
 
-      
-      
      PETSc
      
       
      
 	
      
@@ -664,8 +679,8 @@
 	
      
       
      
 
-      
-      
      SLEPc
      
+      
+      
      SLEPc
      
       
      
 	
      
 	  
       
      
 
-      
-      
      Threading Building Blocks (TBB)
      
       
      
 	
      
@@ -702,8 +717,8 @@
 	
      
       
      
 
-      
-      
      Trilinos
      
+      
+      
      Trilinos
      
       
      
 	
      
 	  Trilinos is a
@@ -727,8 +742,8 @@
 	
      
       
      
 
-      
-      
      UMFPACK
      
       
      
 	
      
@@ -759,7 +774,7 @@
     
      
 
 
-    
+    
     
      More information on configuring and building the library
      
 
     
      
@@ -821,27 +836,31 @@
       information.
     
      
     
      
-      Please note:
+      Please note:
      
       
        
 	
      • 
           Detailed license information can be found following
           this link.
+	
        • 
 	
      • 
           As a contributor to this project, you agree that any of your
           contributions be licensed under the same terms and conditions as
           the license of the deal.II project granted to you.
+	
        • 
         
      • 
           We, the deal.II authors, do not require
           copyright assignments for contributions. This means that the
           copyright for code contributions in the deal.II project is held by
           its respective contributors who have each agreed to release their
           contributed code under the terms of the LGPL v2.1 or later.
+	
        • 
 	
      • 
           In addition to the terms imposed by the LGPL (version 2.1 or
           later), we ask for the courtesy that scientific publications
           presenting results obtained with this libraries acknowledge its
           use. Please cite one of the papers referenced
           here.
+	
        • 
 	
      • 
           deal.II can interface with a number of other packages that you either
@@ -856,13 +875,21 @@
             target="_top">FunctionParser, courtesy of their authors.
           These are also covered by their own licenses; please refer to
           their webpages for more information.
+	
        • 
       
      
 
-      
      
+      
      
 
       
      
 	The deal.II Group
 	$Date$
       
      
+    
      
+      
+	Valid HTML 4.01!
+      
+    
      
   
 
diff --git a/deal.II/doc/screen.css b/deal.II/doc/screen.css
index 65610cf281..f41edfa72e 100644
--- a/deal.II/doc/screen.css
+++ b/deal.II/doc/screen.css
@@ -6,6 +6,7 @@ body {
 	margin: 10px;
 	padding: 0;
 	font-family: sans-serif;
+	counter-reset: section;
 }
 
 body.title {
@@ -28,6 +29,23 @@ body.navbar {
 	font-size: 12px;
 }
 
+div.toc {
+    width: 50%;
+    color: black;
+    border: 1px solid #aaa;
+    background-color: #f9f9f9;
+    padding: 5px;
+    font-size: 88%;
+}
+
+div.toc:before {
+    display: block;
+    content: "Table of contents";
+    text-align: center;
+    font-size: 88%;
+    font-weight: bold;
+}
+
 frameset {
 	background-color: white;
 	border-color: black;
@@ -84,17 +102,20 @@ h1, h2, h3, h4, h5, h6 {
 }
 
 h1 {
-	font-size: 175%;
+	font-size: 175%;    
+	counter-reset: section;
 }
 
 h2 {
 	font-size: 150%;
+	counter-reset: subsection;
 }
 
 h3 { font-size: 140%;
      padding-top: 1.5em;
      padding-bottom: 0.17em;
      border-bottom: 1px dashed #aaaaaa; }
+
 h4 { font-size: 110%; 
      border-bottom: none;
      font-weight: bold;
@@ -108,6 +129,16 @@ h6 { font-size: 80%;
      font-weight: bold;
 }
 
+h2:before {
+    counter-increment: section;
+    content: counter(section) ". ";
+}
+
+h3:before {
+    counter-increment: subsection;
+    content: counter(section) "." counter(subsection) ". ";
+}
+
 div.right {
   text-align: right;
 }
@@ -121,13 +152,18 @@ ul {
     margin: 0.3em 0 0 1.5em;
     padding:0;
 }
+
 ol {
     line-height: 1.5em;
     margin: 0.3em 0 0 3.2em;
     padding:0;
     list-style-image: none;
 }
-li { margin-bottom: 0.1em; }
+
+li {
+    margin-bottom: 0.1em;
+}
+
 dt {
     font-weight: bold;
     margin-bottom: 0.1em;
@@ -180,6 +216,26 @@ pre {
 	line-height: 1.1em;
 }
 	
+pre.cmake {
+	padding: 1em;
+	text-align: left;
+	text-indent: 0;
+	border: 1px solid #d06fab;
+	color: Black;
+	background-color: #f9f9f9;
+	line-height: 1.1em;
+}
+	
+pre.sample {
+	padding: 1em;
+	text-align: left;
+	text-indent: 0;
+	border: 1px dashed #2f6fab;
+	color: Black;
+	background-color: #f9f9f9;
+	line-height: 1.1em;
+}
+	
 table.navbar { }
 table.tutorial {
     color: black;
@@ -239,6 +295,14 @@ td.build {
 	text-decoration: line-through;
 }
 
+.todo:before {
+    content: "TODO: ";
+}
+
+.todo {
+    color: Red;
+}
+
 .deprecated {
 	text-decoration: line-through;
 }
diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc
index 37c7826788..30a8f36706 100644
--- a/deal.II/examples/step-42/step-42.cc
+++ b/deal.II/examples/step-42/step-42.cc
@@ -21,7 +21,6 @@
  *          Timo Heister, Texas A&M University, 2013
  */
 
-
 // @sect3{Include files}
 // We are using the the same
 // include files as in step-41:
@@ -80,8 +79,9 @@
 
 #include 
 
-namespace Step42 {
-using namespace dealii;
+namespace Step42
+{
+  using namespace dealii;
 
 // @sect3{The Input class template}
 
@@ -110,135 +110,158 @@ using namespace dealii;
 // update_solution_and_constraints () of
 // the class PlasticityContactProblem.
 
-template
-class Input {
-public:
-	Input(const char* _name) :
-			name(_name), mpi_communicator(MPI_COMM_WORLD), pcout(std::cout,
-					(Utilities::MPI::this_mpi_process(mpi_communicator) == 0)), obstacle_data(
-					0), hx(0), hy(0), nx(0), ny(0) {
-		read_obstacle(name);
-	}
-
-	double
-	hv(int i, int j);
-
-	double
-	obstacle_function(double x, double y);
-
-	void
-	read_obstacle(const char* name);
-
-private:
-	const char* name;
-	MPI_Comm mpi_communicator;
-	ConditionalOStream pcout;
-	std::vector obstacle_data;
-	double hx, hy;
-	int nx, ny;
-};
+  template 
+    class Input
+    {
+      public:
+        Input (
+            const char* _name)
+            :
+                name(_name),
+                mpi_communicator(MPI_COMM_WORLD),
+                pcout(std::cout,
+                    (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)),
+                obstacle_data(0),
+                hx(0),
+                hy(0),
+                nx(0),
+                ny(0)
+        {
+          read_obstacle(name);
+        }
+
+        double
+        hv (
+            int i, int j);
+
+        double
+        obstacle_function (
+            double x, double y);
+
+        void
+        read_obstacle (
+            const char* name);
+
+      private:
+        const char* name;
+        MPI_Comm mpi_communicator;
+        ConditionalOStream pcout;
+        std::vector obstacle_data;
+        double hx, hy;
+        int nx, ny;
+    };
 
 // This function is used in obstacle_function ()
 // to provide the proper value of the obstacle.
-template
-double Input::hv(int i, int j) {
-	assert(i >= 0 && i < nx);
-	assert(j >= 0 && j < ny);
-	return obstacle_data[nx * (ny - 1 - j) + i]; // i indiziert x-werte, j indiziert y-werte
-}
+  template 
+    double
+    Input::hv (
+        int i, int j)
+    {
+      assert(i >= 0 && i < nx);
+      assert(j >= 0 && j < ny);
+      return obstacle_data[nx * (ny - 1 - j) + i]; // i indiziert x-werte, j indiziert y-werte
+    }
 
 // obstacle_function () calculates the bilinear interpolated
 // value in the point (x,y).
-template
-double Input::obstacle_function(double x, double y) {
-	int ix = (int) (x / hx);
-	int iy = (int) (y / hy);
-
-	if (ix < 0)
-		ix = 0;
-
-	if (iy < 0)
-		iy = 0;
-
-	if (ix >= nx - 1)
-		ix = nx - 2;
-
-	if (iy >= ny - 1)
-		iy = ny - 2;
-
-	double val = 0.0;
-	{
-		FullMatrix H(4, 4);
-		Vector X(4);
-		Vector b(4);
-
-		double xx = 0.0;
-		double yy = 0.0;
-
-		xx = ix * hx;
-		yy = iy * hy;
-		H(0, 0) = xx;
-		H(0, 1) = yy;
-		H(0, 2) = xx * yy;
-		H(0, 3) = 1.0;
-		b(0) = hv(ix, iy);
-
-		xx = (ix + 1) * hx;
-		yy = iy * hy;
-		H(1, 0) = xx;
-		H(1, 1) = yy;
-		H(1, 2) = xx * yy;
-		H(1, 3) = 1.0;
-		b(1) = hv(ix + 1, iy);
-
-		xx = (ix + 1) * hx;
-		yy = (iy + 1) * hy;
-		H(2, 0) = xx;
-		H(2, 1) = yy;
-		H(2, 2) = xx * yy;
-		H(2, 3) = 1.0;
-		b(2) = hv(ix + 1, iy + 1);
-
-		xx = ix * hx;
-		yy = (iy + 1) * hy;
-		H(3, 0) = xx;
-		H(3, 1) = yy;
-		H(3, 2) = xx * yy;
-		H(3, 3) = 1.0;
-		b(3) = hv(ix, iy + 1);
-
-		H.gauss_jordan();
-		H.vmult(X, b);
-
-		val = X(0) * x + X(1) * y + X(2) * x * y + X(3);
-	}
-
-	return val;
-}
+  template 
+    double
+    Input::obstacle_function (
+        double x, double y)
+    {
+      int ix = (int) (x / hx);
+      int iy = (int) (y / hy);
+
+      if (ix < 0)
+        ix = 0;
+
+      if (iy < 0)
+        iy = 0;
+
+      if (ix >= nx - 1)
+        ix = nx - 2;
+
+      if (iy >= ny - 1)
+        iy = ny - 2;
+
+      double val = 0.0;
+        {
+          FullMatrix H(4, 4);
+          Vector X(4);
+          Vector b(4);
+
+          double xx = 0.0;
+          double yy = 0.0;
+
+          xx = ix * hx;
+          yy = iy * hy;
+          H(0, 0) = xx;
+          H(0, 1) = yy;
+          H(0, 2) = xx * yy;
+          H(0, 3) = 1.0;
+          b(0) = hv(ix, iy);
+
+          xx = (ix + 1) * hx;
+          yy = iy * hy;
+          H(1, 0) = xx;
+          H(1, 1) = yy;
+          H(1, 2) = xx * yy;
+          H(1, 3) = 1.0;
+          b(1) = hv(ix + 1, iy);
+
+          xx = (ix + 1) * hx;
+          yy = (iy + 1) * hy;
+          H(2, 0) = xx;
+          H(2, 1) = yy;
+          H(2, 2) = xx * yy;
+          H(2, 3) = 1.0;
+          b(2) = hv(ix + 1, iy + 1);
+
+          xx = ix * hx;
+          yy = (iy + 1) * hy;
+          H(3, 0) = xx;
+          H(3, 1) = yy;
+          H(3, 2) = xx * yy;
+          H(3, 3) = 1.0;
+          b(3) = hv(ix, iy + 1);
+
+          H.gauss_jordan();
+          H.vmult(X, b);
+
+          val = X(0) * x + X(1) * y + X(2) * x * y + X(3);
+        }
+
+      return val;
+    }
 
 // As mentioned above this function reads in the
 // obstacle data and stores them in the std::vector
 // obstacle_data. It will be used only in run ().
-template
-void Input::read_obstacle(const char* name) {
-	std::ifstream f(name);
+  template 
+    void
+    Input::read_obstacle (
+        const char* name)
+    {
+      std::ifstream f(name);
 
-	std::string temp;
-	f >> temp >> nx >> ny;
-	assert(nx > 0 && ny > 0);
+      std::string temp;
+      f >> temp >> nx >> ny;
+      assert(nx > 0 && ny > 0);
 
-	for (int k = 0; k < nx * ny; k++) {
-		double val;
-		f >> val;
-		obstacle_data.push_back(val);
-	}
+      for (int k = 0; k < nx * ny; k++)
+        {
+          double val;
+          f >> val;
+          obstacle_data.push_back(val);
+        }
 
-	hx = 1.0 / (nx - 1);
-	hy = 1.0 / (ny - 1);
+      hx = 1.0 / (nx - 1);
+      hy = 1.0 / (ny - 1);
 
-	pcout << "Resolution of the scanned obstacle picture: " << nx << " x " << ny
-			<< std::endl;
-}
+      pcout << "Resolution of the scanned obstacle picture: " << nx << " x "
+          << ny << std::endl;
+    }
 
 // @sect3{The ConstitutiveLaw class template}
 
@@ -249,41 +272,48 @@ void Input::read_obstacle(const char* name) {
 // isotropic hardening.
 // For gamma = 0 we obtain perfect elastoplastic
 // behavior.
-template
-class ConstitutiveLaw {
-public:
-	ConstitutiveLaw(double _E, double _nu, double _sigma_0, double _gamma,
-			MPI_Comm _mpi_communicator, ConditionalOStream _pcout);
-
-	void
-	plast_linear_hardening(SymmetricTensor<4, dim> &stress_strain_tensor,
-			const SymmetricTensor<2, dim> &strain_tensor,
-			unsigned int &elast_points, unsigned int &plast_points,
-			double &yield);
-	void
-	linearized_plast_linear_hardening(
-			SymmetricTensor<4, dim> &stress_strain_tensor_linearized,
-			SymmetricTensor<4, dim> &stress_strain_tensor,
-			const SymmetricTensor<2, dim> &strain_tensor);
-	inline SymmetricTensor<2, dim>
-	get_strain(const FEValues &fe_values, const unsigned int shape_func,
-			const unsigned int q_point) const;
-	void set_sigma_0(double sigma_hlp) {
-		sigma_0 = sigma_hlp;
-	}
-
-private:
-	SymmetricTensor<4, dim> stress_strain_tensor_mu;
-	SymmetricTensor<4, dim> stress_strain_tensor_kappa;
-	double E;
-	double nu;
-	double sigma_0;
-	double gamma;
-	double mu;
-	double kappa;
-	MPI_Comm mpi_communicator;
-	ConditionalOStream pcout;
-};
+  template 
+    class ConstitutiveLaw
+    {
+      public:
+        ConstitutiveLaw (
+            double _E, double _nu, double _sigma_0, double _gamma,
+            MPI_Comm _mpi_communicator, ConditionalOStream _pcout);
+
+        void
+        plast_linear_hardening (
+            SymmetricTensor<4, dim> &stress_strain_tensor,
+            const SymmetricTensor<2, dim> &strain_tensor,
+            unsigned int &elast_points, unsigned int &plast_points,
+            double &yield);
+        void
+        linearized_plast_linear_hardening (
+            SymmetricTensor<4, dim> &stress_strain_tensor_linearized,
+            SymmetricTensor<4, dim> &stress_strain_tensor,
+            const SymmetricTensor<2, dim> &strain_tensor);
+        inline SymmetricTensor<2, dim>
+        get_strain (
+            const FEValues &fe_values, const unsigned int shape_func,
+            const unsigned int q_point) const;
+        void
+        set_sigma_0 (
+            double sigma_hlp)
+        {
+          sigma_0 = sigma_hlp;
+        }
+
+      private:
+        SymmetricTensor<4, dim> stress_strain_tensor_mu;
+        SymmetricTensor<4, dim> stress_strain_tensor_kappa;
+        double E;
+        double nu;
+        double sigma_0;
+        double gamma;
+        double mu;
+        double kappa;
+        MPI_Comm mpi_communicator;
+        ConditionalOStream pcout;
+    };
 
 // The constructor of the ConstitutiveLaw class sets the
 // required material parameter for our deformable body:
@@ -294,37 +324,46 @@ private:
 // Also it supplies the stress strain tensor of forth order
 // of the volumetric and deviator part. For further details
 // see the documentation above.
-template
-ConstitutiveLaw::ConstitutiveLaw(double _E, double _nu, double _sigma_0,
-		double _gamma, MPI_Comm _mpi_communicator, ConditionalOStream _pcout) :
-		E(_E), nu(_nu), sigma_0(_sigma_0), gamma(_gamma), mpi_communicator(
-				_mpi_communicator), pcout(_pcout) {
-	mu = E / (2 * (1 + nu));
-	kappa = E / (3 * (1 - 2 * nu));
-	stress_strain_tensor_kappa = kappa
-			* outer_product(unit_symmetric_tensor(),
-					unit_symmetric_tensor());
-	stress_strain_tensor_mu = 2 * mu
-			* (identity_tensor()
-					- outer_product(unit_symmetric_tensor(),
-							unit_symmetric_tensor()) / 3.0);
-}
+  template 
+    ConstitutiveLaw::ConstitutiveLaw (
+        double _E, double _nu, double _sigma_0, double _gamma,
+        MPI_Comm _mpi_communicator, ConditionalOStream _pcout)
+        :
+            E(_E),
+            nu(_nu),
+            sigma_0(_sigma_0),
+            gamma(_gamma),
+            mpi_communicator(_mpi_communicator),
+            pcout(_pcout)
+    {
+      mu = E / (2 * (1 + nu));
+      kappa = E / (3 * (1 - 2 * nu));
+      stress_strain_tensor_kappa = kappa
+          * outer_product(unit_symmetric_tensor(),
+              unit_symmetric_tensor());
+      stress_strain_tensor_mu = 2 * mu
+          * (identity_tensor()
+              - outer_product(unit_symmetric_tensor(),
+                  unit_symmetric_tensor()) / 3.0);
+    }
 
 // @sect3{ConstitutiveLaw::ConstitutiveLaw}
 
 // Calculates the strain $\varepsilon(\varphi)=\dfrac{1}{2}\left(\nabla\varphi + \nabla\varphi^T$
 // for the shape functions $\varphi$.
-template
-inline SymmetricTensor<2, dim> ConstitutiveLaw::get_strain(
-		const FEValues &fe_values, const unsigned int shape_func,
-		const unsigned int q_point) const {
-	const FEValuesExtractors::Vector displacement(0);
-	SymmetricTensor < 2, dim > tmp;
+  template 
+    inline SymmetricTensor<2, dim>
+    ConstitutiveLaw::get_strain (
+        const FEValues &fe_values, const unsigned int shape_func,
+        const unsigned int q_point) const
+    {
+      const FEValuesExtractors::Vector displacement(0);
+      SymmetricTensor<2, dim> tmp;
 
-	tmp = fe_values[displacement].symmetric_gradient(shape_func, q_point);
+      tmp = fe_values[displacement].symmetric_gradient(shape_func, q_point);
 
-	return tmp;
-}
+      return tmp;
+    }
 
 // @sect3{ConstitutiveLaw::plast_linear_hardening}
 
@@ -335,35 +374,40 @@ inline SymmetricTensor<2, dim> ConstitutiveLaw::get_strain(
 // points. We need this function to calculate the nonlinear
 // residual in
 // PlasticityContactProblem::residual_nl_system(TrilinosWrappers::MPI::Vector &u).
-template
-void ConstitutiveLaw::plast_linear_hardening(
-		SymmetricTensor<4, dim> &stress_strain_tensor,
-		const SymmetricTensor<2, dim> &strain_tensor,
-		unsigned int &elast_points, unsigned int &plast_points, double &yield) {
-	if (dim == 3) {
-		SymmetricTensor < 2, dim > stress_tensor;
-		stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
-				* strain_tensor;
-
-		SymmetricTensor < 2, dim > deviator_stress_tensor = deviator(
-				stress_tensor);
-
-		double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
-
-		yield = 0;
-		stress_strain_tensor = stress_strain_tensor_mu;
-		double beta = 1.0;
-		if (deviator_stress_tensor_norm > sigma_0) {
-			beta = sigma_0 / deviator_stress_tensor_norm;
-			stress_strain_tensor *= (gamma + (1 - gamma) * beta);
-			yield = 1;
-			plast_points += 1;
-		} else
-			elast_points += 1;
-
-		stress_strain_tensor += stress_strain_tensor_kappa;
-	}
-}
+  template 
+    void
+    ConstitutiveLaw::plast_linear_hardening (
+        SymmetricTensor<4, dim> &stress_strain_tensor,
+        const SymmetricTensor<2, dim> &strain_tensor,
+        unsigned int &elast_points, unsigned int &plast_points, double &yield)
+    {
+      if (dim == 3)
+        {
+          SymmetricTensor<2, dim> stress_tensor;
+          stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
+              * strain_tensor;
+
+          SymmetricTensor<2, dim> deviator_stress_tensor = deviator(
+              stress_tensor);
+
+          double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
+
+          yield = 0;
+          stress_strain_tensor = stress_strain_tensor_mu;
+          double beta = 1.0;
+          if (deviator_stress_tensor_norm > sigma_0)
+            {
+              beta = sigma_0 / deviator_stress_tensor_norm;
+              stress_strain_tensor *= (gamma + (1 - gamma) * beta);
+              yield = 1;
+              plast_points += 1;
+            }
+          else
+            elast_points += 1;
+
+          stress_strain_tensor += stress_strain_tensor_kappa;
+        }
+    }
 
 // @sect3{ConstitutiveLaw::linearized_plast_linear_hardening}
 
@@ -377,116 +421,139 @@ void ConstitutiveLaw::plast_linear_hardening(
 // See
 // PlasticityContactProblem::assemble_nl_system(TrilinosWrappers::MPI::Vector &u)
 // where this function is used.
-template
-void ConstitutiveLaw::linearized_plast_linear_hardening(
-		SymmetricTensor<4, dim> &stress_strain_tensor_linearized,
-		SymmetricTensor<4, dim> &stress_strain_tensor,
-		const SymmetricTensor<2, dim> &strain_tensor) {
-	if (dim == 3) {
-		SymmetricTensor < 2, dim > stress_tensor;
-		stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
-				* strain_tensor;
-
-		SymmetricTensor < 2, dim > deviator_stress_tensor = deviator(
-				stress_tensor);
-
-		double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
-
-		stress_strain_tensor = stress_strain_tensor_mu;
-		stress_strain_tensor_linearized = stress_strain_tensor_mu;
-		double beta = 1.0;
-		if (deviator_stress_tensor_norm > sigma_0) {
-			beta = sigma_0 / deviator_stress_tensor_norm;
-			stress_strain_tensor *= (gamma + (1 - gamma) * beta);
-			stress_strain_tensor_linearized *= (gamma + (1 - gamma) * beta);
-			deviator_stress_tensor /= deviator_stress_tensor_norm;
-			stress_strain_tensor_linearized -= (1 - gamma) * beta * 2 * mu
-					* outer_product(deviator_stress_tensor,
-							deviator_stress_tensor);
-		}
-
-		stress_strain_tensor += stress_strain_tensor_kappa;
-		stress_strain_tensor_linearized += stress_strain_tensor_kappa;
-	}
-}
+  template 
+    void
+    ConstitutiveLaw::linearized_plast_linear_hardening (
+        SymmetricTensor<4, dim> &stress_strain_tensor_linearized,
+        SymmetricTensor<4, dim> &stress_strain_tensor,
+        const SymmetricTensor<2, dim> &strain_tensor)
+    {
+      if (dim == 3)
+        {
+          SymmetricTensor<2, dim> stress_tensor;
+          stress_tensor = (stress_strain_tensor_kappa + stress_strain_tensor_mu)
+              * strain_tensor;
+
+          SymmetricTensor<2, dim> deviator_stress_tensor = deviator(
+              stress_tensor);
+
+          double deviator_stress_tensor_norm = deviator_stress_tensor.norm();
+
+          stress_strain_tensor = stress_strain_tensor_mu;
+          stress_strain_tensor_linearized = stress_strain_tensor_mu;
+          double beta = 1.0;
+          if (deviator_stress_tensor_norm > sigma_0)
+            {
+              beta = sigma_0 / deviator_stress_tensor_norm;
+              stress_strain_tensor *= (gamma + (1 - gamma) * beta);
+              stress_strain_tensor_linearized *= (gamma + (1 - gamma) * beta);
+              deviator_stress_tensor /= deviator_stress_tensor_norm;
+              stress_strain_tensor_linearized -= (1 - gamma) * beta * 2 * mu
+                  * outer_product(deviator_stress_tensor,
+                      deviator_stress_tensor);
+            }
+
+          stress_strain_tensor += stress_strain_tensor_kappa;
+          stress_strain_tensor_linearized += stress_strain_tensor_kappa;
+        }
+    }
 
-namespace EquationData {
+  namespace EquationData
+  {
 // It possible to apply an additional body force
 // but in here it is set to zero.
-template
-class RightHandSide: public Function {
-public:
-	RightHandSide() :
-			Function(dim) {
-	}
-
-	virtual double
-	value(const Point &p, const unsigned int component = 0) const;
-
-	virtual void
-	vector_value(const Point &p, Vector &values) const;
-};
-
-template
-double RightHandSide::value(const Point &p,
-		const unsigned int component) const {
-	double return_value = 0.0;
-
-	if (component == 0)
-		return_value = 0.0;
-	if (component == 1)
-		return_value = 0.0;
-	if (component == 2)
-	  return_value = 0.0;
-
-	return return_value;
-}
+    template 
+      class RightHandSide : public Function
+      {
+        public:
+          RightHandSide ()
+              :
+                  Function(dim)
+          {
+          }
 
-template
-void RightHandSide::vector_value(const Point &p,
-		Vector &values) const {
-	for (unsigned int c = 0; c < this->n_components; ++c)
-		values(c) = RightHandSide::value(p, c);
-}
+          virtual double
+          value (
+              const Point &p, const unsigned int component = 0) const;
+
+          virtual void
+          vector_value (
+              const Point &p, Vector &values) const;
+      };
+
+    template 
+      double
+      RightHandSide::value (
+          const Point &p, const unsigned int component) const
+      {
+        double return_value = 0.0;
+
+        if (component == 0)
+          return_value = 0.0;
+        if (component == 1)
+          return_value = 0.0;
+        if (component == 2)
+          return_value = 0.0;
+
+        return return_value;
+      }
+
+    template 
+      void
+      RightHandSide::vector_value (
+          const Point &p, Vector &values) const
+      {
+        for (unsigned int c = 0; c < this->n_components; ++c)
+          values(c) = RightHandSide::value(p, c);
+      }
 
 // This function class is used to describe the prescribed displacements
 // at the boundary. But again we set this to zero.
-template
-class BoundaryValues: public Function {
-public:
-	BoundaryValues() :
-			Function(dim) {
-	}
-	;
-
-	virtual double
-	value(const Point &p, const unsigned int component = 0) const;
-
-	virtual void
-	vector_value(const Point &p, Vector &values) const;
-};
-
-template
-double BoundaryValues::value(const Point &p,
-		const unsigned int component) const {
-	double return_value = 0;
-
-	if (component == 0)
-		return_value = 0.0;
-	if (component == 1)
-		return_value = 0.0;
-	if (component == 2)
-		return_value = 0.0;
-
-	return return_value;
-}
-
-template
-void BoundaryValues::vector_value(const Point &p,
-		Vector &values) const {
-	for (unsigned int c = 0; c < this->n_components; ++c)
-		values(c) = BoundaryValues::value(p, c);
-}
+    template 
+      class BoundaryValues : public Function
+      {
+        public:
+          BoundaryValues ()
+              :
+                  Function(dim)
+          {
+          }
+          ;
+
+          virtual double
+          value (
+              const Point &p, const unsigned int component = 0) const;
+
+          virtual void
+          vector_value (
+              const Point &p, Vector &values) const;
+      };
+
+    template 
+      double
+      BoundaryValues::value (
+          const Point &p, const unsigned int component) const
+      {
+        double return_value = 0;
+
+        if (component == 0)
+          return_value = 0.0;
+        if (component == 1)
+          return_value = 0.0;
+        if (component == 2)
+          return_value = 0.0;
+
+        return return_value;
+      }
+
+    template 
+      void
+      BoundaryValues::vector_value (
+          const Point &p, Vector &values) const
+      {
+        for (unsigned int c = 0; c < this->n_components; ++c)
+          values(c) = BoundaryValues::value(p, c);
+      }
 
 // This function is obviously implemented to
 // define the obstacle that penetrates our deformable
@@ -494,60 +561,72 @@ void BoundaryValues::vector_value(const Point &p,
 // your obstacle: to read it from a file or to use
 // a function (here a ball).
 // z_max_domain is the z value of the surface of the work piece  
-template
-class Obstacle: public Function {
-public:
-	Obstacle(std_cxx1x::shared_ptr > const &_input,
-		 bool _use_read_obstacle, double z_max_domain) :
-			Function(dim), input_obstacle_copy(_input), use_read_obstacle(
-			  _use_read_obstacle),
-			z_max_domain(z_max_domain){
-	}
-
-	virtual double
-	value(const Point &p, const unsigned int component = 0) const;
-
-	virtual void
-	vector_value(const Point &p, Vector &values) const;
-
-private:
-	std_cxx1x::shared_ptr > const &input_obstacle_copy;
-	bool use_read_obstacle;
-    double z_max_domain;
-};
-
-template
-double Obstacle::value(const Point &p,
-		const unsigned int component) const {
-	if (component == 0)
-	  return p(0);
-	if (component == 1)
-	  return p(1);
-
-	//component==2:
-	if (use_read_obstacle)
-	  {
-	    if (p(0) >= 0.0 && p(0) <= 1.0 && p(1) >= 0.0 && p(1) <= 1.0)
-	      return z_max_domain + 0.999 - input_obstacle_copy->obstacle_function(p(0), p(1));
-	    else
-	      return 10000.0;
-	  }
-	else
-	  {
-	    //sphere:
-	    return -std::sqrt(
-	      0.36 - (p(0) - 0.5) * (p(0) - 0.5)
-	      - (p(1) - 0.5) * (p(1) - 0.5)) + z_max_domain + 0.59;
-	  }
-}
+    template 
+      class Obstacle : public Function
+      {
+        public:
+          Obstacle (
+              std_cxx1x::shared_ptr > const &_input,
+              bool _use_read_obstacle, double z_max_domain)
+              :
+                  Function(dim),
+                  input_obstacle_copy(_input),
+                  use_read_obstacle(_use_read_obstacle),
+                  z_max_domain(z_max_domain)
+          {
+          }
 
-template
-void Obstacle::vector_value(const Point &p,
-		Vector &values) const {
-	for (unsigned int c = 0; c < this->n_components; ++c)
-		values(c) = Obstacle::value(p, c);
-}
-}
+          virtual double
+          value (
+              const Point &p, const unsigned int component = 0) const;
+
+          virtual void
+          vector_value (
+              const Point &p, Vector &values) const;
+
+        private:
+          std_cxx1x::shared_ptr > const &input_obstacle_copy;
+          bool use_read_obstacle;
+          double z_max_domain;
+      };
+
+    template 
+      double
+      Obstacle::value (
+          const Point &p, const unsigned int component) const
+      {
+        if (component == 0)
+          return p(0);
+        if (component == 1)
+          return p(1);
+
+        //component==2:
+        if (use_read_obstacle)
+          {
+            if (p(0) >= 0.0 && p(0) <= 1.0 && p(1) >= 0.0 && p(1) <= 1.0)
+              return z_max_domain + 0.999
+                  - input_obstacle_copy->obstacle_function(p(0), p(1));
+            else
+              return 10000.0;
+          }
+        else
+          {
+            //sphere:
+            return -std::sqrt(
+                0.36 - (p(0) - 0.5) * (p(0) - 0.5)
+                    - (p(1) - 0.5) * (p(1) - 0.5)) + z_max_domain + 0.59;
+          }
+      }
+
+    template 
+      void
+      Obstacle::vector_value (
+          const Point &p, Vector &values) const
+      {
+        for (unsigned int c = 0; c < this->n_components; ++c)
+          values(c) = Obstacle::value(p, c);
+      }
+  }
 
 // @sect3{The PlasticityContactProblem class template}
 
@@ -566,116 +645,124 @@ void Obstacle::vector_value(const Point &p,
 // situation and to handle the nonlinear
 // operator for the constitutive law.
 
-template
-class PlasticityContactProblem {
-public:
-	PlasticityContactProblem(const ParameterHandler &prm);
-	void
-	run();
-
-	static void
-	declare(ParameterHandler &prm);
-
-private:
-	void
-	make_grid();
-	void
-	setup_system();
-	void
-	assemble_nl_system(TrilinosWrappers::MPI::Vector &u);
-	void
-	residual_nl_system(TrilinosWrappers::MPI::Vector &u);
-	void
-	assemble_mass_matrix_diagonal(TrilinosWrappers::SparseMatrix &mass_matrix);
-	void
-	update_solution_and_constraints();
-	void
-	dirichlet_constraints();
-	void
-	solve();
-	void
-	solve_newton();
-	void
-	refine_grid();
-	void
-	move_mesh(
-			const TrilinosWrappers::MPI::Vector &_complete_displacement) const;
-	void
-	output_results(const std::string &title);
-	void
-	output_contact_force(const unsigned int cycle);
-
-	double to_refine_factor;
-	double to_coarsen_factor;
-	unsigned int cycle;
-
-	MPI_Comm mpi_communicator;
-
-	parallel::distributed::Triangulation triangulation;
-
-        FE_Q     u;
-	FESystem fe;
-	DoFHandler dof_handler;
-
-	// We are using the SolutionTransfer class to interpolate the
-	// solution on the new refined mesh. It appears in th refine_grid()
-	// and the run() function.
-	std_cxx1x::shared_ptr<
-			parallel::distributed::SolutionTransfer > soltrans;
-
-	IndexSet locally_owned_dofs;
-	IndexSet locally_relevant_dofs;
-
-	unsigned int number_iterations;
-
-	ConstraintMatrix constraints;
-	ConstraintMatrix constraints_hanging_nodes;
-	ConstraintMatrix constraints_dirichlet_hanging_nodes;
-
-	TrilinosWrappers::SparseMatrix system_matrix_newton;
-
-	TrilinosWrappers::MPI::Vector solution;
-	TrilinosWrappers::MPI::Vector system_rhs_newton;
-	TrilinosWrappers::MPI::Vector system_rhs_lambda;
-	TrilinosWrappers::MPI::Vector resid_vector;
-	TrilinosWrappers::MPI::Vector diag_mass_matrix_vector;
-	Vector cell_constitution;
-	IndexSet active_set;
-
-	ConditionalOStream pcout;
-
-	TrilinosWrappers::PreconditionAMG::AdditionalData additional_data;
-	TrilinosWrappers::PreconditionAMG preconditioner_u;
-
-	std_cxx1x::shared_ptr > input_obstacle;
-	std_cxx1x::shared_ptr > plast_lin_hard;
-
-	double sigma_0; // Yield stress
-	double gamma; // Parameter for the linear isotropic hardening
-	double e_modul; // E-Modul
-	double nu; // Poisson ratio
-
-	TimerOutput computing_timer;
-
-	unsigned int degree;
-	unsigned int n_initial_refinements;
-	struct RefinementStrategy
-	{
-	  enum value
-	  {
-	    refine_global,
-	    refine_percentage,
-	    refine_fix_dofs
-	  };
-	};
-	typename RefinementStrategy::value refinement_strategy;
-	unsigned int n_cycles;
-	std::string obstacle_filename;
-	std::string output_dir;
-    bool transfer_solution;
-    std::string base_mesh;
-};
+  template 
+    class PlasticityContactProblem
+    {
+      public:
+        PlasticityContactProblem (
+            const ParameterHandler &prm);
+        void
+        run ();
+
+        static void
+        declare (
+            ParameterHandler &prm);
+
+      private:
+        void
+        make_grid ();
+        void
+        setup_system ();
+        void
+        assemble_nl_system (
+            TrilinosWrappers::MPI::Vector &u);
+        void
+        residual_nl_system (
+            TrilinosWrappers::MPI::Vector &u);
+        void
+        assemble_mass_matrix_diagonal (
+            TrilinosWrappers::SparseMatrix &mass_matrix);
+        void
+        update_solution_and_constraints ();
+        void
+        dirichlet_constraints ();
+        void
+        solve ();
+        void
+        solve_newton ();
+        void
+        refine_grid ();
+        void
+        move_mesh (
+            const TrilinosWrappers::MPI::Vector &_complete_displacement) const;
+        void
+        output_results (
+            const std::string &title);
+        void
+        output_contact_force (
+            const unsigned int cycle);
+
+        double to_refine_factor;
+        double to_coarsen_factor;
+        unsigned int cycle;
+
+        MPI_Comm mpi_communicator;
+
+        parallel::distributed::Triangulation triangulation;
+
+        FE_Q u;
+        FESystem fe;
+        DoFHandler dof_handler;
+
+        // We are using the SolutionTransfer class to interpolate the
+        // solution on the new refined mesh. It appears in th refine_grid()
+        // and the run() function.
+        std_cxx1x::shared_ptr<
+            parallel::distributed::SolutionTransfer > soltrans;
+
+        IndexSet locally_owned_dofs;
+        IndexSet locally_relevant_dofs;
+
+        unsigned int number_iterations;
+
+        ConstraintMatrix constraints;
+        ConstraintMatrix constraints_hanging_nodes;
+        ConstraintMatrix constraints_dirichlet_hanging_nodes;
+
+        TrilinosWrappers::SparseMatrix system_matrix_newton;
+
+        TrilinosWrappers::MPI::Vector solution;
+        TrilinosWrappers::MPI::Vector system_rhs_newton;
+        TrilinosWrappers::MPI::Vector system_rhs_lambda;
+        TrilinosWrappers::MPI::Vector resid_vector;
+        TrilinosWrappers::MPI::Vector diag_mass_matrix_vector;
+        Vector cell_constitution;
+        IndexSet active_set;
+
+        ConditionalOStream pcout;
+
+        TrilinosWrappers::PreconditionAMG::AdditionalData additional_data;
+        TrilinosWrappers::PreconditionAMG preconditioner_u;
+
+        std_cxx1x::shared_ptr > input_obstacle;
+        std_cxx1x::shared_ptr > plast_lin_hard;
+
+        double sigma_0; // Yield stress
+        double gamma; // Parameter for the linear isotropic hardening
+        double e_modul; // E-Modul
+        double nu; // Poisson ratio
+
+        TimerOutput computing_timer;
+
+        unsigned int degree;
+        unsigned int n_initial_refinements;
+        struct RefinementStrategy
+        {
+            enum value
+            {
+              refine_global,
+              refine_percentage,
+              refine_fix_dofs
+            };
+        };
+        typename RefinementStrategy::value refinement_strategy;
+        unsigned int n_cycles;
+        std::string obstacle_filename;
+        std::string output_dir;
+        bool transfer_solution;
+        std::string base_mesh;
+    };
 
 // @sect3{Implementation of the PlasticityContactProblem class}
 
@@ -683,655 +770,657 @@ private:
 // template that makes use of the functions
 // above. As before, we will write everything
 
-template
-PlasticityContactProblem::PlasticityContactProblem(
-    const ParameterHandler &prm) :
-		mpi_communicator(MPI_COMM_WORLD), triangulation(mpi_communicator),
-		u(QGaussLobatto< 1 > (prm.get_integer("polynomial degree")+1)),
-		fe(u, dim),
-		dof_handler(triangulation), pcout(
-				std::cout,
-				(Utilities::MPI::this_mpi_process(mpi_communicator) == 0)), sigma_0(
-				400.0), gamma(0.01), e_modul(2.0e+5), nu(0.3), computing_timer(
-				MPI_COMM_WORLD, pcout, TimerOutput::never,
-				TimerOutput::wall_times) {
-	// double _E, double _nu, double _sigma_0, double _gamma
-	plast_lin_hard.reset(
-			new ConstitutiveLaw(e_modul, nu, sigma_0, gamma,
-					mpi_communicator, pcout));
-
-	degree = prm.get_integer("polynomial degree");
-	n_initial_refinements = prm.get_integer("number of initial refinements");
-	std::string strat = prm.get("refinement strategy");
-	if (strat == "global")
-	  refinement_strategy = RefinementStrategy::refine_global;
-	else if (strat == "percentage")
-          refinement_strategy = RefinementStrategy::refine_percentage;
-	else if (strat == "fix dofs")
-          refinement_strategy = RefinementStrategy::refine_fix_dofs;
-	else
-	  throw ExcNotImplemented();
-
-	n_cycles = prm.get_integer("number of cycles");
-	obstacle_filename = prm.get("obstacle filename");
-	output_dir = prm.get("output directory");
-	if (output_dir!="" && *(output_dir.rbegin())!='/')
-	  output_dir += "/";
-	mkdir(output_dir.c_str(), 0777);
-
-	transfer_solution = prm.get_bool("transfer solution");
-	base_mesh = prm.get("base mesh");
-	
-        pcout << "    Using output directory '" << output_dir << "'" << std::endl;
-        pcout << "    FE degree " << degree << std::endl;
-        pcout << "    Obstacle '" << obstacle_filename << "'" << std::endl;
-	pcout << "    transfer solution " << (transfer_solution?"true":"false") << std::endl;
-}
+  template 
+    PlasticityContactProblem::PlasticityContactProblem (
+        const ParameterHandler &prm)
+        :
+            mpi_communicator(MPI_COMM_WORLD),
+            triangulation(mpi_communicator),
+            u(QGaussLobatto<1>(prm.get_integer("polynomial degree") + 1)),
+            fe(u, dim),
+            dof_handler(triangulation),
+            pcout(std::cout,
+                (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)),
+            sigma_0(400.0),
+            gamma(0.01),
+            e_modul(2.0e+5),
+            nu(0.3),
+            computing_timer(MPI_COMM_WORLD, pcout, TimerOutput::never,
+                TimerOutput::wall_times)
+    {
+      // double _E, double _nu, double _sigma_0, double _gamma
+      plast_lin_hard.reset(
+          new ConstitutiveLaw(e_modul, nu, sigma_0, gamma,
+              mpi_communicator, pcout));
+
+      degree = prm.get_integer("polynomial degree");
+      n_initial_refinements = prm.get_integer("number of initial refinements");
+      std::string strat = prm.get("refinement strategy");
+      if (strat == "global")
+        refinement_strategy = RefinementStrategy::refine_global;
+      else if (strat == "percentage")
+        refinement_strategy = RefinementStrategy::refine_percentage;
+      else
+        throw ExcNotImplemented();
+
+      n_cycles = prm.get_integer("number of cycles");
+      obstacle_filename = prm.get("obstacle filename");
+      output_dir = prm.get("output directory");
+      if (output_dir != "" && *(output_dir.rbegin()) != '/')
+        output_dir += "/";
+      mkdir(output_dir.c_str(), 0777);
+
+      transfer_solution = prm.get_bool("transfer solution");
+      base_mesh = prm.get("base mesh");
+
+      pcout << "    Using output directory '" << output_dir << "'" << std::endl;
+      pcout << "    FE degree " << degree << std::endl;
+      pcout << "    Obstacle '" << obstacle_filename << "'" << std::endl;
+      pcout << "    transfer solution "
+          << (transfer_solution ? "true" : "false") << std::endl;
+    }
 
 // @sect4{PlasticityContactProblem::declare}
 
-template 
-void
-PlasticityContactProblem::declare(ParameterHandler &prm)
-{
-  prm.declare_entry("polynomial degree","1",Patterns::Integer(),"polynomial degree of the FE_Q finite element space, typically 1 or 2");
-  prm.declare_entry("number of initial refinements","2",Patterns::Integer(),"number of initial global refinements before the first computation");
-  prm.declare_entry("refinement strategy","percentage",Patterns::Selection("global|percentage|fix dofs"),
-      "refinement strategy for each cycle:\n"
-      " global: one global refinement\n"
-      "percentage: fixed percentage gets refined using kelly\n"
-      " fix dofs: tries to achieve 2^initial_refinement*300 dofs after cycle 1 (only use 2 cycles!). Changes the coarse mesh!");
-  prm.declare_entry("number of cycles","5",Patterns::Integer(),"number of adaptive cycles to run");
-  prm.declare_entry("obstacle filename","",Patterns::Anything(),"obstacle file to read, use 'obstacle_file.pbm' or leave empty to use a sphere");
-  prm.declare_entry("output directory","",Patterns::Anything(),"directory to put output files (graphical output and benchmark statistics), leave empty to put into current directory");
-  prm.declare_entry("transfer solution","false",Patterns::Bool(),"decide if the solution should be used as a starting guess for the finer mesh, use 0 otherwise.");
-  prm.declare_entry("base mesh","box",Patterns::Selection("box|half sphere"),
-		    "select the shape of the work piece: 'box' or 'half sphere'");
-  
-}
+  template 
+    void
+    PlasticityContactProblem::declare (
+        ParameterHandler &prm)
+    {
+      prm.declare_entry("polynomial degree", "1", Patterns::Integer(),
+          "polynomial degree of the FE_Q finite element space, typically 1 or 2");
+      prm.declare_entry("number of initial refinements", "2",
+          Patterns::Integer(),
+          "number of initial global refinements before the first computation");
+      prm.declare_entry("refinement strategy", "percentage",
+          Patterns::Selection("global|percentage|fix dofs"),
+          "refinement strategy for each cycle:\n"
+              " global: one global refinement\n"
+              "percentage: fixed percentage gets refined using kelly\n"
+              " fix dofs: tries to achieve 2^initial_refinement*300 dofs after cycle 1 (only use 2 cycles!). Changes the coarse mesh!");
+      prm.declare_entry("number of cycles", "5", Patterns::Integer(),
+          "number of adaptive cycles to run");
+      prm.declare_entry("obstacle filename", "", Patterns::Anything(),
+          "obstacle file to read, use 'obstacle_file.pbm' or leave empty to use a sphere");
+      prm.declare_entry("output directory", "", Patterns::Anything(),
+          "directory to put output files (graphical output and benchmark statistics), leave empty to put into current directory");
+      prm.declare_entry("transfer solution", "false", Patterns::Bool(),
+          "decide if the solution should be used as a starting guess for the finer mesh, use 0 otherwise.");
+      prm.declare_entry("base mesh", "box",
+          Patterns::Selection("box|half sphere"),
+          "select the shape of the work piece: 'box' or 'half sphere'");
+
+    }
 
   Point<3>
-  rotate_half_sphere(const Point<3> &in)
+  rotate_half_sphere (
+      const Point<3> &in)
   {
     return Point<3>(in(2), in(1), -in(0));
   }
 
 // @sect4{PlasticityContactProblem::make_grid}
-  
-template
-void PlasticityContactProblem::make_grid() {
 
-  if (base_mesh == "half sphere")
+  template 
+    void
+    PlasticityContactProblem::make_grid ()
     {
-      Point < dim > center(0, 0, 0);
-      double radius = 0.8;
-      GridGenerator::half_hyper_ball(triangulation, center, radius);
-      GridTools::transform(&rotate_half_sphere, triangulation);
-      Point < dim > shift(0.5, 0.5, 0.5);
-      GridTools::shift(shift, triangulation);
-      static PolarManifold boundary_description(Point(0.5,0.5,0.5));
-      triangulation.set_manifold (0, boundary_description);
-
-      triangulation.refine_global(n_initial_refinements);
+      if (base_mesh == "half sphere")
+        {
+          Point center(0, 0, 0);
+          double radius = 0.8;
+          GridGenerator::half_hyper_ball(triangulation, center, radius);
+          GridTools::transform(&rotate_half_sphere, triangulation);
+          Point shift(0.5, 0.5, 0.5);
+          GridTools::shift(shift, triangulation);
+          static PolarManifold boundary_description(
+              Point(0.5, 0.5, 0.5));
+          triangulation.set_manifold(0, boundary_description);
+
+          triangulation.refine_global(n_initial_refinements);
+
+          to_refine_factor = 0.3;
+          to_coarsen_factor = 0.03;
+          return;
+        }
+
+      Point p1(0, 0, 0);
+      Point p2(1.0, 1.0, 1.0);
 
+      GridGenerator::hyper_rectangle(triangulation, p1, p2);
       to_refine_factor = 0.3;
       to_coarsen_factor = 0.03;
-      return;
+
+      Triangulation<3>::active_cell_iterator cell =
+          triangulation.begin_active(), endc = triangulation.end();
+
+      /* boundary_indicators:
+           _______
+         /  1    /|
+        /______ / |
+       |       | 8|
+       |   8   | /
+       |_______|/
+           6
+
+       The boundary indicators of the sides of the cube are 8.
+       The boundary indicator of the bottom is indicated with 6
+       and the top with 1.
+       */
+
+      for (; cell != endc; ++cell)
+        for (unsigned int face = 0; face < GeometryInfo::faces_per_cell;
+            ++face)
+          {
+            if (cell->face(face)->center()[2] == p2(2))
+              cell->face(face)->set_boundary_indicator(1);
+            if (cell->face(face)->center()[0] == p1(0)
+                || cell->face(face)->center()[0] == p2(0)
+                || cell->face(face)->center()[1] == p1(1)
+                || cell->face(face)->center()[1] == p2(1))
+              cell->face(face)->set_boundary_indicator(8);
+            if (cell->face(face)->center()[2] == p1(2))
+              cell->face(face)->set_boundary_indicator(6);
+          }
+
+      triangulation.refine_global(n_initial_refinements);
     }
-  
-  Point < dim > p1(0, 0, 0);
-  Point < dim > p2(1.0, 1.0, 1.0);
-  unsigned int ref = n_initial_refinements;
 
-  if (refinement_strategy == RefinementStrategy::refine_fix_dofs)
+  template 
+    void
+    PlasticityContactProblem::setup_system ()
     {
-	/**
-	 * This complicated logic creates a mesh and a refinement fraction to_refine_factor,
-	 * so that the resulting mesh after adaptive refinement has approximately
-	 * 2^n_refinements_global*300 dofs. This allows parallel scalability tests.
-	 * About 5%-10% of the cells are being adaptively refined.
-	 * We start with a 3x3,4x4, or 5x5 base mesh (whichever is closed in cell
-	 * count).
-	 */
-	unsigned int ref = (n_initial_refinements + 1) / 3;
-	unsigned int remain = n_initial_refinements + 1 - ref * 3;
-	unsigned int rep = 3;
-	if (remain == 1)
-		rep = 4;
-	else if (remain == 2)
-		rep = 5;
-
-	unsigned int n_cells_x = (1 << ref) * rep;
-	unsigned int goal_dofs = (1 << n_initial_refinements) * 300;
-	double goal_cells = std::pow(std::pow(goal_dofs / 3.0, 1.0 / 3.0) - 1.0,
-			3.0);
-	double n_cells = std::pow(n_cells_x, 3.0);
-	to_refine_factor = (goal_cells - n_cells) / n_cells;
-	//convert from fraction of cells to add to fraction of cells to refine:
-	to_refine_factor /= 7.0;
-	to_coarsen_factor = 0.0;
-
-	std::vector repet(3);
-	repet[0] = rep;
-	repet[1] = rep;
-	repet[2] = rep;
-
-	GridGenerator::subdivided_hyper_rectangle(triangulation, repet, p1, p2);
+      // setup dofs
+        {
+          TimerOutput::Scope t(computing_timer, "Setup: distribute DoFs");
+          dof_handler.distribute_dofs(fe);
+
+          locally_owned_dofs = dof_handler.locally_owned_dofs();
+          locally_relevant_dofs.clear();
+          DoFTools::extract_locally_relevant_dofs(dof_handler,
+              locally_relevant_dofs);
+        }
+
+      // setup hanging nodes and dirichlet constraints
+        {
+          TimerOutput::Scope t(computing_timer, "Setup: constraints");
+          constraints_hanging_nodes.reinit(locally_relevant_dofs);
+          DoFTools::make_hanging_node_constraints(dof_handler,
+              constraints_hanging_nodes);
+          constraints_hanging_nodes.close();
+
+          pcout << "   Number of active cells: "
+              << triangulation.n_global_active_cells() << std::endl
+              << "   Number of degrees of freedom: " << dof_handler.n_dofs()
+              << std::endl;
+
+          dirichlet_constraints();
+        }
+
+      // Initialization for matrices and vectors
+        {
+          TimerOutput::Scope t(computing_timer, "Setup: vectors");
+          solution.reinit(locally_relevant_dofs, mpi_communicator);
+          system_rhs_newton.reinit(locally_owned_dofs, mpi_communicator);
+          system_rhs_lambda.reinit(system_rhs_newton);
+          resid_vector.reinit(system_rhs_newton);
+          diag_mass_matrix_vector.reinit(system_rhs_newton);
+          cell_constitution.reinit(triangulation.n_active_cells());
+          active_set.clear();
+          active_set.set_size(locally_relevant_dofs.size());
+        }
+
+      // setup sparsity pattern
+        {
+          TimerOutput::Scope t(computing_timer, "Setup: matrix");
+          TrilinosWrappers::SparsityPattern sp(locally_owned_dofs,
+              mpi_communicator);
+
+          DoFTools::make_sparsity_pattern(dof_handler, sp,
+              constraints_dirichlet_hanging_nodes, false,
+              Utilities::MPI::this_mpi_process(mpi_communicator));
+
+          sp.compress();
+
+          system_matrix_newton.reinit(sp);
+
+          // we are going to reuse the system
+          // matrix for assembling the diagonal
+          // of the mass matrix so that we do not
+          // need to allocate two sparse matrices
+          // at the same time:
+          TrilinosWrappers::SparseMatrix & mass_matrix = system_matrix_newton;
+          assemble_mass_matrix_diagonal(mass_matrix);
+          const unsigned int start = (system_rhs_newton.local_range().first),
+              end = (system_rhs_newton.local_range().second);
+          for (unsigned int j = start; j < end; j++)
+            diag_mass_matrix_vector(j) = mass_matrix.diag_element(j);
+
+          number_iterations = 0;
+
+          diag_mass_matrix_vector.compress(VectorOperation::insert);
+
+          // remove the mass matrix entries from the matrix:
+          mass_matrix = 0;
+        }
     }
-  else
+
+  template 
+    void
+    PlasticityContactProblem::assemble_nl_system (
+        TrilinosWrappers::MPI::Vector &u)
     {
-      GridGenerator::hyper_rectangle(triangulation, p1, p2);
-      to_refine_factor = 0.3;
-      to_coarsen_factor = 0.03;
+      TimerOutput::Scope t(computing_timer, "Assembling");
+
+      QGauss quadrature_formula(fe.degree + 1);
+      QGauss face_quadrature_formula(fe.degree + 1);
+
+      FEValues fe_values(fe, quadrature_formula,
+          UpdateFlags(
+              update_values | update_gradients | update_q_points
+                  | update_JxW_values));
+
+      FEFaceValues fe_values_face(fe, face_quadrature_formula,
+          update_values | update_quadrature_points | update_JxW_values);
+
+      const unsigned int dofs_per_cell = fe.dofs_per_cell;
+      const unsigned int n_q_points = quadrature_formula.size();
+      const unsigned int n_face_q_points = face_quadrature_formula.size();
+
+      const EquationData::RightHandSide right_hand_side;
+      std::vector > right_hand_side_values(n_q_points,
+          Vector(dim));
+      std::vector > right_hand_side_values_face(n_face_q_points,
+          Vector(dim));
+
+      FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell);
+      Vector cell_rhs(dofs_per_cell);
+
+      std::vector local_dof_indices(dofs_per_cell);
+
+      typename DoFHandler::active_cell_iterator cell =
+          dof_handler.begin_active(), endc = dof_handler.end();
+
+      const FEValuesExtractors::Vector displacement(0);
+
+      const double kappa = 1.0;
+      for (; cell != endc; ++cell)
+        if (cell->is_locally_owned())
+          {
+            fe_values.reinit(cell);
+            cell_matrix = 0;
+            cell_rhs = 0;
+
+            right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
+                right_hand_side_values);
+
+            std::vector > strain_tensor(n_q_points);
+            fe_values[displacement].get_function_symmetric_gradients(u,
+                strain_tensor);
+
+            for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+              {
+                SymmetricTensor<4, dim> stress_strain_tensor_linearized;
+                SymmetricTensor<4, dim> stress_strain_tensor;
+                SymmetricTensor<2, dim> stress_tensor;
+
+                plast_lin_hard->linearized_plast_linear_hardening(
+                    stress_strain_tensor_linearized, stress_strain_tensor,
+                    strain_tensor[q_point]);
+
+                for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                  {
+                    stress_tensor = stress_strain_tensor_linearized
+                        * plast_lin_hard->get_strain(fe_values, i, q_point);
+
+                    for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                      {
+                        cell_matrix(i, j) += (stress_tensor
+                            * plast_lin_hard->get_strain(fe_values, j, q_point)
+                            * fe_values.JxW(q_point));
+                      }
+
+                    // the linearized part a(v^i;v^i,v) of the rhs
+                    cell_rhs(i) += (stress_tensor * strain_tensor[q_point]
+                        * fe_values.JxW(q_point));
+
+                    // the residual part a(v^i;v) of the rhs
+                    cell_rhs(i) -= (strain_tensor[q_point]
+                        * stress_strain_tensor
+                        * plast_lin_hard->get_strain(fe_values, i, q_point)
+                        * fe_values.JxW(q_point));
+
+                    // the residual part F(v) of the rhs
+                    Tensor<1, dim> rhs_values;
+                    rhs_values = 0;
+                    cell_rhs(i) += (fe_values[displacement].value(i, q_point)
+                        * rhs_values * fe_values.JxW(q_point));
+                  }
+              }
+
+            for (unsigned int face = 0;
+                face < GeometryInfo::faces_per_cell; ++face)
+              {
+                if (cell->face(face)->at_boundary()
+                    && cell->face(face)->boundary_indicator() == 1)
+                  {
+                    fe_values_face.reinit(cell, face);
+
+                    right_hand_side.vector_value_list(
+                        fe_values_face.get_quadrature_points(),
+                        right_hand_side_values_face);
+
+                    for (unsigned int q_point = 0; q_point < n_face_q_points;
+                        ++q_point)
+                      {
+                        Tensor<1, dim> rhs_values;
+                        rhs_values[2] = right_hand_side_values[q_point][2];
+                        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                          cell_rhs(i) += (fe_values_face[displacement].value(i,
+                              q_point) * rhs_values
+                              * fe_values_face.JxW(q_point));
+                      }
+                  }
+              }
+
+            cell->get_dof_indices(local_dof_indices);
+            constraints.distribute_local_to_global(cell_matrix, cell_rhs,
+                local_dof_indices, system_matrix_newton, system_rhs_newton,
+                true);
+
+          };
+
+      system_matrix_newton.compress(VectorOperation::add);
+      system_rhs_newton.compress(VectorOperation::add);
     }
 
+  template 
+    void
+    PlasticityContactProblem::residual_nl_system (
+        TrilinosWrappers::MPI::Vector &u)
+    {
+      QGauss quadrature_formula(fe.degree + 1);
+      QGauss face_quadrature_formula(fe.degree + 1);
+
+      FEValues fe_values(fe, quadrature_formula,
+          UpdateFlags(
+              update_values | update_gradients | update_q_points
+                  | update_JxW_values));
+
+      FEFaceValues fe_values_face(fe, face_quadrature_formula,
+          update_values | update_quadrature_points | update_JxW_values);
+
+      const unsigned int dofs_per_cell = fe.dofs_per_cell;
+      const unsigned int n_q_points = quadrature_formula.size();
+      const unsigned int n_face_q_points = face_quadrature_formula.size();
+
+      const EquationData::RightHandSide right_hand_side;
+      std::vector > right_hand_side_values(n_q_points,
+          Vector(dim));
+      std::vector > right_hand_side_values_face(n_face_q_points,
+          Vector(dim));
+
+      Vector cell_rhs(dofs_per_cell);
+
+      std::vector local_dof_indices(dofs_per_cell);
+
+      const FEValuesExtractors::Vector displacement(0);
+
+      typename DoFHandler::active_cell_iterator cell =
+          dof_handler.begin_active(), endc = dof_handler.end();
+
+      unsigned int elast_points = 0;
+      unsigned int plast_points = 0;
+      double yield = 0;
+      unsigned int cell_number = 0;
+      cell_constitution = 0;
+
+      for (; cell != endc; ++cell)
+        if (cell->is_locally_owned())
+          {
+            fe_values.reinit(cell);
+            cell_rhs = 0;
+
+            right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
+                right_hand_side_values);
+
+            std::vector > strain_tensor(n_q_points);
+            fe_values[displacement].get_function_symmetric_gradients(u,
+                strain_tensor);
+
+            for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+              {
+                SymmetricTensor<4, dim> stress_strain_tensor;
+                SymmetricTensor<2, dim> stress_tensor;
+
+                plast_lin_hard->plast_linear_hardening(stress_strain_tensor,
+                    strain_tensor[q_point], elast_points, plast_points, yield);
+
+                cell_constitution(cell_number) += yield;
+                for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                  {
+                    cell_rhs(i) -= (strain_tensor[q_point]
+                        * stress_strain_tensor
+                        * //(stress_tensor) *
+                        plast_lin_hard->get_strain(fe_values, i, q_point)
+                        * fe_values.JxW(q_point));
+
+                    Tensor<1, dim> rhs_values;
+                    rhs_values = 0;
+                    cell_rhs(i) += ((fe_values[displacement].value(i, q_point)
+                        * rhs_values) * fe_values.JxW(q_point));
+                  };
+              };
+
+            for (unsigned int face = 0;
+                face < GeometryInfo::faces_per_cell; ++face)
+              {
+                if (cell->face(face)->at_boundary()
+                    && cell->face(face)->boundary_indicator() == 1)
+                  {
+                    fe_values_face.reinit(cell, face);
+
+                    right_hand_side.vector_value_list(
+                        fe_values_face.get_quadrature_points(),
+                        right_hand_side_values_face);
+
+                    for (unsigned int q_point = 0; q_point < n_face_q_points;
+                        ++q_point)
+                      {
+                        Tensor<1, dim> rhs_values;
+                        rhs_values[2] = right_hand_side_values[q_point][2];
+                        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                          cell_rhs(i) += (fe_values_face[displacement].value(i,
+                              q_point) * rhs_values
+                              * fe_values_face.JxW(q_point));
+                      }
+                  }
+              }
+
+            cell->get_dof_indices(local_dof_indices);
+            constraints_dirichlet_hanging_nodes.distribute_local_to_global(
+                cell_rhs, local_dof_indices, system_rhs_newton);
+
+            for (unsigned int i = 0; i < dofs_per_cell; i++)
+              system_rhs_lambda(local_dof_indices[i]) += cell_rhs(i);
+
+            cell_number += 1;
+          }
+        else
+          {
+            cell_constitution(cell_number) = 0;
+            cell_number += 1;
+          };
 
-	Triangulation<3>::active_cell_iterator cell = triangulation.begin_active(),
-			endc = triangulation.end();
-
-	/* boundary_indicators:
-	 _______
-	 /  1    /|
-	 /______ / |
-	 8|       | 8|
-	 |   8   | /
-	 |_______|/
-	 6
-	 */
-
-	for (; cell != endc; ++cell)
-		for (unsigned int face = 0;
-				face < GeometryInfo < dim > ::faces_per_cell; ++face) {
-			if (cell->face(face)->center()[2] == p2(2))
-				cell->face(face)->set_boundary_indicator(1);
-			if (cell->face(face)->center()[0] == p1(0)
-					|| cell->face(face)->center()[0] == p2(0)
-					|| cell->face(face)->center()[1] == p1(1)
-					|| cell->face(face)->center()[1] == p2(1))
-				cell->face(face)->set_boundary_indicator(8);
-			if (cell->face(face)->center()[2] == p1(2))
-				cell->face(face)->set_boundary_indicator(6);
-		}
-
-	triangulation.refine_global(ref);
-}
-
-template
-void PlasticityContactProblem::setup_system() {
-	// setup dofs
-	{
-		TimerOutput::Scope t(computing_timer, "Setup: distribute DoFs");
-		dof_handler.distribute_dofs(fe);
-
-		locally_owned_dofs = dof_handler.locally_owned_dofs();
-		locally_relevant_dofs.clear();
-		DoFTools::extract_locally_relevant_dofs(dof_handler,
-				locally_relevant_dofs);
-	}
-
-	// setup hanging nodes and dirichlet constraints
-	{
-		TimerOutput::Scope t(computing_timer, "Setup: constraints");
-		constraints_hanging_nodes.reinit(locally_relevant_dofs);
-		DoFTools::make_hanging_node_constraints(dof_handler,
-				constraints_hanging_nodes);
-		constraints_hanging_nodes.close();
-
-		pcout << "   Number of active cells: "
-				<< triangulation.n_global_active_cells() << std::endl
-				<< "   Number of degrees of freedom: " << dof_handler.n_dofs()
-				<< std::endl;
-
-		dirichlet_constraints();
-	}
-
-	// Initialization for matrices and vectors
-	{
-		TimerOutput::Scope t(computing_timer, "Setup: vectors");
-		solution.reinit(locally_relevant_dofs, mpi_communicator);
-		system_rhs_newton.reinit(locally_owned_dofs, mpi_communicator);
-		system_rhs_lambda.reinit(system_rhs_newton);
-		resid_vector.reinit(system_rhs_newton);
-		diag_mass_matrix_vector.reinit(system_rhs_newton);
-		cell_constitution.reinit(triangulation.n_active_cells());
-		active_set.clear();
-		active_set.set_size(locally_relevant_dofs.size());
-	}
-
-	// setup sparsity pattern
-	{
-		TimerOutput::Scope t(computing_timer, "Setup: matrix");
-		TrilinosWrappers::SparsityPattern sp(locally_owned_dofs,
-				mpi_communicator);
-
-		DoFTools::make_sparsity_pattern(dof_handler, sp,
-				constraints_dirichlet_hanging_nodes, false,
-				Utilities::MPI::this_mpi_process(mpi_communicator));
-
-		sp.compress();
-
-		system_matrix_newton.reinit(sp);
-
-		// we are going to reuse the system
-		// matrix for assembling the diagonal
-		// of the mass matrix so that we do not
-		// need to allocate two sparse matrices
-		// at the same time:
-		TrilinosWrappers::SparseMatrix & mass_matrix = system_matrix_newton;
-		assemble_mass_matrix_diagonal(mass_matrix);
-		const unsigned int start = (system_rhs_newton.local_range().first),
-				end = (system_rhs_newton.local_range().second);
-		for (unsigned int j = start; j < end; j++)
-		  diag_mass_matrix_vector(j) = mass_matrix.diag_element(j);
-
-		number_iterations = 0;
-
-		diag_mass_matrix_vector.compress(VectorOperation::insert);
-
-		// remove the mass matrix entries from the matrix:
-		mass_matrix = 0;
-	}
-}
-
-template
-void PlasticityContactProblem::assemble_nl_system(
-		TrilinosWrappers::MPI::Vector &u) {
-	TimerOutput::Scope t(computing_timer, "Assembling");
-
-        QGauss quadrature_formula(fe.degree + 1);
-        QGauss face_quadrature_formula(fe.degree + 1);
-
-	FEValues < dim
-			> fe_values(fe, quadrature_formula,
-					UpdateFlags(
-							update_values | update_gradients | update_q_points
-									| update_JxW_values));
-
-	FEFaceValues < dim
-			> fe_values_face(fe, face_quadrature_formula,
-					update_values | update_quadrature_points
-							| update_JxW_values);
-
-	const unsigned int dofs_per_cell = fe.dofs_per_cell;
-	const unsigned int n_q_points = quadrature_formula.size();
-	const unsigned int n_face_q_points = face_quadrature_formula.size();
-
-	const EquationData::RightHandSide right_hand_side;
-	std::vector < Vector
-			> right_hand_side_values(n_q_points, Vector(dim));
-	std::vector < Vector
-			> right_hand_side_values_face(n_face_q_points, Vector(dim));
-
-	FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell);
-	Vector cell_rhs(dofs_per_cell);
-
-	std::vector local_dof_indices(dofs_per_cell);
-
-	typename DoFHandler::active_cell_iterator cell =
-			dof_handler.begin_active(), endc = dof_handler.end();
-
-	const FEValuesExtractors::Vector displacement(0);
-
-	const double kappa = 1.0;
-	for (; cell != endc; ++cell)
-		if (cell->is_locally_owned()) {
-			fe_values.reinit(cell);
-			cell_matrix = 0;
-			cell_rhs = 0;
-
-			right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
-					right_hand_side_values);
-
-			std::vector < SymmetricTensor<2, dim> > strain_tensor(n_q_points);
-			fe_values[displacement].get_function_symmetric_gradients(u,
-					strain_tensor);
-
-			for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
-				SymmetricTensor < 4, dim > stress_strain_tensor_linearized;
-				SymmetricTensor < 4, dim > stress_strain_tensor;
-				SymmetricTensor < 2, dim > stress_tensor;
-
-				plast_lin_hard->linearized_plast_linear_hardening(
-						stress_strain_tensor_linearized, stress_strain_tensor,
-						strain_tensor[q_point]);
-
-				for (unsigned int i = 0; i < dofs_per_cell; ++i) {
-					stress_tensor = stress_strain_tensor_linearized
-							* plast_lin_hard->get_strain(fe_values, i, q_point);
-
-					for (unsigned int j = 0; j < dofs_per_cell; ++j) {
-						cell_matrix(i, j) += (stress_tensor
-								* plast_lin_hard->get_strain(fe_values, j,
-										q_point) * fe_values.JxW(q_point));
-					}
-
-					// the linearized part a(v^i;v^i,v) of the rhs
-					cell_rhs(i) += (stress_tensor * strain_tensor[q_point]
-							* fe_values.JxW(q_point));
-
-					// the residual part a(v^i;v) of the rhs
-					cell_rhs(i) -= (strain_tensor[q_point]
-							* stress_strain_tensor
-							* plast_lin_hard->get_strain(fe_values, i, q_point)
-							* fe_values.JxW(q_point));
-
-					// the residual part F(v) of the rhs
-					Tensor < 1, dim > rhs_values;
-					rhs_values = 0;
-					cell_rhs(i) += (fe_values[displacement].value(i, q_point)
-							* rhs_values * fe_values.JxW(q_point));
-				}
-			}
-
-			for (unsigned int face = 0;
-					face < GeometryInfo < dim > ::faces_per_cell; ++face) {
-				if (cell->face(face)->at_boundary()
-						&& cell->face(face)->boundary_indicator() == 1) {
-					fe_values_face.reinit(cell, face);
-
-					right_hand_side.vector_value_list(
-							fe_values_face.get_quadrature_points(),
-							right_hand_side_values_face);
-
-					for (unsigned int q_point = 0; q_point < n_face_q_points;
-							++q_point) {
-						Tensor < 1, dim > rhs_values;
-						rhs_values[2] = right_hand_side_values[q_point][2];
-						for (unsigned int i = 0; i < dofs_per_cell; ++i)
-							cell_rhs(i) += (fe_values_face[displacement].value(
-									i, q_point) * rhs_values
-									* fe_values_face.JxW(q_point));
-					}
-				}
-			}
-
-			cell->get_dof_indices(local_dof_indices);
-			constraints.distribute_local_to_global(cell_matrix, cell_rhs,
-					local_dof_indices, system_matrix_newton, system_rhs_newton,
-					true);
-
-
-		};
-
-	system_matrix_newton.compress(VectorOperation::add);
-	system_rhs_newton.compress(VectorOperation::add);
-}
-
-template
-void PlasticityContactProblem::residual_nl_system(
-		TrilinosWrappers::MPI::Vector &u) {
-        QGauss quadrature_formula(fe.degree + 1);
-        QGauss face_quadrature_formula(fe.degree + 1);
-
-	FEValues < dim
-			> fe_values(fe, quadrature_formula,
-					UpdateFlags(
-							update_values | update_gradients | update_q_points
-									| update_JxW_values));
-
-	FEFaceValues < dim
-			> fe_values_face(fe, face_quadrature_formula,
-					update_values | update_quadrature_points
-							| update_JxW_values);
-
-	const unsigned int dofs_per_cell = fe.dofs_per_cell;
-	const unsigned int n_q_points = quadrature_formula.size();
-	const unsigned int n_face_q_points = face_quadrature_formula.size();
-
-	const EquationData::RightHandSide right_hand_side;
-	std::vector < Vector
-			> right_hand_side_values(n_q_points, Vector(dim));
-	std::vector < Vector
-			> right_hand_side_values_face(n_face_q_points, Vector(dim));
-
-	Vector cell_rhs(dofs_per_cell);
-
-	std::vector local_dof_indices(dofs_per_cell);
-
-	const FEValuesExtractors::Vector displacement(0);
-
-	typename DoFHandler::active_cell_iterator cell =
-			dof_handler.begin_active(), endc = dof_handler.end();
-
-	unsigned int elast_points = 0;
-	unsigned int plast_points = 0;
-	double yield = 0;
-	unsigned int cell_number = 0;
-	cell_constitution = 0;
-
-	for (; cell != endc; ++cell)
-		if (cell->is_locally_owned()) {
-			fe_values.reinit(cell);
-			cell_rhs = 0;
-
-			right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
-					right_hand_side_values);
-
-			std::vector < SymmetricTensor<2, dim> > strain_tensor(n_q_points);
-			fe_values[displacement].get_function_symmetric_gradients(u,
-					strain_tensor);
-
-			for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
-				SymmetricTensor < 4, dim > stress_strain_tensor;
-				SymmetricTensor < 2, dim > stress_tensor;
-
-				plast_lin_hard->plast_linear_hardening(stress_strain_tensor,
-						strain_tensor[q_point], elast_points, plast_points,
-						yield);
-
-				cell_constitution(cell_number) += yield;
-				for (unsigned int i = 0; i < dofs_per_cell; ++i) {
-					cell_rhs(i) -= (strain_tensor[q_point]
-							* stress_strain_tensor
-							* //(stress_tensor) *
-							plast_lin_hard->get_strain(fe_values, i, q_point)
-							* fe_values.JxW(q_point));
-
-					Tensor < 1, dim > rhs_values;
-					rhs_values = 0;
-					cell_rhs(i) += ((fe_values[displacement].value(i, q_point)
-							* rhs_values) * fe_values.JxW(q_point));
-				};
-			};
-
-			for (unsigned int face = 0;
-					face < GeometryInfo < dim > ::faces_per_cell; ++face) {
-				if (cell->face(face)->at_boundary()
-						&& cell->face(face)->boundary_indicator() == 1) {
-					fe_values_face.reinit(cell, face);
-
-					right_hand_side.vector_value_list(
-							fe_values_face.get_quadrature_points(),
-							right_hand_side_values_face);
-
-					for (unsigned int q_point = 0; q_point < n_face_q_points;
-							++q_point) {
-						Tensor < 1, dim > rhs_values;
-						rhs_values[2] = right_hand_side_values[q_point][2];
-						for (unsigned int i = 0; i < dofs_per_cell; ++i)
-							cell_rhs(i) += (fe_values_face[displacement].value(
-									i, q_point) * rhs_values
-									* fe_values_face.JxW(q_point));
-					}
-				}
-			}
-
-			cell->get_dof_indices(local_dof_indices);
-			constraints_dirichlet_hanging_nodes.distribute_local_to_global(
-					cell_rhs, local_dof_indices, system_rhs_newton);
-
-			for (unsigned int i=0; i
-void PlasticityContactProblem::assemble_mass_matrix_diagonal(
-		TrilinosWrappers::SparseMatrix &mass_matrix) {
-        QGaussLobatto < dim - 1 > face_quadrature_formula(fe.degree + 1);
+  template 
+    void
+    PlasticityContactProblem::assemble_mass_matrix_diagonal (
+        TrilinosWrappers::SparseMatrix &mass_matrix)
+    {
+      QGaussLobatto face_quadrature_formula(fe.degree + 1);
 
-	FEFaceValues < dim
-			> fe_values_face(fe, face_quadrature_formula,
-					update_values | update_quadrature_points
-							| update_JxW_values);
+      FEFaceValues fe_values_face(fe, face_quadrature_formula,
+          update_values | update_quadrature_points | update_JxW_values);
 
-	const unsigned int dofs_per_cell = fe.dofs_per_cell;
-	const unsigned int n_face_q_points = face_quadrature_formula.size();
+      const unsigned int dofs_per_cell = fe.dofs_per_cell;
+      const unsigned int n_face_q_points = face_quadrature_formula.size();
 
-	FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell);
-	Tensor<1, dim, double> ones(dim);
-	for (unsigned i = 0; i < dim; i++)
-		ones[i] = 1.0;
+      FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell);
+      Tensor<1, dim, double> ones(dim);
+      for (unsigned i = 0; i < dim; i++)
+        ones[i] = 1.0;
 
-	std::vector local_dof_indices(dofs_per_cell);
+      std::vector local_dof_indices(dofs_per_cell);
 
-	const FEValuesExtractors::Vector displacement(0);
+      const FEValuesExtractors::Vector displacement(0);
 
-	typename DoFHandler::active_cell_iterator cell =
-			dof_handler.begin_active(), endc = dof_handler.end();
+      typename DoFHandler::active_cell_iterator cell =
+          dof_handler.begin_active(), endc = dof_handler.end();
 
-	for (; cell != endc; ++cell)
-		if (cell->is_locally_owned())
-			for (unsigned int face = 0;
-					face < GeometryInfo < dim > ::faces_per_cell; ++face)
-				if (cell->face(face)->at_boundary()
-						&& cell->face(face)->boundary_indicator() == 1) {
-					fe_values_face.reinit(cell, face);
-					cell_matrix = 0;
+      for (; cell != endc; ++cell)
+        if (cell->is_locally_owned())
+          for (unsigned int face = 0; face < GeometryInfo::faces_per_cell;
+              ++face)
+            if (cell->face(face)->at_boundary()
+                && cell->face(face)->boundary_indicator() == 1)
+              {
+                fe_values_face.reinit(cell, face);
+                cell_matrix = 0;
 
-					for (unsigned int q_point = 0; q_point < n_face_q_points;
-							++q_point)
-						for (unsigned int i = 0; i < dofs_per_cell; ++i)
-							cell_matrix(i, i) +=
-									(fe_values_face[displacement].value(i,
-											q_point) * ones
-											* fe_values_face.JxW(q_point));
+                for (unsigned int q_point = 0; q_point < n_face_q_points;
+                    ++q_point)
+                  for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                    cell_matrix(i, i) += (fe_values_face[displacement].value(i,
+                        q_point) * ones * fe_values_face.JxW(q_point));
 
-					cell->get_dof_indices(local_dof_indices);
+                cell->get_dof_indices(local_dof_indices);
 
 // 					constraints_dirichlet_hanging_nodes.distribute_local_to_global(
 // 							cell_matrix, local_dof_indices, mass_matrix);
 
-					for (unsigned int i=0; i
-void PlasticityContactProblem::update_solution_and_constraints() {
-	const EquationData::Obstacle obstacle(input_obstacle,
-						   (obstacle_filename!=""),
-						   (base_mesh=="box"?1.0:0.5));
-	std::vector vertex_touched(dof_handler.n_dofs(), false);
-
-	typename DoFHandler::active_cell_iterator cell =
-			dof_handler.begin_active(), endc = dof_handler.end();
-
-	TrilinosWrappers::MPI::Vector distributed_solution(system_rhs_newton);
-	distributed_solution = solution;
-	TrilinosWrappers::MPI::Vector lambda(solution);
-	lambda = resid_vector;
-	TrilinosWrappers::MPI::Vector diag_mass_matrix_vector_relevant(solution);
-	diag_mass_matrix_vector_relevant = diag_mass_matrix_vector;
-
-	constraints.reinit(locally_relevant_dofs);
-	active_set.clear();
-	IndexSet active_set_locally_owned;
-	active_set_locally_owned.set_size(locally_owned_dofs.size());
-	const double c = 100.0 * e_modul;
-
-	Quadrature          face_quadrature (fe.get_unit_face_support_points());
-	FEFaceValues          fe_values_face (fe, face_quadrature, update_quadrature_points);
-
-	const unsigned int         dofs_per_face = fe.dofs_per_face;
-	const unsigned int         n_face_q_points = face_quadrature.size ();
-
-	// pcout<< "dofs_per_face = " << dofs_per_face
-	//      << "n_face_q_points = " << n_face_q_points
-	//      <is_artificial())
-			for (unsigned int face = 0;
-					face < GeometryInfo < dim > ::faces_per_cell; ++face)
-				if (cell->face(face)->at_boundary()
-						&& cell->face(face)->boundary_indicator() == 1)
-				  {
-				        fe_values_face.reinit (cell, face);
-				        std::vector dof_indices (dofs_per_face);
-					cell->face(face)->get_dof_indices (dof_indices);
-
-					for (unsigned int q_point=0; q_point point(fe_values_face.quadrature_point(q_point));
-
-						double obstacle_value = obstacle.value(point, 2);
-						double solution_index_z = solution(index_z);
-						double gap = obstacle_value - point(2);
-
-						if (lambda(index_z)/diag_mass_matrix_vector_relevant(index_z)
-								+ c * (solution_index_z - gap) > 0
-								&& !(constraints_hanging_nodes.is_constrained(
-										index_z))) {
-							constraints.add_line(index_z);
-							constraints.set_inhomogeneity(index_z, gap);
-							distributed_solution(index_z) = gap;
-
-							if (locally_owned_dofs.is_element(index_z)) {
-								active_set_locally_owned.add_index(index_z);
-								if (locally_relevant_dofs.is_element(index_z))
-									active_set.add_index(index_z);
-							}
-
-						}
-						else if (lambda(index_z)/diag_mass_matrix_vector_relevant(index_z) 
-							 + c * (solution_index_z - gap) > 0
-							 && constraints_hanging_nodes.is_constrained(
-												     index_z))
-						  {
-						    if (locally_owned_dofs.is_element(index_z))
-						      {
-							counter_hanging_nodes += 1;
+  template 
+    void
+    PlasticityContactProblem::update_solution_and_constraints ()
+    {
+      const EquationData::Obstacle obstacle(input_obstacle,
+          (obstacle_filename != ""), (base_mesh == "box" ? 1.0 : 0.5));
+      std::vector vertex_touched(dof_handler.n_dofs(), false);
+
+      typename DoFHandler::active_cell_iterator cell =
+          dof_handler.begin_active(), endc = dof_handler.end();
+
+      TrilinosWrappers::MPI::Vector distributed_solution(system_rhs_newton);
+      distributed_solution = solution;
+      TrilinosWrappers::MPI::Vector lambda(solution);
+      lambda = resid_vector;
+      TrilinosWrappers::MPI::Vector diag_mass_matrix_vector_relevant(solution);
+      diag_mass_matrix_vector_relevant = diag_mass_matrix_vector;
+
+      constraints.reinit(locally_relevant_dofs);
+      active_set.clear();
+      IndexSet active_set_locally_owned;
+      active_set_locally_owned.set_size(locally_owned_dofs.size());
+      const double c = 100.0 * e_modul;
+
+      Quadrature face_quadrature(fe.get_unit_face_support_points());
+      FEFaceValues fe_values_face(fe, face_quadrature,
+          update_quadrature_points);
+
+      const unsigned int dofs_per_face = fe.dofs_per_face;
+      const unsigned int n_face_q_points = face_quadrature.size();
+
+      // pcout<< "dofs_per_face = " << dofs_per_face
+      //      << "n_face_q_points = " << n_face_q_points
+      //      <is_artificial())
+          for (unsigned int face = 0; face < GeometryInfo::faces_per_cell;
+              ++face)
+            if (cell->face(face)->at_boundary()
+                && cell->face(face)->boundary_indicator() == 1)
+              {
+                fe_values_face.reinit(cell, face);
+                std::vector dof_indices(dofs_per_face);
+                cell->face(face)->get_dof_indices(dof_indices);
+
+                for (unsigned int q_point = 0; q_point < n_face_q_points;
+                    ++q_point)
+                  {
+                    unsigned int component = fe.face_system_to_component_index(
+                        q_point).first;
+
+                    if (component == 2)
+                      {
+                        unsigned int index_z = dof_indices[q_point];
+
+                        if (vertex_touched[index_z] == false)
+                          vertex_touched[index_z] = true;
+                        else
+                          continue;
+
+                        // the local row where
+                        Point point(
+                            fe_values_face.quadrature_point(q_point));
+
+                        double obstacle_value = obstacle.value(point, 2);
+                        double solution_index_z = solution(index_z);
+                        double gap = obstacle_value - point(2);
+
+                        if (lambda(index_z)
+                            / diag_mass_matrix_vector_relevant(index_z)
+                            + c * (solution_index_z - gap) > 0
+                            && !(constraints_hanging_nodes.is_constrained(
+                                index_z)))
+                          {
+                            constraints.add_line(index_z);
+                            constraints.set_inhomogeneity(index_z, gap);
+                            distributed_solution(index_z) = gap;
+
+                            if (locally_owned_dofs.is_element(index_z))
+                              {
+                                active_set_locally_owned.add_index(index_z);
+                                if (locally_relevant_dofs.is_element(index_z))
+                                  active_set.add_index(index_z);
+                              }
+
+                          }
+                        else if (lambda(index_z)
+                            / diag_mass_matrix_vector_relevant(index_z)
+                            + c * (solution_index_z - gap) > 0
+                            && constraints_hanging_nodes.is_constrained(
+                                index_z))
+                          {
+                            if (locally_owned_dofs.is_element(index_z))
+                              {
+                                counter_hanging_nodes += 1;
 
 // 							std::cout << "index_z = " << index_z
 // 								<< ", lambda = " << lambda (index_z)
@@ -1340,32 +1429,32 @@ void PlasticityContactProblem::update_solution_and_constraints() {
 //  								<< ", x = " << point(0)
 //  								<< ", y = " << point(1)
 //  								<< std::endl;
-						      }
-						  }
-					    }
-					}
-				  }
-	distributed_solution.compress(VectorOperation::insert);
+                              }
+                          }
+                      }
+                  }
+              }
+      distributed_solution.compress(VectorOperation::insert);
 
-	unsigned int sum_contact_constraints = Utilities::MPI::sum(
-			active_set_locally_owned.n_elements(), mpi_communicator);
-	pcout << "         Size of active set: " << sum_contact_constraints
-			<< std::endl;
-	unsigned int sum_contact_hanging_nodes = Utilities::MPI::sum(
-			counter_hanging_nodes, mpi_communicator);
-	pcout << "         Number of hanging nodes in contact: " << sum_contact_hanging_nodes
-			<< std::endl;
+      unsigned int sum_contact_constraints = Utilities::MPI::sum(
+          active_set_locally_owned.n_elements(), mpi_communicator);
+      pcout << "         Size of active set: " << sum_contact_constraints
+          << std::endl;
+      unsigned int sum_contact_hanging_nodes = Utilities::MPI::sum(
+          counter_hanging_nodes, mpi_communicator);
+      pcout << "         Number of hanging nodes in contact: "
+          << sum_contact_hanging_nodes << std::endl;
 
-	solution = distributed_solution;
+      solution = distributed_solution;
 
-	constraints.close();
+      constraints.close();
 
- //  	constraints_dirichlet_hanging_nodes.print (std::cout);
+      //  	constraints_dirichlet_hanging_nodes.print (std::cout);
 
-	constraints.merge(constraints_dirichlet_hanging_nodes);
+      constraints.merge(constraints_dirichlet_hanging_nodes);
 
- 	//constraints.print (std::cout);
-}
+      //constraints.print (std::cout);
+    }
 
 // @sect4{PlasticityContactProblem::dirichlet_constraints}
 
@@ -1373,38 +1462,39 @@ void PlasticityContactProblem::update_solution_and_constraints() {
 // constraints_dirichlet_hanging_nodes. It contains
 // the dirichlet boundary values as well as the
 // hanging nodes constraints.
-template
-void PlasticityContactProblem::dirichlet_constraints() {
-	/* boundary_indicators:
-	 _______
-	 /  1    /|
-	 /______ / |
-	 8|       | 8|
-	 |   8   | /
-	 |_______|/
-	 6
-	 */
-
-	constraints_dirichlet_hanging_nodes.reinit(locally_relevant_dofs);
-	constraints_dirichlet_hanging_nodes.merge(constraints_hanging_nodes);
-
-	// interpolate all components of the solution
-	VectorTools::interpolate_boundary_values(dof_handler, base_mesh=="box"?6:0,
-			EquationData::BoundaryValues(),
-			constraints_dirichlet_hanging_nodes, ComponentMask());
-
-	// interpolate x- and y-components of the
-	// solution (this is a bit mask, so apply
-	// operator| )
-	FEValuesExtractors::Scalar x_displacement(0);
-	FEValuesExtractors::Scalar y_displacement(1);
-	VectorTools::interpolate_boundary_values(dof_handler, 8,
-			EquationData::BoundaryValues(),
-			constraints_dirichlet_hanging_nodes,
-			(fe.component_mask(x_displacement)
-					| fe.component_mask(y_displacement)));
-	constraints_dirichlet_hanging_nodes.close();
-}
+  template 
+    void
+    PlasticityContactProblem::dirichlet_constraints ()
+    {
+      /* boundary_indicators:
+       _______
+       /  1    /|
+       /______ / |
+       8|       | 8|
+       |   8   | /
+       |_______|/
+       6
+       */
+
+      constraints_dirichlet_hanging_nodes.reinit(locally_relevant_dofs);
+      constraints_dirichlet_hanging_nodes.merge(constraints_hanging_nodes);
+
+      // interpolate all components of the solution
+      VectorTools::interpolate_boundary_values(dof_handler,
+          base_mesh == "box" ? 6 : 0, EquationData::BoundaryValues(),
+          constraints_dirichlet_hanging_nodes, ComponentMask());
+
+      // interpolate x- and y-components of the
+      // solution (this is a bit mask, so apply
+      // operator| )
+      FEValuesExtractors::Scalar x_displacement(0);
+      FEValuesExtractors::Scalar y_displacement(1);
+      VectorTools::interpolate_boundary_values(dof_handler, 8,
+          EquationData::BoundaryValues(),
+          constraints_dirichlet_hanging_nodes,
+          (fe.component_mask(x_displacement) | fe.component_mask(y_displacement)));
+      constraints_dirichlet_hanging_nodes.close();
+    }
 
 // @sect4{PlasticityContactProblem::solve}
 
@@ -1432,59 +1522,61 @@ void PlasticityContactProblem::dirichlet_constraints() {
 // instead of CG. For a very small hardening
 // value gamma the linear system becomes
 // almost semi definite but still symmetric.
-template
-void PlasticityContactProblem::solve() {
-	TimerOutput::Scope t(computing_timer, "Solve");
-
-	TrilinosWrappers::MPI::Vector distributed_solution(system_rhs_newton);
-	distributed_solution = solution;
-
-	constraints_hanging_nodes.set_zero(distributed_solution);
-	constraints_hanging_nodes.set_zero(system_rhs_newton);
-	distributed_solution.compress(VectorOperation::insert);
-	system_rhs_newton.compress(VectorOperation::insert);
-
-	{
-		TimerOutput::Scope t(computing_timer, "Solve: setup preconditioner");
-		preconditioner_u.initialize(system_matrix_newton, additional_data);
-	}
-
-	{
-		TimerOutput::Scope t(computing_timer, "Solve: iterate");
-
-		PrimitiveVectorMemory < TrilinosWrappers::MPI::Vector > mem;
-		TrilinosWrappers::MPI::Vector tmp(system_rhs_newton);
-		// 1e-4 seems to be the fasted option altogether, but to get more
-		// reproducible parallel benchmark results, we use a small residual:
-		double relative_accuracy = 1e-8;
-		if (output_dir.compare("its/") == 0)
-		  relative_accuracy = 1e-4;
-
-		const double solver_tolerance = relative_accuracy
-		  * system_matrix_newton.residual(tmp, distributed_solution,
-						  system_rhs_newton);
-
-		SolverControl solver_control(system_matrix_newton.m(),
-				solver_tolerance);
-		SolverBicgstab < TrilinosWrappers::MPI::Vector
-				> solver(solver_control, mem/*,
-				 SolverFGMRES::
-				 AdditionalData(30, true)*/);
-		solver.solve(system_matrix_newton, distributed_solution,
-				system_rhs_newton, preconditioner_u);
-
-		pcout << "         Error: " << solver_control.initial_value() << " -> "
-				<< solver_control.last_value() << " in "
-				<< solver_control.last_step() << " Bicgstab iterations."
-				<< std::endl;
-
-		number_iterations += solver_control.last_step();
-	}
-
-	constraints.distribute(distributed_solution);
-
-	solution = distributed_solution;
-}
+  template 
+    void
+    PlasticityContactProblem::solve ()
+    {
+      TimerOutput::Scope t(computing_timer, "Solve");
+
+      TrilinosWrappers::MPI::Vector distributed_solution(system_rhs_newton);
+      distributed_solution = solution;
+
+      constraints_hanging_nodes.set_zero(distributed_solution);
+      constraints_hanging_nodes.set_zero(system_rhs_newton);
+      distributed_solution.compress(VectorOperation::insert);
+      system_rhs_newton.compress(VectorOperation::insert);
+
+        {
+          TimerOutput::Scope t(computing_timer, "Solve: setup preconditioner");
+          preconditioner_u.initialize(system_matrix_newton, additional_data);
+        }
+
+        {
+          TimerOutput::Scope t(computing_timer, "Solve: iterate");
+
+          PrimitiveVectorMemory mem;
+          TrilinosWrappers::MPI::Vector tmp(system_rhs_newton);
+          // 1e-4 seems to be the fasted option altogether, but to get more
+          // reproducible parallel benchmark results, we use a small residual:
+          double relative_accuracy = 1e-8;
+          if (output_dir.compare("its/") == 0)
+            relative_accuracy = 1e-4;
+
+          const double solver_tolerance = relative_accuracy
+              * system_matrix_newton.residual(tmp, distributed_solution,
+                  system_rhs_newton);
+
+          SolverControl solver_control(system_matrix_newton.m(),
+              solver_tolerance);
+          SolverBicgstab solver(solver_control,
+              mem/*,
+               SolverFGMRES::
+               AdditionalData(30, true)*/);
+          solver.solve(system_matrix_newton, distributed_solution,
+              system_rhs_newton, preconditioner_u);
+
+          pcout << "         Error: " << solver_control.initial_value()
+              << " -> " << solver_control.last_value() << " in "
+              << solver_control.last_step() << " Bicgstab iterations."
+              << std::endl;
+
+          number_iterations += solver_control.last_step();
+        }
+
+      constraints.distribute(distributed_solution);
+
+      solution = distributed_solution;
+    }
 
 // @sect4{PlasticityContactProblem::solve_newton}
 
@@ -1493,281 +1585,296 @@ void PlasticityContactProblem::solve() {
 // iteration and the inner loop for the damping steps which
 // will be used only if necessary. To obtain a good and reasonable
 // starting value we solve an elastic problem in very first step (j=1).
-template
-void PlasticityContactProblem::solve_newton() {
-	TimerOutput::Scope t(computing_timer, "solve newton setup");
-
-	double resid = 0;
-	double resid_old = 100000;
-	TrilinosWrappers::MPI::Vector old_solution(system_rhs_newton);
-	TrilinosWrappers::MPI::Vector res(system_rhs_newton);
-	TrilinosWrappers::MPI::Vector tmp_vector(system_rhs_newton);
-
-	std::vector < std::vector > constant_modes;
-	DoFTools::extract_constant_modes(dof_handler, ComponentMask(),
-			constant_modes);
-
-	double sigma_hlp = sigma_0;
-
-	additional_data.constant_modes = constant_modes;
-	additional_data.elliptic = true;
-	additional_data.n_cycles = 1;
-	additional_data.w_cycle = false;
-	additional_data.output_details = false;
-	additional_data.smoother_sweeps = 2;
-	additional_data.aggregation_threshold = 1e-2;
-
-	IndexSet active_set_old(active_set);
-
-	t.stop(); // stop newton setup timer
-
-	unsigned int j = 1;
-	unsigned int number_assemble_system = 0;
-	for (; j <= 100; j++) {
-	  if (transfer_solution)
-	    {
-	      if (transfer_solution && j == 1 && cycle == 0)
-		plast_lin_hard->set_sigma_0(1e+10);
-	      else if (transfer_solution && (j == 2 || cycle > 0))
-		plast_lin_hard->set_sigma_0(sigma_hlp);
-	    }
-	  else
-	    {
-	      if (j == 1)
-		plast_lin_hard->set_sigma_0(1e+10);
-	      else
-		plast_lin_hard->set_sigma_0(sigma_hlp);
-	    }
-
-		pcout << " " << std::endl;
-		pcout << "   Newton iteration " << j << std::endl;
-		pcout << "      Updating active set..." << std::endl;
-
-		{
-			TimerOutput::Scope t(computing_timer, "update active set");
-			update_solution_and_constraints();
-		}
-
-		pcout << "      Assembling system... " << std::endl;
-		system_matrix_newton = 0;
-		system_rhs_newton = 0;
-		assemble_nl_system(solution); //compute Newton-Matrix
-
-		number_assemble_system += 1;
-
-		pcout << "      Solving system... " << std::endl;
-		solve();
-
-		TrilinosWrappers::MPI::Vector distributed_solution(system_rhs_newton);
-		distributed_solution = solution;
-
-		// We handle a highly nonlinear problem so we have to damp
-		// the Newtons method. We refer that we iterate the new solution
-		// in each Newton step and not only the solution update.
-		// Since the solution set is a convex set and not a space we
-		// compute for the damping a linear combination of the
-		// previous and the current solution to guarantee that the
-		// damped solution is in our solution set again.
-		// At most we apply 10 damping steps.
-		bool damped = false;
-		tmp_vector = old_solution;
-		double a = 0;
-		for (unsigned int i = 0; (i < 5) && (!damped); i++) {
-			a = std::pow(0.5, static_cast(i));
-			old_solution = tmp_vector;
-			old_solution.sadd(1 - a, a, distributed_solution);
-			old_solution.compress(VectorOperation::add);
-
-			TimerOutput::Scope t(computing_timer, "Residual and lambda");
-
-			system_rhs_newton = 0;
-			system_rhs_lambda = 0;
-
-			solution = old_solution;
-			residual_nl_system(solution);
-			res = system_rhs_newton;
-
-			const unsigned int start_res = (res.local_range().first), end_res =
-					(res.local_range().second);
-			for (unsigned int n = start_res; n < end_res; ++n)
-				if (constraints.is_inhomogeneously_constrained(n))
-					res(n) = 0;
-
-			res.compress(VectorOperation::insert);
-
-			resid = res.l2_norm();
-
-			if (resid < resid_old)
-				damped = true;
-
-			pcout << "      Residual of the non-contact part of the system: "
-					<< resid << std::endl
-					<< "         with a damping parameter alpha = " << a
-					<< std::endl;
-
-			// The previous iteration of step 0 is the solution of an elastic problem.
-			// So a linear combination of a plastic and an elastic solution makes no sense
-			// since the elastic solution is not in the convex set of the plastic solution.
-			if (!transfer_solution && j == 2)
-				break;
-			if (transfer_solution && j == 2 && cycle == 0)
-				break;
-		}
-
-		resid_old = resid;
-
-		resid_vector = system_rhs_lambda;
-		resid_vector.compress(VectorOperation::insert);
-
-		int is_my_set_changed = (active_set == active_set_old) ? 0 : 1;
-		int num_changed = Utilities::MPI::sum(is_my_set_changed,
-				MPI_COMM_WORLD);
-		if (num_changed == 0)
-		  {
-		    pcout<< "      Active set did not change!" < > constant_modes;
+      DoFTools::extract_constant_modes(dof_handler, ComponentMask(),
+          constant_modes);
+
+      double sigma_hlp = sigma_0;
+
+      additional_data.constant_modes = constant_modes;
+      additional_data.elliptic = true;
+      additional_data.n_cycles = 1;
+      additional_data.w_cycle = false;
+      additional_data.output_details = false;
+      additional_data.smoother_sweeps = 2;
+      additional_data.aggregation_threshold = 1e-2;
+
+      IndexSet active_set_old(active_set);
+
+      t.stop(); // stop newton setup timer
+
+      unsigned int j = 1;
+      unsigned int number_assemble_system = 0;
+      for (; j <= 100; j++)
+        {
+          if (transfer_solution)
+            {
+              if (transfer_solution && j == 1 && cycle == 0)
+                plast_lin_hard->set_sigma_0(1e+10);
+              else if (transfer_solution && (j == 2 || cycle > 0))
+                plast_lin_hard->set_sigma_0(sigma_hlp);
+            }
+          else
+            {
+              if (j == 1)
+                plast_lin_hard->set_sigma_0(1e+10);
+              else
+                plast_lin_hard->set_sigma_0(sigma_hlp);
+            }
+
+          pcout << " " << std::endl;
+          pcout << "   Newton iteration " << j << std::endl;
+          pcout << "      Updating active set..." << std::endl;
+
+            {
+              TimerOutput::Scope t(computing_timer, "update active set");
+              update_solution_and_constraints();
+            }
+
+          pcout << "      Assembling system... " << std::endl;
+          system_matrix_newton = 0;
+          system_rhs_newton = 0;
+          assemble_nl_system(solution); //compute Newton-Matrix
+
+          number_assemble_system += 1;
+
+          pcout << "      Solving system... " << std::endl;
+          solve();
+
+          TrilinosWrappers::MPI::Vector distributed_solution(system_rhs_newton);
+          distributed_solution = solution;
+
+          // We handle a highly nonlinear problem so we have to damp
+          // the Newtons method. We refer that we iterate the new solution
+          // in each Newton step and not only the solution update.
+          // Since the solution set is a convex set and not a space we
+          // compute for the damping a linear combination of the
+          // previous and the current solution to guarantee that the
+          // damped solution is in our solution set again.
+          // At most we apply 10 damping steps.
+          bool damped = false;
+          tmp_vector = old_solution;
+          double a = 0;
+          for (unsigned int i = 0; (i < 5) && (!damped); i++)
+            {
+              a = std::pow(0.5, static_cast(i));
+              old_solution = tmp_vector;
+              old_solution.sadd(1 - a, a, distributed_solution);
+              old_solution.compress(VectorOperation::add);
+
+              TimerOutput::Scope t(computing_timer, "Residual and lambda");
+
+              system_rhs_newton = 0;
+              system_rhs_lambda = 0;
+
+              solution = old_solution;
+              residual_nl_system(solution);
+              res = system_rhs_newton;
+
+              const unsigned int start_res = (res.local_range().first),
+                  end_res = (res.local_range().second);
+              for (unsigned int n = start_res; n < end_res; ++n)
+                if (constraints.is_inhomogeneously_constrained(n))
+                  res(n) = 0;
+
+              res.compress(VectorOperation::insert);
+
+              resid = res.l2_norm();
+
+              if (resid < resid_old)
+                damped = true;
+
+              pcout << "      Residual of the non-contact part of the system: "
+                  << resid << std::endl
+                  << "         with a damping parameter alpha = " << a
+                  << std::endl;
+
+              // The previous iteration of step 0 is the solution of an elastic problem.
+              // So a linear combination of a plastic and an elastic solution makes no sense
+              // since the elastic solution is not in the convex set of the plastic solution.
+              if (!transfer_solution && j == 2)
+                break;
+              if (transfer_solution && j == 2 && cycle == 0)
+                break;
+            }
+
+          resid_old = resid;
+
+          resid_vector = system_rhs_lambda;
+          resid_vector.compress(VectorOperation::insert);
+
+          int is_my_set_changed = (active_set == active_set_old) ? 0 : 1;
+          int num_changed = Utilities::MPI::sum(is_my_set_changed,
+              MPI_COMM_WORLD);
+          if (num_changed == 0)
+            {
+              pcout << "      Active set did not change!" << std::endl;
+              if (output_dir.compare("its/") != 0 && resid < 1e-7)
+                break;
+              else if (output_dir.compare("its/") == 0 && resid < 1e-10)
+                break;
+            }
+          active_set_old = active_set;
+        }
+
+      pcout << "" << std::endl << "      Number of assembled systems = "
+          << number_assemble_system << std::endl
+          << "      Number of Solver-Iterations = " << number_iterations
+          << std::endl;
     }
-  else
-    {
-	Vector estimated_error_per_cell(triangulation.n_active_cells());
-	KellyErrorEstimator < dim
-			      > ::estimate(dof_handler, QGauss < dim - 1 > (fe.degree + 2),
-					typename FunctionMap::type(), solution,
-					estimated_error_per_cell);
 
-	parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number(
-			triangulation, estimated_error_per_cell, 0.3, 0.03);
-
-	triangulation.prepare_coarsening_and_refinement();
-	if (transfer_solution)
-	  soltrans->prepare_for_coarsening_and_refinement(solution);
+// @sect3{The refine_grid function}
 
-	triangulation.execute_coarsening_and_refinement();
+  template 
+    void
+    PlasticityContactProblem::refine_grid ()
+    {
+      if (refinement_strategy == RefinementStrategy::refine_global)
+        {
+          triangulation.refine_global(1);
+        }
+      else
+        {
+          Vector estimated_error_per_cell(
+              triangulation.n_active_cells());
+          KellyErrorEstimator::estimate(dof_handler,
+              QGauss(fe.degree + 2), typename FunctionMap::type(),
+              solution, estimated_error_per_cell);
+
+          parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number(
+              triangulation, estimated_error_per_cell, 0.3, 0.03);
+
+          triangulation.prepare_coarsening_and_refinement();
+          if (transfer_solution)
+            soltrans->prepare_for_coarsening_and_refinement(solution);
+
+          triangulation.execute_coarsening_and_refinement();
+        }
     }
-}
 
 // @sect3{The move_mesh function}
 
-template
-void PlasticityContactProblem::move_mesh(
-		const TrilinosWrappers::MPI::Vector &_complete_displacement) const {
-	std::vector vertex_touched(triangulation.n_vertices(), false);
-
-	for (typename DoFHandler::active_cell_iterator cell =
-			dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
-		if (cell->is_locally_owned())
-			for (unsigned int v = 0;
-					v < GeometryInfo < dim > ::vertices_per_cell; ++v) {
-				if (vertex_touched[cell->vertex_index(v)] == false) {
-					vertex_touched[cell->vertex_index(v)] = true;
-
-					Point < dim > vertex_displacement;
-					for (unsigned int d = 0; d < dim; ++d) {
-						if (_complete_displacement(cell->vertex_dof_index(v, d))
-								!= 0)
-							vertex_displacement[d] = _complete_displacement(
-									cell->vertex_dof_index(v, d));
-					}
-
-					cell->vertex(v) += vertex_displacement;
-				}
-			}
-}
+  template 
+    void
+    PlasticityContactProblem::move_mesh (
+        const TrilinosWrappers::MPI::Vector &_complete_displacement) const
+    {
+      std::vector vertex_touched(triangulation.n_vertices(), false);
+
+      for (typename DoFHandler::active_cell_iterator cell =
+          dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
+        if (cell->is_locally_owned())
+          for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell;
+              ++v)
+            {
+              if (vertex_touched[cell->vertex_index(v)] == false)
+                {
+                  vertex_touched[cell->vertex_index(v)] = true;
+
+                  Point vertex_displacement;
+                  for (unsigned int d = 0; d < dim; ++d)
+                    {
+                      if (_complete_displacement(cell->vertex_dof_index(v, d))
+                          != 0)
+                        vertex_displacement[d] = _complete_displacement(
+                            cell->vertex_dof_index(v, d));
+                    }
+
+                  cell->vertex(v) += vertex_displacement;
+                }
+            }
+    }
 
 // @sect4{PlasticityContactProblem::output_results}
 
-template
-void PlasticityContactProblem::output_results(
-		const std::string &title) {
-	move_mesh(solution);
-
-	// Calculation of the contact forces
-	TrilinosWrappers::MPI::Vector lambda(solution);
-	TrilinosWrappers::MPI::Vector distributed_lambda(system_rhs_newton);
-	const unsigned int start_res = (resid_vector.local_range().first), end_res =
-	  (resid_vector.local_range().second);
-	for (unsigned int n = start_res; n < end_res; ++n)
-	  if (constraints.is_inhomogeneously_constrained(n))
-	    distributed_lambda(n) = resid_vector(n)/diag_mass_matrix_vector(n);
-	distributed_lambda.compress(VectorOperation::insert);
-	constraints_hanging_nodes.distribute(distributed_lambda);
-	lambda = distributed_lambda;
-	TrilinosWrappers::MPI::Vector resid_vector_relevant(solution);
-	TrilinosWrappers::MPI::Vector distributed_resid_vector(resid_vector);
-	constraints_hanging_nodes.distribute(distributed_resid_vector);
-	resid_vector_relevant = distributed_resid_vector;
-
-	DataOut < dim > data_out;
-
-	data_out.attach_dof_handler(dof_handler);
-
-	const std::vector data_component_interpretation(
-			dim, DataComponentInterpretation::component_is_part_of_vector);
-	data_out.add_data_vector(solution,
-			std::vector < std::string > (dim, "Displacement"),
-			DataOut < dim > ::type_dof_data, data_component_interpretation);
-	data_out.add_data_vector(lambda,
-			std::vector < std::string > (dim, "ContactForce"),
-			DataOut < dim > ::type_dof_data, data_component_interpretation);
-	data_out.add_data_vector(active_set,
-			std::vector < std::string > (dim, "ActiveSet"),
-			DataOut < dim > ::type_dof_data, data_component_interpretation);
-	data_out.add_data_vector(resid_vector_relevant,
-			std::vector < std::string > (dim, "Residual"),
-			DataOut < dim > ::type_dof_data, data_component_interpretation);
-
-	Vector subdomain(triangulation.n_active_cells());
-	for (unsigned int i = 0; i < subdomain.size(); ++i)
-		subdomain(i) = triangulation.locally_owned_subdomain();
-	data_out.add_data_vector(subdomain, "subdomain");
-
-	data_out.add_data_vector(cell_constitution, "CellConstitution");
-
-	data_out.build_patches();
-
-	const std::string filename = (output_dir + title + "-"
-			+ Utilities::int_to_string(triangulation.locally_owned_subdomain(),
-					4));
-
-	std::ofstream output_vtu((filename + ".vtu").c_str());
-	data_out.write_vtu(output_vtu);
-        pcout << output_dir + title << ".pvtu" << std::endl;
-
-	if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) {
-		std::vector filenames;
-		for (unsigned int i = 0;
-				i < Utilities::MPI::n_mpi_processes(mpi_communicator); ++i)
-			filenames.push_back(
-			    title + "-" + Utilities::int_to_string(i, 4) + ".vtu");
-
-		std::ofstream master_output((output_dir + title + ".pvtu").c_str());
-		data_out.write_pvtu_record(master_output, filenames);
-	}
-
-	TrilinosWrappers::MPI::Vector tmp(solution);
-	tmp *= -1;
-	move_mesh(tmp);
-}
+  template 
+    void
+    PlasticityContactProblem::output_results (
+        const std::string &title)
+    {
+      move_mesh(solution);
+
+      // Calculation of the contact forces
+      TrilinosWrappers::MPI::Vector lambda(solution);
+      TrilinosWrappers::MPI::Vector distributed_lambda(system_rhs_newton);
+      const unsigned int start_res = (resid_vector.local_range().first),
+          end_res = (resid_vector.local_range().second);
+      for (unsigned int n = start_res; n < end_res; ++n)
+        if (constraints.is_inhomogeneously_constrained(n))
+          distributed_lambda(n) = resid_vector(n) / diag_mass_matrix_vector(n);
+      distributed_lambda.compress(VectorOperation::insert);
+      constraints_hanging_nodes.distribute(distributed_lambda);
+      lambda = distributed_lambda;
+      TrilinosWrappers::MPI::Vector resid_vector_relevant(solution);
+      TrilinosWrappers::MPI::Vector distributed_resid_vector(resid_vector);
+      constraints_hanging_nodes.distribute(distributed_resid_vector);
+      resid_vector_relevant = distributed_resid_vector;
+
+      DataOut data_out;
+
+      data_out.attach_dof_handler(dof_handler);
+
+      const std::vector data_component_interpretation(
+          dim, DataComponentInterpretation::component_is_part_of_vector);
+      data_out.add_data_vector(solution,
+          std::vector < std::string > (dim, "Displacement"),
+          DataOut::type_dof_data, data_component_interpretation);
+      data_out.add_data_vector(lambda,
+          std::vector < std::string > (dim, "ContactForce"),
+          DataOut::type_dof_data, data_component_interpretation);
+      data_out.add_data_vector(active_set,
+          std::vector < std::string > (dim, "ActiveSet"),
+          DataOut::type_dof_data, data_component_interpretation);
+      data_out.add_data_vector(resid_vector_relevant,
+          std::vector < std::string > (dim, "Residual"),
+          DataOut::type_dof_data, data_component_interpretation);
+
+      Vector subdomain(triangulation.n_active_cells());
+      for (unsigned int i = 0; i < subdomain.size(); ++i)
+        subdomain(i) = triangulation.locally_owned_subdomain();
+      data_out.add_data_vector(subdomain, "subdomain");
+
+      data_out.add_data_vector(cell_constitution, "CellConstitution");
+
+      data_out.build_patches();
+
+      const std::string filename =
+          (output_dir + title + "-"
+              + Utilities::int_to_string(
+                  triangulation.locally_owned_subdomain(), 4));
+
+      std::ofstream output_vtu((filename + ".vtu").c_str());
+      data_out.write_vtu(output_vtu);
+      pcout << output_dir + title << ".pvtu" << std::endl;
+
+      if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
+        {
+          std::vector < std::string > filenames;
+          for (unsigned int i = 0;
+              i < Utilities::MPI::n_mpi_processes(mpi_communicator); ++i)
+            filenames.push_back(
+                title + "-" + Utilities::int_to_string(i, 4) + ".vtu");
+
+          std::ofstream master_output((output_dir + title + ".pvtu").c_str());
+          data_out.write_pvtu_record(master_output, filenames);
+        }
+
+      TrilinosWrappers::MPI::Vector tmp(solution);
+      tmp *= -1;
+      move_mesh(tmp);
+    }
 
 // @sect4{PlasticityContactProblem::output_contact_force}
 
@@ -1782,309 +1889,217 @@ void PlasticityContactProblem::output_results(
 // it is important to apply contraints_hanging_nodes.distribute
 // to the distributed_lambda vector.
 // To calculate the contact pressure in a certain point in the
-// contact area, we have make use of the Functions::FEFieldFunction
-// In parallel this is little tricky because we have to find the
-// process with the right cell which contains this point.
-template
-void PlasticityContactProblem::output_contact_force(
-		const unsigned int cycle) {
-	Functions::FEFieldFunction,
-			TrilinosWrappers::MPI::Vector> solution_function(dof_handler,
-			solution);
-        std::cout.precision(10);
-
-        Vector solution_p1(dim);
-	std::vector < Tensor<1, dim> > solution_gradient_p1(dim);
-
-	const Point p1_of_interest(0.5001, 0.5001, 0.9501);
-	bool point1_found = true;
-	bool point2_found = true;
-
-	// Calculation of the contact forces
-	TrilinosWrappers::MPI::Vector lambda(solution);
-	TrilinosWrappers::MPI::Vector distributed_lambda(system_rhs_newton);
-	const unsigned int start_res = (resid_vector.local_range().first), end_res =
-	  (resid_vector.local_range().second);
-	for (unsigned int n = start_res; n < end_res; ++n)
-	  if (constraints.is_inhomogeneously_constrained(n))
-	    distributed_lambda(n) = resid_vector(n)/diag_mass_matrix_vector(n);
-	  else
-	    distributed_lambda(n) = 0;
-	distributed_lambda.compress(VectorOperation::insert);
-	constraints_hanging_nodes.distribute(distributed_lambda);
-	lambda = distributed_lambda;
-	Functions::FEFieldFunction,
-		        TrilinosWrappers::MPI::Vector> lambda_function(dof_handler,
-		        lambda);
-	const Point p2_of_interest(0.49, 0.5001, 1.0);
-	Vector lambda_p2(dim);
-	
-	MPI_Barrier(MPI_COMM_WORLD);
-	try {
-	  lambda_function.vector_value(p2_of_interest, lambda_p2);
-	} catch (const typename Functions::FEFieldFunction,
-		 TrilinosWrappers::MPI::Vector>::ExcPointNotAvailableHere &) {
-	  point2_found = false;
-	}
-	
-	if (point2_found == true) {
-	  std::cout << "PoI lambda_z: " << lambda_p2(2) << std::endl;
-	}
-
-	// Integral of the contact force in z-direction over the whole contact area.
-	double contact_force = 0.0;
-	{
-		QGauss< dim - 1 > face_quadrature_formula(fe.degree + 1);
-
-		FEFaceValues < dim > fe_values_face(fe, face_quadrature_formula,
-						update_values | update_quadrature_points
-								| update_JxW_values);
-
-		const unsigned int n_face_q_points = face_quadrature_formula.size();
-
-		const FEValuesExtractors::Vector displacement(0);
-
-		typename DoFHandler::active_cell_iterator cell =
-				dof_handler.begin_active(), endc = dof_handler.end();
-		for (; cell != endc; ++cell)
-		  if (cell->is_locally_owned())
-				for (unsigned int face = 0;
-						face < GeometryInfo < dim > ::faces_per_cell; ++face)
-					if (cell->face(face)->at_boundary()
-							&& cell->face(face)->boundary_indicator() == 1) {
-						fe_values_face.reinit(cell, face);
-
-						std::vector < Tensor<1, dim>
-								> lambda_values(n_face_q_points);
-						fe_values_face[displacement].get_function_values(lambda,
-								lambda_values);
-
-						for (unsigned int q_point = 0;
-								q_point < n_face_q_points; ++q_point)
-						{
-							contact_force += lambda_values[q_point][2]
-									* fe_values_face.JxW(q_point);
-						}
-					}
-		contact_force = Utilities::MPI::sum(contact_force,
-				MPI_COMM_WORLD);
-		pcout << "Contact force = " << contact_force << std::endl;
-	}
-
-	// To calculate the contact area between deformable body and obstacle
-	double contact_area = 0.0;
-	{
-	  move_mesh(solution);
-
-	  QGaussLobatto < dim - 1 > face_quadrature_formula(fe.degree + 1);
-	  
-	  FEFaceValues < dim > fe_values_face(fe, face_quadrature_formula,
-					      update_values | update_quadrature_points
-					      | update_JxW_values);
-	  
-	  const unsigned int dofs_per_face = fe.dofs_per_face;
-	  const unsigned int n_face_q_points = face_quadrature_formula.size();
-	  
-	  const FEValuesExtractors::Vector displacement(0);
-	  typename DoFHandler::active_cell_iterator cell =
-	    dof_handler.begin_active(), endc = dof_handler.end();
-	  for (; cell != endc; ++cell)
-	    if (cell->is_locally_owned())
-	      for (unsigned int face = 0;
-		   face < GeometryInfo < dim > ::faces_per_cell; ++face){
-		if (cell->face(face)->at_boundary()
-		    && cell->face(face)->boundary_indicator() == 1) {
-		  fe_values_face.reinit(cell, face);
-
-		  unsigned int contact_counter = 0;
-		  std::vector dof_indices (dofs_per_face);
-		  cell->face(face)->get_dof_indices (dof_indices);
-		  
-		  for (unsigned int q_point=0; q_point
+    void
+    PlasticityContactProblem::output_contact_force (
+        const unsigned int cycle)
+    {
+      Functions::FEFieldFunction,
+          TrilinosWrappers::MPI::Vector> solution_function(dof_handler,
+          solution);
+      std::cout.precision(10);
+
+      Vector solution_p1(dim);
+      std::vector > solution_gradient_p1(dim);
+
+      // Here we calculate the contact pressure as a vector lambda.
+      // If a dof is element of the active set lambda contains the
+      // nonlinear residual this dof divided by the according entry
+      // of the mass matrix. In all other dofs lambda will be set to
+      // zero.
+      TrilinosWrappers::MPI::Vector lambda(solution);
+      TrilinosWrappers::MPI::Vector distributed_lambda(system_rhs_newton);
+      const unsigned int start_res = (resid_vector.local_range().first),
+          end_res = (resid_vector.local_range().second);
+      for (unsigned int n = start_res; n < end_res; ++n)
+        if (constraints.is_inhomogeneously_constrained(n))
+          distributed_lambda(n) = resid_vector(n) / diag_mass_matrix_vector(n);
+        else
+          distributed_lambda(n) = 0;
+      distributed_lambda.compress(VectorOperation::insert);
+      constraints_hanging_nodes.distribute(distributed_lambda);
+      lambda = distributed_lambda;
+      Functions::FEFieldFunction,
+          TrilinosWrappers::MPI::Vector> lambda_function(dof_handler, lambda);
+
+      // Here we try to find the MPI-process which owns the cell
+      // with the point_of_interest. If it is the wrong MPI-process
+      // we catch this case and set point_found to false.
+      const Point point_of_interest(0.49, 0.5001, 1.0);
+      Vector contact_pressure_in_point(dim);
+      bool point_found = true;
+
+      MPI_Barrier(MPI_COMM_WORLD);
+      try
+        {
+          lambda_function.vector_value(point_of_interest,
+              contact_pressure_in_point);
+        }
+      catch (const typename Functions::FEFieldFunction,
+          TrilinosWrappers::MPI::Vector>::ExcPointNotAvailableHere &)
+        {
+          point_found = false;
+        }
+
+      if (point_found == true)
+        {
+          std::cout << "PoI contact pressure: " << contact_pressure_in_point(2)
+              << std::endl;
+        }
+
+      // To obtain the contact force we have to compute an integral of the contact pressure
+      // in z-direction over the whole contact area. To be accurate enough we use the
+      // Gaussian quadrature rule with fe.degree + 1.
+      double contact_force = 0.0;
+        {
+          QGauss face_quadrature_formula(fe.degree + 1);
+
+          FEFaceValues fe_values_face(fe, face_quadrature_formula,
+              update_values | update_quadrature_points | update_JxW_values);
+
+          const unsigned int n_face_q_points = face_quadrature_formula.size();
+
+          const FEValuesExtractors::Vector displacement(0);
+
+          typename DoFHandler::active_cell_iterator cell =
+              dof_handler.begin_active(), endc = dof_handler.end();
+          for (; cell != endc; ++cell)
+            if (cell->is_locally_owned())
+              for (unsigned int face = 0;
+                  face < GeometryInfo::faces_per_cell; ++face)
+                if (cell->face(face)->at_boundary()
+                    && cell->face(face)->boundary_indicator() == 1)
+                  {
+                    fe_values_face.reinit(cell, face);
+
+                    std::vector > lambda_values(n_face_q_points);
+                    fe_values_face[displacement].get_function_values(lambda,
+                        lambda_values);
+
+                    for (unsigned int q_point = 0; q_point < n_face_q_points;
+                        ++q_point)
+                      {
+                        contact_force += lambda_values[q_point][2]
+                            * fe_values_face.JxW(q_point);
+                      }
+                  }
+          contact_force = Utilities::MPI::sum(contact_force, MPI_COMM_WORLD);
+          pcout << "Contact force = " << contact_force << std::endl;
+        }
+      MPI_Barrier(MPI_COMM_WORLD);
+    }
 
 // @sect4{PlasticityContactProblem::run}
 
-template
-void PlasticityContactProblem::run() {
-
-	if (obstacle_filename!="")
-	  {
-                  pcout << "Read the obstacle from '" << obstacle_filename
-                      << "' ... " << std::flush;
-                  input_obstacle.reset(new Input(obstacle_filename.c_str()));
-                  pcout << "done." << std::endl;
-          }
+  template 
+    void
+    PlasticityContactProblem::run ()
+    {
 
-	computing_timer.reset();
-	for (cycle = 0; cycle < n_cycles; ++cycle) {
-		{
-			TimerOutput::Scope t(computing_timer, "Setup");
-
-			pcout << std::endl;
-			pcout << "Cycle " << cycle << ':' << std::endl;
-
-			if (cycle == 0) {
-				make_grid();
-			} else {
-				TimerOutput::Scope t(computing_timer, "Setup: refine mesh");
-				if (transfer_solution)
-				  soltrans.reset (new parallel::distributed::SolutionTransfer(dof_handler));
-				refine_grid();
-			}
-
-			setup_system();
-
-			if (transfer_solution && cycle > 0)
-			{
-				TrilinosWrappers::MPI::Vector distributed_solution(
-						system_rhs_newton);
-				distributed_solution = solution;
-				soltrans->interpolate(distributed_solution);
-				solution = distributed_solution;
-				residual_nl_system(solution);
-				resid_vector = system_rhs_lambda;
-				resid_vector.compress(VectorOperation::insert);
-			}
-
-		}
-
-		solve_newton();
-
-		if (true)    //Utilities::MPI::n_mpi_processes(mpi_communicator) <= 64)
-		{
-			pcout << "      Writing graphical output... " << std::flush;
-
-			TimerOutput::Scope t(computing_timer, "Graphical output");
-
-			std::ostringstream filename_solution;
-			filename_solution << "solution-";
-			filename_solution << Utilities::int_to_string(cycle, 2);
-			output_results(filename_solution.str());
-		}
-
-		computing_timer.print_summary();
-		computing_timer.reset();
-
-		Utilities::System::MemoryStats stats;
-		Utilities::System::get_memory_stats(stats);
-		pcout << "VMPEAK, Resident in kB: " << stats.VmSize << " "
-				<< stats.VmRSS << std::endl;
-
-		if (base_mesh=="box")
-		  output_for_benchmark(cycle);
-	}
-}
+      if (obstacle_filename != "")
+        {
+          pcout << "Read the obstacle from '" << obstacle_filename << "' ... "
+              << std::flush;
+          input_obstacle.reset(new Input(obstacle_filename.c_str()));
+          pcout << "done." << std::endl;
+        }
+
+      computing_timer.reset();
+      for (cycle = 0; cycle < n_cycles; ++cycle)
+        {
+            {
+              TimerOutput::Scope t(computing_timer, "Setup");
+
+              pcout << std::endl;
+              pcout << "Cycle " << cycle << ':' << std::endl;
+
+              if (cycle == 0)
+                {
+                  make_grid();
+                }
+              else
+                {
+                  TimerOutput::Scope t(computing_timer, "Setup: refine mesh");
+                  if (transfer_solution)
+                    soltrans.reset(
+                        new parallel::distributed::SolutionTransfer(dof_handler));
+                  refine_grid();
+                }
+
+              setup_system();
+
+              if (transfer_solution && cycle > 0)
+                {
+                  TrilinosWrappers::MPI::Vector distributed_solution(
+                      system_rhs_newton);
+                  distributed_solution = solution;
+                  soltrans->interpolate(distributed_solution);
+                  solution = distributed_solution;
+                  residual_nl_system(solution);
+                  resid_vector = system_rhs_lambda;
+                  resid_vector.compress(VectorOperation::insert);
+                }
+
+            }
+
+          solve_newton();
+
+          if (true) //Utilities::MPI::n_mpi_processes(mpi_communicator) <= 64)
+            {
+              pcout << "      Writing graphical output... " << std::flush;
+
+              TimerOutput::Scope t(computing_timer, "Graphical output");
+
+              std::ostringstream filename_solution;
+              filename_solution << "solution-";
+              filename_solution << Utilities::int_to_string(cycle, 2);
+              output_results(filename_solution.str());
+            }
+
+          computing_timer.print_summary();
+          computing_timer.reset();
+
+          Utilities::System::MemoryStats stats;
+          Utilities::System::get_memory_stats(stats);
+          pcout << "VMPEAK, Resident in kB: " << stats.VmSize << " "
+              << stats.VmRSS << std::endl;
+
+          if (base_mesh == "box")
+            output_contact_force(cycle);
+        }
+    }
 }
 
 // @sect3{The main function}
 
-int main(int argc, char *argv[]) {
-	using namespace dealii;
-	using namespace Step42;
-
-	deallog.depth_console(0);
-	ParameterHandler prm;
-	PlasticityContactProblem<3>::declare(prm);
-	if (argc!=2)
-	{
-	    prm.print_parameters(std::cout, ParameterHandler::Text);
-	    return 0;
-	}
-
-	prm.read_input(argv[1]);
-	Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
-	{
-	  PlasticityContactProblem<3> problem(prm);
-	  problem.run();
-	}
-
-	return 0;
+int
+main (
+    int argc, char *argv[])
+{
+  using namespace dealii;
+  using namespace Step42;
+
+  deallog.depth_console(0);
+  ParameterHandler prm;
+  PlasticityContactProblem<3>::declare(prm);
+  if (argc != 2)
+    {
+      prm.print_parameters(std::cout, ParameterHandler::Text);
+      return 0;
+    }
+
+  prm.read_input(argv[1]);
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
+    {
+      PlasticityContactProblem<3> problem(prm);
+      problem.run();
+    }
+
+  return 0;
 }
diff --git a/deal.II/examples/step-51/doc/intro.dox b/deal.II/examples/step-51/doc/intro.dox
index 22dc984e08..1c14fa37c9 100644
--- a/deal.II/examples/step-51/doc/intro.dox
+++ b/deal.II/examples/step-51/doc/intro.dox
@@ -37,14 +37,14 @@ The Dirichlet-to-Neumann map concept then permits the solution procedure:
 
        
   
      1.  Solve a Dirichlet problem on each element individually, where the local/interior
 element solutions no longer couple to neighboring elements.  This is known as the
-'local solver', and is an embarrissingly parallel element-by-element 
+'local solver', and is an embarrassingly parallel element-by-element 
 solution process.
-  
      2.   Use the local element solutions to construct the global problem for the 
+  
      3.   Use the local element equations to construct the global problem for the 
 trace values.  These are the only globally coupled unknowns.
-  
      4.   Reconstruct the local solutions from the computed trace values.
+  
      5.   Solve for the element-local solutions using the computed trace values.
 
      
 
-The above procedure also has a linear algebra interpretation and referred to
+The above procedure also has a linear algebra interpretation and is referred to
 as static condensation. Let us write the complete linear system associated to
 the HDG problem as a block system with the discrete DG variables U as
 first block and the skeleton variables $\Lambda$ as the second block:
@@ -73,7 +73,7 @@ degrees of freedom of the interior DG variables which are always only related
 to a single cell. The coupling to other cells is introduced by the matrices
 B and C over the skeleton variable. The block-diagonality of
 A and the structure in B and C allow us to invert the
-matrix A element by element (the local solution of the Dirichelt
+matrix A element by element (the local solution of the Dirichlet
 problem) and subtract $CA^{-1}B$ from $D$.
 
 
       Solution quality and rates of convergence
      
@@ -103,7 +103,7 @@ We consider the convection-diffusion equation:
 	\quad \text{ in } \Omega, \\
 	u &=& g_D, \quad \text{ on } \partial \Omega_D, \\
 	(\mathbf{c} u - \kappa \nabla u)\cdot \mathbf{n} &=& g_N, 
-	\quad \text{ on }  \partial \Omega_D.
+	\quad \text{ on }  \partial \Omega_N.
 @f}
 
 Introduce the auxiliary variable $\mathbf{q}=-\kappa \nabla u$ and rewrite
@@ -113,7 +113,62 @@ the above equation as the first order system:
   \nabla \cdot (\mathbf{c} u + \mathbf{q}) &=& f, \quad \text{ in } \Omega, \\
   u &=& g_D, \quad \text{ on } \partial \Omega_D, \\
   (\mathbf{q} + \mathbf{c}u)\cdot\mathbf{n}  &=& g_N, 
-	\quad \text{ on }  \partial \Omega_D.
+	\quad \text{ on }  \partial \Omega_N.
+@f}
+
+We multiply these equations by the weight functions $\mathbf{v}, w$ 
+and integrate by parts over every element $K$ to obtain:
+@f{eqnarray*}
+  (\kappa^{-1} \mathbf{q}, \mathbf{v})_K - (u, \nabla\cdot\mathbf{v})_K
+    + <\hat{u}, \mathbf{v}\cdot\mathbf{n}>_{\partial K} &=& 0, \\
+  - (\mathbf{c} u + \mathbf{q}, \nabla w)_K 
+    + <(\hat{\mathbf{c} u}+\hat{\mathbf{q}})\cdot\mathbf{n}, w>_{\partial K}
+    &=& (f,w)_K.
+@f}
+
+The terms decorated with a $\hat{}$ denote the numerical traces (also commonly referred
+to as numerical fluxes).   They are approximations
+to the interior values on the boundary of the element.  To ensure conservation, 
+these terms must be single-valued on any given element edge $\partial K$.
+We eliminate the numerical trace $\hat{\mathbf{q}}$ by using traces of the form:
+@f{eqnarray*}
+  \hat{\mathbf{c} u}+\hat{\mathbf{q}} = \mathbf{c}\hat{u} + \mathbf{q}
+  + \tau(u - \hat{u})\mathbf{n} \quad \text{ on } \partial K.
+@f}
+The local stabilization parameter $\tau$ has effects on stability and accuracy of 
+HDG solutions; see the literature for a further discussion.
+
+The trace/skeleton variables in HDG methods are single-valued on element faces.  As such,
+they must strongly represent the Dirichlet data on $\partial\Omega_D$.  We introduce
+a new variable $\lambda$ such that
+@f{equation*}
+  \hat{u} = 
+  \begin{cases}
+     g_D  & \text{ on } \partial \Omega_D, \\
+     \lambda & \text{otherwise}.
+\end{cases}
+@f} 
+
+Eliminating $\hat{u}$ from the weak form in favor of $\lambda$, summing the elemental
+contributions across all elements in the triangulation, and enforcing the normal
+component of the numerical flux, we arrive at the final form of the problem:  
+Find $(\mathbf{q}, u, \lambda)$ such that
+@f{eqnarray*}
+  (\kappa^{-1} \mathbf{q}, \mathbf{v})_{\mathcal{T}}
+    - (u, \nabla\cdot\mathbf{v})_{\mathcal{T}}
+    + <\lambda, \mathbf{v}\cdot\mathbf{n}>_{\partial\mathcal{T}}
+    &=&
+    - _{\partial\Omega_D},
+\\
+  - (\mathbf{c} u + \mathbf{q}, \nabla w)_{\mathcal{T}}
+    + <(\hat{\mathbf{c}u} + \hat{\mathbf{q}})\cdot\mathbf{n}, w>_{\partial \mathcal{T}}
+    &=&
+    (f, w)_{\mathcal{T}},
+\\
+  < \left[ \negthinspace \left[ (\hat{\mathbf{c}u} + \hat{\mathbf{q}})\cdot\mathbf{n}
+    \right] \negthinspace \right], \mu>_{\partial \mathcal{T}}
+    &=&
+    _{\partial\Omega_N}
 @f}
 
 
diff --git a/deal.II/examples/step-51/step-51.cc b/deal.II/examples/step-51/step-51.cc
index 9049e88e1a..717356f1c7 100644
--- a/deal.II/examples/step-51/step-51.cc
+++ b/deal.II/examples/step-51/step-51.cc
@@ -62,12 +62,14 @@
 // on codimension-2 surfaces.
 #include 
 
-// The second new file we include defines a new type of sparse matrix.
-// The regular SparseMatrix type stores indices to all non-zero entries.
-// The ChunkSparseMatrix takes advantage of the coupled nature of
-// DG solutions.  It stores an index to a matrix sub-block of a specified
-// size.  In the HDG context, this sub-block-size is actually the number
-// of degrees of freedom per face defined by the skeleton solution field.
+// The second new file we include defines a new type of sparse matrix.  The
+// regular SparseMatrix type stores indices to all non-zero
+// entries.  The ChunkSparseMatrix takes advantage of the coupled
+// nature of DG solutions.  It stores an index to a matrix sub-block of a
+// specified size.  In the HDG context, this sub-block-size is actually the
+// number of degrees of freedom per face defined by the skeleton solution
+// field. This reduces the memory consumption of the matrix by up to one third
+// and results in similar speedups when using the matrix in solvers.
 #include 
 
 // The final new include for this example deals with data output.  Since
@@ -82,10 +84,10 @@ using namespace dealii;
 
 // @sect3{Equation data}
 //
-// The structure of the analytic solution is the same as in step-7. There
-// are two exceptions. Firstly, we also create a solution for the 3d case,
-// and secondly, we take into account the convection velocity in the right
-// hand side that is variable in this case.
+// The structure of the analytic solution is the same as in step-7. There are
+// two exceptions. Firstly, we also create a solution for the 3d case, and
+// secondly, we scale the solution so its norm is of order unity for all
+// values of the solution width.
 template 
 class SolutionBase
 {
@@ -102,7 +104,7 @@ SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers]
 = { Point<1>(-1.0 / 3.0),
     Point<1>(0.0),
     Point<1>(+1.0 / 3.0)
-};
+  };
 
 
 template <>
@@ -111,60 +113,23 @@ SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers]
 = { Point<2>(-0.5, +0.5),
     Point<2>(-0.5, -0.5),
     Point<2>(+0.5, -0.5)
-};
+  };
 
 template <>
 const Point<3>
 SolutionBase<3>::source_centers[SolutionBase<3>::n_source_centers]
 = { Point<3>(-0.5, +0.5, 0.25),
     Point<3>(-0.6, -0.5, -0.125),
-    Point<3>(+0.5, -0.5, 0.5)   };
+    Point<3>(+0.5, -0.5, 0.5)
+  };
 
 template 
 const double SolutionBase::width = 1./5.;
 
 
-
-template 
-class ConvectionVelocity : public TensorFunction<1,dim>
-{
-public:
-  ConvectionVelocity() : TensorFunction<1,dim>() {}
-
-  virtual Tensor<1,dim> value (const Point &p) const;
-};
-
-
-
-template 
-Tensor<1,dim>
-ConvectionVelocity::value(const Point &p) const
-{
-  Tensor<1,dim> convection;
-  switch (dim)
-    {
-    case 1:
-      convection[0] = 1;
-      break;
-    case 2:
-      convection[0] = p[1];
-      convection[1] = -p[0];
-      break;
-    case 3:
-      convection[0] = p[1];
-      convection[1] = -p[0];
-      convection[2] = 1;
-      break;
-    default:
-      Assert(false, ExcNotImplemented());
-    }
-  return convection;
-}
-
-
 template 
 class Solution : public Function,
-                 protected SolutionBase
+  protected SolutionBase
 {
 public:
   Solution () : Function() {}
@@ -191,7 +156,7 @@ double Solution::value (const Point   &p,
     }
 
   return return_value /
-    Utilities::fixed_power(std::sqrt(2. * numbers::PI) * this->width);
+         Utilities::fixed_power(std::sqrt(2. * numbers::PI) * this->width);
 }
 
 
@@ -218,30 +183,84 @@ Tensor<1,dim> Solution::gradient (const Point   &p,
 
 
 
+// This class implements a function where the scalar solution and its negative
+// gradient are collected together. This function is used when computing the
+// error of the HDG approximation and its implementation is to simply call
+// value and gradient function of the Solution class.
 template 
 class SolutionAndGradient : public Function,
-                            protected SolutionBase
+  protected SolutionBase
 {
 public:
   SolutionAndGradient () : Function(dim) {}
 
   virtual void vector_value (const Point   &p,
-                             Vector     &v) const
-  {
-    AssertDimension(v.size(), dim+1);
-    Solution solution;
-    Tensor<1,dim> grad = solution.gradient(p);
-    for (unsigned int d=0; d     &v) const;
+};
+
+template 
+void SolutionAndGradient::vector_value (const Point &p,
+                                             Vector   &v) const
+{
+  AssertDimension(v.size(), dim+1);
+  Solution solution;
+  Tensor<1,dim> grad = solution.gradient(p);
+  for (unsigned int d=0; d
+class ConvectionVelocity : public TensorFunction<1,dim>
+{
+public:
+  ConvectionVelocity() : TensorFunction<1,dim>() {}
+
+  virtual Tensor<1,dim> value (const Point &p) const;
 };
 
 
 
+template 
+Tensor<1,dim>
+ConvectionVelocity::value(const Point &p) const
+{
+  Tensor<1,dim> convection;
+  switch (dim)
+    {
+    case 1:
+      convection[0] = 1;
+      break;
+    case 2:
+      convection[0] = p[1];
+      convection[1] = -p[0];
+      break;
+    case 3:
+      convection[0] = p[1];
+      convection[1] = -p[0];
+      convection[2] = 1;
+      break;
+    default:
+      Assert(false, ExcNotImplemented());
+    }
+  return convection;
+}
+
+
+
+// The last function we implement is the right hand side for the manufactured
+// solution. It is very similar to step-7, with the exception that we now have
+// a convection term instead of the reaction term. Since the velocity field is
+// incompressible, i.e. $\nabla \cdot \mathbf{c} = 0$, this term simply reads
+// $\mathbf{c} \nabla \ve u$.
 template 
 class RightHandSide : public Function,
-                      protected SolutionBase
+  protected SolutionBase
 {
 public:
   RightHandSide () : Function() {}
@@ -278,99 +297,97 @@ double RightHandSide::value (const Point   &p,
 
 // @sect3{The Step51 HDG solver class}
 
-// The HDG solution procedure follows closely that of step-7.  The major
-// difference is the use of 3 different sets of DoFHandler and FE objects,
-// along with the ChunkSparseMatrix and the corresponding solutions vectors.
-// We also use WorkStream to enable a multi-threaded local solution process which exploits
-// the embarrassingly parallel nature of the local solver.
-
+// The HDG solution procedure follows closely that of step-7. The major
+// difference is the use of 3 different sets of DoFHandler and FE
+// objects, along with the ChunkSparseMatrix and the
+// corresponding solutions vectors. We also use WorkStream to enable a
+// multi-threaded local solution process which exploits the embarrassingly
+// parallel nature of the local solver. For WorkStream, we define the local
+// operations on a cell and a copy function into the global matrix and
+// vector. We do this once for the assembly (which is run twice, once when we
+// generate the system matrix and once when we compute the element-interior
+// solutions from the skeleton values) and once for the postprocessing where
+// we extract a solution that converges at higher order.
 template 
 class Step51
 {
 public:
   enum RefinementMode
-    {
-      global_refinement, adaptive_refinement
-    };
+  {
+    global_refinement, adaptive_refinement
+  };
 
   Step51 (const unsigned int degree,
           const RefinementMode refinement_mode);
   void run ();
 
 private:
-    
+
   struct PerTaskData;
   struct ScratchData;
-  
-  struct EmptyData;
+
   struct PostProcessScratchData;
-  void NullFunction(const EmptyData &data);
-    
+
   void setup_system ();
   void assemble_system (const bool reconstruct_trace = false);
   void assemble_system_one_cell (const typename DoFHandler::active_cell_iterator &cell,
-                                   ScratchData &scratch,
-                                   PerTaskData &task_data);
+                                 ScratchData &scratch,
+                                 PerTaskData &task_data);
   void copy_local_to_global(const PerTaskData &data);
   void solve ();
   void postprocess ();
   void postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell,
-                                   PostProcessScratchData &scratch,
-                                   EmptyData &task_data);
+                             PostProcessScratchData &scratch,
+                             unsigned int &empty_data);
   void refine_grid (const unsigned int cylce);
   void output_results (const unsigned int cycle);
 
   Triangulation   triangulation;
-  
-// The 'local' solutions are interior to each element.  These
-// represent the primal solution field $u$ as well as the auxiliary
-// field $\mathbf{q} = \nabla u$. 
+
+  // The 'local' solutions are interior to each element.  These
+  // represent the primal solution field $u$ as well as the auxiliary
+  // field $\mathbf{q} = -\nabla u$.
   FESystem        fe_local;
   DoFHandler      dof_handler_local;
+  Vector       solution_local;
 
-// The new finite element type and corresponding DoFHandler
-// are used for the global solution that couples the element-level local
-// solution.
+  // The new finite element type and corresponding DoFHandler are
+  // used for the global solution that couples the element-level local
+  // solutions.
   FE_FaceQ        fe;
   DoFHandler      dof_handler;
+  Vector       solution;
+  Vector       system_rhs;
 
-// As stated in the introduction, HDG solutions can be post-processed to
-// attain superconvegence rates of $\mathcal{O}(h^{p+2})$.
-// The post-processed solution is a discontinuous finite element solution
-// representing the primal variable on the interior of each cell.
-// We define a FE type to represent this post-processed solution, which we
-// only use for output after constructing it.
+  // As stated in the introduction, HDG solutions can be post-processed to
+  // attain superconvegence rates of $\mathcal{O}(h^{p+2})$.  The
+  // post-processed solution is a discontinuous finite element solution
+  // representing the primal variable on the interior of each cell.  We define
+  // a FE type of degree $p+1$ to represent this post-processed solution,
+  // which we only use for output after constructing it.
   FE_DGQ          fe_u_post;
   DoFHandler      dof_handler_u_post;
+  Vector       solution_u_post;
 
-// The degrees of freedom corresponding to the skeleton strongly enforce
-// Dirichlet boundary conditions, just as in a continuous Galerkin finite
-// element method.  We can enforce the boundary conditions in an analogous
-// manner through the use of ConstrainMatrix constructs.
+  // The degrees of freedom corresponding to the skeleton strongly enforce
+  // Dirichlet boundary conditions, just as in a continuous Galerkin finite
+  // element method.  We can enforce the boundary conditions in an analogous
+  // manner through the use of ConstrainMatrix constructs. In
+  // addition, hanging nodes where cells of different refinement levels meet
+  // are set as for continuous finite elements: For the face elements which
+  // only define degrees of freedom on the face, this process sets the
+  // solution on the refined to be the one from the coarse side.
   ConstraintMatrix     constraints;
-  
-  // Comment on chunk.
+
+  // The usage of the ChunkSparseMatrix class is similar to the usual sparse
+  // matrices: You need a sparsity pattern of type ChunkSparsityPattern and
+  // the actual matrix object. When creating the sparsity pattern, we just
+  // have to additionally pass the size of local blocks.
   ChunkSparsityPattern sparsity_pattern;
   ChunkSparseMatrix system_matrix;
 
-  // Global/skeleton solution/rhs
-  Vector       solution;
-  Vector       system_rhs;
-
-  // Local elementwise solution
-  Vector       solution_local;
-  
-  // HDG solutions can be post-processed
-  // to gain one order of accuracy.
-  // solution_u_post will be
-  // our post-processed DG solution on the
-  // interior of cells represented by a 
-  // DG solution of order (degree+1)
-  Vector       solution_u_post;
-
   // Same as step-7:
   const RefinementMode refinement_mode;
-
   ConvergenceTable     convergence_table;
 };
 
@@ -389,137 +406,141 @@ Step51::Step51 (const unsigned int degree,
   refinement_mode (refinement_mode)
 {}
 
+
+
+// First come the definition of the local data structures for the parallel
+// assembly. The first structure @p PerTaskData contains the local vector and
+// matrix that are written into the global matrix, whereas the ScratchData
+// contains all data that we need for the local assembly.
 template 
 struct Step51::PerTaskData
 {
-    FullMatrix cell_matrix;
-    Vector     cell_vector;
-    std::vector dof_indices;
-    
-    bool trace_reconstruct;
-    
-    PerTaskData(const unsigned int n_dofs, const bool trace_reconstruct)
+  FullMatrix cell_matrix;
+  Vector     cell_vector;
+  std::vector dof_indices;
+
+  bool trace_reconstruct;
+
+  PerTaskData(const unsigned int n_dofs, const bool trace_reconstruct)
     : cell_matrix(n_dofs, n_dofs),
       cell_vector(n_dofs),
       dof_indices(n_dofs),
       trace_reconstruct(trace_reconstruct)
-    {}
-    
-    void reset(){
-        cell_matrix = 0.0;
-        cell_vector = 0.0;
-    }
+  {}
+
+  void reset()
+  {
+    cell_matrix = 0.0;
+    cell_vector = 0.0;
+  }
 };
 
+
+
 template 
 struct Step51::ScratchData
 {
-    FEValues     fe_values_local;
-    FEFaceValues fe_face_values_local;
-    FEFaceValues fe_face_values;
-    
-    FullMatrix ll_matrix;
-    FullMatrix lf_matrix;
-    FullMatrix fl_matrix;
-    FullMatrix tmp_matrix;
-    Vector     l_rhs;
-    Vector     tmp_rhs;
-    
-    std::vector > q_phi;
-    std::vector         q_phi_div;
-    std::vector         u_phi;
-    std::vector > u_phi_grad;
-    std::vector         tr_phi;
-    std::vector         trace_values;
-    
-    std::vector > fe_local_support_on_face;
-    std::vector > fe_support_on_face;
-    
-    ConvectionVelocity convection_velocity;
-    RightHandSide right_hand_side;
-    const Solution exact_solution;
-    
-        // Full constructor
-    ScratchData(const FiniteElement &fe,
-                const FiniteElement &fe_local,
-                const QGauss   &quadrature_formula,
-                const QGauss &face_quadrature_formula,
-                const UpdateFlags local_flags,
-                const UpdateFlags local_face_flags,
-                const UpdateFlags flags)
-      :
-      fe_values_local (fe_local, quadrature_formula, local_flags),
-      fe_face_values_local (fe_local, face_quadrature_formula, local_face_flags),
-      fe_face_values (fe, face_quadrature_formula, flags),
-      ll_matrix (fe_local.dofs_per_cell, fe_local.dofs_per_cell),
-      lf_matrix (fe_local.dofs_per_cell, fe.dofs_per_cell),
-      fl_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell),
-      tmp_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell),
-      l_rhs (fe_local.dofs_per_cell),
-      tmp_rhs (fe_local.dofs_per_cell),
-      q_phi (fe_local.dofs_per_cell),
-      q_phi_div (fe_local.dofs_per_cell),
-      u_phi (fe_local.dofs_per_cell),
-      u_phi_grad (fe_local.dofs_per_cell),
-      tr_phi (fe.dofs_per_cell),
-      trace_values(face_quadrature_formula.size()),
-      fe_local_support_on_face(GeometryInfo::faces_per_cell),
-      fe_support_on_face(GeometryInfo::faces_per_cell)
-    {
-        for (unsigned int face=0; face::faces_per_cell; ++face)
-            for (unsigned int i=0; i::faces_per_cell; ++face)
-            for (unsigned int i=0; ireset(), but
-// need to have it for the WorkStream interface.   
-    void reset() {}
-    
-};
+  FEValues     fe_values_local;
+  FEFaceValues fe_face_values_local;
+  FEFaceValues fe_face_values;
+
+  FullMatrix ll_matrix;
+  FullMatrix lf_matrix;
+  FullMatrix fl_matrix;
+  FullMatrix tmp_matrix;
+  Vector     l_rhs;
+  Vector     tmp_rhs;
+
+  std::vector > q_phi;
+  std::vector         q_phi_div;
+  std::vector         u_phi;
+  std::vector > u_phi_grad;
+  std::vector         tr_phi;
+  std::vector         trace_values;
+
+  std::vector > fe_local_support_on_face;
+  std::vector > fe_support_on_face;
+
+  ConvectionVelocity convection_velocity;
+  RightHandSide right_hand_side;
+  const Solution exact_solution;
+
+  // Full constructor
+  ScratchData(const FiniteElement &fe,
+              const FiniteElement &fe_local,
+              const QGauss   &quadrature_formula,
+              const QGauss &face_quadrature_formula,
+              const UpdateFlags local_flags,
+              const UpdateFlags local_face_flags,
+              const UpdateFlags flags)
+    :
+    fe_values_local (fe_local, quadrature_formula, local_flags),
+    fe_face_values_local (fe_local, face_quadrature_formula, local_face_flags),
+    fe_face_values (fe, face_quadrature_formula, flags),
+    ll_matrix (fe_local.dofs_per_cell, fe_local.dofs_per_cell),
+    lf_matrix (fe_local.dofs_per_cell, fe.dofs_per_cell),
+    fl_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell),
+    tmp_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell),
+    l_rhs (fe_local.dofs_per_cell),
+    tmp_rhs (fe_local.dofs_per_cell),
+    q_phi (fe_local.dofs_per_cell),
+    q_phi_div (fe_local.dofs_per_cell),
+    u_phi (fe_local.dofs_per_cell),
+    u_phi_grad (fe_local.dofs_per_cell),
+    tr_phi (fe.dofs_per_cell),
+    trace_values(face_quadrature_formula.size()),
+    fe_local_support_on_face(GeometryInfo::faces_per_cell),
+    fe_support_on_face(GeometryInfo::faces_per_cell)
+  {
+    for (unsigned int face=0; face::faces_per_cell; ++face)
+      for (unsigned int i=0; i::faces_per_cell; ++face)
+      for (unsigned int i=0; ireset(), but
+  // need to have it for the WorkStream interface.
+  void reset() {}
 
-template 
-struct Step51::EmptyData
-{
-    EmptyData(){}
-    void reset(){}
 };
 
 template 
@@ -527,65 +548,63 @@ struct Step51::PostProcessScratchData
 {
   FEValues fe_values_local;
   FEValues fe_values;
-    
+
   std::vector u_values;
   std::vector > u_gradients;
   FullMatrix cell_matrix;
-  
+
   Vector cell_rhs;
   Vector cell_sol;
-    
-        // Full constructor
+
+  // Full constructor
   PostProcessScratchData(const FiniteElement &fe,
-                const FiniteElement &fe_local,
-                const QGauss   &quadrature_formula,
-                const UpdateFlags local_flags,
-                const UpdateFlags flags)
-      :
-      fe_values_local (fe_local, quadrature_formula, local_flags),
-      fe_values (fe, quadrature_formula, flags),
-      u_values (quadrature_formula.size()),
-      u_gradients (quadrature_formula.size()),
-      cell_matrix (fe.dofs_per_cell, fe.dofs_per_cell),
-      cell_rhs (fe.dofs_per_cell),
-      cell_sol (fe.dofs_per_cell)
-      {}
-    
-        // Copy constructor
+                         const FiniteElement &fe_local,
+                         const QGauss   &quadrature_formula,
+                         const UpdateFlags local_flags,
+                         const UpdateFlags flags)
+    :
+    fe_values_local (fe_local, quadrature_formula, local_flags),
+    fe_values (fe, quadrature_formula, flags),
+    u_values (quadrature_formula.size()),
+    u_gradients (quadrature_formula.size()),
+    cell_matrix (fe.dofs_per_cell, fe.dofs_per_cell),
+    cell_rhs (fe.dofs_per_cell),
+    cell_sol (fe.dofs_per_cell)
+  {}
+
+  // Copy constructor
   PostProcessScratchData(const PostProcessScratchData &sd)
-      :
-      fe_values_local (sd.fe_values_local.get_fe(),
-                       sd.fe_values_local.get_quadrature(),
-                       sd.fe_values_local.get_update_flags()),
-      fe_values (sd.fe_values.get_fe(),
-                      sd.fe_values.get_quadrature(),
-                      sd.fe_values.get_update_flags()),
-      u_values (sd.u_values),
-      u_gradients (sd.u_gradients),
-      cell_matrix (sd.cell_matrix),
-      cell_rhs (sd.cell_rhs),
-      cell_sol (sd.cell_sol)
-    {}
-       
-    void reset() {
-      cell_matrix = 0.;
-      cell_rhs = 0.;
-      cell_sol = 0.;
-    }
-    
-};
+    :
+    fe_values_local (sd.fe_values_local.get_fe(),
+                     sd.fe_values_local.get_quadrature(),
+                     sd.fe_values_local.get_update_flags()),
+    fe_values (sd.fe_values.get_fe(),
+               sd.fe_values.get_quadrature(),
+               sd.fe_values.get_update_flags()),
+    u_values (sd.u_values),
+    u_gradients (sd.u_gradients),
+    cell_matrix (sd.cell_matrix),
+    cell_rhs (sd.cell_rhs),
+    cell_sol (sd.cell_sol)
+  {}
+
+  void reset()
+  {
+    cell_matrix = 0.;
+    cell_rhs = 0.;
+    cell_sol = 0.;
+  }
 
-template 
-void Step51::NullFunction(const EmptyData &data){}
+};
 
 template 
 void Step51::copy_local_to_global(const PerTaskData &data)
 {
-    if(data.trace_reconstruct == false)
-      constraints.distribute_local_to_global (data.cell_matrix,
-                                              data.cell_vector,
-                                              data.dof_indices,
-                                              system_matrix, system_rhs);
+  if (data.trace_reconstruct == false)
+    constraints.distribute_local_to_global (data.cell_matrix,
+                                            data.cell_vector,
+                                            data.dof_indices,
+                                            system_matrix, system_rhs);
 }
 
 template 
@@ -632,34 +651,34 @@ template 
 void
 Step51::assemble_system (const bool trace_reconstruct)
 {
-    const QGauss   quadrature_formula(fe.degree+1);
-    const QGauss face_quadrature_formula(fe.degree+1);
-    
-    const UpdateFlags local_flags (update_values | update_gradients |
-                                   update_JxW_values | update_quadrature_points);
-    
-    const UpdateFlags local_face_flags (update_values);
-    
-    const UpdateFlags flags ( update_values | update_normal_vectors |
-                             update_quadrature_points |
-                             update_JxW_values);
-    
-    PerTaskData task_data (fe.dofs_per_cell,
-                           trace_reconstruct);
-    ScratchData scratch (fe, fe_local,
-                         quadrature_formula,
-                         face_quadrature_formula,
-                         local_flags,
-                         local_face_flags,
-                         flags);
-    
-    WorkStream::run(dof_handler.begin_active(),
-                    dof_handler.end(),
-                    *this,
-                    &Step51::assemble_system_one_cell,
-                    &Step51::copy_local_to_global,
-                    scratch,
-                    task_data);
+  const QGauss   quadrature_formula(fe.degree+1);
+  const QGauss face_quadrature_formula(fe.degree+1);
+
+  const UpdateFlags local_flags (update_values | update_gradients |
+                                 update_JxW_values | update_quadrature_points);
+
+  const UpdateFlags local_face_flags (update_values);
+
+  const UpdateFlags flags ( update_values | update_normal_vectors |
+                            update_quadrature_points |
+                            update_JxW_values);
+
+  PerTaskData task_data (fe.dofs_per_cell,
+                         trace_reconstruct);
+  ScratchData scratch (fe, fe_local,
+                       quadrature_formula,
+                       face_quadrature_formula,
+                       local_flags,
+                       local_face_flags,
+                       flags);
+
+  WorkStream::run(dof_handler.begin_active(),
+                  dof_handler.end(),
+                  *this,
+                  &Step51::assemble_system_one_cell,
+                  &Step51::copy_local_to_global,
+                  scratch,
+                  task_data);
 }
 
 
@@ -669,17 +688,17 @@ Step51::assemble_system_one_cell (const typename DoFHandler::active_ce
                                        ScratchData &scratch,
                                        PerTaskData &task_data)
 {
-        // Construct iterator for dof_handler_local
+  // Construct iterator for dof_handler_local
   typename DoFHandler::active_cell_iterator
-    loc_cell (&triangulation,
-              cell->level(),
-              cell->index(),
-              &dof_handler_local);
+  loc_cell (&triangulation,
+            cell->level(),
+            cell->index(),
+            &dof_handler_local);
 
   const unsigned int n_q_points    = scratch.fe_values_local.get_quadrature().size();
   const unsigned int n_face_q_points = scratch.fe_face_values_local.get_quadrature().size();
 
-        //  const unsigned int dofs_per_cell = scratch.fe_face_values.get_fe().dofs_per_cell;
+  //  const unsigned int dofs_per_cell = scratch.fe_face_values.get_fe().dofs_per_cell;
   const unsigned int loc_dofs_per_cell = scratch.fe_values_local.get_fe().dofs_per_cell;
 
   // Choose stabilization parameter to be 5 * diffusion = 5
@@ -688,156 +707,156 @@ Step51::assemble_system_one_cell (const typename DoFHandler::active_ce
   const FEValuesExtractors::Vector fluxes (0);
   const FEValuesExtractors::Scalar scalar (dim);
 
-      scratch.ll_matrix = 0;
-      scratch.l_rhs = 0;
-      if (!task_data.trace_reconstruct)
+  scratch.ll_matrix = 0;
+  scratch.l_rhs = 0;
+  if (!task_data.trace_reconstruct)
+    {
+      scratch.lf_matrix = 0;
+      scratch.fl_matrix = 0;
+      task_data.reset();
+    }
+  scratch.fe_values_local.reinit (loc_cell);
+
+  for (unsigned int q=0; q convection
+      = scratch.convection_velocity.value(scratch.fe_values_local.quadrature_point(q));
+      const double JxW = scratch.fe_values_local.JxW(q);
+      for (unsigned int k=0; k convection
-            = scratch.convection_velocity.value(scratch.fe_values_local.quadrature_point(q));
-          const double JxW = scratch.fe_values_local.JxW(q);
-          for (unsigned int k=0; k::faces_per_cell; ++face)
+  for (unsigned int face=0; face::faces_per_cell; ++face)
+    {
+      scratch.fe_face_values_local.reinit(loc_cell, face);
+      scratch.fe_face_values.reinit(cell, face);
+      if (task_data.trace_reconstruct)
+        scratch.fe_face_values.get_function_values (solution, scratch.trace_values);
+
+      for (unsigned int q=0; q quadrature_point =
+            scratch.fe_face_values.quadrature_point(q);
+          const Point normal = scratch.fe_face_values.normal_vector(q);
+          const Tensor<1,dim> convection
+          = scratch.convection_velocity.value(quadrature_point);
+          const double tau_stab = (tau_stab_diffusion +
+                                   std::abs(convection * normal));
 
-          for (unsigned int q=0; q quadrature_point =
-                scratch.fe_face_values.quadrature_point(q);
-              const Point normal = scratch.fe_face_values.normal_vector(q);
-              const Tensor<1,dim> convection
-                = scratch.convection_velocity.value(quadrature_point);
-              const double tau_stab = (tau_stab_diffusion +
-                                       std::abs(convection * normal));
-
-              for (unsigned int k=0; kface(face)->at_boundary()
-                      &&
-                      (cell->face(face)->boundary_indicator() == 1))
-                    {
-                      const double neumann_value =
-                        - scratch.exact_solution.gradient (quadrature_point) * normal
-                        + convection * normal * scratch.exact_solution.value(quadrature_point);
-                      for (unsigned int i=0; iface(face)->at_boundary()
+                  &&
+                  (cell->face(face)->boundary_indicator() == 1))
+                {
+                  const double neumann_value =
+                    - scratch.exact_solution.gradient (quadrature_point) * normal
+                    + convection * normal * scratch.exact_solution.value(quadrature_point);
+                  for (unsigned int i=0; iget_dof_indices(task_data.dof_indices);
-        }
-      else
-        {
-          scratch.ll_matrix.vmult(scratch.tmp_rhs, scratch.l_rhs);
-          loc_cell->set_dof_values(scratch.tmp_rhs, solution_local);
+          for (unsigned int i=0; iget_dof_indices(task_data.dof_indices);
+    }
+  else
+    {
+      scratch.ll_matrix.vmult(scratch.tmp_rhs, scratch.l_rhs);
+      loc_cell->set_dof_values(scratch.tmp_rhs, solution_local);
+    }
 }
 
 
@@ -856,7 +875,7 @@ void Step51::solve ()
 
   system_matrix.clear();
   sparsity_pattern.reinit(0,0,0,1);
-    
+
   constraints.distribute(solution);
 
   // update local values
@@ -875,25 +894,24 @@ Step51::postprocess()
     const QGauss   quadrature_formula(fe_u_post.degree+1);
     const UpdateFlags local_flags (update_values);
     const UpdateFlags flags ( update_values | update_gradients |
-                             update_JxW_values);
-  
-    EmptyData task_data;
-    
+                              update_JxW_values);
+
     PostProcessScratchData scratch (fe_u_post, fe_local,
-                         quadrature_formula,
-                         local_flags,
-                         flags);
-    
+                                    quadrature_formula,
+                                    local_flags,
+                                    flags);
+
     WorkStream::run(dof_handler_u_post.begin_active(),
                     dof_handler_u_post.end(),
-                    *this,
-                    &Step51::postprocess_one_cell,
-                    &Step51::NullFunction,
+                    std_cxx1x::bind (&Step51::postprocess_one_cell,
+                                     std_cxx1x::ref(*this),
+                                     std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3),
+                    std_cxx1x::function(),
                     scratch,
-                    task_data);
-    }
+                    0U);
+  }
 
-// Compute some convergence rates, etc., and add to a table
+  // Compute some convergence rates, etc., and add to a table
   Vector difference_per_cell (triangulation.n_active_cells());
 
   ComponentSelectFunction value_select (dim, dim+1);
@@ -924,31 +942,31 @@ Step51::postprocess()
                                      QGauss(fe.degree+3),
                                      VectorTools::L2_norm);
   const double post_error = difference_per_cell.l2_norm();
-  
-  convergence_table.add_value("cells", 		 triangulation.n_active_cells());
-  convergence_table.add_value("dofs", 	 	 dof_handler.n_dofs());
-  convergence_table.add_value("val L2", 	 L2_error);
-  convergence_table.add_value("grad L2", 	 grad_error);
-  convergence_table.add_value("val L2-post", post_error);  
+
+  convergence_table.add_value("cells",     triangulation.n_active_cells());
+  convergence_table.add_value("dofs",      dof_handler.n_dofs());
+  convergence_table.add_value("val L2",    L2_error);
+  convergence_table.add_value("grad L2",   grad_error);
+  convergence_table.add_value("val L2-post", post_error);
 }
 
 template 
 void
 Step51::postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell,
-                                       PostProcessScratchData &scratch,
-                                       EmptyData &task_data)
+                                   PostProcessScratchData &scratch,
+                                   unsigned int &)
 {
   scratch.reset();
 
   typename DoFHandler::active_cell_iterator
-    loc_cell (&triangulation,
-              cell->level(),
-              cell->index(),
-              &dof_handler_local);
-              
+  loc_cell (&triangulation,
+            cell->level(),
+            cell->index(),
+            &dof_handler_local);
+
   scratch.fe_values_local.reinit (loc_cell);
   scratch.fe_values.reinit(cell);
-  
+
   FEValuesExtractors::Vector fluxes(0);
   FEValuesExtractors::Scalar scalar(dim);
 
@@ -957,38 +975,38 @@ Step51::postprocess_one_cell (const typename DoFHandler::active_cell_i
 
   scratch.fe_values_local[scalar].get_function_values(solution_local, scratch.u_values);
   scratch.fe_values_local[fluxes].get_function_values(solution_local, scratch.u_gradients);
-  
+
   double sum = 0;
   for (unsigned int i=1; i::output_results (const unsigned int cycle)
     default:
       Assert (false, ExcNotImplemented());
     }
-    
+
   std::string face_out(filename);
   face_out += "-face";
 
@@ -1025,24 +1043,24 @@ void Step51::output_results (const unsigned int cycle)
   std::vector names (dim, "gradient");
   names.push_back ("solution");
   std::vector
-    component_interpretation
-    (dim+1, DataComponentInterpretation::component_is_part_of_vector);
+  component_interpretation
+  (dim+1, DataComponentInterpretation::component_is_part_of_vector);
   component_interpretation[dim]
-    = DataComponentInterpretation::component_is_scalar;
+  = DataComponentInterpretation::component_is_scalar;
   data_out.add_data_vector (dof_handler_local, solution_local,
                             names, component_interpretation);
-                            
-  // Post-processed solution:  can now add more than 1 dof_handler to 
+
+  // Post-processed solution:  can now add more than 1 dof_handler to
   // the DataOut object!
   std::vector post_name(1,"u_post");
-  std::vector 
-			 post_comp_type(1, DataComponentInterpretation::component_is_scalar);
+  std::vector
+  post_comp_type(1, DataComponentInterpretation::component_is_scalar);
   data_out.add_data_vector (dof_handler_u_post, solution_u_post,
-						  post_name, post_comp_type);
+                            post_name, post_comp_type);
 
   data_out.build_patches (fe.degree);
   data_out.write_vtk (output);
-    
+
   face_out += "-q" + Utilities::int_to_string(fe.degree,1);
   face_out += "-" + Utilities::int_to_string(cycle,2);
   face_out += ".vtk";
@@ -1050,14 +1068,14 @@ void Step51::output_results (const unsigned int cycle)
 
   DataOutFaces data_out_face(false);
   std::vector face_name(1,"lambda");
-  std::vector 
-			 face_component_type(1, DataComponentInterpretation::component_is_scalar);
-
-  data_out_face.add_data_vector (dof_handler, 
-							  	solution, 
-								face_name,
-								face_component_type);
-								
+  std::vector
+  face_component_type(1, DataComponentInterpretation::component_is_scalar);
+
+  data_out_face.add_data_vector (dof_handler,
+                                 solution,
+                                 face_name,
+                                 face_component_type);
+
   data_out_face.build_patches (fe.degree);
   data_out_face.write_vtk (face_output);
 }
@@ -1076,39 +1094,39 @@ void Step51::refine_grid (const unsigned int cycle)
     switch (refinement_mode)
       {
       case global_refinement:
-        {
-          triangulation.clear();
-          GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1);
-          triangulation.refine_global(3-dim+cycle/2);
-          break;
-        }
+      {
+        triangulation.clear();
+        GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1);
+        triangulation.refine_global(3-dim+cycle/2);
+        break;
+      }
 
       case adaptive_refinement:
-        {
-          Vector estimated_error_per_cell (triangulation.n_active_cells());
+      {
+        Vector estimated_error_per_cell (triangulation.n_active_cells());
 
-          FEValuesExtractors::Scalar scalar(dim);
-          typename FunctionMap::type neumann_boundary;
-          KellyErrorEstimator::estimate (dof_handler_local,
-                                              QGauss(3),
-                                              neumann_boundary,
-                                              solution_local,
-                                              estimated_error_per_cell,
-                                              fe_local.component_mask(scalar));
+        FEValuesExtractors::Scalar scalar(dim);
+        typename FunctionMap::type neumann_boundary;
+        KellyErrorEstimator::estimate (dof_handler_local,
+                                            QGauss(3),
+                                            neumann_boundary,
+                                            solution_local,
+                                            estimated_error_per_cell,
+                                            fe_local.component_mask(scalar));
 
-          GridRefinement::refine_and_coarsen_fixed_number (triangulation,
-                                                           estimated_error_per_cell,
-                                                           0.3, 0.);
+        GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                         estimated_error_per_cell,
+                                                         0.3, 0.);
 
-          triangulation.execute_coarsening_and_refinement ();
+        triangulation.execute_coarsening_and_refinement ();
 
-          break;
-        }
+        break;
+      }
 
       default:
-        {
-          Assert (false, ExcNotImplemented());
-        }
+      {
+        Assert (false, ExcNotImplemented());
+      }
       }
 
   // Just as in step-7, we set the boundary indicator of one of the faces to 1
@@ -1117,8 +1135,8 @@ void Step51::refine_grid (const unsigned int cycle)
   // refinement, the flags are set in every refinement step, not just at the
   // beginning.
   typename Triangulation::cell_iterator
-    cell = triangulation.begin (),
-    endc = triangulation.end();
+  cell = triangulation.begin (),
+  endc = triangulation.end();
   for (; cell!=endc; ++cell)
     for (unsigned int face=0; face::faces_per_cell; ++face)
       if ((std::fabs(cell->face(face)->center()(0) - (-1)) < 1e-12)
@@ -1137,7 +1155,7 @@ void Step51::run ()
   for (unsigned int cycle=0; cycle<10; ++cycle)
     {
       std::cout << "Cycle " << cycle << ':' << std::endl;
-      
+
       refine_grid (cycle);
       setup_system ();
       assemble_system (false);
@@ -1158,17 +1176,17 @@ void Step51::run ()
   if (refinement_mode == global_refinement)
     {
       convergence_table
-        .evaluate_convergence_rates("val L2", "cells", ConvergenceTable::reduction_rate_log2, dim);
+      .evaluate_convergence_rates("val L2", "cells", ConvergenceTable::reduction_rate_log2, dim);
       convergence_table
-        .evaluate_convergence_rates("grad L2", "cells", ConvergenceTable::reduction_rate_log2, dim);
+      .evaluate_convergence_rates("grad L2", "cells", ConvergenceTable::reduction_rate_log2, dim);
       convergence_table
-        .evaluate_convergence_rates("val L2-post", "cells", ConvergenceTable::reduction_rate_log2, dim);
+      .evaluate_convergence_rates("val L2-post", "cells", ConvergenceTable::reduction_rate_log2, dim);
     }
   convergence_table.write_text(std::cout);
 }
 
 
-int main (int argc, char** argv)
+int main (int argc, char **argv)
 {
   const unsigned int dim = 2;
 
diff --git a/deal.II/include/deal.II/base/memory_consumption.h b/deal.II/include/deal.II/base/memory_consumption.h
index f5add851ca..a42f4ca401 100644
--- a/deal.II/include/deal.II/base/memory_consumption.h
+++ b/deal.II/include/deal.II/base/memory_consumption.h
@@ -483,6 +483,14 @@ namespace MemoryConsumption
 
 
 
+  inline
+  std::size_t memory_consumption (const unsigned long int)
+  {
+    return sizeof(unsigned long int);
+  }
+
+
+
   inline
   std::size_t memory_consumption (const unsigned long long int)
   {
diff --git a/deal.II/include/deal.II/base/parameter_handler.h b/deal.II/include/deal.II/base/parameter_handler.h
index 70712ea751..0e6acd4ed0 100644
--- a/deal.II/include/deal.II/base/parameter_handler.h
+++ b/deal.II/include/deal.II/base/parameter_handler.h
@@ -2246,7 +2246,28 @@ public:
    */
   virtual ~MultipleParameterLoop ();
 
-  virtual bool read_input (std::istream &Input);
+  /**
+   * Read input from a stream until the stream returns the eof condition
+   * or error. The second argument can be used to denote the name of the file
+   * (if that's what the input stream represents) we are reading from; this
+   * is only used when creating output for error messages.
+   *
+   * Return whether the read was successful.
+   */
+  virtual bool read_input (std::istream &input,
+                           const std::string &filename = "input file");
+
+  /**
+   * Read input from a file the name of which is given. The PathSearch
+   * class "PARAMETERS" is used to find the file.
+   *
+   * Return whether the read was successful.
+   *
+   * Unless optional is true, this function will
+   * automatically generate the requested file with default values if the
+   * file did not exist. This file will not contain additional comments if
+   * write_stripped_file is true.
+   */
   virtual bool read_input (const std::string &FileName,
                            const bool optional = false,
                            const bool write_stripped_file = false);
diff --git a/deal.II/include/deal.II/base/work_stream.h b/deal.II/include/deal.II/base/work_stream.h
index d4c9f9f619..d06f06cbec 100644
--- a/deal.II/include/deal.II/base/work_stream.h
+++ b/deal.II/include/deal.II/base/work_stream.h
@@ -530,9 +530,10 @@ namespace WorkStream
           {
             try
               {
-                worker (current_item->work_items[i],
-                        *scratch_data,
-                        current_item->copy_datas[i]);
+                if (worker)
+                  worker (current_item->work_items[i],
+                          *scratch_data,
+                          current_item->copy_datas[i]);
               }
             catch (const std::exception &exc)
               {
@@ -633,7 +634,8 @@ namespace WorkStream
           {
             try
               {
-                copier (current_item->copy_datas[i]);
+                if (copier)
+                  copier (current_item->copy_datas[i]);
               }
             catch (const std::exception &exc)
               {
@@ -788,8 +790,13 @@ namespace WorkStream
 
     for (Iterator i=begin; i!=end; ++i)
       {
-        worker (i, scratch_data, copy_data);
-        copier (copy_data);
+        if (static_cast >(worker))
+          worker (i, scratch_data, copy_data);
+        if (static_cast >
+            (copier))
+          copier (copy_data);
       }
 #endif
   }
diff --git a/deal.II/include/deal.II/hp/dof_faces.h b/deal.II/include/deal.II/hp/dof_faces.h
index bc7129518d..2937bf902d 100644
--- a/deal.II/include/deal.II/hp/dof_faces.h
+++ b/deal.II/include/deal.II/hp/dof_faces.h
@@ -40,9 +40,6 @@ namespace internal
   {
 
     /**
-     *
-     * 
      DoFFaces
      
-     *
      * These classes are similar to the internal::hp::DoFLevel classes. We here store
      * information that is associated with faces, rather than cells, as this information is
      * independent of the hierarchical structure of cells, which are organized in levels. In 2D
@@ -66,27 +63,13 @@ namespace internal
      * @ingroup dofs
      * @author Tobias Leicht, 2006
      */
-
-
-
-
-
-
     template
-    class DoFFaces
-    {
-      /**
-       * Make the constrctor private to prevent the use
-       * of this template, only the specializations
-       * should be used
-       */
-    private:
-      DoFFaces();
-    };
+    class DoFFaces;
+
 
     /**
      * Store the indices of degrees of freedom on faces in 1D. As these would be vertices, which
-     * are treted separately, don't do anything.
+     * are treated separately, don't do anything.
      *
      * @ingroup hp
      * @author Tobias Leicht, 2006
@@ -113,6 +96,9 @@ namespace internal
     class DoFFaces<2>
     {
     public:
+      /**
+       * Indices of DoFs on the lines that bound cells.
+       */
       internal::hp::DoFObjects<1> lines;
 
       /**
@@ -134,7 +120,14 @@ namespace internal
     class DoFFaces<3>
     {
     public:
+      /**
+       * Indices of DoFs on the lines that form the edges of cells.
+       */
       internal::hp::DoFObjects<1> lines;
+
+      /**
+       * Indices of DoFs on the quads that bound cells.
+       */
       internal::hp::DoFObjects<2> quads;
 
       /**
diff --git a/deal.II/include/deal.II/hp/dof_handler.h b/deal.II/include/deal.II/hp/dof_handler.h
index 1baba27526..d4ff6b16f7 100644
--- a/deal.II/include/deal.II/hp/dof_handler.h
+++ b/deal.II/include/deal.II/hp/dof_handler.h
@@ -693,7 +693,7 @@ namespace hp
      *  the active_fe_indices in
      *  the
      *  dealii::internal::hp::DoFLevel. They
-     *  are initialized with the a
+     *  are initialized with a
      *  zero indicator, meaning
      *  that fe[0] is going to be
      *  used by default.  This
@@ -704,7 +704,6 @@ namespace hp
      *  cell has a valid
      *  active_fe_index.
      */
-
     void create_active_fe_table ();
 
     /**
@@ -791,7 +790,8 @@ namespace hp
      * levels[] tree of
      * the Triangulation objects.
      */
-    std::vector*>    levels;
+    std::vector*> levels;
+
     /**
      * Space to store the DoF
      * numbers for the faces.
@@ -837,7 +837,7 @@ namespace hp
      * actual data format used to
      * the present class.
      */
-    std::vector      vertex_dofs;
+    std::vector vertex_dofs;
 
     /**
      * For each vertex in the
@@ -866,10 +866,10 @@ namespace hp
 
     /**
      * Array to store the
-     * information, if a cell on
+     * information if a cell on
      * some level has children or
      * not. It is used by the
-     * refinement listeners as a
+     * signal slots as a
      * persistent buffer during the
      * refinement, i.e. from between
      * when pre_refinement_action is
diff --git a/deal.II/include/deal.II/hp/dof_levels.h b/deal.II/include/deal.II/hp/dof_levels.h
index 9bc69d2143..19d9f7e33d 100644
--- a/deal.II/include/deal.II/hp/dof_levels.h
+++ b/deal.II/include/deal.II/hp/dof_levels.h
@@ -29,22 +29,21 @@ namespace internal
 {
   namespace hp
   {
+    template 
+    class DoFLevel;
 
     /**
-     * Store the indices of the degrees of freedom which are located on
-     * objects of dimension @p N. We declare this general template
-     * class, but do not actually use it. Rather, only specializations of
-     * this class are used.
+     * Store the indices of the degrees of freedom that are located on
+     * objects of dimension @p structdim.
      *
-     * The things we store here is very similar to what is stored in the
+     * The things we store here are very similar to what is stored in the
      * internal::DoFHandler::DoFLevel class hierarchy (see there for more
      * information, in particular on the layout of the class hierarchy,
      * and the use of file names). There are two main
      * differences, discussed in the following subsections. In addition to
      * the data already stored by the internal::DoFHandler::DoFLevel
      * classes, we also have to store which finite element each cell
-     * uses. This is done in the DoFLevel<0> class, where for each cell we
-     * have an entry within the active_fe_indices field.
+     * uses.
      *
      *
      * 
      Offset computations
      
@@ -105,35 +104,14 @@ namespace internal
      * @ingroup hp
      * @author Wolfgang Bangerth, 1998, 2006, Oliver Kayser-Herold 2003.
      */
-    template 
+    template 
     class DoFLevel
     {
-    private:
-      /**
-       * Make the constructor private
-       * to avoid that someone uses
-       * this class.
-       */
-      DoFLevel ();
-    };
-
-
-    /**
-     * Storage for degrees of freedom on cells. See the documentation of
-     * the DoFLevel class template for more complete information on the
-     * purpose and layout of this class.
-     *
-     * @ingroup hp
-     * @author Wolfgang Bangerth, 1998, 2006, Oliver Kayser-Herold 2003.
-     */
-    template <>
-    class DoFLevel<0>
-    {
     public:
       /**
        *  Indices specifying the finite
        *  element of hp::FECollection to use
-       *  for the different cells. The
+       *  for the different cells on the current level. The
        *  meaning what a cell is, is
        *  dimension specific, therefore also
        *  the length of this vector depends
@@ -141,35 +119,16 @@ namespace internal
        *  the length of this vector equals
        *  the length of the @p lines vector,
        *  in two dimensions that of the @p
-       *  quads vector, etc.
+       *  quads vector, etc. The vector stores one element per cell
+       *  since the actiev_fe_index is unique for cells.
        */
-
       std::vector active_fe_indices;
-      /**
-       * Determine an estimate for the
-       * memory consumption (in bytes)
-       * of this object.
-       */
-      std::size_t memory_consumption () const;
-    };
-
 
-    /**
-     * Store the indices of the degrees of freedom which are located on
-     * the lines. See the general template DoFLevel for more information.
-     *
-     * @ingroup hp
-     * @author Wolfgang Bangerth, 1998, 2006, Oliver Kayser-Herold 2003.
-     */
-    template <>
-    class DoFLevel<1> : public DoFLevel<0>
-    {
-    public:
       /**
-       *  store the dof-indices and related functions of
-       *  lines
+       *  Store the dof-indices and related data of
+       *  the cells on the current level corresponding to this object.
        */
-      internal::hp::DoFObjects<1> dof_object;
+      internal::hp::DoFObjects dof_object;
 
       /**
        * Determine an estimate for the
@@ -178,60 +137,6 @@ namespace internal
        */
       std::size_t memory_consumption () const;
     };
-
-
-    /**
-     * Store the indices of the degrees of freedom which are located on
-     * quads. See the general template DoFLevel for more information.
-     *
-     * @ingroup hp
-     * @author Wolfgang Bangerth, 1998, 2006, Oliver Kayser-Herold 2003.
-     */
-    template <>
-    class DoFLevel<2> : public DoFLevel<0>
-    {
-    public:
-      /**
-       *  store the dof-indices and related functions of
-       *  quads
-       */
-      internal::hp::DoFObjects<2> dof_object;
-
-      /**
-       * Determine an estimate for the
-       * memory consumption (in bytes)
-       * of this object.
-       */
-      std::size_t memory_consumption () const;
-    };
-
-
-
-    /**
-     * Store the indices of the degrees of freedom which are located on
-     * hexhedra. See the general template DoFLevel for more information.
-     *
-     * @ingroup hp
-     * @author Wolfgang Bangerth, 1998, 2006, Oliver Kayser-Herold 2003.
-     */
-    template <>
-    class DoFLevel<3> : public DoFLevel<0>
-    {
-    public:
-      /**
-       *  store the dof-indices and related functions of
-       *  hexes
-       */
-      internal::hp::DoFObjects<3> dof_object;
-
-      /**
-       * Determine an estimate for the
-       * memory consumption (in bytes)
-       * of this object.
-       */
-      std::size_t memory_consumption () const;
-    };
-
   } // namespace hp
 
 } // namespace internal
diff --git a/deal.II/include/deal.II/hp/dof_objects.h b/deal.II/include/deal.II/hp/dof_objects.h
index 54aee6dc70..15bae95dad 100644
--- a/deal.II/include/deal.II/hp/dof_objects.h
+++ b/deal.II/include/deal.II/hp/dof_objects.h
@@ -96,7 +96,6 @@ namespace internal
      * @ingroup hp
      * @author Tobias Leicht, 2006
      */
-
     template 
     class DoFObjects
     {
@@ -104,9 +103,12 @@ namespace internal
       /**
        * Store the start index for
        * the degrees of freedom of each
-       * hex in the @p hex_dofs array.
+       * object in the @p dofs array.
+       *
+       * The type we store is then obviously the type the @p dofs array
+       * uses for indexing.
        */
-      std::vector dof_offsets;
+      std::vector::size_type> dof_offsets;
 
       /**
        * Store the global indices of
@@ -114,16 +116,15 @@ namespace internal
        * DoFLevel() for detailed
        * information.
        */
-      std::vector dofs;
+      std::vector dofs;
 
       /**
        * Set the global index of
        * the @p local_index-th
        * degree of freedom located
-       * on the hex with number @p
-       * hex_index to the value
-       * given by the last
-       * argument. The @p
+       * on the object with number @p
+       * obj_index to the value
+       * given by @p global_index. The @p
        * dof_handler argument is
        * used to access the finite
        * element that is to be used
@@ -154,8 +155,8 @@ namespace internal
        * Return the global index of
        * the @p local_index-th
        * degree of freedom located
-       * on the hex with number @p
-       * hex_index. The @p
+       * on the object with number @p
+       * obj_index. The @p
        * dof_handler argument is
        * used to access the finite
        * element that is to be used
diff --git a/deal.II/include/deal.II/lac/generic_linear_algebra.h b/deal.II/include/deal.II/lac/generic_linear_algebra.h
index 693d1b420b..765337da13 100644
--- a/deal.II/include/deal.II/lac/generic_linear_algebra.h
+++ b/deal.II/include/deal.II/lac/generic_linear_algebra.h
@@ -39,32 +39,6 @@ namespace LinearAlgebraDealII
 }
 
 
-// Dummy class. This used to check your program
-// to make sure it is compatible with all
-// linear algebra classes. In other words,
-// this is the minimal interface.
-// TODO: should we move this into tests/ only?
-namespace LinearAlgebraDummy
-{
-  class Vector
-  {
-    void compress();
-  };
-  class BlockVector
-  {
-    void compress();
-  };
-
-  class SparseMatrix
-  {
-    void compress();
-  };
-
-  class PreconditionSSOR {};
-
-}
-
-
 DEAL_II_NAMESPACE_CLOSE
 
 
@@ -89,7 +63,6 @@ namespace LinearAlgebraPETSc
   typedef PETScWrappers::BlockVector BlockVector;
 
   typedef PETScWrappers::SparseMatrix SparseMatrix;
-  typedef PETScWrappers::PreconditionSSOR PreconditionSSOR;
 
   typedef PETScWrappers::SolverCG SolverCG;
 
@@ -122,29 +95,30 @@ namespace LinearAlgebraPETSc
     typedef dealii::BlockCompressedSimpleSparsityPattern CompressedBlockSparsityPattern;
 
     /**
-     * Typedef for the AMG preconditioner type used for the
-     * top left block of the Stokes matrix.
+     * Typedef for the AMG preconditioner type.
      */
     typedef PETScWrappers::PreconditionBoomerAMG PreconditionAMG;
 
     /**
-     * Typedef for the Incomplete Cholesky preconditioner used
-     * for other blocks of the system matrix.
+     * Typedef for the Incomplete Cholesky preconditioner.
      */
     typedef PETScWrappers::PreconditionICC PreconditionIC;
 
     /**
-     * Typedef for the Incomplete LU decomposition preconditioner used
-     * for other blocks of the system matrix.
+     * Typedef for the Incomplete LU decomposition preconditioner.
      */
     typedef PETScWrappers::PreconditionILU PreconditionILU;
 
     /**
-     * Typedef for the Incomplete Jacobi decomposition preconditioner used
-     * for other blocks of the system matrix.
+     * Typedef for the Incomplete Jacobi decomposition preconditioner.
      */
     typedef PETScWrappers::PreconditionJacobi PreconditionJacobi;
 
+    /**
+     * Typedef for the SSOR preconditioner.
+     */
+    typedef PETScWrappers::PreconditionSSOR PreconditionSSOR;
+
   }
 
 }
@@ -201,29 +175,30 @@ namespace LinearAlgebraTrilinos
     typedef TrilinosWrappers::BlockSparsityPattern CompressedBlockSparsityPattern;
 
     /**
-     * Typedef for the AMG preconditioner type used for the
-     * top left block of the Stokes matrix.
+     * Typedef for the AMG preconditioner type.
      */
     typedef TrilinosWrappers::PreconditionAMG PreconditionAMG;
 
     /**
-     * Typedef for the Incomplete Cholesky preconditioner used
-     * for other blocks of the system matrix.
+     * Typedef for the Incomplete Cholesky preconditioner.
      */
     typedef TrilinosWrappers::PreconditionIC PreconditionIC;
 
     /**
-     * Typedef for the Incomplete LU decomposition preconditioner used
-     * for other blocks of the system matrix.
+     * Typedef for the Incomplete LU decomposition preconditioner.
      */
     typedef TrilinosWrappers::PreconditionILU PreconditionILU;
 
     /**
-     * Typedef for the Incomplete Jacobi decomposition preconditioner used
-     * for other blocks of the system matrix.
+     * Typedef for the Incomplete Jacobi decomposition preconditioner.
      */
     typedef TrilinosWrappers::PreconditionJacobi PreconditionJacobi;
 
+    /**
+     * Typedef for the SSOR preconditioner
+     */
+    typedef TrilinosWrappers::PreconditionSSOR PreconditionSSOR;
+
 
   }
 
diff --git a/deal.II/include/deal.II/numerics/vector_tools.templates.h b/deal.II/include/deal.II/numerics/vector_tools.templates.h
index 39b16c8cf1..0a906013a5 100644
--- a/deal.II/include/deal.II/numerics/vector_tools.templates.h
+++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h
@@ -4729,8 +4729,9 @@ namespace VectorTools
     // iterate over the list of all vector components we found and see if we
     // can find constrained ones
     unsigned int n_total_constraints_found = 0;
-    for (typename std::set,PointComparator >::
-           const_iterator it=vector_dofs.begin(); it!=vector_dofs.end(); ++it)
+    for (typename std::set,
+           PointComparator >::const_iterator it=vector_dofs.begin(); 
+           it!=vector_dofs.end(); ++it)
       {
         unsigned int n_constraints = 0;
         bool is_constrained[dim];
@@ -4769,7 +4770,7 @@ namespace VectorTools
                   normal[d] = 1.;
                 }
             AssertIndexRange(constrained_index, dim);
-            const std::vector >* constrained
+            const std::vector >* constrained
               = no_normal_flux_constraints.get_constraint_entries((*it)[constrained_index]);
             // find components to which this index is constrained to
             Assert(constrained != 0, ExcInternalError());
diff --git a/deal.II/source/base/parameter_handler.cc b/deal.II/source/base/parameter_handler.cc
index eb0bde28f4..267fef855e 100644
--- a/deal.II/source/base/parameter_handler.cc
+++ b/deal.II/source/base/parameter_handler.cc
@@ -2491,11 +2491,12 @@ MultipleParameterLoop::~MultipleParameterLoop ()
 
 
 
-bool MultipleParameterLoop::read_input (std::istream &input)
+bool MultipleParameterLoop::read_input (std::istream &input,
+                                        const std::string &filename)
 {
   AssertThrow (input, ExcIO());
 
-  bool x = ParameterHandler::read_input (input);
+  bool x = ParameterHandler::read_input (input, filename);
   if (x)
     init_branches ();
   return x;
diff --git a/deal.II/source/grid/tria.cc b/deal.II/source/grid/tria.cc
index 1aa9eda828..80c9b8420a 100644
--- a/deal.II/source/grid/tria.cc
+++ b/deal.II/source/grid/tria.cc
@@ -11327,7 +11327,8 @@ template 
 typename Triangulation::active_cell_iterator
 Triangulation::end_active (const unsigned int level) const
 {
-  return (level == levels.size()-1 ?
+  Assert (level= levels.size()-1 ?
           active_cell_iterator(end()) :
           begin_active (level+1));
 }
diff --git a/deal.II/source/hp/dof_handler.cc b/deal.II/source/hp/dof_handler.cc
index 0b4b865540..205071413e 100644
--- a/deal.II/source/hp/dof_handler.cc
+++ b/deal.II/source/hp/dof_handler.cc
@@ -577,7 +577,7 @@ namespace internal
           for (unsigned int level=0; leveln_levels(); ++level)
             {
               dof_handler.levels[level]->dof_object.dof_offsets
-                = std::vector (
+                = std::vector::size_type> (
                     dof_handler.tria->n_raw_lines(level),
                     DoFHandler::invalid_dof_index);
 
@@ -694,7 +694,7 @@ namespace internal
           for (unsigned int level=0; leveln_levels(); ++level)
             {
               dof_handler.levels[level]->dof_object.dof_offsets
-                = std::vector (
+                = std::vector::size_type> (
                     dof_handler.tria->n_raw_quads(level),
                     DoFHandler::invalid_dof_index);
 
@@ -846,7 +846,8 @@ namespace internal
             // active ones will have a
             // non-invalid value later on
             dof_handler.faces->lines.dof_offsets
-              = std::vector (dof_handler.tria->n_raw_lines(),
+              = std::vector::size_type>
+                                                     (dof_handler.tria->n_raw_lines(),
                                                       DoFHandler::invalid_dof_index);
             dof_handler.faces->lines.dofs
               = std::vector (n_line_slots,
@@ -1059,7 +1060,8 @@ namespace internal
           for (unsigned int level=0; leveln_levels(); ++level)
             {
               dof_handler.levels[level]->dof_object.dof_offsets
-                = std::vector (dof_handler.tria->n_raw_hexs(level),
+                = std::vector::size_type>
+                                                       (dof_handler.tria->n_raw_hexs(level),
                                                         DoFHandler::invalid_dof_index);
 
               types::global_dof_index next_free_dof = 0;
@@ -1214,7 +1216,8 @@ namespace internal
             if (true)
               {
                 dof_handler.faces->quads.dof_offsets
-                  = std::vector (dof_handler.tria->n_raw_quads(),
+                  = std::vector::size_type>
+                                                         (dof_handler.tria->n_raw_quads(),
                                                           DoFHandler::invalid_dof_index);
                 dof_handler.faces->quads.dofs
                   = std::vector (n_quad_slots,
@@ -3119,7 +3122,7 @@ namespace hp
   {
     create_active_fe_table ();
 
-    // Remember if the cells have already
+    // Remember if the cells already have
     // children. That will make the transfer
     // of the active_fe_index to the finer
     // levels easier.
@@ -3130,101 +3133,20 @@ namespace hp
         std::vector *has_children_level =
           new std::vector (cells_on_level);
 
-        // Check for each cell, if it has children.
-        std::transform (tria->levels[i]->cells.refinement_cases.begin (),
-                        tria->levels[i]->cells.refinement_cases.end (),
-                        has_children_level->begin (),
-                        std::bind2nd (std::not_equal_to(),
-                                      static_cast(RefinementCase::no_refinement)));
-
-        has_children.push_back (has_children_level);
-      }
-  }
-
-
-
-  template <>
-  void DoFHandler<1>::pre_refinement_action ()
-  {
-    create_active_fe_table ();
-
-    // Remember if the cells have already
-    // children. That will make the transfer
-    // of the active_fe_index to the finer
-    // levels easier.
-    Assert (has_children.size () == 0, ExcInternalError ());
-    for (unsigned int i=0; in_raw_cells (i);
-        std::vector *has_children_level =
-          new std::vector (cells_on_level);
-
-        // Check for each cell, if it has
-        // children. here we cannot use
-        // refinement_cases, since it is unused in
-        // 1d (as there is only one choice
-        // anyway). use the 'children' vector
-        // instead
-        std::transform (tria->levels[i]->cells.children.begin (),
-                        tria->levels[i]->cells.children.end (),
-                        has_children_level->begin (),
-                        std::bind2nd (std::not_equal_to(), -1));
-
-        has_children.push_back (has_children_level);
-      }
-  }
-
-
-
-  template <>
-  void DoFHandler<1,2>::pre_refinement_action ()
-  {
-    create_active_fe_table ();
-
-    // Remember if the cells have already
-    // children. That will make the transfer
-    // of the active_fe_index to the finer
-    // levels easier.
-    Assert (has_children.size () == 0, ExcInternalError ());
-    for (unsigned int i=0; in_raw_lines(i);
-        std::vector *has_children_level =
-          new std::vector (lines_on_level);
-
-        // Check for each cell, if it has children.
-        std::transform (tria->levels[i]->cells.children.begin (),
-                        tria->levels[i]->cells.children.end (),
-                        has_children_level->begin (),
-                        std::bind2nd (std::not_equal_to(), -1));
-
-        has_children.push_back (has_children_level);
-      }
-  }
-
-
-
-  template <>
-  void DoFHandler<1,3>::pre_refinement_action ()
-  {
-    create_active_fe_table ();
-
-    // Remember if the cells have already
-    // children. That will make the transfer
-    // of the active_fe_index to the finer
-    // levels easier.
-    Assert (has_children.size () == 0, ExcInternalError ());
-    for (unsigned int i=0; in_raw_lines(i);
-        std::vector *has_children_level =
-          new std::vector (lines_on_level);
-
-        // Check for each cell, if it has children.
-        std::transform (tria->levels[i]->cells.children.begin (),
-                        tria->levels[i]->cells.children.end (),
-                        has_children_level->begin (),
-                        std::bind2nd (std::not_equal_to(), -1));
+        // Check for each cell, if it has children. in 1d,
+        // we don't store refinement cases, so use the 'children'
+        // vector instead
+        if (dim == 1)
+          std::transform (tria->levels[i]->cells.children.begin (),
+                          tria->levels[i]->cells.children.end (),
+                          has_children_level->begin (),
+                          std::bind2nd (std::not_equal_to(), -1));
+        else
+          std::transform (tria->levels[i]->cells.refinement_cases.begin (),
+                          tria->levels[i]->cells.refinement_cases.end (),
+                          has_children_level->begin (),
+                          std::bind2nd (std::not_equal_to(),
+                                        static_cast(RefinementCase::no_refinement)));
 
         has_children.push_back (has_children_level);
       }
@@ -3309,11 +3231,6 @@ namespace hp
          children_level != has_children.end ();
          ++children_level)
       delete (*children_level);
-    /*
-      for_each (has_children.begin (),
-      has_children.end (),
-      delete());
-    */
     has_children.clear ();
   }
 
diff --git a/deal.II/source/hp/dof_levels.cc b/deal.II/source/hp/dof_levels.cc
index 6ee1b772fd..ff5dc868b1 100644
--- a/deal.II/source/hp/dof_levels.cc
+++ b/deal.II/source/hp/dof_levels.cc
@@ -23,40 +23,19 @@ namespace internal
 {
   namespace hp
   {
-
-    std::size_t
-    DoFLevel<1>::memory_consumption () const
-    {
-      return (DoFLevel<0>::memory_consumption() +
-              MemoryConsumption::memory_consumption (dof_object));
-    }
-
-
-
-    std::size_t
-    DoFLevel<2>::memory_consumption () const
-    {
-      return (DoFLevel<0>::memory_consumption () +
-              MemoryConsumption::memory_consumption (dof_object));
-    }
-
-
-
+    template 
     std::size_t
-    DoFLevel<3>::memory_consumption () const
+    DoFLevel::memory_consumption () const
     {
-      return (DoFLevel<0>::memory_consumption () +
-              MemoryConsumption::memory_consumption (dof_object));
+      return MemoryConsumption::memory_consumption (active_fe_indices) +
+             MemoryConsumption::memory_consumption (dof_object);
     }
 
 
-
-    std::size_t
-    DoFLevel<0>::memory_consumption () const
-    {
-      return MemoryConsumption::memory_consumption (active_fe_indices);
-    }
-
+    // explicit instantiations
+    template std::size_t DoFLevel<1>::memory_consumption () const;
+    template std::size_t DoFLevel<2>::memory_consumption () const;
+    template std::size_t DoFLevel<3>::memory_consumption () const;
   }
 }
 
diff --git a/deal.II/source/numerics/vector_tools_interpolate.inst.in b/deal.II/source/numerics/vector_tools_interpolate.inst.in
index 5eba56a30f..f877d0f4f0 100644
--- a/deal.II/source/numerics/vector_tools_interpolate.inst.in
+++ b/deal.II/source/numerics/vector_tools_interpolate.inst.in
@@ -70,81 +70,6 @@ for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimens
   }
 
 
-
-for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; DH : DOFHANDLER_TEMPLATES)
-  {
-#if deal_II_dimension <= deal_II_space_dimension
-    namespace VectorTools \{
-      template
-        void interpolate_boundary_values
-        (const Mapping    &,
-         const DH &,
-         const FunctionMap::type &,
-         std::map       &,
-         const ComponentMask    &);
-
-      template
-        void interpolate_boundary_values
-        (const Mapping    &,
-         const DH &,
-         const types::boundary_id,
-         const Function   &,
-         std::map       &,
-         const ComponentMask    &);
-
-      template
-        void interpolate_boundary_values (
-          const DH &,
-          const types::boundary_id,
-          const Function   &,
-          std::map       &,
-          const ComponentMask    &);
-
-      template
-        void interpolate_boundary_values
-        (const DH &,
-         const FunctionMap::type &,
-         std::map       &,
-         const ComponentMask    &);
-
-      template
-        void interpolate_boundary_values (
-          const Mapping    &,
-          const DH &,
-          const FunctionMap::type   &,
-          ConstraintMatrix                    &,
-          const ComponentMask    &);
-
-      template
-        void interpolate_boundary_values
-        (const Mapping    &,
-         const DH &,
-         const types::boundary_id,
-         const Function   &,
-         ConstraintMatrix                    &,
-         const ComponentMask             &);
-
-      template
-        void interpolate_boundary_values (
-          const DH &,
-          const types::boundary_id,
-          const Function   &,
-          ConstraintMatrix                    &,
-          const ComponentMask    &);
-
-      template
-        void interpolate_boundary_values (
-          const DH &,
-          const FunctionMap::type   &,
-          ConstraintMatrix                    &,
-          const ComponentMask    &);
-
-      \}
-#endif
-  }
-
-
-
 for (deal_II_dimension : DIMENSIONS)
 {
   namespace VectorTools \{