Solution voltage strangely set to zero in Poisson equation with dielectrics.

Question

I'm trying to set up a demo that has a dielectric inside of some capacitor plates to show how it alters the electric field when a dielectric body is present. However, it seems like the dielectric body is set to zero voltage instead.

There are two regions, basically a rectangle with an irregular shape in the center. The solution voltage inside of the irregular shape that is in the center always shows up as zero voltage. It doesn't seem to change if I set the voltage boundary conditions differently - or if I set the dielectric constant differently.

Is it that I am not doing the boundary conditions correctly?

  #-------------------------------#
  # Construct Mesh, FunctionSpace #
  #-------------------------------#
  mesh = Mesh("mesh.xml.gz")
  V = FunctionSpace(mesh, "Lagrange", 1)

  #----------------------#
  # Dielectric Materials #
  #----------------------#
  D         = FunctionSpace(mesh, "DG", 0)
  epsilon   = Function(D)
  materials = MeshFunction("size_t", mesh, "material.xml.gz")
  helper    = numpy.asarray(materials.array(), dtype=numpy.int32) - 1
  dm        = D.dofmap()
  for cell in cells(mesh):
    helper[dm.cell_dofs(cell.index())] = materials[cell] - 1
  epsilon.vector()[:] = numpy.choose(helper, (3, 1))
  # if I plot epsilon here, I get what I expect to get ...

  #-------------------------------#
  # Dirichlet boundary conditions #
  #-------------------------------#
  dirichlet_bcs = []
  dirichlet_bcs.append(DirichletBC(V, Constant(1), TopBoundary))
  dirichlet_bcs.append(DirichletBC(V, Constant(-1), BottomBoundary))

  #------------#
  # Define PDE #
  #------------#
  u = TrialFunction(V)
  v = TestFunction(V)
  a = epsilon*inner(grad(u), grad(v))*dx
  L = Constant(0.0)*v*dx
  A, b = assemble_system(a, L, dirichlet_bcs)
  u = Function(V)
  U = u.vector()

  # Compute solution
  solve(A, U, b)

  # Save solution in VTK format
  f = File("problem.pvd")
  f << u

  # Plot solution
  plot(u, interactive=True)

Comments

Hi it is difficult to comment without the mesh but do you mean that you have 0 potential inside the whole dielectric? Also how does the intensity of electric field look if you plot it? Finally could you try with this dielectric that can be easily reproduced by others?

mesh = UnitSquareMesh(20, 20)                                                    
materials = CellFunction("size_t", mesh, 0)                                      
for cell in cells(mesh):
  mp = cell.midpoint()                                                   
  x = mp.x()                                                        
  y = mp.y()                                                        
  if (0.25 < x < 0.75) and (0.25 < y < 0.75):                                                          
    materials[cell] = 1   

With this one and your material settings I observe no zero potential inside the small square.
The effect of dielectric on the potential is almost invisible in the plot of u but check out the
intensity vector. If you want to see bigger effects on u try to increase the ratio between permitivities or a dielectric as long as the capacitor itself.

---

This post keeps being eaten for some reason … so I will cut it in half and try again. It really isn't that large. I made a mesh file and a mesh function file, that when fed into the code above generate the results described. Admittedly, the answer may be related to this question.

Part 1: the mesh function file:

<?xml version="1.0"?>
<dolfin xmlns:dolfin="http://fenicsproject.org">
  <mesh_function>
    <mesh_value_collection name="cells" type="uint" dim="2" size="16">
      <value cell_index="0" local_entity="0" value="1" />
      <value cell_index="1" local_entity="0" value="1" />
      <value cell_index="2" local_entity="0" value="2" />
      <value cell_index="3" local_entity="0" value="2" />
      <value cell_index="4" local_entity="0" value="2" />
      <value cell_index="5" local_entity="0" value="2" />
      <value cell_index="6" local_entity="0" value="2" />
      <value cell_index="7" local_entity="0" value="2" />
      <value cell_index="8" local_entity="0" value="2" />
      <value cell_index="9" local_entity="0" value="2" />
      <value cell_index="10" local_entity="0" value="2" />
      <value cell_index="11" local_entity="0" value="2" />
      <value cell_index="12" local_entity="0" value="2" />
      <value cell_index="13" local_entity="0" value="2" />
      <value cell_index="14" local_entity="0" value="2" />
      <value cell_index="15" local_entity="0" value="2" />
    </mesh_value_collection>
  </mesh_function>
</dolfin>

---

Part 2, the mesh file:

<?xml version="1.0" encoding="UTF-8"?>
<dolfin xmlns:dolfin="http://fenicsproject.org">
  <mesh celltype="triangle" dim="2">
    <cells size="16">
      <triangle index="0" v0="0" v1="2" v2="1" />
      <triangle index="1" v0="0" v1="3" v2="2" />
      <triangle index="2" v0="4" v1="6" v2="5" />
      <triangle index="3" v0="7" v1="9" v2="8" />
      <triangle index="4" v0="10" v1="12" v2="11" />
      <triangle index="5" v0="13" v1="14" v2="11" />
      <triangle index="6" v0="14" v1="13" v2="15" />
      <triangle index="7" v0="5" v1="14" v2="15" />
      <triangle index="8" v0="10" v1="9" v2="7" />
      <triangle index="9" v0="6" v1="9" v2="5" />
      <triangle index="10" v0="4" v1="5" v2="15" />
      <triangle index="11" v0="16" v1="10" v2="11" />
      <triangle index="12" v0="12" v1="10" v2="7" />
      <triangle index="13" v0="8" v1="9" v2="6" />
      <triangle index="14" v0="14" v1="16" v2="11" />
      <triangle index="15" v0="16" v1="14" v2="5" />
    </cells>
    <vertices size="17">
      <vertex index="0" x="1.94062500e+01" y="8.12500000e-01" />
      <vertex index="1" x="2.01250000e+01" y="8.12500000e-01" />
      <vertex index="2" x="2.01250000e+01" y="0.00000000e+00" />
      <vertex index="3" x="1.83750000e+01" y="0.00000000e+00" />
      <vertex index="4" x="1.73750000e+01" y="-1.00000000e+00" />
      <vertex index="5" x="1.83750000e+01" y="0.00000000e+00" />
      <vertex index="6" x="1.92500000e+01" y="-1.00000000e+00" />
      <vertex index="7" x="2.11250000e+01" y="4.06250000e-01" />
      <vertex index="8" x="2.11250000e+01" y="-1.00000000e+00" />
      <vertex index="9" x="2.01250000e+01" y="0.00000000e+00" />
      <vertex index="10" x="2.01250000e+01" y="8.12500000e-01" />
      <vertex index="11" x="1.92500000e+01" y="1.81250000e+00" />
      <vertex index="12" x="2.11250000e+01" y="1.81250000e+00" />
      <vertex index="13" x="1.73750000e+01" y="1.81250000e+00" />
      <vertex index="14" x="1.83562500e+01" y="9.68750000e-01" />
      <vertex index="15" x="1.73750000e+01" y="4.06250000e-01" />
      <vertex index="16" x="1.94062500e+01" y="8.12500000e-01" />
    </vertices>
  </mesh>
</dolfin>

Answer 1

Hi, you are getting 0 voltage because the mesh is wrong. I don't think that the problem is due to ordering that you mention in your linked question but rather the number of vertices. In the above mesh you declare 17 vertices but note that vertex 3 is the same as vertex 5 and similarly for ((2, 9), (1, 10), (0, 16)) so there are only 13 unique vertices. If you fix the mesh, you will get correct solution. Try with

<?xml version="1.0" encoding="UTF-8"?>
<dolfin xmlns:dolfin="http://fenicsproject.org">
  <mesh celltype="triangle" dim="2">
    <cells size="16">
      <triangle index="0" v0="0" v1="2" v2="1" />
      <triangle index="1" v0="0" v1="3" v2="2" />
      <triangle index="2" v0="4" v1="5" v2="3" />
      <triangle index="3" v0="6" v1="2" v2="7" />
      <triangle index="4" v0="1" v1="9" v2="8" />
      <triangle index="5" v0="10" v1="11" v2="8" />
      <triangle index="6" v0="11" v1="10" v2="12" />
      <triangle index="7" v0="3" v1="11" v2="12" />
      <triangle index="8" v0="1" v1="2" v2="6" />
      <triangle index="9" v0="5" v1="2" v2="3" />
      <triangle index="10" v0="4" v1="3" v2="12" />
      <triangle index="11" v0="0" v1="1" v2="8" />
      <triangle index="12" v0="9" v1="1" v2="6" />
      <triangle index="13" v0="7" v1="2" v2="5" />
      <triangle index="14" v0="11" v1="0" v2="8" />
      <triangle index="15" v0="0" v1="11" v2="3" />
    </cells>
    <vertices size="13">
      <vertex index="0" x="1.94062500e+01" y="8.12500000e-01" />
      <vertex index="1" x="2.01250000e+01" y="8.12500000e-01" />
      <vertex index="2" x="2.01250000e+01" y="0.00000000e+00" />
      <vertex index="3" x="1.83750000e+01" y="0.00000000e+00" />
      <vertex index="4" x="1.73750000e+01" y="-1.00000000e+00" />
      <vertex index="5" x="1.92500000e+01" y="-1.00000000e+00" />
      <vertex index="6" x="2.11250000e+01" y="4.06250000e-01" />
      <vertex index="7" x="2.11250000e+01" y="-1.00000000e+00" />
      <vertex index="8" x="1.92500000e+01" y="1.81250000e+00" />
      <vertex index="9" x="2.11250000e+01" y="1.81250000e+00" />
      <vertex index="10" x="1.73750000e+01" y="1.81250000e+00" />
      <vertex index="11" x="1.83562500e+01" y="9.68750000e-01" />
      <vertex index="12" x="1.73750000e+01" y="4.06250000e-01" />
    </vertices>
  </mesh>
</dolfin>

These duplicit vertices all sat in cells that defined the dielectric so perhaps when you stored the mesh they were written to file twice, as part of the whole mesh and as part of the dielectric.

Comments

I think you're right about how the duplicate vertices got there. I am grateful to you for finding that out!