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

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

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

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

**needs the variational form in the problem definition.**

`FreeFem++`

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

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

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.