I'm trying to understand how one approximates functions in $H_0^1(\Omega)$ by piecewise affine functions. The normal construction I have seen begins by breaking up your domain $\Omega$ into cubes $Q_i$ and on each $Q_i$ defining an *averaged gradient* $\frac{1}{|Q_i|} \int_{Q_i} Du =: \xi_i$. I can see why this approximates $Du$ well in the $L^2$ norm, but why can we then conclude that have in fact an $H_0^1(\Omega)$ function?

In other words, how can we "glue" the different affine pieces together? How do we know we can even do this? My first guess was to think in two dimensions and consider a cube and the four cubes adjacent to it. So let $\xi_1$ be the value of the gradient in the 'center cube'. Then if $\xi_2$ is the gradient on the piece on top, then I think we need to make the boundary normal to $\xi_1 - \xi_2$ so that we can have a continuous function across the boundary. I then suppose we need to do the same for the other pieces.

Is there a clearer way of seeing what is going on/how to do this?