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.


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.


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.

  • $\begingroup$ Let's say I can align my mesh with the interface. How can I specify the boundary condition $\varepsilon_1\frac{\partial\varphi_1}{\partial n}=\varepsilon_2\frac{\partial\varphi_2}{\partial n}$ on that interface then? $\endgroup$
    – facetus
    May 26 '13 at 11:12
  • $\begingroup$ As far as I can see, you don't need to care about the interface condition then. From a variational point of view, only one of the conditions has to be satisfied. For instance, when you look at the Hansbo-Paper, you can see that for $g=0$ (which is your case), all interface terms vanish if your trial and test functions are continuous across the interface. $\endgroup$ May 26 '13 at 11:18
  • $\begingroup$ Ah, right! FEM solves an integral form of the Laplace equation. In integral form $\varphi$ doesn't have to be smooth. Thanks! $\endgroup$
    – facetus
    May 26 '13 at 11:48
  • $\begingroup$ I read the tutorial. The tutorial is chapter 1 of the book. I didn't say I read anything that includes Chapter 9. I don't appreciate that your comment makes it look like as if I lied or something. Thanks for pointing me to the chapter that I should read anyway. $\endgroup$
    – facetus
    May 26 '13 at 19:01
  • $\begingroup$ I know about Section 5 of the tutorial. I carefully read it, as I absolutely clearly mentioned it in my initial post. That section does not contain information about the kind of boundary conditions at interfaces I was interested in. I mentioned this in my initial questions as well. I hope this will resolve your doubts about the fact that I can attentively read English texts. $\endgroup$
    – facetus
    May 28 '13 at 1:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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