Skip to content

Poisson's equation

Solving Poisson’s equation#

For a given function f(x,y), find a function u(x,y) satisfying :

\begin{eqnarray} -\Delta u(x,y) &= f(x,y) & \mbox{ for all }(x,y)\mbox{ in }\Omega\label{eqn:Poisson}\\ u(x,y) &= 0 & \mbox{ for all }(x,y)\mbox{ on }\p\Omega\label{eqn:Dirichlet} \end{eqnarray}

Here \p\Omega is the boundary of the bounded open set \Omega\subset\R^2 and \Delta u = \frac{\p^2 u}{\p x^2 } + \frac{\p^2 u}{\p y^2}.

We will compute u with f(x,y)=xy and \Omega the unit disk. The boundary C=\partial\Omega is defined as:

Note

In FreeFem++, the domain \Omega is assumed to be described by the left side of its boundary.

The following is the Freefem++ program which computes u:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
// Define mesh boundary
border C(t=0, 2*pi){x=cos(t); y=sin(t);}

// The triangulated domain Th is on the left side of its boundary
mesh Th = buildmesh(C(50));

// The finite element space defined over Th is called here Vh
fespace Vh(Th, P1);
Vh u, v;// Define u and v as piecewise-P1 continuous functions

// Define a function f
func f= x*y;

// Get the clock in second
real cpu=clock();

// Define the PDE
solve Poisson(u, v, solver=LU)
    = int2d(Th)(    // The bilinear part
          dx(u)*dx(v)
        + dy(u)*dy(v)
    )
    - int2d(Th)(    // The right hand side
          f*v
    )
    + on(C, u=0);   // The Dirichlet boundary condition

// Plot the result
plot(u);

// Display the total computational time
cout << "CPU time = " << (clock()-cpu) << endl;
As illustrated in figure 1.2, we can see the isovalue of u by using Freefem++ plot command (see line 29 above).

Figure 1.1 - Mesh Th by buildmesh(C(50)) Figure 1.2 - isovalue by plot(u)
mesh TH isovalue

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.

Note how close to the mathematics FreeFem++ language is. Lines 19 to 24 correspond to the mathematical variational equation: for all v which are in the finite element space V_h and zero on the boundary C.

Exercise:

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 \Gamma is described analytically by a parametric equation for x and for y. When \Gamma=\sum_{j=0}^J \Gamma_j then each curve \Gamma_j must be specified and crossings of \Gamma_j are not allowed except at end points .

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:

1
2
border Gamma1(t=0, pi){x=cos(t); y=sin(t); label=C};
border 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 \ref{Meshing Examples}).

On the line 5

The triangulation \mathcal{T}_h of \Omega is automatically generated by buildmesh(C(50)) using 50 points on C as in figure 1.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:

1
border C(t=2.*pi, 0){x=0.1+0.3*cos(t); y=0.5*sin(t);};

If by mistake one had written

1
border 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.

Info

Automatic mesh generation is based on the Delaunay-Voronoi algorithm. Refinement of the mesh are done by increasing the number of points on \Gamma, for example 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 \mathcal{T}_h (Th in FreeFem++) refers to the family \{T_k\}_{k=1,\cdots,n_t} of triangles shown in figure 1.1.

Traditionally h refers to the mesh size, n_t to the number of triangles in \mathcal{T}_h and n_v to the number of vertices, but it is seldom that we will have to use them explicitly.

If \Omega is not a polygonal domain, a "skin" remains between the exact domain \Omega and its approximation \Omega_h=\cup_{k=1}^{n_t}T_k. However, we notice that all corners of \Gamma_h = \p\Omega_h are on \Gamma.

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 V_h to be the space of continuous functions which are affine in x,y on each triangle of T_h.

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 \phi_k which are continuous piecewise affine and are equal to 1 on one vertex and 0 on all others. A typical hat function is shown on figure 2.2.

Figure 2.1: mesh Th Figure 2.2: Graph of \phi_1 (left) and \phi_6
mesh Th Typical hat functions

Info

The easiest way to define \phi_k is by making use of the barycentric coordinates \lambda_i(x,y),~i=1,2,3 of a point q=(x,y)\in T, defined by \sum_i\lambda_i=1,~~~\sum_i\lambda_i\vec q^i=\vec q where q^i,~i=1,2,3 are the 3 vertices of T. Then it is easy to see that the restriction of \phi_k on T is precisely \lambda_k.

Then

where M is the dimension of V_h, i.e. the number of vertices. The w_k are called the degree of freedom of w and M the number of degree of freedom.

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 u and v are approximated as above, namely

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 (\ref{eqn:Poisson}) and its Dirichlet boundary conditions (\ref{eqn:Dirichlet}).

This variational formulation is derived by multiplying (\ref{eqn:Poisson}) by v(x,y) and integrating the result over \Omega:

Then, by Green's formula, the problem is converted into finding u such that

\begin{eqnarray} \label{eqn:weakform} &&a(u,v) - \ell(f,v) = 0 \qquad \forall v \hbox{ satisfying $v=0$ on }\p\Omega.\\ &&\hbox{with }a(u,v)=\int_{\Omega}\nabla u\cdot \nabla v \,\d x\d y , \quad \ell(f,v)=\int_{\Omega}fv\, \d x\d y \label{eqn:bilinear} \end{eqnarray}

In FreeFem++ the Poisson problem can be declared only as in

1
Vh u,v; problem Poisson(u,v) = ...
and solved later as in
1
Poisson; //the problem is solved here
or declared and solved at the same time as in
1
Vh u,v; solve Poisson(u,v) = ...
and (\ref{eqn:weakform}) is written with dx(u) =\p u/\p x, dy(u) =\p u/\p y and

\displaystyle{\int_{\Omega}\nabla u\cdot \nabla v\, \d x\d y \longrightarrow} int2d(Th)( dx(u)*dx(v) + dy(u)*dy(v) )

\displaystyle{\int_{\Omega}fv\, \d x\d y \longrightarrow} int2d(Th)( f*v ) (Notice here, u is unused)

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 figure 1.2.

(see \ref{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

\begin{eqnarray} \label{eqn:Equation} \sum_{j=0}^{M-1} A_{ij}u_j - F_i=0 ,\quad i=0,\cdots,M-1;\qquad F_i=\int_{\Omega}f\phi_i\, \d x\d y \end{eqnarray}

which is found by using (\ref{defu}) and replacing v by \phi_i in (\ref{eqn:weakform}). The Dirichlet conditions are implemented by penalty, namely by setting A_{ii}=10^{30} and F_i=10^{30}*0 if i is a boundary degree of freedom.

Info

The number 10^{30} is called tgv (très grande valeur or very high value in english) and it is generally possible to change this value, see the item solve, tgv=

The matrix A=(A_{ij}) is called stiffness matrix. If the user wants to access A directly he can do so by using (see section \ref{matrix-varf} page \pageref{matrix-varf} for details).

1
2
3
4
5
6
7
8
varf a(u,v)
    = int2d(Th)(
          dx(u)*dx(v)
        + dy(u)*dy(v)
    )
    + on(C, u=0)
    ;
matrix A = a(Vh, Vh); //stiffness matrix

The vector F in (\ref{eqn:Equation}) can also be constructed manually

1
2
3
4
5
6
7
8
varf l(unused,v)
    = int2d(Th)(
          f*v
    )
    + on(C, unused=0)
    ;
Vh F;
F[] = l(0,Vh); //F[] is the vector associated to the function F

The problem can then be solved by

1
u[] = A^-1*F[]; //u[] is the vector associated to the function u

Info

Here u and F are finite element function, and u[] and F[] give the array of value associated (u[] \equiv (u_i)_{i=0,\dots,M-1} and F[] \equiv (F_i)_{i=0,\dots,M-1}).

So we have where \phi_i, i=0...,,M-1 are the basis functions of Vh like in equation (\ref{equation3}), and M = \mathtt{Vh.ndof} is the number of degree of freedom (i.e. the dimension of the space Vh).

The linear system (\ref{eqn:Equation}) is solved by UMFPACK unless another option is mentioned specifically as in

1
2
Vh u, v;
problem Poisson(u, v, solver=CG) = int2d(...
meaning that Poisson is declared only here and when it is called (by simply writing Poisson;) then (\ref{eqn:Equation}) will be solved by the Conjugate Gradient method.