Partial derivatives for triangular meshes (in 3D)

Question

A grid offers an obvious definition for the partial derivatives at a grid point, given

$x$ the value of a point $p$ in an $n$ dimensional grid, the forward partial derivative that point for coordinate $i$ would be:

$$ \frac{x - x_{p_{i-1}}}{h}$$

i.e. take the value at $p$ the value at the point "behind" $p$ on the chosen dimension, take the difference, divide by the grid spacing.

This works with grids because there is an intrinsic relationship between a grid and the axial directions.

But let's say I have a triangular mesh embedded in 3D. How do I define differential operators? Like if I have values at the nodes of the mesh, how do I compute a gradient, or a hessian, or any other differential quantity?

Answer 1

I will give an answer for the two dimensional case. The extension to three dimensions is straight forward. Consider a scalar field $f$ as a function of $f(x(\xi,\eta),y(\xi,\eta))$, where each element has its own reference space.

You can start from the physical side

$$ \begin{align} \frac{ \partial f}{ \partial x} = \frac{ \partial f}{ \partial \xi} \frac{ \partial \xi}{ \partial x} + \frac{ \partial f}{ \partial \eta} \frac{ \partial \eta}{ \partial x}, \\ \frac{ \partial f}{ \partial y} = \frac{ \partial f}{ \partial \xi} \frac{ \partial \xi}{ \partial y} + \frac{ \partial f}{ \partial \eta} \frac{ \partial \eta}{ \partial y}, \end{align} \qquad \text{short} \qquad \begin{pmatrix} f_x \\ f_y \end{pmatrix}= \begin{pmatrix} \xi_x & \eta_x \\ \xi_y & \eta_y \end{pmatrix} \begin{pmatrix} f_{\xi} \\ f_{\eta} \end{pmatrix}. $$

Or you start from the computational side with

$$ \begin{align} \frac{ \partial f}{ \partial \xi} = \frac{ \partial f}{ \partial x} \frac{ \partial x}{ \partial \xi} + \frac{ \partial f}{ \partial y} \frac{ \partial y}{ \partial \xi}, \\ \frac{ \partial f}{ \partial \eta} = \frac{ \partial f}{ \partial x} \frac{ \partial x}{ \partial \eta} + \frac{ \partial f}{ \partial y} \frac{ \partial y}{ \partial \eta}, \end{align} \qquad \text{short} \qquad \begin{pmatrix} f_{\xi} \\ f_{\eta} \end{pmatrix}= \begin{pmatrix} x_{\xi} & y_{\xi} \\ x_{\eta} & y_{\eta} \end{pmatrix} \begin{pmatrix} f_{x} \\ f_{y} \end{pmatrix}. $$

Now since both sides are identical it holds

$$ \begin{pmatrix} \xi_x & \eta_x \\ \xi_y & \eta_y \end{pmatrix} = \begin{pmatrix} x_{\xi} & y_{\xi} \\ x_{\eta} & y_{\eta} \end{pmatrix}^{-1} = \frac{1}{J}\begin{pmatrix} y_{\eta} & -y_{\xi} \\ -x_{\eta} & x_{\xi} \end{pmatrix} .$$

Simply replace your derivatives with the relations above, before discretization. On arbitrary elements you also have to discretize the covariant metric terms $ x_{\xi}, y_{\xi}, x_{\eta}, y_{\eta} $ and the Jacobian $J$. Note that, the derivatives are continuous for each element, only.

If you want to calculate derivatives across element boundaries you may rather use a least squares method.

Answer 2

This may or may not help you, but with practice it is often easier to work discretely with integral rather than differential calculus. There is then no need to invoke arguments in which you approach the limit and then back off from it. For instance if you are finding derivatives in order to compute a divergence, note instead that $div(u,v)\simeq(1/S)\oint(vdx-udv)$ for a small control volume.

Answer 3

Usually, a mesh with "values" at the vertices is the result of some Galerkin scheme or interpolation.

I say "values", because those are really degrees of freedom in a basis of functions. It just so happens that a very commonly used basis of functions is the Lagrange basis with the property that the coefficients in the basis are evaluations of the function. This is not true of all bases. Take the Bernstein (Bézier) basis, for instance. In this case, you'd have to compute linear combinations of the DoFs to obtain values at the vertices. Another case is with discontinuous representations, where solutions are described as polynomial over each element but only $L^2$ over the entire domain (e.g. Discontinuous Galerkin schemes); this means values at the vertices are multiple (one per adjacent element). Yet another example is Hermite FE where some DoFs relate to gradients directly.

What this means is that the solution field over the mesh is partial information: the other half is knowing which basis $\{\phi_i\}_{i=1}^{N_P}$, where $N_P$ is the number of DoFs, is used. Then your solution field writes

$f_h = \sum f_i \phi_i$

where $f_i$ are the values supplied by the discrete solution field. Now, to compute partial derivatives, you simply differentiate the basis functions:

$\partial_j f_h = \sum_i f_i \partial_j \phi_i$

By convention, the basis is the Lagrange basis. This is by far the most common one. Only if another basis is used would a person note it explicitely. The answer by @ConvexHull gives you the expression in that case.

Gradients at vertices

These bases are defined element-by-element and, if you only have DoFs at the vertices, are bases of affine functions. As such, derivatives are discontinuous accross elements. This means you cannot compute the derivative at a vertex.

To solve this, there are several approaches. One is to enrich the FE space and resolve the new degrees of freedom by $L^2$ projection. The end result, for the gradient at a point $P$ belonging to elements $\mathcal{B}(P)$, is

$\nabla f_{h^+}(P) = \frac{\sum_{K\in\mathcal{B}(P)} |K| \nabla (f_h)_{|K}}{\sum_{K\in\mathcal{B}(P)} |K|}$

with $|K|$ the element measure. As often, an intimidating theoretical setting leads to the mundane: a weighted average. One thing of note is that the gradient at $P$ then involves values at all points belonging to adjacent elements (degree 1 neighbourhood). You can also reconstruct the Hessian, in which case you'll involve the degree 2 neighbourhood (up to neighbours of adjacent elements).

For a more Finite-Differences-friendly approach, you could also approximate the gradient by least-squares (or other norm) regression. What is the gradient if not the linear component of the Taylor expansion of $f$ at $P$? Then seek $f$ of the form

$f(x) = f(P) + b^T (x-P)$

and minimize the energy

$\sum_{i\ \text{s.t.} \ P_i \ \text{neighbour of} \ P}||f(x_i) - f_i||^2$

You could also fit the more lax $f(x) = a + b^T (x-P)$ to possibly improve the gradient, then add the term $||f_{i_0} - a||^2$ with $f_{i_0}$ the DoF at $P$.