Late this night … I can’t believe it, I found the error. One has to use:
sphere_inner = geom.add_ball([sphere_x1, sphere_y1, sphere_z1], sphere_r1)
sphere_outer = geom.add_ball([sphere_x2, sphere_y2, sphere_z2], sphere_r2)
sphere_final = geom.boolean_difference([sphere_outer],[sphere_inner], delete_first=True, delete_other=True)
geom.add_raw_code("Physical Surface(\"1\") = {1};") # inner
geom.add_raw_code("Physical Surface(\"2\") = {2};") # outer
geom.add_raw_code("Physical Volume(\"3\") = {" + sphere_final.id + "};") # volume
The numbers on the right side of Physical Surface(\"...\") were wrong, they are now correct (= {1} and = {2}).
Here is a complete code, which shall demonstrate, how to calculate the electrostatic potential/field of a spherical capacitor … . Have fun with the example and thanks again for the help!
# ##### BEGIN GPL LICENSE BLOCK #####
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software Foundation,
# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
#
# ##### END GPL LICENSE BLOCK #####
#
# The electrostatic potential and field of a spherical capacitor.
#
#
# Clemens Barth
#
#
# Start : 2020-04-01
# Last modified : 2020-04-01
#
#
import pygmsh
import meshio
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
import scipy.constants as const
# _____________________________________________________ The values to play with
# The potentials on the boundaries (inner and outer sphere)
pot_1 = 1.0 # inner
pot_2 = 0.0 # outer
# The position and size of the two spheres
sphere_x1 = 0.0
sphere_y1 = 0.0
sphere_z1 = 0.0
sphere_r1 = 10.0
sphere_x2 = 0.0
sphere_y2 = 0.0
sphere_z2 = 0.0
sphere_r2 = 100.0
# ____________________________________________________________________ The code
# The file name
filename_base = "Sphere-sphere-capacity_v06"
# Build the mesh first with pygmsh
geom = pygmsh.opencascade.Geometry(characteristic_length_min=10.0,
characteristic_length_max=10.0)
# Create the spheres and ...
sphere_inner = geom.add_ball([sphere_x1, sphere_y1, sphere_z1], sphere_r1)
sphere_outer = geom.add_ball([sphere_x2, sphere_y2, sphere_z2], sphere_r2)
# ... take the difference (boolean).
sphere_final = geom.boolean_difference([sphere_outer],[sphere_inner],
delete_first=True,
delete_other=True)
# The numbers 2, 3, ... 7 are identifiers of subsurfaces. They can be seen in
# Gmsh as follows (open the .geo file):
# 1. Key 'Alt-s' for surface, 'Alt-v' for volume
# 2. Hover over the dashed lines and read the values in the small rectangle.
geom.add_raw_code("Physical Surface(\"1\") = {1};") # inner sphere
geom.add_raw_code("Physical Surface(\"2\") = {2};") # outer sphere
geom.add_raw_code("Physical Volume(\"3\") = {" + sphere_final.id + "};") # vol.
# Build the .geo file.
mesh = pygmsh.generate_mesh(geom, geo_filename=filename_base+".geo")
# What follows is this (strange) mesh conversion stuff ... . ;-)
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(filename_base+"_tetra_mesh.xdmf", tetra_mesh)
meshio.write(filename_base+"_mf.xdmf", triangle_mesh)
# All the latter .xdmf files are read such that they can be used in FEniCS.
mesh = Mesh()
with XDMFFile(filename_base+"_tetra_mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile(filename_base+"_mf.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
# The mesh function.
V = FunctionSpace(mesh, 'CG', 2)
# Set the boundary conditions using the boundary numbers specified in Gmsh.
inner_boundary = DirichletBC(V, Constant(pot_1), mf, 1)
outer_boundary = DirichletBC(V, Constant(pot_2), mf, 2)
bcs =[inner_boundary, outer_boundary]
# Solve the Poisson equation with the source set to 0
# (the Laplace equation) via the L = Constant('0')... line
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v)) * dx
L = Constant('0') * v * dx
u = Function(V)
solve(a == L, u, bcs)
# Find the electric field by taking the negative of the
# gradient of the electric potential. Project this onto
# the function space for evaluation.
electric_field = project(-grad(u))
# output the results to two separate files
potentialFile = File(filename_base+"_potential.pvd")
potentialFile << u
e_fieldfile = File(filename_base+"_Efield.pvd")
e_fieldfile << electric_field
# Obtain a profile with matplotlib.
x = np.linspace(sphere_r1, sphere_r2-1.0, 1000)
y = np.zeros(len(x))
z = np.zeros(len(x))
coords = np.dstack((x, y, z))[0]
# This is the profile of the potential at x=x1..x2, 0, 0
u_line = np.array(list(map(u, coords)))
# We now calculate the analytical solution:
def Analytic_potential(phi_i, phi_o, r_i, r_o, r):
phi_1 = (phi_i * r_i - phi_o * r_o) / (r_i - r_o)
phi_2 = (r_i * r_o * (phi_i - phi_o)) / (r * (r_i - r_o))
phi = phi_1 - phi_2
return phi
# The line along the y axis.
y_analytic = np.linspace(sphere_r1, sphere_r2-1.0, 1000)
# The profile data.
u_line_analytic = [Analytic_potential(pot_1,
pot_2,
sphere_r1,
sphere_r2,
r) for r in y_analytic]
# Plot the potentials and errors
fig, (ax1, ax2) = plt.subplots(2, sharex=True)
fig.suptitle('Results: electrostatic potential of a spherical capacitor')
ax1.plot(y_analytic, u_line , linewidth=2)
ax1.plot(y_analytic, u_line_analytic, linewidth=2)
ax1.legend(['Pot_numeric', 'Pot_analytic'], loc='upper right')
ax2.plot(y_analytic, (u_line-u_line_analytic), linewidth=2)
ax2.legend(['Pot_numeric - Pot_analytic'], loc='upper right')
plt.savefig(filename_base+"_potential_profiles.png")
