* Added Tutorial

doc/tutorial/heat.tex

+\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 &theta; };
+  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}