FreeFEM Documentation on GitHub

stars - forks


Standard types


Integer value (equivalent to long in C++).

1int i = 0;


Boolean value.

1bool b = true;


The result of a comparison is a boolean

bool b = (1 < 2);


Real value (equivalent to double in C++).

1real r = 0.;


Complex value (equivalent to two double or complex<double> in C++).

1complex c = 0. + 1i;

The imaginary number \(i\) is defined as 1i



1complex a = 1i, b = 2 + 3i;
2cout << "a + b = " << a + b << endl;
3cout << "a - b = " << a - b << endl;
4cout << "a*b = " << a*b << endl;
5cout << "a/b = " << a/b << endl;

The output of this script is:

a + b = (2,4)
a - b = (-2,-2)
a*b = (-3,2)
a/b = (0.230769,0.153846)


See Complex example for a detailed example.


String value.

1string s = "this is a string";


string value is enclosed within double quotes.

Other types can be concatenate to a string, like:

1int i = 1;
2real r = 1.;
3string s = "the int i = " + i +", the real r = " + r + ", the complex z = " + (1. + 1i);

To append a string in a string at position 4:

1s(4:3) = "+++";

To copy a substring in an other string:

1string s2 = s1(5:10);

See String Example for a complete example.

Mesh design


Border type.

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

Define the 2D geometrical border in parametric coordinates.



A label can be defined with the border:

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


Inner variable

An inner variable can be defined inside a border:

1border b(t=0., 1.){real tt=2.*pi*t; x=cos(tt); y=sin(tt);}


From vector

A border can be defined from two vectors using P.x and P.y:

1border b(t=0, vectorX.n-1){P.x=vectorX[t]; P.y=vectorY[t];}


2D Mesh type (see Mesh Generation).

1mesh Th;


3D mesh type (see Mesh Generation).

1mesh3 Th;

Finite element space design


Finite element space type (see Finite Element).

1fespace Uh(Th, P1);
2fespace UPh(Th, [P2, P2, P1]);

A finite element space is based on a mesh (Th) with an element definition, scalar (P1) or vector ([P2, P2, P1]).

Available finite element space:


  • P0 / P03d

  • P0Edge

  • P1 / P13d

  • P1dc

  • P1b / P1b3d

  • P1bl / P1bl3d

  • P1nc

  • P2 / P23d

  • P2b

  • P2dc

  • P2h

  • RT0 / RT03d

  • RT0Ortho

  • Edge03d

Using Element_P3:

  • P3

Using Element_P3dc:

  • P3dc

Using Element_P4:

  • P4

Using Element_P4dc:

  • P4dc

Using Element_PkEdge:

  • P1Edge

  • P2Edge

  • P3Edge

  • P4Edge

  • P5Edge

Using Morlay:

  • P2Morley

Using HCT:

  • HCT

Using BernardiRaugel:

  • P2BR

Using Element_Mixte:

  • RT1

  • RT1Ortho

  • RT2

  • RT2Ortho

  • BDM1

  • BDM1Ortho

Using Element_Mixte3d:

  • Edge13d

  • Edge23d

Using Element_QF:

  • FEQF

A finite element function is defined as follow:

1fespace Uh(Th, P1);
2Uh u;
4fespace UPh(Th, [P2, P2, P1]);
5UPh [Ux, Uy, p];

Macro design


Macro type.

1macro vU() [Ux, Uy] //
2macro grad(u) [dx(u), dy(u)] //

Macro ends with //.


Macro concatenation

You can use the C concatenation operator ## inside a macro using #.

If Ux and Uy are defined as finite element function, you can define:

1macro Grad(U) [grad(U#x), grad(U#y)] // End of macro

See Macro example

NewMacro / EndMacro


In developement - Not tested

Set and end a macro

1NewMacro grad(u) [dx(u), dy(u)] EndMacro


Check if a macro exists and check its value.

1IFMACRO(AA) //check if macro AA exists
5IFMACRO(AA, tt) //check if amcro exists and is equall to tt


Functions design


Function type.

Function without parameters (\(x\), \(y\) and \(z\) are implicitly considered):

1func f = x^2 + y^2;


Function’s type is defined by the expression’s type.

Function with parameters:

1func real f (real var){
2    return x^2 + y^2 + var^2;

Elementary functions

Class of basic functions (polynomials, exponential, logarithmic, trigonometric, circular) and the functions obtained from those by the four arithmetic operations

\[f(x) + g(x),\, f(x) - g(x),\, f(x)g(x),\, f(x)/g(x)\]

and by composition \(f(g(x))\), each applied a finite number of times.

In FreeFEM, all elementary functions can thus be created. The derivative of an elementary function is also an elementary function; however, the indefinite integral of an elementary function cannot always be expressed in terms of elementary functions.

See Elementary function example for a complete example.

Random functions

FreeFEM includes the Mersenne Twister random number generator. It is a very fast and accurate random number generator of period \(2^{219937}-1\).

See randint32(), randint31(), randreal1(), randreal2(), randreal3(), randres53(), randinit(seed).

In addition, the ffrandom plugin interface random, srandom and srandomdev functions of the Unix libc library. The range is \(0 -- 2^{31}-1\).


If srandomdev is not defined, a seed based on the current time is used.

gsl plugin equally allows usage of all random functions of the gsllib, see gsl external library.


Finite element functions are also constructed like elementary functions by an arithmetic formula involving elementary functions.

The difference is that they are evaluated at declaration time and FreeFEM stores the array of its values at the places associated with he degree of freedom of the finite element type. By opposition, elementary functions are evaluated only when needed. Hence FE-functions are not defined only by their formula but also by the mesh and the finite element which enter in their definitions.

If the value of a FE-function is requested at a point which is not a degree of freedom, an interpolation is used, leading to an interpolation error, while by contrast, an elementary function can be evaluated at any point exactly.

1func f = x^2*(1+y)^3 + y^2;
2mesh Th = square(20, 20, [-2+4*x, -2+4*y]); // ]-2, 2[^2
3fespace Vh(Th, P1);
4Vh fh=f; //fh is the projection of f to Vh (real value)
5func zf = (x^2*(1+y)^3 + y^2)*exp(x + 1i*y);
6Vh<complex> zh = zf; //zh is the projection of zf to complex value Vh space

The construction of fh = f is explained in Finite Element.


The plot command only works for real or complex FE-functions, not for elementary functions.

Problem design


Problem type.

1problem Laplacian (u, uh) = ...

FreeFEM needs the variational form in the problem definition.

In order to solve the problem, just call:




A solver can be specified in the problem definition:

1problem Laplacian(u, uh, solver=CG) = ...

The default solver is sparsesolver or LU if any direct sparse solver is available.

Solvers are listed in the Global variables section.


Stop test

A criterion \(\varepsilon\) can be defined for iterative methods, like CG for example:

1problem Laplacian(u, uh, solver=CG, eps=1.e-6) = ...

If \(\varepsilon>0\), the stop test is:

\[||Ax-b|| < \varepsilon\]

Else, the stop test is:

\[||Ax-b|| < \frac{|\varepsilon|}{||Ax_0-b||}\]



The keyword init controls the reconstruction of the internal problem matrix.

If init is set to false or 0, the matrix is reconstructed et each problem calls (or after a mesh modification), else the previously constructed matrix is used.

1problem Laplacian(u, uh, init=1) = ...



A preconditioner can be specified in the problem definition:

1problem Laplacian(u, uh, precon=P) = ...

The preconditioning function must have a prototype like:

1func real[int] P(real[int] &xx);


Très grande valeur

The “Très grand valeurtgv (or Terrible giant value) used to implement the Dirichlet conditions can be modified in the problem definition:

1problem Laplacian(u, uh, tgv=1e30) = ...

Refere to Problem definition for a description of the Dirichlet condition implementation.


Pivot tolerance

The tolerance of the pivot in UMFPACK, LU, Crout, Cholesky factorization can be modified in the problem definition:

1problem Laplacian(u, uh, solver=LU, tolpivot=1e-20) = ...



Two specific parameters for the UMFPACK can be modifed:

  • Tolerance of the pivot sym

  • strategy

1problem Laplacian(u, uh, solver=LU, tolpivotsym=1e-1, strategy=0) = ...

Refer to the UMFPACK website for more informations.



Dimension of the Krylov space

Usage of problem is detailled in the tutorials.


Solve type.

Identical to problem but automatically solved.

Usage of solve is detailled in the tutorials.


Variational form type.

1varf vLaplacian (u, uh) = ...

Directly define a variational form.

This is the other way to define a problem in order to directly manage matrix and right hang side.

Usage of varf is detailed in the tutorial.


An array stores multiple objects, and there are 2 kinds of arrays:

  • the first is similar to vector, i.e. array with integer indices

  • the second is array with string indices

In the first case, the size of the array must be known at execution time, and implementation is done with the KN<> class and all the vector operator of KN<> are implemented.

Arrays can be set like in Matlab or Scilab with the operator ::, the array generator of a:c is equivalent to a:1:c, and the array set by a:b:c is set to size \(\lfloor |(b-a)/c|+1 \rfloor\) and the value \(i\) is set by \(a + i (b-a)/c\).

There are int,real, complex array with, in the third case, two operators (.im, .re) to generate the real and imaginary real array from the complex array (without copy).


Quantiles are points taken at regular intervals from the cumulative distribution function of a random variable. Here the array values are random.

This statistical function a.quantile(q) computes \(v\) from an array \(a\) of size \(n\) for a given number \(q\in ]0,1[\) such that:

\[\#\{ i / a[i] < v \} \sim q*n\]

it is equivalent to \(v = a[q*n]\) when the array \(a\) is sorted.

For example, to declare, fill and display an array of real of size n:

1int n = 5;
2real[int] Ai(n);
3for (int i = 0; i < n; i++)
4    Ai[i] = i;
5cout << Ai << endl;

The output of this script is:

      0   1   2   3   4

See the Array example for a complete example.

Array index

Array index can be int or string:

1real[int] Ai = [1, 1, 0, 0];
2real[string] As = [1, 1, 0, 0];

Array size

The size of an array is obtained using the keyword n:

1int ArraySize = Ai.n;

Array sort

To sort an array:


Double array

A double array (matrix) can be defined using two indexes:

1real[int, int] Aii = [[1, 1], [0, 0]];

The two sizes are obtained using the keywords n and m:

1int ArraySize1 = Aii.n;
2int ArraySize2 = Aii.m;

The minimum and maximum values of an array (simple or double) can be obtained using:

1real ArrayMin = Aii.min;
2real ArrayMax = Aii.max;

Th minimum and maximum position of an array can be obtained using:

1int mini = Aii.imin;
2int minj = Aii.jmin;
4int maxi = Aii.imax;
5int maxj = Aii.jmax;


An array can be obtained from a finite element function using:

1real[int] aU = U[];

where U is a finite element function.

Array of FE functions

It is also possible to make an array of FE functions, with the same syntax, and we can treat them as vector valued function if we need them.

The syntax for space or vector finite function is

 1int n = 100; //size of the array.
 2Vh[int] wh(n); //real scalar case
 3Wh[int] [uh,vh](n); //real vectorial case
 4Vh<complex>[int] cwh(n); //complex scalar case
 5Wh<complex>[int] [cuh, cvh](n); //complex vectorial case
 6[cuh[2], cvh[2]] = [x, y]; //set interpolation of index 2
 8// Array of Array
 9real [int][int] V(10);
10matrix[int] B(10);
11real [int, int][int] A(10);



In the following example, Poisson’s equation is solved for 3 different given functions \(f=1,\, \sin(\pi x)\cos(\pi y),\, |x-1||y-1|\), whose solutions are stored in an array of FE function.

 1// Mesh
 2mesh Th = square(20, 20, [2*x, 2*y]);
 4// Fespace
 5fespace Vh(Th, P1);
 6Vh u, v, f;
 8// Problem
 9problem Poisson (u, v)
10   = int2d(Th)(
11        dx(u)*dx(v)
12      + dy(u)*dy(v)
13   )
14   + int2d(Th)(
15      - f*v
16   )
17   + on(1, 2, 3, 4, u=0)
18   ;
20Vh[int] uu(3); //an array of FE function
21// Solve problem 1
22f = 1;
24uu[0] = u;
25// Solve problem 2
26f = sin(pi*x)*cos(pi*y);
28uu[1] = u;
29// Solve problem 3
30f = abs(x-1)*abs(y-1);
32uu[2] = u;
34// Plot
35for (int i = 0; i < 3; i++)
36   plot(uu[i], wait=true);

See FE array example.

Map arrays

1real[string] map; //a dynamic array
3map["1"] = 2.0;
4map[2] = 3.0; //2 is automatically cast to the string "2"
6cout << "map[\"1\"] = " << map["1"] << endl;
7cout << "map[2] = " << map[2] << endl;

It is just a map of the standard template library so no operations on vector are allowed, except the selection of an item.


Defines a sparse matrix.

Matrices can be defined like vectors:

1matrix A = [[1, 2, 3],
2            [4, 5, 6],
3            [7, 8, 9]];

or using a variational form type (see Finite Element):

1matrix Laplacian = vLaplacian(Uh, Uh);

or from block of matrices:

1matrix A1, ..., An;
2matrix A = [[A1, ...], ..., [..., An]];

or using sparse matrix set:

1A = [I, J, C];


I and J are int[int] and C is real[int]. The matrix is defined as:

\[A = \sum_k{C[k]M_{I[k], J[k]}}\]

where \(M_{a, b} = \left(\delta_{ia}\delta_{jb}\right)_{ij}\)

I, J and C can be retrieved using [I, J, C] = A (arrays are automatically resized).

The size of the matrix is n = I.max;, m = J.max;.

Matrices are designed using templates, so they can be real or complex:

1matrix<real> A = ...
2matrix<complex> Ai = ...



See problem.

The default solver is GMRES.

1matrix A = vLaplacian(Uh, Uh, solver=sparsesolver);


1set(A , solver=sparsesolver);



If true, the factorization is done for LU, Cholesky or Crout.

1matrix A = vLaplacian(Uh, Uh, solver=LU, factorize=1);


1set(A , solver=LU, factorize=1);


Stop test

See problem.


Très grande valeur

See problem.



See problem.


Pivot tolerance

See problem.



See problem.



See problem.



Name of the file containing solver parameters, see Parallel sparse solvers



Vector of integer parameters for the solver, see Parallel sparse solvers



Vector of real parameters for the solver, see Parallel sparse solvers



String parameters for the solver, see Parallel sparse solvers


To modify the solver, the stop test,… after the matrix construction, use the set keyword.

Matrix size

The size of a matrix is obtain using:

1int NRows = A.n;
2int NColumns = A.m;

Matrix resize

To resize a matrix, use:

1A.resize(n, m);


When resizing, all new terms are set to zero.

Matrix diagonal

The diagonal of the matrix is obtained using:

1real[int] Aii = A.diag;

Matrix renumbering

1int[int] I(15), J(15);
2matrix B = A;
3B = A(I, J);
4B = A(I^-1, J^-1);

Complex matrix

Use .im and .re to get the imaginary and real part of a complex matrix, respectively:

1matrix<complex> C = ...
2matrix R =;
3matrix I =;

Dot product / Outer product

The dot product of two matrices is realized using:

1real d = A' * B;

The outer product of two matrices is realized using:

1matrix C = A * B'

See Matrix operations example for a complete example.

Matrix inversion

See Matrix inversion example.

Table of content