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 ?


1 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.


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


$$ \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!


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.
  • $\begingroup$ Thank you for the answer. I understand now. I also edited the jump condition where I had missed the normal. $\endgroup$
    – me10240
    Commented Aug 29, 2014 at 12:57
  • $\begingroup$ A further question on the approximation. As I understand it, if the mesh is not fitted to the interface, then the element containing the interface will introduce a discretization error when using quadrature. Is that correct ? $\endgroup$
    – me10240
    Commented Aug 29, 2014 at 13:08
  • $\begingroup$ It is not because of the quadrature, but because of the inability of the finite element space to approximate correctly the exact solution. Even if you could interpolate the exact solution on the finite element space, the error would be quite large because of the "kink" that is located at the interface. $\endgroup$
    – Dr_Sam
    Commented Aug 29, 2014 at 13:37
  • 1
    $\begingroup$ Just think about the interpolation, that already contains the whole problem. Make a draw of a possible solution in 1D: from $a$ to $\Gamma$ a smooth curve, then create a kink in $\Gamma$ (keep the curve continuous but the change a lot the derivative) and go on smoothly from $\Gamma$ to $b$. Now interpolate the curve once with a node in $\Gamma$ (curve A) and once without node in $\Gamma$ (curve B). Do you see the difference in the approximation? $\endgroup$
    – Dr_Sam
    Commented Aug 29, 2014 at 14:01
  • 1
    $\begingroup$ This is because curve A is just two interpolations of two smooth pieces of curves, whereas curve B is the interpolation of a continuous but not $C^1$ curve. So the convergence of the curve B will be slower than with the curve A. And that is about the same thing with the finite element approximation. $\endgroup$
    – Dr_Sam
    Commented Aug 29, 2014 at 14:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.