Gradients on 2D FEM Triangular Mesh

Question

A scalar field f approximated on a triangular mesh using the FEM method can be given as

$f(x,y) \approx \bar{f}(x,y) = \sum_{n=1}^3 \left(\phi_n(x,y) \cdot f_n \right)$

I am of the understanding, then, that the gradient of f can be found inside each triangle with

$\nabla \bar{f}_t = \sum_{n=1}^3 \left(\frac{\partial\phi_n(x,y)}{\partial x} \cdot f_n , \frac{\partial\phi_n(x,y)}{\partial y} \cdot f_n\right)$

Where $f_t$ is the value at the triangle centroids. Since the value at each point $f_n$ is only defined at single points, it acts as a scalar in the derivative.

But, what if we then want to take the gradient of each of the components from the previous solution? One way seems to be to interpolate the gradient results from the triangle centroids back to the given points, and then repeat the process for each component. However, what if this interpolation is unsatisfactory -- is there a way to possibly define the gradient at each point from data at the triangle centroids?

There seems to be a similiar question posed here, but I'm not sure if the solution is relevant -- the answer solution seems to involve going through all the basis functions of adjacent triangles, whereas it seems to me that only one basis function from each adjacent triangle is needed -- the one which has a nonzero value at the observed point.

With this in mind, I've derived the following equation:

$\nabla \bar{f}_i = \sum_{n=1}^{t_i} \left(\left( \frac{\partial\phi_i^n(x,y)}{\partial x}, \frac{\partial\phi_i^n(x,y)}{\partial y} \right) \frac{f_n}{t_i} \right)$

Where i is the observed point, and $t_i$ is the number of triangles which contain the vertex point i.

Unfortunately, it doesn't seem to work, and I'm not sure if averaging over the number of triangles is appropriate. Can someone point me in the right direction?

Answer 1

Let me try a quick and dirty answer. Suppose that you have linear shape functions $\phi_n(x,y)$: if you take the gradient of the shape functions (e.g. for computing the energy norm and the stiffness matrix) you end up with a constant function. (Hence the name CST in continuum mechanics: Constant Strain Triangle, since strain is related to the (deformation) gradient tensor of the displacement field.)

As a consequence if you compute the gradient $\nabla \bar{f}$, this will be piecewise constant on each FE, and is only a crude approximation of $\nabla f$. However Barlow noted super convergence of $\nabla \bar{f}$ in certain interior points of the FE. This observation leads to the following recipe, which is universally adopted in FE post processing programs, when a representation of the gradient $\nabla f$ is required.

1. evaluate the gradient field at numerical quadrature points (since most FE's are numerically integrated): for a CST triangle this is a single point at the centroid.

2. extrapolate the gradient from the quadrature points to the element nodes. For a CST since you have a single quadrature point, no extrapolation is necessary: simply take the centroid value.

3. at each node take the average of the gradient value extrapolated from all elements that have that node in common.

4. use the shape functions to interpolate (inside the FEs) the nodal averaged value.

This recipe is indeed very similar to what I understand from the OP question: so I think his basic idea is correct.

The theory for this technique (nodal averaging over a path of elements extrapolating from the super convergent points for the flux) can be found in O.C. Zienkiewicz and J.Z. Zhu, The Superconvergent patch recovery and a posteriori error estimators. Part 1. The recovery technique, Int. J. Numer. Methods Eng., 33, 1331-1364 (1992).