manual-full.tex

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

The language it
defines is typed, polymorphic and reentrant with macro generation (see
\ref{macro}).  Every variable must be typed and declared in a
statement; each statement separated from the next by a semicolon
`\texttt{;}'.

For purposes of explanation, we used $\mathcal{T}_h$ (\texttt{Th}),\, $V_h$ (\texttt{Vh}), unknown function $u$ (\texttt{u}), test functions $v$ (\texttt{v}) and the problem \textbf{Poisson}, etc. (the term inside the parentheses are symbols in \freefempp programming), but
you can use \emph{any name except reserved words and names already used}.
Reserved words are shown in blue. \texttt{pi, x, y, label, solver} are
reserved variables. It is allowed (although not advisable) to redefine
these variables, so they will not be highlighted again in the
following example programs.\\\\

In each step, the independence in \freefempp programming is very high
as stated below.
\begin{itemize}
  \item
  For example, by changing 1st and 2nd lines as
  following, we can solve (\ref{eqn:Poisson}) and (\ref{eqn:Dirichlet})
  in L-shape domain with $\Gamma=\sum_{j=1}^6\Gamma_j$.
\bT
@border G1(t=0,1){x=t;y=0;@label=1;}  // $\Gamma_1$
@border G2(t=0,0.5){x=1;y=t;@label=1;}  // $\Gamma_2$
@border G3(t=0,0.5){x=1-t;y=0.5;@label=1;}  // $\Gamma_3$
@border G4(t=0.5,1){x=0.5;y=t;@label=1;}  // $\Gamma_4$
@border G5(t=0.5,1){x=1-t;y=1;@label=1;}  // $\Gamma_5$
@border G6(t=0,1){x=0;y=1-t;@label=1;}  // $\Gamma_6$
@mesh Th = @buildmesh ( G1(6)+G2(4)+G3(4)+G4(4)+G5(4)+G6(6));
\eT
  \item
  In Step 3, you can control where the solution will be approximated.
  If you write ``\texttt{Vh(Th,P1);}'' in 3rd line, you can get $P_1$-approximation. The machine time by $P_1$-element will be faster than $P_2$-element
  and the storage is less.
  \item
  In Step 4, you can change the equation and boundary conditions easily.
  For example, if you want solve
  \begin{eqnarray*}
  -\textrm{div}(k(x,y)\nabla u(x,y))&=&f(x,y)\quad\textrm{in }\Omega\\
  u(x,y)&=&0\qquad \textrm{if $(x,y)$ on }\Gamma
  \end{eqnarray*}
  you only write as follows
  \bT
  @func f= $\cdots$;
  @func k= 1+sin(2*pi*x)*cos(2*pi*y);
  @problem Poisson(u,v) =
     @int2d(Th)(k*(@dx(u)*@dx(v) + @dy(u)*@dy(v)))
     - @int2d(Th)( f*v)
     + @on(1,u=0)  ;
  \eT
  \emph{The user can use FE-function as the given function $f$ } (see \refSec{FE-function}), for example, obtained function \texttt{u} in Example \ref{exm:first}.
  In \freefempp programming, the easy reuse of the obtained results
  is important feature.
  \item
  The user can easily compare between mathematical modelling and \freefempp program.
\end{itemize}

\subsection{Projection or Interpolation}
\index{P1-projection}
For a finite element space $V_h$, $P_1$-projection $\Pi_h$ is defined by
\[
\Pi_h f=f(q^1)\phi_1+f(q^2)\phi_2+\cdots +f(q^{n_v})\phi_{n_v}
\]
for all continuous functions $f$. In \freefempp we can easily create the
projection $f_h(=\texttt{fh})$ by
\bT
Vh fh = $f(x,y)$;
\eT
$\Pi_h$ is also called $P_1$-interpolate.

\subsection{\setS{Matrix and Vector}}
Here, we show how to get the stiffness matrix using \freefempp
The first command \texttt{varf} is to define the \key{variational formula}.
\begin{example}\label{exm:second}
Here we solve the same problem (\ref{eqn:Poisson}) and (\ref{eqn:Dirichlet})
using matrix.
For purposes of explanation, we chage mesh size and use $P_1$-element:
\bT
 1: @border C(t=0,2*pi) { x = cos(t); y = sin(t); }
 2: @mesh Th = @buildmesh(C(7));  // changed from Example \ref{exm:first}
 3: @fespace Vh(Th,P1);
 4: Vh u,v,f,F;
 5: @varf a(u,v) = @int2d(Th)( @dx(u)*@dx(v) + @dy(u)*@dy(v))
 6:             + @on(C,u=0) ; // see (\ref{eqn:Stiffness0})
 7: @varf b([v],[f]) = @int2d(Th)(v*f);
 8:
 9: f = x*y;  // interpolate $(x,y) \longleftarrow x*y$ function
10: @matrix A=a(Vh,Vh);  // stiffness matrix, see (\ref{eqn:Equation})
11: @matrix B=b(Vh,Vh);
12: F[]=B*f[];     // load vector, see (\ref{eqn:Equation})
13: @cout << "F=" << F[] << @endl;
14: @cout <<"A="<< A << @endl;
15: u[]=A^-1*F[];  // solve $A\vec{U}_h = \vec{F}$, see (\ref{eqn:Equation})
16: @plot(u);
\eT

We get the mesh $\mathcal{T}_h=\{T_1,\cdots,T_7\}$ (see Fig. \ref{fig:secondTh}).
\begin{note}
\label{note:mesh}
In what follows, we denote the vertices by $q^i,\, i=1,\cdots,8$,
the number of vertices by $n_v$, the number of triangles by $n_t$.
For each triangle $T_k\in \mathcal{T}_h$, we index the vertices by
$q^{k_1},\, q^{k_2},\, q^{k_3}$ and denote the edges by
$[q^{k_1}, q^{k_2}]$, $[q^{k_2}, q^{k_3}]$, $[q^{k_3}, q^{k_i}]$, that is,
$[q^i,q^j]$ is the segment connecting $q^i$ and $q^j$.
We denote the number of edges $[q^i,q^j]$
by $n_e$ for all $q^i,\, q^j\partial\Omega_h$,
$\Omega_h=\sum_{k=1}^7T_k$.
Here $n_v=8$, $n_t=7$, $n_e=7$.
\end{note}
\begin{figure}[htbp]
\begin{minipage}{\textwidth}
\begin{minipage}{0.3\textwidth}
\includegraphics[width=\textwidth]{secondT}%
\caption{mesh \texttt{Th}}
\label{fig:secondTh}
\end{minipage}
\hspace{0.5mm}
\begin{minipage}{0.7\textwidth}
\includegraphics[width=\textwidth]{hat}%
\caption{Graph of $\phi_1$ (left hand side) and $\phi_6$}
\label{fig:hatFunction}
\end{minipage}
\end{minipage}
\end{figure}


The function $v$ in ``\texttt{Vh}'' is expressed
\begin{equation*}
v(x)=v_{1}\varphi _{1}(x) +\cdots +v_{n_{v}}\varphi _{n_{v}}(x)
\end{equation*}%
using the \key{hat functions} $\varphi _{j},\,j=1,\cdots ,n_{v}$ (see Fig.
\ref{fig:hatFunction}).
Here the $j$-th hat function $\varphi _{j}$ associated with $j$-th
vertex $q^j$ is defined in the following way:
\begin{enumerate}
\item $\varphi _{j}$ is continuous function on $\Omega_h$.

\item $\varphi _{j}$ is linear on each triangle $T_k,\, k=1,\cdots,n_{t}$
 of ``\texttt{Th}''.

\item $\varphi _{j}\left( {q^{i}}\right) =\delta _{ji}$ where $q^{i}$
denotes the $i$-th vertex, for all $i=1,\cdots ,n_{v}$.
\end{enumerate}
Here $\delta_{ij}$ is the Kronecker symbol.

\begin{note}
Other finite element spaces in \freefempp are explained in \refSec{Finite Elements}.
\end{note}
\begin{note}
\label{note:FE2VEC}
For an element $v=v_1\phi_1+\cdots+v_M\phi_M$ in
a finite element space $V_h$, we get the \emph{column vector}
$\{v\}$
\[
\{v\}=\left[
\begin{array}{c}
v_1\\
\vdots\\
v_M
\end{array}
\right]\qquad
\{v\}=\texttt{v[]}\quad \textrm{in \freefempp}
\]
\end{note}
Theoretically, it is natural to use the finite element space
\[
H^1_{0h}=
\left\{v\in V_h(\mathcal{T}_h,P_1)\left|
\phi_i(x)=0\quad \textrm{if }q^i\in \partial\Omega_h
\right.\right\}
\]
Let $I_{\Omega}$ be the set of indices $i$ of all internal vertices of the mesh
\texttt{Vh}. In this example, $I_{\Omega}=\{6\}$.
The stiffness matrix $A$ in 10th line is:
\begin{eqnarray}
\label{eqn:StiffnessDirichlet}
\borderarray{[}{]}{1em}{1.2ex}{rrrrrrrrr}{
&1&2&3&4&5&6&7&8\\
1&10^{30}& -0.31& 0& 0&-0.46& -0.51& 0&0\\
2&-0.31& 10^{30}& -0.23& 0& 0&-0.71& 0&0\\
3&0&-0.23&10^{30}&-0.31& 0&-0.71& 0&0\\
4&0& 0&-0.31&10^{30}& 0& -0.51& -0.46&0 \\
5&-0.46& 0& 0& 0&10^{30}& -0.35& 0&-0.54\\
6&-0.51& -0.71& -0.71& -0.51& -0.35& 3.47& -0.35& -0.30\\
7&0& 0& 0&-0.46& &-0.35& 10^{30}& -0.54\\
8&0& 0& 0& 0&-0.54&-0.30&-0.54& 10^{30}
}
\end{eqnarray}
that is
\begin{eqnarray}
A_{ij}&=&\int_{\Omega_h}\nabla u_j\cdot \nabla u_i\quad
\textrm{if }i\neq j,\, i=j\in I_{\Omega}\\
A_{ij}&=&E\qquad (E=10^{30})\quad \textrm{if }j\in I_{\Omega}.
\end{eqnarray}
The load matrix $F^T$ is:
\begin{eqnarray*}
\left(
\begin{array}{cccccccc}
-0.020& -0.037& 0.037& 0.020&0.064& 0& -0.064&1.\times 10^{-17}
\end{array}
\right)
\end{eqnarray*}
For $i\not\in I_{\Omega}$,
\[
Eu_i+\sum_{i\neq j}A_{ij}u_j=b_i
\]
which means that
\[
u_i=(b_i-\sum_{i\neq j}A_{ij}u_j)\times E^{-1}\simeq 10^{-30}\simeq 0
\]
\end{example}

Mathematical results indicate that the Poisson equation
with Neaumann boundary condition has not unique solution,
whose weak form is same to (\ref{eqn:Poisson}) except the boundary condition:
\[
\int_{\Omega}\nabla u\cdot \nabla v = \int_{\Omega}fvdx
\]
without pernarization $E$. Then the stiffness matrix is created by
\bT
@varf a(u,v) = @int2d(Th)( @dx(u)*@dx(v) + @dy(u)*@dy(v))
@matrix A=a(Vh,Vh); // stiffness matrix
\eT
and the obtained stiffness matrix is the following
\[
\borderarray{[}{]}{1em}{1.2ex}{rrrrrrrrr}{
&1&2&3&4&5&6&7&8\\
1&1.29&-0.31&0&0&-0.46&-0.51&0&0\\
2&-0.31&1.26&-0.23&0&0&-0.71&0&0\\
3&0&-0.23&1.26&-0.31&0&-0.71&0&0\\
4&0&0&-0.31&1.29&0&-0.51&-0.46&0\\
5&-0.46&0&0&0&1.35&-0.35&0&-0.54\\
6&-0.51&-0.71&-0.71&-0.51&-0.35&3.47&-0.35&-0.30\\
7&0&0&0&-0.46&0&-0.35&1.35&-0.54\\
8&0&0&0&0&-0.54&-0.30&-0.54&1.38
}
\]
The determinant of this matrix is $-1.7082274230870981\times 10^{-9}\approx 0$
(The matrix here differ from original one by omitting from third decimal decimal point).

\subsubsection{Non-homogeneous Dirichlet Condition}
If we want solve the problem
\[
-\Delta u=f\quad \textrm{in }\Omega;\qquad
u=g\quad \textrm{on }\partial\Omega
\]
We rewrite Example \ref{exm:first} as
\bT
 5: @func f= x*y; @func g = sin(pi*x)*cos(pi*y);
 6: @problem Poisson(u,v,@solver=LU) =
 7:    @int2d(Th)(@dx(u)*@dx(v) + @dy(u)*@dy(v))   //  bilinear part
 8:    - @int2d(Th)( f*v)          // right hand side
 9:    + @on(1,u=g)  ;  // Non-homogeneous Dirichlet
\eT
This make the following linear system, for $i\not\in I_{\Omega}$,
\[
Eu_i+\sum_{i\neq j}A_{ij}u_j=b_i+Eg(q^i)
\]
which means that
\[
u_i=g(q^i)+(b_i-\sum_{i\neq j}A_{ij}u_j)\times E^{-1}\simeq g(q^i)+O(1/E)
\]
\begin{note}
To solve non-homogeneous Dirichlet, we rewrite Example \ref{exm:second} as
\bT
 4: Vh u,v,f,F,g,bc; g = sin(pi*x)*cos(pi*y);
 5: @varf a(u,v) = @int2d(Th)( @dx(u)*@dx(v) + @dy(u)*@dy(v))
 6:             + @on(C,u=1) ; // see (\ref{eqn:Stiffness0})
10: @matrix A=a(Vh,Vh); bc[]=a(0,Vh);
12: F[]=B*f[];  F[] += bc[] .* g[];
\eT
Here ``\texttt{bc[]=a(0,Vh)}'' create the vector
$[bc_1,bc_2,\cdots,bc_M]$, $bc_i=0$ if $i\in I_{\Omega}$ and
$bc_i=E (=10^{30})$ if $i\not\in I_{\Omega}$.
If the finite approximation of $g$ is $g\approx g_1\phi_1+\cdots+g_M\phi_M$
\begin{equation}
\texttt{bc[] .* g[]}=\sum_{j=0}^{M-1} bc_jg_j
\end{equation}
\end{note}

\subsubsection{\setS{Matrix Operations}}

The multiplicative operators *, /, and \% group left to right.

\begin{itemize}
\item  \verb@'@  is unary right transposition of array, matrix \index{transpose}
 \item \verb@.*@ is the term to term multiply operator. \index{.*} \index{\string'} \index{divide!term to term}
 \item \verb@./@ is the term to term divide operator. \index{./} \index{\string'} \index{product!term to term}
\end{itemize}
there are some compound operator:
\begin{itemize}
 \item \verb@^-1@ is for  solving the linear system (example: \verb$ b = A^-1 x$) \index{solve!linear system}
 \item \verb@' *@ is the compound  of transposition and matrix product, so it is the dot product
(example \verb$real DotProduct=a'*b$) \index{dot product}\index{product!dot}
\end{itemize}
\begin{example}~
\bFF
@mesh Th = @square(2,1);
@fespace Vh(Th,P1);
Vh f,g;
f = x*y;
g = sin(pi*x);
Vh<complex> ff,gg; // a complex valued finite element function \index{FE function!complex}
ff= x*(y+1i);
gg = exp(pi*x*i);
@varf mat(u,v) =
  int2d(Th)(1*dx(u)*dx(v)+2*dx(u)*dy(v)+3*dy(u)*dx(v)+4*dy(u)*dy(v))
  + on(1,2,3,4,u=1);
@varf mati(u,v) =
  int2d(Th)(1*dx(u)*dx(v)+2i*dx(u)*dy(v)+3*dy(u)*dx(v)+4*dy(u)*dy(v))
  + on(1,2,3,4,u=1);
@matrix A = mat(Vh,Vh); @matrix<complex> AA = mati(Vh,Vh); // a
complex sparse matrix \index{matrix:complex}

Vh m0; m0[] = A*f[];
Vh m01; m01[] = A'*f[];
Vh m1; m1[] = f[].*g[];
Vh m2; m2[] = f[]./g[];
@cout << "f = " << f[] << @endl;
@cout << "g = " << g[] << @endl;
@cout << "A = " << A << @endl;
@cout << "m0 = " << m0[] << @endl;
@cout << "m01 = " << m01[] << @endl;
@cout << "m1 = "<< m1[] << @endl;
@cout << "m2 = "<< m2[] << @endl;
@cout << "dot Product = "<< f[]'*g[] << @endl;
@cout << "hermitien Product = "<< ff[]'*gg[] << @endl;
\eFF
This produce the following:
\begin{eqnarray*}
A&=&\borderarray{[}{]}{1em}{1.2ex}{rrrrrrrr}{
&1&2&3&4&5&6\\
1&10^{30}& 1 0.5& 0&3 0.& 4 -2.5& 0\\
2&0.&10^{30}&0.5& 0&0.5&-2.5\\
3&0&0.&10^{30}& 0& 0&0.5\\
4&0.5& 0& 0& 10^{30}& 0.& 0\\
5&-2.5&0.5& 0&0.5&10^{30}&0.\\
6&0&-2.5&0.& 0&0.5& 10^{30}
}
\\
\{v\}=\texttt{f[]}&=&
\left(
\begin{array}{rrrrrr}
0 & 0 & 0 & 0 & 0.5 & 1
\end{array}
\right)^T\\
\{w\}=\texttt{g[]}&=&
\left(