Calculating the gradient of a triangular mesh

Question

This question is regarding a homework assignment (want to be upfront because many people on stack exchange don't like helping with homework questions).

I have a triangular mesh, which is defined by a set of vertices and a set of faces. I want to calculate the discrete gradient of the mesh, and I even got instructions on how to implement it. The problem is that even after reading them multiple times, I don't understand them...

The implementation instructions:

My Problem:

I don't understand how to implement the marked sum, which in turn prevents me from constructing the following E matrix.

I would appreciate any help. I'm not looking for specific code examples, rather an explanation/example of the calculation for a single triangle.

Answer 1

A function $f$ on your triangulated surface is an assignment of real numbers on the vertices of the triangulation and therefore $f \in \mathbb{R}^{|\mathcal{V}|}$. According to the continuous theory, the gradient should be a linear (differential) operator that takes a function and spits out a vector field. Hence, in the discrete theory, the gradient should be a linear operator that takes a function $f \in \mathbb{R}^{|\mathcal{V}|}$ and spit out some kind of a discrete version of a vector field. In your case, the vector field is given as an assignment of a vector to each triangle of the triangulation. That is why the discrete gradient operator should be a linear operator, i.e. a linear transformation $$\text{grad} : \mathbb{R}^{\mathcal{|V|}} \to \mathbb{R}^{3|\mathcal{F}|}$$
which is represented by a matrix $\text{grad} \in \mathbb{R}^{3|\mathcal{F}| \times |\mathcal{V}|}$. In the continuous theory, the gradient is determined only by local data, i.e. the behaviou of the function near the point at which the gradient is being calculated and by the geometry near that point. Therefore, in the discrete theory, the gradient is determined only by the geometry of one triangle and the values of the function at the triangle's vertices. Let the triangle $j \in \mathcal{F}$ have vertices denoted (according to the instructions' convention) by $1, 2, 3$. The labels are set up so that the vertices $1,2,3$ define the labels of the opposing edges $\vec{e}_{j_1}, \vec{e}_{j_2}, \vec{e}_{j_3}$, i.e. for the vertex $i$ the edge opposite to the vertex $i$ is denoted by $\vec{e}_{j_i}$, where the latter vector is oriented so that the vectors form a closed chain of vectors summing to zero in a counterclockwise direction.
Hence, the gradient can be defined as $$(\text{grad}\,f)(j) = \frac{1}{2\text{Area}(j)}\, \sum_{i=1}^{3} f_i\, J\vec{e}_{j_i}$$ Or to write it in coordinates, $$\begin{bmatrix} (\text{grad}\,f)(j)_1\\ (\text{grad}\,f)(j)_2\\(\text{grad}\,f)(j)_3\end{bmatrix} = \frac{1}{2\text{Area}(j)}\left( f_1\, \begin{bmatrix} (Je_{j_1})_1\\ (Je_{j_1})_2\\(Je_{j_1})_3\end{bmatrix} + f_2\, \begin{bmatrix} (Je_{j_2})_1\\ (Je_{j_2})_2\\(Je_{j_2})_3\end{bmatrix} + f_3\, \begin{bmatrix} (Je_{j_3})_1\\ (Je_{j_3})_2\\(Je_{j_3})_3\end{bmatrix}\right)$$ and in matrix form $$\begin{bmatrix} (\text{grad}\,f)(j)_1\\ (\text{grad}\,f)(j)_2\\(\text{grad}\,f)(j)_3\end{bmatrix} = \frac{1}{2\text{Area}(j)}\, \begin{bmatrix} (Je_{j_1})_1 & (Je_{j_2})_1 & (Je_{j_3})_1 \\ (Je_{j_1})_2 & (Je_{j_2})_2 & (Je_{j_3})_2\\(Je_{j_1})_3 & (Je_{j_2})_3 & (Je_{j_3})_3 \end{bmatrix} \begin{bmatrix} f_1\\ f_2\\f_3\end{bmatrix} $$ So, for the matrix $E \in \mathbb{R}^{3|\mathcal{F}|\times|\mathcal{V}|}$ is defined as follows. Let $j$ be a triangle and $i$ one of its vertices. Take the edge $e_{j_i}$ of triangle $j$ opposite to the vertex $i$, and consider the vector $\vec{e}_{j_i}$ along the edge, starting from one vertex, different from $i$, to the other vertex, different from $i$, oriented so that the the vector follows counterclockwise orientation of the boundary of the triangle $j$. Then rotate the vector $\vec{e}_{j_i}$ by $90^{\circ}$ so that it points inside the triangle $j$. Denote this rotated vector by $$J\vec{e}_{j_i} = \begin{bmatrix} (Je_{j_i})_1\\ (Je_{j_i})_2\\(Je_{j_i})_3\end{bmatrix}$$ Then fill up the following entries of the matrix $E \in \mathbb{R}^{3|\mathcal{F}|\times|\mathcal{V}|}$:
\begin{align} &E(j,i) = (Je_{j_i})_1\\ &E(j + |\mathcal{F}|, i) = (Je_{j_i})_2\\ &E(j + 2|\mathcal{F}|, i) = (Je_{j_i})_3 \end{align}
You can rotate the vector $\vec{e}_{j_i}$ by using the cross-product operation twice. Take two of the edge-vectors of the triangle $j$, say $\vec{e}_{j_1}$ and $\vec{e}_{j_2}$, and calculate the unit normal vector to the triangle $j$: $$\vec{n}_j = \frac{1}{\|\vec{e}_{j_1} \times \vec{e}_{j_2}\|}\,\Big(\vec{e}_{j_1} \times \vec{e}_{j_2}\Big)$$ It doesn't matter whcih two out of the three edge-vector you would choose, the result is always the same unit normal $\vec{n}_j$. Then $$J\vec{e}_{j_i} = \vec{n}_j \times \vec{e}_{j_i}$$

Answer 2

In a homework we were tasked to calculate the gradient of a triangle respective to a corner vertex. The instructions didn't say what to calculate and mentioned, one could derive it by a derivative of the triangle area using the corner vertex coordinates. I exactly came to the same formula explained by Futurologist.

In case someone wants to know, this is the derivation:

Let $\vec{x},\vec{y}, \vec{z} \in \mathcal{V}$ be the vertex corners of the triangle.

First, the area of a triangle is the length of the cross product between two sides of it divided by 2 (because the determinant calculates the area of a parallelogram which is twice the area of the triangle): $A(\vec{x}) = \left|(\vec{y}-\vec{x})\times(\vec{z}-\vec{x})\right|/2 = \sqrt{\left((\vec{y}-\vec{x})\times(\vec{z}-\vec{x})\right)^2}/2$

The multiplication of vectors should be the dot product whose meaning depends on the operands left and right to it. Using differential calculus rules for vectors, this is what I get:

$\frac{\partial A}{\partial \vec{x}} = \frac{\frac{\partial \left((\vec{y}-\vec{x})\times(\vec{z}-\vec{x})\right)^2}{\partial \vec{x}}}{4\sqrt{\left((\vec{y}-\vec{x})\times(\vec{z}-\vec{x})\right)^2}}$ $= \frac{\frac{\partial \left((\vec{y}-\vec{x})\times(\vec{z}-\vec{x})\right)}{\partial \vec{x}} \cdot 2\left((\vec{y}-\vec{x})\times(\vec{z}-\vec{x})\right) }{4A}$

How to calculate the gradient of a cross product? If we would look for the divergence, we would use the determinant $\text{det}(\nabla_\vec{x},\vec{y}-\vec{x},\vec{z}-\vec{x})$ but that's not the gradient. The gradient $\partial A/\partial \vec{c}$ must be a vector itself so that $\partial(\vec{a}\times\vec{b})/\partial \vec{b}$ must be a matrix! Now, just think of $M \times \vec{b}$ like a matrix product but using $M_i \times \vec{b}$ to replace the dot product $M_i \cdot \vec{b}$ of a normal matrix multiplication (each result component will now be a row vector, not a scalar).

Next step is to split those cross products up into a sum:

$\frac{\partial A}{\partial \vec{x}}= \frac{\frac{\partial \left(\vec{y}\times\vec{z}-\vec{x}\times\vec{z} - \vec{y}\times\vec{x} + \vec{x}\times\vec{x}\right)}{\partial \vec{x}} \cdot \left((\vec{y}-\vec{x})\times(\vec{z}-\vec{x})\right) }{2A}$

$\vec{y}\times\vec{z}$ does not depend on $\vec{x}$, $\vec{x}\times\vec{x}$ is $\vec{0}$. It now is used that $\partial(\vec{a}\times\vec{b})/\partial \vec{a} = I\times\vec{b}$, where $I$ is the eye matrix. It is a fact that $Y_\times = I\times\vec{y}$ is the (left) cross product matrix from $\vec{y}$ so that $\left(I \times \vec{z}\right)^T = \left(I\times(-\vec{z})\right)$ is the right cross product matrix from $\vec{z}$. Hence

$\frac{\partial A}{\partial \vec{x}}= \frac{\frac{\partial \left(\vec{x}\times\vec{y}-\vec{x}\times\vec{z}\right)}{\partial \vec{x}} \cdot \left((\vec{y}-\vec{x})\times(\vec{z}-\vec{x})\right) }{2A}$ $= \frac{\left(I\times\vec{y}-I\times\vec{z}\right) \cdot \left((\vec{y}-\vec{x})\times(\vec{z}-\vec{x})\right) }{2A}$ $= \frac{\left(I\times(\vec{y}-\vec{z})\right) \cdot \left((\vec{y}-\vec{x})\times(\vec{z}-\vec{x})\right) }{2A}$ $= \frac{\left((\vec{y}-\vec{x})\times(\vec{z}-\vec{x})\right)\times(\vec{z}-\vec{y}) }{2A}$

So this result is exactly the same as $-\frac{\mathcal{J}\vec{e_{j_i}}}{2} = \frac{(-\vec{n})\times\vec{e_{j_i}}}{2}$ where $\vec{e_{j_i}} = \vec{z}-\vec{y}$, $i = \vec{x}$ and the direction will always show towards the opposite edge of $\vec{x}$ indepent of which vertex is chosen to be $\vec{y}$ and which is $\vec{z}$. That's the reason, why the $2$ in the original formula for $(\text{grad} f)(j)$ is needed, to divide away the factor of $2$ within the $\vec{e_{j_i}}$.

If you work with parallelograms instead of triangles you can remove the $2$ entirely.

What does the gradient of a triangle mean after all? The triangle gradient is the direction, where -- if you move $\vec{x}$ in that direction -- the triangle experiences the biggest change in area. You can add together the gradient of multiple triangles from $\vec{x}$ to use it for minimizing surface area of a discrete 3D surface.