Problem with plotting subdomains after importing mesh

Question

Hi
I am trying to plot 5 subdomains defined for a mesh imported from GMSH. The imported mesh includes 5 rectangles (meshed separately)corresponding to the 5 subdomains. Here is the GMSH file:

Point(1) = {0, -2e-6, 0, 1.0};
Point(2) = {0, 0, 0, 1.0};
Point(3) = {0, 10e-6, 0, 1.0};
Point(4) = {0, 190e-6, 0, 1.0};
Point(5) = {0, 200e-6, 0, 1.0};
Point(6) = {0, 202e-6, 0, 1.0};
Point(7) = {0.03, 202e-6, 0, 1.0};
Point(8) = {0.03, 200e-6, 0, 1.0};
Point(9) = {0.03, 190e-6, 0, 1.0};
Point(10) = {0.03, 10e-6, 0, 1.0};
Point(11) = {0.03, 0, 0, 1.0};
Point(12) = {0.03, -2e-6, 0, 1.0};
Line(1) = {1, 12};
Line(2) = {12, 11};
Line(3) = {11, 2};
Line(4) = {2, 1};
Line(5) = {2, 3};
Line(6) = {10, 3};
Line(7) = {11, 10};
Line(8) = {10, 9};
Line(9) = {9, 4};
Line(10) = {4, 3};
Line(11) = {4, 5};
Line(12) = {5, 8};
Line(13) = {8, 9};
Line(14) = {8, 7};
Line(15) = {7, 6};
Line(16) = {6, 5};

Line Loop(17) = {1, 2, 3, 4};
Plane Surface(18) = {17};

Transfinite Surface{18};
Transfinite Line{2,4}=2;
Transfinite Line{1,3}=2000;

Line Loop(19) = {3, 5, -6, -7};
Plane Surface(20) = {19};

Transfinite Surface{20};
Transfinite Line{5,7}=4;
Transfinite Line{3,6}=2000;

Line Loop(21) = {10, -6, 8, 9};
Plane Surface(22) = {21};

Transfinite Surface{22};
Transfinite Line{8,10}=20;
Transfinite Line{6,9}=2000;

Line Loop(23) = {11, 12, 13, 9};
Plane Surface(24) = {23};

Transfinite Surface{24};
Transfinite Line{11,13}=4;
Transfinite Line{9,12}=2000;

Line Loop(25) = {16, 12, 14, 15};
Plane Surface(26) = {25};

Transfinite Surface{26};
Transfinite Line{14,16}=2;
Transfinite Line{12,15}=2000;

When I defined and plotted the subdomains like:

from dolfin import *
mesh = Mesh("mesh.xml")

class electrode_1(SubDomain):
  def inside(self, x, on_boundary):
      return True if x[0] >= 0.0 and x[0] <= 0.03 and x[1] >= (-2.)*(10**(-6)) and x[1] <= 0.0 else False

class diffusion_layer_1(SubDomain):
  def inside(self, x, on_boundary):
      return True if x[0] >= 0.0 and x[0] <= 0.03 and x[1] >=  0.0 and x[1] <= 10*(10**(-6)) else False

class polymer(SubDomain):
  def inside(self, x, on_boundary):
      return True if x[0] >= 0.0 and x[0] <= 0.03 and x[1] >=  10.*(10**-6) and x[1] <= 190.*(10**-6) else False

class diffusion_layer_2(SubDomain):
  def inside(self, x, on_boundary):
      return True if x[0] >= 0.0 and x[0] <= 0.03 and x[1] >=  190.*(10**-6) and x[1] <= 200.*(10**-6) else False

class electrode_2(SubDomain):
  def inside(self, x, on_boundary):
      return True if x[0] >= 0.0 and x[0] <= 0.03 and x[1] >= 200.*(10**-6) and x[1] <= 202.*(10**-6) else False

#marking subdomains
domains = CellFunction("size_t", mesh)
domains.set_all(0)


electrode_1 = electrode_1()
electrode_1.mark(domains, 1)
diffusion_layer_1 = diffusion_layer_1()
diffusion_layer_1.mark(domains, 2)
polymer = polymer()
polymer.mark(domains, 3)
diffusion_layer_2 = diffusion_layer_2()
diffusion_layer_2.mark(domains, 4)
electrode_2 = electrode_2()
electrode_2.mark(domains, 5)

plot(domains, interactive=True)    

I expect the defined subdomains to be plotted perfectly like:

But instead this is what I am getting:

Does anybody know what is wrong with that? How can I fix it?

Answer 1

Leo,

The problem might be that you're not using DOLFIN_EPS in the conditional statements which use floats to define the subdomains. For example when you're checking a subdomain from 0 to 1.3 you should use (x[0] > -DOLFIN_EPS) and (x[0] < 1.3 + DOLFIN_EPS) to make sure the boundaries at 0 and 1.3 are included

Comments

I modified the definition of my subdomains according to what you suggested:

    class electrode_1(SubDomain):
  def inside(self, x, on_boundary):
      return True if x[0] >= 0.0 and x[0] <= 0.03 and (x[1] >= (-2.)*(10**(-6))-DOLFIN_EPS)\ and (x[1] <= 0.0+DOLFIN_EPS) else False

class diffusion_layer_1(SubDomain):
  def inside(self, x, on_boundary):
      return True if x[0] >= 0.0 and x[0] <= 0.03 and (x[1] >=  0.0-DOLFIN_EPS)\
                     and (x[1] <= 10*(10**(-6))+DOLFIN_EPS) else False

class polymer(SubDomain):
  def inside(self, x, on_boundary):
      return True if x[0] >= 0.0 and x[0] <= 0.03 and (x[1] >=  10.*(10**-6)-DOLFIN_EPS)\
                     and (x[1] <= 190.*(10**-6)+DOLFIN_EPS) else False

class diffusion_layer_2(SubDomain):
  def inside(self, x, on_boundary):
      return True if x[0] >= 0.0 and x[0] <= 0.03 and (x[1] >=  190.*(10**-6)-DOLFIN_EPS)\
                     and (x[1] <= 200.*(10**-6)+DOLFIN_EPS) else False

class electrode_2(SubDomain):
  def inside(self, x, on_boundary):
      return True if x[0] >= 0.0 and x[0] <= 0.03 and (x[1] >= 200.*(10**-6)- DOLFIN_EPS)\
                     and (x[1] <= 202.*(10**-6)+DOLFIN_EPS) else False

It solved the problem nicely.
Thank you so much!

---

You're welcome. Its one of the first things I remember learning while doing programming in undergrad; if the number is not an integer, you can't use statements like x == 1.3 or the example I gave you, because the values might be slightly off due to machine precision.

So that's why fenics has a small number you can use (DOLFIN_EPS), and a function to do x == 1.3 for floats which is near(x, 1.3). What this function actually does is:

 abs(x - 1.3) < DOLFIN_EPS 

Good luck on your project.