generalities.tex

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

y1:=1.6; 
a:=f(x1,2*y1); 
save('a.dta',a);
\end{Verbatim}
\textbf{Remark:} Recall that ,\verb+a+ being a scalar, its value is \textbf{appended}  to the 
file \texttt{a.dta}.

%% node  Special functions
\subsection{Special functions:dx(), dy(), convect()}

\index{convect@\texttt{convect}}
\index{dx@\texttt{dx}}
\index{dy@\texttt{dy}}

\verb+dx(f)+ is the partial derivative of \verb+f+  with respect to \verb+x+ ; the result
is piecewise constant when \verb+precise+  is set and interpolated with
mass lumping as a the piecewise linear function when precise is not set.

%$$ dx(u)^i = {\frac{1}{\Sigma_i\Sigma_{k,q^i in T_k} u^k |T_k|}}$$

%where $-q^i$ is a vertex,$dx(u)^i$ the value at this vertex, $\Sigma_k$
%is the area of all the triangles which have $q^i$ for vertex. $-T_k$ is
%a triangle, $|T_k|$ its area, $u^k$ the mean of u on this triangle. 


Note that \verb+dx()+ is a non local operator so statements like \verb+f=dx(f)+ 
would give the wrong answer because the new value for \verb+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; \verb+convect+
\index{convect@\texttt{convect}} provides a tool for Lagrangian upwinding. By
\verb+g=convect(f,u,v,dt)+ \textbf{Gfem} construct an approximation of 

$$f(x-u1(x,y)dt,y-u2(x,y)dt)$$

Recall that when 

$$\frac{\partial f}{\partial t} + u\frac{\partial f}{\partial x} + v\frac{\partial f}{\partial y}
    = lim_{dt \to 0}\frac{f(x,y,t) - f(x-u(x,y)dt,y-v(x,y)dt,t-dt)}{dt}$$

Thus  to solve

$$
\frac{\partial f}{\partial t} + u\frac{\partial f}{\partial x}
        + v\frac{\partial f}{\partial y} -div(\mu grad f) = g,
$$

in a much more stable way that if the operator \verb+dx(.)+ and \verb+dy(.)+
were use, the following scheme can be used:
 
\begin{Verbatim}
iter(n) begin
   id(f)/dt - laplace(f)*mu =g + convect(oldf,u,v,dt)/dt; 
   oldf = f
end;
\end{Verbatim}


\textbf{Remark:} Note that \verb+convect+ \index{convect@\texttt{convect}}  is
a nonlocal operator. The statement \texttt{f = convect(f,u,v,dt)}  would give
an incorrect  result because it modifies \verb+f+  at some vertex where the
old value is still needed later. It is necessary to do
\begin{Verbatim}
g=convect(f,u,v,dt); 
f=g;
\end{Verbatim}

%% node  Integrals, Solving an equation, Functions, Generalities

\section{Global operators}
\label{sec:globop}

It is important to understand the difference between
a global operator such as \verb+dx()+  and a local operator
such as \verb+sin()+ .

To compute \verb+dx(f)+ at vertex \verb+q+  we need \texttt{f}  at all neighbors
of \verb+q+.  Therefore evaluation of \texttt{dx(2*f)}  require
the computation of \verb+2*f+  at all neighbor vertices of
\verb+q+ before applying \texttt{dx()} ; but in which memory would the
result be stored?  Therefore Gfem does not allow this
and forces the user to declare a function \verb+g =2*f+ 
before evaluation of  \verb+dx(g)+ ;  Hence in
\begin{Verbatim}
g = 2*f;
h = dx(g) * dy(f);
\end{Verbatim}

the equal sign forces the evaluation of \verb+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
\begin{Verbatim}
dx(), dy(), convect(), intt[], int()[]
\end{Verbatim}

Example of forbidden expressions:
\begin{Verbatim}
intt[f+g], dx(dy(g)), dx(sin(f)), convect(2*u...)
\end{Verbatim}

\section{Integrals}
%%  node-name,  next,  previous,  up

\begin{itemize}
\item 
\textbf{effect}: 

\verb+intt+ returns a complex or real number, an integral with respect to \verb+x,y+ 
\verb+int+  returns a complex or real number, an integral on a curve
\item 
\textbf{syntax}:   

\verb+intt[f]+ or \verb+intt(n1)[f]+ or \verb+intt(n1,n2)[f]+  or \verb+intt(n1,n2,n3)[f]+ 
\verb+int(n1)[f]+ or \verb+int(n1,n2)[f]+  or \verb+int(n1,n2,n3)[f]+ 

where \verb+n1,n2,n3+  are boundary or subdomain identification numbers and where
\verb+f+  is an array function.

\begin{Verbatim}
border(1)...  end;  /* a border has number 1 */
... buildmesh(...);

f = 2 * x;

/* 
 * nx,ny are the components of the boundary normal 
 */
g = f * (nx + ny);  

/*
 * can't do r:= int[2*x] 
 */
r:= int[f];         
s:=int(1)[g];

/*
 * this is the only way to display the result 
 */
save('r.dta',r);
save('s.dta',s);
\end{Verbatim}

\item
\textbf{Restrictions}:

\verb+int+ and \verb+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 \verb+intt(n)[f]+ .

Unfortunately \textbf{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:
\verb+int(n)[g]+ computes the integral of \verb+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 \verb+border+  call or by an
extra argument in \verb+border+ )

\end{itemize}


%% node  Solving an equation, Results, Integrals, Generalities
\section{Solving an equation}
%%  node-name,  next,  previous,  up



%% Onbdy()::                     
%% Pde()::                       
%% Solve()::                     
%% 2-Systems::                   
%% Boundary conditions at corners::  


%% node  Onbdy(), Pde(), Solving an equation, Solving an equation
\subsection{Onbdy()}

Its purpose is to define a boundary condition for a Partial Differential
Equation (PDE).

The general syntax is
\begin{Verbatim}
onbdy(ib1, ib2,...) id(u)+<expression>*dnu(u) = g
onbdy(ib1, ib2,...) u = g
\end{Verbatim}
where \verb+ib+'s are boundary identification numbers, \verb+<expression>+  is a generic 
expression and \texttt{g}  a generic function. 
 
The term \verb+id(u)+ may be absent as in \verb+-dnu(u)=g+ .
\verb+dnu(u)+  represents the conormal of the PDE, i.e. 
  $$\pm\vec\nabla u . n $$

when the PDE operator is 
  $$  a * u - \vec\nabla.(\pm\vec\nabla u)$$

%% node  Pde(), Solve(), Onbdy(), Solving an equation
\subsection{Pde()}
%%  node-name,  next,  previous,  up

The syntax for pde is
\begin{Verbatim}
pde(u) [+-] op1(u)[*/]exp1 [+-] op2(u)[*/]exp2...=exp3
\end{Verbatim}

It defines a partial differential equation with non constant
coefficients where \verb+op+  is one of the operator:
\begin{itemize}
\item 
\verb+id()+ 
\item 
\verb+dx()+ 
\item
\verb+dy()+ 
\item
\verb+laplace()+ 
\item
\verb+dxx()+ 
\item
\verb+dxy()+ 
\item
\verb+dyx()+ 
\item
\verb+dyy()+ 
\end{itemize}
and where \verb+[*/]+ means either a \verb+*+  or a \verb+/+  and similarly for $\pm$.
Note that the expressions are necessarily \textbf{AFTER} the operator while in
practice they are between the 2 differentiations for \verb+laplace...dyy+ .
Thus \texttt{laplace(u)*(x+y)} means 

  $\nabla.((x+y)\nabla u)$

.Similarly
\texttt{dxy(u)*f} means 

  $\frac{\partial f \frac{\partial u}{\partial y}}{\partial x}$.


%% node  Solve(), 2-Systems, Pde(), Solving an equation
\subsection{Solve()}
%%  node-name,  next,  previous,  up<sect>  Solve()
\index{solve@\texttt{solve}}
The syntax for a single unknown function u solution of a PDE is
\begin{Verbatim}
solve(u) 
 begin 
   onbdy()...; 
   onbdy()...; 
   ...; 
   pde(u)... 
 end; 
\end{Verbatim}

For 2-systems and the use of \verb+solve(u,v)+, see the section \verb+2-Systems+ .
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 \verb+solve+  may be called again by
\verb+solve(u,-1)...;+\index{solve@\texttt{solve}} 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
\begin{Verbatim}
solve(u,i)...;
... 
solve(u,-i)...;
\end{Verbatim}
where \verb+i+  is a scalar variable (not an array function).

However matrices must be constructed in the natural order: \verb+i=1+  first 
then \verb+i=2....+  after they can be re-used in any order.  One can also 
re-use an old matrix with a new definition, as in
\begin{Verbatim}
solve(u,i)...;
... 
solve(u,i)...;
solve(u,\pm i)...;
\end{Verbatim}
Notice that \verb+solve(u)+ is equivalent to \verb+solve(u,1)+\index{solve@\texttt{solve}} .

\textbf{Remark:} 2-Systems have their own matrices, so they do not count in the previous
ordering.

%% node  2-Systems, Boundary conditions at corners, Solve(), Solving an equation
\subsection{2-Systems}
%%  node-name,  next,  previous,  up

Before going to systems make sure that your 2 pde's are indeed coupled and
that no simple iteration loop will solve it, because  2-systems are
significantly more computer intensive than scalar equations.

Systems with 2 unknowns can be solved by
\begin{Verbatim} 
solve(u,v) 
begin 
 onbdy(..) ...dnu(u)...=.../* defines a bdy condition for u */
 onbdy(..) u =...          /* defines a bdy conditions for v */
 pde(u)  ...               /* defines PDE for u */
 onbdy(..)<v=... or ...dnu(v)...> /* defines bdy conditions for v */
 pde(v)  ...               /* defines PDE for u */
end; 
\end{Verbatim}

The syntax for \verb+solve+ is the same as for scalar PDEs; so
\verb+solve(u,v,1)+\index{solve@\texttt{solve}}  is ok for instance.  The
equations above can be in any orders; several \verb+onbdy()+  can be used in
case of multiple boundaries... 

The syntax for \verb+onbdy+  is the same as in the scalar case; either Dirichlet 
or mixed-Neumann,  but the later can have more than one \verb+id()+  and 
only one \verb+dnu()+ .

Dirichlet is treated as if it was mixed Neumann with a small coefficient.
For instance u=2 is replaced by  \texttt{dnu(u)+1.e10*u=2.e10} , with quadrature
at the vertices.  
 
Conditions coupling \verb+u,v+  are allowed for mixed Neumann
only, such as \texttt{id(u)+id(v)+dnu(v)=1}. (As said later this is an 
equation for \verb+v+ ).

In \verb+solve(u,v,i) begin .. end;+ when \verb+i>0+  the linear system is built 
factorized and solved.  When \verb+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 \verb+i<0+, \verb+i+  refers to a linear system \verb+i>0+  of
\textbf{SAME TYPE}, so that scalar systems and 2-systems have their own count.

\textbf{Remark:} \verb+saveall('filename',u,v)+   works also.

The syntax for \verb+pde()+  is the same as for the scalar case.  Deciding which
equation is an equation for \verb+u+ or \verb+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.
\begin{itemize}
\item
A boundary condition like \verb+ onbdy(...) ... dnu(u) ... = ...; +  is a
boundary condition associated to \verb+u+, even if it contains \verb+id(v)+ .

\item
Obviously a boundary condition like \verb+ onbdy(...) u...=...;+  is also
associated with \verb+u+ (the syntax forbids any \verb+v+ -operator in this case).

\item
If \verb+u+ is the array function in a \verb+pde(u)+  then what follows is the
PDE associated to \verb+u+ . 
  
\end{itemize}

%% node  Boundary conditions at corners,  , 2-Systems, Solving an equation
\subsection{Boundary conditions at corners}
%%  node-name,  next,  previous,  up


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 
\textbf{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:
\begin{Verbatim}
id(u) a + dnu(u) b = c;
\end{Verbatim}

For Dirichlet conditions a is set to \verb+1.0e12+ and \verb+c+  is multiplied by the same;
for Neumann \verb+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 \verb+a,b,c+  are piecewise linear continuous, or
discontinuous if \verb+precise+  is set.

In case of Dirichlet-Neumann 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 \verb+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.

%% node  Results, Other features, Solving an equation, Generalities

\subsection{Weak formulation}
\label{sec:varform}

The new keyword \texttt{varsolve} allows the user to enter PDEs in
weak form. Syntax:
\begin{Verbatim}
varsolve(<unknown function list>;<test function list>
          ,<<int>>) <<instruction>> : <expression>>;
\end{Verbatim}

where
\begin{itemize}
\item 
\verb+<unknown function list>+ and
\item 
\verb+<test function list>+
are one or two function names separated by a comma.
\item
\verb+<int>+ is a positive or negative integer
\item
instruction is one instruction or more if they are
   enclosed within begin end or \texttt{\{\}}
\item
\verb+<expression>+ is an expression returning a real or
   complex number
\end{itemize}

We have used the notation \verb+<< >>+ whenever the entities
can be omitted.
Examples
\begin{Verbatim}
varsolve(u;w)  /* Dirichlet problem -laplace(u) =x*y */
begin
        onbdy(1) u = 0;
        f = dx(u)*dx(w) + dy(u)*dy(w)
        g = x*y;
end : intt[f] - intt[g,w];

varsolve(u;w,-1)  /* same  with prefactorized matrix */
begin
        onbdy(1) u = 0;
        f = dx(u)*dx(w) + dy(u)*dy(w)
        g = x*y;
end : intt[f] - intt[g,w];

varsolve(u;w)  /* Neuman problem u-laplace(u) = x*y */
begin
        f = dx(u)*dx(w) + dy(u)*dy(w) -x*y;
        g = x;
end : intt[f] + int[u,w] - int(1)[g,w];

varsolve(u,v;w,s) /* Lame's equations */
begin
   onbdy(1) u=0;
   onbdy(1) v=0;
   e11 = dx(u);
   e22 = dy(v);
   e12 = 0.5*(dx(v)+dy(u));
   e21 = e12;
   dive = e11 + e22;
   s11w=2*(lambda*dive+2*mu*e11)*dx(w);
   s22s=2*(lambda*dive+2*mu+e22)*dy(s);
   s12s = 2*mu*e12*(dy(w)+dx(s));
   s21w = s12s;
   a = s11w+s22s+s12s+s21w +0.1*s;
end : intt[a];
\end{Verbatim}