Solution (source code)

= Solution

Integrate the normal component on each face directly, with outward normals. Each unlisted pair of coordinates ranges independently from zero to one:
$$
\begin{array}{c|c|c}
\text{face}&\mathbf v\cdot\mathbf n&\text{flux}\\\hline
x=0&0&0\\
x=1&2(y+z)&2\\
y=0&-x^2&-1/3\\
y=1&x^2-1&-2/3\\
z=0&-y^2&-1/3\\
z=1&y^2-1&-2/3
\end{array}
$$
For example, on $x=1$ the <integral> is $\int_0^1\int_0^1 2(y+z)\,dy\,dz=2$. Adding all six contributions gives
$$
\boxed{\int_{\partial V}\mathbf v\cdot d\mathbf S=0}.
$$
This calculation uses the six explicit <surface integrals>, independently of the <divergence theorem>.