Solution (source code)

= Solution

Take a smooth <variation> $u+\varepsilon\eta$ with $\eta$ vanishing on the boundary; compactly supported test functions suffice for the interior equation. At an extremum, the <first variation> is zero:
$$
0=\left.\frac{d}{d\varepsilon}I[u+\varepsilon\eta]\right|_{\varepsilon=0}
=\int_V\left(F_u\eta+F_{u_x}\eta_x+F_{u_y}\eta_y+F_{u_z}\eta_z\right)\,dV.
$$
<Integration by parts> removes the derivatives from $\eta$. The boundary terms vanish, giving
$$
\int_V\left(F_u-\partial_xF_{u_x}-\partial_yF_{u_y}-\partial_zF_{u_z}\right)\eta\,dV=0.
$$
The <fundamental lemma of the calculus of variations> gives the <Euler-Lagrange equation>
$$
\boxed{\partial_xF_{u_x}+\partial_yF_{u_y}+\partial_zF_{u_z}-F_u=0}.
$$
This is the classical equation for sufficiently smooth admissible extrema; for weaker regularity the first-<variation> identity supplies its weak form.

For the squared-<gradient> integrand, the equation is $2(u_{xx}+u_{yy}+u_{zz})=0$, namely the <Laplace equation>. Put $u_0=ax+by+cz$ and write any admissible <function> as $u=u_0+v$, with $v=0$ on the boundary. The <energy> identity is
$$
I[u]=(a^2+b^2+c^2)\operatorname{vol}V
+2\int_V(a,b,c)\cdot\nabla v\,dV+\int_V|\nabla v|^2\,dV.
$$
By the <divergence theorem>, the middle integral is $\int_{\partial V}v(a,b,c)\cdot n\,dS=0$; equivalently it vanishes by weak <integration by parts> for zero-trace functions. The final term is nonnegative and vanishes at $v=0$. Thus this is a global minimum, not merely a stationary value:
$$
\boxed{\min I=(a^2+b^2+c^2)\operatorname{vol}V}.
$$