manual-full.tex

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

%&LaTeX
\documentclass[a4paper,twoside,12pt]{book}
\usepackage{styles}
\def\key#1{\emph{#1}\index{#1}}
\def\vec#1{\mbox{\boldmath $#1$}}
\def\C{\mathbb{C}}
\def\Z{\mathbb{Z}}
\def\N{\mathbb{N}}

\def\setS#1{#1\label{sec:#1}}
\def\refSec#1{Section \ref{sec:#1}}
\def\refChap#1{Chapter \ref{sec:#1}}
\def\borderarray#1#2#3#4#5#6{%
\setbox0\hbox{$\begin{array}{#5}#6\end{array}$}
\setlength{\dimen1}{\wd0}\addtolength{\dimen1}{-#3}\addtolength{\dimen1}{-\arraycolsep}
\setlength{\dimen2}{\ht0}\addtolength{\dimen2}{-#4}
\setbox1\hbox{$\left#1\rule{\dimen1}{0pt}\rule{0pt}{\dimen2}\right#2$}
\setbox0\hbox{\raisebox{\dp0}{\box0}\kern-\dimen1\kern-5pt\raisebox{\dp1}{\box1}}
\vcenter{\box0}
}

\title{
%\DeclareFixedFont{\TitreFont}{\encodingdefault}{pnc}{r}{\shapedefault}{70pt}
 {\Blue{ \TitreFont FreeFem++ \\ \vglue 1cm  Manual}} \\ \vglue 3cm  ~ \\
      \normalsize  { version 1.46-1  \Red{(Under construction)} }
 \\ \vglue 0.7cm
 \large \url{http://www.freefem.org} \\
\url{http://www.freefem.org/ff++/index.htm}
}
\author{F. Hecht \thanks{\url{mailto:[e-mail]}},
 O. Pironneau \thanks{\url{mailto:[e-mail]}},
\\ Universit\'{e} Pierre et Marie Curie,
\\  Laboratoire Jacques-Louis Lions,
\\175 rue du Chevaleret ,PARIS XIII
\\ ~ \\  K. Ohtsuka  \thanks{\url{mailto:[e-mail]}}, \\
Hiroshima Kokusai Gakuin University, Hiroshima, Japan
}
\makeindex
\begin{document}
\graphicspath{{./}{plots/}{figures/}}
\ifpdf
\DeclareGraphicsExtensions{.pdf ,.jpg , .tif}
\else
\DeclareGraphicsExtensions{.eps ,.ps ,.jpg}
\fi

\maketitle

\tableofcontents
\let\subsubsection\subsection
\let\subsection\section
\let\section\chapter
\section{Introduction}
A partial differential equation is a relation between a function
of several variables and its (partial) derivatives.
 Many problems in physics, engineering, mathematics and even banking
are modeled by one or several partial differential equations.
\\\\

\texttt{Freefem} is a software to solve these equations numerically.
As its name says, it is a free software (see copyright for full detail)
based  on the Finite
Element Method. This software runs on all unix OS (with g++ 2.95.2 or
better and X11R6) , on Window95, 98, 2000, NT, XP,  on MacOS 9 and
X.  \\\\
Many phenomena involve several different fields.  Fluid-structure
interactions, Lorenz forces in liquid aluminium and ocean-atmosphere
problems are three such systems. These require approximations of
different degrees, possibly on different meshes.  Some algorithms such
as Schwarz' domain decomposition method also require data
interpolation on different meshes within one
program. \texttt{freefem++} can handle these difficulties, i.e. {\it
arbitrary finite element spaces on arbitrary unstructured and adapted
meshes}.  \\\\

The characteristics of \freefempp are:
\begin{itemize}
\item A large variety of finite elements : linear and quadratic
 Lagrangian elements, discontinuous P1 and Raviart-Thomas elements,
 elements of a non-scalar type, mini-element, ...
\item  Automatic interpolation of data on different meshes to an over mesh, store the interpolation matrix.
%
\item Linear problems description (real or complex) thanks to a formal variational form,
with access to the vectors and the matrix if needed.
%
\item Includes tools to define discontinuous Galerkin formulations
(please refer to the following keywords: ``jump'', ``average'',
``intalledges'').
%
\item Analytic description of boundaries. When two boundaries
intersect, the user must specify the intersection points.
%
\item Automatic mesh generator, based on the Delaunay-Vorono\"{i}
algorithm. Inner points density is proportional to the density of
points on the boundary \cite{George}.
%
\item Metric-based anisotropic mesh adaptation. The metric can be
computed automatically from the Hessien of a solution \cite{bamg}.
%
%\item Every variable is typed. For instance \texttt{f} is a function, specified
%by the keyword \texttt{func}.
%
\item Solvers : LU, Cholesky, Crout, CG, GMRES, UMFPACK linear solver,
eigenvalue and eigenvector computation.
%
\item Online graphics, C++-like syntax.
%
\item Many examples: Navier-Stokes, elasticity, Fluid structure,
Schwarz's domain decomposition method, Eigen value problem, residual
error indicator, ...
%
\item Parallel version using \texttt{mpi}
\end{itemize}


\subsection{History}

The project has evolved from \texttt{MacFem, PCfem}, written in
Pascal. The first C version was \texttt{freefem 3.4}; it offered mesh
adaptation on a single mesh.  A thorough rewriting in C++ led to
\texttt{freefem+ 1.2.10}, which also included interpolation over
multiple meshes (functions defined on one mesh can be used on any
other mesh).  Implementing the interpolation from one unstructured
mesh to another was not easy because it had to be fast and
non-diffusive; for each point, one had to find the containing
triangle. This is one of the basic problems of computational geometry
(see Preparata \& Shamos\cite{Preparata} for example). Doing it in a
minimum number of operations was a challenge.  Our implementation was
$O(n\log n)$ and based on a quadtree.  \\\\

We are now introducing \freefempp, an entirely new program written
in C++ and based on \texttt{bison} for a more adaptable
freefem language.
\\\\

The freefem language allows for a quick specification of any partial
differential system of equations.  The language syntax of \freefempp
is the result of a new design which makes use of the STL \cite{cpp},
templates and \texttt{bison} for its implementation,
 for more detail see \cite{FHcpp}.
  The outcome is a
versatile software in which any new finite element can be included in
a few hours; but a recompilation is then necessary.  The library of
finite elements available in \freefempp will therefore grow with the
version number.  So far we have linear and quadratic Lagrangian
elements, discontinuous P1 and Raviart-Thomas elements.

\subsection{Getting Started}
\label{sec:example}

We explain how \freefempp solve the problem \textbf{Poisson};
For a given function $f(x,y)$, find a function $u(x,y)$ satisfying
\begin{eqnarray}
\label{eqn:Poisson}
-\Delta u(x,y) &=& f(x,y)\quad \mbox{if $(x,y)$ is in }\Omega,
\qquad \Delta u = \partial^2 u/\partial x^2 + \partial^2 u/\partial y^2,\\
\label{eqn:Dirichlet}
u(x,y) &=& 0\quad \mbox{if $(x,y)$ is on }\partial\Omega
\end{eqnarray}
The following example shows \freefempp program solving $u$ when
$f(x,y)=xy$ (see 5th line) and $\Omega$ is the unit disk. The boundary $C=\partial\Omega$ is
$$
C=\{(x,y)|\; x=\cos(t),\, y=\sin(t),\, 0\le t\le 2\pi\}
\quad \textrm{(see 1st line)}
$$
The domain will
be on the left side of the oriented boundary by the parameter $t$.
As illustrated in Fig. \ref{firstU}, we can see the isovalue of $u$ by using \ttCC{@plot} (see 13th line).
\twoplot[height=5cm]{firstTh}{firstU}{mesh \ttCC{Th} by \ttCC{build(C(50))}}{isovalue by \ttCC{plot(u)}}

\begin{example}\label{exm:first}~
\bFF
 1: @border C(t=0,2*@pi){@x=cos(t); @y=sin(t);label=1;}
 2: @mesh Th = @buildmesh (C(50));
 3; @fespace Vh(Th,@P2);
 4: Vh u,v;
 5: @func f= x*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=0)  ;  // Dirichlet boundary condition
10:
11: @real cpu=clock();
12: Poisson; // SOLVE THE PDE
13: @plot(u);
14: @cout << " CPU time = " << clock()-cpu << @endl;
\eFF
\end{example}

\subsubsection{FEM by \freefempp}
The example shows \freefempp covers  easily all standard step in FEM (finite element method).
We explain how they are done by \freefempp in a step-by-step manner.\\

\textbf{Step1: Mesh Generation}\\

\textbf{1st line:} the boundary $\Gamma$ are described analytically (by opposition to CSG)
as stated before.
In the case $\Gamma=\sum_{j=0}^J \Gamma_j$ with curves $\Gamma_j$, then
the user must specify the
intersection points in case two boundaries intersect.
By the use of the keyword ``label'' such as
\bT
@border Gamma1(t=a1,b1) { x=$\cdots$; y=$\cdots$ ;label=1; }
 $\vdots$        $\vdots$           $\vdots$      $\vdots$
@border GammaJ(t=aJ,bJ) { x=$\cdots$; y=$\cdots$;label=1; }
\eT
the user can refer to $\Gamma$  by the number ``1''.
(examples are in \refSec{Meshing Examples}).

\textbf{2nd line:} the triangulation $\mathcal{T}_h$ of $\Omega$ is
automatically generated  by
``\ttCC{@buildmesh}(C(50))'' giving 50 points on $C$ as in Fig. \ref{firstTh}.
Automatic mesh generation is based on the Delaunay-Voronoi algorithm.
Refinement of the mesh are done by increasing the points on $\Gamma$,
for example, ``\ttCC{@buildmesh}(C(100))'', because inner vertices are
determined by the density of points on the boundary.

The symbol $\mathcal{T}_h$ (\texttt{Th} in \freefempp) shows
a family $\{T_k\}_{k=1,\cdots,n_t}$ of triangles in Fig. \ref{firstTh}
with the size $h$ of the mesh.
Here $n_t$ stands for the number of
triangules in $\mathcal{T}_h$.
If $\Omega$ is not polygonal domain, a ``skin'' remains between
the exact domain $\Omega$ and its approximation
$\Omega_h=\cup_{k=1}^{n_t}T_k$.
However, we notice that all corners of $\Gamma_h = \partial\Omega_h$ are
on $\Gamma$.\\

\textbf{Step2: Making finite element space}\\

\textbf{3rd line:} ``\ttCC{@fespace} Vh(Th,P2)'' makes the continuous Finite Element SPACE
\begin{equation}
V_h(\mathcal{T}_h,P_2)=\left\{w(x,y)\left|\;
w(x,y)=\sum_{k=0}^{M-1}w_k\phi_k(x,y),\, w_k\textrm{ are real numbers}\right.\right\}
\end{equation}
where $P_2$ indicate $\phi_k$ is a polynomial of degree $\le 2$, that is,
in each $T_k$,
$$
\phi_k(x,y)=\alpha_1+\alpha_2x+\alpha_3y+\alpha_4x^2+\alpha_5xy+\alpha_6y^2
$$
and the constants $\alpha_1,\,\cdots,\, \alpha_6$ are defined by its values
at the vertices of $T_k$ and their middle points that continuous in $\Omega$.
Here $w_k$ are called the degree of freedom of $w$ and $M$ the number of
the degree of freedom.
Already \freefempp implemented
P0, P1, P2, RT0, P1nc, P1dc, P2dc, P1b, P2b elements, .
The user can easily add a part of arbitrary degree elements to \freefempp,
so the available finite elements will differ with the version.\\

\textbf{Step3: Setting the problem}\\

\textbf{4th line:} ``\texttt{Vh u}'' declare that $u$ is approximated
through the use of the basis functions $\phi_k$ in $V_h$, that is,
\begin{equation*}
u(x,y)\simeq u_h(x,y)=\sum_{k=0}^{M-1} u_k\phi_k(x,y)
\end{equation*}
\textbf{5th line:} the given function $f$ is defined analytically using the keyword
\ttCC{@func}.

\textbf{6th--9th lines:} the formulation of (\ref{eqn:Poisson}) and
(\ref{eqn:Dirichlet}) are done.

Multiplying (\ref{eqn:Poisson}) by $v(x,y)$ and integrating
the result over $\Omega$, we have
$$
-\int_{\Omega}v\Delta u \,dxdy = \int_{\Omega} vf\, dxdy
$$
Then, by Green's formula, the problem Poisson is translated into finding $u$
such that
\begin{eqnarray}
\label{eqn:weakform}
&&a(u,v) - \ell(f,v) = 0
\qquad \textrm{for all }v\\
&&a(u,v)=\int_{\Omega}\nabla u\cdot \nabla v \,dxdy
\quad \ell(f,v)=\int_{\Omega}fv\, dxdy
\label{eqn:bilinear}
\end{eqnarray}
satisfying $v=0$ on $\partial\Omega$. In \freefempp
the problem \textbf{Poisson} is declared by
\begin{center}
\ttCC{Vh u,v; @problem Poisson(u,v) =}
\end{center}
and (\ref{eqn:weakform}) is expressed by symbols \ttCC{@dx}(u) $=\partial u/\partial x$, \ttCC{@dy}(u) $=\partial u/\partial y$ and
\begin{eqnarray*}
&&\int_{\Omega}\nabla u\cdot \nabla v\, dxdy \longrightarrow
\ttCC{@int@2d(Th)( @dx(u)*@dx(v) + @dy(u)*@dy(v) )}\\
&&\int_{\Omega}fv\, dxdy \longrightarrow
\ttCC{@int@2d(Th)( f*v )}\qquad
\textrm{(notice here, $u$ is unused)}\\
\end{eqnarray*}
In \freefempp the first and second formulas just above must be distinguished each other.
Because, the linear system to be solved are created from substituting $u_h$ for $u$ and $\phi_i$ for $v$ in (\ref{eqn:weakform}),
\begin{eqnarray}
\label{eqn:Equation}
A_{ij}u_j - F_i=0\quad i,j=0,\cdots,M-1;\qquad
F_i=\int_{\Omega}f\phi_i\, dxdy
\end{eqnarray}
and the solution $u_h=\sum_{j=0}^{M-1}u_j\phi_j$  must
satisfy ``$u_h=0$ on $\Gamma_h\simeq C$''.
The matrix $A=(A_{ij})$ is called \emph{stiffness matrix}
\index{matrix!stiffness matrix} and is modified from
\begin{eqnarray}
\label{eqn:Stiffness0}
A^0=\{A^0_{ij}\},\, A^0_{ij}=\int_{\Omega}\nabla \phi_j\cdot \nabla \phi_i \,dxdy\quad i.j=0,\cdots,M-1
\end{eqnarray}
to ensure $u=0$ on $C$ by ``+\ttCC{@on}\{1,u=0\}'' in 9th line.
If you want use the symbol ``C'' such as ``+\ttCC{@on}\{C,u=0\}'' (9th line),
then \emph{the user do not use the keyword ``label''}.

we can create directly the stiffness matrix $A$ in
(\ref{eqn:Equation}) by
\bT
@varf a(u,v) = @int2d(Th)( @dx(u)*@dx(v) + @dy(u)*@dy(v))
               + @on(C,u=0) ;
@matrix A=a(Vh,Vh);  // stiffness matrix, see (\ref{eqn:Stiffness0})
\eT
(see Example \ref{exm:second}).
The vector $B$ in (\ref{eqn:Equation}) is called \emph{load matrix}
\index{matrix!load matrix}, and also we get it:
\bT
@varf b([v],[f]) = @int2d(Th)(v*f);
@matrix B=b(Vh,Vh);
\eT

The linear system (\ref{eqn:Equation}) is solved by a Gauss LU
factorization. You can declare the solver of (\ref{eqn:Equation}),
for example,
\begin{center}
\ttCC{Vh u,v; @problem Poisson(u,v,@solver=CG) =}
\end{center}
means that (\ref{eqn:Equation}) will be solved by Conjugate Gradient method.\\

\textbf{Step4: Solving and visualization}\\

\textbf{11th line:} the current time is stored into the real-valued variable \ttCC{cpu}.

\textbf{12th line:} the problem is solved by calling its name.

\textbf{13th line:} the visualization is done as illustrated in Fig. \ref{firstU}
(see \refSec{Plot} for zoom, postscript and other commands).

\textbf{14th line:} the time in calculation is outputed into your console
(= default of standard output) using C++-like syntax.
The user need not study C++ for using \freefempp, because
C++-like syntax is used for input/output, loops, flow-controls etc.

\subsubsection{Features of \freefempp}