Skip to content

Types

Standard types#

int#

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

1
int i = 0;

bool#

Boolean value.

1
bool b = true;

The result of a comparison is a boolean

1
bool b = (1 < 2);

real#

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

1
real r = 0.;

complex#

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

1
complex c = 0. + 1i;
The imaginary number i is defined as 1i

Example

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

The output of this script is:

1
2
3
4
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.

1
string s = "this is a string";

Note

string value is enclosed within double quotes.

Other types can be concatenate to a string, like:

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

To append a string in a string at position 4:

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

To copy a substring in an other string:

1
string s2 = s1(5:10);

See String Example for a complete example.

Mesh design#

border#

Border type.

1
border b(t=0., 1.){x=cos(2.*pi*t); y=sin(2.*pi*t);}
Define the 2D geometrical border in parametric coordinates.

Label

A label can be defined with the border:

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

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

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

mesh#

2D Mesh type (see Mesh Generation).

1
mesh Th;

mesh3#

3D mesh type (see Mesh Generation).

1
mesh3 Th;

Finite element space design#

fespace#

Finite element space type (see Finite Element).

1
2
fespace Uh(Th, P1);
fespace 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:

1
2
3
4
5
fespace Uh(Th, P1);
Uh u;

fespace UPh(Th, [P2, P2, P1]);
UPh [Ux, Uy, p];

Macro design#

macro#

Macro type.

1
2
macro vU() [Ux, Uy] //
macro 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:

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

See Macro example

NewMacro / EndMacro#

Warning

In developement - Not tested

Set and end a macro

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

IFMACRO#

Check if a macro exists and check its value.

1
2
3
4
5
6
7
IFMACRO(AA) //check if macro AA exists
...
ENDIFMACRO

IFMACRO(AA, tt) //check if amcro exists and is equall to tt
...
ENDIFMACRO

ENDIFMACRO#

Functions design#

func#

Function type.

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

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

Note

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

Function with parameters:

1
2
3
func real f (real var){
    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 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.

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.

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

1
problem Laplacian (u, uh) = ...
FreeFem++ needs the variational form in the problem definition.

In order to solve the problem, just call:

1
Laplacian;

Solver

A solver can be specified in the problem definition:

1
problem 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:

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

If \varepsilon>0, the stop test is: Else, the stop test is:

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.

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

Preconditioning

A preconditioner can be specified in the problem definition:

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

The preconditioning function must have a prototype like:

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

Très grande valeur

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

1
problem 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:

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

UMFPACK

Two specific parameters for the UMFPACK can be modifed:

  • Tolerance of the pivot sym
  • strategy
1
problem Laplacian(u, uh, solver=LU, tolpivotsym=1e-1, strategy=0) = ...

Refer to the UMFPACK website for more informations.

dimKrylov

Dimension of the Krylov space

Usage of problem is detailled in the tutorial.

solve#

Solve type.

Identical to problem but automatically solved.

Usage of solve is detailled in the tutorial.

varf#

Variational form type.

1
varf 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: 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:

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

The output of this script is:

1
2
5
      0   1   2   3   4

See the Array example for a complete example.

Array index#

Array index can be int or string:

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

Array size#

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

1
int ArraySize = Ai.n;

Array sort#

To sort an array:

1
Ai.sort;

Double array#

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

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

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

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

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

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

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

1
2
3
4
5
int mini = Aii.imin;
int minj = Aii.jmin;

int maxi = Aii.imax;
int maxj = Aii.jmax;

Tip

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

1
real[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

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

// Array of Array
real [int][int] V(10);
matrix[int] B(10);
real [int, int][int] A(10);

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
 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
33
34
35
36
// Mesh
mesh Th = square(20, 20, [2*x, 2*y]);

// Fespace
fespace Vh(Th, P1);
Vh u, v, f;

// Problem
problem Poisson (u, v)
    = int2d(Th)(
          dx(u)*dx(v)
        + dy(u)*dy(v)
    )
    + int2d(Th)(
        - f*v
    )
    + on(1, 2, 3, 4, u=0)
    ;

Vh[int] uu(3); //an array of FE function
// Solve problem 1
f = 1;
Poisson;
uu[0] = u;
// Solve problem 2
f = sin(pi*x)*cos(pi*y);
Poisson;
uu[1] = u;
// Solve problem 3
f = abs(x-1)*abs(y-1);
Poisson;
uu[2] = u;

// Plot
for (int i = 0; i < 3; i++)
    plot(uu[i], wait=true);

See FE array example.

Map arrays#

1
2
3
4
5
6
7
real[string] map; //a dynamic array

map["1"] = 2.0;
map[2] = 3.0; //2 is automatically cast to the string "2"

cout << "map[\"1\"] = " << map["1"] << endl;
cout << "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:

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

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

1
matrix Laplacian = vLaplacian(Uh, Uh);

or from block of matrices:

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

or using sparse matrix set:

1
A = [I, J, C];

Note

I and J are int[int] and C is real[int]. The matrix is defined as: 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:

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

Solver

See problem.

The default solver is GMRES.

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

Factorize

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

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

Stop test

See problem.

Très grande valeur

See problem.

Preconditioning

See problem.

Pivot tolerance

See problem.

UMFPACK

See problem.

dimKrylov

See problem.

datafilename

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

lparams

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

dparams

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

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:

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

Matrix resize#

To resize a matrix, use:

1
A.resize(n, m);

Warning

When resizing, all new terms are set to zero.

Matrix diagonal#

The diagonal of the matrix is obtained using:

1
real[int] Aii = A.diag;

Matrix renumbering#

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

Complex matrix#

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

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

Dot product / Outer product#

The dot product of two matrices is realized using:

1
real d = A' * B;

The outer product of two matrices is realized using:

1
matrix C = A * B'

See Matrix operations example for a complete example.

Matrix inversion#

See Matrix inversion example.