Jump condition for elliptic equation in standard finite element method

Question

Consider the elliptic PDE

$ -\mathrm{div}(k(x)\nabla u) = f(x) \in \Omega$ and

$u = 0$ on $\partial \Omega$

If $k(x)$ is piecewise constant across an interface $\Gamma$ in the domain, then we have the jump condition

$k(x)\nabla u\cdot \mathbf{n}|_{\Gamma+} = k(x)\nabla u \cdot \mathbf{n}|_{\Gamma-}$ where $\mathbf{n}$ is the unit outward normal to the interface.

When standard Lagrange finite elements are used for this problem, I understand that $u$ is continuous across the interface provided the mesh is fitted to the interface.

How does the standard formulation enforce the jump condition ?

Answer 1

Preliminary remarks:

1. The natural jump condition is on the normal derivative, not on the full gradient.
2. The fact that the solution is continuous across the interface is part of the formulation of the problem. This is not only due to the finite element formulation.

Formulation

The standard formulation simply enforces it with... no extra term.

Start by writing the formulation on the two subdomains (say $\Omega_1$ and $\Omega_2$, separated by the interface $\Gamma$):

$$ \int_{\Omega_1} k_1 \nabla u_1 \cdot \nabla v_1 - \int_{\Gamma} k_1 (\nabla u_1 \cdot n_1) v_1 = \int_{\Omega_1} f v_1 $$

and

$$ \int_{\Omega_2} k_2 \nabla u_2 \cdot \nabla v_2 - \int_{\Gamma} k_2 (\nabla u_2 \cdot n_2) v_2 = \int_{\Omega_1} f v_2 $$

where $u_i = u|_{\Omega_i}$ and similar for $v$. Now sum the two equations and you get:

$$ \int_{\Omega} k \nabla u \cdot \nabla v - \int_{\Gamma} k_1 (\nabla u_1 \cdot n_1) v_1 - \int_{\Gamma} k_2 (\nabla u_2 \cdot n_2) v_2 = \int_{\Omega} fv$$

Now observe that the two normals to $\Gamma$ are in the opposite direction $n_1 = -n_2$ and that on $\Gamma$, $v_1 = v_2$ (due to the continuity condition) and you finally get

$$ \int_{\Omega} k \nabla u \cdot \nabla v = \int_{\Omega} fv$$

So, no extra term, like for homogeneous Neumann conditions. Remark that if the condition was inhomogeneous, there would be an extra term!

Approximation

Now for the approximation, if you do not pay attention and simply apply the standard finite element method, your approximation won't be good close to the interface, since the elements cannot capture this discontinuity. There are in general two ways to overcome this issue:

• Use a mesh that is fitted with the interface. Then you can use the formulation above and everything should work well. However, this can be tricky if the interface is complicated or if it moves (in a time dependant problem).
• Use enriched finite elements, aka XFEM, that are built with the purpose of capturing this kind of discontinuities. But the formulation is then a bit different and the implementation more involved.