Solution (source code)

= Solution
{c}

Take an admissible variation $u+\epsilon\eta$ with $\eta=0$ on the boundary. Differentiating the functional at an extremum gives
$$
0=\int_V(F_u\eta+F_{u_x}\eta_x+F_{u_y}\eta_y+F_{u_z}\eta_z)dV.
$$
Integration by parts removes the derivatives from $\eta$; the boundary term vanishes. The <fundamental lemma of the calculus of variations> yields
$$
\boxed{\partial_xF_{u_x}+\partial_yF_{u_y}+\partial_zF_{u_z}-F_u=0.}
$$
For the specified squared-<gradient> energy this is $2\Delta u=0$, so $u$ is harmonic. The affine function $\ell=ax+by+cz$ is harmonic and satisfies the boundary data.

To verify a global minimum, write $u=\ell+v$ with zero boundary trace. Then
$$
I[u]=I[\ell]+2\int_V\nabla\ell\cdot\nabla v\,dV+\int_V|\nabla v|^2dV.
$$
The middle term vanishes by integration by parts, since $\Delta\ell=0$ and $v=0$ at the boundary. Hence $I[u]\geq I[\ell]$, and
$$
\boxed{\min I=(a^2+b^2+c^2)\operatorname{vol}V.}
$$
This proves minimality directly, rather than relying only on a necessary Euler-Lagrange equation.