The mixed hybrid finite element method (MHFEM) is based on the mixed finite element method (MFEM). So, I'd recall the implementation of MFEM.
The mixed formulation of Poisson equation reads $$\begin{equation}\begin{split}\mathbf{u}=\nabla h \text{ in } \Omega, \\ \nabla\cdot\mathbf{u}=0 \text{ in } \Omega, \\ \mathbf{u \cdot n} = 0 \text{ on } \partial\Omega_N, \\ h = h_D \text{ on } \partial\Omega_D \end{split}\tag{1} \label{darcy}\end{equation}$$ where $\mathbf{u}$ is velocity, $h$ is hydraulic head, and $\mathbf{n}$ is outer unit normal.
By introducing the $RT0$ basis function, the discrete form of the Eq. (\ref{darcy}) can be written as $$\begin{equation} \begin{split}% \mathbf{u_E}&\int_\Omega \phi_i \cdot \phi_j + h_T\int_\Omega \nabla \cdot \phi_i &= h_D\int_{\partial \Omega_D} \phi_i \cdot \mathbf{n},\\% \mathbf{u_E}&\int_{\Omega} \nabla \cdot\phi_i &= 0, % \end{split} \tag{2} \label{weakform}\end{equation}$$ where $\phi_i$ is a vector (basis function for interpolating normal velocity of each edge of the triangular element), $\mathbf{u_E}$ is the normal components of velocity on each edge, and $h_T$ is the mean head of an element. Therefore, the global matrix can be written as $$\begin{equation}% \left[\begin{matrix}% \mathbf{A} &\mathbf{B^T}\\ % \mathbf{B} \end{matrix}\right]% \left[\begin{matrix}% \mathbf{u_E}\\ % \mathbf{h_T} \end{matrix}\right]=% \left[\begin{matrix}% \mathbf{b_D}\\ % 0 \end{matrix}\right].% \end{equation} \tag{3} \label{matrix_MFEM}$$
Now, let's switch to MHFEM. I have read some papers on that. If I understood correctly, the MHFEM is a kind of method that introduces another variable, i.e. pressure on edge, $\mathbf{h_E}$, and then eliminates $\mathbf{u_E}$ and $\mathbf{h_T}$. The $\mathbf{h_E}$ is also called as the Lagrange multiplier.
Hence, the Eq. (\ref{weakform}) can be rewritten as $$\begin{equation} \begin{split}% \mathbf{u_E}&\int_\Omega \phi_i \cdot \phi_j + h_T\int_\Omega \nabla \cdot \phi_i - \mathbf{h_E}\int_{\partial \Omega} (\phi_i \cdot \mathbf{n})\mu_j&= 0, % \left(h_E = h_D \text{ if $\partial\Omega = \partial\Omega_D$ } \right)\\% \mathbf{u_E}&\int_{\Omega} \nabla \cdot\phi_i &= 0, \\% -\mathbf{u_E} &\int_{\partial \Omega}(\phi_i\cdot \mathbf{n})\mu_j &=0.% \end{split} \tag{4} \label{MHFEM_weakform}\end{equation}$$ Therefore, the matrix now is $$\begin{equation}% \left[\begin{matrix}% \mathbf{A} &\mathbf{B^T} &\mathbf{C^T}\\ % \mathbf{B} \\ % \mathbf{C} % \end{matrix}\right]% \left[\begin{matrix}% \mathbf{u_E}\\ % \mathbf{h_T} \\ % \mathbf{h_E} % \end{matrix}\right]=% \left[\begin{matrix}% \mathbf{b_D}\\ % 0 \\ % 0 % \end{matrix}\right],% \end{equation} \tag{5} \label{matrix_MHFEM}$$ and the $\mathbf{u_E}$ and $\mathbf{h_T}$ can be vanished by some implements.
My questions are
Why the third equation (of the Eq. (\ref{MHFEM_weakform})) is introduced. A lot of papers said, the reason is that the continuity condition on the normal components of the $\mathbf{u}$ across the finite element edges or faces is relaxed, i.e. $\mathbf{u}$ is allowed to be discontinuous across element interfaces; and so, the continuity of the normal flux across the interelement boundaries is enforced by Lagrange multipliers ($\mathbf{h_E}$). But I see the continuity equation (i.e. the 2nd equation in Eq. (\ref{MHFEM_weakform})) still exists in the linear system, why the velocity is discontinuous at element faces?
what is the physical meaning of $\mu$ in the third equation (of the Eq. (\ref{MHFEM_weakform}))? basis function for $\mathbf{h_E}$?
how to impose the Neumann condition in the linear system (MHFEM). In MFEM method, this can be considered in the matrix by, let's say, if normal velocity on edge $i$ is known as zero, then just set elements in the $(:, i)$ and $(i, :)$ of stiffness matrix to be zero, and $(i, i) = 1$, and the $(i, 1)$ of right hand side is set to zero.
how to deduce the exact form of $\mathbf{C}$?
UPDATE
I've achieved MHFEM. If you are interested, the matlab script is here.