I spent a while trying to get this working but wasn’t able to with your mesh. I still don’t know what the problem with the mesh was, but I was able to get it working by using pygmsh to generate another simple solenoid-like coil.
The mesh generation will require gmsh (4.5.6), meshio, and pygmsh. A lot of this code was adapted from this thread describing how to use pygmsh and meshio. I’d take a look there if you want to use pygmsh in the future. Here’s the mesh generation code:
from pygmsh.opencascade import Geometry
from pygmsh import generate_mesh
import meshio
import numpy as np
#Domain Sphere properties (centered)
domain_radius = 20.0
# The position and size of the two cylinders
length = 15.0
rect1_r = 2.0
rect2_r = .5
# Build the mesh first with pygmsh
geom = Geometry()
rect1 = geom.add_box([-rect1_r/2,-rect1_r/2,-length/2], [rect1_r,rect1_r,length], char_length=0.2)
rect2 = geom.add_box([-rect2_r/2,-rect2_r/2,-length/2], [rect2_r,rect2_r,length], char_length=0.2)
solenoid = geom.boolean_difference([rect1], [rect2])
#Create the domain
sphere = geom.add_ball([0,0,0], domain_radius, char_length = 4.0)
sphere_domain = geom.boolean_difference([sphere], [solenoid], delete_other=False)
geom.boolean_fragments([solenoid], [sphere_domain], delete_first = False, delete_other=False)
# 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 Volume(4) = {{{}}};".format(solenoid.id)) # Torus 1
geom.add_raw_code("Physical Volume(2) = {{{}}};".format(sphere_domain.id))
geom.add_raw_code("Physical Surface(3) = {{{}}};".format("17")) # Sphere border
# Build the .geo file.
mesh = generate_mesh(geom, geo_filename="mesh.geo", )
# What follows is this (strange) mesh conversion stuff ... . ;-)
tetra_cells = None
tetra_data = None
triangle_cells = None
triangle_data = None
print(mesh.cell_data_dict["gmsh:physical"].keys())
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.xdmf.write("tetra_mesh.xdmf", tetra_mesh)
meshio.xdmf.write("triangle_mesh.xdmf", triangle_mesh)
For the FEniCs code, I changed the permeability of copper to be non-negative (was there a reason it was negative originally?) and also borrowed from the linked thread above to read in the mesh and mesh boundaries. I also added code to save an XDMF file to read the output in Paraview as opposed to matplotlib, I’ve found it to be much easier to use especially when using the quay.io docker images.
Here’s the modified code to generate A_x:
from fenics import *
from mshr import *
from mpi4py import MPI
import dolfin
import matplotlib.pyplot as plt
domain_radius = 19
# The position and size of the two cylinders
length = 15.0
rect1_r = 2.0
rect2_r = .5
R = domain_radius #Outer radius air sphere in mm
r = 1 #Radius, coil
J = 1 #Current density in A/m^2
#mesh = Mesh("airCoilMesh.xml")
#markers = MeshFunction("size_t", mesh,"airCoilMesh_physical_region.xml")
#boundaries = MeshFunction("size_t", mesh,"airCoilMesh_facet_region.xml")
mesh = Mesh()
with XDMFFile("tetra_mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("tetra_mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
markers = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc)
# Define function space and integration measure
V = FunctionSpace(mesh, "P", 1)
dV = Measure('dx', domain=mesh, subdomain_data=markers)
# Define boundary condition
class Boundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return x[0]**2+x[1]**2+x[2]**2 >= R**2-tol
boundary = Boundary()
bc = DirichletBC(V, Constant(0), boundary)
# define current densities as function of subdomains
class CurrentDensity_x(UserExpression):
def __init__(self, markers, **kwargs):
super().__init__(**kwargs)
self.markers = markers
def eval_cell(self, values, x, cell):
if (self.markers[cell.index] == 4):
values[0] = -x[1]/r*J
else:
values[0] = Constant(0.0)
# permeability
class Permeability(UserExpression):
def __init__(self, markers, **kwargs):
super().__init__(**kwargs)
self.markers = markers
def eval_cell(self, values, x, cell):
if self.markers[cell.index] == 4:
values[0] = 6.4e-6 # copper
else:
values[0] = 4*pi*1e-7 # vacuum
mu = Permeability(markers, degree=1)
# initialise current densities
JX = CurrentDensity_x(markers, degree = 1)
Jx = project(JX, V)
# define variational problems for x, y, z
Ax = TrialFunction(V)
vx = TestFunction(V)
# left hand sides of variational problems
ax = (1/mu)*dot(grad(Ax), grad(vx))*dV
# right hand sides of variational problems
Lx = Jx*vx*dV
# solve variational problems for magnetic Ax (in Vs/mm)
Ax = Function(V)
solve(ax == Lx, Ax, bc)
with XDMFFile(MPI.COMM_WORLD,"output.xdmf") as outfile:
outfile.write(mesh)
outfile.write_checkpoint(Ax, 'Ax')
I’m in a similar boat to you, I could switch to using COMSOL but I want to try to see if an open-source tool could work. Hope this helps!