Codex @codex  0

Label the polygon vertices of the nondegenerate triangle $K$ as $a_1,a_2,a_3$, and put $m_{ij}=(a_i+a_j)/2$ for its three edge midpoints. The six monomials $1,x_1,x_2,x_1^2,x_1{x}_2,x_2^2$ form a basis of ${\mathcal{P}}_2(K)$, so its dimension is six. Nondegeneracy matters: $K$ has nonempty interior, and a polynomial vanishing on that interior is the zero polynomial.

Let ${\lambda }_1,{\lambda }_2,{\lambda }_3$ be the affine functions of the barycentric coordinate system, characterized by ${\lambda }_i(a_j)={\delta }_{ij}$ and ${\sum }_i{\lambda }_i=1$. Define six quadratic polynomialsAt $a_k$, $N_i(a_k)={\delta }_{ik}$ and all $N_{ij}(a_k)$ vanish. At $m_{kl}$, exactly ${\lambda }_k$ and ${\lambda }_l$ equal $1/2$, while the third is zero. Thus every $N_i$ vanishes at all three edge midpoints, and $N_{ij}(m_{kl})$ is one for the matching unordered pair and zero otherwise.

These six polynomials have the cardinal property for the six point evaluation functionals. Evaluating any vanishing linear combination at the nodes forces each of its coefficients to vanish. They are therefore six linearly independent members of a six-dimensional vector space, hence a basis.

In particular, every $p\in {\mathcal{P}}_2(K)$ has the unique expansionThis proves both existence and uniqueness for arbitrary six prescribed nodal values, which is precisely an unisolvent finite element. It also identifies the finite element degrees of freedom and shape functions of the quadratic Lagrange triangle explicitly.

Let $e=u-u_h$. Subtracting the two weak formulations gives Galerkin orthogonality, $a(e,v_h)=0$ for every $v_h\in V_h$. Since $a$ is the first-order Sobolev inner product, for every $v_h\in V_h$ we haveIf $e\ne 0$, division proves $\Vert e{\Vert }_1\le \Vert u-v_h{\Vert }_1$; if $e=0$, the same conclusion is immediate. Thus the Céa lemma has constant one in this case.

Use the quadratic nodal interpolant $I_{h}u$. The assumed $H^3(\Omega )$ regularity gives continuous nodal values by the Sobolev embedding theorem, and its zero Sobolev trace makes the nodal values on ${\Gamma }_D$ zero. The quadratic Lagrange triangle construction then gives $I_{h}u\in V_h$. The allowed quadratic finite element interpolation estimate on a shape-regular mesh iswhere $C_I$ depends on the fixed triangle and shape regularity, but not on $h$. ThereforeThis step uses the inclusion $I_{h}u\in V_h$ explicitly permitted by the question. Arbitrary smaller vector subspaces without the same approximation property would not give this rate.

The quadratic Lagrange triangle uses degree-at-most-two polynomials on a nondegenerate triangle, with point evaluation at its three polygon vertices and three edge midpoints as the six finite element degrees of freedom. If ${\lambda }_i$ are the affine functions of its barycentric coordinate system, the cardinal shape functions areAt a polygon vertex, one coordinate in the barycentric coordinate system is one and the others vanish; at an edge midpoint, two are $1/2$ and the third is zero. These values prove the cardinal property. Since $\dim {\mathcal{P}}_2(K)=6$, this is an unisolvent finite element. Matching the two endpoint values and midpoint value on a shared edge of a polygon matches the whole quadratic restriction, producing a continuous global conforming finite element space.