Hello.
I’m new to FEniCS and would like to use its fantastic capabilities. I use FEniCS under Lubuntu 18.04 in a docker environment and have therefore FEniCS version 2019.1.0. In the container, I have also installed python3-pip and gmsh, and with pip I installed pygmsh meshio h5py lxml.
I followed the example described here (Example), which deals with the electrostatic potential of a sphere-sphere capacitor. The example works very well! So far so good. 
Now, I want to create my personal 3D capacitor (a 3D sphere in front of a plane) with the help of pygmsh and got stuck in the description of the boundary conditions. Here is how I start the code:
import pygmsh
import meshio
import numpy as np
from fenics import *
filename_base = "Create_sphere-plane_v02"
geom = pygmsh.opencascade.Geometry(characteristic_length_min=0.3,
characteristic_length_max=0.3)
sphere_x = 0.0
sphere_y = 0.0
sphere_z = 0.0
sphere_r = 1.0
plane_size = 30.0
plane_thickness = 0.5
plane_z = -4.0
plane_x1 = -plane_size/2.0
plane_y1 = -plane_size/2.0
plane_z1 = plane_z
plane_x2 = plane_size
plane_y2 = plane_size
plane_z2 = plane_thickness
sphere = geom.add_ball([sphere_x, sphere_y, sphere_z], sphere_r)
rectangle = geom.add_box([plane_x1, plane_y1, plane_z1],
[plane_x2, plane_y2, plane_z2])
Now, I want to label the two surfaces, the entire one of the sphere (1) and of the plane (2). I have read that one needs to use geom.add_raw_code, but how?
geom.add_raw_code(‘What comes inside here?’)
geom.add_raw_code(‘What comes inside here?’)
After this code, I have read that one has to well-prepare the mesh, which I got from here:
mesh = pygmsh.generate_mesh(geom)
tetra_cells = None
for cell in mesh.cells:
if cell.type == "tetra":
# Collect the individual meshes
if tetra_cells is None:
tetra_cells = cell.data
else:
tetra_cells = np.vstack((tetra_cells, cell.data))
tetra_data = None
for key in mesh.cell_data_dict["gmsh:physical"].keys():
if key == "tetra":
tetra_data = mesh.cell_data_dict["gmsh:physical"][key]
# Create tetra mesh with cell data attached
tetra_mesh = meshio.Mesh(points=mesh.points, cells={"tetra": tetra_cells},
cell_data={"name_to_read":[tetra_data]})
meshio.write(filename_base+".xdmf", tetra_mesh)
This code works well (without the geom.add_raw_code thing from above), like also the following code needed to read the geometry for FEniCS:
mesh_file = XDMFFile(MPI.comm_world, filename_base+".xdmf")
mesh = Mesh()
mesh_file.read(mesh);
mvc = MeshValueCollection("size_t", mesh, 3)
mesh_file.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
File(filename_base+".pvd").write(mesh)
So, the questions are:
-
What must I insert into
geom.add_raw_code('What comes inside here?')such that I can assign the entire surface of the sphere to label ‘1’ and the surface of the plane to label ‘2’? -
How can I assign these two identifiers, 1 and 2, to two Dirichlet boundaries conditions where the sphere is on potential 1 and the plane on potential 0. So, I would like to something like:
boundaries = … ?
inner_boundary = fn.DirichletBC(V, fn.Constant(1.0), boundaries, 1)
outer_boundary = fn.DirichletBC(V, fn.Constant(0.0), boundaries, 2)
Thanks for some ideas!
