I solved the Poission equation which is given by \begin{equation} \Delta h = 0, \end{equation} where $h$, in my case, is the hydraulic pressure. Because I want to solve steady-state flow, the $\nabla h =[\frac{\partial h}{\partial x}, \frac{\partial h}{\partial y}]$ is the flux velocity.
I used the FEATOOL in MATLAB to solve the Poisson equation over a trapezpodal domain, with three-noded triangular elements. The left and right sides are Dirichlet boundary conditions ($h_L = 100$ and $h_R = 20$), and the top and bottom sides are Neumann boundary condition ($\nabla h \cdot \vec{n} = 0$, $\vec{n}$ is the normal vector of boundary). The solved pressure field is as follows: pressure field
Then I wanted to determine the velocity of each element. Since the pressure value of each node is obtained, I used linear shape function $\vec{N}$ to get velocity: \begin{equation} \frac{\partial h}{\partial x} = \frac{\partial \vec{N}^T}{\partial x}\vec{h_e} \text{ and } \frac{\partial h}{\partial y} = \frac{\partial \vec{N}^T}{\partial y}\vec{h_e}, \end{equation} where $\vec {h_e}$ is the three pressure values of an element, $T$ denotes transpose. The derivates of $\vec{N}$ are \begin{equation} \frac{\partial \vec{N}}{\partial x} = \left(b_1, b_2, b_3\right)^T, b_i = \frac{1}{2A}(y_j - y_k), \end{equation} where $A$ is the area of an element, and \begin{equation} \frac{\partial \vec{N}}{\partial y} = \left(c_1, c_2, c_3\right)^T, c_i = \frac{1}{2A}(x_k - x_j). \end{equation}
Using the above formulas, the velocity vectors are as follows: velocity vectors
Finally, I used the $\nabla h \cdot \vec{n} \cdot L_\Gamma$ ($L_\Gamma$ is the length of boundary edge of an boundary element) to calculate the inlet and outlet flux. I found the difference is around 10% (inlet = 10.1017, outlet = 8.9828, error = 0.1108). Why this happened? Is it due to the length difference between inlet and outlet (I found when they are equal, the fluxes are equal)?
I also used quadratic elements, but the same problem happened when the domain is irregular, particularly the lengths of inlet and outlet are different.
UPDATE
With RT0 elements, the inlet and outlet fluxes are the same now! If you are interested, download a matlab script from github (I just revised the original script a little and I still not very clear about mixed FEM).
This finite element basis (3-noded triangular element) satisfies no discrete conservation properties. In particular, the conservation properties present in the continuous formulation $\Delta h = 0$ do not carry over to the discrete formulation. Moreover, the boundary flux has a quite low convergence rate, if I recall correctly it's $O(h^{1/2})$ for linear elements.
There are other more exotic finite elements and variational formulations that may better satisfy discrete conservation properties. Search for information, e.g., about $H(div)$ conforming finite elements.