ellipt.tex 14.6 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1
2
3
4
5
6
7
8
9
10
11
12
13
\section{Stationary problem with Dirichlet boundary condition}
\label{s:ellipt}
As example for a stationary problem, we choose the Poisson equation
\begin{eqnarray}
\label{eq:ellipt}
-\Delta u &=& f ~~~\mbox{in } \Omega \subset \mathbb{R}^{dim}\\
u &=& g ~~~ \mbox{on } \partial\Omega
\end{eqnarray}
with
\begin{eqnarray}
f(x) &=& -\left( 400x^2 - 20 dow \right) e^{-10x^2}\\
g(x) &=& e^{-10x^2}.
\end{eqnarray}
Thomas Witkowski's avatar
Thomas Witkowski committed
14
15
16
17
18
19
20
21
22
23
24
$dim$ is the problem dimension and thus the dimension of the mesh the
problem will be discretized on. $dow$ is the dimension of the world
the problem lives in. So world coordinates are always real valued
vectors of size $dow$. Note that the problem is defined in a dimension
independent way. Furthermore, $dow$ has not to be equal to $dim$ as
long as $1 \leq dim \leq dow \leq 3$ holds.

Although the implementation described in Section \ref{s:ellipt code}
is dimension independent, we focus on the case $dim = dow = 2$ for the
rest of this section. The analytical solution on $\Omega = [0,1]
\times [0,1]$ is shown in Figure \ref{f:ellipt solution}.
Thomas Witkowski's avatar
Thomas Witkowski committed
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
\begin{figure}
\center
\includegraphics[width=0.3\textwidth]{fig/ellipt.jpg}
\caption{Solution of the Poisson equation on the unit square.}
\label{f:ellipt solution}
\end{figure}

\begin{table}
\center
\begin{tabular}{|cc|c|c|c|}
\hline 
& & {\bf 1d} & {\bf 2d} & {\bf 3d} \\
\hline 
{\bf source code} & \tt src/ & \multicolumn{3}{|c|}{\tt ellipt.cc} \\
\hline 
{\bf parameter file} & \tt init/ & \tt ellipt.dat.1d & \tt ellipt.dat.2d & \tt ellipt.dat.3d \\
\hline 
{\bf macro file} & \tt macro/ & \tt macro.stand.1d & \tt macro.stand.2d & \tt macro.stand.3d \\
\hline 
{\bf output files} & \tt output/ & \multicolumn{3}{|c|}{\tt ellipt.mesh, ellipt.dat} \\
\hline 
\end{tabular}
\caption{Files of the {\tt ellipt} example.}
\label{t:ellipt}
\end{table}

\subsection{Source code}
\label{s:ellipt code}
Thomas Witkowski's avatar
Thomas Witkowski committed
53
54
55
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:
Thomas Witkowski's avatar
Thomas Witkowski committed
56
57
58
59
 
\begin{lstlisting}{}
#include "AMDiS.h"
using namespace AMDiS;
Thomas Witkowski's avatar
Thomas Witkowski committed
60
using namespace std;
Thomas Witkowski's avatar
Thomas Witkowski committed
61
62
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
63
64
65
66
67
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.
Thomas Witkowski's avatar
Thomas Witkowski committed
68

Thomas Witkowski's avatar
Thomas Witkowski committed
69
70
Now, the functions $f$ and $g$ will be defined by the classes \verb+F+
and \verb+G+:
Thomas Witkowski's avatar
Thomas Witkowski committed
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97

\begin{lstlisting}{}
class G : public AbstractFunction<double, WorldVector<double> >
{
public:
  double operator()(const WorldVector<double>& x) const
  {
    return exp(-10.0 * (x * x));
  }
};
\end{lstlisting}

\verb+G+ is a sub class of the templated class
\verb+AbstractFunction<R, T>+ that represents a mapping from type
\verb+T+ to type \verb+R+. Here, we want to define a mapping from
$\mathbb{R}^{dow}$, implemented by the class
\verb+WorldVector<double>+, to $\mathbb{R}$, represented by the data
type \verb+double+. The actual mapping is defined by overloading the
\verb+operator()+. \verb+x*x+ stands for the scalar product of vector
\verb+x+ with itself.

The class \verb+F+ is defined in a similar way:

\begin{lstlisting}{}
class F : public AbstractFunction<double, WorldVector<double> >
{
public:
98
  F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}
Thomas Witkowski's avatar
Thomas Witkowski committed
99

100
101
  double operator()(const WorldVector<double>& x) const 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
102
103
104
105
106
107
108
109
    int dow = Global::getGeo(WORLD);
    double r2 = (x * x);
    double ux = exp(-10.0 * r2);
    return -(400.0 * r2 - 20.0 * dow) * ux;
  }
};
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
110
111
112
113
114
115
\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.
Thomas Witkowski's avatar
Thomas Witkowski committed
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

Now, we start with the main program:

\begin{lstlisting}{}
// ===== main program =====
int main(int argc, char* argv[])
{
  FUNCNAME("main");

  // ===== check for init file =====
  TEST_EXIT(argc == 2)("usage: ellipt initfile\n");

  // ===== init parameters =====
  Parameters::init(true, argv[1]);
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
132
133
134
135
136
137
138
139
140
141
142
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.
Thomas Witkowski's avatar
Thomas Witkowski committed
143
144
145
146
147
148
149
150
151

Now, a scalar problem with name \verb+ellipt+ is created and initialized:  

\begin{lstlisting}{}
  // ===== create and init the scalar problem ===== 
  ProblemScal ellipt("ellipt");
  ellipt.initialize(INIT_ALL);
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
152
153
154
155
156
157
158
159
160
161
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}).
Thomas Witkowski's avatar
Thomas Witkowski committed
162

Thomas Witkowski's avatar
Thomas Witkowski committed
163
164
The next steps are the creation of the adaptation loop and the
corresponding \verb+AdaptInfo+:
Thomas Witkowski's avatar
Thomas Witkowski committed
165
166
167

\begin{lstlisting}{}
  // === create adapt info ===
Thomas Witkowski's avatar
Thomas Witkowski committed
168
  AdaptInfo adaptInfo("ellipt->adapt");
Thomas Witkowski's avatar
Thomas Witkowski committed
169
170

  // === create adapt ===
Thomas Witkowski's avatar
Thomas Witkowski committed
171
  AdaptStationary adapt("ellipt->adapt", ellipt, adaptInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
172
173
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
174
175
176
177
178
179
180
181
182
183
184
185
186
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+.
Thomas Witkowski's avatar
Thomas Witkowski committed
187
188
189
190
191

Now, we define boundary conditions:

\begin{lstlisting}{}
  // ===== add boundary conditions =====
192
  ellipt.addDirichletBC(1, new G);
Thomas Witkowski's avatar
Thomas Witkowski committed
193
194
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
195
196
197
198
199
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.
Thomas Witkowski's avatar
Thomas Witkowski committed
200
201
202
203
204

The operators now are defined as follows:

\begin{lstlisting}{}
  // ===== create matrix operator =====
205
  Operator matrixOperator(ellipt.getFESpace());
206
  matrixOperator.addSecondOrderTerm(new Laplace_SOT);
Thomas Witkowski's avatar
Thomas Witkowski committed
207
208
209
210
  ellipt.addMatrixOperator(&matrixOperator);

  // ===== create rhs operator =====
  int degree = ellipt.getFESpace()->getBasisFcts()->getDegree();
211
  Operator rhsOperator(ellipt.getFESpace());
212
  rhsOperator.addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
Thomas Witkowski's avatar
Thomas Witkowski committed
213
214
215
  ellipt.addVectorOperator(&rhsOperator);
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
216
217
218
219
220
221
222
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.
Thomas Witkowski's avatar
Thomas Witkowski committed
223
224
225
226
227

Finally we start the adaptation loop and afterwards write out the results:

\begin{lstlisting}{}
  // ===== start adaption loop =====
Thomas Witkowski's avatar
Thomas Witkowski committed
228
  adapt.adapt();
Thomas Witkowski's avatar
Thomas Witkowski committed
229
230
231
232
233
234

  // ===== write result =====
  ellipt.writeFiles(adaptInfo, true);
}
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
235
236
237
238
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+.
Thomas Witkowski's avatar
Thomas Witkowski committed
239
240
241

\subsection{Parameter file}
\label{s:ellipt parameter}
Thomas Witkowski's avatar
Thomas Witkowski committed
242
243
The name of the parameter file must be given as command line
argument. In the 2d case we call:
Thomas Witkowski's avatar
Thomas Witkowski committed
244
245
246
247
248

\begin{lstlisting}{}
> ./ellipt init/ellipt.dat.2d
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
249
250
In the following, the content of file \verb+init/ellipt.dat.2d+ is
described:
Thomas Witkowski's avatar
Thomas Witkowski committed
251
252
253
254
255
256
257
258

\begin{lstlisting}{}
dimension of world:                2

elliptMesh->macro file name:       ./macro/macro.stand.2d
elliptMesh->global refinements:    0
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
259
260
261
262
263
264
265
266
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.
Thomas Witkowski's avatar
Thomas Witkowski committed
267
268
269
270

\begin{lstlisting}{}
ellipt->mesh:                      elliptMesh
ellipt->dim:                       2
271
ellipt->polynomial degree:         1
Thomas Witkowski's avatar
Thomas Witkowski committed
272
273
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
274
275
276
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
277
basis functions of degree 1.
Thomas Witkowski's avatar
Thomas Witkowski committed
278
279

\begin{lstlisting}{}
280
ellipt->solver:                    cg   
Thomas Witkowski's avatar
Thomas Witkowski committed
281
282
ellipt->solver->max iteration:     1000
ellipt->solver->tolerance:         1.e-8
283
284
ellipt->solver->left precon:       diag
ellipt->solver->right precon:      no
Thomas Witkowski's avatar
Thomas Witkowski committed
285
286
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
287
288
We use the {\it conjugate gradient method} as iterative solver. The
solving process stops after maximal $1000$ iterations or when a
289
290
tolerance of $10^{-8}$ is reached. Furthermore, we apply diagonal left
preconditioning, and no right preconditioning.
Thomas Witkowski's avatar
Thomas Witkowski committed
291
292

\begin{lstlisting}{}
293
ellipt->estimator:                 residual
Thomas Witkowski's avatar
Thomas Witkowski committed
294
295
296
297
298
ellipt->estimator->error norm:     1
ellipt->estimator->C0:             0.1
ellipt->estimator->C1:             0.1
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
299
300
301
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$.
Thomas Witkowski's avatar
Thomas Witkowski committed
302
303
304
305
306
307

\begin{lstlisting}{}
ellipt->marker->strategy:          2     % 0: no 1: GR 2: MS 3: ES 4:GERS
ellipt->marker->MSGamma:           0.5
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
308
309
After error estimation, elements are marked for refinement and
coarsening. Here, we use the maximum strategy with $\gamma = 0.5$.
Thomas Witkowski's avatar
Thomas Witkowski committed
310
311
312

\begin{lstlisting}{}
ellipt->adapt->tolerance:          1e-4
313
ellipt->adapt->max iteration:      10
Thomas Witkowski's avatar
Thomas Witkowski committed
314
315
316
ellipt->adapt->refine bisections:  2
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
317
The adaptation loop stops, when an error tolerance of $10^{-4}$ is
318
reached, or after maximal $10$ iterations. An element that is marked
Thomas Witkowski's avatar
Thomas Witkowski committed
319
320
321
for refinement, is bisected twice within one iteration. Analog
elements that are marked for coarsening are coarsened twice per
iteration.
Thomas Witkowski's avatar
Thomas Witkowski committed
322
323
324

\begin{lstlisting}{}
ellipt->output->filename:          output/ellipt
325
ellipt->output->ParaView format:   1
Thomas Witkowski's avatar
Thomas Witkowski committed
326
327
\end{lstlisting}

328
329
The result is written in ParaView-format to the file
\verb+output/ellipt.vtu+.
Thomas Witkowski's avatar
Thomas Witkowski committed
330
331
332
333
334
335
336
337
338

\subsection{Macro file}
\label{s:ellipt macro}
\begin{figure}
\center
\includegraphics[width=0.4\textwidth]{fig/ellipt_macro.pdf}
\caption{Two dimensional macro mesh}
\label{f:ellipt macro}
\end{figure}
Thomas Witkowski's avatar
Thomas Witkowski committed
339
340
341
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:
Thomas Witkowski's avatar
Thomas Witkowski committed
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365

\begin{lstlisting}{}
DIM: 2
DIM_OF_WORLD: 2
\end{lstlisting}

Then the total number of elements and vertices are given:

\begin{lstlisting}{}
number of elements: 4  
number of vertices: 5  
\end{lstlisting}

The next block describes the two dimensional coordinates of the five vertices:  

\begin{lstlisting}{}
vertex coordinates:
 0.0 0.0  
 1.0 0.0
 1.0 1.0
 0.0 1.0
 0.5 0.5
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
366
367
The first two numbers are interpreted as the coordinates of vertex 0,
and so on.
Thomas Witkowski's avatar
Thomas Witkowski committed
368

Thomas Witkowski's avatar
Thomas Witkowski committed
369
370
Corresponding to these vertex indices now the four triangles are
given:
Thomas Witkowski's avatar
Thomas Witkowski committed
371
372
373
374
375
376
377
378
379

\begin{lstlisting}{}
element vertices:
0 1 4 
1 2 4  
2 3 4 
3 0 4 
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
380
381
Element 0 consists in the vertices 0, 1 and 4. The numbering is done
anticlockwise starting with the vertices of the longest edge.
Thomas Witkowski's avatar
Thomas Witkowski committed
382
383
384
385
386
387
388
389
390
391
392

It follows the definition of boundary conditions:

\begin{lstlisting}{}
element boundaries:
0 0 1 
0 0 1
0 0 1 
0 0 1 
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
393
394
395
396
397
398
399
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.
Thomas Witkowski's avatar
Thomas Witkowski committed
400

Thomas Witkowski's avatar
Thomas Witkowski committed
401
402
403
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.
Thomas Witkowski's avatar
Thomas Witkowski committed
404
405
406
407
408
409
410
411
412

\begin{lstlisting}{}
element neighbours:
1 3 -1 
2 0 -1 
3 1 -1 
0 2 -1
\end{lstlisting}

Thomas Witkowski's avatar
Thomas Witkowski committed
413
414
This block is optional. If it isn't given in the macro file, element
neighborships are computed automatically.
Thomas Witkowski's avatar
Thomas Witkowski committed
415
416
417
418
419
420
421
422
423
424
425
426
427

\subsection{Output}
\begin{figure}
\center
\begin{tabular}{cc}
\includegraphics[width=0.3\textwidth]{fig/ellipt_dat.jpg}&
\includegraphics[width=0.3\textwidth]{fig/ellipt_mesh.jpg}\\
(a)&(b)
\end{tabular}
\caption{(a): Solution after 9 iterations, (b): corresponding mesh}
\label{f:ellipt solution and mesh}
\end{figure}

Thomas Witkowski's avatar
Thomas Witkowski committed
428
429
430
431
432
433
434
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}.