# implementing higher order derivatives for finite element

I am implementing higher order derivatives for FEM. Example, to solve a Poisson problem, biharmonic or triharmonic PDE one needs first, second or third order derivatives respectively.

As usually done in FEM, there is a mapping from a reference element to a physical element e.g. $F : \hat{K} \rightarrow K$ with the associated coordinates $(\hat{x}, \hat{y})$ and $(x,y).$

If we assume a linear affine map on a quadrilateral mesh, then for a function $u,$ the derivative is easily computed and given by $\nabla u = J^{-T} \nabla \hat{u}.$ Now, my aim is to compute $\nabla^2 u, \nabla^3 u,...$ possibly $\nabla^n u.$

I will appreciate any ideas on how to compute these derivatives.

• You do, however, know that the common elements only have continuous, but not continuously differentiable shape functions, right? You can't compute the second derivative of a shape function -- it has delta functions on cell faces. Commented Oct 14, 2016 at 15:59
• Not all common elements are even continuous, if it comes to that. Commented Oct 14, 2016 at 16:10
• Sure sure. I just wanted to point out that there are problems at cell interfaces even for the simplest and most common cases. Commented Oct 14, 2016 at 21:47
• I am aware of many occasions where you would need to evaluate second, third or fourth order derivatives elementwise. Second derivatives arise in simple adaptive refinement of second order problems based on residual estimators (same thing with stabilized mixed methods). Third order derivatives arise in non-conforming methods for sixth-order problems (i.e. gradient elasticity). Fourth order derivatives arise at least when doing adaptive refinement for plate (biharmonic) problems. So in my opinion it makes sense to support these kinds of things in a FE code.
– knl
Commented Oct 18, 2016 at 12:48
• @kln, you are very right. I will appreciate a formula for these derivatives. any reference will be good. Commented Oct 25, 2016 at 11:31

Changing notation slightly to use $\bf x$ for physical coordinates and $\bf \xi$ for the reference coordinate (because I have a nasty habit of leaving tildes off variables where I shouldn't) you have in component form that

$\frac{\partial}{\partial x_{i}}u\left(\bf{x}\right)=\sum_{j}\frac{\partial\xi_{j}}{\partial x_{i}}\frac{\partial}{\partial\xi_{j}}u\left(\bf{\xi}\right)$

so using the definition

$\frac{\partial^2 f}{\partial x_i^2} := \frac{\partial}{\partial x_i}\frac{\partial f}{\partial x_i}$

we get

$\frac{\partial^2 u}{\partial x_i^2} = \frac{\partial}{\partial x_i}\frac{\partial u}{\partial x_i} = \frac{\partial}{\partial x_i} \frac{\partial}{\partial x_{i}}u\left(\bf{x}\right)= \frac{\partial}{\partial x_i} \sum_{j}\frac{\partial\xi_{j}}{\partial x_{i}}\frac{\partial}{\partial\xi_{j}}u\left(\bf{\xi}\right) = \sum_k \frac{\partial\xi_{k}}{\partial x_{i}}\frac{\partial}{\partial\xi_{k}}\sum_{j}\frac{\partial\xi_{j}}{\partial x_{i}}\frac{\partial}{\partial\xi_{j}}u\left(\bf{\xi}\right)$.

Assuming your mapping is linear, we can take the derivative through the summation and write

$\frac{\partial^2 u}{\partial x_i^2} = \sum_{j,k} \frac{\partial\xi_{k}}{\partial x_{i}}\frac{\partial\xi_{j}}{\partial x_{i}}\frac{\partial^2 }{\partial \xi_k \partial\xi_{j}}u\left(\bf{\xi}\right)$.

Hopefully you see the connection to what you're already doing with first derivatives.

• This is correct for linear mappings (e.g., on triangles or tetrahedra) but the formulas become much more complicated if the mapping is not linear (e.g., on quadrilaterals or hexahedra). Commented Oct 14, 2016 at 15:58
• True. In that case the assumption in the last line fails, and we have to differentiate the terms in the Jacobian, evaluated at the quadrature points, which gets tedious in general. Commented Oct 14, 2016 at 16:26
• Actually, you have to differentiate the elements of the inverse of the Jacobian. That gets even more tedious ;-) Commented Oct 14, 2016 at 21:47
• Thanks Wolfgang and origimo. Truly, I'm using an affine linear map for quadrilateral meshes. i guess i should have added this earlier. Commented Oct 15, 2016 at 10:49

I think we can compute it using the chain rule, similar to how we compute the Jacobian. For example, if we want to compute the second derivative of $\Phi$ w.r.t physical coordinates in 2D, i.e $\dfrac{\partial^2 \Phi}{\partial x^2}$, $\dfrac{\partial^2 \Phi}{\partial y^2}$ and $\dfrac{\partial^2 \Phi}{\partial x \partial y}$, we can first write

\begin{align} \dfrac{\partial^2 \Phi}{\partial \xi^2} = \dfrac{\partial}{\partial \xi} \left( \dfrac{\partial \Phi}{\partial \xi} \right) &= \dfrac{\partial }{\partial \xi} \left( \dfrac{\partial \Phi}{\partial x} \dfrac{\partial x}{\partial \xi} + \dfrac{\partial \Phi}{\partial y} \dfrac{\partial y}{\partial \xi} \right) \nonumber \\ &=\left[ \dfrac{\partial }{\partial \xi} \left( \dfrac{\partial \Phi}{\partial x} \right) \dfrac{\partial x}{\partial \xi} + \dfrac{\partial \Phi}{\partial x} \dfrac{\partial^2 x}{\partial \xi^2} \right] + \left[ \dfrac{\partial }{\partial \xi} \left( \dfrac{\partial \Phi}{\partial y} \right) \dfrac{\partial y}{\partial \xi} + \dfrac{\partial \Phi}{\partial y} \dfrac{\partial^2 y}{\partial \xi^2} \right] \tag{1}\label{eq:d2_phi_d_xi_2-1} \end{align}

The terms $\dfrac{\partial }{\partial \xi} \left( \dfrac{\partial \Phi}{\partial x} \right)$ and $\dfrac{\partial }{\partial \xi} \left( \dfrac{\partial \Phi}{\partial y} \right)$ can be rewritten in physical coordinates, which yields \begin{align} \dfrac{\partial }{\partial \xi} \left( \dfrac{\partial \Phi}{\partial x} \right) = \dfrac{\partial }{\partial x} \left( \dfrac{\partial \Phi}{\partial x} \right) \dfrac{\partial x}{\partial \xi} + \dfrac{\partial }{\partial y} \left( \dfrac{\partial \Phi}{\partial x} \right) \dfrac{\partial y}{\partial \xi} = \dfrac{\partial^2 \Phi}{\partial x^2} \dfrac{\partial x}{\partial \xi} + \dfrac{\partial^2 \Phi}{\partial x \partial y} \dfrac{\partial y}{\partial \xi} \tag{2a}\label{eq:d2_phi_d_x_d_xi}\\ \dfrac{\partial }{\partial \xi} \left( \dfrac{\partial \Phi}{\partial y} \right) = \dfrac{\partial }{\partial x} \left( \dfrac{\partial \Phi}{\partial y} \right) \dfrac{\partial x}{\partial \xi} + \dfrac{\partial }{\partial y} \left( \dfrac{\partial \Phi}{\partial y} \right) \dfrac{\partial y}{\partial \xi} = \dfrac{\partial^2 \Phi}{\partial x \partial y} \dfrac{\partial x}{\partial \xi} + \dfrac{\partial^2 \Phi}{\partial y^2} \dfrac{\partial y}{\partial \xi} \tag{2b}\label{eq:d2_phi_d_y_d_xi} \end{align}

Replace \eqref{eq:d2_phi_d_x_d_xi} and \eqref{eq:d2_phi_d_y_d_xi} to \eqref{eq:d2_phi_d_xi_2-1} yields $$\dfrac{\partial^2 \Phi}{\partial \xi^2} = \left[ \dfrac{\partial^2 \Phi}{\partial x^2} \left( \dfrac{\partial x}{\partial \xi} \right)^2 + \dfrac{\partial^2 \Phi}{\partial x \partial y} \dfrac{\partial x}{\partial \xi} \dfrac{\partial y}{\partial \xi} + \dfrac{\partial \Phi}{\partial x} \dfrac{\partial^2 x}{\partial \xi^2} \right] + \left[ \dfrac{\partial^2 \Phi}{\partial x \partial y} \dfrac{\partial x}{\partial \xi} \dfrac{\partial y}{\partial \xi} + \dfrac{\partial^2 \Phi}{\partial y^2} \left( \dfrac{\partial y}{\partial \xi} \right)^2 + \dfrac{\partial \Phi}{\partial y} \dfrac{\partial^2 y}{\partial \xi^2} \right] \tag{3}\label{eq:d2_phi_d_xi_2-2}$$

Rearrange \eqref{eq:d2_phi_d_xi_2-2} to localize the second derivatives terms w.r.t physical coordinates $$\dfrac{\partial^2 \Phi}{\partial \xi^2} - \dfrac{\partial \Phi}{\partial x} \dfrac{\partial^2 x}{\partial \xi^2} - \dfrac{\partial \Phi}{\partial y} \dfrac{\partial^2 y}{\partial \xi^2} = \dfrac{\partial^2 \Phi}{\partial x^2} \left( \dfrac{\partial x}{\partial \xi} \right)^2 + \dfrac{\partial^2 \Phi}{\partial y^2} \left( \dfrac{\partial y}{\partial \xi} \right)^2 + 2\dfrac{\partial^2 \Phi}{\partial x \partial y} \dfrac{\partial x}{\partial \xi} \dfrac{\partial y}{\partial \xi} \tag{4}\label{eq:local_2d_2nd_derivative-1}$$

In the same manner, we can also derive for $\dfrac{\partial^2 \Phi}{\partial \eta^2}$ and yield $$\dfrac{\partial^2 \Phi}{\partial \eta^2} - \dfrac{\partial \Phi}{\partial x} \dfrac{\partial^2 x}{\partial \eta^2} - \dfrac{\partial \Phi}{\partial y} \dfrac{\partial^2 y}{\partial \eta^2} = \dfrac{\partial^2 \Phi}{\partial x^2} \left( \dfrac{\partial x}{\partial \eta} \right)^2 + \dfrac{\partial^2 \Phi}{\partial y^2} \left( \dfrac{\partial y}{\partial \eta} \right)^2 + 2\dfrac{\partial^2 \Phi}{\partial x \partial y} \dfrac{\partial x}{\partial \eta} \dfrac{\partial y}{\partial \eta}\tag{5} \label{eq:local_2d_2nd_derivative-2}$$

Derivation of $\dfrac{\partial^2 \Phi}{\partial \xi \partial \eta}$ yield \begin{align} \dfrac{\partial^2 \Phi}{\partial \xi \partial \eta} &= \dfrac{\partial}{\partial \xi} \left( \dfrac{\partial \Phi}{\partial \eta} \right) = \dfrac{\partial }{\partial \xi} \left( \dfrac{\partial \Phi}{\partial x} \dfrac{\partial x}{\partial \eta} + \dfrac{\partial \Phi}{\partial y} \dfrac{\partial y}{\partial \eta} \right) \nonumber \\ &=\left[ \dfrac{\partial }{\partial \xi} \left( \dfrac{\partial \Phi}{\partial x} \right) \dfrac{\partial x}{\partial \eta} + \dfrac{\partial \Phi}{\partial x} \dfrac{\partial^2 x}{\partial \xi \partial \eta} \right] + \left[ \dfrac{\partial }{\partial \xi} \left( \dfrac{\partial \Phi}{\partial y} \right) \dfrac{\partial y}{\partial \eta} + \dfrac{\partial \Phi}{\partial y} \dfrac{\partial^2 y}{\partial \xi \partial \eta} \right] \nonumber \\ &= \left[ \dfrac{\partial^2 \Phi}{\partial x^2} \dfrac{\partial x}{\partial \xi} \dfrac{\partial x}{\partial \eta} + \dfrac{\partial^2 \Phi}{\partial x \partial y} \dfrac{\partial y}{\partial \xi} \dfrac{\partial x}{\partial \eta} + \dfrac{\partial \Phi}{\partial x} \dfrac{\partial^2 x}{\partial \xi \partial \eta} \right] + \left[ \dfrac{\partial^2 \Phi}{\partial x \partial y} \dfrac{\partial x}{\partial \xi} \dfrac{\partial y}{\partial \eta} + \dfrac{\partial^2 \Phi}{\partial y^2} \dfrac{\partial y}{\partial \xi} \dfrac{\partial y}{\partial \eta} + \dfrac{\partial \Phi}{\partial y} \dfrac{\partial^2 y}{\partial \xi \partial \eta} \right] \nonumber \end{align}

Or equivalently $$\dfrac{\partial^2 \Phi}{\partial \xi \partial \eta} - \dfrac{\partial \Phi}{\partial x} \dfrac{\partial^2 x}{\partial \xi \partial \eta} - \dfrac{\partial \Phi}{\partial y} \dfrac{\partial^2 y}{\partial \xi \partial \eta} = \dfrac{\partial^2 \Phi}{\partial x^2} \dfrac{\partial x}{\partial \xi} \dfrac{\partial x}{\partial \eta} + \dfrac{\partial^2 \Phi}{\partial y^2} \dfrac{\partial y}{\partial \xi} \dfrac{\partial y}{\partial \eta} + \dfrac{\partial^2 \Phi}{\partial x \partial y} \left( \dfrac{\partial x}{\partial \xi} \dfrac{\partial y}{\partial \eta} + \dfrac{\partial x}{\partial \eta} \dfrac{\partial y}{\partial \xi} \right) \tag{6}\label{eq:local_2d_2nd_derivative-3}$$

Grouping \eqref{eq:local_2d_2nd_derivative-1}, \eqref{eq:local_2d_2nd_derivative-2} and \eqref{eq:local_2d_2nd_derivative-3} into the matrix form, we have $$\begin{bmatrix} \left( \dfrac{\partial x}{\partial \xi} \right)^2 & \left( \dfrac{\partial y}{\partial \xi} \right)^2 & 2 \dfrac{\partial x}{\partial \xi} \dfrac{\partial y}{\partial \xi} \\[1em] \left( \dfrac{\partial x}{\partial \eta} \right)^2 & \left( \dfrac{\partial y}{\partial \eta} \right)^2 & 2 \dfrac{\partial x}{\partial \eta} \dfrac{\partial y}{\partial \eta} \\[1em] \dfrac{\partial x}{\partial \xi} \dfrac{\partial x}{\partial \eta} & \dfrac{\partial y}{\partial \xi} \dfrac{\partial y}{\partial \eta} & \dfrac{\partial x}{\partial \xi} \dfrac{\partial y}{\partial \eta} + \dfrac{\partial x}{\partial \eta} \dfrac{\partial y}{\partial \xi} \end{bmatrix} \begin{bmatrix} \dfrac{\partial^2 \Phi}{\partial x^2} \\[1em] \dfrac{\partial^2 \Phi}{\partial y^2} \\[1em] \dfrac{\partial^2 \Phi}{\partial x \partial y} \end{bmatrix} =\begin{bmatrix} \dfrac{\partial^2 \Phi}{\partial \xi^2} - \dfrac{\partial \Phi}{\partial x} \dfrac{\partial^2 x}{\partial \xi^2} - \dfrac{\partial \Phi}{\partial y} \dfrac{\partial^2 y}{\partial \xi^2} \\[1em] \dfrac{\partial^2 \Phi}{\partial \eta^2} - \dfrac{\partial \Phi}{\partial x} \dfrac{\partial^2 x}{\partial \eta^2} - \dfrac{\partial \Phi}{\partial y} \dfrac{\partial^2 y}{\partial \eta^2} \\[1em] \dfrac{\partial^2 \Phi}{\partial \xi \partial \eta} - \dfrac{\partial \Phi}{\partial x} \dfrac{\partial^2 x}{\partial \xi \partial \eta} - \dfrac{\partial \Phi}{\partial y} \dfrac{\partial^2 y}{\partial \xi \partial \eta} \end{bmatrix}$$

Solve the above system will give us the answer.

The linear system in 3D has the form

$$\begin{bmatrix} \left( \dfrac{\partial x}{\partial \xi} \right)^2 & \left( \dfrac{\partial y}{\partial \xi} \right)^2 & \left( \dfrac{\partial z}{\partial \xi} \right)^2 & 2 \dfrac{\partial x}{\partial \xi} \dfrac{\partial y}{\partial \xi} & 2 \dfrac{\partial y}{\partial \xi} \dfrac{\partial z}{\partial \xi} & 2 \dfrac{\partial x}{\partial \xi} \dfrac{\partial z}{\partial \xi} \\[1em] \left( \dfrac{\partial x}{\partial \eta} \right)^2 & \left( \dfrac{\partial y}{\partial \eta} \right)^2 & \left( \dfrac{\partial z}{\partial \eta} \right)^2 & 2 \dfrac{\partial x}{\partial \eta} \dfrac{\partial y}{\partial \eta} & 2 \dfrac{\partial y}{\partial \eta} \dfrac{\partial z}{\partial \eta} & 2 \dfrac{\partial x}{\partial \eta} \dfrac{\partial z}{\partial \eta} \\[1em] \left( \dfrac{\partial x}{\partial \zeta} \right)^2 & \left( \dfrac{\partial y}{\partial \zeta} \right)^2 & \left( \dfrac{\partial z}{\partial \zeta} \right)^2 & 2 \dfrac{\partial x}{\partial \zeta} \dfrac{\partial y}{\partial \zeta} & 2 \dfrac{\partial y}{\partial \zeta} \dfrac{\partial z}{\partial \zeta} & 2 \dfrac{\partial x}{\partial \zeta} \dfrac{\partial z}{\partial \zeta} \\[1em] \dfrac{\partial x}{\partial \xi} \dfrac{\partial x}{\partial \eta} & \dfrac{\partial y}{\partial \xi} \dfrac{\partial y}{\partial \eta} & \dfrac{\partial z}{\partial \xi} \dfrac{\partial z}{\partial \eta} & \dfrac{\partial x}{\partial \xi} \dfrac{\partial y}{\partial \eta} + \dfrac{\partial x}{\partial \eta} \dfrac{\partial y}{\partial \xi} & \dfrac{\partial y}{\partial \xi} \dfrac{\partial z}{\partial \eta} + \dfrac{\partial y}{\partial \eta} \dfrac{\partial z}{\partial \xi} & \dfrac{\partial x}{\partial \xi} \dfrac{\partial z}{\partial \eta} + \dfrac{\partial x}{\partial \eta} \dfrac{\partial z}{\partial \xi} \\[1em] \dfrac{\partial x}{\partial \eta} \dfrac{\partial x}{\partial \zeta} & \dfrac{\partial y}{\partial \eta} \dfrac{\partial y}{\partial \zeta} & \dfrac{\partial z}{\partial \eta} \dfrac{\partial z}{\partial \zeta} & \dfrac{\partial x}{\partial \eta} \dfrac{\partial y}{\partial \zeta} + \dfrac{\partial x}{\partial \zeta} \dfrac{\partial y}{\partial \eta} & \dfrac{\partial y}{\partial \eta} \dfrac{\partial z}{\partial \zeta} + \dfrac{\partial y}{\partial \zeta} \dfrac{\partial z}{\partial \eta} & \dfrac{\partial x}{\partial \eta} \dfrac{\partial z}{\partial \zeta} + \dfrac{\partial x}{\partial \zeta} \dfrac{\partial z}{\partial \eta} \\[1em] \dfrac{\partial x}{\partial \xi} \dfrac{\partial x}{\partial \zeta} & \dfrac{\partial y}{\partial \xi} \dfrac{\partial y}{\partial \zeta} & \dfrac{\partial z}{\partial \xi} \dfrac{\partial z}{\partial \zeta} & \dfrac{\partial x}{\partial \xi} \dfrac{\partial y}{\partial \zeta} + \dfrac{\partial x}{\partial \zeta} \dfrac{\partial y}{\partial \xi} & \dfrac{\partial y}{\partial \xi} \dfrac{\partial z}{\partial \zeta} + \dfrac{\partial y}{\partial \zeta} \dfrac{\partial z}{\partial \xi} & \dfrac{\partial x}{\partial \xi} \dfrac{\partial z}{\partial \zeta} + \dfrac{\partial x}{\partial \zeta} \dfrac{\partial z}{\partial \xi} \end{bmatrix} \begin{bmatrix} \dfrac{\partial^2 \Phi}{\partial x^2} \\[1em] \dfrac{\partial^2 \Phi}{\partial y^2} \\[1em] \dfrac{\partial^2 \Phi}{\partial z^2} \\[1em] \dfrac{\partial^2 \Phi}{\partial x \partial y} \\[1em] \dfrac{\partial^2 \Phi}{\partial y \partial z} \\[1em] \dfrac{\partial^2 \Phi}{\partial x \partial z} \end{bmatrix} =\begin{bmatrix} \dfrac{\partial^2 \Phi}{\partial \xi^2} - \dfrac{\partial \Phi}{\partial x} \dfrac{\partial^2 x}{\partial \xi^2} - \dfrac{\partial \Phi}{\partial y} \dfrac{\partial^2 y}{\partial \xi^2} - \dfrac{\partial \Phi}{\partial z} \dfrac{\partial^2 z}{\partial \xi^2} \\[1em] \dfrac{\partial^2 \Phi}{\partial \eta^2} - \dfrac{\partial \Phi}{\partial x} \dfrac{\partial^2 x}{\partial \eta^2} - \dfrac{\partial \Phi}{\partial y} \dfrac{\partial^2 y}{\partial \eta^2} - \dfrac{\partial \Phi}{\partial z} \dfrac{\partial^2 z}{\partial \eta^2} \\[1em] \dfrac{\partial^2 \Phi}{\partial \zeta^2} - \dfrac{\partial \Phi}{\partial x} \dfrac{\partial^2 x}{\partial \zeta^2} - \dfrac{\partial \Phi}{\partial y} \dfrac{\partial^2 y}{\partial \zeta^2} - \dfrac{\partial \Phi}{\partial z} \dfrac{\partial^2 z}{\partial \zeta^2} \\[1em] \dfrac{\partial^2 \Phi}{\partial \xi \partial \eta} - \dfrac{\partial \Phi}{\partial x} \dfrac{\partial^2 x}{\partial \xi \partial \eta} - \dfrac{\partial \Phi}{\partial y} \dfrac{\partial^2 y}{\partial \xi \partial \eta} - \dfrac{\partial \Phi}{\partial z} \dfrac{\partial^2 z}{\partial \xi \partial \eta} \\[1em] \dfrac{\partial^2 \Phi}{\partial \eta \partial \zeta} - \dfrac{\partial \Phi}{\partial x} \dfrac{\partial^2 x}{\partial \eta \partial \zeta} - \dfrac{\partial \Phi}{\partial y} \dfrac{\partial^2 y}{\partial \eta \partial \zeta} - \dfrac{\partial \Phi}{\partial z} \dfrac{\partial^2 z}{\partial \eta \partial \zeta} \\[1em] \dfrac{\partial^2 \Phi}{\partial \xi \partial \zeta} - \dfrac{\partial \Phi}{\partial x} \dfrac{\partial^2 x}{\partial \xi \partial \zeta} - \dfrac{\partial \Phi}{\partial y} \dfrac{\partial^2 y}{\partial \xi \partial \zeta} - \dfrac{\partial \Phi}{\partial z} \dfrac{\partial^2 z}{\partial \xi \partial \zeta} \end{bmatrix}$$

(Credit: Thanks C. Clason to edit the label for the equation on Mathjax.)