Gfem is a small language which generally follows the syntax of the language Pascal. See the list below for the reserved words of this language.
The reserved word begin can be replaced by { and end by }.
C programmers: caution the syntax is "...};" while most C constructs use " ...;}"
Example 1: Triangulates a circle and plot f = x*y

Example 2: on the circle solve the Dirichlet problem
Δ(u) = x*y with u=0 on ∂ Ω

Keywords Explanations begin, { Begin a new block end, } End a block if, then, else, or, and,iter Conditionnals and Loops x,y,t,ib,iv,region,nx,ny Mesh variables log, exp, sqrt, abs, min, max sin, asin, tan, atan, cos, acos Mathematical Functions cosh, sinh, tanh I, Re, Im Complex Numbers complex Enter Complex Number Mode buildmesh, savemesh, loadmesh, adaptmesh Mesh related Functions build, save , load and mesh adaptation one, scal, dx, dy, dxx, dyy, dxy, dyx, convect Mathematical operators can be called wherever you want solve Enter Solver Mode pde, id, dnu laplace, div onbdy Solver Functions plot, plot3d graphical functions: plot isolines in 2D and Elevated Surface in 3D save, load, saveall Saving and Loading Functions Change the Wait State: if wait wait, nowait, changewait then the user must click in the window to continue precise, halt, include, evalfct, exec, user Miscellaneous Functions
To create a triangulation you must either
buildmesh
Triangulate
(Macintosh only).In integrated environments, once created, triangulations can be displayed, stored on disk in Gfem format or in text format or even a zoom of its graphic display can be stored in Postscript format (to include it in a TeX file for example).
Gfem stores for each vertex its number and boundary identification number and for each triangle its number and region identification number. Edges number are not stored, but can be recovered by program if needed.
Use it to triangulate domain defined by its boundary. The syntax is

where each line with border
could be replaced by a line with polygon

where f,g
are generic functions and the [...] denotes an optional addition.
The boundary is given in parametric form. The name of the parameter must
be t
and the two coordinates must be x
and y
. When the parameter
goes from t_min
to t_max
the boundary must be scanned so as to have Ω
on its left, meaning counter clockwise if it is the outer boundary and
clockwise if it is the boundary of a hole. Boundaries must be closed but
they may be entered by parts, each part with one instruction border
, and
have inner dividing curves; nb_t
points are put on the boundary with values
t = t_{min} + i * (t_{max} − t_{min}) / (nb_{t}−1) where i
takes integer values
from 0
to nb_t1
.
The triangulation is created by a DelaunayVoronoi method with nb_max
vertices at most. The size of the triangles is monitored by the size of
the nearest boundary edge. Therefore local refinement can be achieved by
adding inner artificial boundaries.
Triangulation may have boundaries with several connected components.
Each connected component is numbered by the integer ib
.
Inner boundaries (i.e. boundaries having the domain on
both sides) can be useful either to separate regions or to
monitor the triangulation by forcing vertices and edges in it. They
must be oriented so that they leave Ω
on their right if they are closed. If they do not carry any boundary
conditions they should be given identification number ib=0
.
The usual if... then ... else
statement can be used with the compound
statement: begin...end
. This allows piecewise parametric definitions of
complicated or polygonal boundaries.
The boundary identification number ib
can be overwritten.For
example:

Recall that begin
and {
is the same and so is end
and }
.
Here one side of the unit square has ib=1. The 3 other sides have ib=2.
The keyword polygon
causes a sequence of boundary points to be read
from the file file_name
which must be in the same directory as the
program. All points must be in sequential order and describing part of
the boundary counter clockwise; the last one should be equal to the first
one for a closed curve.
The format is

each being separated by a tab or a carriage return. The last
parameter nb_t
is optional; it means that each segment will be divided into
nb_t
1+ equal segments (i.e. nb_t
points are added on each segments).
For example

with the file mypoints.dta containing

triangulates the unit square with 4 points on each side and gives ib=1 to its boundary. Note that polygon(1,’mypoints.dta’) is like polygon(1,’mypoints.dta’,0).
There is a problem with buildmesh when doing domain decomposition: by default Gfem swap the diagonals at the corners of the domain if the triangle has two boundary edges. This will lead to bad domain decomposition at the subdomain interfaces.
To solve this, there is a new flag for buildmesh which is optional:
⎧ ⎨ ⎩ 

x,y
refers to the coordinates in the domain
ib
refers to the boundary identification number; it is zero inside
the domain.
nx
and ny
refer to the xy components of the normal on the boundary
vertices; it is zero on all inner vertices.
region
refers to the domain identification number
which is itself based on an internal number, ngt, assigned to each
triangle by the triangulation constructor.
Inner boundaries which divide the domain into several simply connected components are useful to define piecewise discontinuous functions such as dielectric constants, Young modulus...
Inner boundaries may meet other boundaries only at their vertices. Such inner boundaries will split the domain in several subdomains.
A subdomain is a piece of the domain which is simply connected and delimited by boundaries.
Each subdomain has a region number assigned to it. This is done
by Gfem, not by the user. Every time border
is called, an internal
number ngt
is incremented by 1. Then when the key word border
is
invoked the last edge of this portion of boundary assigns this number
to the triangle which lies on its left. Finally all triangles which are in
this subdomain are reassigned this number.
At the end of the triangulation process, each triangle has a well defined
number ngt
. The number region is a piecewise linear continuous interpolation
at the vertices of ngt
. To be exact, the value of region at a
vertex (x_{0},y_{0}) is the value of ngt at (x_{0},y_{0}−10^{−6}),
except if precise is set in which case region is equal to
ngt.
Functions are either read or created.
f=x*y
really means that f(x,y)=x*yfor all x
and y
. Here x
and y
refer
to the coordinates in the domain represented by the triangulation.
g=sin(x*y); f=exp(g);
.x, y, iv, t
:
nx, ny, ib, region
.Most usual functions can be used:

one(x
y<0)+ for instance means 1
if x
y<0+ and 0
otherwise.
Operators:

Functions created by a program are displayed only if the key word
plot()
or plot3d()
is used ( here plot(f)
).
Derivatives of functions can be created by the keywords dx()
and
dy()
.
Unless precise
is set, they are interpolated so the results is also
continuous piecewise linear (or discontinuous when precise is set).
Similarly the convection operator convect(f,u1,u2,dt)
defines a new
function which is approximately
f(x−u1(x,y)dt , y−u2(x,y)dt) 
Scalars are also helpful to create functions. Since no data array is
attached to a scalar the symbol :=
is useful to create them, as in

Here f
is a function, a
is a scalar and pi
is a (predefined)
a scalar.
It is possible to evaluate a function at a point as in
a:=f(1,0)
Here the value of a
will be 1
because f(1,0)
means
f
at x=1
and y=0
.
There are 6 predefined functions: x,y,ib,region, nx, ny
.
x
is horizontal and y
is vertical.
ib = {
0 inside Ω 
> 0 on ∂ Ω 
On ∂ Ω it is equal to the boundary identification number.
The usual if... then ... else statement can be used with an important restriction on the logical expression which must return a scalar value:

The logical expression controls the if by its return being 0 or >0, it is evaluated only once (i.e. with x, y being the coordinates of the first vertex, if there are functions inside the logical expression). Auxiliary variables can be used.
In order to minimize the memory the symbol :=
tells the compiler not
to allocate a data array to this variable. Thus v=sin(a*pi*x);
generates an array for v
but no array is assigned to a
in the
statement a:=2
.
If f
has been defined earlier then it is possible to write a:=f(1.2,3.1);
Then a has the value of f
at x=1.2
and y=3.1
.
It is also allowed to do

Remark: Recall that ,a
being a scalar, its value is appended to the
file a.dta.
dx(f)
is the partial derivative of f
with respect to x
; the result
is piecewise constant when precise
is set and interpolated with
mass lumping as a the piecewise linear function when precise is not set.
Note that dx()
is a non local operator so statements like f=dx(f)
would give the wrong answer because the new value for f
is place before
the end of the use of the old one.
The Finite Element Method does not handle convection terms properly when
they dominate the viscous terms: upwinding is necessary; convect
provides a tool for Lagrangian upwinding. By
g=convect(f,u,v,dt)
Gfem construct an approximation of
f(x−u1(x,y)dt,y−u2(x,y)dt) 
Recall that when
 + u 
 + v 
 = lim_{dt → 0} 

Thus to solve
 + u 
 + v 
 −div(µ grad f) = g, 
in a much more stable way that if the operator dx(.)
and dy(.)
were use, the following scheme can be used:

Remark: Note that convect
is
a nonlocal operator. The statement f = convect(f,u,v,dt) would give
an incorrect result because it modifies f
at some vertex where the
old value is still needed later. It is necessary to do

It is important to understand the difference between
a global operator such as dx()
and a local operator
such as sin()
.
To compute dx(f)
at vertex q
we need f at all neighbors
of q
. Therefore evaluation of dx(2*f) require
the computation of 2*f
at all neighbor vertices of
q
before applying dx() ; but in which memory would the
result be stored? Therefore Gfem does not allow this
and forces the user to declare a function g =2*f
before evaluation of dx(g)
; Hence in

the equal sign forces the evaluation of g
at all
vertices, then when the second equal signs forces the
evaluation of the expression on the right of h at
all vertices , everything is ready for it.
Global operators are

Example of forbidden expressions:

intt
returns a complex or real number, an integral with respect to x,y
int
returns a complex or real number, an integral on a curve
intt[f]
or intt(n1)[f]
or intt(n1,n2)[f]
or intt(n1,n2,n3)[f]
int(n1)[f]
or int(n1,n2)[f]
or int(n1,n2,n3)[f]
where n1,n2,n3
are boundary or subdomain identification numbers and where
f
is an array function.

int
and intt
are global operators, so the values of the integrands are
needed at all vertices at once, therefore you can’t put an expression
for the integrand, it must be a function.
Be careful to check that the region number are correct when you use intt(n)[f]
.
Unfortunately freefem does not store the edges numbers. Hence there are
ambiguities at vertices which are at the intersections of 2 boundaries.
The following convention is used:
int(n)[g]
computes the integral of g
on all segments of the boundary (both ends
have id boundary number !=0) with one vertex boundary id number = n.
(Remember that you can control the boundary id number of the boundary ends
by the order in which you place the corresponding border
call or by an
extra argument in border
)
Its purpose is to define a boundary condition for a Partial Differential Equation (PDE).
The general syntax is

where ib
’s are boundary identification numbers, <expression>
is a generic
expression and g a generic function.
The term id(u)
may be absent as in dnu(u)=g
.
dnu(u)
represents the conormal of the PDE, i.e.
± 
 u . n 
when the PDE operator is
a * u − 
 .(± 
 u) 
The syntax for pde is

It defines a partial differential equation with non constant
coefficients where op
is one of the operator:
id()
dx()
dy()
laplace()
dxx()
dxy()
dyx()
dyy()
and where [*/]
means either a *
or a /
and similarly for ±.
Note that the expressions are necessarily AFTER the operator while in
practice they are between the 2 differentiations for laplace...dyy
.
Thus laplace(u)*(x+y) means
∇.((x+y)∇ u)
.Similarly dxy(u)*f means
∂ f ∂ u/∂ y/∂ x.
The syntax for a single unknown function u solution of a PDE is

For 2systems and the use of solve(u,v)
, see the section 2Systems
.
It defines a PDE and its boundary conditions.
It will be solved by the Finite Element Method of degree 1 on triangles
and a Gauss factorization.
Once the matrix is built and factorized solve
may be called again by
solve(u,1)...;
then the matrix is not rebuilt nor factorized and
only a solution of the linear system is performed by an up and a down
sweep in the Gauss algorithm only. This saves a lot of CPU time whenever
possible. Several matrices can be stored and used simultaneously, in
which case the sequence is

where i
is a scalar variable (not an array function).
However matrices must be constructed in the natural order: i=1
first
then i=2....
after they can be reused in any order. One can also
reuse an old matrix with a new definition, as in

Notice that solve(u)
is equivalent to solve(u,1)
.
Remark: 2Systems have their own matrices, so they do not count in the previous ordering.
Before going to systems make sure that your 2 pde’s are indeed coupled and that no simple iteration loop will solve it, because 2systems are significantly more computer intensive than scalar equations.
Systems with 2 unknowns can be solved by

The syntax for solve
is the same as for scalar PDEs; so
solve(u,v,1)
is ok for instance. The
equations above can be in any orders; several onbdy()
can be used in
case of multiple boundaries...
The syntax for onbdy
is the same as in the scalar case; either Dirichlet
or mixedNeumann, but the later can have more than one id()
and
only one dnu()
.
Dirichlet is treated as if it was mixed Neumann with a small coefficient. For instance u=2 is replaced by dnu(u)+1.e10*u=2.e10 , with quadrature at the vertices.
Conditions coupling u,v
are allowed for mixed Neumann
only, such as id(u)+id(v)+dnu(v)=1. (As said later this is an
equation for v
).
In solve(u,v,i) begin .. end;
when i>0
the linear system is built
factorized and solved. When i<0
, it is only solved; this is useful when
only the right hand side in the boundary conditions or in the equations
have change. When i<0
, i
refers to a linear system i>0
of
SAME TYPE, so that scalar systems and 2systems have their own count.
Remark: saveall('filename',u,v)
works also.
The syntax for pde()
is the same as for the scalar case. Deciding which
equation is an equation for u
or v
is important in the resolution of the
linear system (which is done by Gauss block factorization) because some
block may not be definite matrices if the equations are not well chosen.
onbdy(...) ... dnu(u) ... = ...;
is a
boundary condition associated to u
, even if it contains id(v)
. onbdy(...) u...=...;
is also
associated with u
(the syntax forbids any v
operator in this case).u
is the array function in a pde(u)
then what follows is the
PDE associated to u
. Corners where boundary conditions change from Neumann to Dirichlet are ambiguous because Dirichlet conditions are assigned to vertices while Neumann conditions should be assigned to boundary edges; yet Gfem does not give access to edge numbers. Understanding how these are implemented helps overcome the difficulty.
All boundary conditions are converted to mixed Fourier/Robin conditions:

For Dirichlet conditions a is set to 1.0e12
and c
is multiplied by the same;
for Neumann a=0
. Thus Neumann condition is present even when there
is Dirichlet but the later overrules the former because of the large
penalty number. Functions a,b,c
are piecewise linear continuous, or
discontinuous if precise
is set.
In case of DirichletNeumann corner (with Dirichlet on one side and Neumann
on the other) it is usually better to put a Dirichlet logic at the corner.
But if fine precision is needed then the option precise
can guarantee that
the integral on the edge near the corner on the Neumann side is properly taken
into account because then the corner has a Dirichlet value and a Neumann value
by the fact that functions are discontinuous.
The new keyword varsolve allows the user to enter PDEs in weak form. Syntax:

where
<unknown function list>
and
<test function list>
are one or two function names separated by a comma.
<int>
is a positive or negative integer
<expression>
is an expression returning a real or
complex number
We have used the notation << >>
whenever the entities
can be omitted.
Examples

How does it works The interpreter evaluates the expression after the ":" for each triangle and for each 3 vertices; if there is an instruction prior the ":" it is also evaluated similarly. Each evaluation is done with one of the unknown and one of the test functions being 1 at one vertices and zero at the 2 others. This will give an element of the contribution of the triangle to the linear system of the problem. The right hand side is constructed by having all unknowns equal to zero and one test function equal to one at one vertex. whenever integrals appear they are computed on the current triangle only.
Note that varsolve
takes longer than solve
because derivatives like dx(u)
are evaluated 9 times
instead of once.
Within a program the keyword plot(u)
will display u
.
Instruction save(’filename’,u) will save the data array u
on
disk. If u
is a scalar variable then the (single) value of u
is
appended to the file (this is useful for time dependent problems or any
problem with iteration loop.).
Instruction savemesh(’filename’) will save the triangulation on disk.
Similarly for reading data with load(’filename’,u) and loadmesh(’filename’).
The file must be in the default directory, else it won’t be found. The
file format is best seen by opening them with a text editor. For a data
array f
it is:

(ns
is the number of vertices)
If f
is a constant, its single value is appended to the end of the file;
this is useful for time dependent problems or any problem with iteration
loop.
If precise
is set still the function stored by save
is interpolated on the vertices as the P^{1} continuous function given by mass lumping (see above).
For triangulations the file format is (nt
= number of triangles):

Gfem uses the Fortran standard for me[][]
and
numbers the vertices starting from number 1 instead of 0 as in the
Cstandard. Thus in Cprograms one must use me[][]1
.
Other formats are also recognized by freefem via their file name extensions
for our own internal use we have defined .amdba
and .am_fmt
.
You can do the same if your format is not ours.

There is an optional arguments for the functions load, save, loadmesh, savemesh. This is the 2nd or 3rd argument of these functions. Here are the prototypes:

As an example see nsstepad.pde which use this feature to save the mesh and the solution at each adaptation of the mesh. This special feature allows you to save or load a generic filename with a counter, the final filename is built like this ’<generic filename><counter>’.
The purpose is to solve all the data for a PDE or a 2system with only one instruction. It is meant for those who want to write their own solvers.
The syntax is:

The syntax is exactly the same as that of solve(,)
except that the
first parameter is the file name. The other parameters are used only
to indicate to the interpreter which is/are the unknown function.
The file format for the scalar equation (laplace is decomposed on nuxx, nuyy)

is that each line has all the values for x,y being a vertex: f, g, p, b, c, a1, a2, nuxx, nuxy, nuyx, nuyy.
The actual routine is in C++

The same function is used for complex coefficients, by overloading the operator <<:

For 2systems also the same is used with

where N=2
.
A Dirichlet condition is applied whenever p[k](?). ( Dirichlet
conditions with value 0
are changed to value 1e10
)
Gfem supports other interesting features:

where j refers to the number of loops; j can be the result of an expression (as in iter(i*k)).
Imbedded loops are not allowed. You can use iter
with the adaptation
features of Gfem.
Gfem can handle complex coefficients with 4 dedicated keywords:
complex
: to tell Gfem that complex number will be used.
When it is used it must be located at the beginning of the program
before any function declarations, otherwise the results will be incorrect.
It can appear more than once in the program but only the first occurrence
counts.There is purposely no conjug function but barz=Re(z)I*Im(z) will do.
By default all graphics display the real part. To display the imaginary part do plot(Im(f)).
The functions implemented for complex numbers are:
The linear systems for the PDE are solved by a Gauss complex LU factorization.
WARNING: failure to declare complex
in the program implies all
computation will be done in real, even if I
is used.
The instruction a:=scal(f,g);
does
a =  ∫ 
 f(x,y) 
 (x,y) dxdy 
where Ω is the triangulated domain.
Whenever there is a plot
command, Gfem stops to let the user see
the result. By using nowait
no stop will be made;
Remark under X11: If you click the right button in the window, the next time the solver will give the hand to the plotter the program will stop.
This function is only available under integrated environments.
The last function defined by the keyword plot
is displayed. It can
be visualized in several fashion, one of which being a one dimensional
plot along any segment defined by the mouse. Selection of this menu
brings causes Gfem to waits for the user input which should be the line
segment on which the function is to be displayed. Thus one should
press the mouse at the beginning point then drag the mouse and release
the button at the final point
It is safe to click in the window after to check that the function display is correct. What is seen is a
t → f(x(t),y(t)) 
Plot where [x(t),y(t)]_{t} is the segment drawn by the user and f
is
the last function displayed in the plot window.
The abscissa is the distance with the beginning point.
This keyword warns Gfem that precise quadrature formula must be used. Hence arrayfunctions are discretized as piecewise linear discontinuous functions on the triangulation. Then all integrals are computed with 3 inner Gauss points slightly inside each triangle.
This option consumes more memory (3nt instead of ns per functions, i.e. 9 times more approximately, but still it is nothing compared with the memory that a matrix of the linear system of a PDE requires) because each arrayfunction is stored by 3 values on each triangles.
It is a good idea to use it in conjunction with convect
and/or discontinuous nonlinear terms.
exec('prog_name')
will launch the application prog_name
. It is
useful to execute an external PDE solver for instance especially under
Unix. It is not implemented for Macintosh because there is no simple way
to return to MacGfem after progname has ended.
The same can be achieved manually by a suitable combination of saveall
, wait
and
load
and simultaneous execution under multifinder.
user(what,f)
calls the C++ function in a jloop:

creal
is a scalar (float
) or a complex number if complex
has
been set; nquad=3*nt
if precise
is set and ns
otherwise. what
is intended for users who need several such functions.
Then all can be put in the super function user
and selection is by
an if statement on what
.gfemuser
access to all global variables are of course possible:
the triangulation (ns,nt, me, q, ng,ngt, area...
) ... refer to the file
fem.C
for more details.Remark: An example of such gfemuser function is in fem.C
; if you
wish to put your own you must compile and link it. Under Unix, it is
easy. Under the Macintosh system, either you use freefem, which is
Gfem’s kernel, or you must ask us a library version of MacGfem.
instruction
.The syntax analysis is
driven by the syntax rules because the language is LL(1). Thus there
is Cfunction for each nonterminal.
x * y
would generate a tree with root *
and two branches x
and y
. Trees are Cstruct with four pointers
for the 4 branches (here two would point to NULL) and a symbol for the
operator. The Cfunction which builds trees is called plante
.
In the end the program is transformed into a full tree and
to execute the program there is only one thing to do: evaluate
the operator at the root of the tree.
eval
. It
looks at the root symbol and perform the corresponding operation. Here
it would do: return "value pointed by L1"*"value pointed by L2" if L1
and L2 where the addresses of the two branches.
x
is required, this is also done
by a tree which has the operator "oldvar" as root. The trickiest of all is
the compound instruction {...;....};
.
Here {
is considered as an operator with one branch on the current
instruction and one branch on the next one. Similarly for the
if...then... else
instruction.Suppose one wants to add an instruction to freefem, here is what must be done:
installe
. As there will be a
new Symbol, update the list of symbols (an enum structure). Add the
Cfunctions for the syntax analysis according to the diagrams. Modify
eval by adding to the switch
the new case.
Here is the very simple structure we used for the nodes of the tree:
