integration of a laplacian

Question

Hi, I solved for a Poisson equation with finite elements, using piecewise linear basis functions on 2d triangles.
Now, I want to evaluate the following expressions:
$$ \int_\Omega \Delta u ~dx$$ and $$ \int_\Omega (\Delta u)^2 ~dx$$ I want to evaluate these expressions using my approximated solution $u$ which has been computed on piecewise linear basis functions.

For the first one, I thought of using the identity $\int_\Omega div \nabla u = \int_{\partial\Omega} \nabla u . \vec{n}~ds$ and summing this expression over each triangle.However, as expected, the result is strictly 0 since the basis functions are linear.
I also tried to use a kind of jump formula (like $f'(x)=\tilde{f}'(x) + f^+-f^-$ where $\tilde{f}'$ is the derivative of the smooth part of f) but I'm stuck on how to do that for each triangle in 2D (the outer normal is likely to cancel out when computing the same formula for two adjacent triangles sharing an edge) - and I'm wondering if it is supposed to work.

For the second one, I just have no clue.

Am I forced to use higher order elements ? Any idea ? Thanks!

Answer 1

Since your approximate solution is piecewise linear, it is $H^1$ but not $H^2$. Therefore your calculation is impossible. You can do to things to overcome the difficulty:

• Either use higher-order elements,
• or post-process your approximate solution $u$. This means constructing a smoother $\bar u$ using some convolution by an appropriate $\phi(x/h)$, where $h$ is the typical mesh size. In one space dimension, this amounts to using splines. Then $\bar u\in H^2$ and your calculation is meaningful.

Answer 2

Denis has this exactly right, if your goal is really to calculate these integrals. However, if your real goal (as you say) is to calculate the residual, then this isn't what you want to do at all.

In a weak sense, the Laplacian is a map $ \Delta \colon H^1 (\Omega) \to H^{-1} (\Omega) $, so the PDE $ \Delta u = f $ makes sense when $ u \in H^1 (\Omega)$ and $ f \in H^{-1} (\Omega)$. Denoting the approximate FEM solution by $u_h$, the residual is $ f - \Delta u_h \in H^{-1} (\Omega) $, so it really makes sense to measure the residual in the $ H^{-1} (\Omega) $ norm, not the $ L^1 (\Omega)$ or $L^2 (\Omega)$ norm. That is, $$ \lVert f - \Delta u_h \rVert _{H^{-1}(\Omega)} = \sup _{ \lVert v \rVert _{H^1 (\Omega)}= 1} \langle f - \Delta u_h , v \rangle _{H^{-1} (\Omega) \times H^1 (\Omega) }.$$

On the other hand, maybe you don't really want to measure the residual itself; you want to estimate the a posteriori error $ e _h = u - u_h $. In this case, $ e _h \in H^1 (\Omega)$ solves the residual equation $$ \Delta e _h = \Delta (u - u_h) = f - \Delta u_h .$$ You can measure $ e _h $ a number of ways, e.g., using the energy norm. Typically, of course, you can't actually solve for $e_h$ (since that would mean solving the original PDE exactly!), but you can estimate it by using a more accurate finite-element method for the residual equation (e.g., finer mesh and/or higher-order elements) than you used for $ u_h $.

To learn more about these sorts of things, you should look up residual-based a posteriori error estimation (Google returns lots of hits for this phrase).