FreeFem++ Tutorial – Part 1

First of all, FreeFem is a numerical computing software which allows a fast and automatized treatment of a variety of problems related to partial differential equations. Its name, FreeFem, speaks for itself: it is free and it uses the finite element method. Here are a few reasons for which you may choose to use FreeFem for a certain task:

  1. It allows the user to easily define 2D (and 3D) geometries and it does all the work regarding the construction of meshes on these domains.
  2. The problems you want to solve can be easily written in the program once we know their weak forms.
  3. Once we have variables defined on meshes or solutions to some PDE, we can easily compute all sorts of quantities like integral energies, etc.

Before showing a first example, you need to install FreeFem. If you are not familiar with command line work or you just want to get to work, like me, you can install the visual version of FreeFem which is available here. Of course, you can find example programs in the FreeFem manual or by making a brief search on the internet.

I’ll present some basic stuff, which will allow us in the end to solve the Laplace equation in a circular domain. Once we have the structure of the program, it is possible to change the shape of the domain in no time.

Define the geometry. It is possible to define the geometry of the domain in multiple ways. The one I’ll use below is the parametric method. Let’s say we want to work on the unit circle. Then we simply write

border C(t=0,2*pi){x=cos(t); y=sin(t);label=1}

Note that on each border we can put a label. This may help us later on when we define the problem or when we want to compute an integral over this specific part of the boundary. It is important to note that using labels we can impose different boundary conditions on several parts of the domain.

Construct the mesh. The mesh can be constructed using the simple command

mesh Th = buildmesh (C(450));

You can identify the command buildmesh and the parameter 450 dictates the number of meshuniform discretization points taken on the curve C. Once the mesh is built you can visualize it using the command plot(Th). Note that in FreeFem, like in C, every variable which appears for the first time must be specified by its type. Also, take care not to define a variable’s type a second time, because you’ll get an error.

Construct and solve the problem. Now we have the geometry and the mesh, so we can write the problem. First we define a finite element space on the constructed triangulation by using the command

fespace Vh(Th,P1);

Note that P1 simply means we use P1 finite elements. Once we have a finite element space we can define variables in this space.

Vh uh,vh;

Since we want to be able to impose a general boundary condition we define a function by using the command

func fu = exp(x^2);

Note that whenever the variables x and y are used they represent the usual coordinates in the plane, so you cannot use them for something else. All the basic functions like exp, sin, cos are already implemented in FreeFem so you can use them at will.

We are ready to state the problem in its variational form. Since we want to solve -\Delta u = 0 with the boundary condition u = f, we consider the weak form of the problem \displaystyle \int \nabla u\cdot \nabla v = 0 and we introduce this weak form in FreeFem using the commands

problem laplace(uh,vh) = int2d(Th)( dx(uh)*dx(vh) + dy(uh)*dy(vh) )

We note the use of the keyword problem to denote a weak form of an EDP. There are several other things present in the above command. Note that to integrate over a two dimensional mesh we simply use the int2d command with the first argument the mesh and the second argument the function we wish to integrate on that mesh. Note that we used variables defined on the finite element space on the corresponding mesh. Here are some simple tricks which allow to compute numerically the area and the perimeter of the triangulation. The idea is to use a dummy variable which is equal to 1 everywhere and then integrate it on the domain to find the area and on the boundary to find the perimeter.

real area;
Vh ff = 1;
area = int2d(Th)(ff);
real perim;
perim = int1d(Th,1)(ff);

If you wish to display a variable use a command like cout << perim << endl; Note that if we have more connected components, we can specify over which boundary component we wish to integrate by just putting the corresponding label after the mesh variable Th.

The on operator allows us to put a boundary condition on parts of the boundary specified by the corresponding label. Here, the command on(1,uh = fu) simply means that we impose a condition of the type u = f on the boundary. In order to solve the stated problem simply type its name:


Visualize the results. Now the variable uh contains the numerical solution of our problem. In order to visualize it just type


This will produce a 2D plot with colors corresponding to the values of the function uh. The soloption value=true adds a colorbar which shows the association between the values of the function and the respective colors.

The whole code can be written in a single .edp file, the type associated to FreeFem problems. In order to solve different problems, it suffices to change the problem formulation or the boundary condition. I hope I’ve convinced you that FreeFem is not that hard and that we can solve PDEs in just a few lines of code.


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s

%d bloggers like this: