Skip to content
Snippets Groups Projects
Commit 2f4c9a9a authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Changed tutorial due to remove of MATRIX_OPERATOR and VECTOR_OPERATOR flags.

parent ffd0307d
No related branches found
No related tags found
No related merge requests found
......@@ -184,22 +184,28 @@ Now, we create \verb+problem2+. It should have its own finite element space, sys
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+.
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());
Operator matrixOperator2(problem2.getFESpace());
matrixOperator2.addZeroOrderTerm(new Simple_ZOT);
problem2.addMatrixOperator(&matrixOperator2);
Operator rhsOperator2(Operator::VECTOR_OPERATOR, problem2.getFESpace());
Operator rhsOperator2(problem2.getFESpace());
rhsOperator2.addZeroOrderTerm(new VecAtQP_ZOT(problem1.getSolution(),
new Identity(degree)));
problem2.addVectorOperator(&rhsOperator2);
\end{lstlisting}
At the left hand side, we have just an ordinary \verb+Simple_ZOT+. At the right hand side, we have a zero order term of the form $f(u)$ with $f=I$ the identity. $u$ is given by the solution DOF vector of \verb+problem1+. $I$ maps the values of the DOF vector evaluated at quadrature points to itself. The function degree is used to determine the needed quadrature degree in the assembler.
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:
......@@ -214,9 +220,12 @@ Now, the adaptation loop is created:
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.
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.
The adaptation loop is now started. After it is finished, the
solutions of both problems are written.
\begin{lstlisting}{}
// ===== start adaptation loop =====
......@@ -238,7 +247,9 @@ 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:
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
......@@ -266,7 +277,10 @@ 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:
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
......@@ -288,9 +302,13 @@ 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}.
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}.
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}.
......@@ -202,13 +202,13 @@ The operators now are defined as follows:
\begin{lstlisting}{}
// ===== create matrix operator =====
Operator matrixOperator(Operator::MATRIX_OPERATOR, 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(Operator::VECTOR_OPERATOR, ellipt.getFESpace());
Operator rhsOperator(ellipt.getFESpace());
rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
ellipt.addVectorOperator(&rhsOperator);
\end{lstlisting}
......
\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
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)\\
......@@ -9,18 +10,30 @@ 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
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
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
$\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}
......@@ -114,11 +127,31 @@ 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}$. $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.
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)
......@@ -137,7 +170,16 @@ Now, the implementation of the \verb+ProblemStatBase+ interface begins. As menti
}
\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:
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)
......@@ -165,7 +207,8 @@ private:
};
\end{lstlisting}
The definition of class \verb+Heat+ is now finished. In the following, the main program is described.
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)
......@@ -185,9 +228,15 @@ int main(int argc, char** argv)
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.
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:
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
......@@ -204,7 +253,11 @@ The next step is the creation of the needed \verb+AdaptInfo+ objects and of the
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 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:
......@@ -222,7 +275,10 @@ The definitions of functions $f$ and $g$ are:
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+.
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:
......@@ -232,8 +288,7 @@ Now, we define the operators:
double zero = 0.0;
// create laplace
Operator A(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR,
heatSpace.getFESpace());
Operator A(heatSpace.getFESpace());
A.addSecondOrderTerm(new Laplace_SOT);
A.setUhOld(heat.getOldSolution());
if (*(heat.getThetaPtr()) != 0.0)
......@@ -255,8 +310,7 @@ operator by \verb+setUhOld+.
\begin{lstlisting}{}
// create zero order operator
Operator C(Operator::MATRIX_OPERATOR | Operator::VECTOR_OPERATOR,
heatSpace.getFESpace());
Operator C(heatSpace.getFESpace());
C.addZeroOrderTerm(new Simple_ZOT);
C.setUhOld(heat.getOldSolution());
heatSpace.addMatrixOperator(C, heat.getTau1Ptr(), heat.getTau1Ptr());
......@@ -273,7 +327,7 @@ and the adaptation loop is started:
\begin{lstlisting}{}
// create RHS operator
Operator F(Operator::VECTOR_OPERATOR, heatSpace.getFESpace());
Operator F(heatSpace.getFESpace());
F.addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct));
heatSpace.addVectorOperator(F);
......
......@@ -12,7 +12,12 @@ 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
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)
......@@ -361,9 +366,26 @@ int main(int argc, char* argv[])
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.
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}
......@@ -378,10 +400,7 @@ We have to add operators representing the Newton step equation $-\Delta d + 4u_n
double zero = 0.0;
double minusOne = -1.0;
Operator *nonlinOperator0 = new Operator(Operator::MATRIX_OPERATOR |
Operator::VECTOR_OPERATOR,
nonlin.getFESpace());
Operator *nonlinOperator0 = new Operator(nonlin.getFESpace());
nonlinOperator0->setUhOld(nonlin.getSolution());
nonlinOperator0->addZeroOrderTerm(new VecAtQP_ZOT(nonlin.getSolution(),
new X3));
......@@ -389,10 +408,7 @@ We have to add operators representing the Newton step equation $-\Delta d + 4u_n
nonlin.addMatrixOperator(nonlinOperator0, &four, &one);
nonlin.addVectorOperator(nonlinOperator0, &one, &zero);
Operator *nonlinOperator2 = new Operator(Operator::MATRIX_OPERATOR |
Operator::VECTOR_OPERATOR,
nonlin.getFESpace());
Operator *nonlinOperator2 = new Operator(nonlin.getFESpace());
nonlinOperator2->setUhOld(nonlin.getSolution());
nonlinOperator2->addSecondOrderTerm(new Laplace_SOT);
......@@ -401,8 +417,7 @@ We have to add operators representing the Newton step equation $-\Delta d + 4u_n
int degree = nonlin.getFESpace()->getBasisFcts()->getDegree();
Operator* rhsFunctionOperator = new Operator(Operator::VECTOR_OPERATOR,
nonlin.getFESpace());
Operator* rhsFunctionOperator = new Operator(nonlin.getFESpace());
rhsFunctionOperator->addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
nonlin.addVectorOperator(rhsFunctionOperator, &minusOne, &one);
......
No preview for this file type
......@@ -79,15 +79,12 @@ The operator definitions for the first equation are:
\begin{lstlisting}{}
// ===== create operators =====
Operator matrixOperator00(Operator::MATRIX_OPERATOR,
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(Operator::VECTOR_OPERATOR,
vecellipt.getFESpace(0));
Operator rhsOperator0(vecellipt.getFESpace(0));
rhsOperator0.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
vecellipt.addVectorOperator(&rhsOperator0, 0);
\end{lstlisting}
......@@ -108,13 +105,11 @@ 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 matrixOperator10(vecellipt.getFESpace(1), vecellipt.getFESpace(0));
matrixOperator10.addZeroOrderTerm(new Simple_ZOT);
vecellipt.addMatrixOperator(matrixOperator10, 1, 0);
Operator matrixOperator11(Operator::MATRIX_OPERATOR,
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}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment