FreeFEM Documentation on GitHub

stars - forks

Types

Standard types

int

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

1int i = 0;


bool

Boolean value.

1bool b = true;


Tip

The result of a comparison is a boolean

bool b = (1 < 2);

real

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

1real r = 0.;


complex

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

1complex c = 0. + 1i;


The imaginary number $$i$$ is defined as 1i

Tip

Example

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)


Note

See Complex example for a detailed example.

string

String value.

1string s = "this is a string";


Note

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

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.

Note

Label

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;}


Note

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);}


Note

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];}


mesh

2D Mesh type (see Mesh Generation).

1mesh Th;


mesh3

3D mesh type (see Mesh Generation).

1mesh3 Th;


Finite element space design

fespace

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:

Generic:

• 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;
3
4fespace UPh(Th, [P2, P2, P1]);
5UPh [Ux, Uy, p];


Macro design

macro

Macro type.

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


Macro ends with //.

Note

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


NewMacro / EndMacro

Warning

In developement - Not tested

Set and end a macro

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


IFMACRO

Check if a macro exists and check its value.

1IFMACRO(AA) //check if macro AA exists
2...
3ENDIFMACRO
4
5IFMACRO(AA, tt) //check if amcro exists and is equall to tt
6...
7ENDIFMACRO


Functions design

func

Function type.

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

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


Note

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;
3}


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$$.

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

Note

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.

FE-functions

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.

Warning

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

Problem design

problem

Problem type.

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


FreeFEM needs the variational form in the problem definition.

In order to solve the problem, just call:

1Laplacian;


Note

Solver

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.

Note

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||}$

Note

Reconstruction

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) = ...


Note

Preconditioning

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);


Note

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.

Note

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) = ...


Note

UMFPACK

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.

Note

dimKrylov

Dimension of the Krylov space

Usage of problem is detailled in the tutorials.

solve

Solve type.

Identical to problem but automatically solved.

Usage of solve is detailled in the tutorials.

varf

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.

Array

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).

Note

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:

5
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:

1Ai.sort;


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;
3
4int maxi = Aii.imax;
5int maxj = Aii.jmax;


Tip

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
7
8// Array of Array
9real [int][int] V(10);
10matrix[int] B(10);
11real [int, int][int] A(10);


Tip

Example

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]);
3
4// Fespace
5fespace Vh(Th, P1);
6Vh u, v, f;
7
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   ;
19
20Vh[int] uu(3); //an array of FE function
21// Solve problem 1
22f = 1;
23Poisson;
24uu[0] = u;
25// Solve problem 2
26f = sin(pi*x)*cos(pi*y);
27Poisson;
28uu[1] = u;
29// Solve problem 3
30f = abs(x-1)*abs(y-1);
31Poisson;
32uu[2] = u;
33
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
2
3map["1"] = 2.0;
4map[2] = 3.0; //2 is automatically cast to the string "2"
5
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.

matrix

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];


Note

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 = ...


Note

Solver

See problem.

The default solver is GMRES.

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


or

1set(A , solver=sparsesolver);


Note

Factorize

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

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


or

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


Note

Stop test

See problem.

Note

Très grande valeur

See problem.

Note

Preconditioning

See problem.

Note

Pivot tolerance

See problem.

Note

UMFPACK

See problem.

Note

dimKrylov

See problem.

Note

datafilename

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

Note

lparams

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

Note

dparams

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

Note

sparams

String parameters for the solver, see Parallel sparse solvers

Tip

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);


Warning

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 = C.re;
3matrix I = C.im;


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.