generalities.tex

FreeFEM is an implementation of the GFEM language dedicated to the finite element method. It provid

%
% $Id: generalities.tex,v 1.4 2001/10/20 23:57:17 prudhomm Exp $
%
% SUMMARY:      
% USAGE:        
%
% AUTHOR:       Christophe Prud'homme <[e-mail]>
% ORG:           MIT
% E-MAIL:       [e-mail]
%
% ORIG-DATE:     8-Feb-97 at 16:45:11
% LAST-MOD: 20-Oct-01 at 15:45:01 by Christophe Prud'homme
%
% DESCRIPTION:  
% This is part of the FreeFEM Documentation Manual
% Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001
%   Christophe Prud'homme and Olivier Pironneau
% See the file fdl.tex for copying conditions.
% DESCRIP-END.

\chapter{The Gory Details}
\label{cha:general}

\section{Programs}
\label{sec:programs}

\textbf{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 \texttt{begin} can be replaced by \texttt{\{} and
\texttt{end} by \texttt{\}}.

\textbf{C programmers}: caution the syntax is "...\};" while most C
constructs use " ...;\}"  


\textbf{ Example 1}: Triangulates a circle and plot $f = x*y$

\begin{Verbatim}
border(1,0,6.28,20)  
 begin 
  x:=cos(t);
  y:=sin(t);
 end; 
buildmesh(200); 
f=x*y; plot(f);
\end{Verbatim}


\textbf{ Example 2}: on the circle solve the Dirichlet problem 

-$\Delta(u) = x*y$ with $u=0$ on $\partial \Omega$

\begin{Verbatim}
border(1,0,6.28,20) 
 begin 
  x:=cos(t); 
  y:=sin(t);
 end; 
buildmesh(200);
solve(u) begin 
 onbdy(1) u =0; 
 pde(u) -laplace(u) = x*y ; 
end; 
plot(u);
\end{Verbatim}


\section{List of reserved words}

\begin{table}[htbp]
\begin{center}
\begin{tabular}[c]{|l|l|}
  \hline
  \multicolumn{1}{|c|}{Keywords} &
  \multicolumn{1}{|c|}{Explanations} \\ \hline
  
  begin, \{ & Begin a new block \\ \hline 
  end, \}   & End a block       \\ \hline 
  
  if, then, else, or, and,iter & Conditionnals and Loops \\ \hline
  
  x,y,t,ib,iv,region,nx,ny & Mesh variables \\ \hline
  
  
  log, exp, sqrt, abs, min, max & \\
  sin, asin, tan, atan, cos, acos & Mathematical Functions\\
  cosh, sinh, tanh & \\  \hline
  
  I, Re, Im & Complex Numbers \\ 
  complex & Enter Complex Number Mode \\ \hline
  
  buildmesh, savemesh, loadmesh, adaptmesh & Mesh related Functions \\
  & build, save , load  and mesh adaptation\\ \hline
  
  one, scal, dx, dy, dxx, dyy, dxy, dyx, convect & Mathematical operators \\ 
  & can be called wherever you want \\ \hline
  
  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 \\ \hline
  
  save, load, saveall & Saving and Loading Functions \\ \hline 
  
  & Change the Wait State: if \texttt{wait}\\ 
  wait, nowait, changewait & then the user must click in \\ 
  & the window to continue \\ \hline
  
  precise, halt, include, evalfct, exec, user & Miscellaneous Functions \\ \hline
  
\end{tabular}
  \caption{Gfem Keywords}
  \label{tab:1}
\end{center}
\end{table}
\section{Building a mesh}

\subsection{Triangulations}


To create a triangulation you must either
\begin{itemize}

\item 
Open a project

\item
Read an old triangulation stored in text format

\item
Execute a program which contains the keyword \verb+buildmesh+ 

\item
Create one by hand drawing the boundary of $\Omega$
and activate the menu \verb+Triangulate+  (Macintosh only).

\end{itemize}

In integrated environments, once created, triangulations can be 
displayed, stored on disk in \textbf{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).

\textbf{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.

%% node  Border(), Geometric variables, Triangulations, Building a mesh
\subsection{Border(), buildmesh(), polygon()}
%%  node-name,  next,  previous,  up

\index{border@\texttt{border}}
\index{buildmesh@\texttt{buildmesh}}
\index{polygon@\texttt{polygon}}
Use it to triangulate domain defined by its boundary.
The syntax is
\begin{Verbatim} 
border(ib,t_min,t_max,nb_t)
 begin  
  ...x:=f(t);
  ...y:=g(t)...
 end;
buildmesh(nb_max);
\end{Verbatim}
where each line with \verb+border+ could be replaced by a line with \verb+polygon+ 
\begin{Verbatim}
polygon(ib,'file_name'[,nb_t]); 
\end{Verbatim}
where \verb+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 \verb+t+ and the two coordinates must be \verb+x+  and \verb+y+ .  When the  parameter
goes from \verb+t_min+ to \verb+t_max+  the boundary must be scanned so as to have $\Omega$
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 \verb+border+ ,  and
have inner dividing curves; \verb+nb_t+  points are put on the boundary with values
$t = t_{min} + i * (t_{max} - t_{min}) / (nb_t-1)$ where \verb+i+  takes  integer values
from \verb+0+ to \verb+nb_t-1+ .
 
 The triangulation is created by a Delaunay-Voronoi method with \verb+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 \verb+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 $\Omega$ 
on their right if they are closed. If they do not carry any boundary
conditions they should be given identification number \verb+ib=0+ .
 
The usual \verb+if... then ... else+  statement can be used with the compound
statement: \verb+begin...end+ . This allows piecewise parametric definitions of 
complicated or polygonal boundaries.

The boundary identification number \verb+ib+  can be overwritten.For
example:
\begin{Verbatim}
border(2,0,4,41) begin
   if(t<=1)then  {  x:=t; y:=0 };
   if((t>1)and(t<2))then {  x:=1; y:=t-1; ib=1 };
   if((t>=2)and(t<=3))then  { x:=3-t; y:=1 };
   if(t>3)then { x:=0; y:=4-t }
end;
buildmesh(400);
\end{Verbatim}

Recall that \verb+begin+ and \verb+{+ is the same and so is \verb+end+  and \verb+}+.
Here one side of the unit square has ib=1. The 3 other sides have ib=2.

 The keyword \verb+polygon+  causes a sequence of boundary points to be read
from the file \verb+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
\begin{Verbatim}
x[0]    y[0]
x[1]    y[1]
....
\end{Verbatim}
each being separated by a tab or a carriage return.  The last
parameter \verb+nb_t+  is optional; it means that each segment will be divided into
\verb+nb_t+1+ equal segments (i.e. \verb+nb_t+  points are added on each segments).

For example 
\begin{Verbatim}
polygon(1,'mypoints.dta',2);
buildmesh(100);
\end{Verbatim}
with the file mypoints.dta containing
\begin{Verbatim}
0.      0.
1.      0.
1.      1.
0.      1.
0.      0.
\end{Verbatim}
triangulates the unit square with 4 points on each side and gives \texttt{ib=1} 
to its boundary.
Note that \texttt{polygon(1,'mypoints.dta')} is like \texttt{polygon(1,'mypoints.dta',0)}\index{polygon@\texttt{polygon}}.


\subsubsection{\texttt{buildmesh} and domain decomposition}

There is a problem with
\texttt{buildmesh}\index{buildmesh@\texttt{buildmesh}} when doing domain
decomposition: by 
default \textbf{Gfem} swap the diagonals\index{diagonal swaping} at the corners of the domain if the
triangle has two boundary edges. This will lead to bad domain
decomposition\index{domain decomposition} at
the sub-domain interfaces.
\par
To solve this, there is a new flag for \texttt{buildmesh} which is optional:


\begin{center}
  \texttt{buildmesh(<max\_number\_of\_vertices>, <flag>)}\\
  where \texttt{<flag>} =
  
    $$    \left\{
      \begin{array}[l]{l}
        = 0 \mbox{ classic way: do diagonal swaping}\\
        = 1 \mbox{ domain decomposition: no diagonal swaping}
      \end{array}
    \right.$$

\end{center}

%% node  Geometric variables, Regions, Border(), Building a mesh
\subsection{Geometric variables, inner and region bdy, normal}
%%  node-name,  next,  previous,  up

\begin{itemize}
\item
\verb+x,y+  refers to the coordinates in the domain
\item
\verb+ib+  refers to the boundary identification number; it is zero inside
 the domain.
\item
\verb+nx+ and \verb+ny+  refer to the x-y components of the normal on the boundary
vertices; it is zero on all inner vertices.
\item 
\verb+region+ refers to the domain identification number
which is itself based on an internal number, ngt, assigned to each
triangle by the triangulation  constructor.
\end{itemize}

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 sub-domains.

%% node  Regions,  , Geometric variables, Building a mesh
\subsection{Regions}
%%  node-name,  next,  previous,  up

A sub-domain is a piece of the domain which is simply connected and 
delimited  by boundaries.

Each sub-domain has a \textbf{region} number assigned to it.  This is done
by  \textbf{Gfem}, not by the user. Every time \verb+border+  is called, an internal 
number \verb+ngt+ is incremented by 1. Then when the key word \verb+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 \verb+ngt+.  The number \textbf{region}  is a piecewise linear continuous interpolation
at the vertices of \verb+ngt+.  To be exact, the value of \textbf{region}  at a
vertex $(x_0,y_0)$ is the value of \texttt{ngt}  at $(x_0,y_0-10^{-6})$,
except if \texttt{precise} is set in which case \textbf{region} is equal to
\texttt{ngt}.

\section{Functions}

\subsection{Functions and scalars}

Functions are either read or created.

\begin{itemize}

\item
Functions can be read from a file if its values at the vertices of the 
triangulation are stored in text format. (Open a .dta example with a text 
editor to see the format).

\item
Functions can be created by executing a  program. An instruction like 
\verb+f=x*y+  really means that $f(x,y)=x*y$

for all \verb+x+ and \verb+y+. Here \verb+x+  and \verb+y+  refer 
to the  coordinates in the domain represented by the triangulation. 

\item
Functions can be created with other previously 
defined functions such as in \verb+g=sin(x*y); f=exp(g);+ .

\item
Four other variables can be used besides \verb+x, y, iv, t+ :  
\verb+ nx, ny, ib, region+.

\end{itemize}
 
Most usual functions can be used:
\begin{Verbatim}
max, min, abs, atan, sqrt,
cos, sin, tan, acos, asin, one,
cosh, sinh, tanh, log, exp
\end{Verbatim} 
\verb+one(x+y<0)+ for instance means \verb+1+ if \verb+x+y<0+  and \texttt{0} 
otherwise.

Operators: 
\begin{Verbatim}
and, or, < , <=, < , >=, ==, +, -, *, /, ^
x^2 means x*x
\end{Verbatim}

Functions created by a program are displayed only if the key word 
\verb+plot()+ or \verb+plot3d()+  is used ( here \verb+plot(f)+  ).
 
Derivatives of functions can be created by the keywords \verb+dx()+  and
\verb+dy()+ . 

Unless \verb+precise+  is set, they are interpolated so the results is also 
continuous piecewise linear (or discontinuous when precise is set).
Similarly the convection operator \texttt{convect(f,u1,u2,dt)}
\index{convect@\texttt{convect}}  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 \verb+:=+  is useful to create them, as in
\begin{Verbatim}
 a:= (1+sqrt(5))/2; 
 f= x*cos(a*pi*y);
\end{Verbatim}
Here \verb+f+ is a function, \verb+a+  is a scalar and \verb+pi+  is a (predefined)
a scalar. 

It is possible to evaluate a function at a point as in
\verb+a:=f(1,0)+
Here the value of \verb+a+ will be \verb+1+ because \verb+f(1,0)+ means
\verb+f+  at \verb+x=1+  and \verb+y=0+ .

%% node  Building functions, Value of a function at one point, Functions and scalars, Functions
\subsection{Building functions}
%%  node-name,  next,  previous,  up

There are 6 predefined functions: \verb+x,y,ib,region, nx, ny+ .
\begin{itemize}
\item 
The coordinate \verb+x+ is horizontal and \verb+y+  is vertical.
\item

    $
    \mbox{\texttt{ib}} = \left\{
      \begin{array}[l]{l}
        0 \mbox{ inside } \Omega \\
        > 0 \mbox{ on } \partial \Omega
      \end{array}
    \right.
    $
  
On $\partial \Omega$ it is equal to the boundary identification number.
\end{itemize}

The usual \texttt{if... then ... else} statement can be used with an
important restriction on the logical expression which must
return a \textbf{scalar} value:

\begin{Verbatim}
if( logical  expression) then
 {
   statement;
   ....;
   statement;
 } 
else
 {
   .....
 };
\end{Verbatim}
The \texttt{logical expression} controls the \texttt{if}  by its return being 0 or >0,
it is evaluated only once (i.e. with \texttt{x}, \texttt{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 \verb+:=+  tells the compiler not
to allocate a data array to this variable.  Thus \verb+v=sin(a*pi*x);+ 
generates an array for \verb+v+ but no array is assigned to \verb+a+  in the
statement \verb+a:=2+ .

\subsection{Value of a function at one point}

If \verb+f+ has been defined earlier then it is possible to write \verb+a:=f(1.2,3.1);+ 
Then a has the value of \verb+f+ at \verb+x=1.2+  and \verb+y=3.1+ .

It is also allowed to do
\begin{Verbatim}
x1:=1.2;