1%%%%%%%%%%%%%%%%%%% 2% XLiFE++ is an extended library of finite elements written in C++ 3% Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin 4% 5% This program is free software: you can redistribute it and/or modify 6% it under the terms of the GNU General Public License as published by 7% the Free Software Foundation, either version 3 of the License, or 8% (at your option) any later version. 9% This program is distributed in the hope that it will be useful, 10% but WITHOUT ANY WARRANTY; without even the implied warranty of 11% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12% GNU General Public License for more details. 13% You should have received a copy of the GNU General Public License 14% along with this program. If not, see <http://www.gnu.org/licenses/>. 15%%%%%%%%%%%%%%%%%%% 16 17\section{The variational approach} 18 19Before learning in details what \xlifepp is able to do, let us explain the basics with an example, the {\bf Helmholtz} equation: 20 21\begin{center} 22For a given function f(x,y), find a function u(x,y) satisfying 23\begin{equation} 24\label{eq.helmholtz} 25\begin{cases} 26\displaystyle -\Delta u(x,y) + u(x,y)=f(x,y) & \forall (x,y) \in \Omega \\ 27\displaystyle \frac{\partial u}{\partial n}(x,y)=0 & \forall (x,y) \in \partial\Omega 28\end{cases} 29\end{equation} 30\end{center} 31To solve this problem by a finite element method, \xlifepp is based on its variational formulation : find $u\in H^1(\Omega)$ such that $\forall v\in H^1(\Omega)$ 32 33\vspace{-0.5cm} 34\begin{center} 35\begin{equation} 36\label{fv.helmholtz} 37\displaystyle \int_\Omega\nabla u.\nabla v\,dx\,dy -\int_\Omega u\,v\,dx\,dy=\int_\Omega f\,v\,dx\,dy. 38\end{equation} 39\end{center} 40All the mathematical objects involved in the variational formulation are described in \xlifepp. The following program solves the Helmholtz problem with $f(x,y)=\cos{\pi x} \cos{\pi y}$ and $\Omega$ is the unit square. 41\vspace{.1cm} 42\inputCode[numbers=left, numberstyle=\tiny, numbersep=5pt, deletekeywords={[3]x,y}]{helmholtz.cpp} 43 44Please notice how close to the Mathematics, \xlifepp input language is. 45 46\section{How does it work ?} 47 48This first example shows how \xlifepp executes all the usual steps required by the {\bf Finite Element Method}. Let us walk through them one by one. 49 50\begin{description} 51\item[line 12 : ] every program using \xlifepp begins by a call to the \lstinline{init} function, taking up to 4 key/value arguments: 52\begin{description} 53\item[{\param \_lang}] enum to set the language for print and log messages. Possible values are \emph{en} for English, \emph{fr} for French, \emph{de} for German, or \emph{es} for Spanish. Default value is \emph{en}. 54\item[{\param \_verbose}] integer to set the verbose level. Default value is 1. 55\item[{\param \_trackingMode}] boolean to set if in the log file, you have a backtrace of every call to a \xlifepp routine. Default value is false. 56\item[{\param \_isLogged}] boolean to activate log. Default value is false. 57\end{description} 58Furthermore, the \lstinline{init} function loads functionalities linked to the trace of where such messages come from. If this function is not called, \xlifepp cannot work !!! 59 60\inputCode[firstline=12,lastline=12]{helmholtz.cpp} 61 62\item[lines 13-14 : ] The mesh will be generated on the unit square geometry with 11 nodes per edge. Arguments of a geometry are given with a key/value system. {\param \_origin} is the bottom left front vertex of {\class Square}. Next, we precise the mesh element type (here triangle), the mesh element order (here 1), and an optional description. See \autoref{ch.meshdef} for more examples of mesh definitions. 63 64\inputCode[firstline=13,lastline=14]{helmholtz.cpp} 65 66\item[line 15 : ] The main domain, named "Omega" in the mesh, is defined. 67\inputCode[firstline=15,lastline=15]{helmholtz.cpp} 68\item[line 16 : ] A finite element space is generally a space of polynomial functions on elements, triangles here only. Here \emph{sp} is defined as the space of continuous functions which are affine on each triangle $T_k$ of the domain $\Omega$, usually named $V_h$. The dimension of such a space is finite, so we can define a basis. 69$$ 70sp(\Omega,P1)=\left\{ w(x,y) \mbox{ such that } \exists (w_1,\ldots, w_N) \in \mathbb{R}^N, w(x,y) = \sum_{i=1}^N w_k \varphi_k(x,y)\right\} 71$$ 72 73where $N$ is the space dimension, i.e. the number of nodes, i.e. the number of vertices here. 74 75Currently, \xlifepp implements the following elements : $P_k$ on segment, triangle and tetrahedron, $Q^k$ on quadrangle and hexahedron, $O_k$ on prism and pyramid (see Mesh chapter for more details). 76 77\inputCode[firstline=16,lastline=16]{helmholtz.cpp} 78\item[lines 17-20 : ] The unknown $u$ here is an approximation of the solution of the problem. $v$ is declared as test function. This comes from the variational formulation of \autoref{eq.helmholtz} : multiplying both sides of equation and integrating over $\Omega$, we obtain : 79 80$$ 81-\int_{\Omega} v \Delta u dx dy +\int_{\Omega} v u dx dy = \int_{\Omega} v f dx dy 82$$ 83 84Then, using Green's formula, the problem is converted into finding $u$ such that : 85 86\begin{equation} 87\label{eq.helmholtz.vf} 88a(u,v) = \int_{\Omega} \nabla u \cdot \nabla v dx dy + \int_{\Omega} u v dx dy = \int_{\Omega} f v dx dy = l(v) 89\end{equation} 90 91The 4 next lines in the program declare $u$ and $v$ and define $a$ and $l$. 92 93\inputCode[firstline=17,lastline=20]{helmholtz.cpp} 94Please notice that: 95 96\begin{itemize} 97\item the test function is defined from the unknown. The reason is that the test function is dula to the unknown. Through the unknown, v is also defined on the same space. 98\item the right hand side needs the definition of the function $f$. Such function can be defined as a classical C++ function, but with a particular prototype. In this example, $f$ (i.e. $cosx2$) is a scalar function. So it takes 2 arguments : the first one is a {\class Point}, containing coordinates $x$ and $y$. The second one is optional and contains parameters to use inside the function. Here, the {\class Parameters} object is not used. At last, as a scalar function, it returns a {\class Real}. 99\inputCode[firstline=5,lastline=9, deletekeywords={[3]x,y}]{helmholtz.cpp} 100\end{itemize} 101\item[lines 21-22 : ] The previous definitions are a description of the variational form. Now, we have to define the matrix and the right-hand side vector which are the algebraic representations of the linear forms in the finite element space. This is done by the first 2 following lines. 102\inputCode[firstline=21,lastline=22]{helmholtz.cpp} 103\item[lines 23-24 : ] Matrix and vector being assembled, you can now choose the solver you want. Here, a conjugate gradient solver is used, with an initial guess constant equal to 1. 104 105\xlifepp offers you a various choice of direct or iterative solvers : 106\begin{itemize} 107\item $LU$, $LDU$, $LL^t$, $LDL^t$, $LDL^*$ factorizations 108\item BICG, BiCGStab, CG, CGS, GMRES, QMR, Sor, SSor, solvers 109\item internal eigen solver 110\item interfaces to external packages such as \umfpack, \arpack 111\end{itemize} 112 113See \autoref{ch.solvers} for more details. 114 115\inputCode[firstline=23,lastline=24]{helmholtz.cpp} 116\item[line 25 : ] To save the solution, \xlifepp provides an export to Paraview format file (vtu). 117\inputCode[firstline=25,lastline=25]{helmholtz.cpp} 118\item[line 26 : ] This is the end of the program. A "main" function always ends with this line. 119\inputCode[firstline=26,lastline=26]{helmholtz.cpp} 120\end{description} 121