How to approximate numerically the gradient of the function on a triangular mesh
EDIT: We found a new in my opinion more elegant solution, see this answer.
Another way is first computing the gradient on each triangle: $\nabla f|_T = [a,b]$ where
$$\begin{bmatrix} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ x_3 & y_3 & 1 \end{bmatrix} \begin{bmatrix} a \\ b \\ c \end{bmatrix} = \begin{bmatrix} f(x_1,y_1) \\ f(x_2,y_2)\\ f(x_3,y_3) \end{bmatrix}$$
(we are fitting a plane $p(x,y) = ax+by+c$ on each triangle).
To compute the gradient at each node we just have to use a weighted average of all triangles that contain that node. We can do that by using the angle $\alpha_{n,T}$ of each triangle $T$ at our node $n$.
$$\nabla f_{n} = \frac{1}{2\pi} \sum_{T \in N} \alpha_{n,T}\nabla f|_T$$
Alternatively you can weigh each $\nabla f|_T$ by the area of the corresponding triangle (and normalize via the total area of the surrounding triangles). This has the advantage that we can also reuse the weights. Both methods do have the disadvantage that you also need to know all elements the given node is contained in, not only the edges.
Another method that just uses knowledge about the edges is fitting $p(x,y) = ax+by+c$ (via least squares) to all points that have edges to our node. This again is not very suitable for non-regular meshes.
The key thing you need to do is search "discrete differential geometry", which is essentially all about topics like this.
A quick-and-dirty answer that works pretty well for well-behaved meshes (i.e., ones where the valence of each vertex is pretty near 6, rather than being, say, 14...) and where triangles are more or less equilateral (no angles less than about 20 degrees, say), is this:
$$ \nabla f (v) = \sum_{u\in N(v)} \left( f(u) - f(v) \right), $$ where $N(v)$ is the set of vertices that are adjacent to $v$ (i.e. vertices $u$ with the property that $uv$ is an edge of the mesh).
This will be a somewhat "brute force" approach. The error in this approximation will scale as $\mathcal{O}(h)$ where $h$ is the mesh size.
Take a point $(x_0, y_0)$ with $i$ neighbors. We can compute derivatives in the direction of neighbors by:
$$\vec{f_i'} \approx \frac{f(x_0, y_0) - f(x_i, y_i)}{r_i} \begin{bmatrix} (x_0 - x_i)/r_i \\ (y_0 - y_i)/r_i \end{bmatrix}$$
where $r$ the distance is given by $\sqrt{(x_0 - x_i)^2 + (y_0 -y_i)^2}$.
From here: you have $i$ approximations for the gradient in the point. Take a weighted average of the the approximations (where weights are defined by the distance such that the closes points have the most influence)
$$\vec{f_{\text{final}}'} = \left(\sum_i \frac{1}{r_i}\right)^{-1} \left(\sum_i \frac{\vec{f_i'}}{r_i}\right)$$