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

Question

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.

Answer 1

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.

Answer 2

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

Answer 3

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

Answer 4

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.

Answer 5

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.