implementing higher order derivatives for finite element

Question

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.

Answer 1

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.

Answer 2

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 \begin{equation} \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} \end{equation}

Rearrange \eqref{eq:d2_phi_d_xi_2-2} to localize the second derivatives terms w.r.t physical coordinates \begin{equation} \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} \end{equation}

In the same manner, we can also derive for $\dfrac{\partial^2 \Phi}{\partial \eta^2}$ and yield \begin{equation} \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} \end{equation}

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 \begin{equation} \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} \end{equation}

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{equation} \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} \end{equation}

Solve the above system will give us the answer.

The linear system in 3D has the form

\begin{equation} \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} \end{equation}

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