# How to approximate numerically the gradient of the function on a triangular mesh

Given an arbitrary (lets say 2D) triangular mesh, with known $(x_i,y_i)$ locations of points, and numerical values of a function $f$ on them (either in the nodes, or in the centroids of the triangles, doesn't matter) like this random example, how can I obtain a numerical approximation of the direction and modulus of the gradient in the nodes?

I can get directional gradients to adjacent nodes or elements, but I don't know how to join them to create a general gradient. I have the feeling that this will involve some sort of least squares solution.

PD: if its relevant, we can assume linear functions on the triangles.

EDIT: I am unaware of the terminology, but as a initial test I can check if the gradient approximates linearly and quadratically the real gradient (by testing it against a known plane and parabola). Currently reading more on this.

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.

• Weighting by the inverse of area would be better than weighting by area. The larger the area, the farther away the other corners of the triangle are, which should weaken their influence, not increase it. Jan 30, 2018 at 17:37
• @PaulSinclair Can you think of a concrete example where this produces a better result? I think when we have a large triangle that e.g. has gradient zero, then it would make sense for the approximations of the gradient and the corresponding nodes to be close to zero too. Jan 30, 2018 at 17:51
• The only concrete examples I can think of where it wouldn't produce a better result are those where the triangles are all roughly the same size or the function is close to constant. Try $f(x,y) = \frac 1{xy}$, with a fine-grained mesh for $x, y < 1$, but very course-grained for $x > 1, y > 1$. Now $$\nabla f(1,1) = (1,1)$$. For triangles whose other vertices are $< 1$, the small triangles will be very close to $(1,1)$ in gradient, for the large triangles on the other side, their gradients will be much closer to $(0,0)$. But they will have the large areas, not the small triangles. Jan 30, 2018 at 19:52
• But a better approach would not be to weight by area at all. Gradient is a local property. The farther the other points are from the point in question, the less reliable the estimate of gradient you will get from them will be. But area - even inverse area - doesn't correspond very well with distance. Weighting by the inverse of the max length of the two sides meeting at your target vertex would be better. Jan 30, 2018 at 20:06
• @PaulSinclair I am test the different approaches proposed here and testing them with an arbitrary mesh with different amount of details, by applying a known plane/parabola as $f$. I suggest you add your idea as an answer Jan 31, 2018 at 10:54

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).

• Ah, yes, however I can certainly not assume this. As I am not really doing FEM, just using FEM-type meshes, I am good (and actually in some cases I look for) with triangles with ugly aspect ratios, or nodes with high valences. Thanks for they key words though! Jan 30, 2018 at 13:36

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)$$

• Yes, I'll correct it now. We want the direction vector to have length one. Jan 30, 2018 at 16:31
• This method unfortunatley does not even converge to the actual result if we consider a linear function $f$. To understand why it is probably easiest to consider the function $f(x,y) = x$ and e.g. $(x_0,y_0) = (0,0), (x_1,y_1) = (-1,0), (x_2,y_2) = (0,-1), (x_3,y_3) = (1,1)$. Feb 2, 2018 at 10:50
• Based on your idea we came up with a similar approach that solves this problem, see this answer. Feb 2, 2018 at 11:12

This idea is based on @PiotrHughes' answer, which unfortunately does not provide good results even for the simple case of a linear function.

Idea: The idea is following: We want to compute the gradient at $P_0 = (x_0,y_0)$ which has edges to $P_i = (x_i,y_i)$ for $i=1,2,3,...n$. Let us consider the unit vectors $d_i$ pointing from $P_0$ to $P_i$, so let $r_i = P_i - P_0$ and $d_i = \frac{1}{||r_i||}(r_i)$.

Then the scalar product $\nabla f(x_0,y_0) \cdot d_i$ is the directional derivative in the direction $d_i$. These directional derivatives can be approximated by $$\nabla f(x_0,y_0) \cdot d_i \approx \frac{f(x_i,y_i)-f(x_0,y_0)}{||r_i||}$$

This gives us an equation for each edge, and if we have at least two non-colinear edges we can solve this system for $\nabla f(x_0,y_0)$ using least squares.

Implementation: We can simplify this system to

$$\begin{bmatrix} x_1-x_0 & y_1-y_0 \\ x_2-x_0 & y_2-y_0 \\ \vdots & \vdots \\ x_n - x_0 & y_n - y_0 \end{bmatrix} \nabla f(x_0,y_0) = \begin{bmatrix} f(x_1,y_1) - f(x_0,y_0) \\ f(x_2,y_2) - f(x_0,y_0) \\ \vdots \\ f(x_n,y_n) - f(x_0,y_0) \end{bmatrix}$$

Convergence: The only approximation we made is the use of the finite differences for the approximation of the directional derivatives. Assuming we have a regular and quasi uniform family of meshes these should converge with the known order. And it certainly provides exact results if $f$ is linear.

EDIT: Equivalently we can rewrite the equation from above as

$$\begin{bmatrix} x_1-x_0 & y_1-y_0 & 1\\ x_2-x_0 & y_2-y_0 & 1\\ \vdots & \vdots &\vdots \\ x_n - x_0 & y_n - y_0 & 1\end{bmatrix} \begin{bmatrix} \nabla f(x_0,y_0) \\ f(x_0,y_0)\end{bmatrix}= \begin{bmatrix} f(x_1,y_1) \\ f(x_2,y_2) \\ \vdots \\ f(x_n,y_n) \end{bmatrix}$$

Now it becomes very apparent that we're just fitting a plane through the surrounding points. Which was what I suggested as a comment here.

• This answer is, in a plane, equally as good as your other answer, but this one also works on boundary nodes. Feb 2, 2018 at 11:28
• To add as extra info: The bottom answer in here is equivalent to the Prewitt operator in a regular mesh. If a weigthed least squares is performed, with weigths $1/r^2$, then its equivalent to a Sobel operator. Feb 12, 2018 at 15:01

First, if you extend $f$ linearly across each triangle, then the resulting function cannot be extended to a smooth function on 3-space unless the triangle all lie in a plane. To put it more bluntly: your problem is ill-posed: you want the gradient of a non-differentiable function.

Perhaps you'll say "Well, just give me the gradient of any differentiably function that has those values at those vertices", and again it's ill-posed, for there and infinitely many such functions.

What you really need to do is specify the problem more precisely (which is tough to do). You might, for instance, say "Among all diff'l functions with those values, find the one whose squared gradient magnitude has the smallest integral over the surface". Without some such criterion, there's no reasonable answer to your question.

In particular, there's a smooth function $f$, meeting your vertex values, with the property that its gradient at every vertex is zero.

• Pardon my illiteracy in math, but would't the same apply for example to a mesh that is composed of squares, i.e. an image? We can however approximate the gradient in an image very easily. I don't want to find the gradient of the fucntion in the whole domain, just know its value in the nodes of the mesh Jan 30, 2018 at 19:19
• Yes, you are technically correct, but I assumed it was clear enough when the question said that I needed a "numerical approximate of the [..] gradient in the nodes". You think I should edit the question to add the constrain? To be honest, I'd be even happy to explore other cases Jan 31, 2018 at 9:19
• Most math problems have an answer, and we can check that the answer is right. I'm doing my best to turn yours into one of those. That's actually a valuable skill for those who work with mathematics. I'm sorry I can't just "give you the answer", but... that's life. I mean, I did give you an answer, but I had the decency to point out that it had limitations, and you didn't like those. If you prefer other answers, ones that don't point out their limitations, so be it. Jan 31, 2018 at 12:42
• I liked your answer as much as any of the others, in fact I am testing all of them. I am genuinely asking what limitations I should add to the question/ how to understand the limitation at hand so I can improve the problem definition Jan 31, 2018 at 12:44
• You ask Wouldn't "How do I formulate a question about computing gradients of real-valued functions interpolated from values at mesh vertices?" be too broad for math.se? It might be, but with an explanation of why just asking for a formula doesn't work, it might be a good and interesting question. Jan 31, 2018 at 18:39