Hybrid Discontinue Galerkin (HDG)
This part was treated at the end of the internship, and it was harder for me to understand. There may be some imprecision left on this page. For safer information about HDG method, you can report to the corresponding page in Feel++ mathematics [FeelppMath], or in the article [HDG2020]. 
1. HDG formulation
The HDG method uses discontinuous elements (contrary to usual, the elements are continuous in standard Galerkin methods).
It consists in rewrite a system of equations in a firstorder system.
In the HDG method, we have an unknown, the potential field (for example the temperature) on which we add the associated flux, who becomes unknown too. It allows formulating the problem with a better cost, by introducing the trace.
The HDG method has some attractive features :

It gives an optimal approximation for the potential fiels and the flux

It requires less degrees of freedom than usual Discontinuous Galerkin methods for a comparable accuracy

It allow local postprocessin to obtain better approximations and conservation properties.
See [HDG2020] for more complete details about this method.
2. Application to Bestest heat equation
Those are the equtation, tken from here :
Where the unknowns of the equation are :

\(T_w\) the temperature in the wall

\(T_r\) the temperature in the room (supposed to be constant in space over the room)
We introduce the thermal flux :
Replacing this in the equation \(\ref{hw}\), we get the mixte problem, with firstorder equations :
Let’s partition the border \(\Gamma\) of \(\Omega\) in three disjoint subsets : \(\Gamma = \Gamma_D \cup \Gamma_N \cup \Gamma_{ibc}\)
The quantities \(\q_w\) and \(T_w\) are subject to those boundary conditions :

\(T_w=0\) on \(\Gamma_D\) (Dirichlet condition)

\(\q_w\cdot\n = 0\) on \(\Gamma_N\) (Neumann condition)

\(T_w\) is a constant (unknown) on \(\Gamma_{ibc}\) and \(\displaystyle\int_{\Gamma_{ibc}}\q_w\cdot\n=I_\text{target}\) (Integral boundary condition, or IBC)
With such a condition, \(T_w\) is constant but unknown. It makes sense because the temperature in the wall is supposed to be constant in space, but can vary in time. It it what tells the EDO \(\ref{hr}\)
The value of \(\int_{\Gamma_{ibc}} \q_w \cdot \n\) is equal to the quantity of heat exchanged with the air, so :
We can deduce the value of \(I_\text{target}\) which is therefore equal to \(\displaystyle\int_{\Gamma_{ibc}} h(T_wT_r)\), because \(T_w\) is supposed to be constant over \(\Gamma_{ibc}\).
Let \(\bar{\varphi}\in H^{1/2}(\Gamma)\) such that \(\bar{\varphi}_{\Gamma_D}=0\) and \(\bar{\varphi}_{\Gamma_{ibc}}=1\). The variational problem of our equation system is :
Find \(\q_w\in H(div,\Omega)\), \(T_w\in L^2(\Omega)\) and \(\widehat{T_w} \in \mathrm{span} \left<\bar{\varphi}\right> \oplus H^{1/2}_{00}(\Gamma_N)\) such as for all \(\v\in H(div;\Omega)\), \(w\in L^2(\Omega)\), \(\mu\in\mathrm{span}\left<\bar{\varphi}\right>\oplus H^{1/2}_{00}(\Gamma_N)\) we have :
We use those spaces :

We make the usual approximation \(\dfrac{\partial T_w}{\partial t} = \dfrac{T_wT_w^\text{prev}}{\Delta t}\), where \(T_w^\text{prev}\) is the temperature from the previous step of the simulation. 
We can show (cf. [HDG2020]) that if \(g_N\in L^2(\Gamma_N)\), this problem admits a unique solution \((\q_w,T_w,\widehat{T}_w)\), satisfying the conditions given earlier.
The trace \(\widehat{T}_w\) of \(T_w\) allows reducing the cost of computation. We can show (cf. [Sala2019]) that the HDG method gives the best approximation error on potential and the flux.
3. Discretization
Let \(\T_h\) be a triangulation of \(\Omega\). We denote by \(\F_h\) the set of all faces of \(\T_h\), \(\F_h^\Gamma\) the set of faces belonging to \(\Gamma\) (and a similar notation for the subsets of \(\Gamma\)), and \(\F_h^0\) the set of the faces belonging to the interior of \(\Omega\).
If \(\q\in H(div,K)\), then its normal trace on the boundary of \(K\), \(\q\cdot\n\) belongs to \(H^{1/2}(K)\). In the following, we note \(\q^K := \q_K\). Let’s suppose that \(\q\in \left(L^2(\Omega)\right)^d\), and \(\q^K\in H(div,K)\) for all \(K\in\T_h\). Let \(F := \partial K_1\cap\partial K_2\) the border between two elements \(K_1,K_2\in \T_h\), this face belonging to \(\F_h^0\). We define the jump of the normal trace of \(\q\) across \(F\) as :
We can prove (cf [HDG2020]) that
We introduce those three spaces :

\(V_h = \displaystyle\prod_{K\in \T_h}\mathbf{P}_k(K)\)

\(W_h = \displaystyle\prod_{K\in \T_h} P_k(K)\)

\(M_h = \mathrm{span}\left<\varphi^*\right> \oplus \displaystyle\prod_{F\in\F_h^0 \cup\F_h^{\Gamma_N}} P_k(F)\)
With \(\varphi^*\in L^2(\F_h)\) such as \(\forall F\in \F_h^{\Gamma_{ibc}}, \varphi^*_{F}=1\) and \(\forall F\in \F\setminus\F_h^{\Gamma_{ibc}}, \varphi^*_{F}=0\). And \(\mathbf{P}_k(K)=(\P_k(K))^d\).
We see that the spaces are richer because the functions belonging to them are in general discontinuous. The functions of \(M_h\) play the role of connectors between adjacent elements, in order to reduce the time of calculation.
We define the numerical normal flux on \(\partial K\) :
Where \(\tau_K\) is a nonnegative function that can vary on \(\partial K\), and \(\tau_K > 0\) on at least one face of \(\partial K\). This parameter plays on the continuité of the potential, and to get a stabilized digital flow, we may have to scale this parameter or the flux. It is called stabilization parameter.
The discrete variational problem is therefore : find \(\q_{w,h}\in V_h\), \(T_{w,h}\in W_h\) and \(\widehat{T}_{w,h}\in M_h\) such as \(\forall \v_h\in V_h, \forall w_h\in W_h, \forall \mu_h\in M_h\)
4. Static condensation
The discrete equation \(\ref{eq:pb_var_disc}\) hlod in the interior of each \(K\in \T_h\) and can b solved for each \(K\) to eliminate \(\q_w^K\) and \(T_{w,h}^K\) in favor of the variable \(\widehat{T}_{w,h}^K\). This method of elimination is called static condensation.
With the equation \(\ref{eq:flux}\), it is possible to express the normal numerical flux as a function of \(\widehat{T}_{w,h}^K\).
We are now going to recast the system \(\ref{eq:flux},\ref{eq:pb_var_disc}\) to a global linear system where only the trace of the solution on the boundaries of the mesh appears. After solving the global system, the uknowns wan be recovered locally, in parallel.
First of all, we write the space for the numerical trace \(\widehat{T}_{w,h}\) as a direct sum :
with :

\(\widetilde{M}_h = \left\{\mu\in L^2(\F_h) : \forall F\in\F_h^0\cup F_h^{\Gamma_N} \mu_F\in\P_k(F) \text{ and } \forall F\in \F_h\setminus (\F_h^0\cup\F_h^{\Gamma_N})\mu_F = 0\right\}\). This is the sapce the the standard trace of the elements we find in usual HDG methods.

\(M_h^* = \left\{\mu\in L^2(\F_h) : \mu_{\Gamma_{ubc}}\in\RR \text{ and } \forall F\in \F_h\setminus \F_h^{\Gamma_{ibc}}\mu_F = 0\right\}\). This space handles the condition \(T_w\) is constant over \(\Gamma_{ibc}\).
Let \(\lambda_{1,h} = \widehat{T}_{w,h}_{\widetilde{M}_h}\) and \(\lambda_{2,h} = \widehat{T}_{w,h}_{M_h^*}\). Making integrations by parts, the formulation \(\ref{eq:flux},\ref{eq:pb_var_disc}\) can be rewritten as :
for all \(\v_h\in V_h, w_h\in W_h, \mu_{1,h}\in\widetilde{M}_h, \mu_{2,h}\in M_h^*\). We can notice that \(M_h^*\) has a dimension 1, so \(M_h^*\simeq\RR\). Because of it, the fourth equation of \(\ref{eq:hdg_ibc}\) is a single scalar equation, enforcing the integral boundary condition.
References

[CEN2007] EN 15026, Hygrothermal performance of building components and building elements  Assessment of moisture transfer by numerical simulation, CEN, 2007.

[FeelppMath] Feel++ Mathematics

[FppPicard] Feel, Non linear problems on http://docs.feelpp.org/math/fem/nonlinear/#_picard_strategy[Feel Mathematics]

[HDG2020] A HDG method for elliptic problems with integral boundary condition: Theory and Applications, Silvia Bertoluzza, Giovanna Guidoboni, Romain Hild, Daniele Prada, Christophe Prud’homme, R. Sacco, Lorenzo Sala, Marcela Szopos, In progess, 2020

[Sala2019] Lorenzo Sala. Mathematical modelling and simulation of ocular blood flows and their interactions.Numerical Analysis [math.NA]. Université de Strasbourg, 2019. English. NNT: 2019STRAD021 . tel02284233v2

[Škerget2014] Škerget, L. Tadeu, A., BEM numerical simulation of coupled heat, air and moisture flow through a porous solid, Engineering Analysis with Boundary Elements, 2014, 40, p154161