FEniCS: boundary conditions for electrostatic problems with dielectrics

Question

I carefully read all circa 70 pages of FEniCS tutorial and I still do not understand how to solve electrostatic problems when I have materials with different dielectric constant. The self contained system of equations in electrostatic $$ \mathrm{div}\boldsymbol{\mathcal{D}}=4\pi\rho,\\ \boldsymbol{\mathcal{D}}=\varepsilon\mathbf{E},\\ \mathbf{E}=-\nabla\varphi. $$ works only in regions where the functions are differentiable and smooth. At the interface between two regions with different materials these equations have to be augmented with the following boundary conditions: $$ \varphi_1=\varphi_2,\\ \varepsilon_1\frac{\partial\varphi_1}{\partial n}=\varepsilon_2\frac{\partial\varphi_2}{\partial n}. $$ How can I specify the second boundary condition? The example in the tutorial considers only a simple case of $$ \frac{\partial\varphi}{\partial n}=0 $$ between two regions.

Answer 1

Vatiational formulation of Poisson problem with $\varepsilon \in L^\infty$

Both interfacial conditions are incorporated in following variational problem. Given

• $\Omega$ Lipschitz domain
• $\partial\Omega = \Gamma_\mathrm{D} \cup \Gamma_\mathrm{N}$
• $V=\{v\in H^1(\Omega); v|_{\Gamma_\mathrm{D}}=0 \} $
• $\varepsilon\in L^\infty(\Omega)$
• $g \in L^2(\Gamma_\mathrm{N}) $

find $\varphi\in V$ such that

$$ \int_\Omega \varepsilon \nabla \varphi \cdot \nabla \psi \;\mathrm{d}x = \int_\Omega \rho \psi \;\mathrm{d}x + \int_{\Gamma_\mathrm{N}} g \psi \;\mathrm{d}S \quad \forall \psi\in V$$

Here $g$ represents Neumann condition on $\varepsilon\frac{\partial \varphi}{\partial \mathbf{n}} = g$ on $\Gamma_N$. Both your interfacial conditions are handled automatically by this formulation. You check it by dividing $\Omega$ into subdomains where $\varepsilon$ takes two different values and testing by something close to characteristic functions of these subdomains. Moreover $\varepsilon$ can be much more arbitrary - essential boundedness suffices. Note that I threw away $4\pi$ which corresponds to SI units. If you would like to derive this formulation from physicist point of view, just take Gateaux derivative of field energy minus charge energy in appropriate (accounting to Dirichlet condition) space.

For implementation in FEniCS just pick conforming FE space, i.e. arbitrary order Lagrange/CG elements. Problem is implemented in subdomains-poisson demo assuming discontinuity in $\epsilon$ matches mesh facets.

---

Discontinuity in $\varepsilon$ non-matching with mesh

If you need discontinuity non-matching mesh, you could try DOLFIN-PUM library which implemets XFEM/PUM method. It is not maintained now and Garth said that there are plans for better implementation directly into FEniCS. Library is able to integrate exactly this type of discontinuity in $\epsilon$. But there is additional problem that XFEM/PUM method introduces dicontinuity also to your solution and you need to enforce somehow continuity. It can be solved with Lagrange multiplier, penalty, etc. but I can't rememeber if I succedeed doing this.

You can also solve problem with non-matching mesh using standard formulation and increase quadrature degree or interpolate $\varepsilon$ to some higher-order space prior to computation. I also use $C^k$ approximation to Heaviside function with $k=0,1,2$ to define material discontinuities. You can implement this using conditional. None of these approaches can be integrate $\varepsilon$ exactly but it could enhance solution somehow.

Answer 2

Since the problem is treated in subdomains-poisson for the case, where the mesh is aligned with the interface, I assume in the following that in your case the jump can occur anywhere.

I think the first step here is not to think about implementation FEniCS, but rather to find an adequate variational form for this problem, where the interface conditions are properly incorporated. The problem looks like it can be treated with a variant of the following method:

Hansbo, A. and Hansbo, P. An unfitted finite element method, based on Nitsche’s method, for elliptic interface problems (Comp. Meth. in Appl. Mech. and Eng., Volume 191, Issues 47–48, 22 November 2002, Pages 5537–5552)

Disclaimer: There may be others and maybe also better alternatives for your specific case, but I am mostly familiar with Nitsche-type techniques.

Once you have an idea about how to discretize your problem, you can think about implementation in FEniCS. As far as I know, it is not possible to implement this method straightforwardly in vanilla FEniCS. dolfin-olm may be used to implement this method. I have not used it so far, I don't know if it works with the current version of FEniCS, and I don't know if they merged (parts of) it into dolfin in the meantime.

Again, note that in standard finite element methods the interface condition can be taken into account by aligning your mesh with the interface. If that is an option, you should go for that.