diff --git a/README.md b/README.md
index 622bfec5ca7745d8c9f8736b4c1e71a7c073c0ae..4ffdffce87c0bb44d298705b243be0f8aec6b6f1 100644
--- a/README.md
+++ b/README.md
@@ -13,126 +13,3 @@ Some executable demos are given to detaill fonctionalities in /demo directory.
- Matthieu Aussal (matthieu.aussal@polytechnique.edu)
- Marc Bakry (marc.bakry@polytechnique.edu)
-
-
-
-
-
-
-# Notes on FEM (by Hadrien Montanelli)
-
-## triangledata class
-
-- triangledata - class for representing triangles
- - m_A - 1 x 3 vector - first node $\overrightarrow{A} = (A_0, A_1, A_2)$
- - m_B - 1 x 3 vector - second node $\overrightarrow{B} = (B_0, B_1, B_2)$
- - m_C - 1 x 3 vector - third node $\overrightarrow{C} = (C_0, C_1, C_2)$
- - m_D - 1 x 3 vector - midpoint of $\overrightarrow{AB} = \frac{1}{2}(\overrightarrow{A} + \overrightarrow{B})$
- - m_E - 1 x 3 vector - midpoint of $\overrightarrow{BC} = \frac{1}{2}(\overrightarrow{B} + \overrightarrow{C})$
- - m_F - 1 x 3 vector - midpoint of $\overrightarrow{CA} = \frac{1}{2}(\overrightarrow{C} + \overrightarrow{A})$
- - m_ctr - 1 x 3 vector - center of $ABC = \frac{1}{3}(\overrightarrow{A} + \overrightarrow{B} + \overrightarrow{C})$
- - m_nrm - 1 x 3 vector - unit normal $\overrightarrow{n} = \overrightarrow{AB}\times\overrightarrow{BC}/\Vert\overrightarrow{AB}\times\overrightarrow{BC}\Vert$
- - m_srf - scalar - surface $S = \frac{1}{2}\Vert\overrightarrow{AB}\times\overrightarrow{BC}\Vert$
- - m_lgt - 1 x 3 vector - edge lengths $= (\Vert\overrightarrow{BC}\Vert, \Vert\overrightarrow{CA}\Vert, \Vert\overrightarrow{AB}\Vert)$
- - m_stp - scalar - max of edge lengths $= \max(\Vert\overrightarrow{BC}\Vert, \Vert\overrightarrow{CA}\Vert, \Vert\overrightarrow{AB}\Vert)$
- - m_tau - 3 x 3 vector - unit edge tangents $= (\overrightarrow{\tau_{BC}}; \overrightarrow{\tau_{CA}}; \overrightarrow{\tau_{AB}})$ with:
- - $\overrightarrow{\tau_{BC}} = \overrightarrow{BC}/\Vert\overrightarrow{BC}\Vert$
- - $\overrightarrow{\tau_{CA}} = \overrightarrow{CA}/\Vert\overrightarrow{CA}\Vert$
- - $\overrightarrow{\tau_{AB}} = \overrightarrow{AB}/\Vert\overrightarrow{AB}\Vert$
- - m_nu - 3 x 3 matrix - unit edge normals $= (\overrightarrow{\nu_{BC}}; \overrightarrow{\nu_{CA}}; \overrightarrow{\nu_{AB}})$ with:
- - $\overrightarrow{\nu_{BC}} = \overrightarrow{BC}/\Vert\overrightarrow{BC}\Vert\times\overrightarrow{n}$
- - $\overrightarrow{\nu_{CA}} = \overrightarrow{CA}/\Vert\overrightarrow{CA}\Vert\times\overrightarrow{n}$
- - $\overrightarrow{\nu_{AB}} = \overrightarrow{AB}/\Vert\overrightarrow{AB}\Vert\times\overrightarrow{n}$
- - m_hgt - 1 x 3 vector - heights $= (h_{BC}, h_{CA}, h_{AB})$ with:
- - $h_{BC} = \overrightarrow{AB}\cdot\overrightarrow{BC}/\Vert\overrightarrow{BC}\Vert\times\overrightarrow{n}$
- - $h_{CA} = \overrightarrow{BC}\cdot\overrightarrow{CA}/\Vert\overrightarrow{CA}\Vert\times\overrightarrow{n}$
- - $h_{AB} = \overrightarrow{CA}\cdot\overrightarrow{AB}/\Vert\overrightarrow{AB}\Vert\times\overrightarrow{n}$
-- isfar - check if a point $P$ is close to a triangle $T$ by computing the distance^2 between $P$ and the center of $T$
-- intsing - exact integration of the kernel when isfar is false, that is, when $\overrightarrow{X_0} = (x_0, y_0, z_0)$ is close to $ABC$
- - xX0S - 1 x 3 vector - differences in $x = (A_0-x_0, B_0-x_0, C_0-x_0)$
- - yX0S - 1 x 3 vector - differences in $y = (A_1-y_0, B_1-y_0, C_1-y_0)$
- - zX0S - 1 x 3 vector - differences in $z = (A_2-z_0, B_2-z_0, C_2-z_0)$
- - nrmX0S - 1 x 3 vector - norms of the differences $= (\Vert\overrightarrow{A}-\overrightarrow{X_0}\Vert, \Vert\overrightarrow{B}-\overrightarrow{X_0}\Vert, \Vert\overrightarrow{C}-\overrightarrow{X_0}\Vert)$
- - h - scalar - point height to triangle $h = (\overrightarrow{A}-\overrightarrow{X_0})\cdot\overrightarrow{n} = (\overrightarrow{B}-\overrightarrow{X_0})\cdot\overrightarrow{n} = (\overrightarrow{C}-\overrightarrow{X_0})\cdot\overrightarrow{n}$
- - ca - 1 x 3 vector - cosines $= (\cos(AX_0B), \cos(BX_0C), \cos(CX_0A))$ with:
- - $\cos(AX_0B) = (\overrightarrow{A}-\overrightarrow{X_0})\cdot(\overrightarrow{B}-\overrightarrow{X_0})/(\Vert\overrightarrow{A}-\overrightarrow{X_0}\Vert\Vert\overrightarrow{B}-\overrightarrow{X_0}\Vert)$
- - $\cos(BX_0C) = (\overrightarrow{B}-\overrightarrow{X_0})\cdot(\overrightarrow{C}-\overrightarrow{X_0})/(\Vert\overrightarrow{B}-\overrightarrow{X_0}\Vert\Vert\overrightarrow{C}-\overrightarrow{X_0}\Vert)$
- - $\cos(CX_0A) = (\overrightarrow{C}-\overrightarrow{X_0})\cdot(\overrightarrow{A}-\overrightarrow{X_0})/(\Vert\overrightarrow{C}-\overrightarrow{X_0}\Vert\Vert\overrightarrow{A}-\overrightarrow{X_0}\Vert)$
- - sa - 1 x 3 vector - sines $= (\sin(AX_0B), \sin(BX_0C), \sin(CX_0A))$
- - omega - scalar - solid angle $\Omega$ $= -\pi + \phi_{AB} + \phi_{BC} + \phi_{CA}$ with:
- - $\phi_{AB} = \displaystyle\mathrm{acos}\left(\frac{\cos(AX_0B)-\cos(BX_0C)\cos(CX_0A)}{\sin(BX_0C)\sin(CX_0A)}\right)$
- - $\phi_{BC} = \displaystyle\mathrm{acos}\left(\frac{\cos(BX_0C)-\cos(CX_0A)\cos(AX_0B)}{\sin(CX_0A)\sin(AX_0B)}\right)$
- - $\phi_{CA} = \displaystyle\mathrm{acos}\left(\frac{\cos(CX_0A)-\cos(AX_0B)\cos(BX_0C)}{\sin(AX_0B)\sin(BX_0C)}\right)$
- - ps - scalar - projections along the tangents, its value depends on the edge/iteration count:
- - $(\overrightarrow{B}-\overrightarrow{X_0})\cdot\overrightarrow{\tau_{BC}}$ (i = 0), $(\overrightarrow{C}-\overrightarrow{X_0})\cdot\overrightarrow{\tau_{CA}}$ (i = 1), $(\overrightarrow{A}-\overrightarrow{X_0})\cdot\overrightarrow{\tau_{AB}}$ (i = 2)
- - psp1 - scalar - projections along the tangents, its value depends on the edge:
- - $(\overrightarrow{C}-\overrightarrow{X_0})\cdot\overrightarrow{\tau_{BC}}$ (i = 0), $(\overrightarrow{A}-\overrightarrow{X_0})\cdot\overrightarrow{\tau_{CA}}$ (i = 1), $(\overrightarrow{B}-\overrightarrow{X_0})\cdot\overrightarrow{\tau_{AB}}$ (i = 2)
- - ps2 - scalar - its value depends on the edge: [Not used.]
- - $\displaystyle\frac{(\overrightarrow{B}-\overrightarrow{X_0})\cdot\overrightarrow{\tau_{BC}}}{\Vert \overrightarrow{C}-\overrightarrow{X}_0\Vert\Vert \overrightarrow{B}-\overrightarrow{X}_0\Vert}$ (i = 0), $\displaystyle\frac{(\overrightarrow{C}-\overrightarrow{X_0})\cdot\overrightarrow{\tau_{CA}}}{\Vert \overrightarrow{A}-\overrightarrow{X}_0\Vert\Vert \overrightarrow{C}-\overrightarrow{X}_0\Vert}$ (i = 1), $\displaystyle\frac{(\overrightarrow{A}-\overrightarrow{X_0})\cdot\overrightarrow{\tau_{AB}}}{\Vert \overrightarrow{B}-\overrightarrow{X}_0\Vert\Vert \overrightarrow{A}-\overrightarrow{X}_0\Vert}$ (i = 2)
- - ar - scalar - absolute value of the diffrence between psp1 & ps, its value depends on the edge: [Not used.]
- - $\vert(\overrightarrow{C}-\overrightarrow{B})\cdot\overrightarrow{\tau_{BC}}\vert$ (i = 0), $\vert(\overrightarrow{A}-\overrightarrow{C})\cdot\overrightarrow{\tau_{CA}}\vert$ (i = 1), $\vert(\overrightarrow{B}-\overrightarrow{A})\cdot\overrightarrow{\tau_{AB}}\vert$ (i = 2)
- - xnu - 1 x 3 vector - components along the normals, its value depends on the edge:
- - $(\overrightarrow{B}-\overrightarrow{X}_0)-[(\overrightarrow{B}-\overrightarrow{X}_0)\cdot\overrightarrow{\tau_{BC}}]\overrightarrow{\tau_{BC}}$ (i = 0)
- - $(\overrightarrow{C}-\overrightarrow{X}_0)-[(\overrightarrow{C}-\overrightarrow{X}_0)\cdot\overrightarrow{\tau_{CA}}]\overrightarrow{\tau_{CA}}$ (i = 1)
- - $(\overrightarrow{A}-\overrightarrow{X}_0)-[(\overrightarrow{A}-\overrightarrow{X}_0)\cdot\overrightarrow{\tau_{AB}}]\overrightarrow{\tau_{AB}}$ (i = 2)
- - ah - scalar - norm of xnu
- - intaRm1 - scalar - its value depends on the edge and whether ah is smaller than tol=1e-7:
- - $\displaystyle \pm \log\left(\frac{\Vert\overrightarrow{C}-\overrightarrow{X_0}\Vert}{\Vert\overrightarrow{B}-\overrightarrow{X_0}\Vert}\right)$, sign depends on L/R wrt edge (i = 0, ah<tol)
- - $\displaystyle\mathrm{asinh}\left(\frac{(\overrightarrow{C}-\overrightarrow{X_0})\cdot\overrightarrow{\tau_{BC}}}{\Vert(\overrightarrow{B}-\overrightarrow{X}_0)-[(\overrightarrow{B}-\overrightarrow{X}_0)\cdot\overrightarrow{\tau_{BC}}]\overrightarrow{\tau_{BC}}\Vert}\right)$
-
- $\displaystyle-\,\mathrm{asinh}\left(\frac{(\overrightarrow{B}-\overrightarrow{X_0})\cdot\overrightarrow{\tau_{BC}}}{\Vert(\overrightarrow{B}-\overrightarrow{X}_0)-[(\overrightarrow{B}-\overrightarrow{X}_0)\cdot\overrightarrow{\tau_{BC}}]\overrightarrow{\tau_{BC}}\Vert}\right)$(i = 0, ah>tol)
-
- - $\displaystyle \pm \log\left(\frac{\Vert\overrightarrow{A}-\overrightarrow{X_0}\Vert}{\Vert\overrightarrow{C}-\overrightarrow{X_0}\Vert}\right)$, sign depends on L/R wrt edge (i = 1, ah<tol)
- - $\displaystyle\mathrm{asinh}\left(\frac{(\overrightarrow{A}-\overrightarrow{X_0})\cdot\overrightarrow{\tau_{CA}}}{\Vert(\overrightarrow{C}-\overrightarrow{X}_0)-[(\overrightarrow{C}-\overrightarrow{X}_0)\cdot\overrightarrow{\tau_{CA}}]\overrightarrow{\tau_{CA}}\Vert}\right)$
-
- $\displaystyle-\,\mathrm{asinh}\left(\frac{(\overrightarrow{C}-\overrightarrow{X_0})\cdot\overrightarrow{\tau_{CA}}}{\Vert(\overrightarrow{C}-\overrightarrow{X}_0)-[(\overrightarrow{C}-\overrightarrow{X}_0)\cdot\overrightarrow{\tau_{CA}}]\overrightarrow{\tau_{CA}}\Vert}\right)$(i = 1, ah>tol)
-
- - $\displaystyle \pm \log\left(\frac{\Vert\overrightarrow{B}-\overrightarrow{X_0}\Vert}{\Vert\overrightarrow{A}-\overrightarrow{X_0}\Vert}\right)$, sign depends on L/R wrt edge (i = 2, ah<tol)
- - $\displaystyle\mathrm{asinh}\left(\frac{(\overrightarrow{B}-\overrightarrow{X_0})\cdot\overrightarrow{\tau_{AB}}}{\Vert(\overrightarrow{A}-\overrightarrow{X}_0)-[(\overrightarrow{A}-\overrightarrow{X}_0)\cdot\overrightarrow{\tau_{AB}}]\overrightarrow{\tau_{AB}}\Vert}\right)$
-
- $\displaystyle-\,\mathrm{asinh}\left(\frac{(\overrightarrow{A}-\overrightarrow{X_0})\cdot\overrightarrow{\tau_{AB}}}{\Vert(\overrightarrow{A}-\overrightarrow{X}_0)-[(\overrightarrow{A}-\overrightarrow{X}_0)\cdot\overrightarrow{\tau_{AB}}]\overrightarrow{\tau_{AB}}\Vert}\right)$(i = 2, ah>tol)
-
- - Rm1 - scalar - exact integration of the kernel $= -h\Omega + I_0 + I_1 + I_2$ with:
- - $I_0 =$ intRm1 (i = 0) $\times\,[(\overrightarrow{B}-\overrightarrow{X_0})\cdot\overrightarrow{\nu_{BC}}]$
- - $I_1 =$ intRm1 (i = 1) $\times\,[(\overrightarrow{C}-\overrightarrow{X_0})\cdot\overrightarrow{\nu_{CA}}]$
- - $I_2 =$ intRm1 (i = 2) $\times\,[(\overrightarrow{A}-\overrightarrow{X_0})\cdot\overrightarrow{\nu_{AB}}]$
-- projectP1 - for a given triangle $ABC$ and a point $\overrightarrow{X}$, compute $\overrightarrow{P} = (P_0, P_1, P_2)$ with:
- - $P_0 = (\overrightarrow{B} - \overrightarrow{X})\cdot\overrightarrow{\nu_{BC}}/h_{BC}$
- - $P_1 = (\overrightarrow{C} - \overrightarrow{X})\cdot\overrightarrow{\nu_{CA}}/h_{CA}$
- - $P_2 =(\overrightarrow{A} - \overrightarrow{X})\cdot\overrightarrow{\nu_{AB}}/h_{AB}$
-- update - compute all the triangledata properties from vertex coordinates $\overrightarrow{A}$, $\overrightarrow{B}$ and $\overrightarrow{C}$
- - [x] Add midpoints
-
-## femdata class
-
-- femdata - class for representing triangular finite elements
- - m_vtx - n_vertex x 3 matrix - coordinates of all vertices (input of femdata)
- - m_elt2vtx - n_elem x 3 matrix - indices of the 3 vertices of each elem (input of femdata)
- - m_fet - string - type of finite elem, techincally a femtype, *e.g.*, lagrangeP1 (input of femdata)
- - m_ngs - scalar - number of Gauss points per elem, *e.g.*, 6 (input of femdata)
- - m_dof - n_DoFs x 3 matrix - coordinates of all DoFs (see m_buildDof)
- - m_elt2dof - n_elem x n_DoFs_per_elem matrix - indices of the DoFs of each elem (see m_buildDof)
- - m_dof2elt - n_DoFs_per_elem x n_elem matrix - inverse index mapping (see m_buildDof)
- - m_gss - (n_elem*n_Gauss_per_elem) x 3 matrix - coordinates of Gauss pts (see m_buildQuadrature)
- - m_elt2gss - n_elem x n_Gauss_per elem matrix - indices of the Gauss pts of each elem (see m_buildQuadrature)
- - m_gssnrm - (n_elem*n_Gauss_per_elem) x 3 matrix - unit normal at Gauss points (see m_buildQuadrature)
- - m_wgt - (n_elem*n_Gauss_per_elem) x 1 matrix - Gauss weights (see m_buildQuadrature)
- - m_base - (n_elem*n_Gauss_per_elem) x n_DoFs_per_elem matrix - function at Gauss pts (see m_buildBase)
- - m_grad - 3 x (n_elem*n_Gauss_per_elem) x n_DoFs_per_elem tensor - $\overrightarrow{\nabla}$ at Gauss pts (see m_buildGrad)
- - m_nxgrad - 3 x (n_elem*n_Gauss_per_elem) x n_DoFs_per_elem tensor - $\overrightarrow{n}\times\overrightarrow{\nabla}$ at Gauss pts (see m_buildGrad)
- - m_tridata - n_elem x 1 vector - each entry is an instance of triangledata (see m_buildGeometry)
-- dofsubdata -
-- refQuadrature - Gauss points and weights on the $(0,0)-(1,0)-(0,1)$ unit triangle
-- m_buildBase - definition of basis functions
- - [x] Add P2 functions
-- m_buildDof - degree of freedoms of elements
- - [ ] Add P2 functions
-- m_buildGrad - gradient of basis functions
- - [x] Add P2 functions
-- m_buildGeometry - get vertex coordinates from elements
-- m_buildQuadrature - Gauss points and weights on arbitrary triangle
-
-## Variational formulations
-
-- mass - assemble mass matrix
-- rigidity - assemble stiffness matrix