This page was generated from wta/maxwell.ipynb.

Maxwell Equations

[1]:
from netgen.csg import *
from ngsolve import *
from ngsolve.webgui import Draw

Geometric model using Netgen constructive solid geometry:

[2]:
def MakeGeometry():
    geometry = CSGeometry()
    box = OrthoBrick(Pnt(-1,-1,-1),Pnt(2,1,2)).bc("outer")

    core = OrthoBrick(Pnt(0,-0.05,0),Pnt(0.8,0.05,1))- \
           OrthoBrick(Pnt(0.1,-1,0.1),Pnt(0.7,1,0.9))- \
           OrthoBrick(Pnt(0.5,-1,0.4),Pnt(1,1,0.6)).maxh(0.2).mat("core")

    coil = (Cylinder(Pnt(0.05,0,0), Pnt(0.05,0,1), 0.3) - \
            Cylinder(Pnt(0.05,0,0), Pnt(0.05,0,1), 0.15)) * \
            OrthoBrick (Pnt(-1,-1,0.3),Pnt(1,1,0.7)).maxh(0.2).mat("coil")

    geometry.Add ((box-core-coil).mat("air"), transparent=True)
    geometry.Add (core, col=(0.5,0.5,0))
    geometry.Add (coil, col=(0,1,0))
    return geometry

geo = MakeGeometry()
# Draw (geo)
[3]:
mesh = Mesh(geo.GenerateMesh(maxh=0.5))
mesh.Curve(5)
Draw (mesh, clipping = { "pnt" : (0,0,0), "vec" : (0,1,0) })
 Start Findpoints
 main-solids: 3
 Found points 51
 Analyze spec points
 Find edges
 52 edges found
 Check intersecting edges
 CalcLocalH: 174 Points 0 Elements 0 Surface Elements
 Start Findpoints
 Analyze spec points
 Find edges
 52 edges found
 Check intersecting edges
 Start Findpoints
 Analyze spec points
 Find edges
 52 edges found
 Check intersecting edges
 Surface 1 / 22
 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
 52 elements, 226 points
 Surface 2 / 22
 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
 54 elements, 244 points
 Surface 3 / 22
 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
 90 elements, 278 points
 Surface 4 / 22
 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
 58 elements, 298 points
 Surface 5 / 22
 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
 52 elements, 315 points
 Surface 6 / 22
 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, 346 points
 Surface 7 / 22
 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
 45 elements, 346 points
 Surface 8 / 22
 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
 12 elements, 346 points
 Surface 9 / 22
 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
 45 elements, 346 points
 Surface 10 / 22
 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
 10 elements, 346 points
 Surface 11 / 22
 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
 8 elements, 346 points
 Surface 12 / 22
 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
 2 elements, 346 points
 Surface 13 / 22
 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
 10 elements, 346 points
 Surface 14 / 22
 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
 2 elements, 346 points
 Surface 15 / 22
 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
 10 elements, 346 points
 Surface 16 / 22
 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
 14 elements, 346 points
 Surface 17 / 22
 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
 8 elements, 346 points
 Surface 18 / 22
 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
 18 elements, 346 points
 Surface 19 / 22
 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
 87 elements, 375 points
 Surface 20 / 22
 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
 43 elements, 383 points
 Surface 21 / 22
 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
 156 elements, 448 points
 Surface 22 / 22
 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
 54 elements, 461 points
 SplitSeparateFaces
 Check subdomain 1 / 3
 Check subdomain 2 / 3
 Check subdomain 3 / 3

 Meshing subdomain 1 of 3
 Use internal rules
 Delaunay meshing
 number of points: 460
 blockfill local h
 number of points: 1303
 Points: 1303
 Elements: 7925
 Tree data entries per element: 1.53943
 Tree nodes per element: 0.0306625
 SwapImprove
 SwapImprove
 SwapImprove
 SwapImprove
 0 degenerated elements removed
 Remove intersecting
 Remove outer
 1303 points, 6886 elements
 start tetmeshing
 Use internal rules
 Success !
 1303 points, 6886 elements

 Meshing subdomain 2 of 3
 Use internal rules
 Delaunay meshing
 number of points: 1303
 blockfill local h
 number of points: 1309
 Points: 100
 Elements: 514
 Tree data entries per element: 1.36187
 Tree nodes per element: 0.0252918
 SwapImprove
 SwapImprove
 SwapImprove
 SwapImprove
 0 degenerated elements removed
 Remove intersecting
 Remove outer
 1303 points, 7024 elements
 start tetmeshing
 Use internal rules
 Success !
 1304 points, 7039 elements

 Meshing subdomain 3 of 3
 Use internal rules
 Delaunay meshing
 number of points: 1304
 blockfill local h
 number of points: 1306
 Points: 171
 Elements: 1111
 Tree data entries per element: 1.53015
 Tree nodes per element: 0.029703
 SwapImprove
 SwapImprove
 SwapImprove
 SwapImprove
 0 degenerated elements removed
 Remove intersecting
 Remove outer
 1306 points, 7490 elements
 start tetmeshing
 Use internal rules
 Success !
 1306 points, 7490 elements
 Remove Illegal Elements
 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 = 5
 Update clusters
 Curving edges
 Curving faces
[3]:
BaseWebGuiScene

Magnetostatics:

find \(u \in H(curl)\) such that

\[\int \mu^{-1} \operatorname{curl} u \operatorname{curl} v = \int j v\]

for all \(v \in H(curl)\).

[4]:
fes = HCurl(mesh, order=4, dirichlet="outer", nograds = True)
print ("ndof =", fes.ndof)
u,v = fes.TnT()

mur = { "core" : 1000, "coil" : 1, "air" : 1 }
mu0 = 1.257e-6
nu_coef = [ 1/(mu0*mur[mat]) for mat in mesh.GetMaterials() ]

nu = CoefficientFunction(nu_coef)
a = BilinearForm(fes, symmetric=True)
a += nu*curl(u)*curl(v)*dx + 1e-6*nu*u*v*dx

c = Preconditioner(a, type="bddc")

f = LinearForm(fes)
f += CoefficientFunction((y,0.05-x,0)) * v * dx("coil")

u = GridFunction(fes)
ndof = 182152

Distribute work by task-parallelization:

[5]:
with TaskManager():
    a.Assemble()
    f.Assemble()
    solver = CGSolver(mat=a.mat, pre=c.mat)
    u.vec.data = solver * f.vec
assemble VOL element 5969/5969
call wirebasket inverse ( with 6747 free dofs out of 182152 )
has inverse
assemble VOL element 5969/5969
0 3.53357e-05
1 1.71926e-05
2 1.25056e-05
3 1.04903e-05
4 8.75842e-06
5 7.17513e-06
6 5.19641e-06
7 3.95767e-06
8 2.7472e-06
9 1.91194e-06
10 1.35735e-06
11 9.68446e-07
12 6.65859e-07
13 4.86263e-07
14 3.53038e-07
15 2.43751e-07
16 1.66913e-07
17 1.17587e-07
18 8.61469e-08
19 6.06277e-08
20 4.43328e-08
21 3.19417e-08
22 2.24339e-08
23 1.61129e-08
24 1.17355e-08
25 7.84536e-09
26 5.58601e-09
27 3.85444e-09
28 2.73293e-09
29 1.87317e-09
30 1.29118e-09
31 9.11983e-10
32 6.57486e-10
33 4.59406e-10
34 3.18354e-10
35 2.25671e-10
36 1.54299e-10
37 1.06903e-10
38 7.39977e-11
39 5.16593e-11
40 3.70763e-11
41 2.6818e-11
42 1.8346e-11
43 1.25812e-11
44 8.87334e-12
45 6.00511e-12
46 4.1106e-12
47 2.88359e-12
48 1.98014e-12
49 1.35394e-12
50 9.60483e-13
51 6.67476e-13
52 4.52379e-13
53 3.05917e-13
[6]:
Draw (curl(u), mesh, "B-field", draw_surf=False, \
      clipping = { "pnt" : (0,0,0), "vec" : (0,1,0), "function" : False },
      vectors = { "grid_size" : 100 })
[6]:
BaseWebGuiScene
[ ]: