Implementation of mixed hybrid finite element method

Question

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

1. 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?

2. what is the physical meaning of $\mu$ in the third equation (of the Eq. (\ref{MHFEM_weakform}))? basis function for $\mathbf{h_E}$?

3. 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.

4. how to deduce the exact form of $\mathbf{C}$?

UPDATE

I've achieved MHFEM. If you are interested, the matlab script is here.

Answer 1

In the hybrid method your basis for the flux is not $RT_0$ but it is the vector space of any piecewise polynomial, discontinuous functions. Then you introduce the third equation to get the normal continuity which you normally get for free when using $RT_0$. Physical interpretation is the flux, i.e. essentially the same as for $RT_0$ DOFs.

In the hybrid method you can give a specific value for the new flux variable on the boundary to get a Neumann condition, just like when using $RT_0$.

What do you mean by exact form of $C$? I think you have written the matrix $ij$ entry correctly in your question.