diff --git a/demo/init/heat.dat.1d b/demo/init/heat.dat.1d
index f82aa315b51ab9229eed65ee5250f10752ccc72b..df8c3970f9bae2ef36ac6ba2ab7dbc7823c325c2 100644
--- a/demo/init/heat.dat.1d
+++ b/demo/init/heat.dat.1d
@@ -42,7 +42,7 @@ heat->space->marker->ESThetaC: 0.05
heat->space->output->filename: output/heat
heat->space->output->ParaView format: 1
heat->space->output->ParaView animation: 1
-heat->space->output->write every i-th timestep: 1
+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
diff --git a/demo/init/heat.dat.2d b/demo/init/heat.dat.2d
index f49b16a9a98cd86c60c2d5be838aacb5e7ac1fe7..4c865d0aaa99caecddb2f2516cfa75d394d24131 100644
--- a/demo/init/heat.dat.2d
+++ b/demo/init/heat.dat.2d
@@ -42,7 +42,7 @@ heat->space->marker->MSGammaC: 0.1
heat->space->output->filename: output/heat
heat->space->output->ParaView format: 1
heat->space->output->ParaView animation: 1
-heat->space->output->write every i-th timestep: 1
+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
diff --git a/demo/init/heat.dat.3d b/demo/init/heat.dat.3d
index 77dbfc0f9b0be2e8adbed80f16ed69d7d4dfe415..986c01f6bfb107836444ea16f8c7692e1da06fad 100644
--- a/demo/init/heat.dat.3d
+++ b/demo/init/heat.dat.3d
@@ -43,7 +43,7 @@ heat->adapt->coarsen bisections: 3
heat->space->output->filename: output/heat
heat->space->output->ParaView format: 1
heat->space->output->ParaView animation: 1
-heat->space->output->write every i-th timestep: 1
+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
diff --git a/demo/src/heat.cc b/demo/src/heat.cc
index 04829495670ef480e0dc18c6fd815967bb8efc55..c889d3ae843f99acd41b76e8f3060a2b93bf6339 100644
--- a/demo/src/heat.cc
+++ b/demo/src/heat.cc
@@ -67,9 +67,7 @@ public:
void setTime(AdaptInfo *adaptInfo)
{
ProblemInstat::setTime(adaptInfo);
-
rhsTime = adaptInfo->getTime() - (1 - theta) * adaptInfo->getTimestep();
- boundaryTime = adaptInfo->getTime();
}
void closeTimestep(AdaptInfo *adaptInfo)
@@ -120,17 +118,11 @@ public:
}
/// Returns pointer to \ref rhsTime.
- double *getRHSTimePtr()
+ double *getRhsTimePtr()
{
return &rhsTime;
}
- /// Returns pointer to \ref theta1.
- double *getBoundaryTimePtr()
- {
- return &boundaryTime;
- }
-
private:
/// Used for theta scheme.
double theta;
@@ -141,9 +133,6 @@ private:
/// time for right hand side functions.
double rhsTime;
- /// time for boundary functions.
- double boundaryTime;
-
/// Pointer to boundary function. Needed for initial problem.
AbstractFunction<double, WorldVector<double> > *exactSolution;
};
@@ -187,7 +176,7 @@ int main(int argc, char** argv)
// ===== create boundary functions =====
G *boundaryFct = new G;
- boundaryFct->setTimePtr(heat.getBoundaryTimePtr());
+ boundaryFct->setTimePtr(heat.getTime());
heat.setExactSolution(boundaryFct);
heatSpace.addDirichletBC(1, boundaryFct);
@@ -195,7 +184,7 @@ int main(int argc, char** argv)
// ===== create rhs functions =====
int degree = heatSpace.getFeSpace()->getBasisFcts()->getDegree();
F *rhsFct = new F(degree);
- rhsFct->setTimePtr(heat.getRHSTimePtr());
+ rhsFct->setTimePtr(heat.getRhsTimePtr());
// ===== create operators =====
diff --git a/doc/tutorial/couple.tex b/doc/tutorial/couple.tex
index 48404b289f2959c6777cb8a870ba484e5afb758a..0aedcfc6de8b6ce6e51cc1c91513613d5a2e0a7a 100644
--- a/doc/tutorial/couple.tex
+++ b/doc/tutorial/couple.tex
@@ -190,11 +190,11 @@ The operators for the first problem are defined like in Section
\begin{lstlisting}{}
// ===== create operators for problem2 =====
- Operator matrixOperator2(problem2.getFESpace());
+ Operator matrixOperator2(problem2.getFeSpace());
matrixOperator2.addZeroOrderTerm(new Simple_ZOT);
problem2.addMatrixOperator(&matrixOperator2);
- Operator rhsOperator2(problem2.getFESpace());
+ Operator rhsOperator2(problem2.getFeSpace());
rhsOperator2.addZeroOrderTerm(new VecAtQP_ZOT(problem1.getSolution(),
new Identity(degree)));
problem2.addVectorOperator(&rhsOperator2);
diff --git a/doc/tutorial/ellipt.tex b/doc/tutorial/ellipt.tex
index 8993682a7182aeb17137a771a5013734ac4df62d..218ac096d79a56760ddc92d09d5b0f7d149d37d5 100644
--- a/doc/tutorial/ellipt.tex
+++ b/doc/tutorial/ellipt.tex
@@ -202,13 +202,13 @@ The operators now are defined as follows:
\begin{lstlisting}{}
// ===== create matrix operator =====
- Operator matrixOperator(ellipt.getFESpace());
+ Operator matrixOperator(ellipt.getFeSpace());
matrixOperator.addSecondOrderTerm(new Laplace_SOT);
ellipt.addMatrixOperator(&matrixOperator);
// ===== create rhs operator =====
- int degree = ellipt.getFESpace()->getBasisFcts()->getDegree();
- Operator rhsOperator(ellipt.getFESpace());
+ int degree = ellipt.getFeSpace()->getBasisFcts()->getDegree();
+ Operator rhsOperator(ellipt.getFeSpace());
rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
ellipt.addVectorOperator(&rhsOperator);
\end{lstlisting}
diff --git a/doc/tutorial/heat.tex b/doc/tutorial/heat.tex
index 0728c13787b363470173d4caca299198450d830a..cac2cb35734fa6502bbe694d418dae19ff5755de 100644
--- a/doc/tutorial/heat.tex
+++ b/doc/tutorial/heat.tex
@@ -38,13 +38,6 @@ to the right hand side, the equation reads
\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|}
@@ -64,7 +57,17 @@ to the right hand side, the equation reads
\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:
+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> >,
@@ -78,7 +81,10 @@ public:
};
\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.
+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
@@ -86,11 +92,26 @@ The variable $timePtr$ is a base class member of \verb+TimedObject+. This pointe
\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:
+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}
+\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:
@@ -108,15 +129,21 @@ public:
}
\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 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)
{
+ ProblemInstat::setTime(adaptInfo);
rhsTime = adaptInfo->getTime() - (1 - theta) * adaptInfo->getTimestep();
- boundaryTime = adaptInfo->getTime();
}
void closeTimestep(AdaptInfo *adaptInfo)
@@ -127,10 +154,12 @@ The next lines show the implementation of the time interface.
\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}$, which is the
-current time.
+the problem about the current time. When this function is
+reimplemented, one should always call the function in the parent
+class, such that all time relavant variables in \verb+ProblemInstat+
+will be updated. 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}$, which is the current time.
The method \verb+closeTimestep+ is called at the end of each timestep
by the adaptation loop. In the default implementation of
@@ -190,13 +219,11 @@ Now, we define some getting functions and the private member variables:
double *getThetaPtr() { return θ };
double *getTheta1Ptr() { return &theta1; };
double *getRHSTimePtr() { return &rhsTime; };
- double *getBoundaryTimePtr() { return &boundaryTime; };
private:
double theta;
double theta1;
double rhsTime;
- double boundaryTime;
AbstractFunction<double, WorldVector<double> > *exactSolution;
};
\end{lstlisting}
@@ -212,6 +239,7 @@ int main(int argc, char** argv)
// ===== init parameters =====
Parameters::init(false, argv[1]);
+ Parameters::readArgv(argc, argv);
// ===== create and init stationary problem =====
ProblemScal heatSpace("heat->space");
@@ -257,14 +285,14 @@ The functions $f$ and $g$ are declared in the following way:
\begin{lstlisting}{}
// ===== create boundary functions =====
G *boundaryFct = new G;
- boundaryFct->setTimePtr(heat.getBoundaryTimePtr());
+ boundaryFct->setTimePtr(heat.getTime());
heat.setExactSolution(boundaryFct);
heatSpace.addDirichletBC(1, boundaryFct);
// ===== create rhs functions =====
- int degree = heatSpace.getFESpace()->getBasisFcts()->getDegree();
+ int degree = heatSpace.getFeSpace()->getBasisFcts()->getDegree();
F *rhsFct = new F(degree);
- rhsFct->setTimePtr(heat.getRHSTimePtr());
+ 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
@@ -279,7 +307,7 @@ Now, we define the operators:
double zero = 0.0;
// create laplace
- Operator A(heatSpace.getFESpace());
+ Operator A(heatSpace.getFeSpace());
A.addSecondOrderTerm(new Laplace_SOT);
A.setUhOld(heat.getOldSolution());
if (*(heat.getThetaPtr()) != 0.0)
@@ -301,7 +329,7 @@ operator by \verb+setUhOld+.
\begin{lstlisting}{}
// create zero order operator
- Operator C(heatSpace.getFESpace());
+ Operator C(heatSpace.getFeSpace());
C.addZeroOrderTerm(new Simple_ZOT);
C.setUhOld(heat.getOldSolution());
heatSpace.addMatrixOperator(C, heat.getInvTau(), heat.getInvTau());
@@ -320,7 +348,7 @@ and the adaptation loop is started:
\begin{lstlisting}{}
// create RHS operator
- Operator F(heatSpace.getFESpace());
+ Operator F(heatSpace.getFeSpace());
F.addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct));
heatSpace.addVectorOperator(F);
@@ -356,39 +384,20 @@ heat->adapt->start time: 0.0
heat->adapt->end time: 1.0
\end{lstlisting}
-Now, tolerances are determined:
+Now, tolerances for the space and the time error 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
+heat->adapt->tolerance: 0.05
+heat->adapt->time tolerance: 0.05
\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+.
+strategy.
The following lines determine, whether coarsening is allowed in
regions with sufficient small errors, and how many refinements or
@@ -402,33 +411,34 @@ heat->adapt->coarsen bisections: 2
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
+heat->space->output->filename: output/heat
+heat->space->output->ParaView format: 1
+heat->space->output->ParaView animation: 1
+heat->space->output->write every i-th timestep: 1
+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.
+\verb+output/heat+ and end with the extension \verb+.vtu+. 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 to the file
+\verb+output/heat00.000.vtu+. If the parameter
+\verb+ParaView animation+ is enabled, AMDiS writes for the whole
+simulation one ParaView \verb+pvd+ file (in this case
+\verb+output/heat.pvd+) including the names of all \verb+vtu+ files
+that were created. This file makes it very easy to view and analyze
+all the results of an instationary problem in ParaView.
+
+Finally, we set parameter \verb+WAIT+ to $0$. If the variable is set
+to $1$, 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
+WAIT: 0
\end{lstlisting}
\subsection{Macro file}
@@ -437,8 +447,7 @@ 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
+\verb+output/heat00.000.vtu+. 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.
diff --git a/doc/tutorial/nonlin.tex b/doc/tutorial/nonlin.tex
deleted file mode 100644
index 401ab7672e18c3cb96ae28c46a69f3bb59c0c4ee..0000000000000000000000000000000000000000
--- a/doc/tutorial/nonlin.tex
+++ /dev/null
@@ -1,491 +0,0 @@
-\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:
- double operator()(const WorldVector<double>& x) const
- {
- return 0.0;
- }
-};
-
-class F : public AbstractFunction<double, WorldVector<double> >
-{
-public:
- /// Constructor
- F(int degree)
- : AbstractFunction<double, WorldVector<double> >(degree)
- {}
-
- /// 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:
- X3() : AbstractFunction<double, double>(3) {}
-
- /// 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(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(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(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/parametric.tex b/doc/tutorial/parametric.tex
index 3ee068006ac8f18fc1a2c7b41e82e2a5785cbb8d..f7febe27297bb406282af5bf84a23eaf4002958c 100644
--- a/doc/tutorial/parametric.tex
+++ b/doc/tutorial/parametric.tex
@@ -170,7 +170,7 @@ Before we let the torus rotate, some variables are defined. We set the rotation
DegreeOfFreedom dof;
WorldVector<double> x;
- const FiniteElemSpace *feSpace = torus.getFESpace();
+ const FiniteElemSpace *feSpace = torus.getFeSpace();
const BasisFunction *basFcts = feSpace->getBasisFcts();
int numBasFcts = basFcts->getNumber();
DegreeOfFreedom *localIndices = new double[numBasFcts];
diff --git a/doc/tutorial/tutorial.pdf b/doc/tutorial/tutorial.pdf
index 613e311c722c4b1827ba5aec7526e580f80c33d1..5590fcc562874b813585f67d244f3f09f233b972 100644
Binary files a/doc/tutorial/tutorial.pdf and b/doc/tutorial/tutorial.pdf differ
diff --git a/doc/tutorial/vecellipt.tex b/doc/tutorial/vecellipt.tex
index a03045a2a1e73dffc590003863a76a024e60826e..035c0b603dff07c0e5f44541be1337ffec96c639 100644
--- a/doc/tutorial/vecellipt.tex
+++ b/doc/tutorial/vecellipt.tex
@@ -79,12 +79,12 @@ The operator definitions for the first equation are:
\begin{lstlisting}{}
// ===== create operators =====
- Operator matrixOperator00(vecellipt.getFESpace(0), vecellipt.getFESpace(0));
+ Operator matrixOperator00(vecellipt.getFeSpace(0), vecellipt.getFeSpace(0));
matrixOperator00.addSecondOrderTerm(new Laplace_SOT);
vecellipt.addMatrixOperator(&matrixOperator00, 0, 0);
- int degree = vecellipt.getFESpace(0)->getBasisFcts()->getDegree();
- Operator rhsOperator0(vecellipt.getFESpace(0));
+ int degree = vecellipt.getFeSpace(0)->getBasisFcts()->getDegree();
+ Operator rhsOperator0(vecellipt.getFeSpace(0));
rhsOperator0.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
vecellipt.addVectorOperator(&rhsOperator0, 0);
\end{lstlisting}
@@ -105,11 +105,11 @@ operator vector at position $0$.
Now, the operators for the second equation are defined:
\begin{lstlisting}{}
- Operator matrixOperator10(vecellipt.getFESpace(1), vecellipt.getFESpace(0));
+ Operator matrixOperator10(vecellipt.getFeSpace(1), vecellipt.getFeSpace(0));
matrixOperator10.addZeroOrderTerm(new Simple_ZOT);
vecellipt.addMatrixOperator(matrixOperator10, 1, 0);
- Operator matrixOperator11(vecellipt.getFESpace(1), vecellipt.getFESpace(1));
+ Operator matrixOperator11(vecellipt.getFeSpace(1), vecellipt.getFeSpace(1));
matrixOperator11.addZeroOrderTerm(new Simple_ZOT(-1.0));
vecellipt.addMatrixOperator(matrixOperator11, 1, 1);
\end{lstlisting}