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.