Certainly, a piecewise affine function $f$ (meaning, a function which is affine on each open simplex of some triangulation of the domain) is in $W^{1,\\ p}_{loc}$ for some $1\leq p\leq \infty$ if and only if it is continuous (the discontinuity set would be otherwise too large, a hypersurface). In particular, you really need your approximating function to be continuous, in order to be of class $H^1$. On the other hand, continuous functions that are affine on each cube of a cubic lattice subdivision, are too a rigid class: they are of the form $f(x_1,\dots,x_n)=\phi_1(x_1)+\dots+\phi_n(x_n)$, with $\phi_i$ some continuous piecewise affine functions of one variable. If you instead consider general continuous piecewise affine functions (as defined above), you get a class with more satisfactory density properties. Since $C^1_0(\Omega)$ functions are $H^1$-dense in $H^1_0(\Omega)$, it is sufficient to approximate a function in $C^1_0(\Omega)$ by $C^0$ piecewise affine functions in the $C^1$ sense; this is classically done by affine interpolation on the points of the 0-skeleton of the triangulation (there is a unique such interpolation: that's the advantage of triangulations compared to cubic subdivisions). So in your case you should further subdivide each cube in $n!$ simplices, and take the approximation affine on each of them. It should not be difficult to bound the $H^1$ distance of the approximation; I guess that books of numerical analysis cover completely this point.