An exact polynomial Poisson solution has a zero residual certificate

The factor $x(1-x)y(1-y)$ vanishes on all four edges. Differentiating gives $\mathrm{\Delta }u=-2x(1-x)-2y(1-y)$, so $f=2x(1-x)+2y(1-y)$ is polynomial and [F2] gives $u\in H_0^1$. For every $v\in C_c^{\mathrm{\infty }}(\mathrm{\Omega })$, integration by parts gives $(\mathrm{\nabla }u,\mathrm{\nabla }v)=(f,v)$; both sides are continuous in the H1 norm, so [F3] extends the identity to every $v\in H_0^1$. Thus $u$ is the weak Dirichlet solution. Also $\sigma =-\mathrm{\nabla }u$ and $\mathrm{div}\sigma =f$ are polynomials in $L^2$, so [F4] gives $\sigma \in H(\mathrm{div};\mathrm{\Omega })$.

By definition $\sigma +\mathrm{\nabla }u=0$ and $div\sigma =-\mathrm{\Delta }u=f$, so both components in [F1] vanish pointwise. The objective is nonnegative and therefore this exact zero residual is a global least-squares minimum.