This page was generated from unit-2.4-Maxwell/Maxwell.ipynb.

2.4 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.

[1]:

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.csg import *

Geometric model and meshing of a bar magnet:

[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"), transparent=True)
geo.Add (magnet.mat("magnet").maxh(1), col=(0.3,0.3,0.1))
geo.Draw()
mesh = Mesh(geo.GenerateMesh(maxh=2, curvaturesafety=1))
mesh.Curve(3)

 Calc Triangle Approximation
 Object 0 has 3596 triangles
 Object 1 has 2984 triangles
 Start Findpoints
 main-solids: 2
 Found points 24
 Analyze spec points
 Find edges
 14 edges found
 Check intersecting edges
 CalcLocalH: 44 Points 0 Elements 0 Surface Elements
 Start Findpoints
 Analyze spec points
 Find edges
 14 edges found
 Check intersecting edges
 Start Findpoints
 Analyze spec points
 Find edges
 14 edges found
 Check intersecting edges
 Surface 1 / 9
 load internal triangle rules
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 20 elements, 49 points
 Surface 2 / 9
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 26 elements, 57 points
 Surface 3 / 9
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 26 elements, 65 points
 Surface 4 / 9
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 30 elements, 75 points
 Surface 5 / 9
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 22 elements, 81 points
 Surface 6 / 9
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 22 elements, 87 points
 Surface 7 / 9
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 84 elements, 123 points
 Surface 8 / 9
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 6 elements, 124 points
 Surface 9 / 9
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 6 elements, 125 points
 SplitSeparateFaces
 Check subdomain 1 / 2
 Check subdomain 2 / 2

 Meshing subdomain 1 of 2
 Use internal rules
 Delaunay meshing
 number of points: 122
 blockfill local h
 number of points: 398
 Points: 398
 Elements: 2431
 Tree data entries per element: 1.56314
 Tree nodes per element: 0.0308515
 SwapImprove
 SwapImprove
 SwapImprove
 SwapImprove
 0 degenerated elements removed
 Remove intersecting
 Remove outer
 398 points, 2136 elements
 start tetmeshing
 Use internal rules
 Success !
 398 points, 2136 elements

 Meshing subdomain 2 of 2
 Use internal rules
 Delaunay meshing
 number of points: 398
 blockfill local h
 number of points: 404
 Points: 56
 Elements: 310
 Tree data entries per element: 1.29032
 Tree nodes per element: 0.0225806
 SwapImprove
 SwapImprove
 SwapImprove
 SwapImprove
 0 degenerated elements removed
 Remove intersecting
 Remove outer
 398 points, 2256 elements
 start tetmeshing
 Use internal rules
 Success !
 398 points, 2256 elements
 Remove Illegal Elements
 SplitImprove
 SwapImprove
 SwapImprove2
 SplitImprove
 SwapImprove
 SwapImprove2
 SplitImprove
 SwapImprove
 SwapImprove2
 SplitImprove
 SwapImprove
 SwapImprove2
 SplitImprove
 SwapImprove
 SwapImprove2
 SplitImprove
 SwapImprove
 SwapImprove2
 SplitImprove
 SwapImprove
 SwapImprove2
 SplitImprove
 SwapImprove
 SwapImprove2
 SplitImprove
 SwapImprove
 SwapImprove2
 SplitImprove
 SwapImprove
 SwapImprove2
 SplitImprove
 SwapImprove
 SwapImprove2
 SplitImprove
 SwapImprove
 SwapImprove2
 SplitImprove
 SwapImprove
 SwapImprove2
 SplitImprove
 SwapImprove
 SwapImprove2
 SplitImprove
 SwapImprove
 SwapImprove2
 SplitImprove
 SwapImprove
 SwapImprove2
 Volume Optimization
 CombineImprove
 ImproveMesh
 SplitImprove
 ImproveMesh
 SwapImprove
 SwapImprove2
 ImproveMesh
 CombineImprove
 ImproveMesh
 SplitImprove
 ImproveMesh
 SwapImprove
 SwapImprove2
 ImproveMesh
 CombineImprove
 ImproveMesh
 SplitImprove
 ImproveMesh
 SwapImprove
 SwapImprove2
 ImproveMesh
 Update mesh topology
 Update clusters
 Curve elements, order = 3
 Update clusters
 Curving edges
 Curving faces

[3]:

mesh.GetMaterials(), mesh.GetBoundaries()

[3]:

(('air', 'magnet'),
 ('outer',
  'outer',
  'outer',
  'outer',
  'outer',
  'outer',
  'default',
  'default',
  'default'))

Define space, forms and preconditioner.

To obtain a regular system matrix, we regularize by adding a very small \(L_2\) term.
We solve magnetostatics, so we can gauge by adding and arbitrary gradient field. A cheap possibility is to delete all basis-functions which are gradients (flag ‘nograds’)

[4]:

fes = HCurl(mesh, order=3, dirichlet="outer", nograds=True)
print ("ndof =", fes.ndof)
u,v = fes.TnT()

from math import pi
mu0 = 4*pi*1e-7
mur = mesh.MaterialCF({"magnet" : 1000}, default=1)

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

f = LinearForm(fes)
mag = mesh.MaterialCF({"magnet" : (1,0,0)}, default=(0,0,0))
f += mag*curl(v) * dx("magnet")

ndof = 25068

Assemble system and setup preconditioner using task-parallelization:

[5]:

with TaskManager():
    a.Assemble()
    f.Assemble()

assemble VOL element 1736/1736
call wirebasket inverse ( with 1975 free dofs out of 25068 )
has inverse
assemble VOL element 1736/1736

Finally, declare GridFunction and solve by preconditioned CG iteration:

[6]:

gfu = GridFunction(fes)
with TaskManager():
    solvers.CG(sol=gfu.vec, rhs=f.vec, mat=a.mat, pre=c.mat)

WARNING: maxsteps is deprecated, use maxiter instead!
CG iteration 1, residual = 0.004832539962846939
CG iteration 2, residual = 0.00288019846536514
CG iteration 3, residual = 0.0022636840727236693
CG iteration 4, residual = 0.0018596783600588406
CG iteration 5, residual = 0.0013706814384952773
CG iteration 6, residual = 0.0009629315332397556
CG iteration 7, residual = 0.0006784248173953639
CG iteration 8, residual = 0.0004537347999527588
CG iteration 9, residual = 0.00033739845183924656
CG iteration 10, residual = 0.00026636366952020196
CG iteration 11, residual = 0.0001578586571740973
CG iteration 12, residual = 0.00010596525796416685
CG iteration 13, residual = 7.440611744102437e-05
CG iteration 14, residual = 5.02177621171484e-05
CG iteration 15, residual = 3.5654571375899266e-05
CG iteration 16, residual = 2.1652775596128424e-05
CG iteration 17, residual = 1.6241088605443216e-05
CG iteration 18, residual = 1.0173867259497692e-05
CG iteration 19, residual = 7.247797850199135e-06
CG iteration 20, residual = 4.782087681570744e-06
CG iteration 21, residual = 3.2977746967721143e-06
CG iteration 22, residual = 2.314739047913969e-06
CG iteration 23, residual = 1.4401117530131444e-06
CG iteration 24, residual = 1.0156213478907965e-06
CG iteration 25, residual = 6.934274995206427e-07
CG iteration 26, residual = 5.493246097613419e-07
CG iteration 27, residual = 4.914301547620254e-07
CG iteration 28, residual = 2.841056622455584e-07
CG iteration 29, residual = 1.8374004730767134e-07
CG iteration 30, residual = 1.257610087828429e-07
CG iteration 31, residual = 8.212388417052288e-08
CG iteration 32, residual = 5.655902172431422e-08
CG iteration 33, residual = 3.7186008616013925e-08
CG iteration 34, residual = 2.3956217081010416e-08
CG iteration 35, residual = 1.649825747950453e-08
CG iteration 36, residual = 1.0759738128442831e-08
CG iteration 37, residual = 7.16506361378016e-09
CG iteration 38, residual = 4.497597041625025e-09
CG iteration 39, residual = 2.9156780301960117e-09
CG iteration 40, residual = 1.946000533781764e-09
CG iteration 41, residual = 1.4015945762644607e-09
CG iteration 42, residual = 8.454807217595495e-10
CG iteration 43, residual = 5.885393082899384e-10
CG iteration 44, residual = 3.7150504607750433e-10
CG iteration 45, residual = 2.803666043881062e-10
CG iteration 46, residual = 2.82955532866504e-10
CG iteration 47, residual = 1.5462884796474538e-10
CG iteration 48, residual = 1.088397301258785e-10
CG iteration 49, residual = 8.237486132359566e-11
CG iteration 50, residual = 5.651518020100518e-11
CG iteration 51, residual = 3.6673148626309044e-11
CG iteration 52, residual = 2.366165356457601e-11
CG iteration 53, residual = 1.632212552898413e-11
CG iteration 54, residual = 1.0677638890920086e-11
CG iteration 55, residual = 7.155951527893395e-12
CG iteration 56, residual = 4.796321545484079e-12
CG iteration 57, residual = 3.09894113850954e-12
CG iteration 58, residual = 2.0500215089870894e-12
CG iteration 59, residual = 1.3892675374369965e-12
CG iteration 60, residual = 8.94075998045078e-13
CG iteration 61, residual = 5.967419014352893e-13
CG iteration 62, residual = 3.903708686929611e-13
CG iteration 63, residual = 2.5785823123257414e-13
CG iteration 64, residual = 1.7705720785169274e-13
CG iteration 65, residual = 1.3638568850931996e-13
CG iteration 66, residual = 1.2821182031781367e-13
CG iteration 67, residual = 6.787573804378798e-14
CG iteration 68, residual = 4.557191301691235e-14
CG iteration 69, residual = 3.0159055477099586e-14
CG iteration 70, residual = 2.061343702346506e-14
CG iteration 71, residual = 1.5063560355198106e-14
CG iteration 72, residual = 1.0246751803582759e-14
CG iteration 73, residual = 6.269147514048469e-15
CG iteration 74, residual = 4.1749772564773254e-15

[7]:

# the vector potential is not supposed to look nice
Draw (gfu, mesh, "vector-potential", draw_surf=False, clipping=True)

Draw (curl(gfu), mesh, "B-field", draw_surf=False, clipping=True)
Draw (1/(mu0*mur)*curl(gfu)-mag, mesh, "H-field", draw_surf=False, clipping=True)

[7]:

BaseWebGuiScene

[ ]: