The method of Lagrange multipliers gives the solution, but is indeed not the way that should be chosen for hybrid methods, because it leads to a very large linear system.
Instead one should use reduction, by inversion of the local operators.
In order to show how to implement an hybrid method with inversion of the local problems in FreeFem++, let me explain how to do on this example problem of the Laplace operator by the hybrid high order method of DiPietro and Ern.
1. Discontinuous spaces
Let me recall that Pkdc (with k an integer) is the space of functions that
are polynomial of degree at most k in each triangle, without any continuity requirement from one triangle to another.
Similarly, Pkedge is the space of functions that are defined only on edges of the mesh, and that are polynomial of degree at most k in each edge, without any continuity requirement from one edge to another.
Pkedgedc is similar to Pkedge except that we have two polynomial functions per edge, one for each side of the edge (or only 1 if the edge is a boundary edge), hence two different values when approaching the edge from one neighbouring triangle or the other.
2. The scheme
3. The reduced formulation
4. FreeFem++ code
HHO-reduced.edp (14.5 KB)