Getting started
For a given function
Here
We will compute
Note
In FreeFEM, the domain
The following is the FreeFEM program which computes
1// Define mesh boundary
2border C(t=0, 2*pi){x=cos(t); y=sin(t);}
3
4// The triangulated domain Th is on the left side of its boundary
5mesh Th = buildmesh(C(50));
6
7// The finite element space defined over Th is called here Vh
8fespace Vh(Th, P1);
9Vh u, v;// Define u and v as piecewise-P1 continuous functions
10
11// Define a function f
12func f= x*y;
13
14// Get the clock in second
15real cpu=clock();
16
17// Define the PDE
18solve Poisson(u, v, solver=LU)
19 = int2d(Th)( // The bilinear part
20 dx(u)*dx(v)
21 + dy(u)*dy(v)
22 )
23 - int2d(Th)( // The right hand side
24 f*v
25 )
26 + on(C, u=0); // The Dirichlet boundary condition
27
28// Plot the result
29plot(u);
30
31// Display the total computational time
32cout << "CPU time = " << (clock()-cpu) << endl;
As illustrated in Fig. 2, we can see the isovalue of plot
command (see line 29 above).
Note
The qualifier solver=LU
(line 18) is not required and by default a multi-frontal LU
is used.
The lines containing clock
are equally not required.
Tip
Note how close to the mathematics FreeFEM language is.
Lines 19 to 24 correspond to the mathematical variational equation:
for all
Tip
Change P1
into P2
and run the program.
This first example shows how FreeFEM executes with no effort all the usual steps required by the finite element method (FEM). Let’s go through them one by one.
On the line 2:
The boundary
The keyword label
can be added to define a group of boundaries for later use (boundary conditions for instance).
Hence the circle could also have been described as two half circle with the same label:
1border Gamma1(t=0, pi){x=cos(t); y=sin(t); label=C};
2border Gamma2(t=pi, 2.*pi){x=cos(t); y=sin(t); label=C};
Boundaries can be referred to either by name (Gamma1
for example) or by label (C
here) or even by its internal number here 1 for the first half circle and 2 for the second (more examples are in Meshing Examples).
On the line 5
The triangulation buildmesh(C(50))
using 50 points on C
as in Fig. 1.
The domain is assumed to be on the left side of the boundary which is implicitly oriented by the parametrization. So an elliptic hole can be added by typing:
1border C(t=2.*pi, 0){x=0.1+0.3*cos(t); y=0.5*sin(t);};
If by mistake one had written:
1border C(t=0, 2.*pi){x=0.1+0.3*cos(t); y=0.5*sin(t);};
then the inside of the ellipse would be triangulated as well as the outside.
Note
Automatic mesh generation is based on the Delaunay-Voronoi algorithm.
Refinement of the mesh are done by increasing the number of points on buildmesh(C(100))
, because inner vertices are determined by the density of points on the boundary.
Mesh adaptation can be performed also against a given function f by calling adaptmesh(Th,f)
.
Now the name Th
in FreeFEM) refers to the family
Traditionally
If
On line 8:
A finite element space is, usually, a space of polynomial functions on elements, triangles here only, with certain matching properties at edges, vertices etc. Here fespace Vh(Th, P1)
defines
As it is a linear vector space of finite dimension, basis can be found.
The canonical basis is made of functions, called the hat function
Note
The easiest way to define
Then:
where
It is said also that the nodes of this finite element method are the vertices.
Setting the problem
On line 9, Vh u, v
declares that
On the line 12, the right hand side f
is defined analytically using the keyword func
.
Line 18 to 26 define the bilinear form of equation (1) and its Dirichlet boundary conditions.
This variational formulation is derived by multiplying (1) by
Then, by Green’s formula, the problem is converted into finding
with:
In FreeFEM the Poisson problem can be declared only as in:
1Vh u,v; problem Poisson(u,v) = ...
and solved later as in:
1Poisson; //the problem is solved here
or declared and solved at the same time as in:
1Vh u,v; solve Poisson(u,v) = ...
and (4) is written with dx(u)
dy(u)
int2d(Th)( dx(u)*dx(v) + dy(u)*dy(v) )
int2d(Th)( f*v )
(Notice here,
Warning
In FreeFEM bilinear terms and linear terms should not be under the same integral indeed to construct the linear systems FreeFEM finds out which integral contributes to the bilinear form by checking if both terms, the unknown (here u
) and test functions (here v
) are present.
Solution and visualization
On line 15, the current time in seconds is stored into the real-valued variable cpu
.
Line 18, the problem is solved.
Line 29, the visualization is done as illustrated in Fig. 2.
(see Plot for zoom, postscript and other commands).
Line 32, the computing time (not counting graphics) is written on the console. Notice the C++-like syntax; the user needs not study C++ for using FreeFEM, but it helps to guess what is allowed in the language.
Access to matrices and vectors
Internally FreeFEM will solve a linear system of the type
which is found by using (3) and replacing
Note
The number tgv
(très grande valeur or very high value in english) and it is generally possible to change this value, see the item :freefem
The matrix
1varf a(u,v)
2 = int2d(Th)(
3 dx(u)*dx(v)
4 + dy(u)*dy(v)
5 )
6 + on(C, u=0)
7 ;
8matrix A = a(Vh, Vh); //stiffness matrix
The vector
1varf l(unused,v)
2 = int2d(Th)(
3 f*v
4 )
5 + on(C, unused=0)
6 ;
7Vh F;
8F[] = l(0,Vh); //F[] is the vector associated to the function F
The problem can then be solved by:
1u[] = A^-1*F[]; //u[] is the vector associated to the function u
Note
Here u
and F
are finite element function, and u[]
and F[]
give the array of value associated (u[]
F[]
So we have:
where
The linear system (5) is solved by UMFPACK
unless another option is mentioned specifically as in:
1Vh u, v;
2problem Poisson(u, v, solver=CG) = int2d(...
meaning that Poisson
is declared only here and when it is called (by simply writing Poisson;
) then (5) will be solved by the Conjugate Gradient method.