I am trying to calculate the magnetic field of a permanent magnet in 3D (just for fun), in the region outside of the magnet (because I have read that in general a permanent magnet is ferromagnetic and the relation between its magnetization and the field is non linear, i.e. very complicated, which may matter if more than 1 magnet is considered at once).

I took the magnet as spherical and the surrounding space as being a box. Here is my gmsh file:

SetFactory(“OpenCASCADE”);

Mesh.RandomFactor=1.0e-4;

Sphere(1) = {0, 0, 0, 1};

Box(2) = {-4.9, -4.9, -4.9, 10, 10, 10};

BooleanDifference{ Volume{2}; Delete; }{ Volume{1}; Delete; }

Physical Volume(“testing_vol”) = {2};

Physical Surface(“sphere_surface_magnet”) = {1};

Physical Surface(“box_infinity_surfaces”) = {2, 6, 5, 4, 7, 3};

I generated the mesh with gmsh GUI and refined it there too. So it looks like the following picture

Even though it isn’t visible, there’s really a spherical empty space near the middle of the box, the region of the permanent magnet and where I am not interested to calculate the magnetic field.

I have then followed the magnetostatics FEniCS tutorial (Solving PDEs in Python - <br> The FEniCS Tutorial Volume I) and tried to extend it to 3D. In my case the \vec B field isn’t two dimensional, it is 3D. And the vector potential isn’t along z only (it could have a spherical symmetry because I picked a spherical magnet but I would rather not enforce such symmetry since I plan to modify the geometry. Furthermore it isn’t exactly centered in the box, so the solution shouldn’t be exactly symmetric with respect to rotations).

My FEniCS code (which works, i.e. returns something) is the following

```
import numpy as np
import meshio
msh = meshio.read("sphere_magnet_test.msh")
''' # Uncomment to create the files for the mesh to use with FEniCS.
cells = np.vstack(np.array([cells.data for cells in msh.cells if cells.type == "tetra"]))
tetra_data = msh.cell_data_dict["gmsh:physical"]["tetra"]
tetra_mesh = meshio.Mesh(points=msh.points, cells=[("tetra", cells)], cell_data={"name_to_read": [tetra_data]})
facet_cells = np.vstack([cells.data for cells in msh.cells if cells.type == "triangle"])
facet_data = msh.cell_data_dict["gmsh:physical"]["triangle"]
facet_mesh = meshio.Mesh(points=msh.points,
cells=[("triangle", facet_cells)],
cell_data={"name_to_read": [facet_data]})
# Write mesh
meshio.xdmf.write("mesh_test0.xdmf", tetra_mesh)
meshio.xdmf.write("facet_mesh_test0.xdmf", facet_mesh)
'''
#exit()
from dolfin import *
mesh = Mesh()
# Read mesh.
with XDMFFile("mesh_test0.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("facet_mesh_test0.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc2 = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("mesh_test0.xdmf") as infile:
infile.read(mvc2, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc2)
dx = Measure("dx", domain=mesh,subdomain_data=cf)
ds = Measure("ds", domain=mesh, subdomain_data=mf)
dS = Measure("dS", domain=mesh, subdomain_data=mf)
subdomains = cf
boundary_parts = mf
print('Expected volume:', 10**3-(4.0/3)*np.pi*1**3)
print('Computed volume: ',assemble(Constant(1)*dx(domain=mesh, subdomain_data=cf)))
print('computed volume again: ',assemble(Constant(1)*dx(domain=mesh, subdomain_data=cf, subdomain_id=1)))
print('Expected magnet surface:', 4*np.pi*1**2)
print('Magnet surface: ',assemble(Constant(1)*ds(domain=mesh, subdomain_data=mf, subdomain_id=2)))
print('Box expected surface: ', 10**2*6)
print('Box computed surface: ', assemble(Constant(1)*ds(domain=mesh, subdomain_data=mf, subdomain_id=3)))
# Define function space
V = FunctionSpace(mesh, 'P', 1)
# Define parameters
#magnetic_potential_sphere = Constant(0)
magnetic_potential_walls = Constant(0)
effective_current = Constant(1.0)
# Vacuum permeability.
mu_0 = 4*np.pi*1e-7
#bc0 = DirichletBC(V, magnetic_potential_walls, mf, 2)
bc1 = DirichletBC(V, magnetic_potential_walls, mf, 3)
#bcs = [bc0, bc1]
# Define variational problem
A_z = TrialFunction(V) # This doesn't work when the eq. is non linear.
A_x = TrialFunction(V)
A_y = TrialFunction(V)
v = TestFunction(V)
a = (1 / mu_0)*dot(grad(A_z), grad(v))*dx
L = effective_current * v * dx
a1 = (1 / mu_0)*dot(grad(A_x), grad(v))*dx
L1 = effective_current * v * dx
a2 = (1 / mu_0)*dot(grad(A_y), grad(v))*dx
L2 = effective_current * v * dx
# Solve variational problem
A_z = Function(V)
solve(a == L, A_z, bcs = bc1)
A_x = Function(V)
solve(a1 == L1, A_x, bcs = bc1)
A_y = Function(V)
solve(a2 == L2, A_y, bcs = bc1)
# Compute magnetic field (B = curl A)
W = VectorFunctionSpace(mesh, 'P', 1)
B = project(as_vector( (A_z.dx(1)-A_y.dx(2), -A_z.dx(0)+A_x.dx(2), A_y.dx(0)-A_x.dx(1)) ), W)
# Save solution to file
vtkfile_A_z = File('magnetostatics_test/potential_z.pvd')
vtkfile_A_x = File('magnetostatics_test/potential_x.pvd')
vtkfile_A_y = File('magnetostatics_test/potential_y.pvd')
vtkfile_B = File('magnetostatics_test/field.pvd')
vtkfile_A_z << A_z
vtkfile_A_x << A_x
vtkfile_A_y << A_y
vtkfile_B << B
```

As you can see, since the potential vector is really a 3D vector in my case, I solve 3 variational problems, one for each component of \vec A. **Is it a correct approach?** Or should have I used a mixed element approach instead? Or simply just a single variational problem without ME?

I would have expected the \vec B field to look like a dipole (just like the Earth’s magnetic field), but it doesn’t look like that, to me. Here’s a view of the field with Paraview:

To me, it looks like \vec B is circular around an axis. I probably picked wrong boundary conditions for the vector potential. If I assume that the magnet is polarized in a particular direction then I expect the boundary condition on \vec A to be different according the direction, which is not my case. In my case it’s as if the polarization of the magnet had the same x, y and z components, maybe explaining why the \vec B field “circulates” around an axis. I am not quite sure where I go wrong. I also assumed that the magnet’s polarization is equivalent to an effective current, i.e. as if the \vec B field was produced by a current density rather than by a magnetic material.Edit: By modyfing a little bit the code, more specifically I changed a part for this one instead:

```
magnetic_potential_walls = Constant(0)
magnetic_potential_sphere = Constant(1000000)
effective_current_z = Constant(0.0)
#effective_current_z = Expression('10000*atan(sqrt(x[0]*x[0]+x[1]*x[1])/x[2])',degree=2)
effective_current_x = Constant(0.0)
effective_current_y = Constant(0.0)
# Magnetic permeability.
mu_0 = 4*np.pi*1e-7
bc_walls = DirichletBC(V, magnetic_potential_walls, mf, 3)
bc_sphere_z = DirichletBC(V, magnetic_potential_sphere, mf, 2)
bc_sphere_x = DirichletBC(V, magnetic_potential_walls, mf, 2)
bc_sphere_y = DirichletBC(V, magnetic_potential_walls, mf, 2)
# Define variational problem
A_z = TrialFunction(V) # This doesn't work when the eq. is non linear.
#A_z = Function(V) # works for non linear equations.
A_x = TrialFunction(V)
A_y = TrialFunction(V)
v = TestFunction(V)
a = (1 / mu_0)*dot(grad(A_z), grad(v))*dx
L = effective_current_z * v * dx
a1 = (1 / mu_0)*dot(grad(A_x), grad(v))*dx
L1 = effective_current_x * v * dx
a2 = (1 / mu_0)*dot(grad(A_y), grad(v))*dx
L2 = effective_current_y * v * dx
# Solve variational problem
A_z = Function(V)
solve(a == L, A_z, bcs = [bc_sphere_z, bc_walls])
A_x = Function(V)
solve(a1 == L1, A_x, bcs = [bc_sphere_x, bc_walls])
A_y = Function(V)
solve(a2 == L2, A_y, bcs = [bc_sphere_y, bc_walls])
```

So that there’s no electric current at all anymore, only different values of the potential vector on the sphere vs on the box walls. I specified the anisotropy in A_z via the potential vector boundary condition on the sphere’s surface.

I get the following magnetic field:

Unfortunately it doesn’t look correct to me. The \vec B field should point either towards or outwards the sphere in the region closest to the sphere, not around it.