Maxwell's Equations

[Peter Monk: "Finite Elements for Maxwell's Equations"]

Magnetostatic field generated by a permanent magnet

magnetic flux $B$, magnetic field $H$, given magnetization $M$:

$$ \DeclareMathOperator{\Grad}{grad} \DeclareMathOperator{\Curl}{curl} \DeclareMathOperator{\Div}{div} B = \mu (H + M), \quad \Div B = 0, \quad \Curl H = 0 $$

Introducing a vector-potential $A$ such that $B = \Curl A$, and putting equations together we get

$$ \Curl \mu^{-1} \Curl A = \Curl M $$

In weak form: Find $A \in H(\Curl)$ such that

$$ \int \mu^{-1} \Curl A \Curl v = \int M \Curl v \qquad \forall \, v \in H(\Curl) $$

Usually, the permeability $\mu$ is given as $\mu = \mu_r \mu_0$, with $\mu_0 = 4 \pi 10^{-7}$ the permeability of vacuum.

In [1]:
from ngsolve import *
from netgen.csg import *
import netgen.gui
%gui tk

Geometric model and meshing of a bar magnet:

In [2]:
box = OrthoBrick(Pnt(-3,-3,-3),Pnt(3,3,3)).bc("outer")
magnet = Cylinder(Pnt(-1,0,0),Pnt(1,0,0), 0.3) * OrthoBrick(Pnt(-1,-3,-3),Pnt(1,3,3))
air = box - magnet

geo = CSGeometry()
geo.Add (air.mat("air"))
geo.Add (magnet.mat("magnet").maxh(1))
mesh = Mesh(geo.GenerateMesh(maxh=2))
In [3]:
mesh.GetMaterials(), mesh.GetBoundaries()
Out[3]:
(['air', 'magnet'],
 ['outer',
  'outer',
  'outer',
  'outer',
  'outer',
  'outer',
  'default',
  'default',
  'default'])

Define space, forms and preconditioner. To obtain a regular system matrix, we add a very small $L_2$ term

In [4]:
fes = HCurl(mesh, order=3, dirichlet="outer", flags = { "nograds" : True })
print ("ndof =", fes.ndof)
u = fes.TrialFunction()
v = fes.TestFunction()

from math import pi
mu0 = 4*pi*1e-7
mur = CoefficientFunction( [1000 if mat== "magnet" else 1 
                            for mat in mesh.GetMaterials()]) 

a = BilinearForm(fes)
a += SymbolicBFI(1/(mu0*mur) * curl(u) * curl(v) + 1e-6/(mu0*mur)*u*v)
c = Preconditioner(a, "bddc")

f = LinearForm(fes)
mag = CoefficientFunction((1,0,0))
f += SymbolicLFI(mag*curl(v), definedon=mesh.Materials("magnet"))

ndof = 87940

Assemble system and setup preconditioner using task-parallelization:

In [5]:
with TaskManager():
    a.Assemble()
    f.Assemble()

Finally, declare GridFunction and solve by preconditioned CG iteration:

In [6]:
gfu = GridFunction(fes)

inv = CGSolver(a.mat, c.mat)
gfu.vec.data = inv * f.vec
Draw (curl(gfu), mesh, "B-field")

we forgot to curve the mesh, this can be done by

In [7]:
mesh.Curve(3)

Go back and assemble and solve again ...

In [8]: