This page was generated from wta/maxwell.ipynb.

Maxwell Equations

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

Geometric model using Netgen constructive solid geometry:

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)

mesh = Mesh(geo.GenerateMesh(maxh=0.5))
mesh.Curve(5)
Draw (mesh, clipping = { "pnt" : (0,0,0), "vec" : (0,1,0) })

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)\).

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 = 159940

Distribute work by task-parallelization:

with TaskManager():
    a.Assemble()
    f.Assemble()
    solver = CGSolver(mat=a.mat, pre=c.mat)
    u.vec.data = solver * f.vec

Draw (curl(u), mesh, "B-field", draw_surf=False, \
      clipping = { "pnt" : (0,0,0), "vec" : (0,1,0), "function" : False },
      vectors = { "grid_size" : 100 })

BaseWebGuiScene