diff --git a/doc/tutorial/composite.tex b/doc/tutorial/composite.tex
new file mode 100644
index 0000000000000000000000000000000000000000..d23cbc57edb5ff4fc599efca8a025e8e4a9d0c5d
--- /dev/null
+++ b/doc/tutorial/composite.tex
@@ -0,0 +1,6 @@
+\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
new file mode 100644
index 0000000000000000000000000000000000000000..3dc4c1ba2d2dd0f98217dd5607912d6fd8d97163
--- /dev/null
+++ b/doc/tutorial/couple.tex
@@ -0,0 +1,299 @@
+\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:
+ MEMORY_MANAGED(Identity);
+
+ Identity(int degree) : AbstractFunction<double, double>(degree) {};
+
+ const double& operator()(const double& x) const {
+ static double result;
+ result = x;
+ return result;
+ };
+};
+\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(Operator::MATRIX_OPERATOR,
+ problem2.getFESpace());
+ matrixOperator2.addZeroOrderTerm(NEW Simple_ZOT);
+ problem2.addMatrixOperator(&matrixOperator2);
+
+ Operator rhsOperator2(Operator::VECTOR_OPERATOR, 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
new file mode 100644
index 0000000000000000000000000000000000000000..7c6bbd491fc2c658c3925b0b7efa4dfec3d1536f
--- /dev/null
+++ b/doc/tutorial/ellipt.tex
@@ -0,0 +1,341 @@
+\section{Stationary problem with Dirichlet boundary condition}
+\label{s:ellipt}
+As example for a stationary problem, we choose the Poisson equation
+\begin{eqnarray}
+\label{eq:ellipt}
+-\Delta u &=& f ~~~\mbox{in } \Omega \subset \mathbb{R}^{dim}\\
+u &=& g ~~~ \mbox{on } \partial\Omega
+\end{eqnarray}
+with
+\begin{eqnarray}
+f(x) &=& -\left( 400x^2 - 20 dow \right) e^{-10x^2}\\
+g(x) &=& e^{-10x^2}.
+\end{eqnarray}
+$dim$ is the problem dimension and thus the dimension of the mesh the problem will be discretized on. $dow$ is the dimension of the world the problem lives in. So world coordinates are always real valued vectors of size $dow$. Note that the problem is defined in a dimension independent way. Furthermore, $dow$ has not to be equal to $dim$ as long as $1 \leq dim \leq dow \leq 3$ holds.
+
+Although the implementation described in Section \ref{s:ellipt code} is dimension independent, we focus on the case $dim = dow = 2$ for the rest of this section. The analytical solution on $\Omega = [0,1] \times [0,1]$ is shown in Figure \ref{f:ellipt solution}.
+\begin{figure}
+\center
+\includegraphics[width=0.3\textwidth]{fig/ellipt.jpg}
+\caption{Solution of the Poisson equation on the unit square.}
+\label{f:ellipt solution}
+\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 ellipt.cc} \\
+\hline
+{\bf parameter file} & \tt init/ & \tt ellipt.dat.1d & \tt ellipt.dat.2d & \tt ellipt.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 ellipt.mesh, ellipt.dat} \\
+\hline
+\end{tabular}
+\caption{Files of the {\tt ellipt} example.}
+\label{t:ellipt}
+\end{table}
+
+\subsection{Source code}
+\label{s:ellipt code}
+For this first example, we give the complete source code. But to avoid loosing the overview, we sometimes interrupt the code to give some explaining comment. The first three lines of the application code are:
+
+\begin{lstlisting}{}
+#include "AMDiS.h"
+using namespace std;
+using namespace AMDiS;
+\end{lstlisting}
+
+In the first line, the AMDiS header is included. In line 2 and 3, used namespaces are introduced. \verb+std+ is the C++ standard library namespace, used e.g.\ for the STL classes. AMDiS provides its own namespace \verb+AMDiS+ to avoid potential naming conflicts with other libraries.
+
+Now, the functions $f$ and $g$ will be defined by the classes \verb+F+ and \verb+G+:
+
+\begin{lstlisting}{}
+// ===== function definitions =====
+class G : public AbstractFunction<double, WorldVector<double> >
+{
+public:
+ MEMORY_MANAGED(G);
+
+ double operator()(const WorldVector<double>& x) const
+ {
+ return exp(-10.0 * (x * x));
+ }
+};
+\end{lstlisting}
+
+\verb+G+ is a sub class of the templated class
+\verb+AbstractFunction<R, T>+ that represents a mapping from type
+\verb+T+ to type \verb+R+. Here, we want to define a mapping from
+$\mathbb{R}^{dow}$, implemented by the class
+\verb+WorldVector<double>+, to $\mathbb{R}$, represented by the data
+type \verb+double+. The actual mapping is defined by overloading the
+\verb+operator()+. \verb+x*x+ stands for the scalar product of vector
+\verb+x+ with itself.
+
+Using the macro call \verb+MEMORY_MANAGED(G)+, the class will be
+managed by the memory management of AMDiS. This memory management
+provides memory monitoring to locate memory leaks and block memory
+allocation to accelerate memory access. Objects of a memory managed
+class can now be allocated through the memory manager by the macro
+\verb+NEW+ and deallocated by the macro \verb+DELETE+.
+
+The class \verb+F+ is defined in a similar way:
+
+\begin{lstlisting}{}
+class F : public AbstractFunction<double, WorldVector<double> >
+{
+public:
+ MEMORY_MANAGED(F);
+
+ F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {};
+
+ double operator()(const WorldVector<double>& x) const {
+ int dow = Global::getGeo(WORLD);
+ double r2 = (x * x);
+ double ux = exp(-10.0 * r2);
+ return -(400.0 * r2 - 20.0 * dow) * ux;
+ }
+};
+\end{lstlisting}
+
+\verb+F+ gets the world dimension from the class \verb+Global+ by a call of the static function \verb+getGeo(WORLD)+. The degree handed to the constructor determines the polynomial degree with which the function should be considered in the numerical integration. A higher degree leads to a quadrature of higher order in the assembling process.
+
+Now, we start with the main program:
+
+\begin{lstlisting}{}
+// ===== main program =====
+int main(int argc, char* argv[])
+{
+ FUNCNAME("main");
+
+ // ===== check for init file =====
+ TEST_EXIT(argc == 2)("usage: ellipt initfile\n");
+
+ // ===== init parameters =====
+ Parameters::init(true, argv[1]);
+\end{lstlisting}
+
+The macro \verb+FUNCNAME+ defines the current function name that is used for command line output, e.g. in error messages. The macro \verb+TEST_EXIT+ tests for the condition within the first pair of brackets. If the condition does not hold, an error message given in the second bracket pair is prompted and the program exits. Here the macro is used to check, whether the parameter file was specified by the user as command line argument. If this is the case, the parameters are initialized by \verb+Parameters::init(true, argv[1])+. The first 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:
+
+\begin{lstlisting}{}
+ // ===== create and init the scalar problem =====
+ ProblemScal ellipt("ellipt");
+ ellipt.initialize(INIT_ALL);
+\end{lstlisting}
+
+The name argument of the problem is used to identify parameters in the 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 components can be controlled through the parameter file (see Section \ref{s:ellipt parameter}).
+
+The next steps are the creation of the adaptation loop and the corresponding \verb+AdaptInfo+:
+
+\begin{lstlisting}{}
+ // === create adapt info ===
+ AdaptInfo *adaptInfo = NEW AdaptInfo("ellipt->adapt", 1);
+
+ // === create adapt ===
+ AdaptStationary *adapt = NEW AdaptStationary("ellipt->adapt", &ellipt,
+ adaptInfo);
+\end{lstlisting}
+
+The \verb+AdaptInfo+ object contains information about the current state of the adaptation loop as well as user given parameters concerning the adaptation loop, like desired tolerances or maximal iteration numbers. Using \verb+adaptInfo+, the adaptation loop can be inspected and controlled at runtime.
+Now, a stationary adaptation loop is created, which implements the standard {\it assemble-solve-estimate-adapt} loop. Arguments are the name, again 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 \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(Operator::MATRIX_OPERATOR, ellipt.getFESpace());
+ matrixOperator.addSecondOrderTerm(NEW Laplace_SOT);
+ ellipt.addMatrixOperator(&matrixOperator);
+
+ // ===== create rhs operator =====
+ int degree = ellipt.getFESpace()->getBasisFcts()->getDegree();
+ Operator rhsOperator(Operator::VECTOR_OPERATOR, ellipt.getFESpace());
+ rhsOperator.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree)));
+ ellipt.addVectorOperator(&rhsOperator);
+\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 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.
+
+Finally we start the adaptation loop and afterwards write out the results:
+
+\begin{lstlisting}{}
+ // ===== start adaption loop =====
+ adapt->adapt();
+
+ // ===== write result =====
+ ellipt.writeFiles(adaptInfo, true);
+}
+\end{lstlisting}
+
+The second argument of \verb+writeFiles+ forces the file writer to print out the results. In time dependent problems it can be useful to write the results only every $i$-th timestep. To allow this behavior the second argument has to be \verb+false+.
+
+\subsection{Parameter file}
+\label{s:ellipt parameter}
+The name of the parameter file must be given as command line argument. In the 2d case we call:
+
+\begin{lstlisting}{}
+> ./ellipt init/ellipt.dat.2d
+\end{lstlisting}
+
+In the following, the content of file \verb+init/ellipt.dat.2d+ is described:
+
+\begin{lstlisting}{}
+dimension of world: 2
+
+elliptMesh->macro file name: ./macro/macro.stand.2d
+elliptMesh->global refinements: 0
+\end{lstlisting}
+
+The dimension of the world is 2, the macro mesh with name \verb+elliptMesh+ is defined in file \\\verb+./macro/macro.stand.2d+ (see Section \ref{s:ellipt macro}). The mesh is not globally refined before the adaptation loop. A value of $n$ for \verb+elliptMesh->global refinements+ means $n$ bisections of every macro element. Global refinements before the adaptation loop can be useful to save computation time by starting adaptive computations with a finer mesh.
+
+\begin{lstlisting}{}
+ellipt->mesh: elliptMesh
+ellipt->dim: 2
+ellipt->polynomial degree: 3
+\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 3.
+
+\begin{lstlisting}{}
+ellipt->solver: cg % no bicgstab cg gmres odir ores
+ellipt->solver->max iteration: 1000
+ellipt->solver->tolerance: 1.e-8
+ellipt->solver->left precon: diag % no, diag
+\end{lstlisting}
+
+We use the {\it conjugate gradient method} as iterative solver. The solving process stops after maximal $1000$ iterations or when a tolerance of $10^{-8}$ is reached. Furthermore, we apply diagonal pre-conditioning.
+
+\begin{lstlisting}{}
+ellipt->estimator: residual % residual, recovery
+ellipt->estimator->error norm: 1
+ellipt->estimator->C0: 0.1
+ellipt->estimator->C1: 0.1
+\end{lstlisting}
+
+As error estimator we use the residual method. The used error norm is 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
+\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->max iteration: 100
+ellipt->adapt->refine bisections: 2
+\end{lstlisting}
+
+The adaptation loop stops, when an error tolerance of $10^{-4}$ is reached, or after maximal $100$ iterations. An element that is marked for refinement, is bisected twice within one iteration. Analog elements that are marked for coarsening are coarsened twice per iteration.
+
+\begin{lstlisting}{}
+ellipt->output->filename: output/ellipt
+ellipt->output->AMDiS format: 1
+ellipt->output->AMDiS mesh ext: .mesh
+ellipt->output->AMDiS data ext: .dat
+\end{lstlisting}
+
+The result is written in AMDiS-format to the files \verb+output/ellipt.mesh+ and \verb+output/ellipt.dat+. The first contains the final mesh, the second contains the corresponding solution values.
+
+\subsection{Macro file}
+\label{s:ellipt macro}
+\begin{figure}
+\center
+\includegraphics[width=0.4\textwidth]{fig/ellipt_macro.pdf}
+\caption{Two dimensional macro mesh}
+\label{f:ellipt macro}
+\end{figure}
+In Figure \ref{f:ellipt macro} one can see the macro mesh which is described by the file \verb+macro/macro.stand.2d+.
+First, the dimension of the mesh and of the world are defined:
+
+\begin{lstlisting}{}
+DIM: 2
+DIM_OF_WORLD: 2
+\end{lstlisting}
+
+Then the total number of elements and vertices are given:
+
+\begin{lstlisting}{}
+number of elements: 4
+number of vertices: 5
+\end{lstlisting}
+
+The next block describes the two dimensional coordinates of the five vertices:
+
+\begin{lstlisting}{}
+vertex coordinates:
+ 0.0 0.0
+ 1.0 0.0
+ 1.0 1.0
+ 0.0 1.0
+ 0.5 0.5
+\end{lstlisting}
+
+The first two numbers are interpreted as the coordinates of vertex 0, and so on.
+
+Corresponding to these vertex indices now the four triangles are given:
+
+\begin{lstlisting}{}
+element vertices:
+0 1 4
+1 2 4
+2 3 4
+3 0 4
+\end{lstlisting}
+
+Element 0 consists in the vertices 0, 1 and 4. The numbering is done anticlockwise starting with the vertices of the longest edge.
+
+It follows the definition of boundary conditions:
+
+\begin{lstlisting}{}
+element boundaries:
+0 0 1
+0 0 1
+0 0 1
+0 0 1
+\end{lstlisting}
+
+The first number line means that element 0 has no boundaries at edge 0 and 1, and a boundary with identifier 1 at edge 2. The edge with number $i$ is the edge opposite to vertex number $i$. The boundary identifier 1 corresponds to the identifier 1 we defined within the source code for the Dirichlet boundary. Since all elements of the macro mesh have a Dirichlet boundary at edge 2, the line \verb+0 0 1+ is repeated three times.
+
+The next block defines element neighborships. \verb+-1+ means there is no neighbor at the corresponding edge. A non-negative number determines the index of the neighbor element.
+
+\begin{lstlisting}{}
+element neighbours:
+1 3 -1
+2 0 -1
+3 1 -1
+0 2 -1
+\end{lstlisting}
+
+This block is optional. If it isn't given in the macro file, element neighborships are computed automatically.
+
+\subsection{Output}
+\begin{figure}
+\center
+\begin{tabular}{cc}
+\includegraphics[width=0.3\textwidth]{fig/ellipt_dat.jpg}&
+\includegraphics[width=0.3\textwidth]{fig/ellipt_mesh.jpg}\\
+(a)&(b)
+\end{tabular}
+\caption{(a): Solution after 9 iterations, (b): corresponding mesh}
+\label{f:ellipt solution and mesh}
+\end{figure}
+
+Now, the program is started by the call \verb+./ellipt init/ellipt.dat.2d+. After 9 iterations the tolerance is reached and the files \verb+output/ellipt.mesh+ and \verb+output/ellipt.dat+ are written. In Figure \ref{f:ellipt solution and mesh}(a) the solution is shown and in \ref{f:ellipt solution and mesh}(b) the corresponding mesh. The visualizations was done by the VTK based tool {\bf CrystalClear}.
\ No newline at end of file
diff --git a/doc/tutorial/examples.tex b/doc/tutorial/examples.tex
new file mode 100644
index 0000000000000000000000000000000000000000..6855262e1747e28ae8e6acfd19ad20bdc55fc9b2
--- /dev/null
+++ b/doc/tutorial/examples.tex
@@ -0,0 +1,15 @@
+\chapter{Implementation of example problems}
+\label{c:examples}
+
+\input{ellipt.tex}
+\input{heat.tex}
+\input{vecellipt.tex}
+\input{couple.tex}
+\input{nonlin.tex}
+\input{neumann.tex}
+\input{periodic.tex}
+\input{projections.tex}
+\input{parametric.tex}
+\input{multigrid.tex}
+\input{parallelheat.tex}
+%\input{composite.tex}
diff --git a/doc/tutorial/fig/ball_half.jpg b/doc/tutorial/fig/ball_half.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..3859a4362f55996361e6bb8b5b7007bc177e50a2
Binary files /dev/null and b/doc/tutorial/fig/ball_half.jpg differ
diff --git a/doc/tutorial/fig/boundary_projection.odg b/doc/tutorial/fig/boundary_projection.odg
new file mode 100644
index 0000000000000000000000000000000000000000..bfb51246d99b26781682b89b729542880cd39bd8
Binary files /dev/null and b/doc/tutorial/fig/boundary_projection.odg differ
diff --git a/doc/tutorial/fig/boundary_projection.pdf b/doc/tutorial/fig/boundary_projection.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..bbeefbda38087a730f37b5b0dc6fe95430c596ad
Binary files /dev/null and b/doc/tutorial/fig/boundary_projection.pdf differ
diff --git a/doc/tutorial/fig/coupled_iteration.odg b/doc/tutorial/fig/coupled_iteration.odg
new file mode 100644
index 0000000000000000000000000000000000000000..48dadb789d7e4903e29d0d9376f12697ed0a74a3
Binary files /dev/null and b/doc/tutorial/fig/coupled_iteration.odg differ
diff --git a/doc/tutorial/fig/coupled_iteration.pdf b/doc/tutorial/fig/coupled_iteration.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..b91ec86d4af47087c6b3d50623cc1eb51707c632
Binary files /dev/null and b/doc/tutorial/fig/coupled_iteration.pdf differ
diff --git a/doc/tutorial/fig/ellipt.jpg b/doc/tutorial/fig/ellipt.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..80a29f723085a680a3a2f55aa004803edf44e84f
Binary files /dev/null and b/doc/tutorial/fig/ellipt.jpg differ
diff --git a/doc/tutorial/fig/ellipt_dat.jpg b/doc/tutorial/fig/ellipt_dat.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..21e87f97913eb6cc5d06e404ba7abfa5c4245849
Binary files /dev/null and b/doc/tutorial/fig/ellipt_dat.jpg differ
diff --git a/doc/tutorial/fig/ellipt_macro.odg b/doc/tutorial/fig/ellipt_macro.odg
new file mode 100644
index 0000000000000000000000000000000000000000..e80b9ade14b02e84bda1757eb4cb0f49bf325c0e
Binary files /dev/null and b/doc/tutorial/fig/ellipt_macro.odg differ
diff --git a/doc/tutorial/fig/ellipt_macro.pdf b/doc/tutorial/fig/ellipt_macro.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..f8ede754910fef0d5aa92b34972b689346b87112
Binary files /dev/null and b/doc/tutorial/fig/ellipt_macro.pdf differ
diff --git a/doc/tutorial/fig/ellipt_mesh.jpg b/doc/tutorial/fig/ellipt_mesh.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..5f6e3efb88aeaa34731850512ee4b47867f10a10
Binary files /dev/null and b/doc/tutorial/fig/ellipt_mesh.jpg differ
diff --git a/doc/tutorial/fig/global.pdf b/doc/tutorial/fig/global.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..bacd9fcb782520547f5ae8eeb70f339c96227bef
Binary files /dev/null and b/doc/tutorial/fig/global.pdf differ
diff --git a/doc/tutorial/fig/heat1.jpg b/doc/tutorial/fig/heat1.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..225cd7672e49775514236e6d04b708434a8488f1
Binary files /dev/null and b/doc/tutorial/fig/heat1.jpg differ
diff --git a/doc/tutorial/fig/heat2.jpg b/doc/tutorial/fig/heat2.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..fd8e704665f581118d6ad33e8262c9ff2adbba4f
Binary files /dev/null and b/doc/tutorial/fig/heat2.jpg differ
diff --git a/doc/tutorial/fig/heat3.jpg b/doc/tutorial/fig/heat3.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..3cd9b2b88d93a7824bab237b21481673d9584d5d
Binary files /dev/null and b/doc/tutorial/fig/heat3.jpg differ
diff --git a/doc/tutorial/fig/heat_uml.odg b/doc/tutorial/fig/heat_uml.odg
new file mode 100644
index 0000000000000000000000000000000000000000..ab6ee88127193f5c8d979ea2b65fd3a35141ed60
Binary files /dev/null and b/doc/tutorial/fig/heat_uml.odg differ
diff --git a/doc/tutorial/fig/heat_uml.pdf b/doc/tutorial/fig/heat_uml.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..4dcfa0a10ec6ac103550656eaf357dd1fab4352e
Binary files /dev/null and b/doc/tutorial/fig/heat_uml.pdf differ
diff --git a/doc/tutorial/fig/local1.pdf b/doc/tutorial/fig/local1.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..0bf93765c09a583547d0703fdb4d81106a4ae36b
Binary files /dev/null and b/doc/tutorial/fig/local1.pdf differ
diff --git a/doc/tutorial/fig/local2.pdf b/doc/tutorial/fig/local2.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..21abcdf899710472f1341a10ca9cadc6c5439904
Binary files /dev/null and b/doc/tutorial/fig/local2.pdf differ
diff --git a/doc/tutorial/fig/macro.pdf b/doc/tutorial/fig/macro.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..e7668135917357e8e9b91922d2b14eb9abeb76b5
--- /dev/null
+++ b/doc/tutorial/fig/macro.pdf
@@ -0,0 +1,58 @@
+%PDF-1.3
+%Çì¢
+6 0 obj
+<</Length 7 0 R/Filter /FlateDecode>>
+stream
+xœM‹A
+€0 ïyE^›45õêDD°õà÷Uhìa—aÖ£{“{J°C3(.'0® d¢ð#6Ó"?@¢W
+µÆÏ‘¨B¾ÅÏb2_^ÛOʨ:åõJ=ÜK5 …endstream
+endobj
+7 0 obj
+98
+endobj
+8 0 obj
+<</R4
+4 0 R>>
+endobj
+5 0 obj
+<</Type/Page/MediaBox [0 0 284 284]
+/Parent 3 0 R
+/Resources<</ProcSet[/PDF]
+/ExtGState 8 0 R
+>>
+/Contents 6 0 R
+>>
+endobj
+3 0 obj
+<< /Type /Pages /Kids [
+5 0 R
+] /Count 1
+>>
+endobj
+1 0 obj
+<</Type /Catalog /Pages 3 0 R
+>>
+endobj
+4 0 obj
+<</Type/ExtGState/Name/R4/TR/Identity/OPM 1/SM 0.02>>
+endobj
+2 0 obj
+<</Producer (GNU Ghostscript 6.53)
+>>endobj
+xref
+0 9
+0000000000 65535 f
+0000000420 00000 n
+0000000537 00000 n
+0000000361 00000 n
+0000000468 00000 n
+0000000230 00000 n
+0000000015 00000 n
+0000000183 00000 n
+0000000201 00000 n
+trailer
+<< /Size 9 /Root 1 0 R /Info 2 0 R
+>>
+startxref
+589
+%%EOF
diff --git a/doc/tutorial/fig/neumann_solution.jpg b/doc/tutorial/fig/neumann_solution.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..9a078b9b7833766538906ebcc47fbe43a5a2d66f
Binary files /dev/null and b/doc/tutorial/fig/neumann_solution.jpg differ
diff --git a/doc/tutorial/fig/nonlin_iteration.odg b/doc/tutorial/fig/nonlin_iteration.odg
new file mode 100644
index 0000000000000000000000000000000000000000..516148ee00b000344d6b5023da3496acbba74caa
Binary files /dev/null and b/doc/tutorial/fig/nonlin_iteration.odg differ
diff --git a/doc/tutorial/fig/nonlin_iteration.pdf b/doc/tutorial/fig/nonlin_iteration.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..9d0996f120c0a45b1d264ebda865d983a00d454a
Binary files /dev/null and b/doc/tutorial/fig/nonlin_iteration.pdf differ
diff --git a/doc/tutorial/fig/nonlin_mesh.jpg b/doc/tutorial/fig/nonlin_mesh.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..a47d5108568dd31f807b9ef77ec57a0d5b692be8
Binary files /dev/null and b/doc/tutorial/fig/nonlin_mesh.jpg differ
diff --git a/doc/tutorial/fig/nonlin_operators.odg b/doc/tutorial/fig/nonlin_operators.odg
new file mode 100644
index 0000000000000000000000000000000000000000..ea587a48f3132f5698727f4796d29f755c7dd1c2
Binary files /dev/null and b/doc/tutorial/fig/nonlin_operators.odg differ
diff --git a/doc/tutorial/fig/nonlin_operators.pdf b/doc/tutorial/fig/nonlin_operators.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..349c406c6c194b11e103dc9edda3df925db36a2c
Binary files /dev/null and b/doc/tutorial/fig/nonlin_operators.pdf differ
diff --git a/doc/tutorial/fig/nonlin_solution.jpg b/doc/tutorial/fig/nonlin_solution.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..ddeab59a85ca3afc0547546e6cb3e7cdd7ec6621
Binary files /dev/null and b/doc/tutorial/fig/nonlin_solution.jpg differ
diff --git a/doc/tutorial/fig/parallelheat_01.jpg b/doc/tutorial/fig/parallelheat_01.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..ea0c519abe4e3562469afbff546587b2d4a23042
Binary files /dev/null and b/doc/tutorial/fig/parallelheat_01.jpg differ
diff --git a/doc/tutorial/fig/parallelheat_05.jpg b/doc/tutorial/fig/parallelheat_05.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..7b133d966b56a076f4b5855b686378aa094da0cf
Binary files /dev/null and b/doc/tutorial/fig/parallelheat_05.jpg differ
diff --git a/doc/tutorial/fig/parallelheat_09.jpg b/doc/tutorial/fig/parallelheat_09.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..669818607baa331e9dc0f9431bdaa3a344311dcc
Binary files /dev/null and b/doc/tutorial/fig/parallelheat_09.jpg differ
diff --git a/doc/tutorial/fig/periodic_boundaries.odg b/doc/tutorial/fig/periodic_boundaries.odg
new file mode 100644
index 0000000000000000000000000000000000000000..69991673a6324f00ac814747e944a1ab338dcc2a
Binary files /dev/null and b/doc/tutorial/fig/periodic_boundaries.odg differ
diff --git a/doc/tutorial/fig/periodic_boundaries.pdf b/doc/tutorial/fig/periodic_boundaries.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..882c39c30d4b26a58f6924e631986f2cb9554a4f
Binary files /dev/null and b/doc/tutorial/fig/periodic_boundaries.pdf differ
diff --git a/doc/tutorial/fig/periodic_identification.odg b/doc/tutorial/fig/periodic_identification.odg
new file mode 100644
index 0000000000000000000000000000000000000000..4ad4c98e7598c57fe24d1349e9b26c010fdf70e5
Binary files /dev/null and b/doc/tutorial/fig/periodic_identification.odg differ
diff --git a/doc/tutorial/fig/periodic_identification.pdf b/doc/tutorial/fig/periodic_identification.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..b4002be0f05f5ceed6f3e45009fcc3bc51b1408c
Binary files /dev/null and b/doc/tutorial/fig/periodic_identification.pdf differ
diff --git a/doc/tutorial/fig/periodic_iso.jpg b/doc/tutorial/fig/periodic_iso.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..9ba61156587886568d3c2ec4e5a6b1f23dcc5eb8
Binary files /dev/null and b/doc/tutorial/fig/periodic_iso.jpg differ
diff --git a/doc/tutorial/fig/periodic_macro.odg b/doc/tutorial/fig/periodic_macro.odg
new file mode 100644
index 0000000000000000000000000000000000000000..c5fa88ea13e9cb58c0d8092ff05e5da8d96a9841
Binary files /dev/null and b/doc/tutorial/fig/periodic_macro.odg differ
diff --git a/doc/tutorial/fig/periodic_macro.pdf b/doc/tutorial/fig/periodic_macro.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..2bffa05b7c236384b7ca80ef71337bd7897a803d
Binary files /dev/null and b/doc/tutorial/fig/periodic_macro.pdf differ
diff --git a/doc/tutorial/fig/periodic_solution.jpg b/doc/tutorial/fig/periodic_solution.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..7492b18fc05de86a97402ca16f041598b746eaa0
Binary files /dev/null and b/doc/tutorial/fig/periodic_solution.jpg differ
diff --git a/doc/tutorial/fig/rotated_torus.jpg b/doc/tutorial/fig/rotated_torus.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..67ce92d128a0548015c30b17f3d43167d763f6d3
Binary files /dev/null and b/doc/tutorial/fig/rotated_torus.jpg differ
diff --git a/doc/tutorial/fig/sphere0.jpg b/doc/tutorial/fig/sphere0.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..f55b3a2cde4c4bd53ce0e360eea69e6d5167e68f
Binary files /dev/null and b/doc/tutorial/fig/sphere0.jpg differ
diff --git a/doc/tutorial/fig/sphere2.jpg b/doc/tutorial/fig/sphere2.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..ce1e3e35d8e9dc8d68651c0d26788f4d43f281aa
Binary files /dev/null and b/doc/tutorial/fig/sphere2.jpg differ
diff --git a/doc/tutorial/fig/sphere4.jpg b/doc/tutorial/fig/sphere4.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..b38bd9a5189cc7d26e574f1f673644e9e71fff77
Binary files /dev/null and b/doc/tutorial/fig/sphere4.jpg differ
diff --git a/doc/tutorial/fig/sphere_half.jpg b/doc/tutorial/fig/sphere_half.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..4c3113128f754e379a31ceedffd4d3eb110212e5
Binary files /dev/null and b/doc/tutorial/fig/sphere_half.jpg differ
diff --git a/doc/tutorial/fig/sphere_solution.jpg b/doc/tutorial/fig/sphere_solution.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..1ba721942f3d10c80a60f9a50686cf7fb4634545
Binary files /dev/null and b/doc/tutorial/fig/sphere_solution.jpg differ
diff --git a/doc/tutorial/fig/stationary_loop.odg b/doc/tutorial/fig/stationary_loop.odg
new file mode 100644
index 0000000000000000000000000000000000000000..a51c9dc1c224ee132b6909afd6df913eba4d04c4
Binary files /dev/null and b/doc/tutorial/fig/stationary_loop.odg differ
diff --git a/doc/tutorial/fig/stationary_loop.pdf b/doc/tutorial/fig/stationary_loop.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..e3dd19d9f21958fb3da79a3600f9ba2ac577c9cd
Binary files /dev/null and b/doc/tutorial/fig/stationary_loop.pdf differ
diff --git a/doc/tutorial/fig/theta.odg b/doc/tutorial/fig/theta.odg
new file mode 100644
index 0000000000000000000000000000000000000000..0ed726ac3d369f141d52399ada8029d23c76f9c5
Binary files /dev/null and b/doc/tutorial/fig/theta.odg differ
diff --git a/doc/tutorial/fig/theta.pdf b/doc/tutorial/fig/theta.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..3ffdbbc80765bc5a300dc1e0b05c1dcd3933ea79
Binary files /dev/null and b/doc/tutorial/fig/theta.pdf differ
diff --git a/doc/tutorial/fig/threelevels.pdf b/doc/tutorial/fig/threelevels.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..bd467bbfc5710db6f08dbf554f3e780e66c8efea
Binary files /dev/null and b/doc/tutorial/fig/threelevels.pdf differ
diff --git a/doc/tutorial/fig/torus.odg b/doc/tutorial/fig/torus.odg
new file mode 100644
index 0000000000000000000000000000000000000000..cb1be38985ed6a0a8e18e998162c8d1430b00b69
Binary files /dev/null and b/doc/tutorial/fig/torus.odg differ
diff --git a/doc/tutorial/fig/torus.pdf b/doc/tutorial/fig/torus.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..fe4ec20ff1e3b2b35786918341c09780307cec7b
Binary files /dev/null and b/doc/tutorial/fig/torus.pdf differ
diff --git a/doc/tutorial/fig/torus_macro.jpg b/doc/tutorial/fig/torus_macro.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..f76d8803e7f789be3fc368e84ac6a0a5c47924b0
Binary files /dev/null and b/doc/tutorial/fig/torus_macro.jpg differ
diff --git a/doc/tutorial/fig/vecellipt0.jpg b/doc/tutorial/fig/vecellipt0.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..5725fd0744718d7c164e3a8fd536c5b2d6e8604b
Binary files /dev/null and b/doc/tutorial/fig/vecellipt0.jpg differ
diff --git a/doc/tutorial/fig/vecellipt1.jpg b/doc/tutorial/fig/vecellipt1.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..763baed2dc2980188410e5c3336a332b93d6ab6d
Binary files /dev/null and b/doc/tutorial/fig/vecellipt1.jpg differ
diff --git a/doc/tutorial/fig/volume_projection.odg b/doc/tutorial/fig/volume_projection.odg
new file mode 100644
index 0000000000000000000000000000000000000000..94da8fefe008713ce710712f05e1a5f89b524b3e
Binary files /dev/null and b/doc/tutorial/fig/volume_projection.odg differ
diff --git a/doc/tutorial/fig/volume_projection.pdf b/doc/tutorial/fig/volume_projection.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..dac666411194acd4f6008f127f62be5f1a619576
Binary files /dev/null and b/doc/tutorial/fig/volume_projection.pdf differ
diff --git a/doc/tutorial/heat.tex b/doc/tutorial/heat.tex
new file mode 100644
index 0000000000000000000000000000000000000000..fdf5cc250e8747a2d6c15d181686878f99c09908
--- /dev/null
+++ b/doc/tutorial/heat.tex
@@ -0,0 +1,376 @@
+\section{Time dependent problem}
+\label{s:heat}
+
+This is an example for a time dependent scalar problem. The problem is described by the heat equation
+\begin{eqnarray}
+\label{eq:heat}
+\partial_t u - \Delta u &=& f ~~~~~~~ \mbox{in } \Omega \subset \mathbb{R}^{dim}\times\left({t^{begin}},t^{end}\right)\\
+u &=& g ~~~~~~~ \mbox{on } \partial\Omega \times\left({t^{begin}},t^{end}\right)\\
+u &=& u_0 ~~~~~ \mbox{on } \Omega \times \left({t^{begin}}\right).
+\end{eqnarray}
+
+We solve the problem in the time interval $\left(t^{begin},t^{end}\right)$ with Dirichlet boundary conditions on $\partial\Omega$. The problem is constructed, such that the exact solution is $u(x,t) = \mbox{sin}(\pi t)e^{-10x^2}$. So we set
+\begin{eqnarray}
+f(x,t) &=& \pi \mbox{cos}(\pi t)e^{-10x^2}-\left(400x^2-20dow\right)\mbox{sin}(\pi t)e^{-10x^2}\\
+g(x,t) &=& \mbox{sin}(\pi t)e^{-10x^2} \\
+u_0(x) &=& \mbox{sin}(\pi t^{begin})e^{-10x^2}.
+\end{eqnarray}
+
+We use a variable time discretization scheme. Equation (\ref{eq:heat}) is approximated by
+\begin{equation}
+\frac{u^{new} - u^{old}}{\tau} - \left( \theta\Delta u^{new} + (1-\theta)\Delta u^{old} \right) = f(\cdot, t^{old}+\theta\tau).
+\end{equation}
+$\tau = t^{new} - t^{old}$ is the timestep size between the old and the new problem time. $u^{new}$ is the (searched) solution at $t=t^{new}$. $u^{old}$ is the solution at $t=t^{old}$, which is already known from the last timestep. The parameter $\theta$ determines the implicit and explicit treatment of $\Delta u$. For $\theta = 0$ we have the forward explicit Euler scheme, for $\theta = 1$ the backward implicit Euler scheme. $\theta = 0.5$ results in the Crank-Nicholson scheme. If we bring all terms that depend on $u^{old}$ to the right hand side, the equation reads
+\begin{equation}
+\frac{u^{new}}{\tau} - \theta\Delta u^{new} = \frac{u^{old}}{\tau} + (1-\theta)\Delta u^{old} + f(\cdot, t^{old}+\theta\tau).
+\end{equation}
+
+%\begin{figure}[h]
+%\center
+%\includegraphics[width=0.5\textwidth]{fig/theta.pdf}
+%\caption{Time discretization with the $\theta$ scheme.}
+%\label{f:theta}
+%\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 heat.cc} \\
+\hline
+{\bf parameter file} & \tt init/ & \tt heat.dat.1d & \tt heat.dat.2d & \tt heat.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 heat\_<t>.mesh, heat\_<t>.dat} \\
+\hline
+\end{tabular}
+\caption{Files of the {\tt heat} example. In the output file names, {\tt <t>} is replaced by the time.}
+\end{table}
+
+\subsection{Source code}
+Now, we describe the crucial parts of the source code. First, the functions $f$ and $g$ are defined. In contrast to the ellipt example, the functions now are time dependent. This is implemented by deriving the function classes also from class \verb+TimedObject+. This class provides a pointer to the current time, as well as corresponding setting and getting methods. The usage of a pointer to a real value allows to manage the current time in one location. All objects that deal with the same time, point to the same value. In our example, $f$ is evaluated at $t=t^{old}+\theta\tau$, while $g$ (the Dirichlet boundary function for $u^{new}$) is evaluated at $t=t^{new}$. Function $g$ is implemented as follows:
+
+\begin{lstlisting}{}
+class G : public AbstractFunction<double, WorldVector<double> >,
+ public TimedObject
+{
+public:
+ MEMORY_MANAGED(G);
+
+ double operator()(const WorldVector<double>& x) const
+ {
+ return sin(M_PI * (*timePtr)) * exp(-10.0 *(x * x));
+ }
+};
+\end{lstlisting}
+
+The variable $timePtr$ is a base class member of \verb+TimedObject+. This pointer has to be set once before $g$ is evaluated the first time. Implementation of function $f$ is done in the same way.
+
+\begin{figure}
+\center
+\includegraphics[width=0.5\textwidth]{fig/heat_uml.pdf}
+\caption{UML diagram for class Heat.}
+\label{f:heat_uml}
+\end{figure}
+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:
+\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.
+\item \verb+Heat+ knows another implementation of \verb+ProblemStatBase+: This other implementation represents a stationary problem which is solved within each timestep.\end{description}
+
+The first lines of class \verb+Heat+ are:
+
+\begin{lstlisting}{}
+class Heat : public ProblemInstatScal
+{
+public:
+ MEMORY_MANAGED(Heat);
+
+ Heat(ProblemScal *heatSpace)
+ : ProblemInstatScal("heat", heatSpace)
+ {
+ theta = -1.0;
+ GET_PARAMETER(0, name + "->theta", "%f", &theta);
+ TEST_EXIT(theta >= 0)("theta not set!\n");
+ theta1 = theta - 1;
+ };
+\end{lstlisting}
+
+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, $\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 $\theta$-scheme.
+
+The next lines show the implementation of the time interface.
+
+\begin{lstlisting}{}
+ void setTime(AdaptInfo *adaptInfo) {
+ rhsTime =
+ adaptInfo->getTime()-
+ (1-theta)*adaptInfo->getTimestep();
+ boundaryTime = adaptInfo->getTime();
+ tau1 = 1.0 / adaptInfo->getTimestep();
+ };
+
+ void closeTimestep(AdaptInfo *adaptInfo) {
+ ProblemInstatScal::closeTimestep(adaptInfo);
+ WAIT;
+ };
+\end{lstlisting}
+
+The method \verb+setTime+ is called by the adaptation loop to inform the problem about the current time. The right hand side function $f$ will be evaluated at $t^{old}+\theta\tau = t^{new} - (1-\theta)\tau$, the Dirichlet boundary function $g$ at $t^{new}$. $t^{new}$ is the current time, $\tau$ is the current timestep, both set by the adaptation loop and stored in \verb+adaptInfo+. \verb+tau1+ stores the value of $\frac{1}{\tau}$, which is used later as factor for the zero order time discretization terms.
+
+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 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 until the \verb+return+ key is pressed by the user, if the corresponding entry in the parameter file is set to 1. The macro \verb+WAIT_REALLY+ would wait, independent of parameter settings. If \verb+closeTimestep+ wouldn't be overloaded here, the default implementation without the \verb+WAIT+ statement would be called after each timestep.
+
+Now, the implementation of the \verb+ProblemStatBase+ interface begins. As mentioned above, the instationary problem plays the role of the initial problem by implementing this interface.
+
+\begin{lstlisting}{}
+ void solve(AdaptInfo *adaptInfo)
+ {
+ problemStat->getSolution()->interpol(exactSolution);
+ };
+
+ void estimate(AdaptInfo *adaptInfo)
+ {
+ double errMax, errSum;
+ errSum = Error<double>::L2Err(*exactSolution,
+ *(problemStat->getSolution()),
+ 0, &errMax, false);
+ adaptInfo->setEstSum(errSum, 0);
+ adaptInfo->setEstMax(errMax, 0);
+ };
+\end{lstlisting}
+
+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 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 estimate step, the L2 error is computed. The maximal element error and the sum over all element errors are stored in \verb+adaptInfo+. To make the exact solution known to the problem, we need a setting function:
+
+\begin{lstlisting}{}
+ void setExactSolution(AbstractFunction<double, WorldVector<double> > *fct)
+ {
+ exactSolution = fct;
+ }
+\end{lstlisting}
+
+Now, we define some getting functions and the private member variables:
+
+\begin{lstlisting}{}
+ double *getThetaPtr() { return θ };
+ double *getTheta1Ptr() { return &theta1; };
+ double *getTau1Ptr() { return &tau1; };
+ double *getRHSTimePtr() { return &rhsTime; };
+ double *getBoundaryTimePtr() { return &boundaryTime; };
+
+private:
+ double theta;
+ double theta1;
+ double tau1;
+ double rhsTime;
+ double boundaryTime;
+ AbstractFunction<double, WorldVector<double> > *exactSolution;
+};
+\end{lstlisting}
+
+The definition of class \verb+Heat+ is now finished. In the following, the main program is described.
+
+\begin{lstlisting}{}
+int main(int argc, char** argv)
+{
+ // ===== check for init file =====
+ TEST_EXIT(argc == 2)("usage: heat initfile\n");
+
+ // ===== init parameters =====
+ Parameters::init(false, argv[1]);
+
+ // ===== create and init stationary problem =====
+ ProblemScal *heatSpace = NEW ProblemScal("heat->space");
+ heatSpace->initialize(INIT_ALL);
+
+ // ===== create instationary problem =====
+ Heat *heat = new Heat(heatSpace);
+ heat->initialize(INIT_ALL);
+\end{lstlisting}
+
+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+Heat+ we defined above. \verb+heatSpace+ is given to \verb+heat+ as its stationary problem.
+
+The next step is the creation of the needed \verb+AdaptInfo+ objects and of the instationary adaptation loop:
+
+\begin{lstlisting}{}
+ // create adapt info for heat
+ AdaptInfo *adaptInfo = NEW AdaptInfo("heat->adapt");
+
+ // create initial adapt info
+ AdaptInfo *adaptInfoInitial = NEW AdaptInfo("heat->initial->adapt");
+
+ // create instationary adapt
+ AdaptInstationary *adaptInstat = NEW AdaptInstationary("heat->adapt",
+ heatSpace,
+ adaptInfo,
+ heat,
+ adaptInfoInitial);
+\end{lstlisting}
+
+The object \verb+heatSpace+ is handed as \verb+ProblemIterationInterface+ (implemented by class\\ \verb+ProblemScal+) to the adaptation loop. \verb+heat+ is interpreted as \verb+ProblemTimeInterface+ (implemented by class \verb+ProblemInstatScal+).
+
+The definitions of functions $f$ and $g$ are:
+
+\begin{lstlisting}{}
+ // ===== create boundary functions =====
+ G *boundaryFct = NEW G;
+ boundaryFct->setTimePtr(heat->getBoundaryTimePtr());
+ heat->setExactSolution(boundaryFct);
+
+ heatSpace->addDirichletBC(1, boundaryFct);
+
+ // ===== create rhs functions =====
+ int degree = heatSpace->getFESpace()->getBasisFcts()->getDegree();
+ F *rhsFct = NEW F(degree);
+ rhsFct->setTimePtr(heat->getRHSTimePtr());
+\end{lstlisting}
+
+The functions interpreted as \verb+TimedObject+s are linked with the corresponding time pointers by \verb+setTimePtr+. The boundary function is handed to \verb+heat+ as exact solution and as Dirichlet boundary function with identifier $1$ to \verb+heatSpace+.
+
+Now, we define the operators:
+
+\begin{lstlisting}{}
+ // ===== create operators =====
+ double one = 1.0;
+ double zero = 0.0;
+
+ // create laplace
+ Operator *A = NEW Operator(Operator::MATRIX_OPERATOR |
+ Operator::VECTOR_OPERATOR,
+ heatSpace->getFESpace());
+
+ A->addSecondOrderTerm(new Laplace_SOT);
+
+ A->setUhOld(heat->getOldSolution());
+
+ if((*(heat->getThetaPtr())) != 0.0)
+ heatSpace->addMatrixOperator(A, heat->getThetaPtr(), &one);
+
+ if((*(heat->getTheta1Ptr())) != 0.0)
+ heatSpace->addVectorOperator(A, heat->getTheta1Ptr(), &zero);
+\end{lstlisting}
+
+Operator \verb+A+ represents $-\Delta u$. It is used as matrix operator on the left hand side with factor $\theta$ and as vector operator on the right hand side with factor $-(1-\theta) = \theta - 1$. These assemble factors are the second arguments of \verb+addMatrixOperator+ and \verb+addVectorOperator+. The third argument is the factor used for estimation. In this example, the estimator will consider the operator only on the left hand side with factor $1$. On the right hand side the operator is applied to the solution of the last timestep. So the old solution is handed to the operator by \verb+setUhOld+.
+
+\begin{lstlisting}{}
+ // create zero order operator
+ Operator *C = NEW Operator(Operator::MATRIX_OPERATOR |
+ Operator::VECTOR_OPERATOR,
+ heatSpace->getFESpace());
+
+ C->addZeroOrderTerm(NEW Simple_ZOT);
+
+ C->setUhOld(heat->getOldSolution());
+
+ heatSpace->addMatrixOperator(C, heat->getTau1Ptr(), heat->getTau1Ptr());
+ heatSpace->addVectorOperator(C, heat->getTau1Ptr(), heat->getTau1Ptr());
+\end{lstlisting}
+
+The \verb+Simple_ZOT+ of operator \verb+C+ represents the zero order terms for the time discretization. On both sides of the equation $u$ is added with $\frac{1}{\tau}$ as assemble factor and as estimate factor.
+
+Finally, the operator for the right hand side function $f$ is added and the adaptation loop is started:
+
+\begin{lstlisting}{}
+ // create RHS operator
+ Operator *F = NEW Operator(Operator::VECTOR_OPERATOR,
+ heatSpace->getFESpace());
+
+ F->addZeroOrderTerm(NEW CoordsAtQP_ZOT(rhsFct));
+
+ heatSpace->addVectorOperator(F);
+
+ // ===== start adaption loop =====
+ adaptInstat->adapt();
+}
+\end{lstlisting}
+
+\verb+CoordsAtQP_ZOT+ is a zero order term that evaluates a given function $fct$ at all needed quadrature points. At the left hand side, it would represent the term $fct(x, t)\cdot u$, on the right hand side, just $fct(x, t)$. Note that the old solution isn't given to the operator here. Otherwise the term would represent $fct(x, t) \cdot u^{old}$ on the right hand side.
+
+\subsection{Parameter file}
+\label{s:heat parameter}
+
+In this section, we show only the relevant parts of the parameter file \verb+heat.dat.2d+.
+
+First the parameter $\theta$ for the time discretization is defined:
+
+\begin{lstlisting}{}
+heat->theta: 1.0
+\end{lstlisting}
+
+Then we define the initial timestep and the time interval:
+
+\begin{lstlisting}{}
+heat->adapt->timestep: 0.1
+heat->adapt->start time: 0.0
+heat->adapt->end time: 1.0
+\end{lstlisting}
+
+Now, tolerances are determined:
+
+\begin{lstlisting}{}
+heat->adapt->tolerance: 0.01
+heat->adapt->rel space error: 0.5
+heat->adapt->rel time error: 0.5
+heat->adapt->time theta 1: 1.0
+heat->adapt->time theta 2: 0.3
+\end{lstlisting}
+
+The total tolerance is divided in a space tolerance and a time tolerance. The space tolerance is the maximal allowed space error, given by the product of \verb+tolerance+ and \verb+rel space error+. It is reached by adaptive mesh refinements. The time tolerance is the maximal allowed error, due to the timestep size. It is given by the product of \verb+tolerance+ and \verb+rel time error+ and \verb+time theta 1+. It is relevant, only if an implicit time strategy with adaptive timestep size is used. The parameter \verb+time theta 2+ is used to enlarge the timestep, if the estimated time error falls beneath a given threshold.
+
+\begin{lstlisting}{}
+heat->adapt->strategy: 1
+heat->adapt->time delta 1: 0.7071
+heat->adapt->time delta 2: 1.4142
+\end{lstlisting}
+
+If \verb+strategy+ is $0$, an explicit time strategy with fixed timestep size is used. A value of $1$ stands for the implicit strategy. The time tolerance is reached by successively multiplying the timestep with \verb+time delta 1+. If the estimated timestep error is smaller than the product of \verb+tolerance+ and \verb+rel time error+ and \verb+time theta 2+ at the end of a timestep, the timestep size is multiplied by \verb+time delta 2+.
+
+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
+\end{lstlisting}
+
+Now, the output behavior is determined:
+
+\begin{lstlisting}{}
+heat->space->output->filename: output/heat
+
+heat->space->output->AMDiS format: 1
+heat->space->output->AMDiS mesh ext: .mesh
+heat->space->output->AMDiS data ext: .dat
+
+heat->space->output->write every i-th timestep: 10
+
+heat->space->output->append index: 1
+heat->space->output->index length: 6
+heat->space->output->index decimals: 3
+\end{lstlisting}
+
+In this example, all output filenames start with prefix \verb+output/heat+ and end with the extensions \verb+.mesh+ and \verb+.dat+. Output is written after every $10$th timestep. The time of the single solution is added after the filename prefix with 6 letters, three of them are decimals. The solution for $t=0$ e.g. would be written in the files \verb+output/heat00.000.mesh+ and \verb+output/heat00.000.dat+.
+
+Finally, we set parameter \verb+WAIT+ to $1$. So each call of the macro \verb+WAIT+ in the application will lead to an interruption of the program, until the \verb+return+ key is pressed.
+
+\begin{lstlisting}{}
+WAIT: 1
+\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}
+As mentioned above, the output files look like \verb+output/heat00.000.mesh+ and \verb+output/heat00.000.dat+. Depending on the corresponding value in the parameter file only the solution after every $i$-th timestep is written. In Figure \ref{f:heat}, the solution at three timesteps is visualized.
+
+\begin{figure}
+\center
+\begin{tabular}{ccc}
+\includegraphics[width=0.3\textwidth]{fig/heat1.jpg}&
+\includegraphics[width=0.3\textwidth]{fig/heat2.jpg}&
+\includegraphics[width=0.3\textwidth]{fig/heat3.jpg}\\
+$t \approx 0.2 $ & $t \approx 0.47 $ & $t \approx 0.81$
+\end{tabular}
+\caption{The solution at three different timesteps.}
+\label{f:heat}
+\end{figure}
diff --git a/doc/tutorial/installation.tex b/doc/tutorial/installation.tex
new file mode 100644
index 0000000000000000000000000000000000000000..28ebc3a53c2af05797fbe3caa143d2f03b650752
--- /dev/null
+++ b/doc/tutorial/installation.tex
@@ -0,0 +1,111 @@
+\chapter{Installation}
+\label{c:amdis installation}
+
+\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.tar.gz}:
+ \begin{itemize}
+ \item {\tt > gunzip AMDiS.tar.gz}
+ \item {\tt > tar xvf AMDiS.tar}
+ \end{itemize}
+ \item Checking out a CVS project:
+ \begin{itemize}
+ \item {\tt > export CVSROOT=<cvsroot>}
+ \item {\tt > cvs checkout AMDiS}
+ \end{itemize}
+ \item Checking out a CVS project:
+ \begin{itemize}
+ \item {\tt > svn checkout file://<SVN-REPOSITORY-PATH>/AMDiS}
+ \end{itemize}
+ \end{itemize}
+
+\item Change into the AMDiS directory:
+\begin{verbatim}
+> cd 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+--with-mpi=<MPI-DIR>+: MPI installation path. Used for parallelization support.
+\item \verb+--with-parmetis=<PARMETIS-DIR>+: ParMETIS installation path. ParMETIS is a parallel graph partitioning library, see \\\verb+http://glaros.dtc.umn.edu/gkhome/metis/parmetis/overview+. Used for parallelization support.
+\end{description}
+If AMDiS should be compiled for parallel usage, the MPI and ParMETIS paths must be set.
+
+\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 commands:
+\begin{verbatim}
+> libtoolize --force --copy
+> aclocal
+> autoconf
+> automake --add-missing --copy
+\end{verbatim}
+Then repeat steps 3 and 4.
+For the additional steps \verb+libtool+, \verb+automake+ and \verb+autotools+ must be installed on your system.
+
+Further information about the installation process can be found in the \verb+README+ file in the AMDiS source directory.
+
+\section{Compilation of an example application}
+\label{s:compilation}
+
+For the compilation of the examples, described in this section, the following steps must be executed:
+\begin{enumerate}
+\item Get the sources:
+ \begin{itemize}
+ \item Unpack an archive file: \\
+ \verb+> gunzip demo.tar.gz+\\
+ \verb+> tar xvf demo.tar+
+
+ or
+ \item CVS checkout:\\
+ \verb+> export CVSROOT=<CVSROOT>+\\
+ \verb+> cvs checkout demo+
+
+ or
+ \item SVN checkout:\\
+ \verb+> svn checkout file://<SVN-REPOSITORY-PATH>/demo+
+ \end{itemize}
+
+\item Change into the demo directory:\\
+ \verb+> cd demo+
+
+\item Edit 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 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:
+\begin{itemize}
+\item In the sequential case:\\
+ \verb+> ./<PROG-NAME> <PARAMETER-FILE>+
+\item In the parallel case:\\
+ \verb+> mpirun <MPI-OPTIONS> ./<PROG-NAME> <PARAMETER-FILE>+\\
+
+ The \verb+<MPI-OPTIONS>+ at least should contain the number of used processes, which is given by \verb+-np <NUM-PROCS>+. For further MPI options see \verb+http://www-unix.mcs.anl.gov/mpi/+.
+\end{itemize}
+
+\input{makefile.tex}
+
diff --git a/doc/tutorial/introduction.tex b/doc/tutorial/introduction.tex
new file mode 100644
index 0000000000000000000000000000000000000000..c6e8ae29fddaf22804a8962e233206953e85015b
--- /dev/null
+++ b/doc/tutorial/introduction.tex
@@ -0,0 +1,35 @@
+\chapter{Introduction}
+\label{c:tutorial introduction}
+
+The objective of this tutorial is to introduce the user into the main
+AMDiS features by giving some application examples.
+
+Section \ref{c:amdis installation} describes the installation of the
+AMDiS library and the building of user applications step by step. The
+corresponding application makefile is given in Section
+\ref{s:application makefile}.
+
+In Section \ref{c:examples}, for every example the following aspects
+are described:\begin{itemize}
+\item {\bf Abstract problem description}: In the header of each
+ example section, the abstract problem definition is
+ given. Sometimes, some solution strategies on a high abstraction
+ level are mentioned, also.
+\item {\bf Source code}: In the source code section, the listing of
+ the example source code is explained.
+\item {\bf Parameter file}: In this section, the parameter file is
+ described. The parameter file contains parameters which are read be
+ the application at runtime. The name of the parameter file is
+ usually passed to the the application as a command line argument.
+\item {\bf Macro file}: In the macro file section, the definition of
+ the coarse macro mesh is shown, which is the basis for adaptive
+ refinements.
+\item {\bf Output}: The AMDiS results are written to output files that
+ contain the final mesh and the problem solution on this mesh. The
+ output can be visualized by proper tools ({\it CrystalClear,
+ ParaView, TecPlot}). In the output section, the visualized problem
+ results are shown and discussed.
+\end{itemize}
+To avoid unnecessary repetitions, not every aspect of every example is
+described, but only those aspects that have not appeared in previous
+examples.
diff --git a/doc/tutorial/makefile.tex b/doc/tutorial/makefile.tex
new file mode 100644
index 0000000000000000000000000000000000000000..fcfde4ba432f8343b4c986d2015db6465c685fa4
--- /dev/null
+++ b/doc/tutorial/makefile.tex
@@ -0,0 +1,126 @@
+\chapter{Application makefile}
+\label{s:application makefile}
+
+In this section, the organization of the application makefile is described which is used for the examples in this tutorial. The same organization can be used for other user applications, too.
+
+In the first block, user flags and directories are specified.
+
+\begin{verbatim}
+# ============================================================================
+# ===== flags and directories (to be modified by the user) ===================
+# ============================================================================
+
+USE_PARALLEL_AMDIS = 0 # 0: sequential AMDiS, 1: parallel AMDiS
+DEBUG = 0 # 0: no debug, 1: debug mode
+
+AMDIS_DIR = <AMDIS-DIR> # fill the AMDiS installation path here
+MPI_DIR = <MPI-DIR> # fill the MPI installation path here
+PARMETIS_DIR = <PARMETIS-DIR> # fill the ParMETIS installation path here
+\end{verbatim}
+
+If \verb+USE_PARALLEL_AMDIS+ is set to $1$, parallel applications will be supported. A necessary condition is that the AMDiS library is configured for parallelization, too (see Chapter \ref{c:amdis installation}). The \verb+DEBUG+ entry specifies, whether applications should be compiled in debug mode, or not. This entry is independent of the corresponding AMDiS settings, but if the AMDiS library was not compiled in debug mode, only application code can be debugged.
+
+\verb+AMDIS_DIR+ stores the AMDiS installation path and must be set by the user. This is the path given to the AMDiS configure script by the \verb+--prefix+ option. The values of \verb+MPI_DIR+ and \verb+PARMETIS_DIR+ only are needed if parallelization should be supported. Here, the installation pathes of MPI and ParMETIS are stored.
+
+In the next block, include pathes are defined.
+
+\begin{verbatim}
+# ============================================================================
+# ===== includes pathes ======================================================
+# ============================================================================
+
+AMDIS_INCLUDE = -I$(AMDIS_DIR)/src
+MPI_INCLUDE = -I$(MPI_DIR)/include
+PARMETIS_INCLUDE = -I$(PARMETIS_DIR)
+
+INCLUDES = -I. $(AMDIS_INCLUDE) $(MPI_INCLUDE) $(PARMETIS_INCLUDE)
+\end{verbatim}
+
+Now, we introduce the needed libraries.
+
+\begin{verbatim}
+# ============================================================================
+# ===== libraries ============================================================
+# ============================================================================
+
+AMDIS_LIB = -L$(AMDIS_DIR)/lib -lamdis
+PARMETIS_LIB = -L$(PARMETIS_DIR) -lparmetis -lmetis
+
+LIBS = $(AMDIS_LIB)
+\end{verbatim}
+
+By default, \verb+LIBS+ contains only the AMDiS library. If \verb+USE_PARALLEL_AMDIS+ is $1$, \verb+LIBS+ is extended by the ParMETIS library. In the same way, other libraries can be added. In the sequential case, we use the GNU C++ compiler \verb-g++-, in the parallel case, the MPI C++ compiler \verb+mpiCC+.
+
+\begin{verbatim}
+# ============================================================================
+# ===== parallel or sequential ? =============================================
+# ============================================================================
+
+ifeq ($(USE_PARALLEL_AMDIS), 0)
+ COMPILE = g++
+else
+ COMPILE = $(MPI_DIR)/bin/mpiCC
+ LIBS += $(PARMETIS_LIB)
+endif
+\end{verbatim}
+
+The next block sets the compile flags. In debug mode, we use no optimization (\verb+-O0+) and add symbolic debug information (\verb+-g+). Otherwise, we compile with optimazation level 2 (\verb+-O2+).
+
+\begin{verbatim}
+# ============================================================================
+# ===== compile flags ========================================================
+# ============================================================================
+
+ifeq ($(DEBUG), 0)
+ CPPFLAGS = -O2
+else
+ CPPFLAGS = -g -O0
+endif
+\end{verbatim}
+
+We use the \verb+libtool+ in the AMDiS installation path for linking.
+
+\begin{verbatim}
+# ============================================================================
+# ===== libtool linking ======================================================
+# ============================================================================
+
+LIBTOOL = $(AMDIS_DIR)/libtool
+LINK = $(LIBTOOL) --mode=link $(COMPILE)
+\end{verbatim}
+
+Now, we define rules to create and delete objects files.
+
+\begin{verbatim}
+# ============================================================================
+# ===== rules ================================================================
+# ============================================================================
+
+clean:
+ -rm -rf *.o
+
+.cc.o: $*.cc
+ $(COMPILE) $(INCLUDES) $(CPPFLAGS) -c -o $*.o $^
+
+\end{verbatim}
+
+The second rule creates needed object files automatically using the corresponding C++ files.
+
+Finally, we define rules for the linking of user applications. Here, we present only the rule for the \verb+ellipt+ application. Other applications can be created in an analog way.
+
+\begin{verbatim}
+# ============================================================================
+# ===== user programs ========================================================
+# ============================================================================
+
+VPATH = .:./src:$(AMDIS_DIR)/src
+
+# ===== myprog ===============================================================
+
+ELLIPT_OFILES = ellipt.o
+
+ellipt: $(ELLIPT_OFILES)
+ $(LINK) $(CPPFLAGS) -o ellipt $(ELLIPT_OFILES) $(LIBS)
+\end{verbatim}
+
+The \verb+VPATH+ variable contains all pathes, where sources can be located. The \verb+ellipt+ rule first creates all needed object files defined in \verb+ELLIPT_OFILES+. In this example, only \verb+ellipt.o+ is needed. Then all needed object files and libraries are linked together. The \verb+-o+ option specifies that the executable will be written to the file \verb+ellipt+.
diff --git a/doc/tutorial/multigrid.tex b/doc/tutorial/multigrid.tex
new file mode 100644
index 0000000000000000000000000000000000000000..012926e2aad649f9bb3f1b489c282ece5161b9d4
--- /dev/null
+++ b/doc/tutorial/multigrid.tex
@@ -0,0 +1,87 @@
+\section{Multigrid}
+The multigrid method is an effective way to solve large linear systems of equations. Its complexity is proportional to the number of unknowns in the system of equations. This can be achieved by using the hierarchical mesh structure of AMDiS. The multigrid idea is based on the following two principles:
+\begin{itemize}
+\item {\it Smoothing principle}: Classical iterative methods often have a strong smoothing effect on the error.
+\item {\it Coarse grid principle}: A quantity that is smooth on a given grid can be approximated on a coarser grid without any essential loss of information.
+\end{itemize}
+To solve $L_h u_h = f_h$ on the fine grid, we start with an arbitrary initial guess of the solution $u_h$. Then we apply the following steps, until a given tolerance criterion for the solution is fulfilled:
+\begin{enumerate}
+\item Apply some smoothing steps to $u_h$ on the fine grid (pre-smoothing).
+\item Compute the fine grid residual $d_h = f_h - L_h u_h$.
+\item Restrict the residual $d_h$ to the coarse grid ($d_h \rightarrow d_H$).
+\item Solve the defect equation $L_H v_H = d_H$ on the coarse grid.
+\item Interpolate $v_H$ to the fine grid ($v_H \rightarrow v_h$).
+\item Correct the solution: $u_h = u_h + v_h$.
+\item Apply some smoothing steps to $u_h$ on the fine grid (post-smoothing).
+\end{enumerate}
+To solve the defect equation on the coarse level, we can apply the multigrid method recursively. This can be done once ({\it V-cycle}) or twice ({\it W-cycle}). At the coarsest level, the system of equations can be solved by a direct solver or again some smoothing steps are applied to it (in AMDiS, so far, no direct solver is applied at the coarsest level).
+
+To use the multigrid solver in AMDiS, only the parameter file has to be modified. So, we omit the other sections here.
+
+\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 multigrid.cc} \\
+\hline
+{\bf parameter file} & \tt init/ & \tt multigrid.dat.1d & \tt multigrid.dat.2d & \tt multigrid.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 multigrid.mesh, multigrid.dat} \\
+\hline
+\end{tabular}
+\caption{Files of the {\tt multigrid} example.}
+\label{t:multigrid files}
+\end{table}
+
+\subsection{Parameter file}
+
+First, we have to avoid, that DOFs at coarse levels are freed, if they are not used at finer levels.
+
+\begin{lstlisting}{}
+multigridMesh->preserve coarse dofs: 1
+\end{lstlisting}
+
+Now, we choose \verb+mg+ as solver, which is the key for the multigrid solver in AMDiS.
+
+\begin{lstlisting}{}
+multigrid->solver: mg
+\end{lstlisting}
+
+The next three lines are not multigrid specific. They determine the solver tolerance, the maximal number of solver iterations, and the pre-conditioner.
+
+\begin{lstlisting}{}
+multigrid->solver->tolerance: 1.e-12
+multigrid->solver->max iteration: 100
+multigrid->solver->left precon: diag
+\end{lstlisting}
+
+Now, multigrid specific parameters follow.
+
+\begin{lstlisting}{}
+multigrid->solver->use galerkin operator: 0
+multigrid->solver->mg cycle index: 1 % 1: V, 2: W
+multigrid->solver->smoother: gs
+multigrid->solver->smoother->omega: 1.0
+multigrid->solver->pre smoothing steps: 3
+multigrid->solver->post smoothing steps: 3
+multigrid->solver->coarse level smoothing steps: 1
+multigrid->solver->min level: 0
+multigrid->solver->min level gap: 1
+multigrid->solver->max mg levels: 100
+\end{lstlisting}
+
+The entry \verb+use galerkin operator+ determines, how the system on coarse levels is assembled. If the value is set to 1, the system of the finest level is successively restricted to coarser levels by the galerkin operator (which is only defined for linear Lagrange elements). If the value is set to 0, the system of coarser levels is assembled using the usual problem operators on the coarser meshes.
+
+The \verb+mg cycle index+ determines the number of recursive multigrid calls within each level. A value of $1$ results in V-cycle iterations, a value of $2$ in W-cycle iterations (in principle, also higher values could be chosen).
+
+Currently, as smoother only a relaxed Gauss-Seidel method is implemented (\verb+gs+). The value \verb+smoother->omega+ is the relaxation parameter.
+
+The meaning of \verb+pre smoothing steps+, \verb+post smoothing steps+ and \verb+coarse level smoothing+ \verb+steps+ is explained above.
+
+\verb+min level+ is the coarsest multigrid level. By default it is 0.
+
+\verb+min level gap+ describes the number of mesh levels that at least are skipped between two multigrid levels. \verb+max mg levels+ sets the maximum number of allowed multigrid levels. If necessary, more than \verb+min level gap+ levels are skipped, to fulfill this requirement.
diff --git a/doc/tutorial/neumann.tex b/doc/tutorial/neumann.tex
new file mode 100644
index 0000000000000000000000000000000000000000..e7153b9d1d5d7eeca3a1bc20e86e95097f38623e
--- /dev/null
+++ b/doc/tutorial/neumann.tex
@@ -0,0 +1,96 @@
+\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).
+
+\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 neumann.cc} \\
+\hline
+{\bf parameter file} & \tt init/ & \tt neumann.dat.1d & \tt neumann.dat.2d & \tt neumann.dat.3d \\
+\hline
+{\bf macro file} & \tt macro/ & \tt neumann.macro.1d & \tt neumann.macro.2d & \tt neumann.macro.3d \\
+\hline
+{\bf output files} & \tt output/ & \multicolumn{3}{|c|}{\tt neumann.mesh, neumann.dat} \\
+\hline
+\end{tabular}
+\caption{Files of the {\tt neumann} example.}
+\end{table}
+
+\subsection{Source code}
+Only a few changes in the source code are necessary to apply Neumann boundary conditions. First, we define the function $N=1$.\\
+
+\begin{lstlisting}{}
+class N : public AbstractFunction<double, WorldVector<double> >
+{
+public:
+ MEMORY_MANAGED(N);
+
+ double operator()(const WorldVector<double>& x) const
+ {
+ return 1.0;
+ }
+};
+\end{lstlisting}
+
+In the main program we add the boundary conditions to our problem \verb+neumann+.
+
+\begin{lstlisting}{}
+int main(int argc, char* argv[])
+{
+ ...
+ neumann.addNeumannBC(1, NEW N);
+ neumann.addDirichletBC(2, 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.
+
+\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.
+
+\subsection{Macro file}
+The file \verb+neumann.macro.2d+ is listed below:
+
+\begin{lstlisting}{}
+DIM: 2
+DIM_OF_WORLD: 2
+
+number of vertices: 5
+number of elements: 4
+
+vertex coordinates:
+-0.5 -0.5
+ 0.5 -0.5
+ 0.5 0.5
+-0.5 0.5
+ 0.0 0.0
+
+element vertices:
+0 1 4
+1 2 4
+2 3 4
+3 0 4
+
+element boundaries:
+0 0 2
+0 0 1
+0 0 2
+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.
+
+\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)$.
+\begin{figure}
+\center
+\includegraphics[width=0.5\textwidth]{fig/neumann_solution.jpg}
+\caption{Solution of the problem with Neumann boundary conditions at two sides.}
+\label{f:neumann solution}
+\end{figure}
diff --git a/doc/tutorial/nonlin.tex b/doc/tutorial/nonlin.tex
new file mode 100644
index 0000000000000000000000000000000000000000..3ed697b5903d3aa661f7651f693059bf0f2dcfb3
--- /dev/null
+++ b/doc/tutorial/nonlin.tex
@@ -0,0 +1,479 @@
+\section{Nonlinear problem}
+\label{ss:nonlinear problem}
+We define the nonlinear problem
+\begin{eqnarray}
+-\Delta u + u^4 &=& f~~~\mbox{in } \Omega \subset \mathbb{R}^{dim}\\
+u &=& g~~~ \mbox{on } \partial\Omega.
+\end{eqnarray}
+We choose the functions $f$ and $g$ so that the exact solution again is $u(x)=e^{-10x^2}$. This leads to
+\begin{eqnarray}
+f(x)&=& -\left(400-20dow\right)e^{-10x^2} + \left(e^{-10x^2}\right)^4\\
+g(x)&=&e^{-10x^2},
+\end{eqnarray}
+with $dow$ the world dimension.
+
+We linearize the problem using the Newton method. First, we define an initial guess $u_0$ of the solution which is $0$ for the first adaptation loop iteration. In later iterations we can use the solution of the last iteration interpolated to the current mesh as initial guess. In each Newton step, a correction $d$ for the solution of the last step is computed by solving
+\begin{eqnarray}
+\label{eq:newton step}
+DF(u_n)(d) = F(u_n)
+\end{eqnarray}
+for $d$, where $F(u):=-\Delta u + u^4 - f$ and
+\begin{eqnarray}
+DF(u_n)(d) &=& \lim_{h \rightarrow 0} \frac{F(u_n+hd)-F(u_n)}{h}\\
+&=& \lim_{h \rightarrow 0} \frac{-\Delta u_n - h \Delta d + \Delta u_n}{h} + \lim_{h \rightarrow 0}\frac{(u_n+hd)^4-u_n^4}{h}\\
+&=& -\Delta d + 4 u_n^3 d
+\end{eqnarray}
+the directional derivative of $F$ at $u_n$ along $d$.
+
+Then the solution is updated:
+\begin{equation}
+u_{n+1} := u_n - d.
+\end{equation}
+We repeat this procedure until $||d||_{L^2} < tol$ with $tol$ a given tolerance for the Newton method.
+
+In our example, equation (\ref{eq:newton step}) reads:
+\begin{equation}
+-\Delta d + 4 u_n^3 d = - \Delta u_n + u_n^4 - f.
+\end{equation}
+
+\begin{figure}
+\center
+\includegraphics[width=0.7\textwidth]{fig/nonlin_iteration.pdf}
+\caption{Solve step of the nonlinear problem.}
+\label{f:nonlin iteration}
+\end{figure}
+In Figure \ref{f:nonlin iteration}, the Newton method is illustrated.
+
+\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 nonlin.cc} \\
+\hline
+{\bf parameter file} & \tt init/ & \tt nonlin.dat.1d & \tt nonlin.dat.2d & \tt nonlin.dat.3d \\
+\hline
+{\bf macro file} & \tt macro/ & \tt macro\_big.stand.1d & \tt macro\_big.stand.2d & \tt macro\_big.stand.3d \\
+\hline
+{\bf output files} & \tt output/ & \multicolumn{3}{|c|}{\tt nonlin.mesh, nonlin.dat} \\
+\hline
+\end{tabular}
+\caption{Files of the {\tt nonlin} example.}
+\end{table}
+
+\subsection{Source code}
+The main idea is to realize the Newton method as implementation of the {\it ProblemIterationInterface}, which replaces the standard iteration in the adaptation loop. The Newton method has to know the problem which has to be solved, as well as the implementation of one Newton step. Estimation and adaptation is done by the problem object. The assemble step and solve step are delegated to a Newton-step object.
+
+Now, we describe the code step by step. The function $g$ is defined like in the previous examples. In the following, we define a zero-function which is later used to implement the Dirichlet boundary condition for the Newton-step implementation (at domain boundaries no correction has to be done). The function \verb+F+ implements the right hand side function $f$.
+
+\begin{lstlisting}{}
+class Zero : public AbstractFunction<double, WorldVector<double> >
+{
+public:
+ MEMORY_MANAGED(Zero);
+
+ double operator()(const WorldVector<double>& x) const {
+ return 0.0;
+ }
+};
+
+class F : public AbstractFunction<double, WorldVector<double> >
+{
+public:
+ MEMORY_MANAGED(F);
+
+ /** \brief
+ * Constructor
+ */
+ F(int degree)
+ : AbstractFunction<double, WorldVector<double> >(degree)
+ {};
+
+ /** \brief
+ * Implementation of AbstractFunction::operator().
+ */
+ double operator()(const WorldVector<double>& x) const {
+ int dow = x.getSize();
+ double r2 = x * x;
+ double ux = exp(-10.0 * r2);
+ double ux4 = ux * ux * ux * ux;
+ return ux4 -(400.0 * r2 - 20.0 * dow) * ux;
+ }
+};
+\end{lstlisting}
+
+The class \verb+X3+ implements the function $u^3(x)$ used within the Newton step.
+
+\begin{lstlisting}{}
+class X3 : public AbstractFunction<double, double>
+{
+public:
+ MEMORY_MANAGED(X3);
+
+ X3() : AbstractFunction<double, double>(3) {};
+
+ /** \brief
+ * Implementation of AbstractFunction::operator().
+ */
+ double operator()(const double& x) const {
+ return x * x * x;
+ }
+};
+\end{lstlisting}
+
+In the following, we define an interface which has to be implemented by the Newton-step object.
+
+\begin{lstlisting}{}
+class NewtonStepInterface
+{
+public:
+ virtual void initNewtonStep(AdaptInfo *adaptInfo) = 0;
+ virtual void exitNewtonStep(AdaptInfo *adaptInfo) = 0;
+ virtual void assembleNewtonStep(AdaptInfo *adaptInfo, Flag flag) = 0;
+ virtual void solveNewtonStep(AdaptInfo *adaptInfo) = 0;
+ virtual DOFVector<double> *getCorrection() = 0;
+};
+\end{lstlisting}
+
+The \verb+initNewtonStep+ method is called before each Newton step, the method \verb+exitNewtonStep+ after each Newton step. \verb+assembleNewtonStep+ assembles the linear system of equations needed for the next step, \verb+solveNewtonStep+ solves this system of equations. The solution is the correction $d$. The method \verb+getCorrection+ returns a pointer to the vector storing the correction.
+
+Now, the Newton method will be implemented. Actually, the class \verb+NewtonMethod+ replaces the whole iteration in the adaptation loop, including mesh adaptation and error estimation. The Newton method, which is a loop over Newton steps, is one part of this iteration.
+
+\begin{lstlisting}{}
+class NewtonMethod : public ProblemIterationInterface
+{
+public:
+ NewtonMethod(const char *name,
+ ProblemScal *problem,
+ NewtonStepInterface *step)
+ : problemNonlin(problem),
+ newtonStep(step),
+ newtonTolerance(1e-8),
+ newtonMaxIter(100)
+ {
+ GET_PARAMETER(0, std::string(name) + "->tolerance", "%f",
+ &newtonTolerance);
+ GET_PARAMETER(0, std::string(name) + "->max iteration", "%d",
+ &newtonMaxIter);
+ solution = problemNonlin->getSolution();
+ correction = newtonStep->getCorrection();
+ };
+\end{lstlisting}
+
+In the constructor, pointers to the nonlinear problem and to the Newton-step object are stored to the class members \verb+problemNonlin+ and \verb+newtonStep+. Furthermore, the parameters \verb+newtonTolerance+ and \verb+newtonMaxIter+ are initialized, and pointers to the nonlinear solution and to the correction vector are stored.
+
+The following methods define one iteration in the adaptation loop.
+
+\begin{lstlisting}{}
+ void beginIteration(AdaptInfo *adaptInfo)
+ {
+ FUNCNAME("NewtonMethod::beginIteration()");
+ MSG("\n");
+ MSG("begin of iteration %d\n", adaptInfo->getSpaceIteration()+1);
+ MSG("=============================\n");
+ }
+
+ Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION)
+ {
+ Flag flag = 0, markFlag = 0;
+
+ if(toDo.isSet(MARK)) markFlag = problemNonlin->markElements(adaptInfo);
+ if(toDo.isSet(ADAPT) && markFlag.isSet(MESH_REFINED))
+ flag = problemNonlin->refineMesh(adaptInfo);
+ if(toDo.isSet(ADAPT) && markFlag.isSet(MESH_COARSENED))
+ flag |= problemNonlin->coarsenMesh(adaptInfo);
+
+ if(toDo.isSet(SOLVE)) {
+ newtonStep->initNewtonStep(adaptInfo);
+ int newtonIteration = 0;
+ double res = 0.0;
+ do {
+ newtonIteration++;
+ newtonStep->assembleNewtonStep(adaptInfo, flag);
+ newtonStep->solveNewtonStep(adaptInfo);
+ res = correction->L2Norm();
+ *solution -= *correction;
+ MSG("newton iteration %d: residual %f (tol: %f)\n",
+ newtonIteration, res, newtonTolerance);
+ } while((res > newtonTolerance) && (newtonIteration < newtonMaxIter));
+
+ newtonStep->exitNewtonStep(adaptInfo);
+ }
+
+ if(toDo.isSet(ESTIMATE)) problemNonlin->estimate(adaptInfo);
+ return flag;
+ };
+
+ void endIteration(AdaptInfo *adaptInfo)
+ {
+ FUNCNAME("NewtonMethod::endIteration()");
+ MSG("\n");
+ MSG("end of iteration number: %d\n", adaptInfo->getSpaceIteration()+1);
+ MSG("=============================\n");
+ }
+\end{lstlisting}
+
+The methods \verb+beginIteration+ and \verb+endIteration+ only print some information to the standard output. In \verb+oneIteration+, the iteration, including the loop over the Newton steps, is defined.
+
+Finally, the methods \verb+getNumProblems+ and \verb+getProblem+ are implemented to complete the \\\verb+ProblemIterationInterface+, and the private class members are defined.
+
+\begin{lstlisting}{}
+ int getNumProblems() { return 1; };
+
+ ProblemStatBase *getProblem(int number = 0)
+ {
+ FUNCNAME("NewtonMethod::getProblem()");
+ if(number == 0) return problemNonlin;
+ ERROR_EXIT("invalid problem number\n");
+ return NULL;
+ };
+
+private:
+ ProblemScal *problemNonlin;
+ NewtonStepInterface *newtonStep;
+ double newtonTolerance;
+ int newtonMaxIter;
+ DOFVector<double> *solution;
+ DOFVector<double> *correction;
+};
+\end{lstlisting}
+
+The class \verb+Nonlin+ implements both, the nonlinear problem and the Newton-step. Since the Newton step is accessed over an own interface, it is always clear, in which role a \verb+Nonlin+ instance is called by the Newton method.
+
+\begin{lstlisting}{}
+class Nonlin : public ProblemScal,
+ public NewtonStepInterface
+{
+public:
+ Nonlin(const char *name)
+ : ProblemScal(name)
+ {};
+\end{lstlisting}
+
+In the constructor, the base class constructor of \verb+ProblemScal+ is called and the name is given to it.
+
+In the initialization, the base class intialization is called, the correction vector is created and initialized, and Dirichlet boundary conditions are created.
+
+\begin{lstlisting}{}
+ void initialize(Flag initFlag,
+ ProblemScal *adoptProblem = NULL,
+ Flag adoptFlag = INIT_NOTHING)
+ {
+ ProblemScal::initialize(initFlag, adoptProblem, adoptFlag);
+ correction = NEW DOFVector<double>(this->getFESpace(), "old solution");
+ correction->set(0.0);
+
+ dirichletZero = NEW DirichletBC(1, &zero, feSpace_);
+ dirichletG = NEW DirichletBC(1, &g, feSpace_);
+
+ solution_->getBoundaryManager()->addBoundaryCondition(dirichletG);
+ systemMatrix_->getBoundaryManager()->
+ addBoundaryCondition(dirichletZero);
+ rhs_->getBoundaryManager()->addBoundaryCondition(dirichletZero);
+ correction->getBoundaryManager()->addBoundaryCondition(dirichletZero);
+ };
+\end{lstlisting}
+
+To the solution of the nonlinear problem the function $g$ is applied as Dirichlet boundary function. The system matrix, the correction vector and the right hand side vector build the system for the Newton step. Since no correction has to be done at the domain boundaries, zero Dirichlet conditions are applied to them.
+
+In the destructor, the allocated memory is freed.
+
+\begin{lstlisting}{}
+ ~Nonlin()
+ {
+ DELETE correction;
+ DELETE dirichletZero;
+ DELETE dirichletG;
+ };
+\end{lstlisting}
+
+Now, we implement the Newton step functionality. First, in \verb+initNewtonStep+, we fill the solution vector with boundary values. This will not be done automatically because we let the \verb+solution_+ pointer point to \verb+correction+. The address of \verb+solution_+ is stored in \verb+tmp+. After the Newton method is finished, the \verb+solution_+ pointer is reset to its original value in \verb+exitNewtonStep+.
+
+\begin{lstlisting}{}
+ void initNewtonStep(AdaptInfo *adaptInfo) {
+ solution_->getBoundaryManager()->initVector(solution_);
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(mesh_, -1,
+ Mesh::CALL_LEAF_EL |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_BOUND);
+ while(elInfo) {
+ solution_->getBoundaryManager()->fillBoundaryConditions(elInfo,
+ solution_);
+ elInfo = stack.traverseNext(elInfo);
+ }
+ solution_->getBoundaryManager()->exitVector(solution_);
+
+ tmp = solution_;
+ solution_ = correction;
+ };
+
+ void exitNewtonStep(AdaptInfo *adaptInfo) {
+ solution_ = tmp;
+ };
+\end{lstlisting}
+
+The implementation of \verb+assembleNewtonStep+ and \verb+solveNewtonStep+ just delegates the calls to the base class implementations in \verb+ProblemScal+.
+
+\begin{lstlisting}{}
+ void assembleNewtonStep(AdaptInfo *adaptInfo, Flag flag) {
+ ProblemScal::buildAfterCoarsen(adaptInfo, flag);
+ };
+
+ void solveNewtonStep(AdaptInfo *adaptInfo) {
+ ProblemScal::solve(adaptInfo);
+ };
+\end{lstlisting}
+
+Finally, the \verb+getCorrection+ method is implemented and private class members are defined.
+
+\begin{lstlisting}{}
+ DOFVector<double> *getCorrection() { return correction; };
+
+private:
+ DOFVector<double> *correction;
+ DOFVector<double> *tmp;
+ Zero zero;
+ G g;
+ DirichletBC *dirichletZero;
+ DirichletBC *dirichletG;
+};
+\end{lstlisting}
+
+Now, we start with the main program.
+
+\begin{lstlisting}{}
+int main(int argc, char* argv[])
+{
+ FUNCNAME("main");
+
+ TEST_EXIT(argc == 2)("usage: nonlin initfile\n");
+
+ Parameters::init(false, argv[1]);
+
+ Nonlin nonlin("nonlin");
+ nonlin.initialize(INIT_ALL);
+
+ AdaptInfo *adaptInfo = NEW AdaptInfo("nonlin->adapt", 1);
+
+ NewtonMethod newtonMethod("nonlin->newton", &nonlin, &nonlin);
+
+ AdaptStationary *adapt = NEW AdaptStationary("nonlin->adapt",
+ &newtonMethod,
+ adaptInfo);
+\end{lstlisting}
+
+An instance of class \verb+NewtonMethod+ was created with a pointer to a \verb+Nonlin+ object as nonlinear problem and as Newton-step implementation. Instead of the nonlinear problem, now, the object \verb+newtonMethod+ is handed to the adaptation loop as implementation of \verb+ProblemIterationInterface+.
+
+We have to add operators representing the Newton step equation $-\Delta d + 4u_n^3d=-\Delta u_n + u_n^4 -f$ as well as operators representing the nonlinear problem $-\Delta u_n + u_n^4 = f$. When the operators are given to the problem, one can determine an assemble factor (second argument of \verb+addMatrixOperator+ and \verb+addVectorOperator+) as well as an estimation factor (third argument) for each operator. So, we can manage both equations in one problem instance. Note that in the Newton step we solve for $d$. $u_n$ there is known from the last Newton step. The term $u_n^4$ was implemented as $u_n^3 v$, where $v$ for the Newton step equation is equal to $d$ and in the nonlinear problem equation equal to $u_n$. So, the corresponding operator can be used in both equations just with different factors. In Figure \ref{f:nonlin operators}, the operator factors for the assemble and the estimate step are shown.
+\begin{figure}
+\center
+\includegraphics[width=0.5\textwidth]{fig/nonlin_operators.pdf}
+\caption{Operator factors for the assemble step and for the estimate step.}
+\label{f:nonlin operators}
+\end{figure}
+
+\begin{lstlisting}{}
+ // ===== create operators =====
+ double four = 4.0;
+ double one = 1.0;
+ double zero = 0.0;
+ double minusOne = -1.0;
+
+ Operator *nonlinOperator0 = NEW Operator(Operator::MATRIX_OPERATOR |
+ Operator::VECTOR_OPERATOR,
+ nonlin.getFESpace());
+
+ nonlinOperator0->setUhOld(nonlin.getSolution());
+ nonlinOperator0->addZeroOrderTerm(NEW VecAtQP_ZOT(nonlin.getSolution(),
+ NEW X3));
+
+ nonlin.addMatrixOperator(nonlinOperator0, &four, &one);
+ nonlin.addVectorOperator(nonlinOperator0, &one, &zero);
+
+ Operator *nonlinOperator2 = NEW Operator(Operator::MATRIX_OPERATOR |
+ Operator::VECTOR_OPERATOR,
+ nonlin.getFESpace());
+
+ nonlinOperator2->setUhOld(nonlin.getSolution());
+ nonlinOperator2->addSecondOrderTerm(NEW Laplace_SOT);
+
+ nonlin.addMatrixOperator(nonlinOperator2, &one, &one);
+ nonlin.addVectorOperator(nonlinOperator2, &one, &zero);
+
+ int degree = nonlin.getFESpace()->getBasisFcts()->getDegree();
+
+ Operator* rhsFunctionOperator = NEW Operator(Operator::VECTOR_OPERATOR,
+ nonlin.getFESpace());
+ rhsFunctionOperator->addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree)));
+
+ nonlin.addVectorOperator(rhsFunctionOperator, &minusOne, &one);
+\end{lstlisting}
+
+Finally, the adaptation loop is started and after it is fnished the result is written.
+
+\begin{lstlisting}{}
+ adapt->adapt();
+ nonlin.writeFiles(adaptInfo, true);
+}
+\end{lstlisting}
+
+\subsection{Parameter file}
+The used parameter file \verb+nonlin.dat.2d+ looks like:
+
+\begin{lstlisting}{}
+dimension of world: 2
+
+nonlinMesh->macro file name: ./macro/macro_big.stand.2d
+nonlinMesh->global refinements: 0
+
+nonlin->mesh: nonlinMesh
+
+nonlin->dim: 2
+nonlin->polynomial degree: 1
+
+nonlin->newton->tolerance: 1e-8
+nonlin->newton->max iteration: 100
+
+nonlin->solver: cg
+nonlin->solver->max iteration: 1000
+nonlin->solver->tolerance: 1.e-8
+nonlin->solver->left precon: diag
+
+nonlin->estimator: residual
+nonlin->estimator->C0: 0.1
+nonlin->estimator->C1: 0.1
+
+nonlin->marker->strategy: 2
+nonlin->marker->MSGamma: 0.5
+
+nonlin->adapt->tolerance: 1e-1
+nonlin->adapt->max iteration: 100
+
+nonlin->output->filename: output/nonlin
+nonlin->output->AMDiS format: 1
+nonlin->output->AMDiS mesh ext: .mesh
+nonlin->output->AMDiS data ext: .dat
+\end{lstlisting}
+
+Here, as macro file \verb+macro_big.stand.2d+ is used, which is described in Section \ref{s:nonlin macro}. The parameters \verb+nonlin->newton->tolerance+ and \verb+nonlin->newton->max iteration+ determine the tolerance of the Newton solver and the maximal number of Newton steps within each iteration of the adaptation loop.
+
+\subsection{Macro file}
+\label{s:nonlin macro}
+The used macro file \verb+macro_big.stand.2d+ is very similar to the macro file of the first example described in Section \ref{s:ellipt macro}. Only the domain was changed from $\Omega=[0,1]^2$ to $\Omega=[-1, 1]^2$ by adapting the vertex coordinates correspondingly.
+
+\subsection{Output}
+
+\begin{figure}
+\center
+\begin{tabular}{cc}
+\includegraphics[width=0.4\textwidth]{fig/nonlin_solution.jpg}&
+\includegraphics[width=0.4\textwidth]{fig/nonlin_mesh.jpg}\\
+solution & mesh
+\end{tabular}
+\caption{Solution and final mesh of the nonlinear problem}
+\label{f:nonlin mesh and solution}
+\end{figure}
+
+In Figure \ref{f:nonlin mesh and solution}, the solution and the final mesh written after the adaptation loop are visualized. The solution is shown as height field where the values are interpreted as z-coordinates.
\ No newline at end of file
diff --git a/doc/tutorial/operatorterms.tex b/doc/tutorial/operatorterms.tex
new file mode 100644
index 0000000000000000000000000000000000000000..ddf91c7ded0098226b7cb2901774551c4e859d8c
--- /dev/null
+++ b/doc/tutorial/operatorterms.tex
@@ -0,0 +1,11 @@
+\section{Operator term reference}
+\label{s:operator term reference}
+
+\subsection{Second order terms}
+\label{ss:second order term reference}
+
+\subsection{First order terms}
+\label{ss:first order term reference}
+
+\subsection{Zero order terms}
+\label{ss:zero order term reference}
diff --git a/doc/tutorial/parallelheat.tex b/doc/tutorial/parallelheat.tex
new file mode 100644
index 0000000000000000000000000000000000000000..26faf5e0a6739e7b0fd148c4ae51c5a51459c66f
--- /dev/null
+++ b/doc/tutorial/parallelheat.tex
@@ -0,0 +1,249 @@
+\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:
+ MEMORY_MANAGED(G);
+
+ /** \brief
+ * 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:
+ MEMORY_MANAGED(F);
+
+ F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {};
+
+ /** \brief
+ * 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",
+ ¶llelheat,
+ adaptInfo,
+ ¶llelheat,
+ 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
new file mode 100644
index 0000000000000000000000000000000000000000..15065709b03da7fce2ed8888173facdee161f68f
--- /dev/null
+++ b/doc/tutorial/parameters.tex
@@ -0,0 +1,256 @@
+\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
+\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
+
+\begin{lstlisting}{}
+GET_PARAMETER(<info-level>, "<parameter-name>", "<format-string>", <var-ptr>)+
+\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.
+\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
+
+\begin{lstlisting}{}
+GET_PARAMETER(0, prefix + "->tolerance", "%f", &tolerance);
+\end{lstlisting}
+
+\noindent within the solver constructor reads the init-file line
+
+\begin{lstlisting}{}
+ellipt->solver->tolerance: 1e-4
+\end{lstlisting}
+
+\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.
+
+\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).
+
+\begin{itemize}
+\item \verb+"start time"+ (type: \verb+double+, default: $0.0$): Start time of the simulation.
+\item \verb+"end time"+ (type: \verb+double+, default: $1.0$): End time of the simulation.
+\item \verb+"timestep"+ (type: \verb+double+, default: $0.0$): Inital timestep.
+\item \verb+"min timestep"+ (type: \verb+double+, default: $0.0$): Minimal allowed timestep (for timestep adaptive strategies).
+\item \verb+"max timestep"+ (type: \verb+double+, default: $1.0$): Maximal allowed timestep (for timestep adaptive strategies).
+\item \verb+"max iteration"+ (type: \verb+int+, default: $-1$): Maximal number of space iterations within one timestep (the default value $-1$ stands for infinity).
+\item \verb+"max timestep iteration"+ (type: \verb+int+, default: $30$): Maximal number of timestep adaptations within one timestep.
+%\item \verb+"max time iteration"+ (type: \verb+int+, default: $30$):
+\end{itemize}
+
+An {\it AdaptInfo} instance additionally has parameters that are connected to the problem component. If there is more than one component, \verb+"[<component-number>]"+ is added to the corresponding parameter names (\verb+ellipt->adapt->tolerance[2]+). Othewise the component number can be omitted (\verb+ellipt->adapt->tolerance+). Component dependent parameters are:
+
+\begin{itemize}
+\item \verb+"tolerance"+ (type: \verb+double+, default: $1.0$): Error tolerance of this component in the adaptation loop
+\item \verb+"rel space error"+ (type: \verb+double+, default: $1.0$): The total space tolerance is the product of the total tolerance and the relative space error.
+\item \verb+"rel time error"+ (type: \verb+double+, default: $0.5$): The total time tolerance is the product of the total tolerance, the relative time tolerance and $\theta_1$ (see below).
+\item \verb+"time theta 1"+ (type: \verb+double+, default: $1.0$): see \verb+"rel time error"+.
+\item \verb+"time theta 2"+ (type: \verb+double+, default: $0.3$): The lower bound of the time error is the product of the total tolerance, the relative time tolerance and $\theta_2$.
+\item \verb+"coarsen allowed"+ (type: \verb+int+, default: $0$): Determines whether the component marks for coarsening (0: forbidden, 1: allowed).
+\item \verb+"refinement allowed"+ (type: \verb+int+, default: $1$): Determines whether the component marks for refinement (0: forbidden, 1: allowed).
+\item \verb+"refine bisections"+ (type: \verb+int+, default: $1$): Number of succesively performed bisections on one element when its local error estimate is too large (for this component).
+\item \verb+"coarsen bisections"+ (type: \verb+int+, default: $1$): Number of desired coarsenings of an element with low local error estimate (for this component).
+\end{itemize}
+
+\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.
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\item \verb+""+ (type: \verb++, default: $1$):
+\item \verb+""+ (type: \verb++, default: $1$):
+\item \verb+""+ (type: \verb++, default: $1$):
+\item \verb+""+ (type: \verb++, default: $1$):
+\item \verb+""+ (type: \verb++, default: $1$):
+\item \verb+""+ (type: \verb++, default: $1$):
+\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}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+
+\subsection{Estimator}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{RecoveryEstimator}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{FileWriter}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{ElementFileWriter}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+
+\subsection{MacroReader}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{Marker}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{Mesh}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{Solver}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{GMResSolver}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{MultiGridSolver}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{GSSmoother}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{JacobiSmoother}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{ILUTPreconditioner}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{NonLinSolver}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{NewtonS}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{Problem classes}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{ParallelProblem}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{ProblemInterpolScal}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{ProblemInterpolVec}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{ProblemNonLin}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{ProblemScal}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+\subsection{ProblemVec}
+
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+
+\subsection{Serializer}
+\begin{itemize}
+\item \verb+""+ (type: \verb++, default: $1$):
+\end{itemize}
+
+
diff --git a/doc/tutorial/parametric.tex b/doc/tutorial/parametric.tex
new file mode 100644
index 0000000000000000000000000000000000000000..3487c72663f6686c500752ac1be0934bbb0110f4
--- /dev/null
+++ b/doc/tutorial/parametric.tex
@@ -0,0 +1,321 @@
+\section{Parametric elements}
+\label{s:parametric}
+
+With parametric elements, problems can be solved on meshes which dimensions are not necessarily equal to the world dimension. Therefore, problems on arbitrary manifolds can be solved. Furthermore, the vertex coordinates of the mesh can be flexible. Hence, moving meshes can be implemented.
+
+In this section, we solve equation (\ref{eq:ellipt}) with $f=2x_0$ ($x_0$ is the first component of $x$) on a torus. Then we rotate the torus about the $y$-axis and solve the problem again.
+
+The torus can be created by revolving a circle about an axis coplanar with the circle, which does not touch the circle. We call $r_1$ the radius of the revolved circle and $r_2$ the radius of the revolution, which is the distance of the center of the tube to the center of the torus. In Figure \ref{f:torus}, a torus with its two radii $r_1$ and $r_2$ is shown.
+\begin{figure}
+\center
+\includegraphics[width=0.7\textwidth]{fig/torus.pdf}
+\caption{A torus and a halfed torus with the two radiuses $r_1$ and $r_2$.}
+\label{f:torus}
+\end{figure}
+
+We create a torus with center $(0; 0; 0)$ and the rotation axis in $z$-direction $(0; 0; 1)$. The projection of a point $x_0$ on the torus is implemented by the following steps:
+\begin{enumerate}
+\item $x_1$ is the projection of $x_0$ on the $xy$-plane
+\item $x_2 = x_1 \frac{r_1}{||x_1||}$. Projection of $x_1$ on the sphere with radius $r_1$ with center $0$. Thereby, $x_2$ is used as the center of a sphere with radius $r_2$.
+\item $x_3 = x_0 - x_2$. Move coordinate system into the center of the sphere with center $x_2$. Thereby, $x_3$ contains the coordinates of $x_0$ in this new coordinate system.
+\item $x_4 = x_3 \frac{r_2}{||x_3||}$. Thereby, $x_4$ is the projection of $x_3$ on the sphere with radius $r_2$.
+\item $x_5 = x_4 + x_2$. Thereby, $x_5$ contains the coordinates of $x_4$ in the original coordinate system. It is the projection of $x_0$ on the torus.
+\end{enumerate}
+
+\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 torus.cc} \\
+\hline
+{\bf parameter file} & \tt init/ & \tt - & \tt - & \tt torus.dat.3d \\
+\hline
+{\bf macro file} & \tt macro/ & \tt - & \tt - & \tt torus.macro.3d \\
+\hline
+{\bf output files} & \tt output/ & \multicolumn{3}{|c|}{\tt torus.mesh/.dat, rotation1.mesh/.dat, rotation2.mesh/dat} \\
+\hline
+\end{tabular}
+\caption{Files of the {\tt torus} example.}
+\label{t:torus files}
+\end{table}
+
+\subsection{Source code}
+
+First, we define the rotation about the $y$-axis, which is used later to rotate the whole torus and the right hand side function.
+
+\begin{lstlisting}{}
+class YRotation
+{
+public:
+ static WorldVector<double>& rotate(WorldVector<double> &x, double angle)
+ {
+ double x0 = x[0] * cos(angle) + x[2] * sin(angle);
+ x[2] = -x[0] * sin(angle) + x[2] * cos(angle);
+ x[0] = x0;
+ return x;
+ }
+};
+\end{lstlisting}
+
+The right hand side function $f$ has to follow the rotation of the torus.
+
+\begin{lstlisting}{}
+class F : public AbstractFunction<double, WorldVector<double> >
+{
+public:
+ MEMORY_MANAGED(F);
+
+ F(int degree)
+ : AbstractFunction<double, WorldVector<double> >(degree),
+ rotation(0.0)
+ {};
+
+ double operator()(const WorldVector<double>& x) const {
+ WorldVector<double> myX = x;
+ YRotation::rotate(myX, -rotation);
+ return = -2.0 * myX[0];
+ }
+
+ void rotate(double r) { rotation += r; };
+
+private:
+ double rotation;
+};
+\end{lstlisting}
+
+Every time, the mesh is rotated, the right hand side function will be informed over the method \verb+rotate+.
+
+Now, we implement the projection on the torus.
+
+\begin{lstlisting}{}
+class TorusProject : public Projection
+{
+public:
+ TorusProject(int id,
+ ProjectionType type,
+ double radius1,
+ double radius2)
+ : Projection(id, type),
+ radius1_(radius1),
+ radius2_(radius2)
+ {};
+
+ virtual ~TorusProject() {};
+
+ void project(WorldVector<double> &x) {
+
+ WorldVector<double> xPlane = x;
+ xPlane[2] = 0.0;
+
+ double norm = std::sqrt(xPlane*xPlane);
+ TEST_EXIT(norm != 0.0)("can't project vector x\n");
+
+ WorldVector<double> center = xPlane;
+ center *= radius1_ / norm;
+
+ x -= center;
+
+ norm = std::sqrt(x*x);
+ TEST_EXIT(norm != 0.0)("can't project vector x\n");
+ x *= radius2_/norm;
+
+ x += center;
+ };
+
+protected:
+ double radius1_;
+ double radius2_;
+};
+\end{lstlisting}
+
+In the main program, we create a torus projection as \verb+VOLUME_PROJECTION+ with ID 1. The values of $r_1$ and $r_2$ are chosen, such that the resulting torus is completely surrounded by the macro mesh that is defined later.
+
+\begin{lstlisting}{}
+int main(int argc, char* argv[])
+{
+ FUNCNAME("torus main");
+
+ // ===== check for init file =====
+ TEST_EXIT(argc == 2)("usage: torus initfile\n");
+
+ // ===== init parameters =====
+ Parameters::init(false, argv[1]);
+
+ // ===== create projection =====
+ double r2 = (1.5 - 1.0/std::sqrt(2.0)) / 2.0;
+ double r1 = 1.0/std::sqrt(2.0) + r2;
+
+ NEW TorusProject(1,
+ VOLUME_PROJECTION,
+ r1,
+ r2);
+
+ ...
+
+ adapt->adapt();
+
+ torus.writeFiles(adaptInfo, true);
+\end{lstlisting}
+
+The problem definition and the creation of the adaptation loop are done in the usual way (here, replaced by \verb+...+) . After the adaptation loop has returned, we write the result.
+
+Before we let the torus rotate, some variables are defined. We set the rotation angle to $\frac{\Pi}{3}$.
+
+\begin{lstlisting}{}
+ double rotation = M_PI/3.0;
+ int i, j;
+ int dim = torus.getMesh()->getDim();
+ int dow = Global::getGeo(WORLD);
+
+ DegreeOfFreedom dof;
+ WorldVector<double> x;
+
+ const FiniteElemSpace *feSpace = torus.getFESpace();
+ const BasisFunction *basFcts = feSpace->getBasisFcts();
+ int numBasFcts = basFcts->getNumber();
+ DegreeOfFreedom *localIndices = GET_MEMORY(DegreeOfFreedom, numBasFcts);
+ DOFAdmin *admin = feSpace->getAdmin();
+
+ WorldVector<DOFVector<double>*> parametricCoords;
+ for(i = 0; i < dow; i++) {
+ parametricCoords[i] = NEW DOFVector<double>(feSpace,
+ "parametric coords");
+ }
+\end{lstlisting}
+
+In the next step, we store the rotated vertex coordinates of the mesh in \verb+parametricCoords+,a vector of DOF vectors, where the first vector stores the first component of each vertex coordinate, and so on. In the STL map \verb+visited+, we store which vertices have already been visited, to avoid multiple rotations of the same point.
+
+\begin{lstlisting}{}
+ std::map<DegreeOfFreedom, bool> visited;
+ TraverseStack stack;
+ ElInfo *elInfo = stack.traverseFirst(torus.getMesh(), -1,
+ Mesh::CALL_LEAF_EL |
+ Mesh::FILL_COORDS);
+ while(elInfo) {
+ basFcts->getLocalIndices(elInfo->getElement(), admin, localIndices);
+ for(i = 0; i < dim + 1; i++) {
+ dof = localIndices[i];
+ x = elInfo->getCoord(i);
+ YRotation::rotate(x, rotation);
+ if(!visited[dof]) {
+ for(j = 0; j < dow; j++) {
+ (*(parametricCoords[j]))[dof] = x[j];
+ }
+ visited[dof] = true;
+ }
+ }
+ elInfo = stack.traverseNext(elInfo);
+ }
+\end{lstlisting}
+
+We create an instance of class \verb+ParametricFirstOrder+ which then is handed to the mesh. Now, in all future mesh traverses the vertex coordinates stored in \verb+parametricCoords+ are returned, instead of the original coordinates.
+
+\begin{lstlisting}{}
+ ParametricFirstOrder parametric(¶metricCoords);
+ torus.getMesh()->setParametric(¶metric);
+\end{lstlisting}
+
+We rotate the right hand side function, reset \verb+adaptInfo+ and start the adaptation loop again. Now, we compute the solution on the rotated torus, which then is written to the files \verb+rotation1.mesh+ and \verb+rotation1.dat+.
+
+\begin{lstlisting}{}
+ f.rotate(rotation);
+ adaptInfo->reset();
+ adapt->adapt();
+
+ DataCollector *dc = NEW DataCollector(feSpace, torus.getSolution());
+ MacroWriter::writeMacro(dc, "output/rotation1.mesh");
+ ValueWriter::writeValues(dc, "output/rotation1.dat");
+ DELETE dc;
+\end{lstlisting}
+
+We perform another rotation. All we have to do is to modify the coordinates in \verb+parametricCoords+ and to inform $f$ about the rotation.
+
+\begin{lstlisting}{}
+ visited.clear();
+ elInfo = stack.traverseFirst(torus.getMesh(), -1,
+ Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
+ while(elInfo) {
+ basFcts->getLocalIndices(elInfo->getElement(), admin, localIndices);
+ for(i = 0; i < dim + 1; i++) {
+ dof = localIndices[i];
+ x = elInfo->getCoord(i);
+ YRotation::rotate(x, rotation);
+ if(!visited[dof]) {
+ for(j = 0; j < dow; j++) {
+ (*(parametricCoords[j]))[dof] = x[j];
+ }
+ visited[dof] = true;
+ }
+ }
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ f.rotate(rotation);
+ adaptInfo->reset();
+ adapt->adapt();
+
+ dc = NEW DataCollector(feSpace, torus.getSolution());
+ MacroWriter::writeMacro(dc, "output/rotation1.mesh");
+ ValueWriter::writeValues(dc, "output/rotation1.dat");
+ DELETE dc;
+\end{lstlisting}
+
+The solution is written to \verb+rotation2.mesh+ and \verb+rotation2.dat+.
+
+Finally, we free some memory and finish the main program.
+
+\begin{lstlisting}{}
+ for(i = 0; i < dow; i++)
+ DELETE parametricCoords[i];
+ FREE_MEMORY(localIndices, DegreeOfFreedom, numBasFcts);
+}
+\end{lstlisting}
+
+\subsection{Parameter file}
+
+In the parameter file, we set the macro file to \verb+./macro/torus_macro.3d+. This two dimensional mesh is 8 times globally refined and successively projected on the torus.
+
+\begin{lstlisting}{}
+dimension of world: 3
+
+torusMesh->macro file name: ./macro/torus_macro.3d
+torusMesh->global refinements: 8
+
+torus->mesh: torusMesh
+torus->dim: 2
+torus->polynomial degree: 1
+
+torus->solver: cg
+torus->solver->max iteration: 1000
+torus->solver->tolerance: 1.e-8
+torus->solver->left precon: diag
+torus->estimator: no
+torus->marker: no
+
+torus->output->filename: output/torus
+torus->output->AMDiS format: 1
+torus->output->AMDiS mesh ext: .mesh
+torus->output->AMDiS data ext: .dat
+\end{lstlisting}
+
+\subsection{Macro file}
+
+ The macro mesh defined in \verb+./macro/torus_macro.3d+ is shown in Figure \ref{f:torus macro}.
+\begin{figure}
+\center
+\includegraphics[width=0.5\textwidth]{fig/torus_macro.jpg}
+\caption{Macro mesh of the torus problem.}
+\label{f:torus macro}
+\end{figure}
+
+\subsection{Output}
+
+In Figure \ref{f:rotated torus}, the solutions on the three tori are shown.
+\begin{figure}
+\center
+\includegraphics[width=0.7\textwidth]{fig/rotated_torus.jpg}
+\caption{Solution on the original torus and on the two rotated tori.}
+\label{f:rotated torus}
+\end{figure}
diff --git a/doc/tutorial/periodic.tex b/doc/tutorial/periodic.tex
new file mode 100644
index 0000000000000000000000000000000000000000..edd3b5a77b65e7c9a16602b402ba0962da21b0c6
--- /dev/null
+++ b/doc/tutorial/periodic.tex
@@ -0,0 +1,153 @@
+\section{Periodic boundary conditions}
+\label{s:periodic}
+
+\begin{figure}
+\center
+\includegraphics[width=0.7\textwidth]{fig/periodic_boundaries.pdf}
+\caption{Two dimensional domain with periodic boundary conditions in both dimensions (left) and in only one dimension (right). In the first case, the solution at $\Omega$ can be propagated to the whole plane of $\mathbb{R}^2$. In the second case, the solution only describes a band within $\mathbb{R}^2$}
+\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.
+
+\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.}
+\label{f:periodic identification}
+\end{figure}
+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.
+\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).
+
+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
+\begin{tabular}{|cc|c|c|c|}
+\hline
+& & {\bf 1d} & {\bf 2d} & {\bf 3d} \\
+\hline
+{\bf source code} & \tt src/ & \multicolumn{3}{|c|}{\tt periodic.cc} \\
+\hline
+{\bf parameter file} & \tt init/ & \tt periodic.dat.1d & \tt periodic.dat.2d & \tt periodic.dat.3d \\
+\hline
+{\bf periodic file} & \tt init/ & \tt periodic.per.1d & \tt periodic.per.2d & \tt periodic.per..3d \\
+\hline
+{\bf macro file} & \tt macro/ & \tt periodic.macro.1d & \tt periodic.macro.2d & \tt periodic.macro.3d \\
+\hline
+{\bf output files} & \tt output/ & \multicolumn{3}{|c|}{\tt periodic.mesh, periodic.dat} \\
+\hline
+\end{tabular}
+\caption{Files of the {\tt periodic} example.}
+\label{t:periodic files}
+\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.
+
+\begin{lstlisting}{}
+periodic.addPeriodicBC(-1);
+\end{lstlisting}
+
+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}.
+
+\begin{lstlisting}{}
+periodicMesh->periodic file: ./init/periodic.per.2d
+\end{lstlisting}
+
+The periodic file \verb+periodic.per.2d+ contains the needed periodic information for the mesh.
+
+\begin{lstlisting}{}
+associations: 2
+
+mode bc el1 - local vertices <-> el2 - local vertices
+ 1 -1 4 1 2 7 2 1
+ 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.
+
+\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}.
+\begin{figure}
+\center
+\includegraphics[width=0.7\textwidth]{fig/periodic_macro.pdf}
+\caption{Macro mesh for the two dimensioanl periodic problem.}
+\label{f:periodic macro}
+\end{figure}
+
+The corresponding file \verb+periodic.macro.2d+ looks like:
+
+\begin{lstlisting}{}
+DIM: 2
+DIM_OF_WORLD: 2
+
+number of elements: 8
+number of vertices: 9
+
+element vertices:
+1 3 0
+3 1 4
+2 4 1
+4 2 5
+4 6 3
+6 4 7
+5 7 4
+7 5 8
+
+element boundaries:
+-1 1 0
+ 0 0 0
+ 0 1 0
+-1 0 0
+-1 0 0
+ 0 1 0
+ 0 0 0
+-1 1 0
+
+vertex coordinates:
+-0.2 -0.5
+ 0.3 -0.5
+ 0.8 -0.5
+-0.2 0.0
+ 0.3 0.0
+ 0.8 0.0
+-0.2 0.5
+ 0.3 0.5
+ 0.8 0.5
+
+element neighbours:
+ 3 -1 1
+ 2 4 0
+ 1 -1 3
+ 0 6 2
+ 7 1 5
+ 6 -1 4
+ 5 3 7
+ 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.
+
+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.
+
+\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$.
+
+\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).}
+\label{f:periodic solution}
+\end{figure}
diff --git a/doc/tutorial/projections.tex b/doc/tutorial/projections.tex
new file mode 100644
index 0000000000000000000000000000000000000000..5da9b03df245e193c800d9279d54075b14378e5e
--- /dev/null
+++ b/doc/tutorial/projections.tex
@@ -0,0 +1,337 @@
+\section{Projections}
+\label{ss: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.
+\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.
+\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.}
+\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.}
+\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.
+
+\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 sphere.cc} \\
+\hline
+{\bf parameter file} & \tt init/ & \tt - & \tt - & \tt sphere.dat.3d \\
+\hline
+{\bf macro file} & \tt macro/ & \tt - & \tt - & \tt sphere.macro.3d \\
+\hline
+{\bf output files} & \tt output/ & \multicolumn{3}{|c|}{\tt sphere.mesh, sphere.dat} \\
+\hline
+\end{tabular}
+\caption{Files of the {\tt sphere} example.}
+\label{t:sphere files}
+\end{table}
+
+\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 ball.cc} \\
+\hline
+{\bf parameter file} & \tt init/ & \tt - & \tt - & \tt ball.dat.3d \\
+\hline
+{\bf macro file} & \tt macro/ & \tt - & \tt - & \tt ball.macro.3d \\
+\hline
+{\bf output files} & \tt output/ & \multicolumn{3}{|c|}{\tt ball.mesh, ball.dat} \\
+\hline
+\end{tabular}
+\caption{Files of the {\tt ball} example.}
+\label{t:ball files}
+\end{table}
+
+\subsection{Source code}
+First we define the projection by implementing a sub class \verb+BallProject+ of the base class\\ \verb+Projection+.
+
+\begin{lstlisting}{}
+class BallProject : public Projection
+{
+public:
+ BallProject(int id,
+ ProjectionType type,
+ WorldVector<double> ¢er,
+ double radius)
+ : Projection(id, type),
+ center_(center),
+ radius_(radius)
+ {};
+
+ void project(WorldVector<double> &x) {
+ x -= center_;
+ double norm = sqrt(x*x);
+ TEST_EXIT(norm != 0.0)("can't project vector x\n");
+ x *= radius_/norm;
+ x += center_;
+ };
+
+protected:
+ WorldVector<double> center_;
+ double radius_;
+};
+\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.
+
+If we compute on the surface, we redefine the function $f$.
+
+\begin{lstlisting}{}
+class F : public AbstractFunction<double, WorldVector<double> >
+{
+public:
+ MEMORY_MANAGED(F);
+
+ F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {};
+
+ double operator()(const WorldVector<double>& x) const
+ {
+ return -2.0 * x[0];
+ }
+};
+\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.
+
+\begin{lstlisting}{}
+ // ===== create projection =====
+ WorldVector<double> ballCenter;
+ ballCenter.set(0.0);
+ NEW BallProject(1,
+ BOUNDARY_PROJECTION,
+ ballCenter,
+ 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.
+
+\begin{lstlisting}{}
+ // ===== create projection =====
+ WorldVector<double> ballCenter;
+ ballCenter.set(0.0);
+ NEW BallProject(1,
+ VOLUME_PROJECTION,
+ ballCenter,
+ 1.0);
+\end{lstlisting}
+
+
+\subsection{Parameter file}
+First, we present the parameter file for the volume projection case (two dimensional mesh).
+
+\begin{lstlisting}{}
+dimension of world: 3
+
+sphereMesh->macro file name: ./macro/sphere_macro.3d
+sphereMesh->global refinements: 10
+
+sphere->mesh: sphereMesh
+sphere->dim: 2
+sphere->polynomial degree: 1
+
+sphere->solver: cg
+sphere->solver->max iteration: 100
+sphere->solver->tolerance: 1.e-8
+sphere->solver->left precon: diag
+
+sphere->estimator: no
+sphere->marker->strategy: 0
+
+sphere->output->filename: output/sphere
+sphere->output->AMDiS format: 1
+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.
+
+Now, we show the parameter file for the boudary projection case (three dimensional mesh).
+
+\begin{lstlisting}{}
+dimension of world: 3
+
+ballMesh->macro file name: ./macro/macro.ball.3d
+ballMesh->global refinements: 15
+
+ball->mesh: ballMesh
+ball->dim: 3
+ball->polynomial degree: 1
+
+ball->solver: cg
+ball->solver->max iteration: 1000
+ball->solver->tolerance: 1.e-8
+ball->solver->left precon: diag
+
+ball->estimator: no
+ball->marker->strategy: 0
+
+ball->output->filename: output/ball
+ball->output->AMDiS format: 1
+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.
+
+\subsection{Macro file}
+
+First, the macro file for the two dimensional mesh.
+
+\begin{lstlisting}{}
+DIM: 2
+DIM_OF_WORLD: 3
+
+number of vertices: 8
+number of elements: 12
+
+vertex coordinates:
+ -1.0 1.0 -1.0
+ 1.0 1.0 -1.0
+ 1.0 1.0 1.0
+ -1.0 1.0 1.0
+ -1.0 -1.0 -1.0
+ 1.0 -1.0 -1.0
+ 1.0 -1.0 1.0
+ -1.0 -1.0 1.0
+
+element vertices:
+ 3 1 0
+ 1 3 2
+ 2 5 1
+ 5 2 6
+ 6 4 5
+ 4 6 7
+ 7 0 4
+ 0 7 3
+ 2 7 6
+ 7 2 3
+ 1 4 0
+ 4 1 5
+
+element boundaries:
+ 0 0 0
+ 0 0 0
+ 0 0 0
+ 0 0 0
+ 0 0 0
+ 0 0 0
+ 0 0 0
+ 0 0 0
+ 0 0 0
+ 0 0 0
+ 0 0 0
+ 0 0 0
+
+projections:
+ 1 0 0
+ 1 0 0
+ 1 0 0
+ 1 0 0
+ 1 0 0
+ 1 0 0
+ 1 0 0
+ 1 0 0
+ 1 0 0
+ 1 0 0
+ 1 0 0
+ 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.
+
+Now, we list the macro file for the three dimensional volume mesh.
+
+\begin{lstlisting}{}
+DIM: 3
+DIM_OF_WORLD: 3
+
+number of vertices: 8
+number of elements: 6
+
+vertex coordinates:
+ -1.0 -1.0 0.0
+ 0.0 -1.0 -1.0
+ 0.0 -1.0 1.0
+ 1.0 -1.0 0.0
+ 0.0 1.0 -1.0
+ 1.0 1.0 0.0
+ -1.0 1.0 0.0
+ 0.0 1.0 1.0
+
+element vertices:
+ 0 5 4 1
+ 0 5 3 1
+ 0 5 3 2
+ 0 5 4 6
+ 0 5 7 6
+ 0 5 7 2
+
+element boundaries:
+ 1 1 0 0
+ 1 1 0 0
+ 1 1 0 0
+ 1 1 0 0
+ 1 1 0 0
+ 1 1 0 0
+
+element neighbours:
+ -1 -1 1 3
+ -1 -1 0 2
+ -1 -1 5 1
+ -1 -1 4 0
+ -1 -1 3 5
+ -1 -1 2 4
+
+projections:
+ 1 1 0 0
+ 1 1 0 0
+ 1 1 0 0
+ 1 1 0 0
+ 1 1 0 0
+ 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.
+
+\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.
+
+\begin{figure}
+\center
+\includegraphics[width=0.32\textwidth]{fig/sphere0.jpg}
+\includegraphics[width=0.32\textwidth]{fig/sphere2.jpg}
+\includegraphics[width=0.32\textwidth]{fig/sphere4.jpg}
+\caption{Surface of a cube successively refinend and projected on the sphere.}
+\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).
+\begin{figure}
+\center
+\begin{tabular}{ccc}
+\includegraphics[width=0.32\textwidth]{fig/sphere_solution.jpg}&
+\includegraphics[width=0.32\textwidth]{fig/sphere_half.jpg}&
+\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).}
+\label{f:sphere solution}
+\end{figure}
diff --git a/doc/tutorial/tutorial.pdf b/doc/tutorial/tutorial.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..993af20f76c9a51b3d60fefad84b4f8b2754e022
Binary files /dev/null and b/doc/tutorial/tutorial.pdf differ
diff --git a/doc/tutorial/tutorial.tex b/doc/tutorial/tutorial.tex
new file mode 100644
index 0000000000000000000000000000000000000000..335f1e7e9532db9593076d0d77cc7add078bf055
--- /dev/null
+++ b/doc/tutorial/tutorial.tex
@@ -0,0 +1,46 @@
+\documentclass[a4paper]{book}
+
+\usepackage{graphicx}
+
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{amsfonts}
+
+\usepackage{listings}
+
+%\usepackage{algorithm}
+%\usepackage{algorithmic}
+
+%\usepackage{color}
+%\definecolor{pink}{rgb}{1, 0.5, 0.5}
+
+% use helvtica as font
+\usepackage{helvet}
+\renewcommand{\rmdefault}{phv}
+
+\newcommand{\myline}{\noindent\rule{\textwidth}{0.3mm}\\}
+
+\usepackage{verbatim}
+
+\usepackage{hyperref}
+
+\voffset =0.mm
+\hoffset =0.mm
+\topmargin =5.mm % beyond 25.mm
+\oddsidemargin =5.mm % beyond 25.mm
+\evensidemargin =5.mm % beyond 25.mm
+\headheight =5.mm
+\headsep =5.mm
+\textheight =225.mm
+\textwidth =150.mm
+
+\begin{document}
+\title{
+\fontsize{30}{35}\selectfont\vspace{-2.9cm}{\bf AMDiS tutorial}}
+\author{Simon Vey, Thomas Witkowski}
+\maketitle
+\tableofcontents
+\input{introduction.tex}
+\input{installation.tex}
+\input{examples.tex}
+\end{document}
diff --git a/doc/tutorial/vecellipt.tex b/doc/tutorial/vecellipt.tex
new file mode 100644
index 0000000000000000000000000000000000000000..d71a9fcef966feecbbd9c3b756d7d5e1c1d8cbc9
--- /dev/null
+++ b/doc/tutorial/vecellipt.tex
@@ -0,0 +1,189 @@
+\section{Systems of PDEs}
+\label{s:vecellipt}
+In this example, we show how to implement a system of coupled PDEs. We define
+\begin{eqnarray}
+-\Delta u &=& f\\
+u - v &=& 0.
+\end{eqnarray}
+For the first equation, we use the boundary condition and definition of function $f$ from Section \ref{s:ellipt}. The second equation defines a second solution component $v$, which is coupled to $u$, such that $v=u$. For the second equation, no boundary conditions have to be defined. The system can be written in matrix-vector form as
+\begin{equation}
+\left( \begin{array}{rr}-\Delta & 0 \\I & -I\end{array}\right)
+\left(\begin{array}{r}u\\v\end{array}\right)
+=
+\left( \begin{array}{r}f\\0\end{array}\right),
+\end{equation}
+where $I$ stands for the identity and $0$ for a {\it zero operator} (or for the absence of any operator). This is a very simple example without practical relevance. But it is appropriate to demonstrate the main principles of implementing vector valued problems.
+
+\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 vecellipt.cc} \\
+\hline
+{\bf parameter file} & \tt init/ & \tt vecellipt.dat.1d & \tt vecellipt.dat.2d & \tt vecellipt.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 vecellipt\_comp<c>.mesh, vecellipt\_comp<c>.dat} \\
+\hline
+\end{tabular}
+\caption{Files of the {\tt vecellipt} example. In the output file names, {\tt <c>} is replaced by the component number.}
+\end{table}
+
+\subsection{Source code}
+Instead of a scalar problem, we now create and initialize the vector valued problem \verb+vecellipt+:
+
+\begin{lstlisting}{}
+ ProblemVec vecellipt("vecellipt");
+ vecellipt.initialize(INIT_ALL);
+\end{lstlisting}
+
+The \verb+AdaptInfo+ constructor is called with the number of problem components, which is defined in the parameter file.
+
+\begin{lstlisting}{}
+ // === create adapt info ===
+ AdaptInfo *adaptInfo = NEW AdaptInfo("vecellipt->adapt",
+ vecellipt.getNumComponents());
+
+ // === create adapt ===
+ AdaptStationary *adapt = NEW AdaptStationary("vecellipt->adapt",
+ &vecellipt,
+ adaptInfo);
+
+
+\end{lstlisting}
+
+The adaptation loop doesn't care about the component number. It treats \verb+vecellipt+ only as implementation of the iteration interface.
+
+The Dirichlet boundary condition for the first equation is defined by
+
+\begin{lstlisting}{}
+ // ===== add boundary conditions =====
+ vecellipt.addDirichletBC(1, 0, NEW G);
+\end{lstlisting}
+
+The first argument is the condition identifier, as in the scalar case. The second argument is the component, the boundary condition belongs to.
+
+The operator definitions for the first equation are:
+
+\begin{lstlisting}{}
+ // ===== create operators =====
+ Operator matrixOperator00(Operator::MATRIX_OPERATOR,
+ vecellipt.getFESpace(0),
+ vecellipt.getFESpace(0));
+ matrixOperator00.addSecondOrderTerm(NEW Laplace_SOT);
+ vecellipt.addMatrixOperator(&matrixOperator00, 0, 0);
+
+ Operator rhsOperator0(Operator::VECTOR_OPERATOR,
+ vecellipt.getFESpace(0));
+
+ int degree = vecellipt.getFESpace(0)->getBasisFcts()->getDegree();
+
+ rhsOperator0.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree)));
+
+ vecellipt.addVectorOperator(&rhsOperator0, 0);
+\end{lstlisting}
+
+Operator \verb+matrixOperator00+ represents the $-\Delta$ operator. Each operator belongs to two finite element spaces, the {\it row space} and the {\it column space}. If an operator has the position $(i,j)$ in the operator matrix, the row space is the finite element space of component $i$ and the column space is the finite element space of component $j$. The finite element spaces can differ in the used basis function degree. The underlying meshes must be the same. After \verb+matrixOperator00+ is created, it is handed to the problems operator matrix at position $(0,0)$.
+The right hand side operator \verb+rhsOperator0+ only needs a row space, which is the finite element space of component $0$ ($u$). It is handed to the operator vector at position $0$.
+
+Now, the operators for the second equation are defined:
+
+\begin{lstlisting}{}
+ Operator matrixOperator10(Operator::MATRIX_OPERATOR,
+ vecellipt.getFESpace(1),
+ vecellipt.getFESpace(0));
+
+ Operator matrixOperator11(Operator::MATRIX_OPERATOR,
+ vecellipt.getFESpace(1),
+ vecellipt.getFESpace(1));
+
+ matrixOperator10.addZeroOrderTerm(NEW Simple_ZOT);
+
+ vecellipt.addMatrixOperator(&matrixOperator10, 1, 0);
+
+ matrixOperator11.addZeroOrderTerm(NEW Simple_ZOT(-1.0));
+
+ vecellipt.addMatrixOperator(&matrixOperator11, 1, 1);
+\end{lstlisting}
+
+Note that the operator \verb+matrixOperator10+ can have different finite element spaces, if the spaces of the two components differ. The operators $I$ and $-I$ are implemented by \verb+Simple_ZOT+, once with a fixed factor of $1$ and once with a factor of $-1$.
+
+\subsection{Parameter file}
+First, the number of components and the basis function degrees are given. We use Lagrange polynomials of degree 1 for the first component and of degree 2 for the second component.
+
+\begin{lstlisting}{}
+vecellipt->components: 2
+
+vecellipt->polynomial degree[0]: 1
+vecellipt->polynomial degree[1]: 2
+\end{lstlisting}
+
+In general, the linear system of equations for systems of PDEs is not symmetric. So with the GMRes solver, we use a solver that doesn't assume symmetric matrices.
+
+\begin{lstlisting}{}
+vecellipt->solver: gmres
+\end{lstlisting}
+
+Note that we have only one solver, because the equations of our system are assembled in one linear system of equations.
+
+Each equation can have its own estimator. In this case, adaptivity should be managed only by the first component. So the second equation has no estimator.
+
+\begin{lstlisting}{}
+vecellipt->estimator[0]: residual
+vecellipt->estimator[1]: no
+\end{lstlisting}
+
+Also the marking strategy can differ between the components. Refinement is done, if at least one component has marked an element for refinement. Coarsening only is done, if all components have marked the element for coarsening. In our example, only component $0$ will mark elements.
+
+\begin{lstlisting}{}
+vecellipt->marker[0]->strategy: 2
+vecellipt->marker[1]->strategy: 0
+\end{lstlisting}
+
+We have only one adaptation loop, which does maximal $6$ iterations. The tolerance can be determined for each component. The total tolerance criterion is fulfilled, if all criteria of all components are fulfilled.
+
+\begin{lstlisting}{}
+vecellipt->adapt->max iteration: 6
+
+vecellipt->adapt[0]->tolerance: 1e-4
+vecellipt->adapt[1]->tolerance: 1e-4
+\end{lstlisting}
+
+Also the output can be controlled for each component individually:
+
+\begin{lstlisting}{}
+vecellipt->output[0]->filename: output/vecellipt_comp0
+
+vecellipt->output[0]->AMDiS format: 1
+vecellipt->output[0]->AMDiS mesh ext: .mesh
+vecellipt->output[0]->AMDiS data ext: .dat
+
+vecellipt->output[1]->filename: output/vecellipt_comp1
+
+vecellipt->output[1]->AMDiS format: 1
+vecellipt->output[1]->AMDiS mesh ext: .mesh
+vecellipt->output[1]->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}
+\label{s:vecellipt output}
+Component $0$ of the solution (approximation of $u$) is written to the files \verb+output/vecellipt0.mesh+ and \verb+output/vecellipt0.dat+.
+Component $1$ of the solution (approximation of $v$) is written to the files \verb+output/vecellipt1.mesh+ and \verb+output/vecellipt1.dat+.
+The two components are visualized in Figure \ref{f:vecellipt}.
+
+\begin{figure}
+\center
+\begin{tabular}{cc}
+\includegraphics[width=0.45\textwidth]{fig/vecellipt0.jpg}&
+\includegraphics[width=0.45\textwidth]{fig/vecellipt1.jpg}\\
+component 0 & component 1
+\end{tabular}
+\caption{The two solution components for $u$ and $v$.}
+\label{f:vecellipt}
+\end{figure}