diff --git a/AMDiS/CMakeLists.txt b/AMDiS/CMakeLists.txt
index 10386a6e36944d20dd5a2ab38d753d8ee7a0f787..01afc914418339442c0ff284e963d782f993d2d1 100644
--- a/AMDiS/CMakeLists.txt
+++ b/AMDiS/CMakeLists.txt
@@ -25,8 +25,6 @@ option(USE_PETSC_DEV false)
 set(ENABLE_PARMETIS false)
 option(ENABLE_ZOLTAN false)
 option(ENABLE_UMFPACK "use umfpack solver" false)
-option(ENABLE_MKL "use the mkl" false)
-SET(MKL_DIR "" CACHE PATH "MKL directory")
 option(ENABLE_TESTS "compile the tests" false)
 
 if(ENABLE_INTEL)
@@ -259,11 +257,6 @@ if(ENABLE_UMFPACK)
 	LIST(APPEND AMDiS_LIBS amdis blas amd umfpack)
 endif(ENABLE_UMFPACK)
 
-if(ENABLE_MKL)
-	include_directories(${MKL_DIR})
-	SET(COMPILEFLAGS "${COMPILEFLAGS} -DHAVE_MKL=1")
-endif(ENABLE_MKL)
-
 SET(COMPOSITE_SOURCE_DIR ${SOURCE_DIR}/compositeFEM)
 SET(COMPOSITE_FEM_SRC ${COMPOSITE_SOURCE_DIR}/CFE_Integration.cc 
 		      ${COMPOSITE_SOURCE_DIR}/CFE_NormAndErrorFcts.cc 
diff --git a/demo/src/ellipt.cc b/demo/src/ellipt.cc
index c737ba5f8520bb7a97e18e964abcee91845709a6..098d44c0e4d27b10ca00b0d3ed989553279aabb1 100644
--- a/demo/src/ellipt.cc
+++ b/demo/src/ellipt.cc
@@ -45,7 +45,7 @@ int main(int argc, char* argv[])
   FUNCNAME("main");
 
   // ===== check for init file =====
-  TEST_EXIT(argc == 2)("usage: ellipt initfile\n");
+  TEST_EXIT(argc >= 2)("usage: ellipt initfile\n");
 
 
   // ===== init parameters =====
diff --git a/doc/tutorial/composite.tex b/doc/tutorial/composite.tex
deleted file mode 100644
index d23cbc57edb5ff4fc599efca8a025e8e4a9d0c5d..0000000000000000000000000000000000000000
--- a/doc/tutorial/composite.tex
+++ /dev/null
@@ -1,6 +0,0 @@
-\section{Composite finite elements}
-
-\subsection{Source code}
-\subsection{Parameter file}
-\subsection{Macro file}
-\subsection{Output}
\ No newline at end of file
diff --git a/doc/tutorial/couple.tex b/doc/tutorial/couple.tex
deleted file mode 100644
index 0aedcfc6de8b6ce6e51cc1c91513613d5a2e0a7a..0000000000000000000000000000000000000000
--- a/doc/tutorial/couple.tex
+++ /dev/null
@@ -1,314 +0,0 @@
-\section{Coupled problems}
-\label{s:couple}
-In this example, we solve the same problem as in Section \ref{s:vecellipt}, but here we treat the two equations as two coupled problems. The main difference is that the equations now aren't assembled into the same large system of equations, but into two separated systems of equations, that have to be solved separately. We define the two problems
-\begin{eqnarray}
--\Delta u &=& f ~~~\mbox{in } \Omega \subset \mathbb{R}^{dim}\\
-u &=& g ~~~ \mbox{on } \partial\Omega
-\end{eqnarray}
-and
-\begin{equation}
-v = u.
-\end{equation}
-We first solve the first problem and then use its solution to solve the second problem. This happens in every iteration of the adaptation loop. Both problems should use the same mesh. Mesh adaptation is done by the first problem. So one iteration now looks like illustrated in Figure \ref{f:coupled iteration}.  
-\begin{figure}
-\center
-\includegraphics[width=0.9\textwidth]{fig/coupled_iteration.pdf}
-\caption{State diagram of the coupled iteration.}
-\label{f:coupled iteration}
-\end{figure} 
-
-\begin{table}
-\center
-\begin{tabular}{|cc|c|c|c|}
-\hline 
-& & {\bf 1d} & {\bf 2d} & {\bf 3d} \\
-\hline 
-{\bf source code} & \tt src/ & \multicolumn{3}{|c|}{\tt couple.cc} \\
-\hline 
-{\bf parameter file} & \tt init/ & \tt couple.dat.1d & \tt couple.dat.2d & \tt couple.dat.3d \\
-\hline 
-{\bf macro file} & \tt macro/ & \tt macro.stand.1d & \tt macro.stand.2d & \tt macro.stand.3d  \\
-\hline 
-{\bf output files} & \tt output/ & \multicolumn{3}{|c|}{\tt couple.mesh, couple.dat} \\
-\hline 
-\end{tabular}
-\caption{Files of the {\tt couple} example.}
-\end{table}
-
-\subsection{Source code}
-In the previous examples, the iteration was implemented always by the corresponding problem class. In this example, we want to couple two problems within one iteration, so the default implementation can't be used. For that reason, we define our own coupled iteration class\\ \verb+MyCoupledIteration+ which implements the \verb+ProblemIterationInterface+:
-
-\begin{lstlisting}{}
-class MyCoupledIteration : public ProblemIterationInterface
-{
-public:
-  MyCoupledIteration(ProblemStatBase *prob1,
-		     ProblemStatBase *prob2)
-    : problem1(prob1),
-      problem2(prob2)
-  {}
-\end{lstlisting}
-
-In the constructor pointers to the two problems are assigned to the private members \verb+problem1+ and \verb+problem2+. Note that the pointers point to the interface \verb+ProblemStatBase+ and not to\\ \verb+ProblemScal+. This leads to a more general implementation. If e.g. two vector valued problems should be coupled in the future, we could use our iteration class without modifications. 
-
-Now, we implement the needed interface methods:  
-
-\begin{lstlisting}{}
-  void beginIteration(AdaptInfo *adaptInfo) 
-  {
-    FUNCNAME("StandardProblemIteration::beginIteration()");
-    MSG("\n");
-    MSG("begin of iteration %d\n", adaptInfo->getSpaceIteration()+1);
-    MSG("=============================\n");
-  }
-
-  void endIteration(AdaptInfo *adaptInfo) {
-    FUNCNAME("StandardProblemIteration::endIteration()");
-    MSG("\n");
-    MSG("end of iteration %d\n", adaptInfo->getSpaceIteration()+1);
-    MSG("=============================\n");
-  }
-\end{lstlisting}
-
-These two functions are called at the beginning and at the end of each iteration. Here, we only prompt some output. 
-
-The method \verb+oneIteration+ is the crucial part of our implementation:
-
-\begin{lstlisting}{}
-  Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION) 
-  {
-    Flag flag, markFlag;
-    if (toDo.isSet(MARK)) markFlag = problem1->markElements(adaptInfo);
-    if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_REFINED))
-      flag = problem1->refineMesh(adaptInfo);
-    
-    if (toDo.isSet(BUILD)) problem1->buildAfterCoarsen(adaptInfo, markFlag);
-    if (toDo.isSet(SOLVE)) problem1->solve(adaptInfo);
- 
-    if (toDo.isSet(BUILD)) problem2->buildAfterCoarsen(adaptInfo, markFlag);
-    if (toDo.isSet(SOLVE)) problem2->solve(adaptInfo);
-
-    if (toDo.isSet(ESTIMATE)) problem1->estimate(adaptInfo);    
-    return flag;
-  }
-\end{lstlisting}
-
-The \verb+toDo+ flag is used by the adaptation loop to determine which parts of the iteration should be performed. The first iteration is always an iteration without mesh adaptation (see Figure \ref{f:stationary loop}). So we start our iteration by marking and adapting the mesh. The mesh and its adaptation is managed by the first problem. So we call \verb+markElements+ and \verb+refineMesh+ of \verb+problem1+. Note that no mesh coarsenings have to be performed in our example. Afterwards, \verb+problem1+ assembles its system of equations by \verb+buildAfterCoarsen+. Assemblage and mesh adaptation are nested operations in AMDiS (\verb+buildBeforeRefine+, \verb+refineMesh+, \verb+buildBeforeCoarsen+,\\\verb+coarsenMesh+,\verb+buildAfterCoarsen+). Here, we implement a simplified version.
-
-After \verb+problem1+ has solved its system of equations, \verb+problem2+ can assemble and solve its equations system using the solution of the first problem as right hand side. In the method \verb+oneIteration+, only the order of method calls is determined. The dependency to the solution of the first problem is created later when the operator for the right hand side of \verb+problem2+ is created. 
-
-After also the second problem computed its solution, \verb+problem1+ does the error estimation (remember: mesh adaptation is managed by \verb+problem1+).
-\begin{figure}
-\center
-\includegraphics[width=0.9\textwidth]{fig/stationary_loop.pdf}
-\caption{Stationary adaptation loop.}
-\label{f:stationary loop}
-\end{figure} 
-
-Now, the access to the coupled problems is implemented and the member variables are defined:
-
-\begin{lstlisting}{}
-  int getNumProblems() 
-  {
-    return 2;
-  }
-
-  ProblemStatBase *getProblem(int number = 0) 
-  {
-    FUNCNAME("CoupledIteration::getProblem()");
-    if (number == 0) return problem1;
-    if (number == 1) return problem2;
-    ERROR_EXIT("invalid problem number\n");
-    return NULL;
-  }
-
-private:
-  ProblemStatBase *problem1;
-  ProblemStatBase *problem2;
-};
-\end{lstlisting}
-
-The class \verb+MyCoupledIteration+ is finished now.
-
-The next class, \verb+Identity+, implements the identity $I(x)=x$ for \verb+double+ variables. An arbitrary degree can be given to it. The class is used later to determine the quadrature degree used for the right hand side of \verb+problem2+. 
-
-\begin{lstlisting}{}
-class Identity : public AbstractFunction<double, double>
-{
-public:
-  Identity(int degree) : AbstractFunction<double, double>(degree) {}
-
-  double operator()(const double& x) const 
-  {
-    return x;
-  }
-};
-\end{lstlisting}
-
-Now, we start with the main program:
-
-\begin{lstlisting}{}
-int main(int argc, char* argv[])
-{
-  FUNCNAME("main");
-  TEST_EXIT(argc == 2)("usage: couple initfile\n");
-  Parameters::init(true, argv[1]);
-
-  // ===== create and init the first problem ===== 
-  ProblemScal problem1("problem1");
-  problem1.initialize(INIT_ALL);
-
-  // ===== add boundary conditions for problem1 =====
-  problem1.addDirichletBC(1, new G);
-\end{lstlisting}
-
-So far, we created and initialized \verb+problem1+ and its boundary conditions.
-
-Now, we create \verb+problem2+. It should have its own finite element space, system, solver and file writer, but the mesh should be adopted from \verb+problem1+.  
-
-\begin{lstlisting}{}
-  // ===== create and init the second problem ===== 
-  Flag initFlag = 
-    INIT_FE_SPACE |
-    INIT_SYSTEM |
-    INIT_SOLVER |
-    INIT_FILEWRITER;
-
-  Flag adoptFlag =
-    CREATE_MESH |
-    INIT_MESH;
-
-  ProblemScal problem2("problem2");
-  problem2.initialize(initFlag,
-		      &problem1,
-		      adoptFlag);
-\end{lstlisting}
-
-The operators for the first problem are defined like in Section
-\ref{s:ellipt code}. Here, we only show the operators of
-\verb+problem2+.
-
-\begin{lstlisting}{}
-  // ===== create operators for problem2 =====
-  Operator matrixOperator2(problem2.getFeSpace());
-  matrixOperator2.addZeroOrderTerm(new Simple_ZOT);
-  problem2.addMatrixOperator(&matrixOperator2);
-  
-  Operator rhsOperator2(problem2.getFeSpace());
-  rhsOperator2.addZeroOrderTerm(new VecAtQP_ZOT(problem1.getSolution(),
-                                                new Identity(degree)));
-  problem2.addVectorOperator(&rhsOperator2);
-\end{lstlisting}
-
-At the left hand side, we have just an ordinary \verb+Simple_ZOT+. At
-the right hand side, we have a zero order term of the form $f(u)$ with
-$f=I$ the identity. $u$ is given by the solution DOF vector of
-\verb+problem1+. $I$ maps the values of the DOF vector evaluated at
-quadrature points to itself. The function degree is used to determine
-the needed quadrature degree in the assembler.
- 
-Now, the adaptation loop is created:
-
-\begin{lstlisting}{}
-  // ===== create adaptation loop and iteration interface =====
-  AdaptInfo *adaptInfo = new AdaptInfo("couple->adapt", 1);
-
-  MyCoupledIteration coupledIteration(&problem1, &problem2);
-
-  AdaptStationary *adapt = new AdaptStationary("couple->adapt",
-                                               &coupledIteration,
-                                               adaptInfo);
-\end{lstlisting}
-
-Note that not a pointer to one of the problems is passed to the
-adaptation loop, but a pointer to the \verb+coupledIteration+ object,
-which in turn knows both problems.
-
-The adaptation loop is now started. After it is finished, the
-solutions of both problems are written.
-
-\begin{lstlisting}{}
-  // ===== start adaptation loop =====
-  adapt->adapt();
-
-  // ===== write solution ===== 
-  problem1.writeFiles(adaptInfo, true);
-  problem2.writeFiles(adaptInfo, true);
-}
-\end{lstlisting}
-
-\subsection{Parameter file}
-
-We have one adaptation loop called \verb+couple->adapt+:
-
-\begin{lstlisting}{}
-couple->adapt->tolerance:         1e-8
-couple->adapt->max iteration:     10
-couple->adapt->refine bisections: 2
-\end{lstlisting}
-
-The coupled problem consits of two sub problems. The first problem
-creates the mesh, solves its linear system of equations, estimates the
-error, adapts the mesh, and finally writes its output:
-
-\begin{lstlisting}{}
-coupleMesh->macro file name:    ./macro/macro.stand.2d
-coupleMesh->global refinements: 0
-
-problem1->mesh:                 coupleMesh
-problem1->dim:                  2
-problem1->polynomial degree:    1
-
-problem1->solver:                 cg % no, bicgstab, cg, gmres, odir, ores
-problem1->solver->max iteration:  1000
-problem1->solver->tolerance:      1.e-8
-problem1->solver->left precon:    diag
-
-problem1->estimator:              residual
-problem1->estimator->C0:          0.1 % constant of element residual
-problem1->estimator->C1:          0.1 % constant of jump residual
-
-problem1->marker->strategy:       2 % 0: no 1: GR 2: MS 3: ES 4:GERS
-problem1->marker->MSGamma:        0.5
-
-problem1->output->filename:       output/problem1
-problem1->output->AMDiS format:   1
-problem1->output->AMDiS mesh ext: .mesh
-problem1->output->AMDiS data ext: .dat
-\end{lstlisting}
-
-The second problem uses the mesh of \verb+problem1+. So it creates no
-mesh, no estimator, and no marker. But a solver is needed to solve
-\verb+problem2+s linear system of equations, and a file writer to
-write the solution:
-
-\begin{lstlisting}{}
-problem2->dim:                  2
-problem2->polynomial degree:    1
-
-problem2->solver:                 cg % no, bicgstab, cg, gmres, odir, ores
-problem2->solver->max iteration:  1000
-problem2->solver->tolerance:      1.e-8
-problem2->solver->left precon:    diag
-
-problem2->estimator:              no
-
-problem2->marker->strategy:       0
-
-problem2->output->filename:       output/problem2
-problem2->output->AMDiS format:   1
-problem2->output->AMDiS mesh ext: .mesh
-problem2->output->AMDiS data ext: .dat
-\end{lstlisting}
-
-\subsection{Macro file}
-We again use the macro file \verb+macro/macro.stand.2d+, which was
-described in Section \ref{s:ellipt macro}.
-
-\subsection{Output}
-The solution of the first problem is written to the files
-\verb+output/problem1.mesh+ and\\ \verb+output/problem1.dat+.  The
-solution of the second problem is written to the files\\
-\verb+output/problem2.mesh+ and \verb+output/problem2.dat+.  We don't
-visualize the results here, because they conform with the results
-showed in Section \ref{s:vecellipt output}.
diff --git a/doc/tutorial/ellipt.tex b/doc/tutorial/ellipt.tex
index 218ac096d79a56760ddc92d09d5b0f7d149d37d5..79d7c01277aeac7a89a092200af8a9db351920eb 100644
--- a/doc/tutorial/ellipt.tex
+++ b/doc/tutorial/ellipt.tex
@@ -123,7 +123,7 @@ int main(int argc, char* argv[])
   FUNCNAME("main");
 
   // ===== check for init file =====
-  TEST_EXIT(argc == 2)("usage: ellipt initfile\n");
+  TEST_EXIT(argc >= 2)("usage: ellipt initfile\n");
 
   // ===== init parameters =====
   Parameters::init(true, argv[1]);
@@ -141,11 +141,11 @@ argument specifies, whether the initialized parameters should be
 printed after initialization for debug reasons. The second argument is
 the name of the parameter file.
 
-Now, a scalar problem with name \verb+ellipt+ is created and initialized:  
+Now, a stationary problem with name \verb+ellipt+ is created and initialized:  
 
 \begin{lstlisting}{}
   // ===== create and init the scalar problem ===== 
-  ProblemScal ellipt("ellipt");
+  ProblemStat ellipt("ellipt");
   ellipt.initialize(INIT_ALL);
 \end{lstlisting}
 
@@ -154,9 +154,9 @@ parameter file that belong to this problem. In this case, all
 parameters with prefix \verb+ellipt->+ are associated to this
 problem. The initialization argument \verb+INIT_ALL+ means that all
 problem modules are created in a standard way. Those are: The finite
-element space including the corresponding mesh, needed system matrices
-and vectors, an iterative solver, an estimator, a marker, and a file
-writer for the computed solution. The initialization of these
+element space including the corresponding mesh, the required system
+matrices and vectors, an iterative solver, an estimator, a marker, and
+a file writer for the computed solution. The initialization of these
 components can be controlled through the parameter file (see Section
 \ref{s:ellipt parameter}).
 
@@ -182,47 +182,49 @@ used as parameter prefix, the problem as implementation of an
 iteration interface, and the \verb+AdaptInfo+ object. The adaptation
 loop only knows when to perform which part of an iteration. The
 implementation and execution of the single steps is delegated to an
-iteration interface, here implemented by the scalar problem
+iteration interface, here implemented by the stationary problem
 \verb+ellipt+.
 
-Now, we define boundary conditions:
-
-\begin{lstlisting}{}
-  // ===== add boundary conditions =====
-  ellipt.addDirichletBC(1, new G);
-\end{lstlisting}
-
-We have one Dirichlet boundary condition associated with identifier
-$1$. All nodes belonging to this boundary are set to the value of
-function \verb+G+ at the corresponding coordinates. In the macro file
-(see Section \ref{s:ellipt macro}) the Dirichlet boundary is marked
-with identifier $1$, too. So the nodes can be uniquely determined.
-
 The operators now are defined as follows:
-
 \begin{lstlisting}{}
   // ===== create matrix operator =====
   Operator matrixOperator(ellipt.getFeSpace());
   matrixOperator.addSecondOrderTerm(new Laplace_SOT);
-  ellipt.addMatrixOperator(&matrixOperator);
+  ellipt.addMatrixOperator(matrixOperator, 0, 0);
 
   // ===== create rhs operator =====
   int degree = ellipt.getFeSpace()->getBasisFcts()->getDegree();
   Operator rhsOperator(ellipt.getFeSpace());
   rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
-  ellipt.addVectorOperator(&rhsOperator);
+  ellipt.addVectorOperator(rhsOperator, 0);
 \end{lstlisting}
-
-First, we define a matrix operator (left hand side operator) on the
-finite element space of the problem. Now, we add the term $-\Delta u$
-to it. Note that the minus sign isn't explicitly given, but implicitly
+We define a matrix operator (left hand side operator) on the finite
+element space of the problem. The term $-\Delta u$ is added to
+it. Note that the minus sign isn't explicitly given, but implicitly
 contained in \verb+Laplace_SOT+. With \verb+addMatrixOperator+ we add
-the operator to the problem. The definition of the right hand side is
-done in a similar way. We choose the degree of our function to be
-equal to the current basis function degree.
+the operator to the stationary problem definition. The both zeros
+represent the position of the operator in the operator matrix. As we
+are about to define a scalar equation, there is only the 0/0 position
+in the operator matrix. The definition of the right hand side is done
+in a similar way. We choose the degree of our function to be equal to
+the current basis function degree.
 
-Finally we start the adaptation loop and afterwards write out the results:
+Now, we define boundary conditions:
+\begin{lstlisting}{}
+  // ===== add boundary conditions =====
+  ellipt.addDirichletBC(1, 0, 0, new G);
+\end{lstlisting}
+We have one Dirichlet boundary condition associated with identifier
+$1$. All nodes belonging to this boundary are set to the value of
+function \verb+G+ at the corresponding coordinates. In the macro file
+(see Section \ref{s:ellipt macro}) the Dirichlet boundary is marked
+with identifier $1$, too. So the nodes can be uniquely determined. As
+with adding operators to the operator matrix, we have to define the
+operator, on which the boundary condition will be applied. Thus we
+have to provide the matrix position indices after the boundary
+identifier.
 
+Finally we start the adaptation loop and afterwards write out the results:
 \begin{lstlisting}{}
   // ===== start adaption loop =====
   adapt.adapt();
@@ -269,12 +271,15 @@ a finer mesh.
 ellipt->mesh:                      elliptMesh
 ellipt->dim:                       2
 ellipt->polynomial degree:         1
+ellipt->components:                1
 \end{lstlisting}
 
 Now, we construct the finite element space for the problem
 \verb+ellipt+ (see Section \ref{s:ellipt code}). We use the mesh
 \verb+elliptMesh+, set the problem dimension to 2, and choose Lagrange
-basis functions of degree 1.
+basis functions of degree 1. The number of components, i.e.\
+variables, in the equation is set to 1, since we are about to define a
+scalar PDE.
 
 \begin{lstlisting}{}
 ellipt->solver:                    cg   
@@ -290,10 +295,10 @@ tolerance of $10^{-8}$ is reached. Furthermore, we apply diagonal left
 preconditioning, and no right preconditioning.
 
 \begin{lstlisting}{}
-ellipt->estimator:                 residual
-ellipt->estimator->error norm:     1
-ellipt->estimator->C0:             0.1
-ellipt->estimator->C1:             0.1
+ellipt->estimator[0]:                 residual
+ellipt->estimator[0]->error norm:     1
+ellipt->estimator[0]->C0:             0.1
+ellipt->estimator[0]->C1:             0.1
 \end{lstlisting}
 
 As error estimator we use the residual method. The used error norm is
@@ -301,17 +306,18 @@ the H1-norm (instead of the L2-norm: 2). Element residuals (C0) and
 jump residuals (C1) both are weighted by factor $0.1$.
 
 \begin{lstlisting}{}
-ellipt->marker->strategy:          2     % 0: no 1: GR 2: MS 3: ES 4:GERS
-ellipt->marker->MSGamma:           0.5
+ellipt->marker[0]->strategy:          2     % 0: no 1: GR 2: MS 3: ES 4:GERS
+ellipt->marker[0]->MSGamma:           0.5
 \end{lstlisting}
 
 After error estimation, elements are marked for refinement and
 coarsening. Here, we use the maximum strategy with $\gamma = 0.5$.
 
 \begin{lstlisting}{}
-ellipt->adapt->tolerance:          1e-4
+ellipt->adapt[0]->tolerance:          1e-4
+ellipt->adapt[0]->refine bisections:  2
+
 ellipt->adapt->max iteration:      10
-ellipt->adapt->refine bisections:  2
 \end{lstlisting}
 
 The adaptation loop stops, when an error tolerance of $10^{-4}$ is
diff --git a/doc/tutorial/examples.tex b/doc/tutorial/examples.tex
index 679bb75a8523ab9b6a6fb278a5de6074cbd69aaa..d1351822e3770ae5676072140f34d2d839735a6b 100644
--- a/doc/tutorial/examples.tex
+++ b/doc/tutorial/examples.tex
@@ -4,7 +4,6 @@
 \input{ellipt.tex}
 \input{heat.tex}
 \input{vecellipt.tex}
-\input{couple.tex}
 \input{neumann.tex}
 \input{periodic.tex}
 \input{projections.tex}
diff --git a/doc/tutorial/heat.tex b/doc/tutorial/heat.tex
index b7ef5118642d4db8a8a7a27855c8f7250b962c09..11253f1d874c8296266f8fe7c6f67352c9eefc54 100644
--- a/doc/tutorial/heat.tex
+++ b/doc/tutorial/heat.tex
@@ -95,19 +95,18 @@ the same way.
 Now, we begin with the implementation of class \verb+Heat+, that
 represents the instationary problem. In Figure \ref{f:heat_uml}, its
 class diagram is shown. \verb+Heat+ is derived from class
-\verb+ProblemInstatScal+ which leads to following properties:
+\verb+ProblemInstat+ which leads to following properties:
 \begin{description}
 \item \verb+Heat+ implements the \verb+ProblemTimeInterface+, so the
   adaptation loop can set the current time and schedule timesteps.
 \item \verb+Heat+ implements \verb+ProblemStatBase+ in the role as
   initial (stationary) problem. The adaptation loop can compute the
   initial solution through this interface. The single iteration steps
-  can be overloaded by sub classes of
-  \verb+ProblemInstatScal+. Actually, the initial solution is computed
-  through the method \verb+solveInitialProblem+ of
-  \verb+ProblemTimeInterface+. But this method is implemented by
-  \verb+ProblemInstatScal+ interpreting itself as initial stationary
-  problem.
+  can be overloaded by sub classes of \verb+ProblemInstat+. Actually,
+  the initial solution is computed through the method
+  \verb+solveInitialProblem+ of \verb+ProblemTimeInterface+. But this
+  method is implemented by \verb+ProblemInstat+ interpreting
+  itself as initial stationary problem.
 \item \verb+Heat+ knows another implementation of
   \verb+ProblemStatBase+: This other implementation represents a
   stationary problem which is solved within each
@@ -116,14 +115,14 @@ class diagram is shown. \verb+Heat+ is derived from class
 The first lines of class \verb+Heat+ are:
 
 \begin{lstlisting}{}
-class Heat : public ProblemInstatScal
+class Heat : public ProblemInstat
 {
 public:
-  Heat(ProblemScal &heatSpace)
-    : ProblemInstatScal("heat", heatSpace)
+  Heat(ProblemStat &heatSpace)
+    : ProblemInstat("heat", heatSpace)
   {
     theta = -1.0;
-    GET_PARAMETER(0, name + "->theta", "%f", &theta);
+    Parameters::get(name + "->theta", theta);
     TEST_EXIT(theta >= 0)("theta not set!\n");
     theta1 = theta - 1.0;
   }
@@ -131,7 +130,7 @@ public:
 
 The argument \verb+heatSpace+ is a pointer to the stationary problem
 which is solved each timestep. It is directly handed to the base class
-constructor of \verb+ProblemInstatScal+. In the body of the
+constructor of \verb+ProblemInstat+. In the body of the
 constructor, $\theta$ is read from the parameter file and stored in a
 member variable. The member variable \verb+theta1+ stores the value of
 $\theta - 1$. A pointer to this value is used later as factor in the
@@ -148,7 +147,7 @@ The next lines show the implementation of the time interface.
 
   void closeTimestep(AdaptInfo *adaptInfo) 
   {
-    ProblemInstatScal::closeTimestep(adaptInfo);
+    ProblemInstat::closeTimestep(adaptInfo);
     WAIT;
   }
 \end{lstlisting}
@@ -163,7 +162,7 @@ boundary function $g$ at $t^{new}$, which is the current time.
 
 The method \verb+closeTimestep+ is called at the end of each timestep
 by the adaptation loop. In the default implementation of
-\verb+ProblemInstatScal::closeTimestep+, the solution is written to
+\verb+ProblemInstat::closeTimestep+, the solution is written to
 output files, if specified in the parameter file. Note that the base
 class implementation of a method must be explicitly called, if the
 method is overwritten in a sub class. The macro \verb+WAIT+ waits
@@ -197,7 +196,7 @@ the initial problem by implementing this interface.
 
 Here, only the solve and the estimate step are overloaded. For the
 other steps, there are empty default implementations in
-\verb+ProblemInstatScal+. Since the mesh is not adapted in the initial
+\verb+ProblemInstat+. Since the mesh is not adapted in the initial
 problem, the initial adaptation loop will stop after one iteration. In
 the solve step, the exact solution is interpolated on the macro mesh
 and stored in the solution vector of the stationary problem. In the
@@ -242,7 +241,7 @@ int main(int argc, char** argv)
   Parameters::readArgv(argc, argv);
 
   // ===== create and init stationary problem =====
-  ProblemScal heatSpace("heat->space");
+  ProblemStat heatSpace("heat->space");
   heatSpace.initialize(INIT_ALL);
 
   // ===== create instationary problem =====
@@ -253,7 +252,7 @@ int main(int argc, char** argv)
 So far, the stationary space problem \verb+heatSpace+ and the
 instationary problem \verb+heat+ were created and
 initialized. \verb+heatSpace+ is an instance of
-\verb+ProblemScal+. \verb+heat+ is an instance of the class
+\verb+ProblemStat+. \verb+heat+ is an instance of the class
 \verb+Heat+ we defined above. \verb+heatSpace+ is given to \verb+heat+
 as its stationary problem.
 
@@ -277,9 +276,9 @@ and of the instationary adaptation loop:
 
 The object \verb+heatSpace+ is handed as
 \verb+ProblemIterationInterface+ (implemented by class\\
-\verb+ProblemScal+) to the adaptation loop. \verb+heat+ is interpreted
+\verb+ProblemStat+) to the adaptation loop. \verb+heat+ is interpreted
 as \verb+ProblemTimeInterface+ (implemented by class
-\verb+ProblemInstatScal+).
+\verb+ProblemInstat+).
 
 The function $g$ is declared in the following way:
 \begin{lstlisting}{}
@@ -309,11 +308,11 @@ Now, we define the operators:
   // create laplace
   Operator A(heatSpace.getFeSpace());
   A.addSecondOrderTerm(new Laplace_SOT);
-  A.setUhOld(heat.getOldSolution());
+  A.setUhOld(heat.getOldSolution(0));
   if (*(heat.getThetaPtr()) != 0.0)
-    heatSpace.addMatrixOperator(A, heat.getThetaPtr(), &one);
+    heatSpace.addMatrixOperator(A, 0, 0, heat.getThetaPtr(), &one);
   if (*(heat.getTheta1Ptr()) != 0.0)
-    heatSpace.addVectorOperator(A, heat.getTheta1Ptr(), &zero);
+    heatSpace.addVectorOperator(A, 0, heat.getTheta1Ptr(), &zero);
 \end{lstlisting}
 
 Operator \verb+A+ represents $-\Delta u$. It is used as matrix
@@ -331,9 +330,9 @@ operator by \verb+setUhOld+.
   // create zero order operator
   Operator C(heatSpace.getFeSpace());
   C.addZeroOrderTerm(new Simple_ZOT);
-  C.setUhOld(heat.getOldSolution());
-  heatSpace.addMatrixOperator(C, heat.getInvTau(), heat.getInvTau());
-  heatSpace.addVectorOperator(C, heat.getInvTau(), heat.getInvTau());
+  C.setUhOld(heat.getOldSolution(0));
+  heatSpace.addMatrixOperator(C, 0, 0, heat.getInvTau(), heat.getInvTau());
+  heatSpace.addVectorOperator(C, 0, heat.getInvTau(), heat.getInvTau());
 \end{lstlisting}
 
 The \verb+Simple_ZOT+ of operator \verb+C+ represents the zero order
@@ -347,7 +346,7 @@ the function \verb+getInvTau()+, which is a member of the class
   // create RHS operator
   Operator F(heatSpace.getFeSpace());
   F.addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct));
-  heatSpace.addVectorOperator(F);
+  heatSpace.addVectorOperator(F, 0);
 \end{lstlisting}
 
 \verb+CoordsAtQP_ZOT+ is a zero order term that evaluates a given
@@ -365,7 +364,7 @@ boundary function. To the last, the adaption loop is started:
   G *boundaryFct = new G;
   boundaryFct->setTimePtr(heat.getTime());
   heat.setExactSolution(boundaryFct);
-  heatSpace.addDirichletBC(1, boundaryFct);
+  heatSpace.addDirichletBC(1, 0, 0, boundaryFct);
 
   // ===== start adaption loop =====
   int errorCode = adaptInstat.adapt();
@@ -397,8 +396,8 @@ heat->adapt->end time:               1.0
 Now, tolerances for the space and the time error are determined:
 
 \begin{lstlisting}{}
-heat->adapt->tolerance:              0.05
-heat->adapt->time tolerance:         0.05
+heat->adapt[0]->tolerance:              0.05
+heat->adapt[0]->time tolerance:         0.05
 \end{lstlisting}
 
 \begin{lstlisting}{}
@@ -413,9 +412,9 @@ The following lines determine, whether coarsening is allowed in
 regions with sufficient small errors, and how many refinements or
 coarsenings are performed for marked elements.
 \begin{lstlisting}{}
-heat->adapt->coarsen allowed:        1
-heat->adapt->refine bisections:      2
-heat->adapt->coarsen bisections:     2
+heat->adapt[0]->coarsen allowed:        1
+heat->adapt[0]->refine bisections:      2
+heat->adapt[0]->coarsen bisections:     2
 \end{lstlisting}
 
 Now, the output behavior is determined:
diff --git a/doc/tutorial/installation.tex b/doc/tutorial/installation.tex
index 5b27e05f350184764a282967dfe36bd27664d71a..61a1f94ab0f8b4476bb60734c9f3208b9efe6acf 100644
--- a/doc/tutorial/installation.tex
+++ b/doc/tutorial/installation.tex
@@ -4,56 +4,12 @@
 \section{Installation of the AMDiS library}
 \label{s:installation}
 
-To install the AMDiS library, the following steps must be performed: 
-\begin{enumerate}
-
-\item Create the AMDiS source directory. This can be done in different ways. Here, we show three possibilities:
-  \begin{itemize}
-  \item Unpacking the archive file {\tt AMDiS.tgz}:
-    \begin{itemize}
-    \item {\tt > tar xvfz AMDiS.tgz}
-    \end{itemize}
-  \item Checking out an SVN project:
-    \begin{itemize}
-    \item {\tt > svn checkout --username developername \\}
-      {\tt https://gforge.zih.tu-dresden.de/svn/amdis/trunk amdis}
-    \end{itemize}
-  \end{itemize}
-
-\item Change into the AMDiS directory:
-\begin{verbatim}
-> cd amdis/AMDiS
-\end{verbatim}
-
-\item Create the makefiles for your system using the \verb+configure+ script:
-\begin{verbatim}
-> ./configure <CONFIGURE-OPTIONS>
-\end{verbatim}
-  The \verb+configure+ script creates the needed makefiles. It can be
-  called with the following options:
-\begin{description}
-\item \verb+--prefix=<AMDIS-DIR>+: Installation path of the AMDiS
-  library. The library file will be stored in
-  \verb+<AMDIS-DIR>/lib+. With \verb+--prefix=`pwd`+ AMDiS will be
-  installed in the current working directory, which mostly is a good
-  choice.
-\item \verb+--enable-debug+: If this option is used, AMDiS will be
-  compiled with debug support. By default, AMDiS is compiled in an
-  optimized mode without debug support.
-\item \verb+--enable-umfpack+: If this option is used, AMDiS is
-  compiled with support for the UMFPACK solver library.
-\end{description}
-
-\item Make the library:
-\begin{verbatim}
-> make install
-\end{verbatim}
-\end{enumerate}
-If you have added a new source file or you want to change something on the automake-system, you have to rerun the following command:
-\begin{verbatim}
-./autogen <CONFIGURE-OPTION>
-\end{verbatim}
-This will generate files for the configure-make-system.
+There are two different ways to install AMDiS, depending if you have a
+binary packages downloaded from www.amdis-fem.org, or you have a
+source package. In the first case, you just need to run the package
+manager of your Linux distribution and install this package. To
+install AMDiS from the sources, you need CMake with version at least
+2.8. Follow the hints in the file AMDiS/Howto\_cmake.html.
 
 \section{Compilation of an example application}
 \label{s:compilation}
@@ -63,21 +19,17 @@ For the compilation of the examples, described in this section, the following st
 \item Change into the demo directory:\\
   \verb+> cd demo+
 
-\item It may be required to change the \verb+Makefile+:
-  \begin{itemize}
-  \item Set the AMDiS path and paths of other needed libraries.
-  \item Set user flags.
-  \end{itemize}  
-  The makefile is described in Section \ref{s:application makefile} in detail.
+\item call cmake:\\
+\verb+> ccmake .+
 
+\item if AMDiS is not found automatically, set \verb+AMDIS_DIR+ to
+  the directory of \verb+AMDiS/share/amdis+
 \item Make the application example:\\
   \verb+> make <PROG-NAME>+\\
 
   \verb+<PROG-NAME>+ is the name of the application example.
-
 \end{enumerate}
 To run the example, call:\\
   \verb+> ./<PROG-NAME> <PARAMETER-FILE>+
-
 \input{makefile.tex}
 
diff --git a/doc/tutorial/neumann.tex b/doc/tutorial/neumann.tex
index 2b20165ebed863bf4f5c76d119566ee6fb39778a..c888ba13791bb1d792c1fab1665db14b03dff1c8 100644
--- a/doc/tutorial/neumann.tex
+++ b/doc/tutorial/neumann.tex
@@ -1,7 +1,13 @@
 \section{Neumann boundary conditions}
 \label{s:neumann}
-In this example, we solve the problem defined in Section \ref{s:ellipt}. But now, we set the domain $\Omega$ to $[-0.5;0.5]^2$, so the source $f$ is located in the middle of $\Omega$. Furthermore, we use Neumann boundary conditions on the left and on the right side of $\Omega$. We set $A\nabla u \cdot \nu = 1$ at the Neumann boundary.
-So, the derivative in direction of the surface normal is set to $1$ at these points. The rest of the boundary keeps unchanged (Dirichlet boundary, set to the true solution).
+In this example, we solve the problem defined in Section
+\ref{s:ellipt}. But now, we set the domain $\Omega$ to $[-0.5;0.5]^2$,
+so the source $f$ is located in the middle of $\Omega$. Furthermore,
+we use Neumann boundary conditions on the left and on the right side
+of $\Omega$. We set $A\nabla u \cdot \nu = 1$ at the Neumann boundary.
+So, the derivative in direction of the surface normal is set to $1$ at
+these points. The rest of the boundary keeps unchanged (Dirichlet
+boundary, set to the true solution).
 
 \begin{table}
 \center
@@ -41,16 +47,23 @@ In the main program we add the boundary conditions to our problem \verb+neumann+
 int main(int argc, char* argv[])
 {
   ...
-  neumann.addNeumannBC(1, new N);
-  neumann.addDirichletBC(2, new G);
+  neumann.addNeumannBC(1, 0, 0, new N);
+  neumann.addDirichletBC(2, 0, 0, new G);
   ...
 }
 \end{lstlisting}
 
-Since the Dirichlet condition has a higher ID, it has a higher priority against the Neumann boundary condition. This is important, where different conditions meet each other in some points. In this example, these are the corner points of $\Omega$. If Dirichlet boundary conditions are used together with boundary conditions of other type, the Dirichlet conditions should always have the higher priority. 
+Since the Dirichlet condition has a higher ID, it has a higher
+priority against the Neumann boundary condition. This is important,
+where different conditions meet each other in some points. In this
+example, these are the corner points of $\Omega$. If Dirichlet
+boundary conditions are used together with boundary conditions of
+other type, the Dirichlet conditions should always have the higher
+priority.
 
 \subsection{Parameter file}
-In the parameter file, we use the file \verb+./macro/neumann.macro.2d+ as macro mesh file, described in the next section.
+In the parameter file, we use the file \verb+./macro/neumann.macro.2d+
+as macro mesh file, described in the next section.
 
 \subsection{Macro file}
 The file \verb+neumann.macro.2d+ is listed below:
@@ -82,10 +95,17 @@ element boundaries:
 0 0 1
 \end{lstlisting}
 
-In contrast to the standard file \verb+macro.stand.2d+, here the vertex coordinates are shifted to describe the domain $[-0.5;0.5]^2$. Furthermore, the boundary block changed. The elements 0 and 2 have the Dirichlet boundary with ID 2 at edge 2. Elements 1 and 3 have the Neumann boundary condition with ID 1 applied to their local edge 2.
+In contrast to the standard file \verb+macro.stand.2d+, here the
+vertex coordinates are shifted to describe the domain
+$[-0.5;0.5]^2$. Furthermore, the boundary block changed. The elements
+0 and 2 have the Dirichlet boundary with ID 2 at edge 2. Elements 1
+and 3 have the Neumann boundary condition with ID 1 applied to their
+local edge 2.
 
 \subsection{Output}
-In Figure \ref{f:neumann solution}, the solution is shown. At the Neumann boundaries, one can see the positive slope. At Dirichlet boundaries, the solution is set to $g(x)$. 
+In Figure \ref{f:neumann solution}, the solution is shown. At the
+Neumann boundaries, one can see the positive slope. At Dirichlet
+boundaries, the solution is set to $g(x)$.
 \begin{figure}
 \center
 \includegraphics[width=0.5\textwidth]{fig/neumann_solution.jpg}
diff --git a/doc/tutorial/parallelheat.tex b/doc/tutorial/parallelheat.tex
deleted file mode 100644
index dd03fd21c2ea3220bbf4e107e4c361436dce7d13..0000000000000000000000000000000000000000
--- a/doc/tutorial/parallelheat.tex
+++ /dev/null
@@ -1,241 +0,0 @@
-\section{Parallelization}
-
-\begin{figure}
-\center
-\begin{tabular}{cccc}
-\includegraphics[width=0.21\textwidth]{fig/macro.pdf} &
-\includegraphics[width=0.21\textwidth]{fig/local1.pdf} &
-\includegraphics[width=0.21\textwidth]{fig/global.pdf} &
-\includegraphics[width=0.21\textwidth]{fig/local2.pdf} \\
-(a)&(b)&(c)&(d)\\
-\end{tabular}
-\caption{(a) A triangular macro mesh, (b) domain decomposition after six global refinements and overlap for partition 0, (c) composite mesh after adaptation loop, (d) local mesh of process 0 after adaptation loop}
-\label{f:full domain covering meshes2}
-\end{figure}
-
-Before we start with tha application example, we give a short overview over the parallelization approach of full domain covering meshes used in AMDiS. The approach can be summerized by the following steps:
-
-\begin{itemize}
-\item Create a partitioning level by some adaptive or global refinements of the macro mesh.
-\item Create a partitioning with approximately the same number of elements in each partition.
-\item Every process computes on the whole domain, but adaptive mesh refinements are only allowed within the local partition including some overlap.
-\item After each adaptive iteration, a new partitioning can be computed if the parallel load balance becomes to bad. Partitionings are always computed at the same mesh level. The partitioning elements are weighted by the number of corresponding leaf elements (elements with no children).
-\item At the end of each timestep or at the end of parallel computations, a global solution is constructed by a partititon of unity method (weighted sum over all process solutions). 
-\end{itemize}
-Figure \ref{f:full domain covering meshes2} illustrates the concept of full domain covering meshes.
-
-The so called three level approach allows to decouple the following levels:
-\begin{enumerate}
-\item Partitioning level: Basis for the partitioning. The partitioning level is built by the current leaf elements of the mesh, when the parallelization is initialized.
-\item Global coarse grid level: Number of global refinements starting from the partitioning level for the whole domain. The global coarse grid level builds the coarsest level for all future computations. 
-\item Local coarse grid level: Number  of global refinements starting from the partitioning level within the local domain including neighbor elements. The local coarse grid level reduces the size of partition overlap.
-\end{enumerate}
-In Figure \ref{f:threelevels2}, an example for the three level approach is given. The dashed line shows the overlap of size 1 for the local partition computed at partitining level.
-
-The parallelization in AMDiS uses the {\it Message Passing Interface} MPI, see
-\begin{lstlisting}{}
-http://www-unix.mcs.anl.gov/mpi/ .
-\end{lstlisting}
-
-
-For the partitioning the parallel graph partitioning library ParMETIS is used, see
-\begin{lstlisting}{}
-http://glaros.dtc.umn.edu/gkhome/metis/parmetis/overview .
-\end{lstlisting}
-
-To use AMDiS in parallel, it must be configured like
-\begin{lstlisting}{}
-> ./configure --prefix=<AMDIS-DIR> 
-              --with-mpi=<MPI-DIR>
-              --with-parmetis=<PARMETIS-DIR> .
-\end{lstlisting}
-
-The parallel application (in this example \verb+parallelheat+) must be compiled with the MPI C++ compiler \verb+mpiCC+ and linked against the ParMETIS library. Then the application can be called with \verb+mpirun+:
-
-\begin{lstlisting}{}
-> mpirun -np <num-procs> ./parallelheat init/parallelheat.dat.2d
-\end{lstlisting}
-
-Now we start with the example. We want to solve a time dependent problem as described in \ref{s:heat}. As right hand side function $f$ we choose the {\it moving source} 
-\begin{equation}
-f(x,t)=\mbox{sin}(\pi t) e^{-10(x-\vec{t})^2}
-\end{equation}
-with $0 \leq t \leq 1$ and $\vec{t}$ a vector with all components set to $t$. The maximum of this function moves from $(0,0)$ to $(1,1)$ within the specified time intervall. 
-We use 4 parallel processes.
-
-\begin{figure}
-\center
-\includegraphics[width=0.95\textwidth]{fig/threelevels.pdf}
-\caption{Example for the three levels of partitioning. The partitioning mesh is globally refined twice to get the global coarse grid level. Then two further global refinement steps applied on $\Omega_i$ and its direct neighbor elements (in each step) result in the local coarse grid level.}
-\label{f:threelevels2}
-\end{figure}
-
-\begin{table}
-\center
-\begin{tabular}{|cc|c|c|c|}
-\hline 
-&& {\bf 1d} & {\bf 2d} & {\bf 3d} \\
-\hline 
-{\bf source code} & {\tt src/} & \multicolumn{3}{|c|}{\tt parallelheat.cc} \\
-\hline 
-{\bf parameter file} & {\tt init/} & \tt - & \tt parallelheat.dat.2d & \tt - \\
-\hline 
-{\bf macro file} & {\tt macro/} & \tt - & \tt macro.stand.2d & \tt - \\
-\hline 
-{\bf output files} & {\tt output/} & \multicolumn{3}{|c|}{\tt parallelheat\_proc<n>\_<t>.mesh/.dat} \\
-\hline 
-\end{tabular}
-\caption{Files of the {\tt parallelheat} example. In the output file names, {\tt <n>} is replaced by the process number and {\tt <t>} is replaced by the time.}
-\label{t:parallelheat}
-\end{table}
-
-
-\subsection{Source code}
-
-In this section the interesting aspects of file \verb+parallelheat.cc+ are described. First, the moving functions $f$ and $g$ are defined.
-
-\begin{lstlisting}{}
-class G : public AbstractFunction<double, WorldVector<double> >,
-	  public TimedObject
-{
-public:
-  /// Implementation of AbstractFunction::operator().
-  double operator()(const WorldVector<double>& argX) const
-  {
-    WorldVector<double> x = argX;
-    int dim = x.getSize();
-    for (int i = 0; i < dim; i++)
-      x[i] -= *timePtr;
-
-    return sin(M_PI * (*timePtr)) * exp(-10.0 * (x * x));
-  }
-};
-
-class F : public AbstractFunction<double, WorldVector<double> >,
-	  public TimedObject
-{
-public:
-  F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}
-
-  /// Implementation of AbstractFunction::operator().
-  double operator()(const WorldVector<double>& argX) const 
-  {
-    WorldVector<double> x = argX;
-    int dim = x.getSize();
-    for (int i = 0; i < dim; i++)
-      x[i] -= *timePtr;
-
-    double r2 = (x * x);
-    double ux = sin(M_PI * (*timePtr)) * exp(-10.0 * r2);
-    double ut = M_PI * cos(M_PI*(*timePtr)) * exp(-10.0 * r2);
-    return ut - (400.0 * r2 - 20.0 * dim) * ux;
-  }
-};
-\end{lstlisting}
-
-The main program starts with \verb+MPI::Init+ to initialize MPI, and ends with \verb+MPI::Finalize+ to finalize MPI.
-
-\begin{lstlisting}{}
-int main(int argc, char** argv)
-{
-  MPI::Init(argc, argv);
-
-  ...
-
-  std::vector<DOFVector<double>*> vectors;
-  ParallelProblemScal parallelheat("heat->parallel", heatSpace, heat, 
-                                   vectors);
-
-  AdaptInstationary *adaptInstat =
-    new AdaptInstationary("heat->adapt",
-			  &parallelheat,
-			  adaptInfo,
-			  &parallelheat,
-			  adaptInfoInitial);
-
-  ... 
-
-  parallelheat.initParallelization(adaptInfo);
-  adaptInstat->adapt();
-  parallelheat.exitParallelization(adaptInfo);
-
-  MPI::Finalize();
-}
-\end{lstlisting}
-
-The parallel problem \verb+parallelheat+ has to know the sequential instationary problem \verb+heat+ and the space problem \verb+heatSpace+. The vector \verb+vectors+ stores pointers to DOF vectors which are used by the operator terms. The values of these vectors have to be exchanged during repartitioning. The solution of the last time step is considered automatically. Hence, \verb+vectors+ here is empty.
-
-The adaptation loop \verb+adaptInstat+ uses \verb+parallelheat+ as iteration interface and as time interface. Before the adaptation loop starts, the parallelization is initialized by \verb+initParallelization+. After the adaptation loop has finished, the parallelization is finalized by \verb+exitParallelization+.
-
-\subsection{Parameter file}
-
-We use the parameter file described in Section \ref{s:heat parameter} as basis for the file \verb+parallelheat.dat.2d+ and describe the relevant changes for the parallel case.
-
-\begin{lstlisting}{}
-heatMesh->global refinements:         6
-\end{lstlisting}
-
-The macro mesh is globally refined 6 times to create the partitioning mesh which is used as basis for the domain decomposition. 
-
-In the following, one can see the parameters for the parallel problem \verb+heat->parallel+.
-
-\begin{lstlisting}{}
-heat->parallel->upper part threshold: 1.5
-heat->parallel->lower part threshold: 0.66
-\end{lstlisting}
-
-These two parameters determine the upper and lower partitioning thresholds $rt_{high}$ and $rt_{low}$. If at least one process exceeds the number of $rt_{high} \cdot sum_{av}$ elements or falls below $rt_{low} \cdot sum_{av}$ elements, a repartitioning must be done ($sum_{av}$ is the average number of elements per partition). 
-
-%If adaptive thresholds are used, these values are also the limits for the threshold adaptation. The start value for the upper threshold is also the minimum for the upper threshold, the start value for the lower threshold is also the maximum for the lower threshold.
-
-%\begin{lstlisting}{}
-%heat->parallel->adaptive thresholds:  1
-%heat->parallel->threshold inc factor: 2.0
-%heat->parallel->threshold dec factor: 0.5
-%heat->parallel->repart time factor:   10.0
-%\end{lstlisting}
-
-%If \verb+adaptive thresholds+ is set to $1$, the \verb+threshold inc factor+ and \verb+threshold dec factor+ determine the incremenent and decrement factors for the thresholds. The value of \verb+repart time factor+ is the desired relation between repartitioning time and computation time.
-
-\begin{lstlisting}{}
-heat->parallel->global coarse grid level: 0
-heat->parallel->local coarse grid level:  4
-\end{lstlisting}
-
-Here, the parameters used for the three level approach are given. The partitioning level has already been determined by \verb+heatMesh->global refinements+. The value of \verb+global coarse grid level+ is set to $0$, and \verb+local coarse grid level+ is set to $4$. This means that 4 global refinements are done within the local partititon at each process (incl. neighoring elements), to reduce the overlap size. No further refinements are done outside of the local partitions. 
-
-\begin{lstlisting}{}
-heat->adapt->timestep:               0.1
-heat->adapt->min timestep:           0.1
-heat->adapt->max timestep:           0.1
-heat->adapt->start time:             0.0
-heat->adapt->end time:               1.0
-\end{lstlisting}
-
-We use a fixed timestep of $0.1$. The start time is set to $0.0$, the end time is set to $1.0$. Thus, we have 10 timesteps (without the initial solution at $0.0$). 
-
-\begin{lstlisting}{}
-heat->space->output->filename:       output/heat
-\end{lstlisting}
-
-For the output the prefix \verb+output/heat+ is used. Each process adds its process ID. Then the time is added and finally the file extensions. The result of the initial problem of process 1 will be written to \verb+output/heat_proc0_00.000.mesh+ and \verb+output/heat_proc0_00.000.dat+.
-
-\subsection{Macro file}
-We again use the macro file \verb+macro/macro.stand.2d+, which was described in Section \ref{s:ellipt macro}.
-
-\subsection{Output}
-For each timestep and for each process, one mesh file and one value file is written (total file number: $11 \cdot 4 \cdot 2 = 88$). In Figure \ref{f:parallel heat}, the local process solutions for $t=0.1$, $t=0.5$ and $t=0.9$ are visualized after the partititon of unity was applied. One can see that the partitioning considers the moving source.
-
-\begin{figure}
-\center
-\begin{tabular}{ccc}
-\includegraphics[width=0.3\textwidth]{fig/parallelheat_01.jpg}&
-\includegraphics[width=0.3\textwidth]{fig/parallelheat_05.jpg}&
-\includegraphics[width=0.3\textwidth]{fig/parallelheat_09.jpg}\\
-(a)&(b)&(c)
-\end{tabular}
-\caption{Local process solutions for $t=0.1$ (a), $t=0.5$ (b), and $t=0.9$ (c).}
-\label{f:parallel heat}
-\end{figure}
-
-
diff --git a/doc/tutorial/parameters.tex b/doc/tutorial/parameters.tex
index 15065709b03da7fce2ed8888173facdee161f68f..639b6f22329e69fb24840900583fc5b1f57bcf27 100644
--- a/doc/tutorial/parameters.tex
+++ b/doc/tutorial/parameters.tex
@@ -1,47 +1,41 @@
 \section{Parameter reference}
 \label{s:parameter reference}
 
-This section gives an overview over all available AMDiS parameters. The parameters are defined by the user in an init-file which is usually handed to the AMDiS application as command line parameter. 
-A parameter in the init-file is defined by
+This section gives an overview over all available AMDiS
+parameters. The parameters are defined by the user in an init-file
+which is usually handed to the AMDiS application as command line
+parameter.  A parameter in the init-file is defined by
 \begin{lstlisting}{}
 <parameter-name>: <parameter-value> % comment
 \end{lstlisting}
 
 \noindent The order of parameters in the init file is arbitrary. 
 
-Within the AMDiS library parameters are read by the macro
-
+Within the AMDiS library parameters are read by the function
 \begin{lstlisting}{}
-GET_PARAMETER(<info-level>, "<parameter-name>", "<format-string>", <var-ptr>)+
+Parameters::get("<parameter-name>", <var>)+
 \end{lstlisting}
 
 \noindent The parameters are
 \begin{itemize}
-\item \verb+<info-level>+: Level of information printed during parameter initialization ($0-4$).
-\item \verb+"<parameter-name>"+: Name of the parameter in the input file (without the \verb+:+).
-\item \verb+"<format-string>"+: Format string. Determines how the \verb+<parameter-value>+ string in the input file is interpreted. It can be
-  \begin{itemize}
-  \item \verb+"%c"+: Character
-  \item \verb+"%d"+: Signed decimal integer
-  \item \verb+"%e"+: Scientific notation (mantise/exponent) using \verb+e+ character
-  \item \verb+"%E"+: Scientific notation (mantise/exponent) using \verb+E+ character
-  \item \verb+"%f"+: Decimal floating point
-  \item \verb+"%u"+: Unsigned decimal integer
-  \end{itemize}
-
-\item \verb+<var-ptr>+: Address of a variable the value should be written to. If the corresponding \verb+"<parameter-name>"+ wasn't found in the init-file, the old value of the variable isn't changed and it keeps it default value, which was set before.
+\item \verb+"<parameter-name>"+: Name of the parameter in the input
+  file (without the \verb+:+).
+\item \verb+<var>+: Variable the value should be written to. If the
+  corresponding \verb+"<parameter-name>"+ wasn't found in the
+  init-file, the old value of the variable isn't changed and it keeps
+  it default value, which was set before.
 \end{itemize}
 
-To read a STL string from the init file another version of \verb+GET_PARAMETER+ can be used where the format string is left out:
-
-\begin{lstlisting}{}
-GET_PARAMETER(<info-level>, "<parameter-name>", &myString);
-\end{lstlisting}
-
-Names of parameters that belong to a specific object usually start with a prefix that uniquely identifies the object followed by \verb+"->"+. The prefix should be handed to the objects constructor, because there mostly also the parameters are initialized. Assume that we have a problem with name \verb+"ellipt"+ and we want to set its solver tolerance to $10^{-4}$. The solver is initialized with the prefix \verb+"ellipt->solver"+. The statement 
+Names of parameters that belong to a specific object usually start
+with a prefix that uniquely identifies the object followed by
+\verb+"->"+. The prefix should be handed to the objects constructor,
+because there mostly also the parameters are initialized. Assume that
+we have a problem with name \verb+"ellipt"+ and we want to set its
+solver tolerance to $10^{-4}$. The solver is initialized with the
+prefix \verb+"ellipt->solver"+. The statement
 
 \begin{lstlisting}{}
-GET_PARAMETER(0, prefix + "->tolerance", "%f", &tolerance);
+Parameters::get(prefix + "->tolerance", tolerance);
 \end{lstlisting}
 
 \noindent within the solver constructor reads the init-file line
@@ -50,15 +44,22 @@ GET_PARAMETER(0, prefix + "->tolerance", "%f", &tolerance);
 ellipt->solver->tolerance: 1e-4
 \end{lstlisting}
 
-\noindent and writes the value $10^{-4}$ to the floating point variable \verb+tolerance+.
+\noindent and writes the value $10^{-4}$ to the floating point
+variable \verb+tolerance+.
 
-In the following sections all AMDiS parameters are listed categorized by the classes they belong to. The user can also create its own parameters in the way described above.
+In the following sections all AMDiS parameters are listed categorized
+by the classes they belong to. The user can also create its own
+parameters in the way described above.
 
 \subsection{Adapt classes}
 
 \subsection{AdaptInfo}
 
-The prefix is defined in the main program (constructor of {\it AdaptInfo}). Usually it is set to \verb#<problem-name> + "->adapt->"#, where \verb#<problem-name># is the name of the problem the \verb+AdaptInfo+ belongs to (the + sign here means string concatenation).
+The prefix is defined in the main program (constructor of {\it
+  AdaptInfo}). Usually it is set to
+\verb#<problem-name> + "->adapt->"#, where \verb#<problem-name># is
+the name of the problem the \verb+AdaptInfo+ belongs to (the + sign
+here means string concatenation).
 
 \begin{itemize}
 \item \verb+"start time"+ (type: \verb+double+, default: $0.0$): Start time of the simulation.
@@ -87,7 +88,10 @@ An {\it AdaptInfo} instance additionally has parameters that are connected to th
 
 \subsection{AdaptStationary}
 
-The prefix is defined in the main program (constructor of {\it AdaptStationary}). Usually it is set to \verb#<problem-name> + "->adapt->"#, where \verb#<problem-name># is the name of the problem the {\it AdaptStationary} instance belongs to.
+The prefix is defined in the main program (constructor of {\it
+  AdaptStationary}). Usually it is set to
+\verb#<problem-name> + "->adapt->"#, where \verb#<problem-name># is
+the name of the problem the {\it AdaptStationary} instance belongs to.
 
 \begin{itemize}
 \item \verb+""+ (type: \verb++, default: $1$):
@@ -100,13 +104,6 @@ The prefix is defined in the main program (constructor of {\it AdaptStationary})
 \item \verb+""+ (type: \verb++, default: $1$):
 \end{itemize}
 
-%AdaptInstationary.cc:170:    GET_PARAMETER(0, aName + "->strategy", "%d", &strategy);
-%AdaptInstationary.cc:171:    GET_PARAMETER(0, aName + "->time delta 1", "%f", &time_delta_1);
-%AdaptInstationary.cc:172:    GET_PARAMETER(0, aName + "->time delta 2", "%f", &time_delta_2);
-%AdaptInstationary.cc:173:    GET_PARAMETER(0, aName + "->info", "%d", &info_);
-%AdaptInstationary.cc:174:    GET_PARAMETER(0, aName + "->break when stable", "%d", &breakWhenStable);
-%AdaptStationary.cc:49:  GET_PARAMETER(0, name_ + "->info", "%d", &info_);
-
 \subsection{AdaptInstationary}
 
 \begin{itemize}
diff --git a/doc/tutorial/periodic.tex b/doc/tutorial/periodic.tex
index edd3b5a77b65e7c9a16602b402ba0962da21b0c6..727c108504c27c3ce7087780bf89f76198916713 100644
--- a/doc/tutorial/periodic.tex
+++ b/doc/tutorial/periodic.tex
@@ -8,23 +8,64 @@
 \label{f:periodic boundaries}
 \end{figure} 
 
-Periodic boundary conditions allow to simulate an effectively infinite tiled domain, where the finite domain $\Omega$ is interpreted as one tile of the infinte problem domain. The solution outside of $\Omega$ can be constructed by periodically continue the solution within $\Omega$. In Figure \ref{f:periodic boundaries} two examples for periodic boundary conditions on a two dimensional domain are illustrated. On the left hand side example, the upper and the lower part of the boundary as well as the left and the right part of the boundary are assigned to each other as periodic boundary. This results in a solution, which tiles the infinte plane. On the right hand side example, only the left and the right part of the domain boundary are assigned to each other, which results in a infinte band.   
+Periodic boundary conditions allow to simulate an effectively infinite
+tiled domain, where the finite domain $\Omega$ is interpreted as one
+tile of the infinte problem domain. The solution outside of $\Omega$
+can be constructed by periodically continue the solution within
+$\Omega$. In Figure \ref{f:periodic boundaries} two examples for
+periodic boundary conditions on a two dimensional domain are
+illustrated. On the left hand side example, the upper and the lower
+part of the boundary as well as the left and the right part of the
+boundary are assigned to each other as periodic boundary. This results
+in a solution, which tiles the infinte plane. On the right hand side
+example, only the left and the right part of the domain boundary are
+assigned to each other, which results in a infinte band.
 
 \begin{figure}
 \center
 \includegraphics[width=0.5\textwidth]{fig/periodic_identification.pdf}
-\caption{A one dimensional mesh with vertex coordinates stored at the elements and a corresponding periodic mesh with changed mesh topology. Note that the geometric data are not changed. The coordinates of the first vertex depend on the element it belongs to.}
+\caption{A one dimensional mesh with vertex coordinates stored at the
+  elements and a corresponding periodic mesh with changed mesh
+  topology. Note that the geometric data are not changed. The
+  coordinates of the first vertex depend on the element it belongs
+  to.}
 \label{f:periodic identification}
 \end{figure} 
-In AMDiS, there are two ways to implement periodic boundary conditions:
+In AMDiS, there are two ways to implement periodic boundary
+conditions:
 \begin{enumerate}
-\item Changing the mesh topology ({\it mode 0}): Before the computation is started, the topology of the macro mesh is changed. Two vertices that are assigned to each other by a periodic boundary condition, are replaced by one single vertex, which is now treated as an inner vertex of the mesh (if it is not part of any other boundary). Since geometric data like coordinates are stored at elements and not at vertices, this modification does not change the geometry of the problem. Topological information, like element neighborhood, does change. The method is illustrated in Figure \ref{f:periodic identification}.
-\item Modify the linear system of equations in each iteration ({\it mode 1}): Sometimes it is necessary to store geometric information at vertices. E.g., if moving meshes are implemented with parametric elements, a DOF vector may store the coordinates. In this case, the mesh topology keeps unchanged, and the periodic boundary conditions are applied, like any other boundary condition, after the assemblage of the linear system of equations. In the application source code a boundary condition object has to be created, and in the macro file the periodic boundary must be specified.
+\item Changing the mesh topology ({\it mode 0}): Before the
+  computation is started, the topology of the macro mesh is
+  changed. Two vertices that are assigned to each other by a periodic
+  boundary condition, are replaced by one single vertex, which is now
+  treated as an inner vertex of the mesh (if it is not part of any
+  other boundary). Since geometric data like coordinates are stored at
+  elements and not at vertices, this modification does not change the
+  geometry of the problem. Topological information, like element
+  neighborhood, does change. The method is illustrated in Figure
+  \ref{f:periodic identification}.
+\item Modify the linear system of equations in each iteration ({\it
+    mode 1}): Sometimes it is necessary to store geometric information
+  at vertices. E.g., if moving meshes are implemented with parametric
+  elements, a DOF vector may store the coordinates. In this case, the
+  mesh topology keeps unchanged, and the periodic boundary conditions
+  are applied, like any other boundary condition, after the assemblage
+  of the linear system of equations. In the application source code a
+  boundary condition object has to be created, and in the macro file
+  the periodic boundary must be specified.
 \end{enumerate}
 
-In this section, we show how to use both variants of periodic boundary conditions. Again, we use the problem defined in Section \ref{s:ellipt}. We choose $\Omega = [-0.2; 0.8]\times[-0.5;0.5]$ (we do not use $\Omega = [-0.5; 0.5]^2$, because in this example periodic boundary conditions would then be equal to the trivial zero flux conditions).
+In this section, we show how to use both variants of periodic boundary
+conditions. Again, we use the problem defined in Section
+\ref{s:ellipt}. We choose $\Omega = [-0.2; 0.8]\times[-0.5;0.5]$ (we
+do not use $\Omega = [-0.5; 0.5]^2$, because in this example periodic
+boundary conditions would then be equal to the trivial zero flux
+conditions).
 
-We apply a periodic boundary condition which connects the left and the right edge of $\Omega$. Since we do not know the exact solution of this periodic problem, we apply zero Dirichlet conditions at the lower and upper edge of the domain. 
+We apply a periodic boundary condition which connects the left and the
+right edge of $\Omega$. Since we do not know the exact solution of
+this periodic problem, we apply zero Dirichlet conditions at the lower
+and upper edge of the domain.
 
 \begin{table}
 \center
@@ -48,13 +89,16 @@ We apply a periodic boundary condition which connects the left and the right edg
 \end{table}
 
 \subsection{Source code}
-If we use {\it mode 0}, no modifications in the source code have to be made. For {\it mode 1}, we have to add a periodic boundary condition object to the problem.
+If we use {\it mode 0}, no modifications in the source code have to be
+made. For {\it mode 1}, we have to add a periodic boundary condition
+object to the problem.
 
 \begin{lstlisting}{}
-periodic.addPeriodicBC(-1);
+periodic.addPeriodicBC(-1, 0, 0);
 \end{lstlisting}
 
-Note that periodic boundary conditions must be described by negative numbers.
+Note that periodic boundary conditions must be described by negative
+numbers.
 
 \subsection{Parameter file}
 In the parameter file, we add an link to the {\it periodic file}.
@@ -63,7 +107,8 @@ In the parameter file, we add an link to the {\it periodic file}.
 periodicMesh->periodic file:      ./init/periodic.per.2d
 \end{lstlisting}
 
-The periodic file \verb+periodic.per.2d+ contains the needed periodic information for the mesh.
+The periodic file \verb+periodic.per.2d+ contains the needed periodic
+information for the mesh.
 
 \begin{lstlisting}{}
 associations: 2
@@ -73,10 +118,25 @@ mode  bc  el1 - local vertices <->  el2 - local vertices
  1    -1   0        1 2              3        2 1
 \end{lstlisting}
 
-First, the number of edge associations (point associations in 2d, face associations in 3d) is given. Then each association is described in one line. The first entry is the mode which should be used for this periodic association. If the mode is 1, the next entry specifies the identifier of the used boundary condition. This identifier also must be used in the source code and in the macro file. If the mode is 0, the identifier is ignored. The rest of the line describes, which (local) vertices of which elements are associated with each other. The first association in this example is interpreted as follows: The local vertices 1 and 2 of element 4 are associated with the vertices 2 and 1 of element 7. Or more precisely, vertex 1 of element 4 is associated with vertex 2 of element 7, and vertex 2 of element 4 with vertex 1 of element 7.  
+First, the number of edge associations (point associations in 2d, face
+associations in 3d) is given. Then each association is described in
+one line. The first entry is the mode which should be used for this
+periodic association. If the mode is 1, the next entry specifies the
+identifier of the used boundary condition. This identifier also must
+be used in the source code and in the macro file. If the mode is 0,
+the identifier is ignored. The rest of the line describes, which
+(local) vertices of which elements are associated with each other. The
+first association in this example is interpreted as follows: The local
+vertices 1 and 2 of element 4 are associated with the vertices 2 and 1
+of element 7. Or more precisely, vertex 1 of element 4 is associated
+with vertex 2 of element 7, and vertex 2 of element 4 with vertex 1 of
+element 7.
 
 \subsection{Macro file}
-To avoid degenerated elements, one macro element must not contain two vertices which are associated with each other. Therefore, we choose a macro mesh with a few more elements, showed in Figure \ref{f:periodic macro}.
+To avoid degenerated elements, one macro element must not contain two
+vertices which are associated with each other. Therefore, we choose a
+macro mesh with a few more elements, showed in Figure \ref{f:periodic
+  macro}.
 \begin{figure}
 \center
 \includegraphics[width=0.7\textwidth]{fig/periodic_macro.pdf}
@@ -135,19 +195,31 @@ element neighbours:
  4 -1  6
 \end{lstlisting}
 
-Compared to the macro file of Section \ref{s:ellipt macro}, the vertex coordinates are shifted by -0.2 in x-direction. 
+Compared to the macro file of Section \ref{s:ellipt macro}, the vertex
+coordinates are shifted by -0.2 in x-direction.
 
-In the boundary block $-1$ specifies the periodic boundary. If we use {\it mode 0}, this boundaries are ignored (Here, the minus sign becomes important! Only boundary conditions with negative IDs are recognized as periodic boundaries and can be ignored if they are not used).
+In the boundary block $-1$ specifies the periodic boundary. If we use
+{\it mode 0}, this boundaries are ignored (Here, the minus sign
+becomes important! Only boundary conditions with negative IDs are
+recognized as periodic boundaries and can be ignored if they are not
+used).
 
-In the neighbors block, neighborships between elements that are connected by a periodic edge (point/face) are added. Note that this must also be done for {\it mode 1} periodic boundaries.
+In the neighbors block, neighborships between elements that are
+connected by a periodic edge (point/face) are added. Note that this
+must also be done for {\it mode 1} periodic boundaries.
 
 \subsection{Output}
-In Figure \ref{f:periodic solution}, the solution of our periodic problem is shown as height field at the left hand side. At the right hand side, one can see iso lines for the values $0.1, 0.2, \dots, 1.1$.
+In Figure \ref{f:periodic solution}, the solution of our periodic
+problem is shown as height field at the left hand side. At the right
+hand side, one can see iso lines for the values $0.1, 0.2, \dots,
+1.1$.
 
 \begin{figure}
 \center
 \includegraphics[width=0.45\textwidth]{fig/periodic_solution.jpg}
 \includegraphics[width=0.45\textwidth]{fig/periodic_iso.jpg}
-\caption{Solution of the problem with periodic boundary conditions at two sides (left) and iso lines of the solution for the values $0.1, 0.2, \dots ,1.1$ (right).}
+\caption{Solution of the problem with periodic boundary conditions at
+  two sides (left) and iso lines of the solution for the values $0.1,
+  0.2, \dots ,1.1$ (right).}
 \label{f:periodic solution}
 \end{figure} 
diff --git a/doc/tutorial/projections.tex b/doc/tutorial/projections.tex
index ced1cdb2b48f3f8fa399ffcc41a5320306e9c019..72ff443e6dd8ce566ee110bee05bb2b4e9a26755 100644
--- a/doc/tutorial/projections.tex
+++ b/doc/tutorial/projections.tex
@@ -1,26 +1,46 @@
 \section{Projections}
 \label{ss:projections}
 
-In AMDiS, projections can be applied to the vertex coordinates of a mesh. There are two types of projections:
+In AMDiS, projections can be applied to the vertex coordinates of a
+mesh. There are two types of projections:
 \begin{enumerate}
-\item {\it Boundary projections}: Only vertices at the domain boundary are projected.
-\item {\it Volume projections}: All vertices of the mesh are projected.
+\item {\it Boundary projections}: Only vertices at the domain boundary
+  are projected.
+\item {\it Volume projections}: All vertices of the mesh are
+  projected.
 \end{enumerate}
-Projections are applied to all (boundary) vertices of the macro mesh and to each new (boundary) vertex, created during adaptive refinements. In Figure \ref{f:boundary projection}, this is illustrated for a two dimensional mesh which boundary vertices are successively projected to a circle. In Figure \ref{f:volume projection} the vertices of a one dimensional mesh are successively projected on the circle.
+Projections are applied to all (boundary) vertices of the macro mesh
+and to each new (boundary) vertex, created during adaptive
+refinements. In Figure \ref{f:boundary projection}, this is
+illustrated for a two dimensional mesh which boundary vertices are
+successively projected to a circle. In Figure \ref{f:volume
+  projection} the vertices of a one dimensional mesh are successively
+projected on the circle.
 \begin{figure}
 \center
 \includegraphics[width=0.6\textwidth]{fig/boundary_projection.pdf}
-\caption{Boundary projection for a two dimensional mesh. The boundary vertices of the mesh are projected on the circle.}
+\caption{Boundary projection for a two dimensional mesh. The boundary
+  vertices of the mesh are projected on the circle.}
 \label{f:boundary projection}
 \end{figure}
 \begin{figure}
 \center
 \includegraphics[width=0.6\textwidth]{fig/volume_projection.pdf}
-\caption{Volume projection for the one dimensional mesh. All vertices of this mesh are projected on the circle.}
+\caption{Volume projection for the one dimensional mesh. All vertices
+  of this mesh are projected on the circle.}
 \label{f:volume projection}
 \end{figure}
 
-In this section, we give an example for both projection types. As projection we choose the projection to the unit sphere in 3d. In the first example, we start with the three dimensional cube $[-1,1]^3$ and solve the three dimensional version of problem \ref{eq:ellipt} in it. Furthermore, we apply a boundary projection to the unit sphere. In the second example, we set the right hand side $f$ of equation (\ref{eq:ellipt}) to $2x_0$ ($x_0$ is the first component of x), and solve on the two dimensional surface of the sphere. Here, we start with a macro mesh that defines the surface of a cube and apply a volume projection to it. 
+In this section, we give an example for both projection types. As
+projection we choose the projection to the unit sphere in 3d. In the
+first example, we start with the three dimensional cube $[-1,1]^3$ and
+solve the three dimensional version of problem \ref{eq:ellipt} in
+it. Furthermore, we apply a boundary projection to the unit sphere. In
+the second example, we set the right hand side $f$ of equation
+(\ref{eq:ellipt}) to $2x_0$ ($x_0$ is the first component of x), and
+solve on the two dimensional surface of the sphere. Here, we start
+with a macro mesh that defines the surface of a cube and apply a
+volume projection to it.
 
 \begin{table}
 \center
@@ -61,7 +81,8 @@ In this section, we give an example for both projection types. As projection we
 \end{table}
 
 \subsection{Source code}
-First we define the projection by implementing a sub class \verb+BallProject+ of the base class\\ \verb+Projection+.
+First we define the projection by implementing a sub class
+\verb+BallProject+ of the base class\\ \verb+Projection+.
 
 \begin{lstlisting}{}
 class BallProject : public Projection
@@ -91,9 +112,14 @@ protected:
 };
 \end{lstlisting}
 
-First, in the constructor, the base class constructor is called with a projection identifier and the projection type which can be \verb+BOUNDARY_PROJECTION+ or \verb+VOLUME_PROJECTION+. The projection identifier is used to associated a projection instance to projections defined in the macro file. The method \verb+project+ implements the concrete projection of a point $x$ in world coordinates.
+First, in the constructor, the base class constructor is called with a
+projection identifier and the projection type which can be
+\verb+BOUNDARY_PROJECTION+ or \verb+VOLUME_PROJECTION+. The projection
+identifier is used to associated a projection instance to projections
+defined in the macro file. The method \verb+project+ implements the
+concrete projection of a point $x$ in world coordinates.
 
-If we compute on the surface, we redefine the function $f$.  
+If we compute on the surface, we redefine the function $f$.
 
 \begin{lstlisting}{}
 class F : public AbstractFunction<double, WorldVector<double> >
@@ -108,7 +134,11 @@ public:
 };
 \end{lstlisting}
 
-In the main program, we create an instance of \verb+BallProject+ with ID $1$, center $0$ and radius $1$. If we solve in the three dimensional volume of the sphere, the projection type is \verb+BOUNDARY_PROJECTION+, because we project only boundary vertices to the sphere. 
+In the main program, we create an instance of \verb+BallProject+ with
+ID $1$, center $0$ and radius $1$. If we solve in the three
+dimensional volume of the sphere, the projection type is
+\verb+BOUNDARY_PROJECTION+, because we project only boundary vertices
+to the sphere.
 
 \begin{lstlisting}{}
   // ===== create projection =====
@@ -120,7 +150,9 @@ In the main program, we create an instance of \verb+BallProject+ with ID $1$, ce
                   1.0);
 \end{lstlisting}
 
-If we solve on the two dimensional surface of the sphere, the projection type is \verb+VOLUME_PROJECTION+, because all vertices of the mesh are projected.
+If we solve on the two dimensional surface of the sphere, the
+projection type is \verb+VOLUME_PROJECTION+, because all vertices of
+the mesh are projected.
 
 \begin{lstlisting}{}
   // ===== create projection =====
@@ -160,9 +192,14 @@ sphere->output->AMDiS mesh ext: .mesh
 sphere->output->AMDiS data ext: .dat
 \end{lstlisting}
 
-The world dimension is 3, whereas the mesh dimension is set to 2. We use a macro mesh which defines the surface of a cube, defined in \verb+./macro/sphere_macro.3d+, and apply 10 global refinements to it. In this example we do not use adaptivity. Thus, no estimator and no marker is used.
+The world dimension is 3, whereas the mesh dimension is set to 2. We
+use a macro mesh which defines the surface of a cube, defined in
+\verb+./macro/sphere_macro.3d+, and apply 10 global refinements to
+it. In this example we do not use adaptivity. Thus, no estimator and
+no marker is used.
 
-Now, we show the parameter file for the boudary projection case (three dimensional mesh).
+Now, we show the parameter file for the boudary projection case (three
+dimensional mesh).
 
 \begin{lstlisting}{}
 dimension of world:             3
@@ -188,7 +225,8 @@ ball->output->AMDiS mesh ext: .mesh
 ball->output->AMDiS data ext: .dat
 \end{lstlisting}
 
-The macro mesh is a three dimensional cube, defined in \verb+./macro/macro.ball.3d+, and 15 times globally refined.
+The macro mesh is a three dimensional cube, defined in
+\verb+./macro/macro.ball.3d+, and 15 times globally refined.
 
 \subsection{Macro file}
 
@@ -254,7 +292,9 @@ projections:
  1 0 0
 \end{lstlisting}
 
-In the projections block, the projection IDs for each element are listed. There is one entry for each side of each element. Since we use volume projection, here, only the first entry of a line is used.
+In the projections block, the projection IDs for each element are
+listed. There is one entry for each side of each element. Since we use
+volume projection, here, only the first entry of a line is used.
 
 Now, we list the macro file for the three dimensional volume mesh.
 
@@ -308,10 +348,14 @@ projections:
   1    1    0    0
 \end{lstlisting}
 
-Here, we use boundary projections. In the boundary block for each boundary side of an element the projection ID is given.
+Here, we use boundary projections. In the boundary block for each
+boundary side of an element the projection ID is given.
 
 \subsection{Output}
-In Figure \ref{f:sphere}, the solution of the two dimensional problem is shown on a successively refined mesh whose vertices are projected on the sphere. The finer the mesh, the better is the approximation to the sphere. 
+In Figure \ref{f:sphere}, the solution of the two dimensional problem
+is shown on a successively refined mesh whose vertices are projected
+on the sphere. The finer the mesh, the better is the approximation to
+the sphere.
 
 \begin{figure}
 \center
@@ -322,7 +366,11 @@ In Figure \ref{f:sphere}, the solution of the two dimensional problem is shown o
 \label{f:sphere}
 \end{figure}
 
-In Figure \ref{f:sphere solution} (a), the final solution of the two dimensional problem is shown, Figure \ref{f:sphere solution} (b) shows the halfed sphere to demonstrate that the solution is really defined on the sphere. The solution of the three dimensional problem is shown in Figure \ref{f:sphere solution} (c).
+In Figure \ref{f:sphere solution} (a), the final solution of the two
+dimensional problem is shown, Figure \ref{f:sphere solution} (b) shows
+the halfed sphere to demonstrate that the solution is really defined
+on the sphere. The solution of the three dimensional problem is shown
+in Figure \ref{f:sphere solution} (c).
 \begin{figure}
 \center
 \begin{tabular}{ccc}
@@ -331,6 +379,8 @@ In Figure \ref{f:sphere solution} (a), the final solution of the two dimensional
 \includegraphics[width=0.32\textwidth]{fig/ball_half.jpg}\\
 (a)&(b)&(c)
 \end{tabular}
-\caption{(a): Solution of the two dimensional problem on the surface of the sphere, (b): Halfed sphere, (c): Solution of the three dimensional problem (halfed ball).}
+\caption{(a): Solution of the two dimensional problem on the surface
+  of the sphere, (b): Halfed sphere, (c): Solution of the three
+  dimensional problem (halfed ball).}
 \label{f:sphere solution}
 \end{figure}
diff --git a/doc/tutorial/tutorial.pdf b/doc/tutorial/tutorial.pdf
index 3fa0a9c2990bd457d5f2afa242ce031fef647b0c..420f82cb5f025e655bd4c493eb2561b68c866fa3 100644
Binary files a/doc/tutorial/tutorial.pdf and b/doc/tutorial/tutorial.pdf differ
diff --git a/doc/tutorial/vecellipt.tex b/doc/tutorial/vecellipt.tex
index 035c0b603dff07c0e5f44541be1337ffec96c639..e7a950bc6e5164149da0da68c0db7aacf23ab2f9 100644
--- a/doc/tutorial/vecellipt.tex
+++ b/doc/tutorial/vecellipt.tex
@@ -45,7 +45,7 @@ Instead of a scalar problem, we now create and initialize the vector
 valued problem \verb+vecellipt+:
 
 \begin{lstlisting}{}
-  ProblemVec vecellipt("vecellipt");
+  ProblemStat vecellipt("vecellipt");
   vecellipt.initialize(INIT_ALL);
 \end{lstlisting}