$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 \leq3$ 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}.
$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 \leq3$ 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}.
@@ -42,20 +50,26 @@ Although the implementation described in Section \ref{s:ellipt code} is dimensi
\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:
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;
using namespace std;
\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.
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+:
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:
...
...
@@ -93,7 +107,12 @@ public:
};
\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.
\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:
...
...
@@ -110,7 +129,17 @@ int main(int argc, char* argv[])
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.
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:
...
...
@@ -120,21 +149,41 @@ Now, a scalar problem with name \verb+ellipt+ is created and initialized:
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 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+:
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);
AdaptInfo adaptInfo("ellipt->adapt");
// === create adapt ===
AdaptStationary *adapt = new AdaptStationary("ellipt->adapt", &ellipt,
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+.
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:
...
...
@@ -143,7 +192,11 @@ Now, we define 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.
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:
...
...
@@ -160,30 +213,41 @@ The operators now are defined as follows:
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.
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();
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+.
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:
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:
In the following, the content of file \verb+init/ellipt.dat.2d+ is
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.
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
...
...
@@ -200,7 +271,10 @@ 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.
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
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.
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
...
...
@@ -218,14 +295,17 @@ 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$.
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$.
After error estimation, elements are marked for refinement and
coarsening. Here, we use the maximum strategy with $\gamma=0.5$.
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.
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
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.
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}
...
...
@@ -252,8 +339,9 @@ The result is written in AMDiS-format to the files \verb+output/ellipt.mesh+ and
\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:
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
...
...
@@ -278,9 +366,11 @@ vertex coordinates:
0.5 0.5
\end{lstlisting}
The first two numbers are interpreted as the coordinates of vertex 0, and so on.
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:
Corresponding to these vertex indices now the four triangles are
given:
\begin{lstlisting}{}
element vertices:
...
...
@@ -290,7 +380,8 @@ element vertices:
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.
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:
...
...
@@ -302,9 +393,17 @@ element boundaries:
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 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.
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:
...
...
@@ -314,7 +413,8 @@ element neighbours:
0 2 -1
\end{lstlisting}
This block is optional. If it isn't given in the macro file, element neighborships are computed automatically.
This block is optional. If it isn't given in the macro file, element
neighborships are computed automatically.
\subsection{Output}
\begin{figure}
...
...
@@ -328,4 +428,10 @@ This block is optional. If it isn't given in the macro file, element neighborshi
\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
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
