Magnetic field of a cylindrical coil with an imported mesh

Hi FEniCS-Community,

I’m very new to FEniCS and I’m trying to simulate the magnetic flux of a cylindrical air coil to get started. I choose this setup because it is possible to compare the results with a Biot-Savart simulation. Later iron will be added to the simulation. The code is similar to the ft11-magnetostatics.py example from the tutorial. I tried to calculate every component of the vector potential independently:

-\nabla \cdot (\mu^{-1} \nabla A_x) = J_x
-\nabla \cdot (\mu^{-1} \nabla A_y) = J_y
-\nabla \cdot (\mu^{-1} \nabla A_z) = J_z

For simplicity I reduced my code here to the calculation of one component. Unfortunately the results are not as expected. The magnetic flux is zero everywhere except the domains where there is a current density. Do I need to add an additional boundary condition? Or is there another problem with my code or my mesh?

I would be very happy if someone finds the mistake/s!

Here you can find the mesh:

The domains were numbered with Gmsh. (Cell index inside the copper of the coil: 4)

I’m using:
Ubuntu 18.04.3
Python 3.6.9
dolfin version 2019.1.0

from fenics import *
from mshr import *

import matplotlib.pyplot as plt

R = 200   #Outer radius air sphere in mm
r = 50    #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")

# Define function space and integration measure
V = FunctionSpace(mesh, "P", 1)
dx = Measure('dx', domain=mesh, subdomain_data=markers)
dV = Measure("dx")(mesh)

# 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 = dot(grad(Ax), grad(vx))*dV

# right hand sides of variational problems
Lx = mu*Jx*vx*dV

# solve variational problems for magnetic Ax (in Vs/mm)
Ax = Function(V)
solve(ax == Lx, Ax, bc)

# plot solutions
plt.figure()
plot(Ax)
plt.show()

I recently spent some time getting the magnetostatics demo working on dolfinx, and was happy with the results. I think I’d be able to help, would you be able to post a link to the mesh that is accessible currently (right now it leads to a 404)? Or even better do you have a simplified geometry that you are able to test using pygmsh and meshio?

Hi awarru,
Thank you for yor reply! Here you can finde my mesh now:


I created this mesh using Gmsh. I hope this is simple enough. Or do you suggest to simplify it even further?
I wasn’t able to solve my problems with FEniCs, so I siwtched to COMSOL instead. Nevertheless I would be very interested in the solution.

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!

Thank you very much! This helps a lot! I’m considering now to go back to FEniCs again.
The reason for the negativ permeability in the copper domain was a comment in the FEniCs Tutorial. On page 104 of the Tutorial they define the domains for the magnetostatic example. In an older Version of the tutorial the permeability of the copper was negative and there was written: “yes, it’s really negative!”. It always confused me, but I copied it anyway.