Dear all.
My problem is probably due to a mesh thing, this is why I post here. Well, I don’t know … we will see. By the way, the exemplifying code is at the end of this message!
I model 3D capacitors with Gmsh (pymsh OPENCASCADE, meshio, etc.) and calculate different electrostatic properties of the capacitors with FEniCS, and all this works very well (small example see here)! I could even implement MPI such that I can run the calculations on several cores. So far so good. The objective now is to cut the scene such that only a small percentage of the otherwise large domain is included in the calculations. This massively reduces computing time.
In the following image, a typical sphere(top)-dielectric film(middle)-plane(bottom) capacitor is shown:
This capacitor is inside a large domain:
Here is the potential of the full domain (2D cut, through the center, ParaView):
The voltages of -2 and +5V are applied at the plane and sphere, respectively, with a Dirichlet BC. As one can see, the electric field is somewhat entering into the dielectric film (see weak contrast change in the dielectric just below the sphere). Everything is alright!
Now, I cut the scene, which is possible thanks to the high symmetry around z. And thanks to pygmsh and OPENCASCADE I can boolean cut this large domain onto a 4th quite easily:
However, the potential is not anymore well calculated, see the following image:
The dielectric film appears as a uniform sharp box, as I would apply a Dirichlet boundary condition. However, I do not apply a DBC!
I though that FEniCS automatically applies Neumann boundary conditions onto facets etc. if one does not explicitly assign a condition for the facets. In the first example (no domain cutting) it is indeed the case! However, as soon as I cut the scene, this does not work anymore for the dielectric film.
Comments
(a) The labels for the facets and volumes (physical identifiers) are all correct for the gmsh, I checked that a dozen of times. (b) I tried to explicitly apply a Neumann boundary condition on the facets of the dielectric film. No luck, well, may be I did it wrong! (c ) It is not the eps assignment, this works.
Do you have an idea what is happening here and what might be the solution? Thanks for some comments in advance.
Here the code:
import pygmsh
import meshio
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
import scipy.constants as const
import os
# __________________________________________________________________ Parameters
z = 5.0
V_sphere = 5.0
V_plane = -2.0
eps_domain = 1.0
eps_dielec = 10.0
FLAG_cut_vertic_4 = False
# Domain
domain_size = 500.0
domain_eps = 1.0
# Sphere
sphere_r = 15.0
sphere_x = 0.0
sphere_y = 0.0
sphere_z = 0.0
# Dielectric film
dielectric_size = 80.0
dielectric_height = 40.0
dielectric_eps = 15.0
dielectric_x = 0.0
dielectric_y = 0.0
dielectric_z = 0.0
dielectric_vx = 0.0
dielectric_vy = 0.0
dielectric_vz = 0.0
# Conducting plane
plane_size = 80.0
plane_x = 0.0
plane_y = 0.0
plane_z = 0.0
plane_vx = 0.0
plane_vy = 0.0
plane_vz = 10.0
# ___________________________________________________________________ Directory
text_spaces = " "
dir_base = "./" + os.path.basename(__file__)[:-3]
if not os.path.isdir(dir_base):
os.mkdir(dir_base)
# ____________________________________________________ Mesh creation via pygmsh
sphere_z = sphere_r + z
plane_z = -(plane_vz + dielectric_height + 0.001)
# The position and size of the dielectric film
dielectric_vz = -dielectric_height
# Build the meshs first with pygmsh and ...
geom = pygmsh.opencascade.Geometry()
domain = geom.add_ball([0.0, 0.0, 0.0], domain_size)
sphere = geom.add_ball([sphere_x, sphere_y, sphere_z], sphere_r)
plane = geom.add_cylinder([plane_x, plane_y, plane_z],
[plane_vx, plane_vy, plane_vz],
plane_size, 2.0*const.pi)
dielect = geom.add_cylinder([dielectric_x, dielectric_y, dielectric_z],
[dielectric_vx, dielectric_vy, dielectric_vz],
dielectric_size, 2.0*const.pi)
# ... take the difference (boolean).
space_1 = geom.boolean_difference([domain],[sphere],
delete_first=True,
delete_other=True)
space_2 = geom.boolean_difference([space_1],[plane],
delete_first=True,
delete_other=True)
# We want the dielectricum inside a 'hole' of exact same size inside the domain,
# therefore keep the dielectricum.
space_3 = geom.boolean_difference([space_2],[dielect],
delete_first=True,
delete_other=False)
# We vertically cut out a 4th of the scene.
if FLAG_cut_vertic_4:
d = domain_size*1.5
box_cut_1 = geom.add_box([-0.0 , -d, -d],
[2.0*d, 2.0*d, 2.0*d],
char_length=1.0)
box_cut_2 = geom.add_box([-d, 0.0, -d],
[2.0*d, 2.0*d, 2.0*d],
char_length=1.0)
# Domain
space_4 = geom.boolean_difference([space_3],[box_cut_1, box_cut_2],
delete_first=True,
delete_other=False)
# Dielectric
dielect_final = geom.boolean_difference([dielect],[box_cut_1, box_cut_2],
delete_first=True,
delete_other=True)
## The .geo file for Gmsh => find the identifiers ...
#mesh = pygmsh.generate_mesh(geom,
# geo_filename=os.path.join(dir_base,
# "base.geo"))
#exit(-1)
if FLAG_cut_vertic_4:
# Sphere
geom.add_raw_code("Physical Surface(\"1\") = {29};")
# Plane
geom.add_raw_code("Physical Surface(\"2\") = {33,34,35};")
# Domain
geom.add_raw_code("Physical Surface(\"3\") = {26};") # Sides: 27, 28
# Dielectricum
geom.add_raw_code("Physical Surface(\"4\") = {36,38,40};") # Sides: 37, 39
# Dielectricum interface
geom.add_raw_code("Physical Surface(\"5\") = {30,31,32};")
# Volume space
geom.add_raw_code("Physical Volume(\"6\") = {" + space_4.id + "};")
# Volume dielectricum
geom.add_raw_code("Physical Volume(\"7\") = {" + dielect_final.id + "};")
else:
# Sphere
geom.add_raw_code("Physical Surface(\"1\") = {10};")
# Plane
geom.add_raw_code("Physical Surface(\"2\") = {11,12,13};")
# Domain
geom.add_raw_code("Physical Surface(\"3\") = {9};")
# Dielectricum
geom.add_raw_code("Physical Surface(\"4\") = {6,7,8};")
# Dielectricum interface
geom.add_raw_code("Physical Surface(\"5\") = {6,7,8};")
# Volume space
geom.add_raw_code("Physical Volume(\"6\") = {" + space_3.id + "};")
# Volume dielectricum
geom.add_raw_code("Physical Volume(\"7\") = {" + dielect.id + "};")
# Some quick mesh refinement:
geom.add_raw_code("Field[1] = Box;")
geom.add_raw_code("Field[1].VIn = " + str(5.0) + ";")
geom.add_raw_code("Field[1].VOut = " + str(50.0) + ";")
geom.add_raw_code("Field[1].XMax = " + str(dielectric_size) + ";")
geom.add_raw_code("Field[1].XMin = " + str(-dielectric_size) + ";")
geom.add_raw_code("Field[1].YMax = " + str(dielectric_size) + ";")
geom.add_raw_code("Field[1].YMin = " + str(-dielectric_size) + ";")
geom.add_raw_code("Field[1].ZMax = " + str(sphere_r*3.0) + ";")
geom.add_raw_code("Field[1].ZMin = " + str(0.0) + ";")
geom.add_raw_code("Background Field = 1;")
# Build the .geo file.
file_geo = os.path.join(dir_base, "Geometry.geo")
mesh = pygmsh.generate_mesh(geom, geo_filename=file_geo)
# ______________________________________________ Mesh conversion and re-reading
tetra_cells = None
tetra_data = None
triangle_cells = None
triangle_data = None
for key in mesh.cell_data_dict["gmsh:physical"].keys():
if key == "triangle":
triangle_data = mesh.cell_data_dict["gmsh:physical"][key]
elif key == "tetra":
tetra_data = mesh.cell_data_dict["gmsh:physical"][key]
for cell in mesh.cells:
if cell.type == "tetra":
if tetra_cells is None:
tetra_cells = cell.data
else:
tetra_cells = np.vstack([tetra_cells, cell.data])
elif cell.type == "triangle":
if triangle_cells is None:
triangle_cells = cell.data
else:
triangle_cells = np.vstack([triangle_cells, cell.data])
tetra_mesh = meshio.Mesh(points=mesh.points,
cells={"tetra": tetra_cells},
cell_data={"name_to_read":[tetra_data]})
triangle_mesh =meshio.Mesh(points=mesh.points,
cells=[("triangle", triangle_cells)],
cell_data={"name_to_read":[triangle_data]})
meshio.write(os.path.join(dir_base,"tetra_mesh.xdmf"), tetra_mesh)
meshio.write(os.path.join(dir_base,"mf.xdmf") , triangle_mesh)
# All the latter .xdmf files are read such that they can be used in FEniCS.
mesh = Mesh()
with XDMFFile(os.path.join(dir_base,"tetra_mesh.xdmf")) as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile(os.path.join(dir_base,"mf.xdmf")) as infile:
infile.read(mvc, "name_to_read")
boundaries = cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc2 = MeshValueCollection("size_t", mesh, 3)
with XDMFFile(os.path.join(dir_base,"tetra_mesh.xdmf")) as infile:
infile.read(mvc2, "name_to_read")
subdomains = cpp.mesh.MeshFunctionSizet(mesh, mvc2)
# ______________________________________ (Sub)domains, boundary conditions, etc.
# The epsilon in divers considered media (vacuum, dielectric, etc.)
def get_epsilon(subdomains):
# Here, we define the eps for the thin film (dielectric_eps) and the vacuum.
class permittivity(UserExpression):
def __init__(self, markers, **kwargs):
super(permittivity, self).__init__(**kwargs)
self.markers = markers
def eval_cell(self, values, x, cell):
if self.markers[cell.index] == 7:
values[0] = dielectric_eps
else:
values[0] = domain_eps
return permittivity(subdomains, degree=1)
# Permittivity
eps = get_epsilon(subdomains)
# Mesh function
V0 = FunctionSpace(mesh, 'DG', 0)
V = FunctionSpace(mesh, 'P', 2)
# The Dirichlet boundary conditions for the sphere and plane:
b_cond = []
b_cond.append(DirichletBC(V, Constant(V_sphere), boundaries, 1))
b_cond.append(DirichletBC(V, Constant(V_plane), boundaries, 2))
# We put the potential at the domain boundary onto 0V:
b_cond.append(DirichletBC(V, Constant(0.0), boundaries, 3))
# Test and trial functions
u = TrialFunction(V)
v = TestFunction(V)
# Create a and L.
a = dot(grad(u), grad(v)) * eps * dx
L = Constant(0.0) * v * dx
u = Function(V)
solve(a == L, u, b_cond,
solver_parameters = {"linear_solver" : "gmres",
"preconditioner" : "amg"})
if FLAG_cut_vertic_4:
file_pot = os.path.join(dir_base, "Potential_cut.pvd")
else:
file_pot = os.path.join(dir_base, "Potential_uncut.pvd")
potentialFile = File(MPI.comm_self, file_pot)
potentialFile << u
, the issue was a mesh thing: if one has two or more ‘connected’ domains, which shall be created by ‘cutting’, then one has to look forward that the interface has only one well-defined facet (before I always had 2 !) To do so, one needs to use 