This page was generated from unit-11.1-fmm/BiotSavart.ipynb.

11.1.4 Biot Savart

Let \(G(x,y) = \frac{1}{4 \pi} \frac{\exp (i k |x-y|)}{|x-y|}\) be Green’s function for the Helmholtz equation.

For a given current path \(j\) along a Curve \(C\), the magnetic field in vacuum on the full space \({\mathbb R}^3\) is

\[H(x) = \int_C j(y) \times \nabla_y G(x,y) dl_y\]

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

vismesh = Mesh(OCCGeometry(Box((-5,-5,-5), (5,5,5))).GenerateMesh(maxh=1))
for l in range(0):
    vismesh.Refine()
Draw (vismesh)

WebGLScene

from ngsolve.bem import BiotSavartCF, BiotSavartSingularMLCF, BiotSavartRegularMLCF

We evaluate the single layer integral using numerical integration on the surface mesh. Thus, we get a sum of many Green’s functions, which is compressed using a multilevel-multipole.

kappa = 0.01*pi
mp = BiotSavartSingularMLCF((0,0,0), r=5, kappa=kappa)
mp.expansion.AddCurrent( (1,0,-1), (1,0,1), 1, num=100)
mp.expansion.AddCurrent( (1,0, 1), (-1,0,1), 1, num=100)
mp.expansion.AddCurrent( (-1,0, 1), (-1,0,-1), 1, num=10)
mp.expansion.AddCurrent( (-1,0, -1), (1,0,-1), 1, num=100)

regmp = mp.CreateRegularExpansion((0,0,0),r=5)

clipping = { "function" : False,  "pnt" : (0,0,0), "vec" : (0,0,-1) }
Draw (regmp.real, vismesh, min=0, max=1, order=2,  vectors={"grid_size" : 40, "offset" : 0 }, clipping=clipping);

visplane = WorkPlane(Axes( (0,0.2,0), Z, X)).RectangleC(5,5).Face()
vismesh2 = Mesh(OCCGeometry(visplane).GenerateMesh(maxh=0.5))
Draw (regmp.real[0], vismesh2, min=-0.1, max=0.1, order=10);

kappa = 0.01*pi
mp = BiotSavartSingularMLCF((0,0,0), r=5, kappa=kappa)

coil = Cylinder((0,-1,0), Y, r=1, h=1, mantle="outer") - Cylinder((0,-1,0), Y, r=0.5, h=1)
coilmesh = Mesh(OCCGeometry(coil).GenerateMesh(maxh=0.3)).Curve(3)
Draw (coilmesh)

current = CF((z,0,-x))
current /= Norm(current)
mp.expansion.AddCurrentDensity(current, coilmesh.Materials(".*"))
# mp.expansion.AddCurrentDensity(current, coilmesh.Boundaries("outer"))

regmp = mp.CreateRegularExpansion((0,0,0),r=5)

clipping = { "function" : False,  "pnt" : (0,0,0), "vec" : (0,0,-1) }
Draw (regmp.real, vismesh, min=0, max=0.2, order=2,  vectors={"grid_size" : 40, "offset" : 0 }, clipping=clipping);

Draw (regmp.real[0], vismesh2, min=-0.1, max=0.1, order=10);

the spikes inside the coil domain stem from numerical integration of the current source.

from ngsolve.webgui import FieldLines
N=14
fl = FieldLines(mp.real, vismesh.Materials('.*'), num_lines=N**3/20, length=1)

clipping = { "function" : False,  "pnt" : (0,0,0), "vec" : (0,0,-1) }
Draw (regmp.real, vismesh, min=0, max=0.5, order=2,  vectors={"grid_size" : 40, "offset" : 0 }, clipping=clipping);

settings = {"objects": { "Surface":False, "Wireframe":False}}
clipping = {"function" : False, "pnt" : (0,0,0), "vec" : (0,0,-1) }
vectors = {"grid_size" : 100, "offset": 0.2}
Draw (mp.real, vismesh, "X", objects=[fl], min=0, max=0.5, autoscale=False, settings=settings, vectors=vectors, clipping=clipping);

Biot-Savart with MaxwellDL integral operator

from ngsolve.bem import *

fescur = HDiv(coilmesh)
u,v = fescur.TnT()
gfcur = GridFunction(fescur)
gfcur.Set(current)

potCF = MaxwellDL(u*dx, 1e-4)(gfcur)
potCF.BuildLocalExpansion(vismesh2.Boundaries(".*"))

Draw(potCF[0], vismesh2, min=-0.1, max=0.1, order=10);

Biot-Savart with a line current

wirebox = Box((-2.5,-2.5,-2.5), (2.5,2.5,2.5))
wp = WorkPlane(Axes((0,0,0), Y, X))
wp.Circle(1)
wire = wp.Wires()[0]
for e in wire.edges:
    e.name = "wire"
    e.maxh = 0.05
wiremesh = Mesh(OCCGeometry(Glue([wirebox, wire])).GenerateMesh(maxh=1))
Draw(wiremesh)

WebGLScene

fesline = HCurl(wiremesh, order=2)
uline = fesline.TrialFunction()
gfline = GridFunction(fesline)
gfline.Set(current, definedon=wiremesh.BBoundaries("wire"))

potline = MaxwellDL(uline*dl(definedon=wiremesh.BBoundaries("wire")), 1e-4)(gfline)

# H_y = I R^2 / (2 (R^2+y^2)^(3/2)) (radius R, current I)
for yy in [0, 0.5, 1, 2]:
    Hnum = potline(wiremesh(0,yy,0))[1].real
    Hex = 1/(2*(1+yy**2)**1.5)
    print (f"y = {yy:3.1f}:   H_y = {Hnum:.5f}   exact: {Hex:.5f}")

y = 0.0:   H_y = 0.50005   exact: 0.50000
y = 0.5:   H_y = 0.35776   exact: 0.35777
y = 1.0:   H_y = 0.17674   exact: 0.17678
y = 2.0:   H_y = 0.04470   exact: 0.04472

potline.BuildLocalExpansion(vismesh2.Boundaries(".*"))
Draw (potline[0].real, vismesh2, min=-1, max=1, order=10);
Draw (potline[1].real, vismesh2, min=-1, max=1, order=10);