Why the inlet and outlet fluxes are different when I solve the Poisson equation by FEM over a trapezoidal domain?
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.