Strange results from simple cylinder mesh and Gauss's Law

General
1

I am trying to run a simple simulation of Gauss’s law as it applies to a cylinder with the top and bottom ends held at a constant charge and the sides of the cylinder kept at 0 charge. I am using gmsh to create the simple mesh which I then convert into a DolfinX mesh using model_to_mesh. I then visualise the solution using Paraview. Things appear to work till this point even though the electric field looks strange to me. However, when I try to add some objects to the simulation that have different permittivitys it looks like the simulation is completely ignoring these objects. I’m pretty sure the boundary conditions are set correctly as the electric potential looks right in Paraview but the electric field lines do not make sense and seem to ignore the objects of very high and very low permittivity. This is my example code.

from mpi4py import MPI
from dolfinx.io import gmshio
from dolfinx.io.utils import XDMFFile
import numpy as np
import gmsh

from ufl import TestFunction, TrialFunction, dot, dx, grad
from dolfinx.fem import (Constant, Function, FunctionSpace, Expression, VectorFunctionSpace,
                         dirichletbc, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from petsc4py.PETSc import ScalarType
import matplotlib.pyplot as plt


# specified in mm
Tube_hight = 200
Tube_radius = 10
tube_perm = 1  # 4 * np.pi*1e-7 # Vacuum
obj1_perm = 100 # 1.26e-4 # Copper
obj2_perm = 0.001  # 1.26e-10 # plastic

gmsh.initialize()
# gmsh.model.add("cylinder_mesh")
gmsh.option.setNumber("Mesh.MeshSizeMin", 2)
gmsh.option.setNumber("Mesh.MeshSizeMax", 2)

# Create model
main_model = gmsh.model()
main_model.add('main_mesh')
main_model.setCurrent('main_mesh')

# Add the cylinder to the model
# disks for the bottom and top of the cylinder
bottom_disk = main_model.occ.addDisk(0, 0, 0, Tube_radius, Tube_radius, tag=1)
gmsh.model.occ.synchronize()
main_model.addPhysicalGroup(2, [bottom_disk], tag = 100, name="bottom_disk")
# extrude to make a cylinder
subdivision = [Tube_hight]
extrusion = gmsh.model.occ.extrude([(2, bottom_disk)], 0, 0, Tube_hight, subdivision)
# Generate the 3D geometry
main_model.occ.synchronize()
# We need to add physical surfaces for Fenicx
main_model.addPhysicalGroup(3, [extrusion[1][1]], tag=10, name="volume")
main_model.addPhysicalGroup(2, [extrusion[2][1]], tag = 101, name="side_surface")
main_model.addPhysicalGroup(2, [extrusion[0][1]], tag = 102, name="top_disk")

# Add an object with large permitivity
obj1 = main_model.occ.addSphere(0, 0, 150, 5, tag=103)
main_model.occ.synchronize()
main_model.addPhysicalGroup(3, [obj1], name="obj1")

# Add an object with small permitivity
obj2 = main_model.occ.addSphere(0, 0, 50, 5, tag=104)
main_model.occ.synchronize()
main_model.addPhysicalGroup(3, [obj2], name="obj2")

main_model.mesh.generate(3)
# main_model.mesh.refine()
gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
domain, ct, ft = gmshio.model_to_mesh(main_model, mesh_comm, gmsh_model_rank, gdim=3)
gmsh.finalize()
# Tags: 
# 10 = tube interior
# 101 = tube sides
# 102 = top disk
# 100 = bottom disk
# 103 = obj1 (large perm)
# 104 = obj2 (small perm)

# Set up the functions
V = FunctionSpace(domain, ("CG", 1))
u = TrialFunction(V)
v = TestFunction(V)

# Define the material parameters
Q = FunctionSpace(domain, ("DG", 0))
material_tags = np.unique(ct.values)
eps = Function(Q)
# As we only set some values in eps, initialize all as vacuum
eps.x.array[:] = tube_perm
for tag in material_tags:
    cells = ct.find(tag)
    if tag == 103:
        eps_ = obj1_perm
    elif tag == 104:
        eps_ = obj2_perm
    else:
        eps_ = tube_perm
    eps.x.array[cells] = np.full_like(cells, eps_, dtype=ScalarType)


# Set the boundary conditions
tdim = domain.topology.dim
fdim = tdim - 1
# The top boundary is set to 10V
top_disk_facets = ft.indices[ft.values == 102]
top_disk_dofs = locate_dofs_topological(V, fdim, top_disk_facets)
u_T = Constant(domain, ScalarType(10))
top_bc = dirichletbc(u_T, top_disk_dofs, V)
# The bottom boundary is set to -10V
bottom_disk_facets = ft.indices[ft.values == 100]
bottom_disk_dofs = locate_dofs_topological(V, fdim, bottom_disk_facets)
u_B = Constant(domain, ScalarType(-10))
bottom_bc = dirichletbc(u_B, bottom_disk_dofs, V)
# The side boundary is set to 0V
side_disk_facets = ft.indices[ft.values == 101]
side_disk_dofs = locate_dofs_topological(V, fdim, side_disk_facets)
u_S = Constant(domain, ScalarType(0))
sides_bc = dirichletbc(u_S, side_disk_dofs, V)
bcs = [top_bc, bottom_bc, sides_bc]

# Construct the PDE and solve
a = dot(grad(u), grad(v)) * eps * dx
L = Constant(domain, ScalarType(1.e-16)) * v * dx 
# Solve the linear problem
problem = LinearProblem(a, L, bcs=bcs)
uh = problem.solve()

# calculate the electric field
E_field = -grad(uh)
# project to mesh
V_E_field = VectorFunctionSpace(domain, ("CG", 2))
E_field_expr = Expression(E_field, V_E_field.element.interpolation_points())
E_field_projected = Function(V_E_field)
E_field_projected.interpolate(E_field_expr)

# Export solution
with XDMFFile(domain.comm, "output.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    uh.name = "Electric potential"
    xdmf.write_function(uh)
    E_field_projected.name = "Electric field"
    xdmf.write_function(E_field_projected)

This is the electric potential which looks correct.

potentail
potentail1299×907 106 KB

This is the field electric field. Note how the field abruptly changes in the center of the cylinder and how the objects inside the cylinder (highlighted) are ignored.
field.PNG
field.PNG1719×1157 89.2 KB

2

You should never do synchronize after adding physical tags, as gmsh removes them.

Also you are not using boolean operations with gmsh, such as union, cut and fragment, usually causing the geometry to be disconnected.

3

Thanks for the advice. When should I call synchronize? Also, do mean that I should be using boolean operations?

4

Before usage of any boolean operation, and before marking physical entities.

Boolean difference is covered here: Using the GMSH Python API to generate complex meshes | Jørgen S. Dokken
And other operations and references here: Using subdomains in diffusion equation - #6 by dokken

5

Thanks for the help. I will look into using boolean operations based on your tutorial. I am still usnsure of when to use synchronize though. If I only synchronize once before marking the physical groups and never after adding any physical tags I get the following errors

Warning : Unknown entity of dimension 3 and tag 1 in physical group 10
Warning : Unknown entity of dimension 2 and tag 2 in physical group 101
Warning : Unknown entity of dimension 2 and tag 3 in physical group 102
Warning : Unknown entity of dimension 3 and tag 103 in physical group 103
Warning : Unknown entity of dimension 3 and tag 104 in physical group 104

I think this is because I cam not calling synchronize before adding the physical groups. But I can’t call it after adding the groups as then, as you mentioned, gmsh removes them? It feels like I’m missing something fundamental.

6

I modified the code so synchronize is only done once before the physical groups are added. The resulting simulation is identical to my last result so I suspect the problem lies with the mesh and the way the two spheres intersect with the main cylinder. I tried to use the cut method to extract the spheres from the cylinder but this completely removes the cells of those spheres from the simulation domain. I need those cells to be part of the domain but with different properties compared to the rest of the cylinder. This is the code I use to to set the material properties. The eps function then gets used directly in the PDE as seen in my original post.

# Define the material parameters
Q = FunctionSpace(domain, ("DG", 0))
material_tags = np.unique(ct.values)
eps = Function(Q)
# As we only set some values in eps, initialize all as vacuum
eps.x.array[:] = tube_perm
for tag in material_tags:
    cells = ct.find(tag)
    if tag == obj1_tag:
        eps_ = obj1_perm
    elif tag == obj2_tag:
        eps_ = obj2_perm
    else:
        eps_ = tube_perm
    eps.x.array[cells] = np.full_like(cells, eps_, dtype=ScalarType)

My code is based on the Fenicx tutorial Defining subdomains for different materials

7

Thanks for your help and advice @dokken. I managed to fix the problem with some experimentation and using the occ.fragment method.

# extrude to make a cylinder
subdivision = [Tube_hight]
extrusion = gmsh.model.occ.extrude([(2, bottom_disk)], 0, 0, Tube_hight, subdivision)
# Add an object with large permitivity
obj1 = main_model.occ.addSphere(0, 0, 150, 5)
# Add an object with small permitivity
obj2 = main_model.occ.addSphere(0, 0, 50, 5)
# Need to filter this volumne to exclude the extra objects 
# so cut out the spheres from the total volumne 
tube_volume = extrusion[1][1]
# empty_space = gmsh.model.occ.cut([(3, tube_volume)], [(3, obj1), (3, obj2)])
cell_fragments = gmsh.model.occ.fragment([(3, tube_volume)], [(3, obj1), (3, obj2)])

# Tag the physical groups
main_model.occ.synchronize()

# Find the volumes based on COM
for vol in gmsh.model.getEntities(dim=3):
        com = main_model.occ.getCenterOfMass(vol[0], vol[1])
        if np.allclose(com, [0, 0, 150]):
            main_model.addPhysicalGroup(3, [vol[1]], tag=obj1_tag, name="obj1")
        elif np.allclose(com, [0, 0, 50]):
            main_model.addPhysicalGroup(3, [vol[1]], tag=obj2_tag, name="obj2")
        else:
            main_model.addPhysicalGroup(3, [vol[1]], tag=tube_interior_tag, name="volume")
# Find the facets based on COM
for vol in gmsh.model.getEntities(dim=2):
        com = main_model.occ.getCenterOfMass(vol[0], vol[1])
        if np.allclose(com, [0, 0, 0]):
            main_model.addPhysicalGroup(2, [vol[1]], tag=bottom_disk_tag, name="bottom_disk")
        elif np.allclose(com, [0, 0, Tube_hight]):
            main_model.addPhysicalGroup(2, [vol[1]], tag=top_disk_tag, name="top_disk")
        elif np.allclose(com, [0, 0, Tube_hight / 2]):
            main_model.addPhysicalGroup(2, [vol[1]], tag=tude_sides_tag, name="side_surface")