Dear all.
Yep, and here is the solution 
, 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 geom.boolean_fragments, see the code at the end.
Cheers all.
Here is an image …
… and here is the code. Note that one can now cut under a specific angle … 0.0° < angle < 90.0°. If the angle is larger, the tags change and need to be adapted.
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 = 3.0
eps_domain = 1.0
FLAG_cut_vertic_4 = True
# The cutting angle, valid for: 0.0° < cut_angle < 90.0°
cut_angle = 15.0
# Domain
domain_size = 500.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_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 = 100.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 + 8.0)
# The height of the dielectric film
dielectric_vz = -dielectric_height
# We use opencascade for the 3D objects.
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 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([0.0 , -d, -d],
[2.0*d, 2.0*d, 2.0*d], char_length=1.0)
# Rotate box No. 2
geom.rotate(box_cut_2,
point=(0.0, 0.0, 0.0),
axis =(0.0, 0.0, 1.0),
angle=(180.0-cut_angle)/360.0 * (2.0*const.pi))
# First, we cut the dielectric film
dielect_final = geom.boolean_difference([dielect],[box_cut_1, box_cut_2],
delete_first=True,
delete_other=False)
# Now we cut the domain.
space_3 = geom.boolean_difference([space_2],[box_cut_1, box_cut_2],
delete_first=True,
delete_other=True)
# And here it comes: we have to use fractions such that we don't have
# multiple interfaces. This was the problem before!
geom.boolean_fragments([space_3],[dielect_final],
delete_first=True,
delete_other=True)
else:
# See just above.
geom.boolean_fragments([space_2],[dielect],
delete_first=True,
delete_other=True)
## This is for finding 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\") = {2};")
# Plane
geom.add_raw_code("Physical Surface(\"2\") = {7,8,9};")
# Domain
geom.add_raw_code("Physical Surface(\"3\") = {10};") # Sides: 1,3
# Dielectricum
geom.add_raw_code("Physical Surface(\"4\") = {4,5,6};") # Sides: 11,12
# Dielectricum interface
geom.add_raw_code("Physical Surface(\"5\") = {4,5,6};")
# Volume space
geom.add_raw_code("Physical Volume(\"6\") = {33};") # This needs to be determined manually!
# Volume dielectricum
geom.add_raw_code("Physical Volume(\"7\") = {32};") # This needs to be determined manually!
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\") = {11};")
# Volume dielectricum
geom.add_raw_code("Physical Volume(\"7\") = {10};")
# 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(-dielectric_height) + ";")
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])
print("aha")
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)
dx = Measure("dx", domain=mesh, subdomain_data=subdomains)
ds = Measure("ds", domain=mesh, subdomain_data=boundaries)
dS = Measure("dS", domain=mesh, subdomain_data=boundaries)
# ______________________________________ (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] = eps_dielec
else:
values[0] = eps_domain
return permittivity(subdomains, degree=1)
# Permittivity
eps = get_epsilon(subdomains)
# Mesh function
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(file_pot)
potentialFile << u