Solution (source code)

= Solution
{c}

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 have
$$
\begin{aligned}
\|e\|_1^2
&=a(e,u-v_h)+a(e,v_h-u_h)\\
&=a(e,u-v_h)\leq\|e\|_1\|u-v_h\|_1.
\end{aligned}
$$
If $e\ne0$, division proves $\|e\|_1\leq\|u-v_h\|_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_hu$. 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_hu\in V_h$. The allowed quadratic <finite element interpolation estimate> on a <shape-regular mesh> is
$$
\|u-I_hu\|_1\leq C_Ih^2\|u\|_{H^3(\Omega)},
$$
where $C_I$ depends on the fixed <triangle> and <shape regularity>, but not on $h$. Therefore
$$
\boxed{\|u-u_h\|_{H^1(\Omega)}\leq\|u-I_hu\|_{H^1(\Omega)}\leq C_1h^2\|u\|_{H^3(\Omega)},\qquad C_1=C_I.}
$$
This step uses the inclusion $I_hu\in V_h$ explicitly permitted by the question. Arbitrary smaller <vector subspaces> without the same approximation property would not give this rate.