For a continuous function $f:[a,b]\to R$ there is a natural and obvious procedure to approximate it with a sequence of continuous, piecewise linear functions: take $N$ equally spaced points in $[a,b]$ and join the corresponding points in the graph of $f$ with segments. If $f$ is $C^1$, this method has the additional benefit that it converges in the Lipschitz norm, without any tweaking.
Is there an equally simple procedure for functions of several variables? What is the most natural definition of a continuous, piecewise linear approximation for a continuous function $f:[a,b]^2\to R$, having the 'best possible' properties concerning $Df$ when $f$ is assumed to be $C^1$? The last part of the question is vague since I do not know exactly what to expect; ideally, I would like the approximating sequence to converge also in the Lipschitz norm. Note that piecewise linear functions, even in 1D, are not dense in Lip, but I want to approximate a $C^1$ function.
Background: I am writing some notes for an advanced calculus course, and a simple piecewise linear approximation procedure would greatly simplify some proofs. Since this must be a very familiar problem to people working in more applied fields, I thought it might be a good idea to ask them before rediscovering the wheel.