diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..db899980c4c3e089f1ef08478d1a079ef5bc1d9f
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,29 @@
+# list of ignored files/folders
+
+build/
+*~
+*.o
+*.obj
+*.bak
+.vscode/
+*.msh
+gmsh-sdk/
+eigen/
+*.tgz
+*.zip
+*.swp
+a.out
+
+# latex
+*.aux
+*.log
+*.nav
+*.out
+#*.pdf
+*.snm
+*.synctex*
+*.toc
+*.fls
+*.fdb_latexmk
+
+Backup_of*
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
new file mode 100644
index 0000000000000000000000000000000000000000..d0212f44044034adb745c31677465c14b844acf6
--- /dev/null
+++ b/.gitlab-ci.yml
@@ -0,0 +1,13 @@
+stages:
+ - build
+
+build_and_test:
+ image: rboman/waves
+ stage: build
+ script:
+ - source ./envs/linux-macos.sh
+ - cd lib && ./get_gmsh.sh && cd ..
+ - mkdir build
+ - cd build
+ - cmake ..
+ - make -j 4
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..90d81b21640dd301feb9b5b20f303fddf9c9b7c4
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,6 @@
+PROJECT(MATH0471 CXX)
+CMAKE_MINIMUM_REQUIRED(VERSION 3.12)
+
+ADD_SUBDIRECTORY( examples/gmsh_api )
+ADD_SUBDIRECTORY( examples/sparse_solve )
+ADD_SUBDIRECTORY( srcs )
diff --git a/doc/README.md b/doc/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..d1ccf35319f9a33ff93c6074ae9e905bf3afebb1
--- /dev/null
+++ b/doc/README.md
@@ -0,0 +1,3 @@
+# math0471/doc
+
+This folder contains the latest version of the presentation about thermoelasticity.
\ No newline at end of file
diff --git a/doc/bem/bem.cdr b/doc/bem/bem.cdr
new file mode 100644
index 0000000000000000000000000000000000000000..c2d0f1c1695211913860e67752883388f79dcbff
Binary files /dev/null and b/doc/bem/bem.cdr differ
diff --git a/doc/bem/bem.pdf b/doc/bem/bem.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..c8251a3b316f830c83a73ded2b959f49c0227f28
Binary files /dev/null and b/doc/bem/bem.pdf differ
diff --git a/doc/bem/bem.tex b/doc/bem/bem.tex
new file mode 100644
index 0000000000000000000000000000000000000000..9ba97c65614ab0294c7c175051bca2ef0e66fe7a
--- /dev/null
+++ b/doc/bem/bem.tex
@@ -0,0 +1,471 @@
+% !TeX spellcheck = en_GB
+
+\documentclass[10pt,xcolor=pdftex,dvipsnames,table]{beamer}
+
+\usetheme{Boadilla}
+
+\usepackage{lmodern}
+\usepackage[T1]{fontenc} % caracteres accentues
+\usepackage[utf8]{inputenc} % encodage unicode des sources .tex
+\usepackage{amsmath,amssymb,amsbsy,dsfont}
+
+% pour \mybox...
+\usepackage{xcolor}
+\setlength{\fboxrule}{0.1pt}
+\newcommand{\mybox}[2]{\fcolorbox{#1}{white}{$\displaystyle#2$}}
+
+
+\title{Boundary Element Method}
+\author{ MATH-0471 }
+%\institute{University of Li\`ege}
+\date{\today}
+
+%\AtBeginSection[]
+%{
+% \begin{frame}
+% \frametitle{Outline}
+% \tableofcontents[currentsection]
+% \end{frame}
+%}
+
+
+\begin{document}
+
+\frame{\titlepage}
+
+%\section*{Outline}
+
+%\begin{frame}
+% \frametitle{Outline}
+% \tableofcontents
+%\end{frame}
+
+% ---------------------------------------------------------------------------------------
+%\section{ FEM for 2D Linear Elasticity }
+% ---------------------------------------------------------------------------------------
+\begin{frame}
+ \frametitle{ Boundary Element Method }
+ \framesubtitle{ Potential problem }
+
+ \textbf{Potential problem} in $\Omega$:
+
+ \begin{equation*}
+ \left\{
+ \begin{aligned}
+ &\Delta u = 0 &&\text{ in } \Omega && \\
+ &u = \bar{u} && \text{ on } \Gamma_D && \text{ (Dirichlet)}\\
+ &\frac{\partial u}{\partial \boldsymbol{n}} = \bar{u}_n && \text{ on } \Gamma_N && \text{ (Neumann)} \\
+ \end{aligned}
+ \right.
+ \end{equation*}
+
+where
+\begin{itemize}
+ \item $u$ is the scalar unknown,
+ \item $\Gamma$ is the boundary of $\Omega$ and $\Gamma_D \cup \Gamma_N = \Gamma$ and $\Gamma_D \cap \Gamma_N = \emptyset$,
+ \item $\frac{\partial u}{\partial \boldsymbol{n}} = \boldsymbol{\nabla} u \cdot \boldsymbol{n}$, where $\boldsymbol{n}$ is the outward normal to $\Omega$.
+\end{itemize}
+ \vspace{\stretch{1}}
+
+Note that the PDE is homogeneous but the boundary conditions are not.
+ \vspace{\stretch{1}}
+
+ We follow here the notations from the book of \cite{Katsikadelis2016}, which is freely available from the \href{https://explore.lib.uliege.be/discovery/fulldisplay?docid=alma9919484598002321&context=L&vid=32ULG_INST:ULIEGE&lang=fr&search_scope=ULIEGE&adaptor=Local}{\color{blue}ULiege library website} (use VPN).
+
+\end{frame}
+
+\begin{frame}
+ \frametitle{ Boundary Element Method }
+ \framesubtitle{ Potential problem }
+
+ \centerline{\includegraphics[width=\textwidth]{fig_domain}}
+ \vspace{\stretch{1}}
+ The domain can be finite (left) or infinite (right).
+\end{frame}
+
+\begin{frame}
+ \frametitle{ Boundary Element Method }
+ \framesubtitle{ Fundamental solution / Free space Green's function in 2D}
+
+ If we ignore the boundary conditions, a \textbf{fundamental solution} $v$ of the Laplace problem defined in $\mathds{R}^2$ can be found by solving the PDE with a Dirac impulse as right-hand side.
+
+ \begin{equation*}
+ \Delta v = \delta(\boldsymbol{Q}-\boldsymbol{P})
+ \end{equation*}
+ where $\boldsymbol{P} = \boldsymbol{P}(x,y)$ is the position of the Dirac impulse and $\boldsymbol{Q} = \boldsymbol{Q}(\xi,\eta)$ is the position where $v$ is evaluated
+
+ \vspace{\stretch{1}}
+ We obtain
+ \begin{equation*}
+ v(\boldsymbol{Q}, \boldsymbol{P}) = \frac{1}{2\pi}\ln r
+ \end{equation*}
+
+ with $r = \sqrt{(\xi-x)^2+(\eta-y)^2} = || \boldsymbol{Q}-\boldsymbol{P} ||$
+
+ \vspace{\stretch{1}}
+
+ The fundamental solution is singular at $\boldsymbol{Q}=\boldsymbol{P}$.
+
+\end{frame}
+
+
+
+
+
+
+
+
+
+
+\begin{frame}
+ \frametitle{ Boundary Element Method }
+ \framesubtitle{ Representation formula in the domain $\Omega$ }
+
+ \textbf{Green's second identity} (also called Green's reciprocal identity):
+
+ \begin{equation*}
+ \int_{\Omega} \left( v\,\Delta u - u\,\Delta v\, \right) d\Omega
+ =
+ \int_{\Gamma} \left( v \frac{\partial u}{\partial \boldsymbol{n}} - u \frac{\partial v}{\partial \boldsymbol{n}} \right)\, ds
+ \end{equation*}
+
+
+
+ \vspace{\stretch{1}}
+
+ From this equality, if we choose $v$ as the fundamental solution introduced earlier, we can write a \textbf{representation formula} for the solution $u$ of the Laplace equation in $\Omega$ (boundary excluded):
+
+ \begin{equation*}
+ u(\boldsymbol{P}) = - \int_{\Gamma} \left( v(\boldsymbol{P},\boldsymbol{q}) \frac{\partial u(\boldsymbol{q})}{\partial \boldsymbol{n}_q} - u(\boldsymbol{q}) \frac{\partial v(\boldsymbol{P},\boldsymbol{q})}{\partial \boldsymbol{n}_q} \right)\, ds_q
+ \end{equation*}
+
+ \vspace{\stretch{1}}
+
+ Important remark concerning the notations:
+ \begin{itemize}
+ \item A capital letter $(\boldsymbol{P}, \boldsymbol{Q})$ represents a point inside the domain $\Omega$,
+ \item A lowercase letter $(\boldsymbol{p},\boldsymbol{q})$ represents a point on the boundary of $\Omega$.
+ \end{itemize}
+
+ \vspace{\stretch{1}}
+ This formula can be used to compute $u$ in $\Omega$ if we know both $u(\boldsymbol{q})$ and $\frac{\partial u(\boldsymbol{q})}{\partial \boldsymbol{n}_q}$ on the boundary (but we don't: we have either one or the other).
+
+\end{frame}
+
+
+\begin{frame}
+ \frametitle{ Boundary Element Method }
+ \framesubtitle{ Representation formula on the boundary $\Gamma$ }
+
+ From the Green's second identity, it is also possible to deduce a representation formula for the solution $u$ lying on $\Gamma$:
+
+ \begin{equation*}
+ \frac{\alpha}{2\pi} u(\boldsymbol{p}) = - \int_{\Gamma} \left( v(\boldsymbol{p},\boldsymbol{q}) \frac{\partial u(\boldsymbol{q})}{\partial \boldsymbol{n}_q} - u(\boldsymbol{q}) \frac{\partial v(\boldsymbol{p},\boldsymbol{q})}{\partial \boldsymbol{n}_q} \right)\, ds_q
+ \end{equation*}
+
+ where $\alpha$ is the angle of the corner made by the boundary at $\boldsymbol{p}$.
+
+ \vspace{\stretch{1}}
+ \centerline{\includegraphics[width=0.5\textwidth]{fig_alpha}}
+
+ If the boundary is smooth, $\alpha=\pi$ and the coefficient in front of $u(\boldsymbol{p})$ is $\frac{1}{2}$.
+
+ \vspace{\stretch{1}}
+ This \textbf{boundary integral equation} can be used to calculate $u$ and its derivative on the boundary.
+
+\end{frame}
+
+
+\begin{frame}
+ \frametitle{ Boundary Element Method }
+ \framesubtitle{ Principle of the BEM }
+
+ \begin{itemize}
+ \item The Boundary Element Method consists in meshing the boundary of the problem with a series of $N$ ``boundary elements''.
+ \item On each element, the solution $u$ is approximated by simple polynomials (constant, linear, etc.).
+ \item Injecting this approximation into the boundary integral equation and taking into account the boundary conditions allows us to calculate the solution and its derivative on the whole boundary.
+ \item Eventually, the solution inside the domain can be computed from the values on the boundary with the help of the representation formula.
+ \end{itemize}
+
+\end{frame}
+
+
+
+\begin{frame}
+ \frametitle{ Boundary Element Method }
+ \framesubtitle{ Constant boundary elements }
+
+ \begin{itemize}
+ \item The simplest way to approximate the solution is to use a piecewise-linear approximation of the boundary (a series of oriented straight segments).
+
+ \item On each element the solution $u$ and its derivative $u_n$ is considered constant (thus discontinuous across elements).
+
+ \item We also define a node at the centre of each segment.
+ \end{itemize}
+
+ \vspace{\stretch{1}}
+ \centerline{\includegraphics[width=0.4\textwidth]{fig_mesh}}
+
+\end{frame}
+
+
+
+\begin{frame}
+ \frametitle{ Boundary Element Method }
+ \framesubtitle{ Solving the boundary }
+
+ Discretized Boundary Integral Equation:
+ \begin{equation*}
+ \frac{1}{2} u_i = - \sum_{j=1}^{N} \int_{\Gamma_j} v(\boldsymbol{p}_i, \boldsymbol{q})
+ \underbrace{ \frac{\partial u(\boldsymbol{q})}{\partial \boldsymbol{n}_q} }_{u_n^j}
+ \,ds_q
+ + \sum_{j=1}^{N} \int_{\Gamma_j}
+ \underbrace{u(\boldsymbol{q})}_{u^j}
+ \frac{\partial v(\boldsymbol{p}_i,\boldsymbol{q})}{\partial \boldsymbol{n}_q}\,ds_q
+ \end{equation*}
+
+ Thanks to the node position centred in the middle of each segment, the boundary is locally smooth and $\alpha=\frac{1}{2}$.
+
+ \vspace{\stretch{1}}
+
+ Reordering the terms:
+ \begin{equation*}
+ -\frac{1}{2} u^i
+ + \sum_{j=1}^{N}
+ \underbrace{
+ \left[\int_{\Gamma_j} \frac{\partial v(\boldsymbol{p}_i,\boldsymbol{q})}{\partial \boldsymbol{n}_q}\,ds_q\right]
+ }_{\hat{H}_{ij}}
+ u^j
+ = \sum_{j=1}^{N}
+ \underbrace{
+ \left[ \int_{\Gamma_j} v(\boldsymbol{p}_i, \boldsymbol{q}) \,ds_q \right]
+ }_{G_{ij}}
+ u_n^j
+ \end{equation*}
+
+\end{frame}
+
+
+
+\begin{frame}
+ \frametitle{ Boundary Element Method }
+ \framesubtitle{ Solving the boundary }
+
+ \begin{equation*}
+ \sum_{j=1}^{N}
+ H_{ij}
+ u^j
+ = \sum_{j=1}^{N}
+ G_{ij}
+ u_n^j
+ \qquad
+ \text{ with }
+ \left\{
+ \begin{aligned}
+ &H_{ij} = \int_{\Gamma_j} \frac{\partial v(\boldsymbol{p}_i,\boldsymbol{q})}{\partial \boldsymbol{n}_q}\,ds_q -\frac{1}{2}\delta_{ij} \\
+ &G_{ij} = \int_{\Gamma_j} v(\boldsymbol{p}_i, \boldsymbol{q}) \,ds_q
+ \end{aligned}
+ \right.
+ \end{equation*}
+
+ %\vspace{\stretch{1}}
+ \centerline{\includegraphics[width=0.75\textwidth]{fig_integration}}
+
+ Difficulties are expected when $i=j$ (singular integrand).
+\end{frame}
+
+
+\begin{frame}
+ \frametitle{ Boundary Element Method }
+ \framesubtitle{ Solving the boundary }
+
+ In matrix form
+ \begin{equation*}
+ [H] \{u\} = [G] \{u_n\}
+ \end{equation*}
+
+ Equations can be sorted according to the type of boundary condition prescribed at the corresponding node (e.g. Dirichlet first, then Neumann):
+
+ \begin{equation*}
+ \begin{bmatrix}
+ \qquad\\
+ [H_{D}] & [H_N] \\
+ \qquad\\
+ \end{bmatrix}
+ \begin{Bmatrix}
+ \{u\}_D \\
+ \{u\}_N \\
+ \end{Bmatrix}
+ =
+ \begin{bmatrix}
+ \qquad\\
+ [G_{D}] & [G_N] \\
+ \qquad\\
+ \end{bmatrix}
+ \begin{Bmatrix}
+ \{u_n\}_D \\
+ \{u_n\}_N \\
+ \end{Bmatrix}
+ \end{equation*}
+ Since $\{u\}_D = \{\bar{u}\}$ and $\{u_n\}_N = \{\bar{u}_n\}$, we can gather the unknowns in the left-hand side:
+
+ \begin{equation*}
+ \underbrace{
+ \begin{bmatrix}
+ \qquad\\
+ [H_N] & -[G_D] \\
+ \qquad\\
+ \end{bmatrix}
+}_{[A]}
+ \underbrace{
+ \begin{Bmatrix}
+ \{u\}_N \\
+ \{u_n\}_D \\
+ \end{Bmatrix}
+}_{\{x\}}
+ =
+ \underbrace{
+ \begin{bmatrix}
+ \qquad\\
+ -[H_{D}] & [G_N] \\
+ \qquad\\
+ \end{bmatrix}
+ \begin{Bmatrix}
+ \{\bar{u}\} \\
+ \{\bar{u}_n\} \\
+ \end{Bmatrix}
+}_{ \{b\}}
+\end{equation*}
+This system is then solved using a linear solver.
+
+\end{frame}
+
+
+
+\begin{frame}
+ \frametitle{ Boundary Element Method }
+ \framesubtitle{ Calculation of $[G]$ }
+
+ $\bullet$ $G_{ij}$ with $i\neq j$ (Gauss quadrature):
+ \begin{equation*}
+ G_{ij} = \int_{\Gamma_j} v(\boldsymbol{p}_i, \boldsymbol{q}) \,ds_q = \int_{-1}^1 \frac{1}{2\pi} \ln (r(\xi))\,\frac{l_j}{2}\,d\xi
+ = \frac{l_j}{4\pi} \sum_{k=1}^{N^{GP}} \ln (r(\xi_k))\, w_k
+ \end{equation*}
+ where $l_j$ is the length of element $j$ and $N^{GP}$ is the number of Gauss points with are located at $\xi_k$ with weight $w_k$.
+
+ \vspace{\stretch{1}}
+
+ $\bullet$ $G_{ii}$ (analytical integration):
+ \begin{equation*}
+ G_{ii} = \frac{l_i}{2\pi} \left( \ln\frac{l_i}{2}-1\right)
+\end{equation*}
+
+\end{frame}
+
+
+\begin{frame}
+ \frametitle{ Boundary Element Method }
+ \framesubtitle{ Calculation of $[H]$ }
+
+ $\bullet$ $H_{ij}$ with $i\neq j$ (analytical integration):
+ \begin{equation*}
+ H_{ij} = \int_{\Gamma_j} \frac{\partial v(\boldsymbol{p}_i,\boldsymbol{q})}{\partial \boldsymbol{n}_q}\,ds_q
+ = \frac{a_{j+1}-a_j}{2\pi}
+ \end{equation*}
+
+ %\vspace{\stretch{1}}
+ \centerline{\includegraphics[width=0.45\textwidth]{fig_hij}}
+
+ \vspace{\stretch{1}}
+
+ $\bullet$ $H_{ii}$ (analytical integration):
+ \begin{equation*}
+ \hat{H}_{ii} = 0 \qquad \Rightarrow
+ H_{ii} = -\frac{1}{2}
+ \end{equation*}
+
+\end{frame}
+
+
+\begin{frame}
+ \frametitle{ Boundary Element Method }
+ \framesubtitle{ Calculation of the solution inside $\Omega$ }
+
+ From the representation formula:
+ \begin{equation*}
+ u(\boldsymbol{P}) = - \int_{\Gamma} \left( v(\boldsymbol{P},\boldsymbol{q}) \frac{\partial u(\boldsymbol{q})}{\partial \boldsymbol{n}_q} - u(\boldsymbol{q}) \frac{\partial v(\boldsymbol{P},\boldsymbol{q})}{\partial \boldsymbol{n}_q} \right)\, ds_q
+ \end{equation*}
+
+ we can write:
+ \begin{equation*}
+ u(\boldsymbol{P}) = -\sum_{j=1}^{N}
+ \left[ \int_{\Gamma_j} v(\boldsymbol{P}, \boldsymbol{q}) \,ds_q \right]
+ u_n^j
+ + \sum_{j=1}^{N}
+ \left[\int_{\Gamma_j} \frac{\partial v(\boldsymbol{P},\boldsymbol{q})}{\partial \boldsymbol{n}_q}\,ds_q\right]
+ u^j
+\end{equation*}
+
+\begin{itemize}
+ \item The solution inside $\Omega$ requires the computation of boundary integrals with regular integrands (the singularity is located at $\boldsymbol{P}$). Refer to the expressions of $H_{ij}$ and $G_{ij}$ for $i\neq j$.
+ \item If the solution at several points need to be computed, it can be easily performed in parallel.
+\end{itemize}
+
+
+\end{frame}
+
+
+
+\begin{frame}
+ \frametitle{ Boundary Element Method }
+ \framesubtitle{ Summary }
+
+ Advantages of the BEM:
+ \begin{itemize}
+ \item Only the surface of the problem should be meshed (the problem dimension is decreased by 1 compared to FEM).
+ \item The domain can be infinite.
+ \item The solution and its derivative can be easily evaluated anywhere in the domain.
+ \item Works well for stationary, linear and homogeneous PDEs with non-homogeneous boundary conditions.
+ \end{itemize}
+
+ \vspace{\stretch{1}}
+
+ Typical applications:
+ \begin{itemize}
+ \item Potential problems: stationary heat conduction, electrostatics, irrotational flows, Darcy flows,
+ \item Acoustics (Helmholtz equation),
+ \item Slow viscous flows (Stokes equations),
+ \item Elasticity (Lamé equations).
+ \end{itemize}
+
+\end{frame}
+
+
+
+\begin{frame}
+ \frametitle{ Boundary Element Method }
+ \framesubtitle{ Summary }
+
+ Issues of the BEM:
+ \begin{itemize}
+ \item The method requires an \textbf{integral representation} of the solution and a fundamental solution to be calculated, which is sometimes impossible (e.g. nonlinear equations).
+ \item Integrals of \textbf{singular functions} must be evaluated very precisely with special techniques for obtaining accurate results.
+ \item The method relies on \textbf{non-sparse, non-symmetric matrices} which lead large memory requirement and CPU time. Thus, complicated acceleration techniques (fast multipole method, hierarchical matrices) are needed for 3D problems.
+ \end{itemize}
+\end{frame}
+
+%------------------------------------------------
+
+\begin{frame}
+ \frametitle{References}
+ \footnotesize{
+ \begin{thebibliography}{99} % Beamer does not support BibTeX so references must be inserted manually as below
+ \bibitem[Katsikadelis, 2016]{Katsikadelis2016} John T. Katsikadelis (2016)
+ \newblock The Boundary Element Method for Engineers and Scientists -- Theory and Applications
+ \newblock Elsevier, \emph{Second edition}, ISBN: 9780128044933
+ \end{thebibliography}
+ }
+\end{frame}
+
+
+\end{document}
diff --git a/doc/bem/fig_alpha.pdf b/doc/bem/fig_alpha.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..f493a1fe2bf4a165259fe1144107b566a975243e
Binary files /dev/null and b/doc/bem/fig_alpha.pdf differ
diff --git a/doc/bem/fig_domain.pdf b/doc/bem/fig_domain.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..0fe59acfe29bfc7f6a455cf9fd935fd411ff088d
Binary files /dev/null and b/doc/bem/fig_domain.pdf differ
diff --git a/doc/bem/fig_hij.pdf b/doc/bem/fig_hij.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..1f2935ba1abade3e3854a8cf63ea8a5d2ba027ac
Binary files /dev/null and b/doc/bem/fig_hij.pdf differ
diff --git a/doc/bem/fig_integration.pdf b/doc/bem/fig_integration.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..39cbc4431419d48d3cef66348051ebfb5fac70e8
Binary files /dev/null and b/doc/bem/fig_integration.pdf differ
diff --git a/doc/bem/fig_mesh.pdf b/doc/bem/fig_mesh.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..d5375e1c218e6c5a2aa1f9d1ee0e42a916eae7a7
Binary files /dev/null and b/doc/bem/fig_mesh.pdf differ
diff --git a/doc/elasticity/elasticity.cdr b/doc/elasticity/elasticity.cdr
new file mode 100644
index 0000000000000000000000000000000000000000..b8ef3c0e9e2fbe6ddbea49726ab4d868a5fdce40
Binary files /dev/null and b/doc/elasticity/elasticity.cdr differ
diff --git a/doc/elasticity/elasticity.pdf b/doc/elasticity/elasticity.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..1892d920ed3ea1e632990f7f9740f6d6321f3638
Binary files /dev/null and b/doc/elasticity/elasticity.pdf differ
diff --git a/doc/elasticity/elasticity.tex b/doc/elasticity/elasticity.tex
new file mode 100644
index 0000000000000000000000000000000000000000..4a79b995b9e91963381f9684e5e42a94ebc5d7ea
--- /dev/null
+++ b/doc/elasticity/elasticity.tex
@@ -0,0 +1,625 @@
+% !TeX spellcheck = en_GB
+
+\documentclass[10pt,xcolor=pdftex,dvipsnames,table]{beamer}
+
+\usetheme{Boadilla}
+
+\usepackage{lmodern}
+\usepackage[T1]{fontenc} % caracteres accentues
+\usepackage[utf8]{inputenc} % encodage unicode des sources .tex
+\usepackage{amsmath,amssymb,amsbsy,dsfont}
+
+% pour \mybox...
+\usepackage{xcolor}
+\setlength{\fboxrule}{0.1pt}
+\newcommand{\mybox}[2]{\fcolorbox{#1}{white}{$\displaystyle#2$}}
+
+
+\title{FEM for 2D Linear Elasticity}
+\author{ MATH-0471 }
+%\institute{University of Li\`ege}
+\date{\today}
+
+%\AtBeginSection[]
+%{
+% \begin{frame}
+% \frametitle{Outline}
+% \tableofcontents[currentsection]
+% \end{frame}
+%}
+
+
+\begin{document}
+
+\frame{\titlepage}
+
+%\section*{Outline}
+
+%\begin{frame}
+% \frametitle{Outline}
+% \tableofcontents
+%\end{frame}
+
+% ---------------------------------------------------------------------------------------
+%\section{ FEM for 2D Linear Elasticity }
+% ---------------------------------------------------------------------------------------
+\begin{frame}
+ \frametitle{ FEM for 2D Linear Elasticity }
+ \framesubtitle{ Summary of equations }
+
+ \textbf{Mechanical equilibrium} in $\Omega$:
+
+ \begin{equation*}
+ \left\{
+ \begin{aligned}
+ \boldsymbol{\nabla}\cdot\boldsymbol{\sigma} + \rho\boldsymbol{b} &= 0
+ \qquad &&\text{ (translation) } \\
+ \boldsymbol{\sigma} &= \boldsymbol{\sigma}^T &&\text{ (rotation) }
+ \end{aligned}
+ \right.
+ \end{equation*}
+
+ where $\boldsymbol{\sigma}$ is the Cauchy stress tensor, $\rho$ is the density and $\boldsymbol{b}$ represents the mechanical volumic forces.
+
+ \vspace{\stretch{1}}
+
+
+and the \textbf{boundary conditions} on $\Gamma$, the boundary of $\Omega$:
+
+\begin{equation*}
+ \left\{
+ \begin{aligned}
+ & \boldsymbol{u} = \bar{\boldsymbol{u}} & \qquad \text { on } \Gamma_{\bar{\boldsymbol{u}}} & \qquad \text{ : prescribed displacement (Dirichlet)} \\
+ & \boldsymbol{\sigma}\cdot\boldsymbol{n} = \bar{\boldsymbol{t}} & \text { on } \Gamma_{\bar{\boldsymbol{t}}} & \qquad \text{ : prescribed surface traction (Neumann)} \\
+ \end{aligned}
+ \right.
+\end{equation*}
+where $\boldsymbol{u}$ is the displacement, $\boldsymbol{n}$ is the outward unit normal to the boundary, $\Gamma_{\bar{\boldsymbol{u}}} \cup \Gamma_{\bar{\boldsymbol{t}}} = \Gamma$ and $\Gamma_{\bar{\boldsymbol{u}}} \cap \Gamma_{\bar{\boldsymbol{t}}} = \emptyset$
+
+
+\end{frame}
+
+\begin{frame}
+ \frametitle{ FEM for 2D Linear Elasticity }
+ \framesubtitle{ Summary of equations }
+
+ \textbf{Strain-displacement relationship} (infinitesimal strains):
+ \begin{equation*}
+ \boldsymbol{\varepsilon}
+ = \tfrac{1}{2} \left( \boldsymbol{\nabla}\boldsymbol{u} + (\boldsymbol{\nabla}\boldsymbol{u})^T \right)
+ \end{equation*}
+ where small strains and small displacements are assumed.
+
+ \vspace{\stretch{1}}
+
+ \textbf{Constitutive law} (Hooke's law):
+ \begin{equation*}
+ \boldsymbol{\sigma} = \mathbb{H}:\boldsymbol{\varepsilon}
+ \end{equation*}
+ with the 4$^{th}$ order tensor
+ \begin{equation*}
+ \mathbb{H}_{ijkl} = \frac{E}{2(1+\nu)} \delta_{ik}\delta_{jl}
+ + \frac{\nu\,E}{(1+\nu)(1-2\nu)}\delta_{ij}\delta_{kl}
+ \quad\text{ with }
+ \left\{
+ \begin{aligned}
+ &E : \text{Young's modulus} \\
+ &\nu : \text{Poisson's ratio} \\
+ \end{aligned}
+ \right.
+ \end{equation*}
+ for an isotropic material.
+\end{frame}
+
+
+
+
+\begin{frame}
+ \frametitle{ FEM for 2D Linear Elasticity }
+ \framesubtitle{ Weak formulation }
+
+ A \textbf{weak formulation} is obtained by multiplying the PDE by a function $\boldsymbol{w}(\boldsymbol{x})$ and by integrating over $\Omega$:
+
+ \begin{equation*}
+ \int_{\Omega} \boldsymbol{w}\cdot \left(\boldsymbol{\nabla}\cdot\boldsymbol{\sigma}\right)\, dV
+ + \int_{\Omega} \rho \boldsymbol{w}\cdot \boldsymbol{b}\, dV = 0
+ \end{equation*}
+
+ \vspace{\stretch{1}}
+
+ The first term of the left-hand side can be integrated by parts:
+ \begin{equation*}
+ \int_{\Omega} \boldsymbol{w}\cdot \left(\boldsymbol{\nabla}\cdot\boldsymbol{\sigma}\right)\, dV
+ = \int_{\Gamma_{\bar{\boldsymbol{u}}}} \boldsymbol{w}\cdot \boldsymbol{\sigma} \cdot \boldsymbol{n}\, dS
+ + \int_{\Gamma_{\bar{\boldsymbol{t}}}} \boldsymbol{w}\cdot \underbrace{\boldsymbol{\sigma} \cdot \boldsymbol{n}}_{\boldsymbol{\bar{t}}}\, dS
+ - \int_{\Omega} \boldsymbol{\nabla}\boldsymbol{w} : \boldsymbol{\sigma}\, dV
+ \end{equation*}
+
+ The integral on $\Gamma_{\bar{u}}$ can disappear if we choose test functions $\boldsymbol{w}(\boldsymbol{x})$ that vanish on this boundary.
+\end{frame}
+
+
+\begin{frame}
+ \frametitle{ FEM for 2D Linear Elasticity }
+ \framesubtitle{ Weak formulation }
+
+ Moreover, since $\boldsymbol{\sigma} = \boldsymbol{\sigma}^T$, we can write
+
+ \begin{equation*}
+ \boldsymbol{\nabla}\boldsymbol{w} : \boldsymbol{\sigma} = \underbrace{\tfrac{1}{2} \left( \boldsymbol{\nabla}\boldsymbol{w} + (\boldsymbol{\nabla}\boldsymbol{w})^T \right)}_{\text{similar to } \boldsymbol{\varepsilon}(\boldsymbol{u})} : \boldsymbol{\sigma}
+ \end{equation*}
+
+
+ \vspace{\stretch{1}}
+
+ The weak formulation becomes\footnote{This is also called ``the principle of virtual work''. \\
+ Notations in \cite{ponthot2020}: $\boldsymbol{w} \leftrightarrow \boldsymbol{\delta u}$ and $\boldsymbol{\delta \varepsilon} \leftrightarrow \tfrac{1}{2} \left( \boldsymbol{\nabla}\boldsymbol{w} + (\boldsymbol{\nabla}\boldsymbol{w})^T \right)$) \\
+ and $\boldsymbol{u}+\boldsymbol{\delta u}=\boldsymbol{\bar{u}}$ on $\Gamma_{\bar{u}}$ (``kinematically admissible virtual displacement'').\\}:
+
+ \begin{block}{Weak formulation}
+ Find $\boldsymbol{u}$ with $\boldsymbol{u} = \bar{\boldsymbol{u}}$ on $\Gamma_{\bar{u}}$ such that
+ \begin{equation*}
+ \int_{\Gamma_{\bar{\boldsymbol{t}}}} \boldsymbol{w}\cdot \boldsymbol{\bar{t}}\, dS
+ - \int_{\Omega} \tfrac{1}{2} \left( \boldsymbol{\nabla}\boldsymbol{w} + (\boldsymbol{\nabla}\boldsymbol{w})^T \right) : \boldsymbol{\sigma}\, dV
+ + \int_{\Omega} \rho\boldsymbol{w} \cdot \boldsymbol{b}\, dV = 0
+ \end{equation*}
+ for all test functions $\boldsymbol{w}$ which are 0 on $\Gamma_{\bar{u}}$.
+ \end{block}
+
+\end{frame}
+
+
+\begin{frame}
+ \frametitle{ FEM for 2D Linear Elasticity }
+ \framesubtitle{ Spatial discretization - Finite Element approximation }
+
+ The discretization procedure is similar to the one followed for the Laplace/Poisson equation in MATH0024, except that the unknown field $\boldsymbol{u}(\boldsymbol{x}, t)$ is now a vector instead of a scalar. The equations are detailed here for the 2D case.
+ \vspace{\stretch{1}}
+
+ Using Cartesian coordinates, each component of $\boldsymbol{u}(\boldsymbol{x}, t)$ and $\boldsymbol{w}(\boldsymbol{x})$ are approximated using well-chosen shape functions $N^k$:
+ \begin{equation*}
+ \left\{
+ \begin{aligned}
+ &u_1(\boldsymbol{x}, t) = \sum_k N^k(\boldsymbol{x})\,u_1^k(t) \\
+ &u_2(\boldsymbol{x}, t) = \sum_k N^k(\boldsymbol{x})\,u_2^k(t) \\
+ \end{aligned}
+ \right.
+ \quad
+ \text{ and }
+ \left\{
+ \begin{aligned}
+ &w_1(\boldsymbol{x}) = \sum_k N^k(\boldsymbol{x})\,w_1^k \\
+ &w_2(\boldsymbol{x}) = \sum_k N^k(\boldsymbol{x})\,w_2^k \\
+ \end{aligned}
+ \right.
+ \end{equation*}
+ where $u_i^k$ is the $i^{\text{th}}$ coordinate of node $k$ (also called ``degree of freedom'').
+\end{frame}
+
+
+\begin{frame}
+ \frametitle{ FEM for 2D Linear Elasticity }
+ \framesubtitle{ Spatial discretization - Finite Element approximation }
+
+ In matrix form (assuming a mesh of $n$ nodes in a 2D space):
+ \begin{equation*}
+ \boldsymbol{u}(\boldsymbol{x},t) = \mathbf{N}(\boldsymbol{x})\,\boldsymbol{d}(t)
+ \end{equation*}
+ with
+ \begin{equation*}
+ \mathbf{N}(\boldsymbol{x}) =
+ \begin{bmatrix}
+ N_1 & 0 & N_2 & 0 & \dots & N_n & 0 \\
+ 0 & N_1 & 0 & N_2 & \dots & 0 & N_n \\
+ \end{bmatrix}
+ \end{equation*}
+ \begin{equation*}
+ \boldsymbol{d}^T =
+ \begin{bmatrix}
+ u_1^1 & u_2^1 & u_1^2 & u_2^2 & \dots & u_1^n & u_2^n
+ \end{bmatrix}
+ \end{equation*}
+
+ \vspace{\stretch{1}}
+
+ Similarly, for the test functions:
+ \begin{equation*}
+ \boldsymbol{w}(\boldsymbol{x}) = \mathbf{N}(\boldsymbol{x})\,\boldsymbol{h}
+ \end{equation*}
+\end{frame}
+
+
+
+\begin{frame}
+ \frametitle{ FEM for 2D Linear Elasticity }
+ \framesubtitle{ Isoparametric Finite Elements }
+
+
+ \centerline{\includegraphics[width=0.7\textwidth]{fig_isoparametric}}
+
+ \vspace{\stretch{1}}
+
+ \textbf{Isoparametric finite elements} ($m$ nodes) are a classical choice. Their geometry is interpolated with the same shape functions\footnote{$N^k$: superscript $k$ is related nodes where $N^k=1$}
+ as the unknown field $\boldsymbol{u}(\boldsymbol{x}, t)$ and the test functions $\boldsymbol{w}(\boldsymbol{x})$.
+ \begin{equation*}
+ \boldsymbol{x}(\boldsymbol{\xi}) = \sum_{k=1}^m N^k(\boldsymbol{\xi})\,\boldsymbol{x}^k
+ \end{equation*}
+ where $\boldsymbol{\xi}$ are reduced coordinates (usually $\xi_i\in[-1,1]$)
+\end{frame}
+
+
+
+\begin{frame}
+ \frametitle{ FEM for 2D Linear Elasticity }
+ \framesubtitle{ 3D Hooke's law (matrix form) }
+
+ \textbf{Voigt's notation}: in a practical implementation, $4^{\text{th}}$ and $2^{\text{nd}}$-order tensors are replaced by matrices and vectors respectively.
+ \begin{equation*}
+ \begin{bmatrix}
+ \sigma_{11} \\
+ \sigma_{22} \\
+ \sigma_{33} \\
+ \sigma_{12} \\
+ \sigma_{13} \\
+ \sigma_{23} \\
+ \end{bmatrix}
+ = \frac{E}{(1+\nu)(1-2\nu)}
+ \begin{bmatrix}
+ 1-\nu & \nu & \nu & 0 & 0 & 0 \\
+ \nu & 1-\nu & \nu & 0 & 0 & 0 \\
+ \nu & \nu & 1-\nu & 0 & 0 & 0 \\
+ 0 &0 & 0 & \frac{1-2\nu}{2} & 0 & 0 \\
+ 0 &0 & 0 & 0 &\frac{1-2\nu}{2} & 0 \\
+ 0 &0 & 0 & 0 & 0 & \frac{1-2\nu}{2} \\
+ \end{bmatrix}
+ \begin{bmatrix}
+ \varepsilon_{11} \\
+ \varepsilon_{22} \\
+ \varepsilon_{33} \\
+ 2\varepsilon_{12} \\
+ 2\varepsilon_{13} \\
+ 2\varepsilon_{23} \\
+ \end{bmatrix}
+ \end{equation*}
+ such that
+ \begin{equation*}
+ \boldsymbol{\varepsilon}:\boldsymbol{\sigma} = \boldsymbol{\varepsilon}^v\cdot \boldsymbol{\sigma}^v
+ \end{equation*}
+
+ Hooke's law becomes:
+ \begin{equation*}
+ \boldsymbol{\sigma}^v = \mathbf{H}\,\boldsymbol{\varepsilon}^v
+ \end{equation*}
+
+ \vspace{\stretch{1}}
+
+ Remark: $\boldsymbol{\varepsilon}^v$ involves ${\color{red}2\times}\varepsilon_{ij}$ = $\gamma_{ij}$ (shear strain/angle) for $i\neq j$.
+
+\end{frame}
+
+
+
+\begin{frame}
+ \frametitle{ FEM for 2D Linear Elasticity }
+ \framesubtitle{ 2D Hooke's law (matrix form) }
+
+ A very common 2D assumption is the \textbf{plane-stress hypothesis} which is valid for thin structures ($\sigma_{33}=\sigma_{13}=\sigma_{23}=0$):
+
+ \begin{equation*}
+ \begin{bmatrix}
+ \varepsilon_{11} \\
+ \varepsilon_{22} \\
+ \varepsilon_{33} \\
+ 2\varepsilon_{12} \\
+ 2\varepsilon_{13} \\
+ 2\varepsilon_{23} \\
+ \end{bmatrix}
+ = \underbrace{\frac{1}{E}
+ \begin{bmatrix}
+ 1 & -\nu & -\nu & 0 & 0 & 0 \\
+ -\nu & 1 & -\nu & 0 & 0 & 0 \\
+ -\nu & -\nu & 1 & 0 & 0 & 0 \\
+ 0 &0 & 0 & \frac{1+\nu}{2} & 0 & 0 \\
+ 0 &0 & 0 & 0 &\frac{1+\nu}{2} & 0 \\
+ 0 &0 & 0 & 0 & 0 & \frac{1+\nu}{2} \\
+ \end{bmatrix}}_{\boldsymbol{H}^{-1}}
+ \begin{bmatrix}
+ \sigma_{11} \\
+ \sigma_{22} \\
+ \color{red}0 \\
+ \sigma_{12} \\
+ \color{red}0 \\
+ \color{red}0 \\
+ \end{bmatrix}
+ \end{equation*}
+
+ \vspace{\stretch{1}}
+
+ The 2D Hooke's law is obtained by inverting this relationship:
+ \begin{equation*}
+ \begin{bmatrix}
+ \sigma_{11} \\
+ \sigma_{22} \\
+ \sigma_{12} \\
+ \end{bmatrix}
+ = \underbrace{\frac{E}{1-\nu^2}
+ \begin{bmatrix}
+ 1 & \nu & 0 \\
+ \nu & 1 & 0 \\
+ 0 & 0 & \frac{1-\nu}{2} \\
+ \end{bmatrix} }_{\boldsymbol{H}_{2D}}
+ \begin{bmatrix}
+ \varepsilon_{11} \\
+ \varepsilon_{22} \\
+ 2\varepsilon_{12} \\
+ \end{bmatrix}
+ \end{equation*}
+ This $3\times 3$ matrix is used as $\boldsymbol{H}$ in the following equations.
+
+\end{frame}
+
+
+\begin{frame}
+ \frametitle{ FEM for 2D Linear Elasticity }
+ \framesubtitle{ Strain-displacement relationship }
+
+ Using Voigt's notation, the strain-displacement relationship becomes:
+
+ \begin{equation*}
+ \begin{bmatrix}
+ \varepsilon_{11} \\
+ \varepsilon_{22} \\
+ 2\varepsilon_{12} \\
+ \end{bmatrix}
+ =
+ \begin{bmatrix}
+ \frac{\partial}{\partial x_1} & 0 \\
+ 0 &\frac{\partial}{\partial x_2} \\
+ \frac{\partial}{\partial x_2} & \frac{\partial}{\partial x_1} \\
+ \end{bmatrix}
+ \begin{bmatrix}
+ u_1 \\
+ u_2 \\
+ \end{bmatrix}
+ \end{equation*}
+
+ or, symbolically
+ \begin{equation*}
+ \boldsymbol{\varepsilon}^v = \boldsymbol{\partial}\,\boldsymbol{u}
+ \end{equation*}
+
+ \vspace{\stretch{1}}
+
+ Remarks:
+ \begin{itemize}
+ \item $\varepsilon_{33}$ does not appear in the 2D equations. However, it can be computed, if needed, from the stresses:
+ \begin{equation*}
+ \varepsilon_{33}=-\frac{\nu}{E}\left(\sigma_{11}+\sigma_{22}\right)
+ \end{equation*}
+ \item Strains (and stresses) are discontinuous across finite elements boundaries unless $\boldsymbol{u}(\boldsymbol{x})$ is $C^1$, which is almost never the case in practise.
+ \end{itemize}
+\end{frame}
+
+
+
+\begin{frame}
+ \frametitle{ FEM for 2D Linear Elasticity }
+ \framesubtitle{ Strains/stresses as a function of the unknowns}
+
+ Computation of \textbf{strains} from the nodal displacements $\boldsymbol{d}$:
+ \begin{equation*}
+ \boldsymbol{\varepsilon}^v = \boldsymbol{\partial}\,\boldsymbol{u} = \boldsymbol{\partial}\,\mathbf{N}\,\boldsymbol{d} = \mathbf{B}\,\boldsymbol{d}
+ \end{equation*}
+ with
+ \begin{equation*}
+ \mathbf{B} =
+ \begin{bmatrix}
+ \frac{\partial}{\partial x_1} & 0 \\
+ 0 &\frac{\partial}{\partial x_2} \\
+ \frac{\partial}{\partial x_2} & \frac{\partial}{\partial x_1} \\
+ \end{bmatrix}
+ \begin{bmatrix}
+ N^1 & 0 & \dots & N^m & 0 \\
+ 0 & N^1 & \dots & 0 & N^m \\
+ \end{bmatrix}
+ \end{equation*}
+
+
+ \vspace{\stretch{1}}
+
+ The \textbf{stresses} can also be obtained from the nodal displacement $\boldsymbol{d}$:
+
+ \begin{equation*}
+ \boldsymbol{\sigma}^v = \mathbf{H}\,\mathbf{B}\,\boldsymbol{d}
+ \end{equation*}
+
+\end{frame}
+
+
+\begin{frame}
+ \frametitle{ FEM for 2D Linear Elasticity }
+ \framesubtitle{ Discretized weak formulation }
+
+ Back to the \textbf{weak formulation}
+ \begin{equation*}
+ \int_{\Gamma_{\bar{\boldsymbol{t}}}} \boldsymbol{w}\cdot \boldsymbol{\bar{t}}\, dS
+ - \int_{\Omega}
+ \underbrace{\frac{1}{2} \left( \boldsymbol{\nabla}\boldsymbol{w} + (\boldsymbol{\nabla}\boldsymbol{w})^T \right)
+ : \boldsymbol{\sigma}}_{(\mathbf{B}\boldsymbol{h})\cdot\boldsymbol{\sigma}^v }\, dV
+ + \int_{\Omega} \rho \boldsymbol{w} \cdot \boldsymbol{b}\, dV = 0
+ \end{equation*}
+
+
+ \vspace{\stretch{1}}
+
+ Replacing $\boldsymbol{w}=\mathbf{N}\boldsymbol{h}$ and $\boldsymbol{u}=\mathbf{N}\boldsymbol{d}$ and $\boldsymbol{\sigma}^v = \mathbf{H}\,\mathbf{B}\,\boldsymbol{d}$
+
+
+ \begin{equation*}
+ \boldsymbol{h}^T \left( \int_{\Gamma_{\bar{\boldsymbol{t}}}} \mathbf{N}^T\boldsymbol{\bar{t}}\, dS \right) - \boldsymbol{h}^T \left(\int_{\Omega}\mathbf{B}^T\mathbf{H}\mathbf{B} \, dV \right) \boldsymbol{d}
+ + \boldsymbol{h}^T \left( \int_{\Omega} \rho\mathbf{N}^T\boldsymbol{b}\, dV \right) = 0
+ \end{equation*}
+
+ This relationship should be satisfied for any $\boldsymbol{h}$.
+
+\end{frame}
+
+
+
+
+\begin{frame}
+ \frametitle{ FEM for 2D Linear Elasticity }
+ \framesubtitle{ Discretized FE equations }
+
+ This leads to the system of discretized equations:
+ \begin{equation*}
+ \underbrace{\left(\int_{\Omega}\mathbf{B}^T\mathbf{H}\mathbf{B} \, dV \right)}_{\mathbf{K}} \boldsymbol{d}
+ =
+ \underbrace{\left( \int_{\Gamma_{\bar{\boldsymbol{t}}}} \mathbf{N}^T\boldsymbol{\bar{t}}\, dS \right)
+ + \left( \int_{\Omega} \rho\,\mathbf{N}^T\boldsymbol{b}\, dV \right)}_{\boldsymbol{f}}
+ \end{equation*}
+
+ \vspace{\stretch{1}}
+
+ \begin{block}{Discretized equations}
+ \begin{equation*}
+ \mathbf{K}\boldsymbol{d}=\boldsymbol{f}
+ \end{equation*}
+ \end{block}
+
+
+\end{frame}
+
+\begin{frame}
+ \frametitle{ FEM for 2D Linear Elasticity }
+ \framesubtitle{ Finite Elements }
+
+ \textbf{Assembly}:
+ \begin{itemize}
+ \item The integrals involved in the calculation of $\mathbf{K}$ and $\boldsymbol{f}$ are expressed as a sum of integrals over each finite element. These elemental contributions are then assembled (summed) in a large and structural matrix and vector.
+ \item $\mathbf{K}$ is a symmetric sparse matrix and this property should be exploited to reduce the required storage and to efficiently solve the linear system of equations.
+ \end{itemize}
+
+ \vspace{\stretch{1}}
+
+ \textbf{Boundary conditions}:
+ \begin{itemize}
+ \item The equations of the system corresponding to nodes where Dirichlet boundary conditions are prescribed should be discarded and replaced by equations enforcing these conditions ($\boldsymbol{u}=\bar{\boldsymbol{u}}$) at these nodes.
+ \item Homogeneous Dirichlet boundary conditions ($\boldsymbol{u}=0$) lead to remove the lines and columns of the matrices corresponding to the components of the displacement where these conditions are prescribed.
+ \end{itemize}
+\end{frame}
+
+
+\begin{frame}
+ \frametitle{ FEM for 2D Linear Elasticity }
+ \framesubtitle{ Practical calculation of the elemental stiffness matrices }
+
+ \begin{equation*}
+ K_{ij} = \int_{\Omega} B_{ki}(\boldsymbol{x})\,H_{kl}\,B_{lj}(\boldsymbol{x}) \, dV
+ \end{equation*}
+ \vspace{\stretch{1}}
+
+ Isoparametric finite elements:
+ \begin{equation*}
+ K_{ij} = \int_{-1}^1\int_{-1}^1\int_{-1}^1 B_{ki}(\boldsymbol{\xi})\,H_{kl}\,B_{lj}(\boldsymbol{\xi})\, \underbrace{\det (\frac{\partial{\boldsymbol{x}} }{\partial\boldsymbol{\xi}})}_{\text{jacobian}} \, d\boldsymbol{\xi}
+ \end{equation*}
+
+ \vspace{\stretch{1}}
+
+ \centerline{\includegraphics[width=0.8\textwidth]{fig_isoparametric2}}
+
+\end{frame}
+
+
+\begin{frame}
+ \frametitle{ FEM for 2D Linear Elasticity }
+ \framesubtitle{ Gauss quadrature }
+
+ This integral is calculated using a Gauss quadrature: ($N^{\text{GP}}$ chosen $\rightarrow$ $\boldsymbol{\xi}^p$, $w_p$)
+
+ \begin{equation*}
+ K_{ij} \approx \sum_{p=1}^{N^{\text{GP}}} \underbrace{
+ B_{ki}(\boldsymbol{\xi})\,H_{kl}\,B_{lj}(\boldsymbol{\xi})
+ \det \mathbf{J}
+ }_{\text{all the factors evaluated at } \boldsymbol{\xi}=\boldsymbol{\xi}^p} \,w_p
+ \end{equation*}
+ with $N^{\text{GP}}$, the number of Gauss points, $\boldsymbol{\xi}^p$ the positions and $w_p$, the weights.
+
+ \vspace{\stretch{1}}
+ \centerline{\includegraphics[width=0.4\textwidth]{fig_gauss}}
+
+\end{frame}
+
+
+\begin{frame}
+ \frametitle{ FEM for 2D Linear Elasticity }
+ \framesubtitle{ Jacobian of isoparametric elements at each Gauss points}
+
+ \textbf{Jacobian matrix}\footnote{transposed in \cite{ponthot2020}}: let $\boldsymbol{x}=(x_1, x_2, x_3)$ and $\boldsymbol{\xi}=(\xi_1, \xi_2, \xi_3)$
+ \begin{equation*}
+ J_{ij} (\boldsymbol{\xi}) =\frac{\partial x_i}{\partial \xi_j} = \frac{\partial N^k (\boldsymbol{\xi})}{\partial \xi_j} x_i^k
+ \end{equation*}
+ where $x_i^k$ is the $i^{\text{th}}$ coordinate of the $k^{\text{th}}$ node of the finite element.
+
+ \vspace{\stretch{1}}
+ The Jacobian matrix and its determinant must be evaluated at each Gauss point ($\boldsymbol{\xi}=\boldsymbol{\xi}^p$) of each finite element.
+
+ \vspace{\stretch{1}}
+
+ Note: the derivatives of the shape functions evaluated at the Gauss points $\frac{\partial N^k}{\partial \xi_j} (\boldsymbol{\xi}^p)$ are the same for all elements and can be computed once.
+
+\end{frame}
+
+
+
+
+\begin{frame}
+ \frametitle{ FEM for 2D Linear Elasticity }
+ \framesubtitle{ Derivatives of the shape functions }
+
+ \begin{equation*}
+ \mathbf{B} =
+ \begin{bmatrix}
+ \frac{\partial}{\partial x_1} & 0 \\
+ 0 &\frac{\partial}{\partial x_2} \\
+ \frac{\partial}{\partial x_2} & \frac{\partial}{\partial x_1} \\
+ \end{bmatrix}
+ \begin{bmatrix}
+ N^1 & 0 & \dots & N^m & 0 \\
+ 0 & N^1 & \dots & 0 & N^m \\
+ \end{bmatrix}
+ \end{equation*}
+
+ The matrix $\mathbf{B}$ contains the derivatives of the shape functions with respect to $\boldsymbol{x}$ which should be transformed into derivatives of the shape functions with respect to $\boldsymbol{\xi}$ using the Jacobian matrix:
+
+ \begin{equation*}
+ \frac{\partial}{\partial\xi_i} = \left( \frac{\partial x_j}{\partial \xi_i} \right) \frac{\partial}{\partial x_j} = J_{ji} \frac{\partial}{\partial x_j} \qquad \Rightarrow \boldsymbol{\nabla}_{\color{red}\boldsymbol{x}} = \mathbf{J}^{-T} \boldsymbol{\nabla}_{\color{red}\boldsymbol{\xi}}
+ \end{equation*}
+
+ \vspace{\stretch{1}}
+
+ Thus, the derivative of $N^i$ with respect to $\boldsymbol{x}$ can be computed by this formula:
+
+ \begin{equation*}
+ \boldsymbol{\nabla}_{\boldsymbol{x}} N^i = \mathbf{J}^{-T} \boldsymbol{\nabla}_{\boldsymbol{\xi}} N^i
+ \end{equation*}
+
+\end{frame}
+
+
+
+%------------------------------------------------
+
+\begin{frame}
+ \frametitle{References}
+ \footnotesize{
+ \begin{thebibliography}{99} % Beamer does not support BibTeX so references must be inserted manually as below
+ \bibitem[Ponthot, 2020]{ponthot2020} J.-P. Ponthot (2020)
+ \newblock An Introduction to the Finite Element Method
+ \newblock \emph{Lecture notes}, chapters 10 -- 11.
+ \end{thebibliography}
+ }
+\end{frame}
+
+
+\end{document}
diff --git a/doc/elasticity/fig_gauss.pdf b/doc/elasticity/fig_gauss.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..c8b0c11008c4258deaba565da6267664a502b6b3
Binary files /dev/null and b/doc/elasticity/fig_gauss.pdf differ
diff --git a/doc/elasticity/fig_isoparametric.pdf b/doc/elasticity/fig_isoparametric.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..d901cf183acf0eee3f9676e67afcc835d8f104a2
Binary files /dev/null and b/doc/elasticity/fig_isoparametric.pdf differ
diff --git a/doc/elasticity/fig_isoparametric2.pdf b/doc/elasticity/fig_isoparametric2.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..5386a4424cfb7dd2c4c6b89864d1c51bc620ef1f
Binary files /dev/null and b/doc/elasticity/fig_isoparametric2.pdf differ
diff --git a/envs/README.md b/envs/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..fbcc703be3add5fc2189f220f7045ddb46b4249e
--- /dev/null
+++ b/envs/README.md
@@ -0,0 +1,14 @@
+# math0471/envs
+
+This folder contains 2 scripts for the correct definition of environment variables of your system so that
+* Your c++ compiler (g++) is available from the command line,
+* Eigen and gmsh libraries and header files can be found by cmake,
+* The gmsh executable can be used in the command line (by just typing `gmsh`).
+* `cmake` and `make` can be used on Windows as if you were on linux/mac (see `cmake.cmd` and `make.cmd` for more information)
+
+Use `windows.cmd` on Windows (run it or double-click on it)
+
+Use `linux-macos.sh` on linux or macOS with the command:
+```
+source linux-macos.sh
+```
diff --git a/envs/cmake.cmd b/envs/cmake.cmd
new file mode 100644
index 0000000000000000000000000000000000000000..662f965de4ed89c8dc9e7a0c9e1d1b411cb4e274
--- /dev/null
+++ b/envs/cmake.cmd
@@ -0,0 +1,2 @@
+@echo off
+cmake.exe -G"MinGW Makefiles" %*
\ No newline at end of file
diff --git a/envs/compiler b/envs/compiler
new file mode 100644
index 0000000000000000000000000000000000000000..820ba5b04e70545315afb8673c465dfeda418112
--- /dev/null
+++ b/envs/compiler
@@ -0,0 +1 @@
+ = NOT found!
diff --git a/envs/example_of_launch_mac.json b/envs/example_of_launch_mac.json
new file mode 100644
index 0000000000000000000000000000000000000000..e4079da034c6ccbdff1f5a9d505ca5df90233b00
--- /dev/null
+++ b/envs/example_of_launch_mac.json
@@ -0,0 +1,26 @@
+{
+"version": "0.2.0",
+"configurations": [
+ {
+ "name": "clang++ - Build and debug active file",
+ "type": "cppdbg",
+ "request": "launch",
+ "program": "${workspaceFolder}/build/bin/advection",
+ "args": [],
+ "stopAtEntry": true,
+ "cwd": "${workspaceFolder}",
+ "environment": [
+ {
+ "name": "PATH",
+ "value": "${env:PATH}:${workspaceFolder}/gmsh-sdk/bin:${workspaceFolder}/gmsh-sdk/lib"
+ },
+ {
+ "name": "DYLD_LIBRARY_PATH",
+ "value": "${env:DYLD_LIBRARY_PATH}:${workspaceFolder}/gmsh-sdk/bin:${workspaceFolder}/gmsh-sdk/lib"
+ }
+ ],
+ "externalConsole": false,
+ "MIMode": "lldb"
+ }
+ ]
+}
\ No newline at end of file
diff --git a/envs/example_of_launch_win.json b/envs/example_of_launch_win.json
new file mode 100644
index 0000000000000000000000000000000000000000..55b3cb34cedd671c557462bd8cea714a66950c1d
--- /dev/null
+++ b/envs/example_of_launch_win.json
@@ -0,0 +1,33 @@
+{
+ // Use IntelliSense to learn about possible attributes.
+ // Hover to view descriptions of existing attributes.
+ // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387
+ "version": "0.2.0",
+ "configurations": [
+ {
+ "name": "(gdb) Launch edges.exe",
+ "type": "cppdbg",
+ "request": "launch",
+ "program": "${workspaceFolder}/build/bin/edges.exe",
+ "args": [],
+ "stopAtEntry": true,
+ "cwd": "${workspaceFolder}",
+ "environment": [
+ {
+ "name": "PATH",
+ "value": "${env:PATH};C:\\mingw-w64\\mingw64\\bin;${workspaceFolder}\\gmsh-sdk\\bin;${workspaceFolder}\\gmsh-sdk\\lib"
+ }
+ ],
+ "externalConsole": false,
+ "MIMode": "gdb",
+ "miDebuggerPath": "C:/mingw-w64/mingw64/bin/gdb.exe",
+ "setupCommands": [
+ {
+ "description": "Enable pretty-printing for gdb",
+ "text": "-enable-pretty-printing",
+ "ignoreFailures": true
+ }
+ ]
+ }
+ ]
+}
\ No newline at end of file
diff --git a/envs/linux-macos.sh b/envs/linux-macos.sh
new file mode 100644
index 0000000000000000000000000000000000000000..9970673d782322100e11d4a020740a14a792be0c
--- /dev/null
+++ b/envs/linux-macos.sh
@@ -0,0 +1,32 @@
+# This file setup an environment which allows you to compile the code
+# on Linux or macOS using the default compiler (gcc or clang)
+#
+# HOW TO USE THIS FILE?
+# open a terminal
+# source ./envs/linux-macos.sh
+# mkdir build
+# cd build
+# cmake ..
+# make
+# [executables are built in the bin/ folder]
+# ctest
+
+# try to get the absolute path of root folder
+DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd && echo x)"; DIR="${DIR%x}"
+abspath=$(dirname $DIR )
+#echo "abspath=${abspath}"
+
+# set the location of gmsh SDK ( **MODIFY THIS LINE FOR YOUR SYSTEM** )
+GMSHSDK=${abspath}/lib/gmsh-sdk
+EIGEN=${abspath}/lib/eigen
+
+# where are gmsh and gmsh-**.so ?
+export PATH=${GMSHSDK}/bin:${GMSHSDK}/lib:${PATH}
+# where is gmsh.h ?
+export INCLUDE=${EIGEN}:${GMSHSDK}/include:${INCLUDE}
+# where is gmsh.lib ?
+export LIB=${GMSHSDK}/lib:${LIB}
+# where is gmsh.py ? (required only if you want to use the python API)
+export PYTHONPATH=${GMSHSDK}/lib:${PYTHONPATH}
+# the following command is only useful for macOS
+export DYLD_LIBRARY_PATH=${GMSHSDK}/lib:${DYLD_LIBRARY_PATH}
diff --git a/envs/make.cmd b/envs/make.cmd
new file mode 100644
index 0000000000000000000000000000000000000000..cd6247a535b99591f68b00b81cc7d587239cb7ab
--- /dev/null
+++ b/envs/make.cmd
@@ -0,0 +1,2 @@
+@echo off
+mingw32-make.exe %*
\ No newline at end of file
diff --git a/envs/windows.cmd b/envs/windows.cmd
new file mode 100644
index 0000000000000000000000000000000000000000..f3de67ad128b07f806c1a8bcb10d319faa00d5fa
--- /dev/null
+++ b/envs/windows.cmd
@@ -0,0 +1,82 @@
+@echo off
+:: This file opens a terminal which allows you to compile the code with
+:: a 64-bit mingw compiler
+:: https://mingw-w64.org/
+:: https://www.msys2.org/
+::
+:: HOW TO USE THIS FILE?
+:: [check the PATHs below]
+:: [run this file]
+:: mkdir build
+:: cd build
+:: cmake ..
+:: make
+:: ctest
+::
+:: How to clean the "build" folder using cmd line?
+:: cd build
+:: rd /q /s .
+
+echo setting MinGW64 environment...
+
+:: set the location of gmsh SDK / mingw compiler
+set GMSHSDK=%~dp0..\lib\gmsh-sdk
+set EIGEN=%~dp0..\lib\eigen
+set COMPILERPATHS=C:\msys64\mingw64\bin;C:\mingw-w64\mingw64\bin
+
+:: perform some tests...
+IF NOT EXIST %GMSHSDK% (
+ ECHO - gmsh-sdk NOT found in %GMSHSDK%!!
+ @REM PAUSE
+ @REM EXIT /B
+) ELSE (
+ ECHO - gmsh-sdk found in %GMSHSDK%.
+)
+
+IF NOT EXIST %EIGEN% (
+ ECHO - eigen NOT found in %EIGEN%!!
+ @REM PAUSE
+ @REM EXIT /B
+) ELSE (
+ ECHO - eigen found in %EIGEN%.
+)
+
+FOR %%a IN (%COMPILERPATHS%) DO (
+ IF EXIST %%a (
+ ECHO - compiler found in %%a.
+ SET COMPILERPATH=%%a
+ ) ELSE (
+ IF NOT DEFINED COMPILERPATH (
+ ECHO - compiler not found in %%a.
+ )
+ )
+)
+IF NOT DEFINED COMPILERPATH (
+ ECHO => compiler NOT found!
+ PAUSE
+ EXIT /B
+)
+
+:: where is gmsh.exe and gmsh-**.dll ?
+set PATH=%GMSHSDK%\bin;%GMSHSDK%\lib;%PATH%
+:: where is gmsh.h ?
+set INCLUDE=%EIGEN%;%GMSHSDK%\include;%INCLUDE%
+:: where is gmsh.lib ?
+set LIB=%GMSHSDK%\lib;%LIB%
+:: where is gmsh.py ? (required only if you want to use the python API)
+set PYTHONPATH=%GMSHSDK%\lib;%PYTHONPATH%
+:: add path to the compiler to the PATH
+set PATH=%COMPILERPATH%;%PATH%
+:: add current folder to PATH for cmake/make aliases
+set PATH=%~dp0;%PATH%
+
+:: clear local vars
+set GMSHSDK=
+set EIGEN=
+set COMPILERPATHS=
+set COMPILERPATH=
+
+:: open terminal
+CD /d "%~dp0"
+CD ..
+%comspec% /K
diff --git a/examples/README.md b/examples/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..c59b02e0e589b89e2bee185063832e7eab654225
--- /dev/null
+++ b/examples/README.md
@@ -0,0 +1,3 @@
+# math0471/examples
+
+This folder contains code examples that are useful for the project.
diff --git a/examples/gmsh_api/CMakeLists.txt b/examples/gmsh_api/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..7b47598dc342cdb1e17dd1cdc74c6802cc36c7cb
--- /dev/null
+++ b/examples/gmsh_api/CMakeLists.txt
@@ -0,0 +1,56 @@
+PROJECT(MATH0471 CXX)
+CMAKE_MINIMUM_REQUIRED(VERSION 3.12)
+
+# ------------------------------------------------------------------------------
+# Find libraries and setup compiler
+# ------------------------------------------------------------------------------
+
+# build type is "" by default in Linux
+IF(NOT CMAKE_BUILD_TYPE)
+ SET( CMAKE_BUILD_TYPE "Release" CACHE STRING "" FORCE)
+ENDIF()
+
+# enable C++11
+SET(CMAKE_CXX_STANDARD 11)
+SET(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+IF(APPLE)
+ # on macOS, do not give priority to frameworks/apps
+ SET(CMAKE_FIND_APPBUNDLE LAST)
+ SET(CMAKE_FIND_FRAMEWORK LAST)
+ENDIF()
+
+# find gmsh-sdk
+# gmsh.h
+FIND_PATH(GMSH_INCLUDE_DIRS NAMES "gmsh.h")
+MESSAGE(STATUS "GMSH_INCLUDE_DIRS=" ${GMSH_INCLUDE_DIRS})
+if(NOT GMSH_INCLUDE_DIRS)
+ MESSAGE(FATAL_ERROR "gmsh.h not found!")
+ENDIF()
+INCLUDE_DIRECTORIES(${GMSH_INCLUDE_DIRS})
+
+# libgmsh.so
+FIND_LIBRARY(GMSH_LIBRARIES gmsh)
+MESSAGE(STATUS "GMSH_LIBRARIES=" ${GMSH_LIBRARIES})
+IF(NOT GMSH_LIBRARIES)
+ MESSAGE(FATAL_ERROR "gmsh library not found!")
+ENDIF()
+
+# find Eigen
+find_path(EIGEN_INCLUDE_DIRS "Eigen/Dense"
+ PATHS "${PROJECT_SOURCE_DIR}/lib/eigen" "/usr/include/eigen3")
+MESSAGE(STATUS "EIGEN_INCLUDE_DIRS=" ${EIGEN_INCLUDE_DIRS})
+IF(NOT EIGEN_INCLUDE_DIRS)
+ MESSAGE(FATAL_ERROR "Eigen include dir not found!")
+ENDIF()
+INCLUDE_DIRECTORIES(${EIGEN_INCLUDE_DIRS})
+
+# ------------------------------------------------------------------------------
+
+# compile each .cpp into an exe
+FILE(GLOB SRCS *.cpp)
+FOREACH(SFILE ${SRCS})
+ GET_FILENAME_COMPONENT(PROG_NAME ${SFILE} NAME_WE)
+ ADD_EXECUTABLE(${PROG_NAME} ${SFILE})
+ TARGET_LINK_LIBRARIES(${PROG_NAME} ${GMSH_LIBRARIES})
+ENDFOREACH()
diff --git a/examples/gmsh_api/README.md b/examples/gmsh_api/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..e1283e8ad6d2f103f973c3c8f3ffc6ae65f34988
--- /dev/null
+++ b/examples/gmsh_api/README.md
@@ -0,0 +1,9 @@
+# gmsh_api
+
+This folder contains programs related to the presentation named "First use of Gmsh and Gmsh SDK".
+
+You can copy any example/tutorial from gmsh into this folder so that it gets included in the build process when you run cmake and make.
+
+Examples can be found in the following folders:
+* `gmsh-sdk/share/doc/gmsh/demos/api/`
+* `gmsh-sdk/share/doc/gmsh/tutorial/c++/`
\ No newline at end of file
diff --git a/examples/gmsh_api/code1_loadgeo.cpp b/examples/gmsh_api/code1_loadgeo.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..adbac2853fcd37b735bed8e5a9261c1117e88df8
--- /dev/null
+++ b/examples/gmsh_api/code1_loadgeo.cpp
@@ -0,0 +1,36 @@
+// This program takes a .geo file as argument, generates the mesh of the
+// surfaces and displays the result in the Gmsh window.
+//
+// This example introduces:
+// - the gmsh.h header
+// - initialize()/finalize() functions
+// - C++ namespaces / standard library (std::cout)
+//
+// How to build/run? (windows)
+// [in examples/gmsh_api]
+// mkdir build
+// cd build
+// cmake .. && make && code1_loadgeo.exe ..\rectangle.geo
+
+#include <gmsh.h> // gmsh main header (READ THIS FILE)
+#include <iostream> // for std::cout
+
+int main(int argc, char **argv)
+{
+ // the program requires 1 argument: a .geo file.
+ if (argc < 2)
+ {
+ std::cout << "Usage: " << argv[0] << " <geo_file>\n";
+ return 0;
+ }
+
+ gmsh::initialize(); // first thing to do before using any gmsh SDK function
+
+ gmsh::open(argv[1]); // similar to "File / Open" in the GUI
+ gmsh::model::mesh::generate(2); // similar to "Mesh / 2D" in the GUI
+
+ gmsh::fltk::run(); // opens the gmsh window
+
+ gmsh::finalize(); // clean-up things before ending the program
+ return 0;
+}
diff --git a/examples/gmsh_api/code2_getgroups.cpp b/examples/gmsh_api/code2_getgroups.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a917e1aaa72dd9a16a407abd5ffb79a614ece299
--- /dev/null
+++ b/examples/gmsh_api/code2_getgroups.cpp
@@ -0,0 +1,57 @@
+// This program loops over the physical groups defined in a .geo file
+// and prints their names, dimensions and tags.
+//
+// This example introduces:
+// - (dim,tags)
+// - std::pair<T1,T2>
+// - std::vector<T>
+// - gmsh::vectorpair == std::vector<std::pair<int,int>>
+// - references / const references
+// - how to call a gmsh function (input vs output arguments)
+// - std::string
+//
+// How to build/run? (windows)
+// [in examples/gmsh_api]
+// mkdir build
+// cd build
+// cmake .. && make && code2_getgroups.exe ..\rectangle.geo
+
+#include <gmsh.h>
+#include <iostream>
+
+int main(int argc, char **argv)
+{
+ if (argc < 2)
+ {
+ std::cout << "Usage: " << argv[0] << " <geo_file>\n";
+ return 0;
+ }
+
+ gmsh::initialize();
+
+ gmsh::open(argv[1]);
+ gmsh::model::mesh::generate(2);
+
+ // -- code added to code1.cpp (begin) --------------------------------------
+
+ gmsh::vectorpair dimTags; // output parameter of "getPhysicalGroups", defining an empty vector of pairs
+ gmsh::model::getPhysicalGroups(dimTags); // vector is filled
+
+ for(int i=0; i<dimTags.size(); ++i) // size gives size of the vector
+ {
+ int dim = dimTags[i].first; // dimTags[i] is a std::pair<int,int>
+ int tag = dimTags[i].second;
+
+ std::string name; // output parameter of "getPhysicalName"
+ gmsh::model::getPhysicalName(dim, tag, name);
+
+ // prints the data to the console
+ std::cout << "Physical group (" << dim << "D) "
+ << "named '"<< name << "' tag=" << tag << '\n';
+ }
+
+ // -- code added to code1.cpp (end) ----------------------------------------
+
+ gmsh::finalize();
+ return 0;
+}
diff --git a/examples/gmsh_api/code3_buildphymap.cpp b/examples/gmsh_api/code3_buildphymap.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f6671dad41e428b32a2fe467d3e764a085b7884e
--- /dev/null
+++ b/examples/gmsh_api/code3_buildphymap.cpp
@@ -0,0 +1,69 @@
+// This program loops over the physical groups defined in a .geo file
+// and prints their names, dimensions and tags.
+//
+// This example introduces:
+// - std::map
+// - Physical Group
+// => Entities
+// => Element types
+// => Elements/nodes hierarchy
+//
+// How to build/run? (windows)
+// [in example/gmsh_api]
+// mkdir build
+// cd build
+// cmake .. && make && code3_buildphymap.exe ..\rectangle.geo
+
+#include <gmsh.h>
+#include <iostream>
+#include <map>
+
+int main(int argc, char **argv)
+{
+ if (argc < 2)
+ {
+ std::cout << "Usage: " << argv[0] << " <geo_file>\n";
+ return 0;
+ }
+
+ gmsh::initialize();
+
+ gmsh::open(argv[1]);
+ gmsh::model::mesh::generate(2);
+
+ // -- code added to previous code (begin) ----------------------------------
+ // definition of the map: "group_name" => (dim,tag)
+ std::map< std::string, std::pair<int,int> > groups;
+ // -- code added to previous code (end) ------------------------------------
+
+ gmsh::vectorpair dimTags;
+ gmsh::model::getPhysicalGroups(dimTags);
+
+ for(int i=0; i<dimTags.size(); ++i)
+ {
+ int dim = dimTags[i].first;
+ int tag = dimTags[i].second;
+
+ std::string name;
+ gmsh::model::getPhysicalName(dim, tag, name);
+
+ std::cout << "Physical group (" << dim << "D) "
+ << "named '"<< name << "' tag=" << tag << '\n';
+
+ // -- code added to previous code (begin) ------------------------------
+ // set value (dim, tag) for key 'name'
+ groups[name] = {dim, tag};
+ // -- code added to previous code (end) --------------------------------
+ }
+
+ // -- code added to previous code (begin) ----------------------------------
+ // print map
+ std::cout << "Physical Groups:\n";
+ for (auto &x : groups)
+ std::cout << "\t '" << x.first << "' => (" << x.second.first << ','
+ << x.second.second << ")\n";
+ // -- code added to previous code (end) ------------------------------------
+
+ gmsh::finalize();
+ return 0;
+}
diff --git a/examples/gmsh_api/code4_hierarchy.cpp b/examples/gmsh_api/code4_hierarchy.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5247ad261e6cedb042d0973045f8e38d1f32030b
--- /dev/null
+++ b/examples/gmsh_api/code4_hierarchy.cpp
@@ -0,0 +1,94 @@
+// This program loops over the physical groups defined in a .geo file
+// and prints their names, dimensions and tags.
+//
+// This example loops over the gmsh data structure hierarchy:
+// - Physical Group
+// => contains 1 or more "Entities"
+// => contains 1 or more "Element types"
+// => access to Elements
+//
+// How to build/run? (windows)
+// [in example/gmsh_api]
+// mkdir build
+// cd build
+// cmake .. && make && code4_hierarchy.exe ..\rectangle.geo
+
+#include <gmsh.h>
+#include <iostream>
+#include <map>
+
+int main(int argc, char **argv)
+{
+ if (argc < 2)
+ {
+ std::cout << "Usage: " << argv[0] << " <geo_file>\n";
+ return 0;
+ }
+
+ gmsh::initialize();
+
+ gmsh::open(argv[1]);
+ gmsh::model::mesh::generate(2);
+
+ std::map<std::string, std::pair<int, int>> groups;
+ gmsh::vectorpair dimTags;
+ gmsh::model::getPhysicalGroups(dimTags);
+ for (int i = 0; i < dimTags.size(); ++i)
+ {
+ int dim = dimTags[i].first;
+ int tag = dimTags[i].second;
+ std::string name;
+ gmsh::model::getPhysicalName(dim, tag, name);
+ groups[name] = {dim, tag};
+ }
+
+ // -- code added to previous code (begin) ----------------------------------
+
+ std::string groupname = "left_edge";
+
+ int dim = groups[groupname].first;
+ int tag = groups[groupname].second;
+ if (tag == 0)
+ {
+ std::cerr << "Group '" << groupname << "' does not exist!\n";
+ return 1;
+ }
+
+ // get entities of the chosen physical group
+ std::vector<int> tags;
+ gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
+
+ std::cout << "Entities in group named '" << groupname << "': ";
+ for (int i = 0; i < tags.size(); ++i)
+ std::cout << tags[i] << " ";
+ std::cout << '\n';
+
+ // loop over entities
+ for (int k = 0; k < tags.size(); ++k)
+ {
+ std::cout << "Entity (" << dim << "," << tags[k] << "):\n";
+ std::vector<int> elementTypes;
+ std::vector<std::vector<std::size_t>> elementTags;
+ std::vector<std::vector<std::size_t>> nodeTags;
+ gmsh::model::mesh::getElements(elementTypes, elementTags,
+ nodeTags, dim, tags[k]);
+
+ // loop over element types
+ for (int i = 0; i < elementTypes.size(); ++i)
+ {
+ std::cout << "\telement type " << elementTypes[i] << ":\n";
+ std::cout << "\t\telements: ";
+ for (int j = 0; j < elementTags[i].size(); ++j)
+ std::cout << elementTags[i][j] << " ";
+ std::cout << '\n';
+ std::cout << "\t\tnodes: ";
+ for (int j = 0; j < nodeTags[i].size(); ++j)
+ std::cout << nodeTags[i][j] << " ";
+ std::cout << '\n';
+ }
+ }
+ // -- code added to previous code (end) ------------------------------------
+
+ gmsh::finalize();
+ return 0;
+}
diff --git a/examples/gmsh_api/code5_elprops.cpp b/examples/gmsh_api/code5_elprops.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..218001b7c5731aa258bb57bc159127cdf842452b
--- /dev/null
+++ b/examples/gmsh_api/code5_elprops.cpp
@@ -0,0 +1,45 @@
+// This program displays info about "element type 1".
+//
+// How to build/run? (windows)
+// [in example/gmsh_api]
+// mkdir build
+// cd build
+// cmake .. && make && code5_elprops.exe
+
+#include <gmsh.h>
+#include <iostream>
+#include <map>
+
+int main(int argc, char **argv)
+{
+ gmsh::initialize();
+
+ // input parameter
+ int elementType = 1;
+
+ // outputs
+ std::string elementName;
+ int dim, order, numNodes, numPrimaryNodes;
+ std::vector<double> localNodeCoord;
+
+ // gmsh call
+ gmsh::model::mesh::getElementProperties(elementType,
+ elementName, dim, order,
+ numNodes, localNodeCoord,
+ numPrimaryNodes);
+
+ // display data returned by gmsh
+ std::cout << "Element type " << elementType << ":\n";
+ std::cout << "\tname = '" << elementName << "'\n";
+ std::cout << "\tdim = " << dim << '\n';
+ std::cout << "\torder = " << order << '\n';
+ std::cout << "\tnumNodes = " << numNodes << '\n';
+ std::cout << "\tlocalNodeCoord = ";
+ for(int i=0; i<localNodeCoord.size(); ++i)
+ std::cout << localNodeCoord[i] << " ";
+ std::cout << '\n';
+ std::cout << "\tnumPrimaryNodes = " << numPrimaryNodes << '\n';
+
+ gmsh::finalize();
+ return 0;
+}
diff --git a/examples/gmsh_api/code6_allnodes.cpp b/examples/gmsh_api/code6_allnodes.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c79187886f0af7d7137a950ee733a865ae87af48
--- /dev/null
+++ b/examples/gmsh_api/code6_allnodes.cpp
@@ -0,0 +1,50 @@
+// This program gets all the node tags and the coordinates of the model.
+//
+// How to build/run? (windows)
+// [in example/gmsh_api]
+// mkdir build
+// cd build
+// cmake .. && make && code6_allnodes.exe ..\rectangle.geo
+
+#include <gmsh.h>
+#include <iostream>
+
+int main(int argc, char **argv)
+{
+ if (argc < 2)
+ {
+ std::cout << "Usage: " << argv[0] << " <geo_file>\n";
+ return 0;
+ }
+
+ gmsh::initialize();
+
+ gmsh::open(argv[1]);
+ gmsh::model::mesh::generate(2);
+
+ // Renumber the node tags in a continuous sequence.
+ gmsh::model::mesh::renumberNodes();
+
+ // output parameters
+ std::vector<std::size_t> nodeTags;
+ std::vector<double> coord;
+ std::vector<double> parametricCoord;
+
+ // gmsh call
+ gmsh::model::mesh::getNodes(nodeTags, coord, parametricCoord);
+
+
+ // display the vectors
+ std::cout << "\tnodeTags = ";
+ for(int i=0; i<nodeTags.size(); ++i)
+ std::cout << nodeTags[i] << " ";
+ std::cout << '\n';
+
+ std::cout << "\tcoord = ";
+ for(int i=0; i<coord.size(); ++i)
+ std::cout << coord[i] << " ";
+ std::cout << '\n';
+
+ gmsh::finalize();
+ return 0;
+}
diff --git a/examples/gmsh_api/code7_parameters.cpp b/examples/gmsh_api/code7_parameters.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6fb94c42c09565875dc7fdf7fc0a2d4a5e829c24
--- /dev/null
+++ b/examples/gmsh_api/code7_parameters.cpp
@@ -0,0 +1,62 @@
+// This program takes a .geo file as argument, and reads the onelab parameters
+// defined by the "SetNumber" command.
+// This technique can be used to define additional parameters for the solver
+// such as the values for the boundary conditions or material parameters.
+//
+// How to build/run? (windows)
+// [in example/gmsh_api]
+// mkdir build
+// cd build
+// cmake .. && make && code7_parameters.exe ..\rectangle.geo
+
+#include <gmsh.h>
+#include <iostream>
+#include <sstream> // for std::stringstream
+
+int main(int argc, char **argv)
+{
+ // the program requires 1 argument: a .geo file.
+ if (argc < 2)
+ {
+ std::cout << "Usage: " << argv[0] << " <geo_file>\n";
+ return 0;
+ }
+
+ gmsh::initialize();
+
+ gmsh::open(argv[1]);
+
+ // get selected keys from the onelab database
+ std::vector<std::string> keys;
+ gmsh::onelab::getNames(keys, "(Boundary Conditions|Materials).+"); //see rectangle.geo file, retrieves boundary conditions and mat
+ for (auto &key : keys)
+ {
+ // get corresponding value
+ std::vector<double> value;
+ gmsh::onelab::getNumber(key, value);
+
+ // print received data
+ std::cout << key << ":";
+ for (auto &v : value)
+ std::cout << " " << v;
+ std::cout << '\n';
+
+ // expected key structure is "type/group_name/field"
+ // => split the key string into components
+ std::stringstream ss(key);
+ std::vector<std::string> words;
+ std::string word;
+ while(std::getline(ss, word, '/')) // read string until '/'
+ words.push_back(word);
+ if(words.size()==3)
+ {
+ std::cout << "\t.type of data = '" << words[0] << "'\n";
+ std::cout << "\t.physical group name = '" << words[1] << "'\n";
+ std::cout << "\t.field = '" << words[2] << "'\n";
+ std::cout << "\t.value = " << value[0] << '\n';
+ }
+ }
+
+ gmsh::finalize();
+ return 0;
+}
diff --git a/examples/gmsh_api/code8_jacobians.cpp b/examples/gmsh_api/code8_jacobians.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1add737d73b727c851b275af29ad33042ac8c92a
--- /dev/null
+++ b/examples/gmsh_api/code8_jacobians.cpp
@@ -0,0 +1,167 @@
+// This program demonstrates how to convert gmsh vectors to Eigen matrices.
+// It also shows:
+// - how to select a quadrature type and get the positions and weights of
+// the Gauss points.
+// - how to compute the inverse of the Jacobians.
+// - how to get the values of the shape functions and their derivatives at
+// the Gauss points.
+//
+// How to build/run? (windows)
+// [in example/gmsh_api]
+// mkdir build
+// cd build
+// cmake .. && make && code8_jacobians.exe ..\mono.geo
+
+#include <gmsh.h>
+#include <iostream>
+#include <map>
+#include <Eigen/Dense> // <= new header
+
+int main(int argc, char **argv)
+{
+ if (argc < 2)
+ {
+ std::cout << "Usage: " << argv[0] << " <geo_file>\n";
+ return 0;
+ }
+
+ gmsh::initialize();
+
+ gmsh::open(argv[1]);
+ gmsh::model::mesh::generate(2);
+
+ std::map<std::string, std::pair<int, int>> groups;
+ gmsh::vectorpair dimTags;
+ gmsh::model::getPhysicalGroups(dimTags);
+ for (int i = 0; i < dimTags.size(); ++i)
+ {
+ int dim = dimTags[i].first;
+ int tag = dimTags[i].second;
+ std::string name;
+ gmsh::model::getPhysicalName(dim, tag, name);
+ groups[name] = {dim, tag};
+ }
+
+ std::string groupname = "domain"; // <= we choose a physical group here!
+
+ int dim = groups[groupname].first;
+ int tag = groups[groupname].second;
+ if (tag == 0)
+ {
+ std::cerr << "Group '" << groupname << "' does not exist!\n";
+ return 1;
+ }
+
+ // get entities of the chosen physical group
+ std::vector<int> tags;
+ gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
+
+ std::cout << "Entities in group named '" << groupname << "': ";
+ for (int i = 0; i < tags.size(); ++i)
+ std::cout << tags[i] << " ";
+ std::cout << '\n';
+
+ // loop over entities
+ for (int k = 0; k < tags.size(); ++k)
+ {
+ std::cout << "Entity (" << dim << "," << tags[k] << "):\n";
+ std::vector<int> elementTypes;
+ std::vector<std::vector<std::size_t>> elementTags;
+ std::vector<std::vector<std::size_t>> nodeTags;
+ gmsh::model::mesh::getElements(elementTypes, elementTags,
+ nodeTags, dim, tags[k]);
+
+ // loop over element types
+ for (int i = 0; i < elementTypes.size(); ++i)
+ {
+ // -- CODE ADDED TO PREVIOUS EXAMPLE (begin) -----------------------
+
+ // get element name as a string
+ std::string elementName;
+ int edim, eorder, enumNodes, enumPrimaryNodes;
+ std::vector<double> localNodeCoord;
+ gmsh::model::mesh::getElementProperties(elementTypes[i],
+ elementName, dim, eorder,
+ enumNodes, localNodeCoord,
+ enumPrimaryNodes);
+ std::cout << '\t' << elementName << " order=" << eorder
+ << " nodes=" << enumNodes << '\n';
+
+ // get Gauss points coordinates and weights
+ std::vector<double> localCoords, weights;
+ gmsh::model::mesh::getIntegrationPoints(elementTypes[i],
+ "CompositeGauss2", // "Gauss2"
+ localCoords, weights);
+ std::cout << "\t" << weights.size() << " integration points\n";
+
+ std::cout << "\tweights: ";
+ for (int j = 0; j < weights.size(); ++j)
+ std::cout << weights[j] << " ";
+ std::cout << '\n';
+
+ std::cout << "\tlocalCoords (gmsh): ";
+ for (int j = 0; j < localCoords.size(); ++j)
+ std::cout << localCoords[j] << " ";
+ std::cout << '\n';
+
+ // convert to Eigen format
+ Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>>
+ localCoords_E(&localCoords[0], localCoords.size()/3, 3);
+ std::cout << "\tlocalCoords (Eigen):\n";
+ Eigen::IOFormat fmt(4, 0, ", ", "\n", "\t\t[", "]");
+ std::cout << localCoords_E.format(fmt) << '\n';
+
+ // get Jacobians
+ std::vector<double> jacobians, determinants, coords;
+ gmsh::model::mesh::getJacobians(elementTypes[i],
+ localCoords, jacobians,
+ determinants, coords);
+ std::cout << "\tGot " << determinants.size()
+ << " Jacobians for type " << elementTypes[i] << "\n";
+
+ std::cout << "\tjacobians (gmsh): ";
+ for (int j = 0; j < jacobians.size(); ++j)
+ std::cout << jacobians[j] << " ";
+ std::cout << '\n';
+
+ // convert first jacobian to Eigen format
+ Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> >
+ jacob(&jacobians[0], 3, 3);
+ std::cout << "\tfirst jacobian (Eigen):\n"
+ << jacob.format(fmt) << '\n';
+
+ // compute the inverse with Eigen
+ Eigen::Matrix<double, 3, 3> jacobinv = jacob.inverse();
+ std::cout << "\tinverse of the first jacobian (Eigen):\n"
+ << jacobinv.format(fmt) << "\n";
+
+ // Compute J*Jinv with Eigen "*" operator
+ std::cout << "\tVerification (J*Jinv):\n"
+ << (jacob*jacobinv).format(fmt) << '\n';
+
+ // shape functions and derivatives at the Gauss points
+ int numComponents, numOrientations;
+ std::vector<double> basisFunctions;
+ gmsh::model::mesh::getBasisFunctions(elementTypes[i], localCoords,
+ "Lagrange", numComponents,
+ basisFunctions,
+ numOrientations);
+ std::cout << "\tGot " << basisFunctions.size()
+ << " basis functions for type "
+ << elementTypes[i] << '\n';
+
+ gmsh::model::mesh::getBasisFunctions(elementTypes[i], localCoords,
+ "GradLagrange", numComponents,
+ basisFunctions,
+ numOrientations);
+ std::cout << "\tGot " << basisFunctions.size()
+ << " basis function gradients for type "
+ << elementTypes[i] << '\n';
+
+ // -- CODE ADDED TO PREVIOUS EXAMPLE (end) -------------------------
+ }
+ }
+
+ gmsh::finalize();
+ return 0;
+}
diff --git a/examples/gmsh_api/code9_views.cpp b/examples/gmsh_api/code9_views.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..20531988035b46a7e350d79cbb60c1d9bcd79bb2
--- /dev/null
+++ b/examples/gmsh_api/code9_views.cpp
@@ -0,0 +1,164 @@
+// This example generate some results over a mesh
+//
+// How to use this example?
+//
+// 1.build the code:
+// cd build
+// cmake ..
+// make
+// 2.generate a mesh from a geo file:
+// gmsh -2 -order 3 ..\src\mymesh.geo
+// 3.run the program with the msh as argument
+// bin\myview.exe ..\src\mymesh.msh
+// 4.display the generated msh data
+// gmsh results.msh
+
+#include <gmsh.h>
+#include <iostream>
+
+int main(int argc, char **argv)
+{
+ if (argc < 2)
+ {
+ std::cout << "Usage: " << argv[0] << " <geo_file>\n";
+ return 0;
+ }
+
+ gmsh::initialize();
+
+ gmsh::open(argv[1]);
+ gmsh::model::mesh::generate(2);
+
+ std::vector<std::string> names;
+ gmsh::model::list(names);
+ std::cout << "the file contains " << names.size() << " model names\n";
+
+ // get all the nodes
+ std::vector<std::size_t> nodeTags;
+ std::vector<double> coords;
+ std::vector<double> parametricCoord;
+ gmsh::model::mesh::getNodes(nodeTags, coords, parametricCoord);
+
+ std::cout << "the mesh has " << nodeTags.size() << " nodes\n";
+
+ // get all the elements for dim=2
+ std::vector<int> elementTypes;
+ std::vector<std::vector<std::size_t>> elementTags;
+ std::vector<std::vector<std::size_t>> elnodeTags;
+ gmsh::model::mesh::getElements(elementTypes, elementTags, elnodeTags, 2);
+
+ // compute total number of 2D elements
+ int nbelems=0;
+ for(std::size_t i=0; i< elementTags.size(); ++i)
+ nbelems+=elementTags[i].size();
+ std::cout << "the mesh has " << nbelems << " elements\n";
+
+ // Create a new post-processing view
+ int viewtag1 = gmsh::view::add("NodeData");
+ int viewtag2 = gmsh::view::add("ElementNodeData");
+
+ std::size_t nstep = 20;
+ double Lx = 20.;
+ double v = 10.0;
+ double tend = Lx / v;
+ for (std::size_t step = 0; step < nstep; step++)
+ {
+ double time = tend * ((double)step / nstep);
+
+ // ** NodeData
+ // loop over the nodes and compute some values from their coordinates
+ std::vector<std::vector<double>> data1(nodeTags.size()); // check this size!!
+ for (std::size_t i = 0; i < nodeTags.size(); ++i)
+ {
+ std::size_t tag = nodeTags[i];
+ double x = coords[i * 3 + 0];
+ double y = coords[i * 3 + 1];
+ double z = coords[i * 3 + 2];
+
+ double val = sin(2 * M_PI * (x - v * time) / Lx * 2) + y;
+ data1[i].resize(1);
+ data1[i][0] = val;
+ }
+ gmsh::view::addModelData(viewtag1, step, names[0], "NodeData",
+ nodeTags, data1, time);
+
+ // ** ElementNodeData
+ // loop over the elements and compute some values at the coords of its nodes
+
+ std::vector<std::vector<double>> data2(nbelems);
+ std::vector<std::size_t> elems2D(nbelems);
+
+ // for each element type
+ //std::cout << "there are " << elementTypes.size() << " element type(s)\n";
+ int k=0;
+ for (std::size_t ityp = 0; ityp < elementTypes.size(); ++ityp)
+ {
+ // get info about this type of element
+ std::string name;
+ int dim, order, numNodes, numPrimaryNodes;
+ std::vector<double> paramCoord;
+ gmsh::model::mesh::getElementProperties(elementTypes[ityp], name, dim, order,
+ numNodes, paramCoord, numPrimaryNodes);
+ //std::cout << "element type has " << numNodes << " nodes\n";
+
+ // select element/node tags for this type of element
+ auto &etags = elementTags[ityp];
+ auto &enods = elnodeTags[ityp];
+
+ // loop over elements of type "ityp"
+ //std::cout << "loop over the " << etags.size() << " elements\n";
+ for (std::size_t i = 0; i < etags.size(); ++i)
+ {
+ elems2D[k] = etags[i];
+
+ double offset = (elems2D[k]%2) *0.5; // so that fields are discontinious across elements
+
+ data2[k].resize(numNodes);
+ // loop over nodes of the element
+ for (std::size_t j = 0; j < numNodes; ++j)
+ {
+ int nodeTag = enods[numNodes * i + j];
+
+ std::vector<double> nodecoord;
+ std::vector<double> nodepcoord;
+ int dim, tag;
+ gmsh::model::mesh::getNode(nodeTag, nodecoord, nodepcoord, dim, tag);
+ double x = nodecoord[0];
+ double y = nodecoord[1];
+ double z = nodecoord[2];
+/*
+ double x = coords[nod * 3 + 0];
+ double y = coords[nod * 3 + 1];
+ double z = coords[nod * 3 + 2];
+*/
+ double val = sin(2 * M_PI * (x - v * time) / Lx * 2) + y;
+ data2[k][j] = val + offset;
+ }
+ k++;
+ }
+ }
+ gmsh::view::addModelData(viewtag2, step, names[0], "ElementNodeData",
+ elems2D, data2, time, 1); // the last ,1 is important!
+
+ }
+
+ // save the results to disk
+ gmsh::view::write(viewtag1, "results.msh");
+ gmsh::view::write(viewtag2, "results.msh", true);
+
+ // set display options & view results
+ gmsh::option::setNumber("View[0].IntervalsType", 3);
+ gmsh::option::setNumber("View[0].AdaptVisualizationGrid", 1);
+ gmsh::option::setNumber("View[0].MaxRecursionLevel", 3);
+ gmsh::option::setNumber("View[0].TargetError", -0.0001);
+
+ gmsh::option::setNumber("View[1].IntervalsType", 3);
+ gmsh::option::setNumber("View[1].AdaptVisualizationGrid", 1);
+ gmsh::option::setNumber("View[1].MaxRecursionLevel", 3);
+ gmsh::option::setNumber("View[1].TargetError", -0.0001);
+
+ gmsh::fltk::run();
+
+ gmsh::finalize();
+ return 0;
+}
diff --git a/examples/gmsh_api/mono.geo b/examples/gmsh_api/mono.geo
new file mode 100644
index 0000000000000000000000000000000000000000..aa4b82811cfaccd441476cb733ea10af6508a25d
--- /dev/null
+++ b/examples/gmsh_api/mono.geo
@@ -0,0 +1,22 @@
+Lx = 10;
+Ly = 10;
+nx = 1;
+ny = 1;
+
+Point(1) = {0, 0, 0, 1.0};
+Point(2) = {Lx, 0, 0, 1.0};
+Point(3) = {Lx, Ly, 0, 1.0};
+Point(4) = {0, Ly, 0, 1.0};
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+Curve Loop(1) = {4, 1, 2, 3};
+Plane Surface(1) = {1};
+Transfinite Curve {3, 1} = nx+1 Using Progression 1;
+Transfinite Curve {4, 2} = ny+1 Using Progression 1;
+Transfinite Surface {1};
+
+Recombine Surface {1}; // quads instead of triangles
+
+Physical Surface("domain", 6) = {1};
diff --git a/examples/gmsh_api/my_code.cpp b/examples/gmsh_api/my_code.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1cf01494c7ef97de5fcc074e0b094abf2913cfab
--- /dev/null
+++ b/examples/gmsh_api/my_code.cpp
@@ -0,0 +1,216 @@
+// This program takes a .geo file as argument, generates the mesh of the
+// surfaces and displays the result in the Gmsh window.
+//
+// This example introduces:
+// - the gmsh.h header
+// - initialize()/finalize() functions
+// - C++ namespaces / standard library (std::cout)
+//
+// How to build/run? (windows)
+// [in examples/gmsh_api]
+// mkdir build
+// cd build
+// cmake .. && make && code1_loadgeo.exe ..\rectangle.geo
+
+#include <gmsh.h> // gmsh main header (READ THIS FILE)
+#include <iostream> // for std::cout
+#include <map>
+#include <sstream> // for std::stringstream
+#include <Eigen/Dense> // Eigen library
+
+int main(int argc, char **argv)
+{
+ /*--------------INIT----------------*/
+ // the program requires 1 argument: a .geo file.
+ if (argc < 2)
+ {
+ std::cout << "Usage: " << argv[0] << " <geo_file>\n";
+ return 0;
+ }
+
+ gmsh::initialize(); // first thing to do before using any gmsh SDK function
+
+ gmsh::open(argv[1]); // similar to "File / Open" in the GUI
+ gmsh::model::mesh::generate(2); // similar to "Mesh / 2D" in the GUI
+
+ /*--------get physical groups (more particularly: the domain)--------------*/
+ std::map<std::string, std::pair<int, int>> groups;
+ gmsh::vectorpair dimTags;
+ gmsh::model::getPhysicalGroups(dimTags);
+ for (int i = 0; i < dimTags.size(); ++i)
+ {
+ int dim = dimTags[i].first;
+ int tag = dimTags[i].second;
+ std::string name;
+ gmsh::model::getPhysicalName(dim, tag, name);
+ groups[name] = {dim, tag};
+ }
+
+ /*-----------get physical properties of domain----------------------*/
+ std::vector<std::string> keys;
+ gmsh::onelab::getNames(keys, "(Materials).+");
+ double E;
+ double nu;
+ for (auto &key : keys)
+ {
+ // get corresponding value
+ std::vector<double> value;
+ gmsh::onelab::getNumber(key, value);
+ // expected key structure is "type/group_name/field"
+ // => split the key string into components
+ std::stringstream ss(key);
+ std::vector<std::string> words;
+ std::string word;
+ while(std::getline(ss, word, '/')) // read string until '/'
+ words.push_back(word);
+ if(words.size()==3)
+ {
+ if(words[2].compare("Young") == 0){
+ E = value[0];
+ }
+ if(words[2].compare("Poisson") == 0){
+ nu = value[0];
+ }
+ }
+ }
+ // std::cout << "E = " << E << "\n";
+ // std::cout << "nu = " << nu << "\n";
+
+ /*-----------------H matrix------------------*/
+ Eigen::Matrix<double, 3, 3> H_matrix;
+ H_matrix(0,0) = 1.;
+ H_matrix(1,0) = nu;
+ H_matrix(2,0) = 0.;
+ H_matrix(0,1) = nu;
+ H_matrix(1,1) = 1.;
+ H_matrix(2,1) = 0.;
+ H_matrix(0,2) = 0.;
+ H_matrix(1,2) = 0.;
+ H_matrix(2,2) = (1-nu)/2;
+ H_matrix = E/(1-nu*nu)*H_matrix;
+ // std::cout << H_matrix << '\n';
+
+ /*------------loop over elements of domain------------------------*/
+ std::string groupname = "domain";
+
+ int dim = groups[groupname].first;
+ int tag = groups[groupname].second;
+ if (tag == 0)
+ {
+ std::cerr << "Group '" << groupname << "' does not exist!\n";
+ return 1;
+ }
+
+ // get entities of the domain
+ std::vector<int> tags;
+ gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
+
+ std::cout << "Entities in group named '" << groupname << "': ";
+ for (int i = 0; i < tags.size(); ++i)
+ std::cout << tags[i] << " ";
+ std::cout << '\n';
+
+ // loop over entities
+ for (int k = 0; k < tags.size(); ++k)
+ {
+ std::cout << "Entity (" << dim << "," << tags[k] << "):\n";
+ std::vector<int> elementTypes;
+ std::vector<std::vector<std::size_t>> elementTags;
+ std::vector<std::vector<std::size_t>> nodeTags;
+ gmsh::model::mesh::getElements(elementTypes, elementTags,
+ nodeTags, dim, tags[k]);
+
+ // loop over element types
+ for (int i = 0; i < elementTypes.size(); ++i)
+ {
+ std::cout << "\telement type " << elementTypes[i] << ":\n";
+ std::cout << "\t\telements: ";
+ for (int j = 0; j < elementTags[i].size(); ++j)
+ std::cout << elementTags[i][j] << " ";
+ std::cout << '\n';
+
+ // get Gauss points coordinates and weights
+ std::vector<double> localCoords, weights;
+ gmsh::model::mesh::getIntegrationPoints(elementTypes[i],
+ "CompositeGauss2", // "Gauss2"
+ localCoords, weights);
+ std::cout << "\t" << weights.size() << " integration points\n";
+
+ std::cout << "\tweights: ";
+ for (int j = 0; j < weights.size(); ++j)
+ std::cout << weights[j] << " ";
+ std::cout << '\n';
+
+ std::cout << "\tlocalCoords (gmsh): ";
+ for (int j = 0; j < localCoords.size(); ++j)
+ std::cout << localCoords[j] << " ";
+ std::cout << '\n';
+
+ // convert to Eigen format
+ Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>>
+ localCoords_E(&localCoords[0], localCoords.size()/3, 3);
+ std::cout << "\tlocalCoords (Eigen):\n";
+ Eigen::IOFormat fmt(4, 0, ", ", "\n", "\t\t[", "]");
+ std::cout << localCoords_E.format(fmt) << '\n';
+
+ // get determinants and Jacobians
+ std::vector<double> jacobians, determinants, coords;
+ gmsh::model::mesh::getJacobians(elementTypes[i],
+ localCoords, jacobians,
+ determinants, coords);
+ std::cout << "\tGot " << determinants.size()
+ << " Jacobians for type " << elementTypes[i] << "\n";
+
+ std::cout << "\tdeterminants (gmsh): ";
+ for (int j = 0; j < determinants.size(); ++j)
+ std::cout << determinants[j] << " ";
+ std::cout << '\n';
+
+ std::cout << "\tjacobians (gmsh): ";
+ for (int j = 0; j < jacobians.size(); ++j)
+ std::cout << jacobians[j] << " ";
+ std::cout << '\n';
+
+ // derivatives of the shape functions at the Gauss points of the reference element
+ // we are here working with derivatives in the reference space, same for all elements
+ int numComponents, numOrientations;
+ std::vector<double> basisFunctions;
+ gmsh::model::mesh::getBasisFunctions(elementTypes[i], localCoords,
+ "GradLagrange", numComponents,
+ basisFunctions,
+ numOrientations);
+ std::cout << "\tbasis functions (gmsh): ";
+ for (int j = 0; j < basisFunctions.size(); ++j)
+ std::cout << basisFunctions[j] << " ";
+ std::cout << '\n';
+
+ /*---------computing the B matrix of the first element's (i.e elementTags[i][0]) first GP:------*/
+ // ne serait-ça pas plutôt utile de directement convertir les Jacobiens en 2D ? travailler en full 2D
+ // et faire pareil pour l'évaluation des gradients des basis functions ci-dessus ?
+ std::vector<double> b_matrix;
+
+ // convert first jacobian to Eigen format
+ Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> >
+ jacob(&jacobians[0], 3, 3);
+ std::cout << "\tfirst jacobian (Eigen):\n"
+ << jacob.format(fmt) << '\n';
+
+ // compute the inverse with Eigen
+ Eigen::Matrix<double, 3, 3> jacobinv = jacob.inverse();
+ std::cout << "\tinverse of the first jacobian (Eigen):\n"
+ << jacobinv.format(fmt) << "\n";
+
+ // compute the transpose of the inverse with Eigen
+ Eigen::Matrix<double, 3,3> jacobinvtrans = jacobinv.transpose();
+ std::cout << "\ttranspose of the inverse of the first jacobian (Eigen):\n"
+ << jacobinvtrans.format(fmt) << "\n";
+
+ for (int j = 0; j < weights.size(); ++j){ // looping over the number of gauss points ( = number of shape functions (A VERIF))
+
+ }
+ }
+ }
+
+ gmsh::finalize(); // clean-up things before ending the program
+ return 0;
+}
diff --git a/examples/gmsh_api/my_geo.geo b/examples/gmsh_api/my_geo.geo
new file mode 100644
index 0000000000000000000000000000000000000000..24307d2d5b28048c8ec38faf100ebb7faff4aef3
--- /dev/null
+++ b/examples/gmsh_api/my_geo.geo
@@ -0,0 +1,32 @@
+Lx = 5;
+Ly = 2;
+nx = 5;
+ny = 1;
+
+Point(1) = {0, 0, 0, 1.0};
+Point(2) = {Lx, 0, 0, 1.0};
+Point(3) = {Lx, Ly, 0, 1.0};
+Point(4) = {0, Ly, 0, 1.0};
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+Curve Loop(1) = {4, 1, 2, 3};
+Plane Surface(1) = {1};
+Transfinite Curve {3, 1} = nx+1 Using Progression 1;
+Transfinite Curve {4, 2} = ny+1 Using Progression 1;
+Transfinite Surface {1};
+
+Recombine Surface {1}; // quads instead of triangles
+
+Physical Curve("left_edge", 5) = {4};
+Physical Surface("domain", 6) = {1};
+Physical Curve("top_edge",7) = {3};
+
+// additional parameters given to the solver
+SetNumber("Boundary Conditions/left_edge/ux", 0.);
+SetNumber("Boundary Conditions/left_edge/uy", 0.);
+SetNumber("Materials/domain/Young", 210e3);
+SetNumber("Materials/domain/Poisson", 0.3);
+SetNumber("Boundary Conditions/top_edge/tx", 0.);
+SetNumber("Boundary Conditions/top_edge/ty", -100.);
\ No newline at end of file
diff --git a/examples/gmsh_api/rectangle.db b/examples/gmsh_api/rectangle.db
new file mode 100644
index 0000000000000000000000000000000000000000..2c3d83363d0063aefb0053171b65346027c353ad
Binary files /dev/null and b/examples/gmsh_api/rectangle.db differ
diff --git a/examples/gmsh_api/rectangle.geo b/examples/gmsh_api/rectangle.geo
new file mode 100644
index 0000000000000000000000000000000000000000..a0666de562275a53648095e91cc179a9ca36d4fa
--- /dev/null
+++ b/examples/gmsh_api/rectangle.geo
@@ -0,0 +1,32 @@
+Lx = 10;
+Ly = 2;
+nx = 5;
+ny = 1;
+
+Point(1) = {0, 0, 0, 1.0};
+Point(2) = {Lx, 0, 0, 1.0};
+Point(3) = {Lx, Ly, 0, 1.0};
+Point(4) = {0, Ly, 0, 1.0};
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+Curve Loop(1) = {4, 1, 2, 3};
+Plane Surface(1) = {1};
+Transfinite Curve {3, 1} = nx+1 Using Progression 1;
+Transfinite Curve {4, 2} = ny+1 Using Progression 1;
+Transfinite Surface {1};
+
+Recombine Surface {1}; // quads instead of triangles
+
+Physical Curve("left_edge", 5) = {4};
+Physical Surface("domain", 6) = {1};
+Physical Curve("top_edge",7) = {3};
+
+// additional parameters given to the solver
+SetNumber("Boundary Conditions/left_edge/ux", 0.);
+SetNumber("Boundary Conditions/left_edge/uy", 0.5);
+SetNumber("Materials/domain/Young", 210e3);
+SetNumber("Materials/domain/Poisson", 0.3);
+SetNumber("Boundary Conditions/top_edge/tx", 0.);
+SetNumber("Boundary Conditions/top_edge/ty", -100.);
\ No newline at end of file
diff --git a/examples/sparse_solve/CMakeLists.txt b/examples/sparse_solve/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c0c4e0c802b630d3aea205cbcee1b970c9775aba
--- /dev/null
+++ b/examples/sparse_solve/CMakeLists.txt
@@ -0,0 +1,52 @@
+PROJECT(SPARSE_SOLVE CXX)
+CMAKE_MINIMUM_REQUIRED(VERSION 3.12)
+
+# ------------------------------------------------------------------------------
+# Find libraries and setup compiler
+# ------------------------------------------------------------------------------
+
+# build type is "" by default in Linux
+IF(NOT CMAKE_BUILD_TYPE)
+ SET( CMAKE_BUILD_TYPE "Release" CACHE STRING "" FORCE)
+ENDIF()
+
+# enable C++11
+SET(CMAKE_CXX_STANDARD 11)
+SET(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+IF(APPLE)
+ # on macOS, do not give priority to frameworks/apps
+ SET(CMAKE_FIND_APPBUNDLE LAST)
+ SET(CMAKE_FIND_FRAMEWORK LAST)
+ENDIF()
+
+# find gmsh-sdk
+# gmsh.h
+FIND_PATH(GMSH_INCLUDE_DIRS NAMES "gmsh.h")
+MESSAGE(STATUS "GMSH_INCLUDE_DIRS=" ${GMSH_INCLUDE_DIRS})
+if(NOT GMSH_INCLUDE_DIRS)
+ MESSAGE(FATAL_ERROR "gmsh.h not found!")
+ENDIF()
+INCLUDE_DIRECTORIES(${GMSH_INCLUDE_DIRS})
+
+# libgmsh.so
+FIND_LIBRARY(GMSH_LIBRARIES gmsh)
+MESSAGE(STATUS "GMSH_LIBRARIES=" ${GMSH_LIBRARIES})
+IF(NOT GMSH_LIBRARIES)
+ MESSAGE(FATAL_ERROR "gmsh library not found!")
+ENDIF()
+
+# find Eigen
+find_path(EIGEN_INCLUDE_DIRS "Eigen/Dense"
+ PATHS "${PROJECT_SOURCE_DIR}/lib/eigen" "/usr/include/eigen3")
+MESSAGE(STATUS "EIGEN_INCLUDE_DIRS=" ${EIGEN_INCLUDE_DIRS})
+IF(NOT EIGEN_INCLUDE_DIRS)
+ MESSAGE(FATAL_ERROR "Eigen include dir not found!")
+ENDIF()
+INCLUDE_DIRECTORIES(${EIGEN_INCLUDE_DIRS})
+
+# ------------------------------------------------------------------------------
+
+SET(SRCS Tutorial_sparse_example.cpp Tutorial_sparse_example_details.cpp)
+ADD_EXECUTABLE(solve_sparse ${SRCS})
+TARGET_LINK_LIBRARIES(solve_sparse ${GMSH_LIBRARIES})
diff --git a/examples/sparse_solve/README.md b/examples/sparse_solve/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..3531cbb942b9fa3cd97b229fbade080bbd5ae45f
--- /dev/null
+++ b/examples/sparse_solve/README.md
@@ -0,0 +1,5 @@
+# sparse_solve
+
+This program solves a sparse linear system with Eigen.
+
+The files are directly taken from the Eigen source tree and put here without any modification.
diff --git a/examples/sparse_solve/Tutorial_sparse_example.cpp b/examples/sparse_solve/Tutorial_sparse_example.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f3587e82e3c8f4f2adbb72434e726c312a230eae
--- /dev/null
+++ b/examples/sparse_solve/Tutorial_sparse_example.cpp
@@ -0,0 +1,38 @@
+#include <Eigen/Sparse>
+#include <vector>
+#include <iostream>
+
+typedef Eigen::SparseMatrix<double> SpMat; // declares a column-major sparse matrix type of double
+typedef Eigen::Triplet<double> T;
+
+void buildProblem(std::vector<T>& coefficients, Eigen::VectorXd& b, int n);
+void saveAsBitmap(const Eigen::VectorXd& x, int n, const char* filename);
+
+int main(int argc, char** argv)
+{
+ // if(argc!=2) {
+ // std::cerr << "Error: expected one and only one argument.\n";
+ // return -1;
+ // }
+
+ int n = 300; // size of the image
+ int m = n*n; // number of unknows (=number of pixels)
+
+ // Assembly:
+ std::vector<T> coefficients; // list of non-zeros coefficients
+ Eigen::VectorXd b(m); // the right hand side-vector resulting from the constraints
+ buildProblem(coefficients, b, n);
+
+ SpMat A(m,m);
+ A.setFromTriplets(coefficients.begin(), coefficients.end());
+
+ // Solving:
+ Eigen::SimplicialCholesky<SpMat> chol(A); // performs a Cholesky factorization of A
+ Eigen::VectorXd x = chol.solve(b); // use the factorization to solve for the given right hand side
+
+ // Export the result to a file:
+ // saveAsBitmap(x, n, argv[1]);
+
+ return 0;
+}
+
diff --git a/examples/sparse_solve/Tutorial_sparse_example_details.cpp b/examples/sparse_solve/Tutorial_sparse_example_details.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2e628363b65f4a57eecbb7c137cfd263e255a8c8
--- /dev/null
+++ b/examples/sparse_solve/Tutorial_sparse_example_details.cpp
@@ -0,0 +1,44 @@
+#include <Eigen/Sparse>
+#include <vector>
+// #include <QImage>
+
+typedef Eigen::SparseMatrix<double> SpMat; // declares a column-major sparse matrix type of double
+typedef Eigen::Triplet<double> T;
+
+void insertCoefficient(int id, int i, int j, double w, std::vector<T>& coeffs,
+ Eigen::VectorXd& b, const Eigen::VectorXd& boundary)
+{
+ int n = int(boundary.size());
+ int id1 = i+j*n;
+
+ if(i==-1 || i==n) b(id) -= w * boundary(j); // constrained coefficient
+ else if(j==-1 || j==n) b(id) -= w * boundary(i); // constrained coefficient
+ else coeffs.push_back(T(id,id1,w)); // unknown coefficient
+}
+
+void buildProblem(std::vector<T>& coefficients, Eigen::VectorXd& b, int n)
+{
+ b.setZero();
+ Eigen::ArrayXd boundary = Eigen::ArrayXd::LinSpaced(n, 0,M_PI).sin().pow(2);
+ for(int j=0; j<n; ++j)
+ {
+ for(int i=0; i<n; ++i)
+ {
+ int id = i+j*n;
+ insertCoefficient(id, i-1,j, -1, coefficients, b, boundary);
+ insertCoefficient(id, i+1,j, -1, coefficients, b, boundary);
+ insertCoefficient(id, i,j-1, -1, coefficients, b, boundary);
+ insertCoefficient(id, i,j+1, -1, coefficients, b, boundary);
+ insertCoefficient(id, i,j, 4, coefficients, b, boundary);
+ }
+ }
+}
+
+// void saveAsBitmap(const Eigen::VectorXd& x, int n, const char* filename)
+// {
+// Eigen::Array<unsigned char,Eigen::Dynamic,Eigen::Dynamic> bits = (x*255).cast<unsigned char>();
+// QImage img(bits.data(), n,n,QImage::Format_Indexed8);
+// img.setColorCount(256);
+// for(int i=0;i<256;i++) img.setColor(i,qRgb(i,i,i));
+// img.save(filename);
+// }
diff --git a/lib/README.md b/lib/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..0cf64c7f9b1e252d8cdc426b5f4385e8b7dbfcbd
--- /dev/null
+++ b/lib/README.md
@@ -0,0 +1,3 @@
+# math0471/lib
+
+This folder contains scripts for downloading the 2 external libraries that are required for this project.
\ No newline at end of file
diff --git a/lib/get_eigen.cmd b/lib/get_eigen.cmd
new file mode 100644
index 0000000000000000000000000000000000000000..7930122a0a8ff970b4abed92e447d4f41ce4ad55
--- /dev/null
+++ b/lib/get_eigen.cmd
@@ -0,0 +1,18 @@
+::@echo off
+:: download eigen & extract to current folder
+
+set version=3.4.0
+set file=eigen-%version%.zip
+
+:: download file
+del %file% >nul 2>&1
+bitsadmin /transfer get_eigen /dynamic /download /priority foreground "https://gitlab.com/libeigen/eigen/-/archive/%version%/%file%" "%CD%\%file%"
+
+:: unzip
+rd /Q/S eigen-%version% >nul 2>&1
+powershell.exe -nologo -noprofile -command "& { Add-Type -A 'System.IO.Compression.FileSystem'; [IO.Compression.ZipFile]::ExtractToDirectory('%CD%\%file%', '%CD%'); }"
+del %file% >nul 2>&1
+
+:: rename folder
+rd /Q/S eigen >nul 2>&1
+ren eigen-%version% eigen
diff --git a/lib/get_eigen.sh b/lib/get_eigen.sh
new file mode 100644
index 0000000000000000000000000000000000000000..5f6e9e75f1fa310fd0e42eae2c757e5ba5f3c1ba
--- /dev/null
+++ b/lib/get_eigen.sh
@@ -0,0 +1,11 @@
+#!/bin/bash
+
+set -e
+
+VERSION=3.4.0
+
+wget -q https://gitlab.com/libeigen/eigen/-/archive/${VERSION}/eigen-${VERSION}.tar.gz
+tar xzf eigen-${VERSION}.tar.gz
+rm eigen-${VERSION}.tar.gz
+rm -rf eigen
+mv eigen-${VERSION} eigen
diff --git a/lib/get_gmsh.cmd b/lib/get_gmsh.cmd
new file mode 100644
index 0000000000000000000000000000000000000000..e960844e2b7b81e7b989979b1fcfb22c6e04a926
--- /dev/null
+++ b/lib/get_gmsh.cmd
@@ -0,0 +1,21 @@
+::@echo off
+:: download gmsh-sdk & extract to current folder
+
+set gmsh=gmsh-4.9.3-Windows64-sdk
+set file=%gmsh%.zip
+
+:: download file
+del %file% >nul 2>&1
+bitsadmin /transfer get_gmsh /dynamic /download /priority foreground "https://gmsh.info/bin/Windows/%file%" "%CD%\%file%"
+
+:: unzip
+rd /Q/S %gmsh% >nul 2>&1
+powershell.exe -nologo -noprofile -command "& { Add-Type -A 'System.IO.Compression.FileSystem'; [IO.Compression.ZipFile]::ExtractToDirectory('%CD%\%file%', '%CD%'); }"
+del %file% >nul 2>&1
+
+:: rename folder
+rd /Q/S gmsh-sdk >nul 2>&1
+ren %gmsh% gmsh-sdk
+
+:: copy dll to bin folder so that double-clic on gmsh.exe works!
+copy /Y gmsh-sdk\lib\gmsh-*.dll gmsh-sdk\bin\
diff --git a/lib/get_gmsh.sh b/lib/get_gmsh.sh
new file mode 100644
index 0000000000000000000000000000000000000000..0fb81224d41b54b17d0aed8b01f45ce743934347
--- /dev/null
+++ b/lib/get_gmsh.sh
@@ -0,0 +1,21 @@
+#!/bin/bash
+
+set -e
+
+VERSION=4.9.3
+
+if [[ "$OSTYPE" == "linux-gnu" ]]; then
+ wget -q https://gmsh.info/bin/Linux/gmsh-${VERSION}-Linux64-sdk.tgz
+ tar xzf gmsh-${VERSION}-Linux64-sdk.tgz
+ rm gmsh-${VERSION}-Linux64-sdk.tgz
+ rm -rf gmsh-sdk
+ mv gmsh-${VERSION}-Linux64-sdk gmsh-sdk
+elif [[ "$OSTYPE" == "darwin"* ]]; then
+ wget -q https://gmsh.info/bin/MacOSX/gmsh-${VERSION}-MacOSX-sdk.tgz
+ tar xzf gmsh-${VERSION}-MacOSX-sdk.tgz
+ rm gmsh-${VERSION}-MacOSX-sdk.tgz
+ rm -rf gmsh-sdk
+ mv gmsh-${VERSION}-MacOSX-sdk gmsh-sdk
+else
+ echo "unknown system"
+fi
diff --git a/srcs/CMakeLists.txt b/srcs/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e71010cdbebb28ac31889be22a7b02481935e378
--- /dev/null
+++ b/srcs/CMakeLists.txt
@@ -0,0 +1,52 @@
+PROJECT(MATH0471 CXX)
+CMAKE_MINIMUM_REQUIRED(VERSION 3.12)
+
+# ------------------------------------------------------------------------------
+# Find libraries and setup compiler
+# ------------------------------------------------------------------------------
+
+# build type is "" by default in Linux
+IF(NOT CMAKE_BUILD_TYPE)
+ SET( CMAKE_BUILD_TYPE "Release" CACHE STRING "" FORCE)
+ENDIF()
+
+# enable C++11
+SET(CMAKE_CXX_STANDARD 11)
+SET(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+IF(APPLE)
+ # on macOS, do not give priority to frameworks/apps
+ SET(CMAKE_FIND_APPBUNDLE LAST)
+ SET(CMAKE_FIND_FRAMEWORK LAST)
+ENDIF()
+
+# find gmsh-sdk
+# gmsh.h
+FIND_PATH(GMSH_INCLUDE_DIRS NAMES "gmsh.h")
+MESSAGE(STATUS "GMSH_INCLUDE_DIRS=" ${GMSH_INCLUDE_DIRS})
+if(NOT GMSH_INCLUDE_DIRS)
+ MESSAGE(FATAL_ERROR "gmsh.h not found!")
+ENDIF()
+INCLUDE_DIRECTORIES(${GMSH_INCLUDE_DIRS})
+
+# libgmsh.so
+FIND_LIBRARY(GMSH_LIBRARIES gmsh)
+MESSAGE(STATUS "GMSH_LIBRARIES=" ${GMSH_LIBRARIES})
+IF(NOT GMSH_LIBRARIES)
+ MESSAGE(FATAL_ERROR "gmsh library not found!")
+ENDIF()
+
+# find Eigen
+find_path(EIGEN_INCLUDE_DIRS "Eigen/Dense"
+ PATHS "${PROJECT_SOURCE_DIR}/lib/eigen" "/usr/include/eigen3")
+MESSAGE(STATUS "EIGEN_INCLUDE_DIRS=" ${EIGEN_INCLUDE_DIRS})
+IF(NOT EIGEN_INCLUDE_DIRS)
+ MESSAGE(FATAL_ERROR "Eigen include dir not found!")
+ENDIF()
+INCLUDE_DIRECTORIES(${EIGEN_INCLUDE_DIRS})
+
+# ------------------------------------------------------------------------------
+
+FILE(GLOB SRCS *.cpp)
+ADD_EXECUTABLE(solver ${SRCS})
+TARGET_LINK_LIBRARIES(solver ${GMSH_LIBRARIES})
diff --git a/srcs/README.md b/srcs/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..c265be7582c451b283e289136a4088a4474670cd
--- /dev/null
+++ b/srcs/README.md
@@ -0,0 +1,3 @@
+# srcs
+
+You can start your project here.
diff --git a/srcs/main.cpp b/srcs/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..05f39ef006504a19e39564482cee715b115b8940
--- /dev/null
+++ b/srcs/main.cpp
@@ -0,0 +1,8 @@
+#include <iostream>
+
+int main(int argc, char **argv)
+{
+ std::cout << "Hello world!\n";
+
+ return 0;
+}
\ No newline at end of file