calculating the Laplacian of the field variable in estimating the local residual error in the finite element method

Question

to perform adaptive refinement in the finite element method according to the explicit residual method, the quantity $$\eta_K^2=h_K^2\left\lVert r\right\rVert_{L_2(K)}^2+h_K\left\lVert R\right\rVert_{L_2(\partial K)}^2$$ must be calculated for all elements of the mesh. For the equation $-\nabla^2u=f$, the norm $\left\lVert r\right\rVert_{L_2(K)}^2$ is calculated by $$\left\lVert r\right\rVert_{L_2(K)}^2=\int_{\Omega^K}\left[\nabla^2u_h+f\right]^2 d\Omega^K$$ in which the Laplacian of the field variable, i.e. $\nabla^2u_h$, appears. Calculating $\nabla^2u_h$ in each element requires second derivates of shape functions and hence is so expensive, particularly for higher order hierarchic shape functions.

How $\nabla^2u_h$ could be calculated (or estimated) for hierarchic shape functions?

Answer 1

Let $\mathbf{x} = \mathbf{F}(\hat{\mathbf{x}})$ be the mapping from the reference cell to the real cell. Then your solution at a point $ \mathbf{x}$ is given by $$ u_h( \mathbf{x}) = \sum_j U_j \varphi_j( \mathbf{x}) = \sum_j U_j \hat\varphi_j( \mathbf{F}^{-1} \mathbf{x})= \sum_j U_j \hat\varphi_j(\hat{\mathbf{x}}). $$ If you want to evaluate its derivatives, you get $$ \nabla_x u_h( \mathbf{x}) = \sum_j U_j \nabla_x \varphi_j( \mathbf{x}) = \sum_j U_j \nabla_x \hat\varphi_j( \mathbf{F}^{-1} \mathbf{x}) = \sum_j U_j J^{-1} \nabla_{\hat x} \hat\varphi_j(\hat{\mathbf{x}}). $$ Here, $J^{-1} = [\nabla_{\hat x} \mathbf F(\hat{\mathbf x})]^{-1}$ is the inverse of the Jacobian of the mapping.

By the same token, you get $$ \nabla_x^2 u_h( \mathbf{x}) = \sum_j U_j J^{-1} \nabla_{\hat x} \left(J^{-1} \nabla_{\hat x} \hat\varphi_j(\hat{\mathbf{x}})\right). $$ Applying the product rule, you have $$ \nabla_x^2 u_h( \mathbf{x}) = \sum_j U_j \left\{ J^{-1} J^{-1} \nabla_{\hat x}^2 \hat\varphi_j(\hat{\mathbf{x}}) + J^{-1} \left[\nabla_{\hat x} J^{-1}\right] \nabla_{\hat x} \hat\varphi_j(\hat{\mathbf{x}}) \right\}. $$ In this formula, you need to pay a bit of attention to what indices you are summing over in all of these tensor contractions, details that I will skip here. The key step that is left to recognize that $$ \nabla_{\hat x} J^{-1} = -J^{-1} [\nabla_{\hat x} J] J^{-1}, $$ a formula you obtain by recognizing that $\nabla_{\hat x}[JJ^{-1}] = [\nabla_{\hat x}J] J^{-1} + J [\nabla_{\hat x} J^{-1}] = 0$.

In other words, what you need to compute in order to evaluate $\nabla_x u_h(\mathbf x)$ are:

• $J$, the Jacobian of the mapping. This is easy because the mapping $\mathbf F$ is generally a low-order polynomial
• $J^{-1}$ its inverse. Computing the inverse of a $2\times 2$ or $3\times 3$ matrix is also relatively easy.
• $\nabla_{\hat x}J$, the second derivatives of the mapping. This is again easy because the mapping $\mathbf F$ is generally a low-order polynomial
• $\nabla_{\hat x} J^{-1} = -J^{-1} [\nabla_{\hat x} J] J^{-1}$. Just some tensor contractions.
• $\nabla_{\hat x}^2 \hat\varphi_j(\hat{\mathbf x})$. This is again simple because you just have to compute second derivatives of a polynomial.

All of this is not exactly cheap, but it's also not super expensive. A lot of it can be pretabulated on the reference cell, for example.

Answer 2

assume $(r,s)$ and $(x,y)$ are the coordinates in the reference and actual element, respectively. the shape functions for mapping and interpolation are defined as $$ \mathbf{M}= \begin{bmatrix} M_1(r,s) & M_2(r,s) & \cdots & M_m(r,s) \end{bmatrix}^T $$$$ \mathbf{I}= \begin{bmatrix} I_1(r,s) & I_2(r,s) & \cdots & I_i(r,s) \end{bmatrix}^T $$ matrix of first derivatives is $$ \begin{bmatrix}\dfrac{\partial}{\partial r}\\ \dfrac{\partial}{\partial s}\end{bmatrix}= \begin{bmatrix} \dfrac{\partial x}{\partial r}&\dfrac{\partial y}{\partial r}\\ \dfrac{\partial x}{\partial s}&\dfrac{\partial y}{\partial s} \end{bmatrix} \begin{bmatrix}\dfrac{\partial}{\partial x}\\ \dfrac{\partial}{\partial y}\end{bmatrix} $$ from which the matrix of second derivatives can be obtained $$ \begin{bmatrix}\dfrac{\partial^2}{\partial r^2}\\ \dfrac{\partial^2}{\partial s^2}\\ \dfrac{\partial^2}{\partial r\partial s}\end{bmatrix} = \begin{bmatrix}\dfrac{\partial^2 x}{\partial r^2}&\dfrac{\partial^2 y}{\partial r^2}\\ \dfrac{\partial^2 x}{\partial s^2}&\dfrac{\partial^2 y}{\partial s^2}\\ \dfrac{\partial^2 x}{\partial r\partial s}&\dfrac{\partial^2 y}{\partial r\partial s}\end{bmatrix} \begin{bmatrix}\dfrac{\partial}{\partial x}\\ \dfrac{\partial}{\partial y}\end{bmatrix} + \begin{bmatrix}\left(\dfrac{\partial x}{\partial r}\right)^2& \left(\dfrac{\partial y}{\partial r}\right)^2& 2\dfrac{\partial x}{\partial r}\dfrac{\partial y}{\partial r}\\ \left(\dfrac{\partial x}{\partial s}\right)^2& \left(\dfrac{\partial y}{\partial s}\right)^2& 2\dfrac{\partial x}{\partial s}\dfrac{\partial y}{\partial s}\\ \dfrac{\partial x}{\partial r}\dfrac{\partial x}{\partial s}& \dfrac{\partial y}{\partial r}\dfrac{\partial y}{\partial s}& \dfrac{\partial x}{\partial s}\dfrac{\partial y}{\partial r}+ \dfrac{\partial x}{\partial r}\dfrac{\partial y}{\partial s}\end{bmatrix} \begin{bmatrix}\dfrac{\partial^2}{\partial x^2}\\ \dfrac{\partial^2}{\partial y^2}\\ \dfrac{\partial^2}{\partial x\partial y}\end{bmatrix} $$ by mapping $x$ and $y$ with mapping shape functions we have $$ x=\mathbf{M}^T\mathbf{X}^e\quad\Rightarrow\quad \frac{\partial x}{\partial u}=\frac{\partial \mathbf{M}^T}{\partial u}\mathbf{X}^e,\quad \frac{\partial^2 x}{\partial u^2}=\frac{\partial^2\mathbf{M}^T}{\partial u^2}\mathbf{X}^e,\quad \frac{\partial^2 x}{\partial u\partial v}=\frac{\partial^2\mathbf{M}^T}{\partial u\partial v}\mathbf{X}^e,\quad u,v=r,s $$$$ y=\mathbf{M}^T\mathbf{Y}^e\quad\Rightarrow\quad \frac{\partial y}{\partial u}=\frac{\partial \mathbf{M}^T}{\partial u}\mathbf{Y}^e,\quad \frac{\partial^2 y}{\partial u^2}=\frac{\partial^2\mathbf{M}^T}{\partial u^2}\mathbf{Y}^e,\quad \frac{\partial^2 y}{\partial u\partial v}=\frac{\partial^2\mathbf{M}^T}{\partial u\partial v}\mathbf{Y}^e,\quad u,v=r,s $$ where $\mathbf{X}^e$ and $\mathbf{Y}^e$ are column matrices of actual coordinates of the element $e$.

using mapped $x$ and $y$ and their derivatives in the matrix of first and second derivatives yields $$ \begin{bmatrix}\dfrac{\partial^2\mathbf{I}^T}{\partial x^2}\\ \dfrac{\partial^2\mathbf{I}^T}{\partial y^2}\\ \dfrac{\partial^2\mathbf{I}^T}{\partial x\partial y}\end{bmatrix} = \begin{bmatrix}\left(\dfrac{\partial\mathbf{M}^T}{\partial r}\mathbf{X}^e\right)^2& \left(\dfrac{\partial\mathbf{M}^T}{\partial r}\mathbf{Y}^e\right)^2& 2\left(\dfrac{\partial\mathbf{M}^T}{\partial r}\mathbf{X}^e\right)\left(\dfrac{\partial\mathbf{M}^T}{\partial r}\mathbf{Y}^e\right)\\ \left(\dfrac{\partial\mathbf{M}^T}{\partial s}\mathbf{X}^e\right)^2& \left(\dfrac{\partial\mathbf{M}^T}{\partial s}\mathbf{Y}^e\right)^2& 2\left(\dfrac{\partial\mathbf{M}^T}{\partial s}\mathbf{X}^e\right)\left(\dfrac{\partial\mathbf{M}^T}{\partial s}\mathbf{Y}^e\right)\\ \left(\dfrac{\partial\mathbf{M}^T}{\partial r}\mathbf{X}^e\right)\left(\dfrac{\partial\mathbf{M}^T}{\partial s}\mathbf{X}^e\right)& \left(\dfrac{\partial\mathbf{M}^T}{\partial r}\mathbf{Y}^e\right)\left(\dfrac{\partial\mathbf{M}^T}{\partial s}\mathbf{Y}^e\right)& \left(\dfrac{\partial\mathbf{M}^T}{\partial s}\mathbf{X}^e\right)\left(\dfrac{\partial\mathbf{M}^T}{\partial r}\mathbf{Y}^e\right)+ \left(\dfrac{\partial\mathbf{M}^T}{\partial r}\mathbf{X}^e\right)\left(\dfrac{\partial\mathbf{M}^T}{\partial s}\mathbf{Y}^e\right) \end{bmatrix}^{-1} \left( \begin{bmatrix}\dfrac{\partial^2\mathbf{I}^T}{\partial r^2}\\ \dfrac{\partial^2\mathbf{I}^T}{\partial s^2}\\ \dfrac{\partial^2\mathbf{I}^T}{\partial r\partial s}\end{bmatrix} - \begin{bmatrix}\dfrac{\partial^2 \mathbf{M}^T}{\partial r^2}\mathbf{X}^e&\dfrac{\partial^2 \mathbf{M}^T}{\partial r^2}\mathbf{Y}^e\\ \dfrac{\partial^2 \mathbf{M}^T}{\partial s^2}\mathbf{X}^e&\dfrac{\partial^2 \mathbf{M}^T}{\partial s^2}\mathbf{Y}^e\\ \dfrac{\partial^2 \mathbf{M}^T}{\partial r\partial s}\mathbf{X}^e&\dfrac{\partial^2 \mathbf{M}^T}{\partial r\partial s}\mathbf{Y}^e\end{bmatrix} \begin{bmatrix} \dfrac{\partial \mathbf{M}^T}{\partial r}\mathbf{X}^e&\dfrac{\partial \mathbf{M}^T}{\partial r}\mathbf{Y}^e\\ \dfrac{\partial \mathbf{M}^T}{\partial s}\mathbf{X}^e&\dfrac{\partial \mathbf{M}^T}{\partial s}\mathbf{Y}^e \end{bmatrix}^{-1} \begin{bmatrix}\dfrac{\partial\mathbf{I}^T}{\partial r}\\\dfrac{\partial\mathbf{I}^T}{\partial s}\end{bmatrix} \right) $$ the Laplacian of the field variable could be calculated from the above formula in each element.

update

for example, for the four-node quadrilateral element with linear interpolation shape functions $$ \mathbf{M}=\mathbf{I}= \begin{bmatrix} (1-r)(1-s)/4&(1+r)(1-s)/4&(1+r)(1+s)/4&(1-r)(1+s)/4 \end{bmatrix}^T\equiv\mathbf{N} $$ and coordinates $$ \mathbf{X}^e=\begin{bmatrix}1&4&2.5&1.2\end{bmatrix}^T $$ $$ \mathbf{Y}^e=\begin{bmatrix}2&1&4&3\end{bmatrix}^T $$ the second derivative matrix would be $$ \begin{bmatrix}\dfrac{\partial^2}{\partial x^2}\\ \dfrac{\partial^2}{\partial y^2}\\ \dfrac{\partial^2}{\partial x\partial y}\end{bmatrix} \mathbf{N}^T = \begin{bmatrix}\left(\dfrac{17}{40}s-\dfrac{43}{40}\right)^2& \dfrac{1}{4}s^2& -\dfrac{17}{40}s^2+\dfrac{43}{40}s\\ \left(\dfrac{17}{40}r+\dfrac{13}{40}\right)^2& \left(\dfrac{1}{2}r+1\right)^2& -\left(\dfrac{1}{2}r+1\right)\left(\dfrac{17}{20}r+\dfrac{13}{20}\right)\\ \left(\dfrac{17}{40}r+\dfrac{13}{40}\right)\left(\dfrac{17}{40}s-\dfrac{43}{40}\right)& \dfrac{1}{2}s\left(\dfrac{1}{2}r+1\right)& -\dfrac{1}{2}s\left(\dfrac{17}{40}r+\dfrac{13}{40}\right)+\left(\dfrac{1}{2}r+1\right)\left(\dfrac{17}{40}s-\dfrac{43}{40}\right) \end{bmatrix}^{-1} \left( \begin{bmatrix}0&0&0&0\\0&0&0&0\\\dfrac{1}{4}&\dfrac{-1}{4}&\dfrac{1}{4}&\dfrac{-1}{4}\end{bmatrix} - \begin{bmatrix}0&0\\0&0\\-\dfrac{17}{40}&\dfrac{1}{2}\end{bmatrix} \begin{bmatrix} \dfrac{43}{40}-\dfrac{17}{40}s&\dfrac{1}{2}s\\ -\dfrac{17}{40}r-\dfrac{13}{40}&\dfrac{1}{2}r+1 \end{bmatrix}^{-1} \begin{bmatrix} \dfrac{1}{4}s-\dfrac{1}{4}&\dfrac{1}{4}-\dfrac{1}{4}s&\dfrac{1}{4}s+\dfrac{1}{4}&-\dfrac{1}{4}s-\dfrac{1}{4}\\ \dfrac{1}{4}r-\dfrac{1}{4}&-\dfrac{1}{4}r-\dfrac{1}{4}&\dfrac{1}{4}r+\dfrac{1}{4}&\dfrac{1}{4}-\dfrac{1}{4}r \end{bmatrix} \right) $$