Magnetic field in 3D of a "permanent magnet"

[raboynics]

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 (https://fenicsproject.org/pub/tutorial/html/._ftut1015.html) 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.

[kamensky]

The linked example’s formulation for A_z is based on the assumption that fields are constant in the z-direction, and that the magnetic field is entirely in-plane. This is not the case in your example, and it is not valid to make this assumption for each component of \mathbf{A}, then superpose the results.

The standard approach for solving 3D magnetostatic problems is to solve the “curl-curl” problem for the vector potential \mathbf{A}, discretizing it with a special type of vector-valued finite element, called a Nédélec element. Nédélec function spaces are subsets of the space H(\text{curl}) of functions whose curls are square-integrable. However, unlike standard C^0 finite element spaces, Nédélec spaces are not subsets of H^1 (i.e., their full gradients are not square-integrable, only their curls). The practical implication of this is that tangential components of Nédélec functions are continuous across element faces, while normal components are not.

Clearly H^1-conforming functions are a subset of H(\text{curl}), but are not in general suitable for the curl-curl problem, because this subset is too restricted, as illustrated using FEniCS here, in the context of eigenvalue problems:

https://fenics.readthedocs.io/projects/dolfin/en/2017.2.0/demos/maxwell-eigenvalues/python/demo_maxwell-eigenvalues.py.html

For further discussion on subtleties of setting up curl-curl problems in FEniCS, you can take a look at this old discussion:

https://fenicsproject.org/qa/11451/nedelec-conforming-elements-converge-simple-test-problem/

[raboynics]

Thank you very much for replying!
This is above my head at the moment, I will try to see how to proceed.

However before I dive into it, I do not understand why the problem isn’t decoupled as I thought it is.
From Maxwell’s equation \nabla \times \vec H = \frac{\partial \vec D}{\partial t}+\vec J with \vec D=\vec 0 and assuming a linear relation between \vec H and \vec B (valid in vacuum or air, which is what matters to me, since this is where I solve the problem), one finds that \nabla \times \left( \frac{1}{\mu} \nabla \times \vec A \right)=\vec J . If \nabla \mu = \vec 0 (usually assumed to hold), then we get \nabla (\nabla \cdot \vec A)-\nabla ^2 \vec A=\mu \vec J where \nabla^2 is the “vector Laplacian” operator (this is eq.5.92 in Jackson’s Classical Electrodynamics textbook), so that \nabla ^2 \vec A= (\nabla^2 A_x, \nabla^2 A_y, \nabla^2 A_z). So this is a vector equation, to me it looks like every component is decoupled to any other. So solving this equation, to me, looks like to solve the equation for each component.

To me this leads to solve (assuming the Coulomb gauge \nabla \cdot \vec A =0)
\nabla^2 A_x=-\mu J_x for A_x (this solves the first component of the vector equation),

\nabla^2 A_y=-\mu J_y for A_y (this solves the second component of the vector equation),

and \nabla^2 A_z=-\mu J_z for A_z (this solves the third and last component of the vector equation).

Here, in these 3 equations, \nabla^2 is the (scalar) Laplacian.

I do not see where I go wrong. Could you please comment?

[kamensky]

If I understand your argument correctly, the logic is that \nabla\cdot\mathbf{A} = 0 implies that each component of \mathbf{A} satisfies an independent Poisson equation. However, if you solve three independent Poisson equations, I don’t see why the resulting vector field would be guaranteed to be solenoidal. Thus, the components would still be coupled through the constraint that \nabla\cdot\mathbf{A} = 0.

[nate]

For details consider Peter Monk’s excellent book, or from an engineering perspective, Jin’s excellent book.

[raboynics]

Thanks for the references, but I am not sure what I should look for, exactly. For instance in Jin’s textbook, it is mentioned that for 3D magnetostatics problems, it is indeed question of solving for the vector potential (and not the “curl curl” way mentioned by Kamensky). He mentions which boundary conditions to implement, but I have yet to get familiar with his notation. I suspect that if I could fix/tweak the boundary conditions in my code, it should work (i.e. return the correct result). But overall, I am still missing what I do wrong, exactly.

Regarding the first reference, I do not see anything related to magnetostatics in particular in the book. It has a part dealing with vector potentials in general, but I fail to see the relation with Maxwell’s equations. The “Coulomb” gauge for example is not mentioned (in neither book for that matter, though it is at least described in Jin’s).

So, yeah, I still do not see what is wrong with my approach. I suspect it’s just the boundary conditions, but that’s a guess.

[raboynics]

Essentially yes. This is written explicitely in Jackson’s textbook (p.181 of the 3rd edition). It should be Poisson equation (since the current density is not necessarily null), but in the case of no electric current this is indeed a Laplace equation for each of the components of the vector potential. I am not sure yet how to answer to your remark regarding the solenoidal condition. Good point to ponder on.

[Ricardo_Adriano]

Hello,
First of all, I would like to congratulate everyone involved in the project. My experience with Fenics has been very enjoyable. I am trying to create a model of a cylindrical magnet in 3D following the paper “Finite Element Analysis of Magnetic Fields with Permanent Magnet” (DOI: 10.1109/TMAG.1984.1063462). The 2D axisymmetric version worked perfectly and was validated with FEMM. However, in the 3D version, the results are correct only in the region outside the magnet. Inside the magnet, B tends to zero. An interesting fact is that if I assume the relative permeability (mu_r) of the magnet is equal to mu_0, removing the (1/mu) term from the bilinear form, both the 2D and 3D problems converge to the same solution.

It seems like I’m encountering an issue when assembling the matrix system, specifically when performing the product (1/mu_r) of type DG of order zero with the term curl(u)⋅curl(v), where u and v are of Nedelec type (first order). Any ideas on how to work around this problem?

The Strong Form of the Problem in \mathbb{R}^{n}, n = 2 or 3 is given by:

Given \mu_r, \mu_0, \textbf{J}_M, find A:\bar{\Omega}\rightarrow \mathbb{R}^{n} such that:

\nabla\times(\frac{1}{\mu_r}\nabla\times\textbf{A}) = \mu_0\textbf{J}_M; \; \forall x \in \Omega
\textbf{A} = 0; \; \forall x \in \Gamma
\textbf{J}_M = \nabla\times\textbf{M}

Where the magnetization vector \textbf{M} = H_c \hat{i}_p is defined by the coercivity of the magnetic material H_c and \hat{i}_p, which represents the direction of the intrinsic remanent magnetization.

Applying the weighted residual method and after some simplifications, the
n-dimensional weak form of the problem can be written as:

\int_\Omega \nabla \times \textbf{w} \cdot\frac{1}{\mu_r}\nabla \times \textbf{A} d\Omega = \mu_0\int_\Omega \textbf{M} \cdot \nabla\times\textbf{w} d\Omega

As presented by Kamensky, the standard approach for solving 3D magnetostatic problems for the vector potential A involves using Nédélec elements. As a consequence, I defined the variational form as:

#ldim = 1 # element order.
#N = functionspace(mesh, (“N1curl”, ldim))
#u = TrialFunction(N), v = TestFunction(N)

#dx = Measure(“dx”, domain=mesh, subdomain_data=cell_tags);

#a = (1/mu_r)dot(curl(u),curl(v)) * dx #
#L = mu0inner(M_,curl(v)) * dx(1)

where mu_r = mr inside the magnet and 1.0 outside. To represent the material, I used the Q = functionspace(mesh, (“DG”, 0)) and the physical groups defined in Gmsh.

To validate the problem, I used a cilindrical magnet in free space. Considering that the magnet and the magnetization vector are oriented along the z-direction, the problem can be solved by exploiting axial symmetry. Thus, the axisymmetric weak form of the problem can be written as:

\int_\Omega \nabla \times \textbf{w} \cdot\frac{1}{\mu_r}\nabla \times \textbf{A} d\Omega = 2\pi \int_{\Omega_{2D}}\frac{1}{\mu_r r} (\nabla w'\cdot \nabla A')\partial r \partial z

and

\mu_0\int_\Omega \textbf{M} \cdot \nabla\times\textbf{w} d\Omega = 2\pi\mu_0 \int_{\Omega_{2D}} M \cdot \nabla\times w' \partial r \partial z

In this case, the vector quantity \textbf{A} and \textbf{w} reduces to the scalars w' = rw_\phi and A'= rA_\phi and the 2D variational form is defined as:

#ldim = 2 # Lagrange element order
#V = functionspace(mesh, (“Lagrange”, ldim))

#x = SpatialCoordinate(mesh)
#r = x[0]

#u = TrialFunction(V), v = TestFunction(V)

#a = (1 / mu)(1 / r) * dot(grad(u), grad(v)) * dx
#L = mu0M_* v.dx(0) * dx

To compare the results, the magnetic flux along the z-axis (r=0) was plotted using 101 points. The results were compared with each other and with the FEMM software (HomePage:Finite Element Method Magnetics).

It can be observed that the 2D formulation and FEMM yield similar results. However, the 3D formulation provides good results only in the region outside the magnet. By exporting the results to ParaView, it is possible to see that this behavior occurs not only along the line (r=0) but also throughout the entire domain.

Fig. 1 Magnetic Induction (B) for mr =1200, (xy) 2d result, (xz) 3d result

However, if I set mr = 1, the 2D and 3D problems converge to the same solution.

Fig. 2 Magnetic Induction (B) for mr =1, (xy) 2d result, (xz) 3d result

The problem I need to solve is neither axisymmetric nor linear (though I am not addressing that part at the moment). Therefore, I cannot use the 2D approximation, and I need a straightforward way to update \mu_r in the future to incorporate the nonlinearities. I have definitely exhausted my ideas on how to work around this problem, and any help is welcome.

See below my MWE for the 2D and 3D cases. I am running the programs via Docker on Windows using the dolfinx/lab:nightly image.

# Geometry parameters
#air box
rair = 80e-3
zair = 80e-3 
# Magnet
rmag = 30e-3
zmag = 30e-3
 
# magnet permeability
#mi0 = 4*np.pi*1e-7
mr = 1200 # if mr=1 -> 3d problem == 2d problem !?
# coercivity 
Hc  = 1000

import numpy as np # For vector and matrix manipulation
import gmsh # Interface with the software GMSH

# Visualization 
import pyvista
from mpi4py import MPI
import matplotlib.pyplot as plt

# FEniCS
from ufl import (SpatialCoordinate, TestFunction, TrialFunction, as_vector,
                 dx, dot, grad, curl, inner, Measure)
from petsc4py import PETSc

from dolfinx import default_scalar_type, geometry
from dolfinx.mesh import create_mesh, compute_midpoints, locate_entities_boundary
from dolfinx.io import gmshio, XDMFFile, VTKFile, VTXWriter
from dolfinx.io.gmshio import model_to_mesh
from dolfinx.fem import (Constant, dirichletbc, Function, Expression, functionspace, assemble_scalar,
                         form, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.plot import vtk_mesh
from dolfinx import fem

2D Problem

2D Mesh

rank = MPI.COMM_WORLD.rank

gmsh.initialize()

gdim = 2  # Geometric dimension of the mesh
model_rank = 0
mesh_comm = MPI.COMM_WORLD

if mesh_comm.rank == model_rank:

    # Creates the air rectangle
    air_box = gmsh.model.occ.addRectangle(0, -zair/2, 0, rair, zair)
    print(air_box)
    # Creates the magnet
    magnet = gmsh.model.occ.addRectangle(0, -zmag/2, 0, rmag, zmag)
    print(magnet)

    whole_domain = gmsh.model.occ.fragment([(gdim, air_box)], [(gdim, magnet)])
    gmsh.model.occ.synchronize()
    print(whole_domain)
    
    # Creating physical groups
    gmsh.model.addPhysicalGroup(gdim, [magnet], tag = 1, name = "mag_tag") 
    gmsh.model.addPhysicalGroup(gdim, [3], tag = 2, name = "air_tag") 
    gmsh.model.occ.synchronize()
    
    #gmsh.option.setNumber("Mesh.CharacteristicLengthMax", hcil/8)
    
    #set mesh size of the magnet
    gmsh.model.mesh.field.add("Box", 1)
    gmsh.model.mesh.field.setNumber(1, "VIn", zmag/40)
    gmsh.model.mesh.field.setNumber(1, "VOut", zmag/5)
    gmsh.model.mesh.field.setNumber(1, "XMin", 0)
    gmsh.model.mesh.field.setNumber(1, "YMin", -zmag/2)
    gmsh.model.mesh.field.setNumber(1, "XMax", rmag)
    gmsh.model.mesh.field.setNumber(1, "YMax", zmag/2)
    gmsh.model.mesh.field.setNumber(1, "Thickness", zmag)
    gmsh.model.mesh.field.setAsBackgroundMesh(1)
    gmsh.model.occ.synchronize()

    gmsh.model.mesh.generate(3)  
    gmsh.write("MWE2D.msh")

# Converts the mesh created in GMSH to the format used by Dolfinx.
mesh,cell_tags,facet_tags = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
# Close GMSH
gmsh.finalize()

# Sets up the visualization with PyVista
pyvista.start_xvfb()

plotter = pyvista.Plotter()
mesh.topology.create_connectivity(gdim, gdim)
grid = pyvista.UnstructuredGrid(*vtk_mesh(mesh, gdim))
num_local_cells = mesh.topology.index_map(gdim).size_local
grid.cell_data["Marker"] = cell_tags.values[cell_tags.indices < num_local_cells]
grid.set_active_scalars("Marker")

# Applies the "viridis" colormap to color the mesh based on "Cell Markers"
actor = plotter.add_mesh(grid, show_edges=True, cmap="viridis", edge_color="black")
actor = plotter.add_title(
    'Fig. 1 Geometria do problema axissimétrico', font='courier', color='k', font_size=7)
plotter.view_xy()
plotter.show();
print('Fig. 3 Geometria do problema axissimétrico')

2D Weak form

Dirichlet boundary condition (A=0)

ldim = 2 # Lagrange element order
V = functionspace(mesh, ("Lagrange", ldim))
tdim = mesh.topology.dim

# Define the radial coordinate r from x[0]
x = SpatialCoordinate(mesh)
r = x[0]

# Dirichlet boundary condition (A=0)
a0 = fem.Function(V)
a0.x.array[:] = 0

# Find the boundary nodes and define that all of them belong to the Dirichlet boundary (A=0).
facets = locate_entities_boundary(mesh, tdim - 1, lambda x: np.full(x.shape[1], True))
dofs = locate_dofs_topological(V, tdim - 1, facets)
bc = dirichletbc(a0, dofs)
del a0

Magnet relative permeability (\mu_r)

Q = functionspace(mesh, ("DG", 0))
mu_r = Function(Q) #

# Permeability
mu0 = Constant(mesh, PETSc.ScalarType(4e-7*np.pi))
muf = Constant(mesh, PETSc.ScalarType(mr))
muair = Constant(mesh, PETSc.ScalarType(1.0))

# List of elements in each domain
mag_elements = cell_tags.find(1) 
air_elements = cell_tags.find(2)

# Defining the permeability values
mu_r.x.array[mag_elements] = np.full_like(mag_elements, muf, dtype=default_scalar_type)
mu_r.x.array[air_elements] = np.full_like(air_elements, muair, dtype=default_scalar_type)

Magnetization vector

Magnetization vector in \hat{z}.

# Magnitude of the magnetization vector
M = Function(Q)
M.x.array[:] = 0
M.x.array[mag_elements] = np.full_like(mag_elements, Hc, dtype=default_scalar_type)

M_ = Function(V)
M_expr = Expression(M, V.element.interpolation_points())
M_.interpolate(M_expr)
M.x.scatter_forward()
del M

Variational equation with cylindrical coordinates (r,\phi,z)

# Trial and test functions
u = TrialFunction(V)
v = TestFunction(V)

a =  (1 / mu_r)*(1 / r) * dot(grad(u), grad(v)) * dx #(1 / mu) *
L = mu0*M_* v.dx(0) * dx  # Derivando v

Solving the problem

A_ = Function(V) # A_ = r*A_alpha
problem = LinearProblem(a, L, u=A_, bcs=[bc])
problem.solve()

2D Results

FEMM result

Reference result using FEMM using mr = 1200, r=0 and -zair/2 < z < zair/2

# FEMM solution
B_FEMM = np.array([1.17233668e-008,1.49526914e-005,2.84316238e-005,4.20871235e-005,5.95092533e-005,7.23426298e-005
,8.50584574e-005,9.61686499e-005,1.07590921e-004,1.21075141e-004,1.34489307e-004,1.47598483e-004,1.60755857e-004
,1.74069628e-004,1.86817760e-004,1.98332071e-004,2.09553259e-004,2.20286245e-004,2.30171191e-004,2.38957291e-004
,2.48202146e-004,2.57906187e-004,2.66789677e-004,2.75030398e-004,2.81140276e-004,2.85939748e-004,2.90481422e-004
,2.94900139e-004,2.98646986e-004,3.02158925e-004,3.04542986e-004,3.05146192e-004,3.54558549e-004,4.10396674e-004
,4.59355861e-004,5.10392938e-004,5.69194487e-004,6.20014955e-004,6.70664157e-004,7.10190363e-004,7.43641753e-004
,7.79413698e-004,8.19134381e-004,8.53232677e-004,8.75003397e-004,8.96628941e-004,9.13815713e-004,9.28699512e-004
,9.40257373e-004,9.45045110e-004,9.49829434e-004,9.44744232e-004,9.39657871e-004,9.28635948e-004,9.14692392e-004
,8.97228223e-004,8.72980645e-004,8.48687213e-004,8.20222470e-004,7.91757726e-004,7.55895231e-004,7.16599272e-004
,6.71897791e-004,6.17161012e-004,5.62153858e-004,5.14125346e-004,4.61298812e-004,4.11124740e-004,3.54136261e-004
,3.05068464e-004,3.04472052e-004,3.02266048e-004,2.99541504e-004,2.96585109e-004,2.92967587e-004,2.88217399e-004
,2.81056377e-004,2.73254521e-004,2.64400152e-004,2.56069486e-004,2.48404135e-004,2.40325597e-004,2.31841239e-004
,2.22486555e-004,2.12460656e-004,1.99948830e-004,1.85945822e-004,1.73391468e-004,1.61496379e-004,1.47364804e-004
,1.32474841e-004,1.19206232e-004,1.06323795e-004,9.51739291e-005,8.42862728e-005,7.11248706e-005,5.77956739e-005
,4.20408863e-005,2.78460894e-005,1.40591593e-005,1.12339595e-008])

Magnetic Induction (B)

W = functionspace(mesh, ("DG", 0, (mesh.geometry.dim, )))
B_2d = Function(W)
B_expr = Expression(as_vector((-(1/r)*A_.dx(1), (1/r)*A_.dx(0))), W.element.interpolation_points())
B_2d.interpolate(B_expr)

#Creating a balanced tree with the elements of the mesh
bb_tree = geometry.bb_tree(mesh, mesh.topology.dim) 

# sample points
tol = 1e-6
z_points = np.linspace(-zair/2 + tol, zair/2-tol, 101)
r_points = np.full_like(z_points, 0)
points = np.zeros((3, 101))
points[0] = r_points
points[1] = z_points


B_2d_values = []
cells = []
points_on_proc = []

# Find cells whose bounding-box collide with the the points
cell_candidates = geometry.compute_collisions_points(bb_tree, points.T)
# Choose one of the cells that contains the point
colliding_cells = geometry.compute_colliding_cells(mesh, cell_candidates, points.T)
for i, point in enumerate(points.T):
    if len(colliding_cells.links(i)) > 0:
        points_on_proc.append(point)
        cells.append(colliding_cells.links(i)[0])

points_on_proc = np.array(points_on_proc, dtype=np.float64)
B_2d_values = B_2d.eval(points_on_proc, cells)

# plot horizontal line
fig = plt.figure(figsize=(8,5))
plt.plot(z_points, B_2d_values[:,0], 'y', linewidth=2); # magnify w
plt.plot(z_points, B_2d_values[:,1], 'g', linewidth=2); # magnify w
plt.plot(z_points, B_FEMM, '-.b', linewidth=2); # magnify w
plt.grid(True);
plt.xlabel('$z$');
plt.legend(['Br_2d (T)','Bz_2d (T)','|B_femm| (T)' ], loc='upper right');
fig.suptitle("2D results $B(z)$");

# Saiving xdmf
geo = mesh.geometry

xdmf = XDMFFile(MPI.COMM_WORLD, "B2d.xdmf", "w")
xdmf.write_mesh(mesh)
# Write the MeshTags object to the file
xdmf.write_meshtags(cell_tags, geo)
xdmf.write_function(B_2d)
#xdmf.write_function(mu)
xdmf.close()

del (geo, mesh, cell_tags, facet_tags, mu_r, A_,Q ,V , B_2d, M_, x, r, 
     bb_tree, points_on_proc, cell_candidates, dx, a, L, bc, dofs)

Fig. 3 Comparing the 2D result with FEMM for mr =1200.

3D Problem

3D Mesh

gdim = 3  # Geometric dimension of the mesh

gmsh.initialize()

if mesh_comm.rank == model_rank:

    # Creates the air box
    air_box = gmsh.model.occ.addBox(-rair, -rair, -zair/2, 2*rair, 2*rair, zair)
    
    # Creates the magnet
    magnet = gmsh.model.occ.addCylinder(0, 0, -zmag/2, 0, 0, zmag, rmag)

    whole_domain = gmsh.model.occ.fragment([(gdim, air_box)], [(gdim, magnet)])
    gmsh.model.occ.synchronize()
    
    # Creating physical groups
    gmsh.model.addPhysicalGroup(gdim, [magnet], tag = 1, name = "mag_tag") 
    gmsh.model.addPhysicalGroup(gdim, [3], tag = 2, name = "air_tag") 
    gmsh.model.occ.synchronize()
    
    #gmsh.option.setNumber("Mesh.CharacteristicLengthMax", hcil/8)
    
    #set mesh size of the magnet
    gmsh.model.mesh.field.add("Box", 1)
    gmsh.model.mesh.field.setNumber(1, "VIn", zmag/20)
    gmsh.model.mesh.field.setNumber(1, "VOut", zmag/5)
    gmsh.model.mesh.field.setNumber(1, "XMin", -rmag)
    gmsh.model.mesh.field.setNumber(1, "YMin", -rmag)
    gmsh.model.mesh.field.setNumber(1, "ZMin", -zmag/2)
    gmsh.model.mesh.field.setNumber(1, "XMax", rmag)
    gmsh.model.mesh.field.setNumber(1, "YMax", rmag)
    gmsh.model.mesh.field.setNumber(1, "ZMax", zmag/2)
    gmsh.model.mesh.field.setNumber(1, "Thickness", zmag)
    gmsh.model.mesh.field.setAsBackgroundMesh(1)
    gmsh.model.occ.synchronize()

    gmsh.model.mesh.generate(3)  
    gmsh.write("MWE3D.msh")

# Converts the mesh created in GMSH to the format used by Dolfinx.
mesh,cell_tags,facet_tags = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)

# Close GMSH
gmsh.finalize()
del gmsh

Dirichlet boundary condition (A=0)

ldim = 1 # element order
N = functionspace(mesh, ("N1curl", ldim))
tdim = mesh.topology.dim

x = SpatialCoordinate(mesh) 

a0 = fem.Function(N)
a0.x.array[:] = 0

# Find the boundary nodes and define that all of them belong to the Dirichlet boundary (A=0).
facets = locate_entities_boundary(mesh, tdim - 1, lambda x: np.full(x.shape[1], True))
dofs = locate_dofs_topological(N, tdim - 1, facets)
bc = dirichletbc(a0, dofs)
del a0

Magnet relative permeability (\mu_r)

Q = functionspace(mesh, ("DG", 0))
mu_r = Function(Q) #

# Permeability
mu0 = Constant(mesh, PETSc.ScalarType(4e-7*np.pi))
mur_iron = Constant(mesh, PETSc.ScalarType(mr))
mur_air = Constant(mesh, PETSc.ScalarType(1.0))

# List of elements in each domain
mag_elements = cell_tags.find(1) 
air_elements = cell_tags.find(2)

# Defining the permeability values
mu_r.x.array[mag_elements] = np.full_like(mag_elements, muf, dtype=default_scalar_type)
mu_r.x.array[air_elements] = np.full_like(air_elements, muair, dtype=default_scalar_type)
mu_r.x.scatter_forward()

Magnetization vector

# Magnitude of the magnetization vector
M = Function(Q)
M.x.array[air_elements] = np.full_like(air_elements, 0.0, dtype=default_scalar_type)
M.x.array[mag_elements] = np.full_like(mag_elements, Hc, dtype=default_scalar_type)
M.x.scatter_forward()

M_ = fem.Function(N)
M_expr = Expression(as_vector((0.0, 0.0, M)), N.element.interpolation_points())
M_.interpolate(M_expr)
M_.x.scatter_forward()

del M

3d Variational equation

# Trial and test functions
u = TrialFunction(N)
v = TestFunction(N)

dx = Measure("dx", domain=mesh, subdomain_data=cell_tags);
 
a =  (1/mu_r)*inner(curl(u),curl(v)) * dx #
L = mu0*inner(M_,curl(v)) * dx(1)

Solving the problem

A = Function(N)
problem = LinearProblem(a, L, u = A, bcs = [bc], petsc_options = {"ksp_type": "gmres", "ksp_rtol":1e-9, "ksp_atol":1e-13, "ksp_max_it": 4000, "pc_type": "none"})
problem.solve()
del M_

3D results

W = functionspace(mesh, ("DG", 0, (mesh.geometry.dim, )))

B_3d = Function(W)
B_exp = Expression(curl(A), W.element.interpolation_points())
B_3d.interpolate(B_exp)

#Creating a balanced tree with the elements of the mesh
bb_tree = geometry.bb_tree(mesh, mesh.topology.dim) 

# sample points
tol = 1e-6
z_points = np.linspace(-zair/2 + tol, zair/2-tol, 101)
r_points = np.full_like(z_points, 0)
points = np.zeros((3, 101))
points[0] = r_points
points[1] = r_points
points[2] = z_points

B_3d_values = []
cells = []
points_on_proc = []

# Find cells whose bounding-box collide with the the points
cell_candidates = geometry.compute_collisions_points(bb_tree, points.T)
# Choose one of the cells that contains the point
colliding_cells = geometry.compute_colliding_cells(mesh, cell_candidates, points.T)
for i, point in enumerate(points.T):
    if len(colliding_cells.links(i)) > 0:
        points_on_proc.append(point)
        cells.append(colliding_cells.links(i)[0])

points_on_proc = np.array(points_on_proc, dtype=np.float64)
B_3d_values = B_3d.eval(points_on_proc, cells)

# plot horizontal line
fig = plt.figure(figsize=(8,5))
plt.plot(z_points, B_3d_values[:,0], 'y', linewidth=2); # magnify w
plt.plot(z_points, B_3d_values[:,1], 'b', linewidth=2); # magnify w
plt.plot(z_points, B_3d_values[:,2], 'g', linewidth=2); # magnify w
plt.plot(z_points, B_2d_values[:,1], '-.b', linewidth=2); # magnify w
plt.grid(True);
plt.xlabel('$z$');
plt.legend(['Bx_3d (T)','By_3d (T)', 'Bz_3d (T)','Bz_2d'], loc='upper right');
fig.suptitle("$B(z)$");

# Saiving xdmf
geo = mesh.geometry

xdmf = XDMFFile(MPI.COMM_WORLD, "B3d.xdmf", "w")
xdmf.write_mesh(mesh)
# Write the MeshTags object to the file
xdmf.write_meshtags(cell_tags, geo)
xdmf.write_function(B_3d)
#xdmf.write_function(mu)
xdmf.close()

Fig. 4 Comparing the 2D and 3D results for mr =1200.