Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It only takes a minute to sign up.
Sign up to join this community
I have a 2-D mesh of triangles and I have a scalar function $f(x,y)$ defined at all the vertices of this mesh.
I want to accurately estimate the values of $\frac{\partial f}{\partial x}$ and $\frac{\partial f}{\partial y}$.
I am aware how to do this when dealing with finite difference methods on the standard orthogonal grid, using the standard formulas. How will I do it on the 2-D mesh?
If you have a triangular mesh with data defined at each vertex, you can of course associate a linear function with each element that interpolates the data values at the vertices. This also provides you with an estimate of the gradient on each cell by just computing the gradient of the linear interpolant. Of course, it is only an estimate because you are compute the gradient of the interpolant of $f$, not of $f$ itself.
However, it is possible to improve the accuracy of this estimate. For example, you could use a process such as the ZZ estimator (short, for Zienkiewicz and Zhou). If you google for "ZZ estimator" you will find quite a number of links to pages that describe how this is done. (Typically, these pages talk about stresses, strains and displacements; in your context, imagine stress=strain=$\nabla f$, and displacement=$f$.)
You can compute differential operators for this meshes using information around the neighborhood of each vertex (discrete in this case). In this reference the authors compute differential-geometry operators in the 1-star neighborhood, i.e., the cells around each vertex.
Another option is to fit a paraboloid around each vertex and analytically calculate the derivatives.
The most conceptually simple approach is to select an interpolation scheme, interpolate through the data points that you have, then compute the derivative of the interpolant. This is a natural generalization of the typical finite-difference differentiation formulas.
While you don't mention your exact requirements (in particular: do you need the derivatives at the mesh points, or within the mesh triangles?), the simplest approach is to interpolate the function values linearly in each triangle, and estimate the derivative as the derivative of the linear interpolant; the interpolant will have a piecewise constant derivative, and be discontinuous at the mesh points and edges, where you can take linear combinations of adjacent interpolants.
You can also look at other inteprolation methods, like radial basis functions.
Also note that many software packages support 2d interpolation themselves, in which case it might not be necessary to do the interpolation step yourself. If a library can interpolate a function, but doesn't calculate the interpolant's derivatives, you can recover the derivatives by using a finite-difference approximation on the interpolant.
I usually use gradient averaging since it has some superconvergence properties under certain regularity assumptions and I always tend to have some sort of finite element code lying around when I want to do this. Let $V_h$ be for example the piecewise linear finite element space defined on the mesh of the domain $\Omega$. The idea is to consider the projection problem: find $w \in V_h$ such that $$\int_\Omega w v \,\mathrm{d}x = \int_\Omega \frac{\partial f}{\partial x} v\,\mathrm{d}x, \quad \forall v \in V_h.$$ The matrix formulation reads $$M\boldsymbol{w}=C\boldsymbol{f},$$ where $M$ is the typical mass matrix, $\boldsymbol{f}$ is the vector of nodal values and $C_{ij}=\int_\Omega \frac{\partial \varphi_j}{\partial x} \varphi_i \,\mathrm{d}x$ where $\varphi_i$ is the basis of $V_h$. Then you simply solve $$\boldsymbol{w}=M^{-1}C\boldsymbol{f}.$$
If you don't have a finite element code lying around, you can implement the equivalent projection by computing the piecewise constant derivative in each neighboring element of a node, taking their weighted average with the areas of the triangles as weights and let that average value represent the derivative at the node.
