Suppose we have a piecewise linear real-valued function $g$ on a triangle mesh. $g$ is completely determined by its values $g_i$ at vertices. Although it is not really important we imagine that the mesh is embedded in $\mathbb{R}^3$. Let $\sigma$ be a triangle with vertices $i,j,k$ and $a,b,c$ the edge vectors of the embedding.

The restriction of $g$ to the triangle is a linear function. Its gradient is constant vector in the plane of the triangle. It is characterized by the fact that

\begin{align*}
\langle \textrm{grad}\,g,e_{jk} \rangle &= g_k-g_j \\
\langle \textrm{grad}\,g, e_{ki} \rangle &= g_i-g_k.
\end{align*}

Because of $c=-a-b$, the third identity

$\langle \textrm{grad}\,g, e_{ij} \rangle = g_j-g_i$

is a consequence of the two above.

The normal vector $N$ to the triangle is given by

$\displaystyle N=\frac{1}{2A}\,e_{jk}\times e_{ki}=\frac{1}{2A}\,e_{ki}\times e_{ij}=\frac{1}{2A}\,e_{ij}\times e_{jk}$

where $A = \frac{1}{2} |a\times b|$ denotes the area of $\sigma$. Then one can verify that the vector

$\displaystyle v:= – \frac{1}{2A} N\times (g_i e_{jk}+g_j e_{ki}+g_k e_{ij})$

solves the above equations. For example,

\begin{eqnarray*}
\langle v , e_{jk} \rangle &=& -\frac{1}{2A} \langle N \times (g_i e_{jk}+g_j  e_{ki}+g_k e_{ij}),e_{jk} \rangle\\
&=& -\frac{1}{2A} (g_j \langle e_{ij} \times e_{jk},N\rangle +g_k \langle e_{ij} \times e_{jk},N \rangle)\\
&=& g_k-g_j.
\end{eqnarray*}

Therefore $\textrm{grad}\,g=v$. The gradient of $g$ rotated by $90^\circ$ in the plane of the triangle is obtained by an even nicer formula than $\textrm{grad}\,g$ itself:

$\displaystyle N \times \textrm{grad}\,g =\frac{1}{2A}( g_i e_{jk}+g_j e_{ki}+g_k e_{ij}).$

Here we just take a linear combination of the edge vectors with the value of $g$ at the opposing vertex as the respective coefficient.