manual-full.tex

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

stored in an array of FE function.
\bT
@mesh Th=@square(20,20,[2*x,2*y]);
@fespace Vh(Th,P1);
Vh u, v, f;
@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
f=1;   // problem1
Poisson; uu[0] = u;
f=sin(pi*x)*cos(pi*y);  // problem2
Poisson; uu[1] = u;
f=abs(x-1)*abs(y-1);    // problem3
Poisson; uu[2] = u;
@for (int i=0; i<3; i++)  // plots all solutions
  @plot(uu[i], @wait=true);
\eT
\end{example}

 For the second case, it is just
 a map of the STL\footnote{Standard template Library, now part of standard \Cpp}\cite{cpp}
 so no vector operation except the
 selection of an item is allowed . \index{dot product}\index{transpose}

The transpose operator is \texttt{\string'} like  MathLab or SciLab, so the way
to compute the dot product of two array \ttCC{a,b} is  \ttCC{@real\ ab=\ a'*b}\index{product!dot}.
\index{min}\index{array!min}
\bFF
@int i;
@real [int] tab(10), tab1(10); // 2 array of 10 real
  \it  real [int] tab2;    //  bug array with no size
tab = 1;                //  set all the array to 1
tab[1]=2;
@cout << tab[1] << " " << tab[9] << " size of tab = "
     << tab.n << " " << tab.min << " " << tab.max << " " <<  @endl;
tab1=tab; tab=tab+tab1; tab=2*tab+tab1*5; tab1=2*tab-tab1*5;
tab+=tab; @cout << " dot product " << tab'*tab << @endl; //
${}^{t}{tab}\,{tab} $ @cout << tab << @endl; @cout << tab[1] << " "
<< tab[9] <<  @endl; @real[string] map;        //  a dynamic array
@for (i=0;i<10;i=i+1)
  {
    tab[i] = i*i;
    @cout << i << " " << tab[i] << "\n";
  };

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;
\eFF


\subsection{\setS{Loops}}

The \texttt{for} and \texttt{while}  loops are implemented
with \texttt{break} and \texttt{continue} keywords.
\index{for}\index{while}
\index{break}\index{continue}

In for-loop, there are three parameters;
the INITIALIZATION of a control variable,
the CONDITION to continue,
the CHANGE of the control variable.
While \texttt{CONDITION} is true, for-loop continue.
\bT
@for (INITIALIZATION; CONDITION; CHANGE)
     { BLOCK of calculations }
\eT
The sum from 1 to 10 is calculated by (the result is in \ttCC{sum}),
\bT
@int sum=0;
@for (@int i=1; i<=10; i++)
   sum += i;
\eT
The while-loop
\bT
@while (CONDITION) {
   BLOCK of calculations or change of control variables
}
\eT
is executed repeatedly until CONDITION become false.
The sum from 1 to 10 is also written by \ttCC{@while} as follows,
\bT
@int i=1, sum=0;
@while (i<=10) {
  sum += i; i++;
}
\eT

We can exit from a loop in midstream by \ttCC{@break}.
The \ttCC{@continue} statement will pass the part from
\emph{continue} to the end of the loop.

\begin{example}~
\bFF
@for (@int i=0;i<10;i=i+1)
    @cout << i << "\n";
@real eps=1;
@while (eps>1e-5)
 { eps = eps/2;
   @if( i++ <100) @break;
   @cout << eps << @endl;}

@for (int j=0; j<20; j++) {
   @if (j<10) @continue;
   @cout << "j = " << j << @endl;
}
\eFF
\end{example}

\subsection{\setS{Input/Output}}
\index{cout}\index{cint}\index{ifstream}\index{ofstream}\index{endl}

The syntax of input/output statements is similar  to \Cpp syntax. It
uses \ttCC{@cout}, \ttCC{@cin}, \ttCC{@endl}, \ttCC{<<,>>}.

To write  to (resp. read from)  a file, \index{$<<$} \index{$>>$}\index{append}\index{ofstream!append}
declare a new variable \texttt{ofstream ofile("filename");} or \texttt{ofstream ofile("filename",append);} (resp.
\texttt{ifstream ifile("filename");} ) and use \texttt{ofile}  (resp. \texttt{ifile})
as \texttt{cout} (resp. \texttt{cin}). The word \texttt{append} in  \texttt{ofstream ofile("filename",append);}
 means openning a file in append mode.

\begin{note} The file is closed
at the exit of the enclosing block,
\end{note}
\begin{example}~
\label{exm:io}
\bFF
@int i;
@cout << " std-out" << @endl;
@cout << " enter i= ? ";
@cin >> i ;
{
  @ofstream f("toto.txt");
  f << i << "coucou'\n";
}; //  close the file f because the variable f is delete

{
  @ifstream f("toto.txt");
   f >> i;
}
{
  @ofstream f("toto.txt",append);
     // to append to the existing file "toto.txt"
  f << i << "coucou'\n";
}; //  close the file f because the variable f is delete

  @cout << i << @endl;
\eFF
\end{example}

\section{\setS{Mesh Generation}}
\subsection{Commands for Mesh Generation}\label{sec:InitialMesh}
In Step1 in Section \ref{sec:example}, the keywords
\texttt{\bf border, buildmesh}
are explained.
%square,
%,  movemesh ,
% adaptmesh, readmesh, trunc, triangulate, splitmesh, emptymesh}
%
%The following keywords are discussed in this section:
%
%\texttt{\bf square, border, buildmesh,  movemesh ,
% adaptmesh, readmesh, trunc, triangulate, splitmesh, emptymesh}
%
% \index{square}\index{border}\index{buildmesh}\index{movemesh}\index{adaptmesh}
% \index{readmesh}\index{triangulate}\index{splitmesh}\index{emptymesh}

All the examples in this section come from the files \texttt{mesh.edp}
and \texttt{tablefunction.edp}.

\subsubsection{\setS{Square}}

For easy and simple testing, there is the command
``\texttt{\bf square}''.
The following
\bT
@mesh Th = @square(4,5);
\eT
generate a $4\times 5$ grid in the unit squre $[0,1]^2$ whose labels
are shown in Fig. \ref{fig:square}.
\begin{figure}[htbp]
\begin{center}
  \includegraphics[height=5cm]{figures/square}
\end{center}
  \caption{Boundary labels of the mesh by \texttt{square(10,10)}}
  \label{fig:square} \index{label}
\end{figure}
If you want constructs a
$n\times m$ grid in the rectangle $[x_0,x_1]\times [y_0,y_1]$, you can
write
\bFF
  @real x0=1.2,x1=1.8;
  @real y0=0,y1=1;
  @int n=5,m=20;
  @mesh Th=@square(n,m,[x0+(x1-x0)*x,y0+(y1-y0)*y]);
\eFF
 \begin{note}
You must notice that if you adding the name parameter \texttt{flags=1}, you  get a
Union Jack fags mesh pattern. \index{square!flags=}
\bFF
  @mesh Th=@square(n,m,[x0+(x1-x0)*x,y0+(y1-y0)*y],flags=1);
\eFF
\end{note}



\subsubsection{\setS{Border}}\index{border}\index{label}

A domain is defined as being on the left (resp right) of its
parameterized boundary
$$
\Gamma_j=\{(x,y)\left|\; x=\varphi_x(t),\, y=\varphi_y(t),\, a_j\le t\le b_j\right.\}
$$
We can easily check the orientation by drawing the curve
$t\mapsto (\varphi_x(t),\varphi_y(t)),\, t_0\le t\le t_1$.
If the figure become like to Fig. \ref{fig:border}, then
the domain lie on the shaded area, otherwise it lie on opposite side
(see also the examples enclosed with the box).
The boundaries $\Gamma_j$ can only intersect at their end points.

\begin{figure}[htbp]
\begin{center}
  \includegraphics[height=4cm]{figures/border}
\end{center}
  \caption{Orientation of the boundary defined by $(\phi_x(t),\phi_y(t))$}
  \label{fig:border} \index{border}
\end{figure}
The general expression of the triangulation is
\[
\ttCC{@mesh  ~~Mesh\_Name = @buildmesh$\left(\Gamma_1(m_1)+\cdots+\Gamma_J(m_j)\right)$;}
\]
where $m_j$ are numbers of marked points on $\Gamma_j,\,
\Gamma=\cup_{j=1}^J \Gamma_J$.
We can change the orientation of boundaries by changing the sign of $m_j$.
The following example shows how to change the orientation.
The example generates the unit disk
with a small circular hole, and assign ``1'' to the unit disk
(``2'' to the circle inside).
The boundary label must be non-zero, however we can omit the label
if we want use only the symbol.

\bFF
1: @border a(t=0,2*pi){ x=cos(t); y=sin(t);label=1;}
2: @border b(t=0,2*pi){ x=0.3+0.3*cos(t); y=0.3*sin(t);label=2;}
3: @plot(a(50)+b(+30)) ; // to see a plot of the border mesh \index{plot!border}
4: @mesh Thwithouthole= @buildmesh(a(50)+b(+30));
5: @mesh Thwithhole   = @buildmesh(a(50)+b(-30));
6: @plot(Thwithouthole,wait=1,ps="Thwithouthole.eps"); //figure \ref{Thwithouthole}\index{plot!mesh}
7: @plot(Thwithhole,wait=1,ps="Thwithhole.eps"); // figure \ref{Thwithhole}
\eFF
\begin{note}
You must notice that the orientation is changed by ``\texttt{b(-30)}'' in 5th line. In 7th line, \texttt{ps="fileName"} is used to generate a postscript file
identification that is shown on screen.
\end{note}

\twoplot[height=6cm]{Thwithouthole}{Thwithhole}{mesh without hole}{mesh with hole}

\subsubsection{Data Structure of Mesh and Reading/Writing a Mesh}

\index{readmesh}\index{savemesh}
Some user asked us that they want use the triangulation made from other tools
or hand-made mesh.
The example
\bT
@border C(t=0,2*pi) { x=cos(t); y=sin(t); }
@mesh Th = @buildmesh(C(10));
@savemesh("mesh_sample.msh");
\eT
make the mesh as in Fig. \ref{fig:meshSample}.
The informations about \texttt{Th} are save in the file ``mesh\_sample.msh''.

We can read from Fig. \ref{fig:meshSample} and ``mesh\_sample.msh''
as in Table \ref{tab:meshSample} where $n_v$ is
the number of vertices, $n_t$ number of triangles
and $n_s$ the number of edges on boundary.
For each vertex $q^i,\, i=1,\cdots,n_v$, we denote by $(q^i_x,q^i_y)$
the $x$-coordinate and $y$-coordinate.

Each triangle $T_k, k=1,\cdots,10$ have three vertices $q^{k_1},\, q^{k_2},\,q^{k_3}$ that are oriented in counterclockwise.
The boundary consists of 10 lines $L_i,\, i=1,\cdots,10$ whose tips are
$q^{i_1},\, q^{i_2}$.

\begin{figure}[htbp]
\begin{minipage}{\textwidth}
\begin{minipage}{0.5\textwidth}
\includegraphics[width=\textwidth]{figures/mesh_sample}%
\caption{mesh by \texttt{buildmesh(C(10))}}
\label{fig:meshSample}
\end{minipage}
\hspace{0.5mm}
\begin{minipage}{0.5\textwidth}
In the left figure, we have the following.\\\\
$n_v=14,\, n_t=16,\, n_s=10$\\\\
$q^1=(-0.309016994375,\, 0.951056516295)$\\
$\vdots\qquad \vdots\qquad \vdots$\\
$q^{14}=(-0.309016994375,\, -0.951056516295)$\\\\
The vertices of $T_1$ are $q^9,\, q^{12},\, q^{10}$.\\
$\vdots\qquad \vdots\qquad \vdots$\\
The vertices of $T_{16}$ are $q^9,\, q^{10},\, q^{6}$.\\\\
The edge of 1st side $L_1$ are $q^6,\, q^5$.\\
$\vdots\qquad \vdots\qquad \vdots$\\
The edge of 10th side $L_{10}$ are $q^{10},\, q^6$.\\
\end{minipage}
\end{minipage}
\end{figure}

\begin{table}[htbp]
\begin{tabular}{|l|l|}
\hline
Contents of file&Explanation\\
\hline
14 16 10& $n_v$\qquad  $n_t$\qquad $n_e$\\
-0.309016994375 0.951056516295 1& $q^1_x$\qquad $q^1_y$\qquad boundary label=1\\
0.309016994375 0.951056516295 1& $q^2_x$\qquad $q^2_y$\qquad boundary label=1\\
$\cdots$  $\cdots$ $\vdots$& \\
-0.309016994375 -0.951056516295 1& $q^{14}_x$\qquad $q^{14}_y$\qquad boundary label=1\\
\hline
9 12 10 0&$1_1$\quad $1_2$\quad $1_3$\quad region label=0 \\
5 9 6 0&$2_1$\quad $2_2$\quad $2_3$\quad region label=0  \\
$\cdots$& \\
9 10 6 0&$16_1$\quad $16_2$\quad $16_3$\quad region label=0 \\
\hline
6 5 1&$1_1\quad 1_2$\quad boundary label=1\\
5 2 1&$2_1\quad 2_2$\quad boundary label=1\\
$\cdots$& \\
10 6 1&$10_1\quad 10_2$\quad boundary label=1\\
\hline
\end{tabular}
  \caption{The structure of ``mesh\_sample.msh''}
  \label{tab:meshSample}
\end{table}

There are many mesh file formats available for communication with
other tools such as emc2, modulef.. (see \refSec{Mesh Files}),
The extension of a file gives the chosen
type.\index{bamg} More details can be found in the article by F. Hecht
"bamg : a bidimentional anisotropic mesh generator" available from the
FreeFem web site.  \\

The following give the example write and read files of generated mesh,
Freefem can read and write files. A
mesh file can be read back into \freefempp but the names of the
borders are lost. So these borders have to be referenced by the number
which corr