FreeFEM Documentation on GitHub

stars - forks


Summary : Here we shall learn how to solve a Dirichlet and/or mixed Dirichlet Neumann problem for the Laplace operator with application to the equilibrium of a membrane under load. We shall also check the accuracy of the method and interface with other graphics packages

An elastic membrane \(\Omega\) is attached to a planar rigid support \(\Gamma\), and a force \(f(x) dx\) is exerted on each surface element \(\text{d}{x}=\text{d}{x}_1 \text{d}{x}_2\). The vertical membrane displacement, \(\varphi(x)\), is obtained by solving Laplace’s equation:

\[-\Delta \varphi =f ~\hbox{in}~ \Omega\]

As the membrane is fixed to its planar support, one has:

\[\varphi |_{\Gamma }=0\]

If the support wasn’t planar but had an elevation \(z(x_1,x_2)\) then the boundary conditions would be of non-homogeneous Dirichlet type.


If a part \(\Gamma_2\) of the membrane border \(\Gamma\) is not fixed to the support but is left hanging, then due to the membrane’s rigidity the angle with the normal vector \(\boldsymbol{n}\) is zero; thus the boundary conditions are:


where \(\Gamma_1=\Gamma-\Gamma_2\); recall that \(\frac{\partial\varphi}{\partial\boldsymbol{n}}=\nabla\varphi\cdot \boldsymbol{n}\) Let us recall also that the Laplace operator \(\Delta\) is defined by:

\[\Delta \varphi = {\partial ^{2}\varphi \over \partial x^{2}_{1} } + {\partial ^{2}\varphi \over \partial x_{2}^{2} }\]


Check references

With such “mixed boundary conditions” the problem has a unique solution (see Dautray-Lions (1988), Strang (1986) and Raviart-Thomas (1983)). The easiest proof is to notice that \(\varphi\) is the state of least energy, i.e.

\[E(\phi) =\min_{\varphi-z\in V} E(v) ,\quad \mbox{with} \quad E(v)=\int_\Omega(\frac12|\nabla v|^2-fv )\]

and where \(V\) is the subspace of the Sobolev space \(H^1(\Omega)\) of functions which have zero trace on \(\Gamma_1\). Recall that (\(x\in\mathbb{R}^d,~d=2\) here):

\[H^1(\Omega)=\{u\in L^2(\Omega)~:~\nabla u\in (L^2(\Omega))^d\}\]

Calculus of variation shows that the minimum must satisfy, what is known as the weak form of the PDE or its variational formulation (also known here as the theorem of virtual work)

\[\int_\Omega \nabla\varphi\cdot\nabla w = \int_\Omega f w\quad\forall w\in V\]

Next an integration by parts (Green’s formula) will show that this is equivalent to the PDE when second derivatives exist.


Unlike the previous version Freefem+ which had both weak and strong forms, FreeFEM implements only weak formulations. It is not possible to go further in using this software if you don’t know the weak form (i.e. variational formulation) of your problem: either you read a book, or ask help form a colleague or drop the matter. Now if you want to solve a system of PDE like \(A(u,v)=0,~ B(u,v)=0\) don’t close this manual, because in weak form it is

\[\int_\Omega(A(u,v)w_1+B(u,v)w_2)=0~~\forall w_1,w_2...\]


Let an ellipse have the length of the semimajor axis \(a=2\), and unitary the semiminor axis. Let the surface force be \(f=1\). Programming this case with FreeFEM gives:

 1// Parameters
 2real theta = 4.*pi/3.;
 3real a = 2.; //The length of the semimajor axis
 4real b = 1.; //The length of the semiminor axis
 5func z = x;
 7// Mesh
 8border Gamma1(t=0., theta){x=a*cos(t); y=b*sin(t);}
 9border Gamma2(t=theta, 2.*pi){x=a*cos(t); y=b*sin(t);}
10mesh Th = buildmesh(Gamma1(100) + Gamma2(50));
12// Fespace
13fespace Vh(Th, P2); //P2 conforming triangular FEM
14Vh phi, w, f=1;
16// Solve
17solve Laplace(phi, w)
18    = int2d(Th)(
19          dx(phi)*dx(w)
20        + dy(phi)*dy(w)
21    )
22    - int2d(Th)(
23          f*w
24    )
25    + on(Gamma1, phi=z)
26    ;
28// Plot
29plot(phi, wait=true, ps="membrane.eps"); //Plot phi
30plot(Th, wait=true, ps="membraneTh.eps"); //Plot Th
32// Save mesh

A triangulation is built by the keyword buildmesh. This keyword calls a triangulation subroutine based on the Delaunay test, which first triangulates with only the boundary points, then adds internal points by subdividing the edges. How fine the triangulation becomes is controlled by the size of the closest boundary edges.

The PDE is then discretized using the triangular second order finite element method on the triangulation; as was briefly indicated in the previous chapter, a linear system is derived from the discrete formulation whose size is the number of vertices plus the number of mid-edges in the triangulation.

The system is solved by a multi-frontal Gauss LU factorization implemented in the package UMFPACK.

The keyword plot will display both \(\T_h\) and \(\varphi\) (remove Th if \(\varphi\) only is desired) and the qualifier fill=true replaces the default option (colored level lines) by a full color display.

1plot(phi,wait=true,fill=true); //Plot phi with full color display

Results are on Fig. 5 and Fig. 6.


Fig. 5 Mesh of the ellipse


Fig. 6 Level lines of the membrane deformation


Next we would like to check the results !

One simple way is to adjust the parameters so as to know the solutions. For instance on the unit circle a=1, \(\varphi_e=\sin(x^2+y^2-1)\) solves the problem when:

\[z=0, f=-4(\cos(x^2+y^2-1)-(x^2+y^2)\sin(x^2+y^2-1))\]

except that on \(\Gamma_2\) \(\partial_n\varphi=2\) instead of zero. So we will consider a non-homogeneous Neumann condition and solve:

\[\int_\Omega\nabla\varphi\cdot\nabla w = \int_\Omega f w+\int_{\Gamma_2}2w\quad\forall w\in V\]

We will do that with two triangulations, compute the \(L^2\) error:

\[\epsilon = \int_\Omega|\varphi-\varphi_e|^2\]

and print the error in both cases as well as the log of their ratio an indication of the rate of convergence.

 1//  Parameters
 2verbosity = 0; //to remove all default output
 3real theta = 4.*pi/3.;
 4real a=1.; //the length of the semimajor axis
 5real b=1.; //the length of the semiminor axis
 6func f = -4*(cos(x^2+y^2-1) - (x^2+y^2)*sin(x^2+y^2-1));
 7func phiexact = sin(x^2 + y^2 - 1);
 9// Mesh
10border Gamma1(t=0., theta){x=a*cos(t); y=b*sin(t);}
11border Gamma2(t=theta, 2.*pi){x=a*cos(t); y=b*sin(t);}
13// Error loop
14real[int] L2error(2); //an array of two values
15for(int n = 0; n < 2; n++){
16    // Mesh
17    mesh Th = buildmesh(Gamma1(20*(n+1)) + Gamma2(10*(n+1)));
19    // Fespace
20    fespace Vh(Th, P2);
21    Vh phi, w;
23    // Solve
24    solve Laplace(phi, w)
25        = int2d(Th)(
26              dx(phi)*dx(w)
27            + dy(phi)*dy(w)
28        )
29        - int2d(Th)(
30              f*w
31        )
32        - int1d(Th, Gamma2)(
33              2*w
34        )
35        + on(Gamma1,phi=0)
36        ;
38    // Plot
39    plot(Th, phi, wait=true, ps="membrane.eps");
41    // Error
42    L2error[n] = sqrt(int2d(Th)((phi-phiexact)^2));
45// Display loop
46for(int n = 0; n < 2; n++)
47    cout << "L2error " << n << " = " << L2error[n] << endl;
49// Convergence rate
50cout << "convergence rate = "<< log(L2error[0]/L2error[1])/log(2.) << endl;

The output is:

1L2error 0 = 0.00462991
2L2error 1 = 0.00117128
3convergence rate = 1.9829
4times: compile 0.02s, execution 6.94s

We find a rate of 1.98 , which is not close enough to the 3 predicted by the theory.

The Geometry is always a polygon so we lose one order due to the geometry approximation in \(O(h^2)\).

Now if you are not satisfied with the .eps plot generated by FreeFEM and you want to use other graphic facilities, then you must store the solution in a file very much like in C++. It will be useless if you don’t save the triangulation as well, consequently you must do

2    ofstream ff("phi.txt");
3    ff << phi[];

For the triangulation the name is important: the extension determines the format.


Fig. 7 The 3D version drawn by gnuplot from a file generated by FreeFEM

Still that may not take you where you want. Here is an interface with gnuplot (see : web site link ) to produce the Fig. 7.

 1//to build a gnuplot data file
 3    ofstream ff("graph.txt");
 4    for (int i = 0; i < Th.nt; i++)
 5    {
 6        for (int j = 0; j < 3; j++)
 7            ff << Th[i][j].x << " "<< Th[i][j].y << " " << phi[][Vh(i,j)] << endl;
 9        ff << Th[i][0].x << " " << Th[i][0].y << " " << phi[][Vh(i,0)] << "\n\n\n"
10    }

We use the finite element numbering, where Wh(i,j) is the global index of \(j^{Th}\) degrees of freedom of triangle number \(i\).

Then open gnuplot and do:

1set palette rgbformulae 30,31,32
2splot "graph.txt" w l pal

This works with P2 and P1, but not with P1nc because the three first degrees of freedom of P1 or P2 are on vertices and not with P1nc.

Table of content