1
\chapter{Getting started with \Rheolef}
3
For obtaining and installing \Rheolef, see the installation instructions on the \Rheolef\ home page:
5
\url{http://www-ljk.imag.fr/membres/Pierre.Saramito/rheolef/}
7
All examples presented along the present book
8
are available in the \code{example/} directory
9
of the \Rheolef\ distribution.
10
This directory is given by the following unix command:
11
\cindex{directory of example files}%
12
\pindex{rheolef-config}%
13
\pindex{rheolef-config!--example-dir}%
15
rheolef-config --exampledir
17
This command returns you a path, something like \code{/usr/share/doc/rheolef-doc/examples/}
18
and you should make a copy of these files:
20
cp -a /usr/share/doc/rheolef-doc/examples/ .
23
Before to run examples, please check your \Rheolef\ installation with:
24
\pindex{rheolef-config!--check}%
26
rheolef-config --check
29
% --------------------------------------
30
\section{The model problem}
31
% --------------------------------------
32
\label{sec-dirichlet}%
35
Let us consider the classical Poisson problem
36
with homogeneous Dirichlet boundary conditions
37
in a domain bounded $\Omega \subset R^d$,
40
\cindex{operator!Laplace}%
41
{\em (P): find $u$, defined in $\Omega$, such that:}
43
-\Delta u &=& 1 \ {\rm in}\ \Omega
44
\label{eq-dirichlet-1}
46
u &=& 0 \ {\rm on}\ \partial \Omega
47
\label{eq-dirichlet-2}
49
where $\Delta$ denotes the Laplace operator.
50
The variational formulation of this problem expresses
51
(see appendix~\ref{sec-green-scalar} for details):
53
{\em (VF): find $u\in H^1_0(\Omega)$ such that:}
56
\ \forall v \in H^1_0(\Omega)
57
\label{eq-dirichlet-fv}
59
where the bilinear form $a(.,.)$ and the linear form $l(.)$ are defined by
61
a(u,v) &=& \int_\Omega \nabla u . \nabla v \, {\rm d}x,
62
\ \ \forall u,v \in H^1_0(\Omega)
64
l(v) &=& \int_\Omega v \, {\rm d}x,
65
\ \ \forall v \in L^2(\Omega)
67
The bilinear form $a(.,.)$ defines a scalar product
69
in $H^1_0(\Omega)$ and is related to the {\em energy} form.
70
This form is associated to the $-\Delta$ operator.
71
% The $m(.,.)$ is here the classical scalar product on $L^2(\Omega)$,
72
% and is related to the {\em mass} form.
74
% --------------------------------------
75
\section{Approximation}
76
% --------------------------------------
77
\cindex{approximation} %
78
Let us introduce a mesh ${\cal T}_h$ of $\Omega$
79
and the finite dimensional space $X_h$ of continuous
80
piecewise polynomial functions.
82
X_h = \{ v \in H^1(\Omega); \
84
\forall K \in {\cal T}_h \}
87
Let $V_h=X_h \cap H^1_0(\Omega)$ be the functions of $X_h$ that
88
vanishes on the boundary of $\Omega$.
89
The approximate problem expresses:
91
{\em $(VF)_h$: find $u_h\in V_h$ such that:}
96
By developing $u_h$ on a basis of $V_h$, this
97
problem reduces to a linear system.
98
The following C++ code
99
implement this problem in the \Rheolef\ environment.
102
% --------------------------------------
103
\section{File \file{dirichlet.cc}}
104
% --------------------------------------
106
\exindex{dirichlet.cc}
107
\myexample{dirichlet.cc}
109
% --------------------------------------
111
% --------------------------------------
117
This code applies for both one, two or three
118
dimensional meshes and
119
for both piecewise linear or quadratic
120
finite element approximations.
121
Four major classes are involved, namely:
122
{\color{rheoleftypcolor} \code{geo}},
123
{\color{rheoleftypcolor} \code{space}},
124
{\color{rheoleftypcolor} \code{form}} and
125
{\color{rheoleftypcolor} \code{field}}.
127
Let us now comment the code, line by line.
128
\begin{lstlisting}[numbers=none,frame=none]
131
The first line includes the \Rheolef\ header file \file{rheolef.h}.
132
\cindex{namespace!std}%
133
\cindex{namespace!rheolef}%
134
\begin{lstlisting}[numbers=none,frame=none]
135
using namespace rheolef;
138
By default, in order to avoid possible name conflicts when using another library,
139
all classe and function names are prefixed by \code{rheolef::}, as in~\code{rheolef::space}.
140
This feature is called the name space.
141
Here, since there is no possible conflict, and in order to simplify the syntax, we drop all the
142
\code{rheolef::} prefixes, and do the same with the standard \code{c++} library classes and
143
variables, that are also prefixed by \code{std::}.
144
\cindex{argc, argv, command line arguments}%
145
\begin{lstlisting}[numbers=none,frame=none]
146
int main(int argc, char**argv) {
148
The entry function of the program is always called \code{main} and accepts arguments
149
from the unix command line: \code{argc} is the counter of command line arguments
150
and \code{argv} is the table of values. The character string \code{argv[0]}
151
is the program name and \code{argv[i]}, for $i=1$ to \code{argc-1}, are the
152
additional command line arguments.
153
\begin{lstlisting}[numbers=none,frame=none]
154
environment rheolef (argc, argv);
156
\toindex{MPI, message passing interface}%
157
\cindex{distributed computation}%
158
\cindex{parallel computation}%
159
\toindex{library!boost}%
160
These two command line parameters are immediatly furnished to the
161
distributed environment initializer of the \code{boost::mpi} library, that is
162
a \code{c++} library based on the usual message passing interface ({\sc mpi}) library.
163
Notice that this initialization is requiered, even when you run with only one processor.
164
\begin{lstlisting}[numbers=none,frame=none]
167
This command get the first unix command-line argument {\tt argv[1]}
168
as a mesh file name and store the corresponding mesh in the variable {\tt omega}.
169
\begin{lstlisting}[numbers=none,frame=none]
170
space Xh (omega, argv[2]);
172
Build the finite element space {\tt Xh} contains all the piecewise polynomial
173
continuous functions. The polynomial type
174
is the second command-line arguments
178
\apindex{high-order}%
179
{\tt argv[2]}, and could be either \code{P1}, \code{P2} or any \code{P}$k$, where $k\geq 1$.
181
\begin{lstlisting}[numbers=none,frame=none]
182
Xh.block ("boundary");
184
The homogeneous Dirichlet conditions are declared on the boundary.
185
\begin{lstlisting}[numbers=none,frame=none]
186
form a (Xh, Xh, "grad_grad");
189
\foindex{grad\_grad}%
190
The form $a(.,.)$ is the energy form.
191
\begin{lstlisting}[numbers=none,frame=none]
192
field lh = riesz (Xh, 1);
194
\cindex{Riesz representer}
195
Here $lh(.)$ is the Riesz representer of the constant right-hand side $f=1$ of the problem.
196
\begin{lstlisting}[numbers=none,frame=none]
199
The field {\tt uh} contains the the degrees of freedom.
200
\begin{lstlisting}[numbers=none,frame=none]
203
Some degrees of freedom are prescribed as zero on the boundary.
205
\cindex{Lagrange!node}%
206
Let $(\varphi_i)_{0\leq i< {\rm dim}(X_h)}$ be the basis
207
of $X_h$ associated to the Lagrange nodes, e.g. the vertices
208
of the mesh for the $P_1$ approximation and the vertices and the
209
middle of the edges for the $P_2$ approximation.
210
The approximate solution $u_h$ expresses as a linear combination
211
of the continuous piecewise polynomial functions $(\varphi_i)$:
213
u_h = \sum_{i} u_i \varphi_i
215
Thus, the field $u_h$ is completely represented by its coefficients $(u_i)$.
216
The coefficients $(u_i)$ of this combination are grouped into
217
to sets: some have zero values, from the boundary condition
218
and are related to {\em blocked} coefficients,
219
and some others are {\em unknown}.
220
Blocked coefficients are stored into the \code{uh.b} array
221
while unknown one are stored into \code{uh.u}.
222
Thus, the restriction of the bilinear form $a(.,.)$ to $X_h\times X_h$
223
can be conveniently represented by a block-matrix structure:
224
\cindex{matrix!block structure}%
227
\left ( \begin{array}{cc}
228
{\tt vh.u} & {\tt vh.b}
230
\left ( \begin{array}{cc}
231
{\tt a.uu} & {\tt a.ub} \\
232
{\tt a.bu} & {\tt a.bb}
234
\left ( \begin{array}{c}
239
This representation also applies for the
243
\left ( \begin{array}{cc}
244
{\tt vh.u} & {\tt vh.b}
246
\left ( \begin{array}{c}
251
Thus, the problem $(VF)_h$ writes now:
253
\left ( \begin{array}{cc}
254
{\tt vh.u} & {\tt vh.b}
256
\left ( \begin{array}{cc}
257
{\tt a.uu} & {\tt a.ub} \\
258
{\tt a.bu} & {\tt a.bb}
260
\left ( \begin{array}{c}
265
\left ( \begin{array}{cc}
266
{\tt vh.u} & {\tt vh.b}
268
\left ( \begin{array}{c}
273
for any {\tt vh.u} and where {\tt vh.b} = 0.
274
After expansion, the problem reduces to
275
{\em find} {\tt uh.u} {\em such that}:
277
{\tt a.uu*uh.u} = {\tt l.u - a.ub*uh.b}
279
\cindex{matrix!factorization!Choleski}%
282
The resolution of this linear system for the {\tt a.uu} matrix is
283
then performed. A preliminary step build the $LDL^T$ factorization:
284
\begin{lstlisting}[numbers=none,frame=none]
287
Then, the second step solves the {\em unknown part}:
288
\begin{lstlisting}[numbers=none,frame=none]
289
uh.set_u() = sa.solve (lh.u() - a.ub()*uh.b());
291
\cindex{algorithm!conjugate gradient}%
292
\cindex{preconditioner!Choleski incomplete factorization}%
293
When $d>3$, a faster iterative strategy is automatically preferred
294
by the {\tt solver} class for solving
295
the linear system: in that case, the preliminary step build
296
an incomplete Choleski factorization preconditioner,
297
while the second step runs an iterative method:
298
the preconditioned conjugate gradient algorithm.
299
Finally, the field is printed to standard output:
300
\begin{lstlisting}[numbers=none,frame=none]
303
The {\color{rheolefvarcolor}\tt dout} stream is a specific variable defined in the \Rheolef\ library:
304
it is a distributed and parallel extension
305
of the usual {\tt cout} stream in {\tt C++}
306
% --------------------------------------
307
\section{How to compile the code}
308
% --------------------------------------
309
\cindex{compilation}%
313
First, create a \code{Makefile} as follow:
314
\exindex{Makefile.demo}%
316
\verbatiminput{Makefile.demo}
322
Now, your program, linked with \Rheolef,
323
is ready to run on a mesh.
325
% --------------------------------------
326
\section{How to run the program}
327
% --------------------------------------
331
\includegraphics[height=6.5cm]{dirichlet-2d-P1-gnuplot.png} &
332
\includegraphics[height=6.5cm]{dirichlet-2d-P2-gnuplot.png}
335
\caption{Solution of the model problem for $d=2$:
336
(left) $P_1$ element;
337
(right) $P_2$ element.}
338
\label{fig-dirichlet-2d}
340
\cindex{geometry!square}%
342
\pindex{mkgeo\_grid}%
343
\pindex{mkgeo\_grid!-t}%
344
\fiindex{\file{.geo} mesh}%
345
\toindex{visualization!mesh}%
348
mkgeo_grid -t 10 > square.geo
351
The first command generates a simple 10x10 bidimensional mesh
353
and stores it in the file \code{square.geo}.
354
The second command shows the mesh.
356
It uses \code{gnuplot} visualization program by default.
359
\fiindex{\file{.field} field}%
360
The next command performs the computation:
362
./dirichlet square.geo P1 > square.field
365
% --------------------------------------
366
\section{Distributed and parallel runs}
367
% --------------------------------------
369
\cindex{distributed computation}%
370
\cindex{parallel computation}%
371
Alternatively, a computation in a distributed and parallel environment writes:
373
mpirun -np 4 ./dirichlet square.geo P1 > square.field
380
\includegraphics[height=6.5cm]{dirichlet-2d-bw-gnuplot.pdf} &
381
%\includegraphics[height=6.5cm]{dirichlet-2d-elevation-gnuplot.png}
382
\includegraphics[height=6.5cm]{dirichlet-2d-elevation.png}
385
\caption{Alternative representations
386
of the solution of the model problem ($d=2$ and the
388
(left) in black-and-white;
389
(right) in elevation and stereoscopic anaglyph mode.}
390
\label{fig-dirichlet-2d-alt}
392
% ----------------------------------------------------------------
393
\section{Stereo visualization}
394
% ----------------------------------------------------------------
395
Also explore some graphic rendering modes
396
(see Fig.~\ref{fig-dirichlet-2d-alt}):
398
\pindex{field!-gray}%
399
\pindex{field!-stereo}%
400
\pindex{field!-nofill}%
401
\pindex{field!-elevation}%
402
\pindex{field!-mayavi}%
404
\cindex{visualization!elevation view}%
406
field square.field -bw
407
field square.field -gray
408
field square.field -mayavi
409
field square.field -elevation -nofill -stereo
411
\cindex{visualization!stereoscopic anaglyph}%
412
The last command shows the solution in elevation
413
and in stereoscopic anaglyph
414
mode (see Fig.~\ref{fig-dirichlet-3d}, left).
415
The anaglyph mode requires red-cyan glasses:
416
red for the left eye and cyan for the right one,
417
as shown on Fig.~\ref{fig-anaglyph-glasses}.
420
\includegraphics{red-cyan-glasses.png}
422
\caption{Red-cyan anaglyph glasses for the stereoscopic visualization.}
423
\label{fig-anaglyph-glasses}
425
See~\code{http://en.wikipedia.org/wiki/Anaglyph\_image}
427
\code{http://www.alpes-stereo.com/lunettes.html}
428
for how to find anaglyph red-cyan glasses.
429
Please, consults the corresponding unix manual page for more
430
on \code{field}, \code{geo} and \code{mkgeo\_grid}:
437
% ----------------------------------------------------------------
438
\section{High-order finite element methods}
439
% ----------------------------------------------------------------
442
\apindex{high-order}%
443
Turning to the \code{P2} or \code{P3} approximations simply writes:
445
./dirichlet square.geo P2 > square-P2.field
446
field square-P2.field
448
Fig.~\ref{fig-dirichlet-2d}.right shows the result.
449
You can replace the \code{P2} command-line argument by any \code{P}$k$, where $k\geq 1$.
451
\cindex{geometry!line}%
452
Now, let us consider a mono-dimensional problem
455
\pindex{mkgeo\_grid!-e}%
457
mkgeo_grid -e 10 > line.geo
459
./dirichlet line.geo P1 | field -
461
The first command generates a subdivision containing ten edge elements.
463
The last two lines show the mesh and the solution
464
via \code{gnuplot} visualization, respectively.
466
Conversely, the \code{P2} case writes:
468
./dirichlet line.geo P2 | field -
470
% ----------------------------------------------------------------
471
\section{Tridimensional computations}
472
% ----------------------------------------------------------------
473
Let us consider a three-dimensional problem $\Omega=]0,1[^3$.
474
First, let us generate a mesh:
475
\cindex{geometry!cube}%
476
\pindex{geo!-stereo}%
479
\pindex{mkgeo\_grid!-T}%
481
mkgeo_grid -T 10 > cube.geo
483
geo cube.geo -stereo -full
484
geo cube.geo -stereo -cut
486
The previous commands draw the mesh with all internal edges (\code{-full}),
487
stereoscopic anaglyph (\code{-stereo})and then with a cut (\code{-cut}
488
inside the internal structure:
489
a simple click on the central arrow draws the cut plane normal vector
490
or its origin, while the red square allows a translation.
494
\includegraphics[scale=0.30]{cube-3d-stereo-fig.png} &
495
\includegraphics[scale=0.35]{dirichlet-3d-mayavi-mix-fig.png}
497
\caption{Solution of the model problem for $d=3$ and the $P_1$ element~:
499
(right) isovalue, cut planes and stereo anaglyph renderings.}
500
\label{fig-dirichlet-3d}
503
Then, we perform the computation and the visualization:
505
./dirichlet cube.geo P1 > cube.field
509
The visualization presents an isosurface.
510
Also here, you can interact with the cutting plane.
511
Click on \code{IsoSurface} in the left menu and
512
change the value of the isosurface.
513
Finally exit from the visualization and explore the
514
stereoscopic anaglyph mode
515
(see Fig.~\ref{fig-dirichlet-3d}, right):
516
\pindex{field!-stereo}%
518
field cube.field -stereo
520
It is also possible to add a second \code{IsoSurface}
521
or \code{ScalarCutPlane} module to this scene
522
by using the \code{Visualize} menu.
524
%TODO 3d volume rendering ?
525
% Exit from the visualization and enter the command:
527
% field -volume cube.field
529
% Click on the \fbox{configure} button of the "Modules" menu at the bottom of the control
530
% panel. Then, select "Use TextureMapper" and inscreases the luminosity of the scene
531
% by editing the color table
532
% You can drag the mouse with different buttons to change the colors.
533
% The following are the mouse buttons and key combinations that can be used to drag
534
% and edit the color table and the luminosity:
536
% \item luminosity~: \fbox{Shift-Button-1}
537
% \item red~: \fbox{Button-1}
538
% \item green~: \fbox{Button-1}
539
% \item blue~: \fbox{Button-3} or \fbox{Ctrl-Button-1}
541
% \fiindex{\file{.vtk} mesh and field data}
542
After this exploration of the 3D visualisation capacities of our environment,
543
let us go back to the Dirichlet problem and
544
perform the \code{P2} approximation:
546
./dirichlet cube.geo P2 | field -
548
% ----------------------------------------------------------------
549
\section{Quandrangles, prisms and hexahedra}
550
% ----------------------------------------------------------------
551
Quadrangles and hexahedra are also supported in meshes:
552
\pindex{mkgeo\_grid!-q}%
553
\pindex{mkgeo\_grid!-H}%
555
mkgeo_grid -q 10 > square.geo
557
mkgeo_grid -H 10 > cube.geo
560
Notices also that the one-dimensional exact solution writes:
562
u(x) = \Frac{x(1-x)}{2}
564
while the two-and three dimensional ones
565
support a Fourier expansion (see e.g. \cite{saramito-roquet-00}, annex).
added
removed
doc/pusrman/dirichlet.tex
