Which finite element to use?

Hi everyone!

I need your help: Which finite element do I use for a Laplace problem of the form

\nabla \sigma \nabla \phi = 0,

where \sigma represents conductivity and \phi a potential, if

  1. \sigma is a step function
  2. I want to calculate the flux \vec{j} = \sigma \nabla \phi at the interface of the two regions with the different sigmas later.

The boundary conditions are dirichlet boundary conditions for two sides.

The setup:

Which finite element / function space do you recommend to use? I’ve thought: the potential \phi has to be differentiable at the interface in the normal direction and continuous.

Any help is very much appreciated! Thanks, Baltasar

If your discretization is based on the standard weak form of: Find \phi\in\mathcal{S} such that

\int_\Omega(\sigma\nabla\phi)\cdot\nabla w\,d\Omega = 0\quad\forall w\in\mathcal{V}\text{ ,}

where \mathcal{S} is a trial space that satisfies Dirichlet boundary conditions and \mathcal{V} is a test space that is zero on the Dirichlet part of the boundary, then it should be fine to use an H^1-conforming discrete space, like

V = FunctionSpace(mesh,"CG",k)

i.e., piecewise polynomial functions of degree k that satisfy C^0 continuity between elements.

I’ve thought: the potential ϕ has to be differentiable at the interface in the normal direction and continuous.

You do want to use a continuous function space for \phi, to ensure that it is in H^1, but you actually expect a jump in \nabla\phi\cdot\mathbf{n} at the interface (where \mathbf{n} is the unit normal to the interface), because, if the normal flux \mathbf{j}\cdot\mathbf{n} is continuous (to satisfy the conservation law) but \sigma is different on the two sides, then \nabla\phi\cdot\mathbf{n} must be discontinuous to compensate.

Of course, in the approximate solution using this method, normal flux continuity is only enforced weakly, so the fluxes on either side of the interface differ slightly on a finite mesh. If you just want to integrate normal flux over the material interface, you could average the value between both sides (using the UFL avg function in a dS integral over interior facets). However, it’s worth pointing out that controlling the orientation of normal vectors on interior facets is tricky; a good example can be found here.

If you want flux as a vector field over the full domain that exactly maintains normal continuity on boundaries between elements, you can project the flux from the H^1-conforming formulation onto an H(\text{div})-conforming space like Raviart–Thomas, e.g.,

j = project(sigma*grad(phi),FunctionSpace(mesh,"RT",k))

Alternatively, you could also use a mixed formulation that enforces flux continuity strongly, by solving for an H(\text{div})-conforming flux as one of the unknowns. The mixed Poisson formulation is demonstrated in this FEniCS demo.


Hi David,
Thank you a lot for the fantastic answer. I guess you’ve answered all my questions and even more :bouquet:.

Have a nice weekend (it’s Friday!)