Types

Standard types

int

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

int i = 0;

bool

Boolean value.

bool b = true;

Tip

The result of a comparison is a boolean

bool b = (1 < 2);

real

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

real r = 0.;

complex

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

complex c = 0. + 1i;

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

Tip

Example

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:

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.

string s = "this is a string";

Note

string value is enclosed within double quotes.

Other types can be concatenate to a string, like:

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:

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

To copy a substring in an other string:

string s2 = s1(5:10);

See String Example for a complete example.

Mesh design

border

Border type.

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

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

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

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

mesh

2D Mesh type (see Mesh Generation).

mesh Th;

mesh3

3D mesh type (see Mesh Generation).

mesh3 Th;

Finite element space design

fespace

Finite element space type (see Finite Element).

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:

fespace Uh(Th, P1);
Uh u;

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

Macro design

macro

Macro type.

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

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

See Macro example

NewMacro / EndMacro

Warning

In developement - Not tested

Set and end a macro

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

IFMACRO

Check if a macro exists and check its value.

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

func f = x^2 + y^2;

Note

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

Function with parameters:

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

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

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.

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.

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

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.

Note

Stop test

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

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

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

Note

Preconditioning

A preconditioner can be specified in the problem definition:

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

The preconditioning function must have a prototype like:

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

Note

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

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

problem 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

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

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:

\[\#\{ 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:

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:

5
      0   1   2   3   4

See the Array example for a complete example.

Array index

Array index can be int or string:

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:

int ArraySize = Ai.n;

Array sort

To sort an array:

Ai.sort;

Double array

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

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

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

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

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

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

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

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:

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

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

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.

// 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

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:

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

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

matrix Laplacian = vLaplacian(Uh, Uh);

or from block of matrices:

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

or using sparse matrix set:

A = [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:

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

Note

Solver

See problem.

The default solver is GMRES.

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

set(A , solver=sparsesolver);

Note

Factorize

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

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

set(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:

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

Matrix resize

To resize a matrix, use:

A.resize(n, m);

Warning

When resizing, all new terms are set to zero.

Matrix diagonal

The diagonal of the matrix is obtained using:

real[int] Aii = A.diag;

Matrix renumbering

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:

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

Dot product / Outer product

The dot product of two matrices is realized using:

real d = A' * B;

The outer product of two matrices is realized using:

matrix C = A * B'

See Matrix operations example for a complete example.

Matrix inversion

See Matrix inversion example.

Here you can use the search bar to search in all FreeFEM examples instead

The list of all matching examples on the right is updated for each new query

Clicking on a search hit or on an example on the right brings up the corresponding script from GitHub

You can also filter the list on the right by selecting one or more intersecting categories on the left tree

Types

Standard types

int

bool

real

complex

string

Mesh design

border

mesh

mesh3

Finite element space design

fespace

Macro design

macro

NewMacro / EndMacro

IFMACRO

ENDIFMACRO

Functions design

func

Elementary functions

Random functions

FE-functions

Problem design

problem

solve

varf

Array

Array index

Array size

Array sort

Double array

Array of FE functions

Map arrays

matrix

Matrix size

Matrix resize

Matrix diagonal

Matrix renumbering

Complex matrix

Dot product / Outer product

Matrix inversion