Which finite element to use?

mathematics
2

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.

2 Likes