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