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

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?

• $u_h$ is just a polynomial on the reference cell mapped to an actual cell. Its second derivative requires you to compute derivatives of the mapping, but it's not "so expensive" that it becomes a problem. In fact, the polynomial degree of the shape functions doesn't matter very much (taking derivatives of polynomials is cheap), it's the degree of the mapping that matters. So my question is why this ends up being of concern to you? Mar 2 at 21:45
• I'm developing a p-adaptive finite element code. To decide which elements to enrich after solving the problem with the current mesh, the element-wise errors $\eta_K$ must be calculated. Here I don't know how to calculate $\nabla^2 u_h$ in a 2D quadrilateral master element $K$ when my interpolating shape functions are high order (my mapping shape functions are four-node linear).
– Masa
Mar 3 at 6:17
• I don't understand the problem. You're just trying to compute second derivatives of a polynomial. You know how to do that. Mar 3 at 14:52
• yes, the Laplacian $\nabla^2=\partial^2/\partial x^2+\partial^2/\partial y^2$ of $u_h=u_h(x,y)$ is simple, but here we have $u_h=u_h(r,s)$. Therefore, Laplacian in the integral must first be expressed in terms of $(r,s)$ instead of $(x,y)$, then it would become a simple differentiation.
– Masa
Mar 3 at 16:04
• It's just a contraction with the gradient of the mapping (a rank-2 tensor) and the second derivatives of the mapping (a rank-3 tensor). But the mapping is generally low-order polynomial, so its derivatives are easy to compute. Mar 3 at 20:17

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.

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)$$