From d79f94355dc983c148557817d68e51712bfea31a Mon Sep 17 00:00:00 2001
From: blegouix <stilynx51@gmail.com>
Date: Mon, 13 Jan 2025 20:55:46 +0100
Subject: [PATCH 01/13] init

---
 docs/CMakeLists.txt           |  1 +
 docs/finite_element_module.md | 52 +++++++++++++++++++++++++++++++++++
 2 files changed, 53 insertions(+)
 create mode 100644 docs/finite_element_module.md

diff --git a/docs/CMakeLists.txt b/docs/CMakeLists.txt
index e8c5d0f0..a5c93c49 100644
--- a/docs/CMakeLists.txt
+++ b/docs/CMakeLists.txt
@@ -71,5 +71,6 @@ doxygen_add_docs(doc
 	"${CMAKE_CURRENT_SOURCE_DIR}/Running.md"
 	"${CMAKE_CURRENT_SOURCE_DIR}/tensor_module.md"
 	"${CMAKE_CURRENT_SOURCE_DIR}/exterior_module.md"
+	"${CMAKE_CURRENT_SOURCE_DIR}/finite_element_module.md"
         "${SimiLie_SOURCE_DIR}/include/similie/"
         ALL)
diff --git a/docs/finite_element_module.md b/docs/finite_element_module.md
new file mode 100644
index 00000000..59a4cbda
--- /dev/null
+++ b/docs/finite_element_module.md
@@ -0,0 +1,52 @@
+# Finite Element module {#finite_element_module}
+<!--
+SPDX-FileCopyrightText: 2025 Baptiste Legouix
+SPDX-License-Identifier: GPL-3.0-or-later
+-->
+
+## Why partial differential equations are importants for physics ?
+
+Given a multidimensional function \f$u:x\rightarrow \; .\f$, any [partial differential equation](https://en.wikipedia.org/wiki/Partial_differential_equation) can be written:
+
+\f\[
+\mathcal{F}(D^k u, D^{k-1} u,... ,D u, u, x) = f
+\f\]
+
+Where \f$D^k\f$ is a \f$k-\f$order differential operator. Any partial differential equation can be splitted into a Hamiltonian and non-Hamiltonian parts:
+
+\f\[
+\frac{\partial u}{\partial t} + \{\mathcal{H}, u\} + \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) = f
+\f\]
+
+Where \f$\{\mathcal{H}, u\}\f$ is the Poisson bracket between the Hamiltonian \f$\mathcal{H}\f$ of the system and \f$u\f$. It can always be written \f$\{\mathcal{H}, u\} = \mathcal{J}\frac{\partial \mathcal{H}}{\partial x}\f$ in term of the symplectic matrix \f$\mathcal{J}\f$. The Hamiltonian part of the PDE is associated to conserved quantities such as mass, momentum or energy.
+
+\note If a PDE involves a second-order temporal derivative \f$\frac{\partial^2 u}{\partial t^2}\f$ (like in wave equation) we can always define \f$U = (\frac{\partial u}{\partial t}, u)\f$ and write the PDE in term of  \f$U\f$.
+
+Partial differential equations are the main mathematical tool to describe rules of classical physics.
+
+## What is Finite Element Method ?
+
+Finite Element Method aims to discretize PDE to produce numerical schemes to solve it. It relies on the weak formulation of the PDE:
+
+\f\[
+\int \frac{\partial u}{\partial t}v\; d\Omega + \int \mathcal{J}\mathcal{H}\frac{\partial v}{\partial x} \; d\Omega + \int \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) v \; d\Omega = \int f v \; d\Omega
+\f\]
+
+\note \f$\int \mathcal{J}\frac{\partial \mathcal{H}}{\partial x}v \; d\Omega\f$ has been transformed into \f$\int \mathcal{J}\mathcal{H}\frac{\partial v}{\partial x} \; d\Omega\f$ by part integration. 
+
+Where \f$v\f$ is a test function belonging to a basis of test functions.
+
+We decompose \f$u\f$ and \f$v\f$ in a discrete basis function \f$\phi\f$ (we consider both \f$u\f$ and \f$v\f$ to be approximated in the same discrete basis function but in principle we could distinguish two different basis function). 
+
+\f\[
+\left\{\begin{matrix}
+u = \sum_i \tilde{u}_i \phi_i\\
+v = \sum_i \tilde{v}_i \phi_i
+\end{matrix}\right.
+\f\]
+
+Leading to:
+
+\f\[
+\int \frac{\partial \tilde{u}_i}{\partial t}\tilde{v}_j\phi_i\phi_j\; d\Omega + \int \mathcal{J}\mathcal{H}\tilde{v}_j\frac{\partial \phi_j}{\partial x} \; d\Omega + \int \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) v \; d\Omega = \int \tilde{f}_i \tilde{v}_j \phi_i\phi_j\; d\Omega
+\f\]

From 3958431088ce210a863d4134593c496d95c76c15 Mon Sep 17 00:00:00 2001
From: blegouix <stilynx51@gmail.com>
Date: Tue, 14 Jan 2025 14:53:31 +0100
Subject: [PATCH 02/13] wip

---
 docs/finite_element_module.md | 25 +++++++++++++++++++------
 1 file changed, 19 insertions(+), 6 deletions(-)

diff --git a/docs/finite_element_module.md b/docs/finite_element_module.md
index 59a4cbda..ad14bd37 100644
--- a/docs/finite_element_module.md
+++ b/docs/finite_element_module.md
@@ -18,7 +18,7 @@ Where \f$D^k\f$ is a \f$k-\f$order differential operator. Any partial differenti
 \frac{\partial u}{\partial t} + \{\mathcal{H}, u\} + \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) = f
 \f\]
 
-Where \f$\{\mathcal{H}, u\}\f$ is the Poisson bracket between the Hamiltonian \f$\mathcal{H}\f$ of the system and \f$u\f$. It can always be written \f$\{\mathcal{H}, u\} = \mathcal{J}\frac{\partial \mathcal{H}}{\partial x}\f$ in term of the symplectic matrix \f$\mathcal{J}\f$. The Hamiltonian part of the PDE is associated to conserved quantities such as mass, momentum or energy.
+Where \f$\{\mathcal{H}, u\}\f$ is the Poisson bracket between the Hamiltonian \f$\mathcal{H}\f$ of the system and \f$u\f$. It can always be written \f$\{\mathcal{H}, u\} = \mathcal{J}\frac{d\mathcal{H}(u, t)}{d x}\f$ in term of the symplectic matrix \f$\mathcal{J}\f$. The Hamiltonian part of the PDE is associated to conserved quantities such as mass, momentum or energy.
 
 \note If a PDE involves a second-order temporal derivative \f$\frac{\partial^2 u}{\partial t^2}\f$ (like in wave equation) we can always define \f$U = (\frac{\partial u}{\partial t}, u)\f$ and write the PDE in term of  \f$U\f$.
 
@@ -29,10 +29,23 @@ Partial differential equations are the main mathematical tool to describe rules
 Finite Element Method aims to discretize PDE to produce numerical schemes to solve it. It relies on the weak formulation of the PDE:
 
 \f\[
-\int \frac{\partial u}{\partial t}v\; d\Omega + \int \mathcal{J}\mathcal{H}\frac{\partial v}{\partial x} \; d\Omega + \int \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) v \; d\Omega = \int f v \; d\Omega
+\int \frac{\partial u}{\partial t}v\; d\Omega + \int \mathcal{J}\frac{d \mathcal{H}(u, t)}{dx}v\; d\Omega + \int \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) v \; d\Omega = \int f v \; d\Omega
 \f\]
 
-\note \f$\int \mathcal{J}\frac{\partial \mathcal{H}}{\partial x}v \; d\Omega\f$ has been transformed into \f$\int \mathcal{J}\mathcal{H}\frac{\partial v}{\partial x} \; d\Omega\f$ by part integration. 
+In which the total derivative of the Hamiltonian can be expanded using the [functional derivative](https://en.wikipedia.org/wiki/Functional_derivative):
+
+\f\[
+\int \frac{\partial u}{\partial t}v\; d\Omega + \int \mathcal{J}\left(\frac{\partial \mathcal{H}(u, t)}{\partial x} - \nabla\cdot \frac{\partial \mathcal{H}(u, t)}{\partial (\nabla x)}\right)v\; d\Omega + \int \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) v \; d\Omega = \int f v \; d\Omega
+\f\]
+
+Leading to:
+
+\f\[
+\int \frac{\partial u}{\partial t}v\; d\Omega + \int \mathcal{J}\frac{\partial \mathcal{H}(u, t)}{\partial x} v\; d\Omega - \int \mathcal{J}\frac{\partial \mathcal{H}(u, t)}{\partial (\nabla x)}\nabla v\; d\Omega + \int \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) v \; d\Omega = \int f v \; d\Omega
+\f\]
+
+
+\note \f$\int \mathcal{J}\frac{\partial \mathcal{H}(u, t)}{\partial x}v \; d\Omega\f$ has been transformed into \f$\int \mathcal{J}\mathcal{H}(u,t)\frac{\partial v}{\partial x} \; d\Omega\f$ by part integration. 
 
 Where \f$v\f$ is a test function belonging to a basis of test functions.
 
@@ -40,13 +53,13 @@ We decompose \f$u\f$ and \f$v\f$ in a discrete basis function \f$\phi\f$ (we con
 
 \f\[
 \left\{\begin{matrix}
-u = \sum_i \tilde{u}_i \phi_i\\
-v = \sum_i \tilde{v}_i \phi_i
+u(x, t) = \sum_i \tilde{u}_i \phi_i(x, t)\\
+v(x, t) = \sum_j \tilde{v}_j \psi_j(x, t)
 \end{matrix}\right.
 \f\]
 
 Leading to:
 
 \f\[
-\int \frac{\partial \tilde{u}_i}{\partial t}\tilde{v}_j\phi_i\phi_j\; d\Omega + \int \mathcal{J}\mathcal{H}\tilde{v}_j\frac{\partial \phi_j}{\partial x} \; d\Omega + \int \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) v \; d\Omega = \int \tilde{f}_i \tilde{v}_j \phi_i\phi_j\; d\Omega
+\frac{\partial \tilde{u}_i}{\partial t}\tilde{v}_j \int \phi_i\psi_j\; d\Omega + \tilde{v}_j\int \mathcal{J}\mathcal{H}(u, t)\frac{\partial \psi_j}{\partial x} \; d\Omega + \tilde{v}_j\int \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) \psi_j \; d\Omega = \tilde{f}_i \tilde{v}_j \int \phi_i\psi_j\; d\Omega
 \f\]

From db4707fbd9c5163c9e50229a2aca81ed587b34bc Mon Sep 17 00:00:00 2001
From: blegouix <stilynx51@gmail.com>
Date: Tue, 14 Jan 2025 15:48:20 +0100
Subject: [PATCH 03/13] wip

---
 docs/finite_element_module.md | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/docs/finite_element_module.md b/docs/finite_element_module.md
index ad14bd37..7b5198ba 100644
--- a/docs/finite_element_module.md
+++ b/docs/finite_element_module.md
@@ -26,22 +26,22 @@ Partial differential equations are the main mathematical tool to describe rules
 
 ## What is Finite Element Method ?
 
-Finite Element Method aims to discretize PDE to produce numerical schemes to solve it. It relies on the weak formulation of the PDE:
+Finite Element Method aims to discretize PDE to produce numerical schemes to solve it. It relies on the weak formulation of the PDE (\f$\left<\cdot, \cdot\right>\f$ is the inner product between functions):
 
 \f\[
-\int \frac{\partial u}{\partial t}v\; d\Omega + \int \mathcal{J}\frac{d \mathcal{H}(u, t)}{dx}v\; d\Omega + \int \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) v \; d\Omega = \int f v \; d\Omega
+\left<\frac{\partial u}{\partial t}, v\right> + \left< \mathcal{J}\frac{d \mathcal{H}(u, t)}{dx}, v\right> + \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), v \right> = \left< f, v\right>
 \f\]
 
 In which the total derivative of the Hamiltonian can be expanded using the [functional derivative](https://en.wikipedia.org/wiki/Functional_derivative):
 
 \f\[
-\int \frac{\partial u}{\partial t}v\; d\Omega + \int \mathcal{J}\left(\frac{\partial \mathcal{H}(u, t)}{\partial x} - \nabla\cdot \frac{\partial \mathcal{H}(u, t)}{\partial (\nabla x)}\right)v\; d\Omega + \int \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) v \; d\Omega = \int f v \; d\Omega
+\left<\frac{\partial u}{\partial t}, v\right> + \left< \mathcal{J}\left(\frac{\partial \mathcal{H}(u, t)}{\partial x} - \nabla\cdot \frac{\partial \mathcal{H}(u, t)}{\partial (\nabla x)}\right), v\right> + \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), v \right> = \left< f, v\right>
 \f\]
 
 Leading to:
 
 \f\[
-\int \frac{\partial u}{\partial t}v\; d\Omega + \int \mathcal{J}\frac{\partial \mathcal{H}(u, t)}{\partial x} v\; d\Omega - \int \mathcal{J}\frac{\partial \mathcal{H}(u, t)}{\partial (\nabla x)}\nabla v\; d\Omega + \int \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) v \; d\Omega = \int f v \; d\Omega
+\left<\frac{\partial u}{\partial t}, v\right> + \left< \mathcal{J}\frac{\partial \mathcal{H}(u, t)}{\partial x}, v\right> - \left< \mathcal{J}\frac{\partial \mathcal{H}(u, t)}{\partial (\nabla x)}, \nabla v\right> + \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), v \right> = \left< f, v\right>
 \f\]
 
 

From 609789ce687eadd16d535aeab7ed8dc8e08fac23 Mon Sep 17 00:00:00 2001
From: blegouix <stilynx51@gmail.com>
Date: Wed, 15 Jan 2025 12:46:57 +0100
Subject: [PATCH 04/13] wip

---
 docs/finite_element_module.md | 29 +++++++++++++++--------------
 1 file changed, 15 insertions(+), 14 deletions(-)

diff --git a/docs/finite_element_module.md b/docs/finite_element_module.md
index 7b5198ba..4152d4ee 100644
--- a/docs/finite_element_module.md
+++ b/docs/finite_element_module.md
@@ -9,13 +9,13 @@ SPDX-License-Identifier: GPL-3.0-or-later
 Given a multidimensional function \f$u:x\rightarrow \; .\f$, any [partial differential equation](https://en.wikipedia.org/wiki/Partial_differential_equation) can be written:
 
 \f\[
-\mathcal{F}(D^k u, D^{k-1} u,... ,D u, u, x) = f
+\mathcal{F}(D^k u, D^{k-1} u,... ,D u, u, x) = 0
 \f\]
 
 Where \f$D^k\f$ is a \f$k-\f$order differential operator. Any partial differential equation can be splitted into a Hamiltonian and non-Hamiltonian parts:
 
 \f\[
-\frac{\partial u}{\partial t} + \{\mathcal{H}, u\} + \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) = f
+\frac{\partial u}{\partial t} + \{\mathcal{H}, u\} = \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) 
 \f\]
 
 Where \f$\{\mathcal{H}, u\}\f$ is the Poisson bracket between the Hamiltonian \f$\mathcal{H}\f$ of the system and \f$u\f$. It can always be written \f$\{\mathcal{H}, u\} = \mathcal{J}\frac{d\mathcal{H}(u, t)}{d x}\f$ in term of the symplectic matrix \f$\mathcal{J}\f$. The Hamiltonian part of the PDE is associated to conserved quantities such as mass, momentum or energy.
@@ -29,37 +29,38 @@ Partial differential equations are the main mathematical tool to describe rules
 Finite Element Method aims to discretize PDE to produce numerical schemes to solve it. It relies on the weak formulation of the PDE (\f$\left<\cdot, \cdot\right>\f$ is the inner product between functions):
 
 \f\[
-\left<\frac{\partial u}{\partial t}, v\right> + \left< \mathcal{J}\frac{d \mathcal{H}(u, t)}{dx}, v\right> + \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), v \right> = \left< f, v\right>
+\left<\frac{\partial u}{\partial t}, \psi_j\right> + \left< \mathcal{J}\frac{d \mathcal{H}(u, t)}{dx}, \psi_j\right> = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>
 \f\]
 
 In which the total derivative of the Hamiltonian can be expanded using the [functional derivative](https://en.wikipedia.org/wiki/Functional_derivative):
 
 \f\[
-\left<\frac{\partial u}{\partial t}, v\right> + \left< \mathcal{J}\left(\frac{\partial \mathcal{H}(u, t)}{\partial x} - \nabla\cdot \frac{\partial \mathcal{H}(u, t)}{\partial (\nabla x)}\right), v\right> + \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), v \right> = \left< f, v\right>
+\left<\frac{\partial u}{\partial t}, \psi_j\right> + \left< \mathcal{J}\left(\frac{\partial \mathcal{H}(u, t)}{\partial x} - \nabla\cdot \frac{\partial \mathcal{H}(u, t)}{\partial (\nabla x)}\right), \psi_j\right> = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>
 \f\]
 
 Leading to:
 
 \f\[
-\left<\frac{\partial u}{\partial t}, v\right> + \left< \mathcal{J}\frac{\partial \mathcal{H}(u, t)}{\partial x}, v\right> - \left< \mathcal{J}\frac{\partial \mathcal{H}(u, t)}{\partial (\nabla x)}, \nabla v\right> + \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), v \right> = \left< f, v\right>
+\left<\frac{\partial u}{\partial t}, \psi_j\right> + \left< \mathcal{J}\frac{\partial \mathcal{H}(u, t)}{\partial x}, \psi_j\right> - \left< \mathcal{J}\frac{\partial \mathcal{H}(u, t)}{\partial (\nabla x)}, \nabla \psi_j\right> = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>
 \f\]
 
 
-\note \f$\int \mathcal{J}\frac{\partial \mathcal{H}(u, t)}{\partial x}v \; d\Omega\f$ has been transformed into \f$\int \mathcal{J}\mathcal{H}(u,t)\frac{\partial v}{\partial x} \; d\Omega\f$ by part integration. 
+\note \f$\left< \mathcal{J}\nabla\cdot\frac{\partial \mathcal{H}(u, t)}{\partial (\nabla x)}, \psi_j\right>\f$ has been transformed into \f$\left< \mathcal{J}\frac{\partial \mathcal{H}(u, t)}{\partial (\nabla x)}, \nabla \psi_j\right>\f$ by part integration. 
 
-Where \f$v\f$ is a test function belonging to a basis of test functions.
+Where \f$\psi\f$ is a basis of test functions.
 
-We decompose \f$u\f$ and \f$v\f$ in a discrete basis function \f$\phi\f$ (we consider both \f$u\f$ and \f$v\f$ to be approximated in the same discrete basis function but in principle we could distinguish two different basis function). 
+We decompose \f$u\f$ and \f$\psi_j\f$ in a discrete basis function \f$\phi\f$. 
 
 \f\[
-\left\{\begin{matrix}
-u(x, t) = \sum_i \tilde{u}_i \phi_i(x, t)\\
-v(x, t) = \sum_j \tilde{v}_j \psi_j(x, t)
-\end{matrix}\right.
+u(x, t) = \sum_i \tilde{u}_i \phi_i(x, t)
 \f\]
 
-Leading to:
+\attention \f$\phi\f$ and \f$\psi\f$ may or may not be the same.
+
+The weak formulation is thus rewritten (Einstein summation convention applies):
 
 \f\[
-\frac{\partial \tilde{u}_i}{\partial t}\tilde{v}_j \int \phi_i\psi_j\; d\Omega + \tilde{v}_j\int \mathcal{J}\mathcal{H}(u, t)\frac{\partial \psi_j}{\partial x} \; d\Omega + \tilde{v}_j\int \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) \psi_j \; d\Omega = \tilde{f}_i \tilde{v}_j \int \phi_i\psi_j\; d\Omega
+\frac{\partial \tilde{u}_i}{\partial t}\left<\phi_i, \psi_j\right> + \left< \mathcal{J}\frac{\partial \mathcal{H}(\tilde{u}_i\phi_i, t)}{\partial x}, \psi_j\right> - \left< \mathcal{J}\frac{\partial \mathcal{H}(\tilde{u}_i\phi_i, t)}{\partial (\nabla x)}, \nabla \psi_j\right> = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>
 \f\]
+
+\note As an example, the Poisson equation \f$\Delta\phi = \rho\f$ is associated to the Hamiltonian \f$\mathcal{H} = \frac{1}{2}|\nabla \phi|^2 - \rho\phi\f$.

From 8947326c3ff96f0d39e0b253f9a0a2cf1e38ebd9 Mon Sep 17 00:00:00 2001
From: blegouix <stilynx51@gmail.com>
Date: Wed, 15 Jan 2025 14:00:01 +0100
Subject: [PATCH 05/13] wip

---
 docs/finite_element_module.md | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/docs/finite_element_module.md b/docs/finite_element_module.md
index 4152d4ee..eb92eef9 100644
--- a/docs/finite_element_module.md
+++ b/docs/finite_element_module.md
@@ -18,7 +18,7 @@ Where \f$D^k\f$ is a \f$k-\f$order differential operator. Any partial differenti
 \frac{\partial u}{\partial t} + \{\mathcal{H}, u\} = \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) 
 \f\]
 
-Where \f$\{\mathcal{H}, u\}\f$ is the Poisson bracket between the Hamiltonian \f$\mathcal{H}\f$ of the system and \f$u\f$. It can always be written \f$\{\mathcal{H}, u\} = \mathcal{J}\frac{d\mathcal{H}(u, t)}{d x}\f$ in term of the symplectic matrix \f$\mathcal{J}\f$. The Hamiltonian part of the PDE is associated to conserved quantities such as mass, momentum or energy.
+Where \f$\{\mathcal{H}, u\}\f$ is the Poisson bracket between the Hamiltonian \f$\mathcal{H}(u, \nabla{u}, t)\f$ of the system and \f$u\f$. It can always be written \f$\{\mathcal{H}, u\} = \mathcal{J}\frac{\delta \mathcal{H}}{\delta u}\f$ in term of the symplectic matrix \f$\mathcal{J}\f$ (where \f$\frac{\delta}{\delta u}\f$ is a [variational derivative](https://en.wikipedia.org/wiki/Functional_derivative)). The Hamiltonian part of the PDE is associated to conserved quantities such as mass, momentum or energy.
 
 \note If a PDE involves a second-order temporal derivative \f$\frac{\partial^2 u}{\partial t^2}\f$ (like in wave equation) we can always define \f$U = (\frac{\partial u}{\partial t}, u)\f$ and write the PDE in term of  \f$U\f$.
 
@@ -26,26 +26,26 @@ Partial differential equations are the main mathematical tool to describe rules
 
 ## What is Finite Element Method ?
 
-Finite Element Method aims to discretize PDE to produce numerical schemes to solve it. It relies on the weak formulation of the PDE (\f$\left<\cdot, \cdot\right>\f$ is the inner product between functions):
+Finite Element Method aims to discretize PDEs to produce numerical schemes able to solve them. It relies on the weak formulation of the PDEs (\f$\left<\cdot, \cdot\right>_\Omega\f$ is the inner product between functions on the domain \f$\Omega\f$):
 
 \f\[
-\left<\frac{\partial u}{\partial t}, \psi_j\right> + \left< \mathcal{J}\frac{d \mathcal{H}(u, t)}{dx}, \psi_j\right> = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>
+\left<\frac{\partial u}{\partial t}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\delta \mathcal{H}}{\delta u}, \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
 \f\]
 
-In which the total derivative of the Hamiltonian can be expanded using the [functional derivative](https://en.wikipedia.org/wiki/Functional_derivative):
+In which the variational derivative of the Hamiltonian can be expanded as: 
 
 \f\[
-\left<\frac{\partial u}{\partial t}, \psi_j\right> + \left< \mathcal{J}\left(\frac{\partial \mathcal{H}(u, t)}{\partial x} - \nabla\cdot \frac{\partial \mathcal{H}(u, t)}{\partial (\nabla x)}\right), \psi_j\right> = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>
+ \frac{\delta \mathcal{H}}{\delta u} = \frac{\partial \mathcal{H}}{\partial u} - \nabla\cdot \frac{\partial \mathcal{H}}{\partial (\nabla u)}
 \f\]
 
 Leading to:
 
 \f\[
-\left<\frac{\partial u}{\partial t}, \psi_j\right> + \left< \mathcal{J}\frac{\partial \mathcal{H}(u, t)}{\partial x}, \psi_j\right> - \left< \mathcal{J}\frac{\partial \mathcal{H}(u, t)}{\partial (\nabla x)}, \nabla \psi_j\right> = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>
+\left<\frac{\partial u}{\partial t}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial u}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \psi_j\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \nabla \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
 \f\]
 
 
-\note \f$\left< \mathcal{J}\nabla\cdot\frac{\partial \mathcal{H}(u, t)}{\partial (\nabla x)}, \psi_j\right>\f$ has been transformed into \f$\left< \mathcal{J}\frac{\partial \mathcal{H}(u, t)}{\partial (\nabla x)}, \nabla \psi_j\right>\f$ by part integration. 
+\note \f$\left< \mathcal{J}\nabla\cdot\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \psi_j\right>_\Omega\f$ has been transformed into \f$\left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \psi_j\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \nabla \psi_j\right>_\Omega\f$ using integration by parts. 
 
 Where \f$\psi\f$ is a basis of test functions.
 
@@ -60,7 +60,7 @@ u(x, t) = \sum_i \tilde{u}_i \phi_i(x, t)
 The weak formulation is thus rewritten (Einstein summation convention applies):
 
 \f\[
-\frac{\partial \tilde{u}_i}{\partial t}\left<\phi_i, \psi_j\right> + \left< \mathcal{J}\frac{\partial \mathcal{H}(\tilde{u}_i\phi_i, t)}{\partial x}, \psi_j\right> - \left< \mathcal{J}\frac{\partial \mathcal{H}(\tilde{u}_i\phi_i, t)}{\partial (\nabla x)}, \nabla \psi_j\right> = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>
+\frac{\partial \tilde{u}_i}{\partial t}\left<\phi_i, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial u}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \psi_j\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \nabla \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
 \f\]
 
-\note As an example, the Poisson equation \f$\Delta\phi = \rho\f$ is associated to the Hamiltonian \f$\mathcal{H} = \frac{1}{2}|\nabla \phi|^2 - \rho\phi\f$.
+\note As an example, the Poisson equation \f$\Delta u = \rho\f$ is associated to the Hamiltonian \f$\mathcal{H} = \frac{1}{2}|\nabla u|^2 - \rho u\f$.

From 1c8532ebafcdbafc6d8fe7ffc0968287f507510b Mon Sep 17 00:00:00 2001
From: blegouix <stilynx51@gmail.com>
Date: Wed, 15 Jan 2025 14:48:19 +0100
Subject: [PATCH 06/13] wip

---
 docs/finite_element_module.md | 17 ++++++++++++++---
 1 file changed, 14 insertions(+), 3 deletions(-)

diff --git a/docs/finite_element_module.md b/docs/finite_element_module.md
index eb92eef9..38726a54 100644
--- a/docs/finite_element_module.md
+++ b/docs/finite_element_module.md
@@ -41,14 +41,20 @@ In which the variational derivative of the Hamiltonian can be expanded as:
 Leading to:
 
 \f\[
-\left<\frac{\partial u}{\partial t}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial u}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \psi_j\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \nabla \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
+\left<\frac{\partial u}{\partial t}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial u}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \psi_j\vec{n}\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \nabla \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
 \f\]
 
 
-\note \f$\left< \mathcal{J}\nabla\cdot\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \psi_j\right>_\Omega\f$ has been transformed into \f$\left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \psi_j\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \nabla \psi_j\right>_\Omega\f$ using integration by parts. 
+\note \f$\left< \mathcal{J}\nabla\cdot\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \psi_j\right>_\Omega\f$ has been transformed into \f$\left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \psi_j\vec{n}\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \nabla \psi_j\right>_\Omega\f$ using integration by parts. 
 
 Where \f$\psi\f$ is a basis of test functions.
 
+\note As an example, the heat equation \f$\frac{\partial u}{\partial t} + \Delta u = s\f$ is associated to the Hamiltonian \f$\mathcal{H} = \frac{1}{2}|\nabla u|^2\f$ and the heat source term \f$\mathcal{D} = s\f$, leading to:
+
+\note \f\[
+\left<\frac{\partial u}{\partial t}, \psi_j\right>_\Omega + \left< \mathcal{J}\nabla u, \psi_j\vec{n}\right>_{\partial\Omega} - \left< \mathcal{J}\nabla u, \nabla \psi_j\right>_\Omega = \left<s, \psi_j \right>_\Omega
+\f\]
+
 We decompose \f$u\f$ and \f$\psi_j\f$ in a discrete basis function \f$\phi\f$. 
 
 \f\[
@@ -63,4 +69,9 @@ The weak formulation is thus rewritten (Einstein summation convention applies):
 \frac{\partial \tilde{u}_i}{\partial t}\left<\phi_i, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial u}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \psi_j\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \nabla \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
 \f\]
 
-\note As an example, the Poisson equation \f$\Delta u = \rho\f$ is associated to the Hamiltonian \f$\mathcal{H} = \frac{1}{2}|\nabla u|^2 - \rho u\f$.
+\note As an example, the Poisson equation \f$\Delta u = \rho\f$ is associated to the Hamiltonian \f$\mathcal{H} = \frac{1}{2}|\nabla u|^2 - \rho u\f$, leading to:
+
+\note \f\[
+\frac{\partial \tilde{u}_i}{\partial t}\left<\phi_i, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial u}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \psi_j\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \nabla \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
+\f\]
+

From 193d81c779225f47a75360c1bfce914e6c51e2a1 Mon Sep 17 00:00:00 2001
From: blegouix <stilynx51@gmail.com>
Date: Wed, 15 Jan 2025 19:50:56 +0100
Subject: [PATCH 07/13] wip

---
 docs/finite_element_module.md | 37 +++++++++++++++++++++++++++++++++--
 1 file changed, 35 insertions(+), 2 deletions(-)

diff --git a/docs/finite_element_module.md b/docs/finite_element_module.md
index 38726a54..4b73e64a 100644
--- a/docs/finite_element_module.md
+++ b/docs/finite_element_module.md
@@ -69,9 +69,42 @@ The weak formulation is thus rewritten (Einstein summation convention applies):
 \frac{\partial \tilde{u}_i}{\partial t}\left<\phi_i, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial u}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \psi_j\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \nabla \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
 \f\]
 
-\note As an example, the Poisson equation \f$\Delta u = \rho\f$ is associated to the Hamiltonian \f$\mathcal{H} = \frac{1}{2}|\nabla u|^2 - \rho u\f$, leading to:
+To continue the computation, we need to assume a quadratic form for the Hamiltonian:
+
+\f\[
+\mathcal{H} \sim \frac{\alpha}{2} u^2 + \frac{\beta}{2} |\nabla u|^2 + \gamma u|\nabla u| + \delta u + \epsilon \nabla u + \zeta
+\f\]
+
+\attention If the Hamiltonian does not have this form, implementing a non-linear solver is required (it would rely on a quadratic approximation of \f$\mathcal{H}\f$, thus being able to solve the quadratic case is still necessary). 
+
+\f\[
+\frac{\partial \tilde{u}_i}{\partial t}\left<\phi_i, \psi_j\right>_\Omega + \tilde{u}_i\left(\left< \mathcal{J}\left(\alpha \phi_i+\gamma |\nabla \phi_i| + \delta\right), \psi_j\right>_\Omega + \left< \mathcal{J}\left(\beta |\nabla \phi_i|+ \gamma \phi + \epsilon\right), \psi_j\right>_{\partial\Omega}\\ - \left< \mathcal{J}\left(\beta |\nabla \phi_i|+ \gamma \phi_i + \epsilon\right), \nabla \psi_j\right>_\Omega\right) = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
+\f\]
+
+Or in other words:
+
+\f[
+M\frac{\partial \tilde{u}}{\partial t} + K\tilde{u} = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
+\f]
+
+With:
+
+\f\[
+\left\{\begin{matrix}
+M = \left<\phi_i, \psi_j\right>_\Omega\\
+K = \left< \mathcal{J}\left(\alpha \phi_i+\gamma |\nabla \phi_i| + \delta\right), \psi_j\right>_\Omega\\ + \left< \mathcal{J}\left(\beta |\nabla \phi_i|+ \gamma \phi + \epsilon\right), \psi_j\right>_{\partial\Omega}\\ - \left< \mathcal{J}\left(\beta |\nabla \phi_i|+ \gamma \phi_i + \epsilon\right), \nabla \psi_j\right>_\Omega
+\end{matrix}\right.
+\f\]
+
+(\f$M\f$ is called the <em>mass matrix</em> and \f$S\f$ the <em>stiffness matrix</em>)
+
+\note As an example, the heat equation is associated to: 
 
 \note \f\[
-\frac{\partial \tilde{u}_i}{\partial t}\left<\phi_i, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial u}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \psi_j\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \nabla \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
+\left\{\begin{matrix}
+M = \left<\phi_i, \psi_j\right>_\Omega\\
+K = \left< \mathcal{J}\nabla \phi_i, \psi_j\right>_{\partial\Omega}\\ - \left< \mathcal{J} |\nabla \phi_i|, \nabla \psi_j\right>_\Omega
+\end{matrix}\right.
 \f\]
 
+## How is Finite Element Method implemented in Similie ?

From 7c5a0cc8e96858137a1cc17d8d7a84dace65a53d Mon Sep 17 00:00:00 2001
From: blegouix <stilynx51@gmail.com>
Date: Wed, 15 Jan 2025 20:10:31 +0100
Subject: [PATCH 08/13] wip

---
 docs/finite_element_module.md | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/docs/finite_element_module.md b/docs/finite_element_module.md
index 4b73e64a..174ee952 100644
--- a/docs/finite_element_module.md
+++ b/docs/finite_element_module.md
@@ -18,7 +18,7 @@ Where \f$D^k\f$ is a \f$k-\f$order differential operator. Any partial differenti
 \frac{\partial u}{\partial t} + \{\mathcal{H}, u\} = \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) 
 \f\]
 
-Where \f$\{\mathcal{H}, u\}\f$ is the Poisson bracket between the Hamiltonian \f$\mathcal{H}(u, \nabla{u}, t)\f$ of the system and \f$u\f$. It can always be written \f$\{\mathcal{H}, u\} = \mathcal{J}\frac{\delta \mathcal{H}}{\delta u}\f$ in term of the symplectic matrix \f$\mathcal{J}\f$ (where \f$\frac{\delta}{\delta u}\f$ is a [variational derivative](https://en.wikipedia.org/wiki/Functional_derivative)). The Hamiltonian part of the PDE is associated to conserved quantities such as mass, momentum or energy.
+Where \f$\{\mathcal{H}, u\}\f$ is the Poisson bracket between the Hamiltonian \f$\mathcal{H}(u, \vec{\nabla}u, t)\f$ of the system and \f$u\f$. It can always be written \f$\{\mathcal{H}, u\} = \mathcal{J}\frac{\delta \mathcal{H}}{\delta u}\f$ in term of the symplectic matrix \f$\mathcal{J}\f$ (where \f$\frac{\delta}{\delta u}\f$ is a [variational derivative](https://en.wikipedia.org/wiki/Functional_derivative)). The Hamiltonian part of the PDE is associated to conserved quantities such as mass, momentum or energy.
 
 \note If a PDE involves a second-order temporal derivative \f$\frac{\partial^2 u}{\partial t^2}\f$ (like in wave equation) we can always define \f$U = (\frac{\partial u}{\partial t}, u)\f$ and write the PDE in term of  \f$U\f$.
 
@@ -35,24 +35,24 @@ Finite Element Method aims to discretize PDEs to produce numerical schemes able
 In which the variational derivative of the Hamiltonian can be expanded as: 
 
 \f\[
- \frac{\delta \mathcal{H}}{\delta u} = \frac{\partial \mathcal{H}}{\partial u} - \nabla\cdot \frac{\partial \mathcal{H}}{\partial (\nabla u)}
+ \frac{\delta \mathcal{H}}{\delta u} = \frac{\partial \mathcal{H}}{\partial u} - \vec{\nabla}\cdot \frac{\partial \mathcal{H}}{\partial (\vec{\nabla} u)}
 \f\]
 
 Leading to:
 
 \f\[
-\left<\frac{\partial u}{\partial t}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial u}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \psi_j\vec{n}\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \nabla \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
+\left<\frac{\partial u}{\partial t}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial u}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\vec{\nabla} u)}, \psi_j\vec{n}\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\vec{\nabla} u)}, \vec{\nabla} \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
 \f\]
 
 
-\note \f$\left< \mathcal{J}\nabla\cdot\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \psi_j\right>_\Omega\f$ has been transformed into \f$\left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \psi_j\vec{n}\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \nabla \psi_j\right>_\Omega\f$ using integration by parts. 
+\note \f$\left< \mathcal{J}\vec{\nabla}\cdot\frac{\partial \mathcal{H}}{\partial (\vec{\nabla} u)}, \psi_j\right>_\Omega\f$ has been transformed into \f$\left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\vec{\nabla} u)}, \psi_j\vec{n}\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\vec{\nabla} u)}, \vec{\nabla} \psi_j\right>_\Omega\f$ using integration by parts. 
 
 Where \f$\psi\f$ is a basis of test functions.
 
-\note As an example, the heat equation \f$\frac{\partial u}{\partial t} + \Delta u = s\f$ is associated to the Hamiltonian \f$\mathcal{H} = \frac{1}{2}|\nabla u|^2\f$ and the heat source term \f$\mathcal{D} = s\f$, leading to:
+\note As an example, the heat equation \f$\frac{\partial u}{\partial t} + \Delta u = s\f$ is associated to the Hamiltonian \f$\mathcal{H} = \frac{1}{2}|\vec{\nabla} u|^2\f$ and the heat source term \f$\mathcal{D} = s\f$, leading to:
 
 \note \f\[
-\left<\frac{\partial u}{\partial t}, \psi_j\right>_\Omega + \left< \mathcal{J}\nabla u, \psi_j\vec{n}\right>_{\partial\Omega} - \left< \mathcal{J}\nabla u, \nabla \psi_j\right>_\Omega = \left<s, \psi_j \right>_\Omega
+\left<\frac{\partial u}{\partial t}, \psi_j\right>_\Omega + \left< \mathcal{J}\vec{\nabla} u, \psi_j\vec{n}\right>_{\partial\Omega} - \left< \mathcal{J}\vec{\nabla} u, \vec{\nabla} \psi_j\right>_\Omega = \left<s, \psi_j \right>_\Omega
 \f\]
 
 We decompose \f$u\f$ and \f$\psi_j\f$ in a discrete basis function \f$\phi\f$. 
@@ -66,19 +66,19 @@ u(x, t) = \sum_i \tilde{u}_i \phi_i(x, t)
 The weak formulation is thus rewritten (Einstein summation convention applies):
 
 \f\[
-\frac{\partial \tilde{u}_i}{\partial t}\left<\phi_i, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial u}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \psi_j\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\nabla u)}, \nabla \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
+\frac{\partial \tilde{u}_i}{\partial t}\left<\phi_i, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial u}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\vec{\nabla} u)}, \psi_j\vec{n}\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\vec{\nabla} u)}, \vec{\nabla} \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
 \f\]
 
 To continue the computation, we need to assume a quadratic form for the Hamiltonian:
 
 \f\[
-\mathcal{H} \sim \frac{\alpha}{2} u^2 + \frac{\beta}{2} |\nabla u|^2 + \gamma u|\nabla u| + \delta u + \epsilon \nabla u + \zeta
+\mathcal{H} \sim \frac{\alpha}{2} u^2 + \frac{\beta}{2} ||\vec{\nabla} u||^2 + \gamma u||\vec{\nabla} u|| + \delta u + \epsilon ||\vec{\nabla} u|| + \zeta
 \f\]
 
 \attention If the Hamiltonian does not have this form, implementing a non-linear solver is required (it would rely on a quadratic approximation of \f$\mathcal{H}\f$, thus being able to solve the quadratic case is still necessary). 
 
 \f\[
-\frac{\partial \tilde{u}_i}{\partial t}\left<\phi_i, \psi_j\right>_\Omega + \tilde{u}_i\left(\left< \mathcal{J}\left(\alpha \phi_i+\gamma |\nabla \phi_i| + \delta\right), \psi_j\right>_\Omega + \left< \mathcal{J}\left(\beta |\nabla \phi_i|+ \gamma \phi + \epsilon\right), \psi_j\right>_{\partial\Omega}\\ - \left< \mathcal{J}\left(\beta |\nabla \phi_i|+ \gamma \phi_i + \epsilon\right), \nabla \psi_j\right>_\Omega\right) = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
+\frac{\partial \tilde{u}_i}{\partial t}\left<\phi_i, \psi_j\right>_\Omega + \tilde{u}_i\left(\left< \mathcal{J}\left(\alpha \phi_i+\gamma ||\vec{\nabla} \phi_i|| + \delta\right), \psi_j\right>_\Omega + \left< \mathcal{J}\left(\beta ||\vec{\nabla} \phi_i||+ \gamma \phi + \epsilon\right), \psi_j\vec{n}\right>_{\partial\Omega}\\ - \left< \mathcal{J}\left(\beta ||\vec{\nabla} \phi_i||+ \gamma \phi_i + \epsilon\right), \vec{\nabla} \psi_j\right>_\Omega\right) = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
 \f\]
 
 Or in other words:
@@ -92,7 +92,7 @@ With:
 \f\[
 \left\{\begin{matrix}
 M = \left<\phi_i, \psi_j\right>_\Omega\\
-K = \left< \mathcal{J}\left(\alpha \phi_i+\gamma |\nabla \phi_i| + \delta\right), \psi_j\right>_\Omega\\ + \left< \mathcal{J}\left(\beta |\nabla \phi_i|+ \gamma \phi + \epsilon\right), \psi_j\right>_{\partial\Omega}\\ - \left< \mathcal{J}\left(\beta |\nabla \phi_i|+ \gamma \phi_i + \epsilon\right), \nabla \psi_j\right>_\Omega
+K = \left< \mathcal{J}\left(\alpha \phi_i+\gamma ||\vec{\nabla} \phi_i|| + \delta\right), \psi_j\right>_\Omega\\ + \left< \mathcal{J}\left(\beta ||\vec{\nabla} \phi_i||+ \gamma \phi + \epsilon\right), \psi_j\vec{n}\right>_{\partial\Omega}\\ - \left< \mathcal{J}\left(\beta ||\vec{\nabla} \phi_i||+ \gamma \phi_i + \epsilon\right), \vec{\nabla} \psi_j\right>_\Omega
 \end{matrix}\right.
 \f\]
 
@@ -103,7 +103,7 @@ K = \left< \mathcal{J}\left(\alpha \phi_i+\gamma |\nabla \phi_i| + \delta\right)
 \note \f\[
 \left\{\begin{matrix}
 M = \left<\phi_i, \psi_j\right>_\Omega\\
-K = \left< \mathcal{J}\nabla \phi_i, \psi_j\right>_{\partial\Omega}\\ - \left< \mathcal{J} |\nabla \phi_i|, \nabla \psi_j\right>_\Omega
+K = \left< \mathcal{J}\vec{\nabla} \phi_i, \psi_j\vec{n}\right>_{\partial\Omega}\\ - \left< \mathcal{J} |\vec{\nabla} \phi_i|, \vec{\nabla} \psi_j\right>_\Omega
 \end{matrix}\right.
 \f\]
 

From 306e6e1b467d5ea41a43aad150fad750b8aa91aa Mon Sep 17 00:00:00 2001
From: blegouix <stilynx51@gmail.com>
Date: Thu, 16 Jan 2025 15:39:40 +0100
Subject: [PATCH 09/13] wip

---
 docs/finite_element_module.md | 32 ++++++++++++++++----------------
 1 file changed, 16 insertions(+), 16 deletions(-)

diff --git a/docs/finite_element_module.md b/docs/finite_element_module.md
index 174ee952..bc4833f4 100644
--- a/docs/finite_element_module.md
+++ b/docs/finite_element_module.md
@@ -1,4 +1,4 @@
-# Finite Element module {#finite_element_module}
+# Finite element module {#finite_element_module}
 <!--
 SPDX-FileCopyrightText: 2025 Baptiste Legouix
 SPDX-License-Identifier: GPL-3.0-or-later
@@ -6,7 +6,7 @@ SPDX-License-Identifier: GPL-3.0-or-later
 
 ## Why partial differential equations are importants for physics ?
 
-Given a multidimensional function \f$u:x\rightarrow \; .\f$, any [partial differential equation](https://en.wikipedia.org/wiki/Partial_differential_equation) can be written:
+Given a differential form \f$u\f$ parametrized by coordinates \f$x\f$, any [partial differential equation](https://en.wikipedia.org/wiki/Partial_differential_equation) can be written:
 
 \f\[
 \mathcal{F}(D^k u, D^{k-1} u,... ,D u, u, x) = 0
@@ -18,7 +18,7 @@ Where \f$D^k\f$ is a \f$k-\f$order differential operator. Any partial differenti
 \frac{\partial u}{\partial t} + \{\mathcal{H}, u\} = \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) 
 \f\]
 
-Where \f$\{\mathcal{H}, u\}\f$ is the Poisson bracket between the Hamiltonian \f$\mathcal{H}(u, \vec{\nabla}u, t)\f$ of the system and \f$u\f$. It can always be written \f$\{\mathcal{H}, u\} = \mathcal{J}\frac{\delta \mathcal{H}}{\delta u}\f$ in term of the symplectic matrix \f$\mathcal{J}\f$ (where \f$\frac{\delta}{\delta u}\f$ is a [variational derivative](https://en.wikipedia.org/wiki/Functional_derivative)). The Hamiltonian part of the PDE is associated to conserved quantities such as mass, momentum or energy.
+Where \f$\{\mathcal{H}, u\}\f$ is the Poisson bracket between the Hamiltonian \f$\mathcal{H}(u, du, t)\f$ of the system and \f$u\f$. It can always be written \f$\{\mathcal{H}, u\} = \mathcal{J}\frac{\delta \mathcal{H}}{\delta u}\f$ in term of the symplectic matrix \f$\mathcal{J}\f$ (where \f$\frac{\delta}{\delta u}\f$ is a [variational derivative](https://en.wikipedia.org/wiki/Functional_derivative)). The Hamiltonian part of the PDE is associated to conserved quantities such as mass, momentum or energy.
 
 \note If a PDE involves a second-order temporal derivative \f$\frac{\partial^2 u}{\partial t^2}\f$ (like in wave equation) we can always define \f$U = (\frac{\partial u}{\partial t}, u)\f$ and write the PDE in term of  \f$U\f$.
 
@@ -26,7 +26,7 @@ Partial differential equations are the main mathematical tool to describe rules
 
 ## What is Finite Element Method ?
 
-Finite Element Method aims to discretize PDEs to produce numerical schemes able to solve them. It relies on the weak formulation of the PDEs (\f$\left<\cdot, \cdot\right>_\Omega\f$ is the inner product between functions on the domain \f$\Omega\f$):
+Finite Element Method aims to discretize PDEs to produce numerical schemes able to solve them. It relies on the weak formulation of the PDEs (\f$\left<\cdot, \cdot\right>_\Omega = \int_\Omega \cdot \wedge \star \cdot\;\f$ is the inner product between forms on the domain \f$\Omega\f$):
 
 \f\[
 \left<\frac{\partial u}{\partial t}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\delta \mathcal{H}}{\delta u}, \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
@@ -35,24 +35,23 @@ Finite Element Method aims to discretize PDEs to produce numerical schemes able
 In which the variational derivative of the Hamiltonian can be expanded as: 
 
 \f\[
- \frac{\delta \mathcal{H}}{\delta u} = \frac{\partial \mathcal{H}}{\partial u} - \vec{\nabla}\cdot \frac{\partial \mathcal{H}}{\partial (\vec{\nabla} u)}
+\frac{\delta \mathcal{H}}{\delta u} = \frac{\partial \mathcal{H}}{\partial u} - d\frac{\partial \mathcal{H}}{\partial (du)}
 \f\]
 
 Leading to:
 
 \f\[
-\left<\frac{\partial u}{\partial t}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial u}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\vec{\nabla} u)}, \psi_j\vec{n}\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\vec{\nabla} u)}, \vec{\nabla} \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
+\left<\frac{\partial u}{\partial t}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial u}, \psi_j\right>_\Omega + \int_{\partial\Omega} \mathrm{tr} \; \mathcal{J}\frac{\partial \mathcal{H}}{\partial (du)} \wedge \mathrm{tr} \star \psi_j + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (du)}, d \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
 \f\]
 
-
-\note \f$\left< \mathcal{J}\vec{\nabla}\cdot\frac{\partial \mathcal{H}}{\partial (\vec{\nabla} u)}, \psi_j\right>_\Omega\f$ has been transformed into \f$\left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\vec{\nabla} u)}, \psi_j\vec{n}\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\vec{\nabla} u)}, \vec{\nabla} \psi_j\right>_\Omega\f$ using integration by parts. 
+\note \f$\left< \mathcal{J}d\frac{\partial \mathcal{H}}{\partial (d u)}, \psi_j\right>_\Omega\f$ has been transformed into \f$\int_{\partial\Omega} \mathrm{tr} \; \mathcal{J}\frac{\partial \mathcal{H}}{\partial (du)} \wedge \mathrm{tr} \star \psi_j + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (du)}, d\psi_j\right>_\Omega\f$ using integration by parts. 
 
 Where \f$\psi\f$ is a basis of test functions.
 
-\note As an example, the heat equation \f$\frac{\partial u}{\partial t} + \Delta u = s\f$ is associated to the Hamiltonian \f$\mathcal{H} = \frac{1}{2}|\vec{\nabla} u|^2\f$ and the heat source term \f$\mathcal{D} = s\f$, leading to:
+\note As an example, the heat equation \f$\frac{\partial u}{\partial t} + \Delta u = s\f$ is associated to the Hamiltonian \f$\mathcal{H} = \frac{1}{2}\left<du,du\right>\f$ and the heat source term \f$\mathcal{D} = s\f$, leading to:
 
 \note \f\[
-\left<\frac{\partial u}{\partial t}, \psi_j\right>_\Omega + \left< \mathcal{J}\vec{\nabla} u, \psi_j\vec{n}\right>_{\partial\Omega} - \left< \mathcal{J}\vec{\nabla} u, \vec{\nabla} \psi_j\right>_\Omega = \left<s, \psi_j \right>_\Omega
+\left<\frac{\partial u}{\partial t}, \psi_j\right>_\Omega + \int_{\partial\Omega} \mathrm{tr} \; \mathcal{J}du \wedge \mathrm{tr} \star \psi_j + \left< \mathcal{J}du, d \psi_j\right>_\Omega = \left<s, \psi_j \right>_\Omega
 \f\]
 
 We decompose \f$u\f$ and \f$\psi_j\f$ in a discrete basis function \f$\phi\f$. 
@@ -66,19 +65,20 @@ u(x, t) = \sum_i \tilde{u}_i \phi_i(x, t)
 The weak formulation is thus rewritten (Einstein summation convention applies):
 
 \f\[
-\frac{\partial \tilde{u}_i}{\partial t}\left<\phi_i, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial u}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\vec{\nabla} u)}, \psi_j\vec{n}\right>_{\partial\Omega} - \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\vec{\nabla} u)}, \vec{\nabla} \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
+\frac{\partial \tilde{u}_i}{\partial t}\left<\phi_i, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\tilde{u}_i\phi_i)}, \psi_j\right>_\Omega + \int_{\partial\Omega} \mathrm{tr} \; \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\tilde{u}_id\phi_i)} \wedge \mathrm{tr} \star \psi_j\\ + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (\tilde{u}_id\phi_i)}, d \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
 \f\]
 
 To continue the computation, we need to assume a quadratic form for the Hamiltonian:
 
 \f\[
-\mathcal{H} \sim \frac{\alpha}{2} u^2 + \frac{\beta}{2} ||\vec{\nabla} u||^2 + \gamma u||\vec{\nabla} u|| + \delta u + \epsilon ||\vec{\nabla} u|| + \zeta
+\mathcal{H} = \frac{\alpha}{2} \left<u|u\right> + \frac{\beta}{2} \left<du|du\right> = \tilde{u}_i^2\left(\frac{\alpha}{2} \left<\phi_i|\phi_i\right> + \frac{\beta}{2} \left<d\phi_i|d\phi_i\right>\right)
+
 \f\]
 
-\attention If the Hamiltonian does not have this form, implementing a non-linear solver is required (it would rely on a quadratic approximation of \f$\mathcal{H}\f$, thus being able to solve the quadratic case is still necessary). 
+\important If the Hamiltonian does not have this form, implementing a non-linear solver is required (it would rely on a quadratic approximation of \f$\mathcal{H}\f$, thus being able to solve the quadratic case is still necessary). 
 
 \f\[
-\frac{\partial \tilde{u}_i}{\partial t}\left<\phi_i, \psi_j\right>_\Omega + \tilde{u}_i\left(\left< \mathcal{J}\left(\alpha \phi_i+\gamma ||\vec{\nabla} \phi_i|| + \delta\right), \psi_j\right>_\Omega + \left< \mathcal{J}\left(\beta ||\vec{\nabla} \phi_i||+ \gamma \phi + \epsilon\right), \psi_j\vec{n}\right>_{\partial\Omega}\\ - \left< \mathcal{J}\left(\beta ||\vec{\nabla} \phi_i||+ \gamma \phi_i + \epsilon\right), \vec{\nabla} \psi_j\right>_\Omega\right) = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
+\frac{\partial \tilde{u}_i}{\partial t}\left<\phi_i, \psi_j\right>_\Omega + \tilde{u}_i\left(\alpha\left< \mathcal{J}\phi_i, \psi_j\right>_\Omega + \beta \int_{\partial\Omega} \mathrm{tr} \; \mathcal{J}d\phi_i \wedge \mathrm{tr} \star \psi_j\\ + \beta \left< \mathcal{J}d\phi_i, d \psi_j\right>_\Omega\right) = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
 \f\]
 
 Or in other words:
@@ -92,7 +92,7 @@ With:
 \f\[
 \left\{\begin{matrix}
 M = \left<\phi_i, \psi_j\right>_\Omega\\
-K = \left< \mathcal{J}\left(\alpha \phi_i+\gamma ||\vec{\nabla} \phi_i|| + \delta\right), \psi_j\right>_\Omega\\ + \left< \mathcal{J}\left(\beta ||\vec{\nabla} \phi_i||+ \gamma \phi + \epsilon\right), \psi_j\vec{n}\right>_{\partial\Omega}\\ - \left< \mathcal{J}\left(\beta ||\vec{\nabla} \phi_i||+ \gamma \phi_i + \epsilon\right), \vec{\nabla} \psi_j\right>_\Omega
+K = \alpha\left< \mathcal{J}\phi_i, \psi_j\right>_\Omega + \beta \int_{\partial\Omega} \mathrm{tr} \; \mathcal{J}d\phi_i \wedge \mathrm{tr} \star \psi_j + \beta \left< \mathcal{J}d\phi_i, d \psi_j\right>_\Omega
 \end{matrix}\right.
 \f\]
 
@@ -103,7 +103,7 @@ K = \left< \mathcal{J}\left(\alpha \phi_i+\gamma ||\vec{\nabla} \phi_i|| + \delt
 \note \f\[
 \left\{\begin{matrix}
 M = \left<\phi_i, \psi_j\right>_\Omega\\
-K = \left< \mathcal{J}\vec{\nabla} \phi_i, \psi_j\vec{n}\right>_{\partial\Omega}\\ - \left< \mathcal{J} |\vec{\nabla} \phi_i|, \vec{\nabla} \psi_j\right>_\Omega
+K = \int_{\partial\Omega} \mathrm{tr} \; \mathcal{J}d\phi_i \wedge \mathrm{tr} \star \psi_j + \left< \mathcal{J}d\phi_i, d \psi_j\right>_\Omega
 \end{matrix}\right.
 \f\]
 

From 84bc68b058e4a6a9ec93cc935a7c14aa0b6a5ab3 Mon Sep 17 00:00:00 2001
From: blegouix <stilynx51@gmail.com>
Date: Thu, 16 Jan 2025 15:57:22 +0100
Subject: [PATCH 10/13] wip

---
 docs/finite_element_module.md | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/docs/finite_element_module.md b/docs/finite_element_module.md
index bc4833f4..9f317a31 100644
--- a/docs/finite_element_module.md
+++ b/docs/finite_element_module.md
@@ -32,10 +32,10 @@ Finite Element Method aims to discretize PDEs to produce numerical schemes able
 \left<\frac{\partial u}{\partial t}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\delta \mathcal{H}}{\delta u}, \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
 \f\]
 
-In which the variational derivative of the Hamiltonian can be expanded as: 
+In which the variational derivative of the Hamiltonian can be expanded as (\f$d\f$ is the exterior derivative and \f$\delta\f$ the codifferential): 
 
 \f\[
-\frac{\delta \mathcal{H}}{\delta u} = \frac{\partial \mathcal{H}}{\partial u} - d\frac{\partial \mathcal{H}}{\partial (du)}
+\frac{\delta \mathcal{H}}{\delta u} = \frac{\partial \mathcal{H}}{\partial u} - \delta\frac{\partial \mathcal{H}}{\partial (du)}
 \f\]
 
 Leading to:
@@ -44,7 +44,7 @@ Leading to:
 \left<\frac{\partial u}{\partial t}, \psi_j\right>_\Omega + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial u}, \psi_j\right>_\Omega + \int_{\partial\Omega} \mathrm{tr} \; \mathcal{J}\frac{\partial \mathcal{H}}{\partial (du)} \wedge \mathrm{tr} \star \psi_j + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (du)}, d \psi_j\right>_\Omega = \left< \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x), \psi_j \right>_\Omega
 \f\]
 
-\note \f$\left< \mathcal{J}d\frac{\partial \mathcal{H}}{\partial (d u)}, \psi_j\right>_\Omega\f$ has been transformed into \f$\int_{\partial\Omega} \mathrm{tr} \; \mathcal{J}\frac{\partial \mathcal{H}}{\partial (du)} \wedge \mathrm{tr} \star \psi_j + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (du)}, d\psi_j\right>_\Omega\f$ using integration by parts. 
+\note \f$\left< \mathcal{J}\delta\frac{\partial \mathcal{H}}{\partial (d u)}, \psi_j\right>_\Omega\f$ has been transformed into \f$\int_{\partial\Omega} \mathrm{tr} \; \mathcal{J}\frac{\partial \mathcal{H}}{\partial (du)} \wedge \mathrm{tr} \star \psi_j + \left< \mathcal{J}\frac{\partial \mathcal{H}}{\partial (du)}, d\psi_j\right>_\Omega\f$ using integration by parts. 
 
 Where \f$\psi\f$ is a basis of test functions.
 

From 2a2011eb4d5a3f172d5ce430d91345b39fb8751d Mon Sep 17 00:00:00 2001
From: blegouix <stilynx51@gmail.com>
Date: Mon, 20 Jan 2025 16:11:42 +0100
Subject: [PATCH 11/13] wip

---
 docs/finite_element_module.md | 60 ++++++++++++++++++++++++++++++++---
 1 file changed, 55 insertions(+), 5 deletions(-)

diff --git a/docs/finite_element_module.md b/docs/finite_element_module.md
index 9f317a31..d813e710 100644
--- a/docs/finite_element_module.md
+++ b/docs/finite_element_module.md
@@ -6,19 +6,69 @@ SPDX-License-Identifier: GPL-3.0-or-later
 
 ## Why partial differential equations are importants for physics ?
 
-Given a differential form \f$u\f$ parametrized by coordinates \f$x\f$, any [partial differential equation](https://en.wikipedia.org/wiki/Partial_differential_equation) can be written:
+### Splitting between Hamiltonian and non-Hamiltonian parts 
+
+Given a differential form \f$u\f$ parametrized by coordinates \f$q^i\f$ and conjugate momenta \f$p_i\f$, any [partial differential equation](https://en.wikipedia.org/wiki/Partial_differential_equation) can be written:
+
+\f\[
+\mathcal{F}(D^k u, D^{k-1} u,... ,D u, u, p, q) = 0
+\f\]
+
+Where \f$D^k\f$ is a \f$k-\f$order differential operator. Defining the [conjugate momentum field](https://en.wikipedia.org/wiki/Hamiltonian_field_theory) \f$v\f$ of \f$u\f$, let's prove that any partial differential equation can be splitted into a Hamiltonian and non-Hamiltonian parts:
+
+\f\[
+\iota_X\Omega = d\mathcal{H}(u, v, q, p)  + \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, q, p) 
+\f\]
+
+Where the vector field \f$X\f$ is the infinitesimal generator of the [vector flow](https://en.wikipedia.org/wiki/Vector_flow) describing the dynamic of the field \f$u\f$. In other words, \f$u\f$ is transported along the integral curves \f$\gamma\f$ such that:
+
+\f\[
+\frac{d\gamma}{dt} = X(\gamma(t))
+\f\]
+
+Let's define the [Hamiltonian vector field](https://en.wikipedia.org/wiki/Hamiltonian_vector_field) \f$X_\mathcal{H}\f$ associated to the [Hamiltonian](https://en.wikipedia.org/wiki/Hamiltonian_mechanics) \f$\mathcal{H}\f$ such that:
+
+\f\[
+\iota_{X_\mathcal{H}} \Omega = d\mathcal{H}
+\f\]
+ 
+Where \f$\iota_{X_\mathcal{H}} \Omega\f$ is the interior product of the Hamiltonian vector field \f$X_\mathcal{H}\f$ with the inverse \f$\Omega\f$ of the [symplectic \f$2-\f$form](https://en.wikipedia.org/wiki/Symplectic_vector_space) \f$\omega = dp_i\wedge dq^i\f$. Such a field is [symplectic](https://en.wikipedia.org/wiki/Symplectic_vector_field), that is to say it obeys the [Liouville equation](https://en.wikipedia.org/wiki/Liouville%27s_theorem_(Hamiltonian)):
 
 \f\[
-\mathcal{F}(D^k u, D^{k-1} u,... ,D u, u, x) = 0
+\mathcal{L}_{X_\mathcal{H}}\Omega = \iota_{X_\mathcal{H}} d\Omega + d(\iota_{X_\mathcal{H}}\Omega) = \iota_{X_\mathcal{H}} d\Omega - dd\mathcal{H} = 0
 \f\]
 
-Where \f$D^k\f$ is a \f$k-\f$order differential operator. Any partial differential equation can be splitted into a Hamiltonian and non-Hamiltonian parts:
+Because \f$\Omega\f$ is a closed form and [Poincaré lemma](https://en.wikipedia.org/wiki/Poincar%C3%A9_lemma) applies. This means \f$\Omega\f$ is preserved along the flow. In other words, along an integral curve \f$\gamma_\mathcal{H}\f$ defined such that:
+
+\f\[
+\frac{d\gamma_\mathcal{H}}{dt} = X_\mathcal{H}(\gamma(t))
+\f\]
+
+The symplectic form \f$\Omega\f$ is preserved.
+
+By defining the non-Hamiltonian part of the vector field \f$X_\mathcal{D} = X - X_\mathcal{D}\f$, obeying:
+
+\f\[
+\iota_{X_\mathcal{D}}\Omega = \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, q, p) 
+\f\]
+
+The linearity of the interior product immediatly leads to:
+
+\f\[
+\iota_X\Omega = d\mathcal{H}(u, v, q, p)  + \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, q, p) 
+\f\]
+
+
+### Splitting between Hamiltonian and non-Hamiltonian parts 
+
+
+Any partial differential equation can be splitted into a Hamiltonian and non-Hamiltonian parts:
 
 \f\[
-\frac{\partial u}{\partial t} + \{\mathcal{H}, u\} = \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) 
+\iota_X \Omega = d\mathcal{H}(D^k u, D^{k-1} u,... ,D u, u, x)  + \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, x) 
 \f\]
 
-Where \f$\{\mathcal{H}, u\}\f$ is the Poisson bracket between the Hamiltonian \f$\mathcal{H}(u, du, t)\f$ of the system and \f$u\f$. It can always be written \f$\{\mathcal{H}, u\} = \mathcal{J}\frac{\delta \mathcal{H}}{\delta u}\f$ in term of the symplectic matrix \f$\mathcal{J}\f$ (where \f$\frac{\delta}{\delta u}\f$ is a [variational derivative](https://en.wikipedia.org/wiki/Functional_derivative)). The Hamiltonian part of the PDE is associated to conserved quantities such as mass, momentum or energy.
+It can always be written \f$\{\mathcal{H}, u\} = \mathcal{J}\frac{\delta \mathcal{H}}{\delta u}\f$ in term of the symplectic matrix \f$\mathcal{J}\f$ (where \f$\frac{\delta}{\delta u}\f$ is a [variational derivative](https://en.wikipedia.org/wiki/Functional_derivative)). The Hamiltonian part of the PDE is associated to conserved quantities such as mass, momentum or energy.
 
 \note If a PDE involves a second-order temporal derivative \f$\frac{\partial^2 u}{\partial t^2}\f$ (like in wave equation) we can always define \f$U = (\frac{\partial u}{\partial t}, u)\f$ and write the PDE in term of  \f$U\f$.
 

From 09e3823c9449931048f49e03cbe32bbce7b8083a Mon Sep 17 00:00:00 2001
From: blegouix <stilynx51@gmail.com>
Date: Mon, 20 Jan 2025 16:39:11 +0100
Subject: [PATCH 12/13] wip

---
 docs/finite_element_module.md | 17 +++++++++++++----
 1 file changed, 13 insertions(+), 4 deletions(-)

diff --git a/docs/finite_element_module.md b/docs/finite_element_module.md
index d813e710..06afefe3 100644
--- a/docs/finite_element_module.md
+++ b/docs/finite_element_module.md
@@ -26,11 +26,20 @@ Where the vector field \f$X\f$ is the infinitesimal generator of the [vector flo
 \frac{d\gamma}{dt} = X(\gamma(t))
 \f\]
 
-Let's define the [Hamiltonian vector field](https://en.wikipedia.org/wiki/Hamiltonian_vector_field) \f$X_\mathcal{H}\f$ associated to the [Hamiltonian](https://en.wikipedia.org/wiki/Hamiltonian_mechanics) \f$\mathcal{H}\f$ such that:
+Let's define the [Hamiltonian vector field](https://en.wikipedia.org/wiki/Hamiltonian_vector_field) \f$X_\mathcal{H}\f$ associated to the [Hamiltonian](https://en.wikipedia.org/wiki/Hamiltonian_mechanics) \f$\mathcal{H}\f$:
 
 \f\[
 \iota_{X_\mathcal{H}} \Omega = d\mathcal{H}
 \f\]
+
+\note This key definition is the exterior calculus counterpart of the [Hamilton's field equations](https://en.wikipedia.org/wiki/Hamiltonian_field_theory), usually written in term of the variational derivative \f$\frac{\delta}{\delta f}\f$:
+
+\note \f\[
+\left\{\begin{matrix}
+\frac{\partial u}{\partial t} = \frac{\delta \mathcal{H}}{\delta v}\\
+\frac{\partial v}{\partial t} = -\frac{\delta \mathcal{H}}{\delta u}
+\end{matrix}\right.
+\f\]
  
 Where \f$\iota_{X_\mathcal{H}} \Omega\f$ is the interior product of the Hamiltonian vector field \f$X_\mathcal{H}\f$ with the inverse \f$\Omega\f$ of the [symplectic \f$2-\f$form](https://en.wikipedia.org/wiki/Symplectic_vector_space) \f$\omega = dp_i\wedge dq^i\f$. Such a field is [symplectic](https://en.wikipedia.org/wiki/Symplectic_vector_field), that is to say it obeys the [Liouville equation](https://en.wikipedia.org/wiki/Liouville%27s_theorem_(Hamiltonian)):
 
@@ -44,15 +53,15 @@ Because \f$\Omega\f$ is a closed form and [Poincaré lemma](https://en.wikipedia
 \frac{d\gamma_\mathcal{H}}{dt} = X_\mathcal{H}(\gamma(t))
 \f\]
 
-The symplectic form \f$\Omega\f$ is preserved.
+The symplectic form \f$\Omega\f$ is invariant.
 
-By defining the non-Hamiltonian part of the vector field \f$X_\mathcal{D} = X - X_\mathcal{D}\f$, obeying:
+By defining the non-Hamiltonian part of the vector field \f$X_\mathcal{D} = X - X_\mathcal{H}\f$, obeying:
 
 \f\[
 \iota_{X_\mathcal{D}}\Omega = \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, q, p) 
 \f\]
 
-The linearity of the interior product immediatly leads to:
+The linearity of the interior product immediately leads to:
 
 \f\[
 \iota_X\Omega = d\mathcal{H}(u, v, q, p)  + \mathcal{D}(D^k u, D^{k-1} u,... ,D u, u, q, p) 

From 4ef73107149eb56af487193266479a83552b7d35 Mon Sep 17 00:00:00 2001
From: Baptiste Legouix <stilynx51@gmail.com>
Date: Fri, 31 Jan 2025 14:39:19 +0100
Subject: [PATCH 13/13] Update docs/finite_element_module.md

---
 docs/finite_element_module.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/finite_element_module.md b/docs/finite_element_module.md
index 06afefe3..d0a684f6 100644
--- a/docs/finite_element_module.md
+++ b/docs/finite_element_module.md
@@ -1,7 +1,7 @@
 # Finite element module {#finite_element_module}
 <!--
 SPDX-FileCopyrightText: 2025 Baptiste Legouix
-SPDX-License-Identifier: GPL-3.0-or-later
+SPDX-License-Identifier: MIT
 -->
 
 ## Why partial differential equations are importants for physics ?