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. – Wolfgang Bangerth Oct 14 '16 at 15:59
• Not all common elements are even continuous, if it comes to that. – origimbo Oct 14 '16 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. – Wolfgang Bangerth Oct 14 '16 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 Oct 18 '16 at 12:48
• @kln, you are very right. I will appreciate a formula for these derivatives. any reference will be good. – uli.xu Oct 25 '16 at 11:31

2 Answers

You appear to have pretty much answered your own question.

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). – Wolfgang Bangerth Oct 14 '16 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. – origimbo Oct 14 '16 at 16:26
• Actually, you have to differentiate the elements of the inverse of the Jacobian. That gets even more tedious ;-) – Wolfgang Bangerth Oct 14 '16 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. – uli.xu Oct 15 '16 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.)