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

2.4.1 Maxwell eigenvalue problem

We solve the Maxwell eigenvalue problem

\[\int \operatorname{curl} u \, \operatorname{curl} v = \lambda \int u v\]

for \(u, v \; \bot \; \nabla H^1\) using a PINVIT solver from the ngsolve solvers module.

The orthogonality to gradient fields is important to eliminate the huge number of zero eigenvalues. The orthogonal sub-space is implemented using a Poisson projection:

\[P u = u - \nabla \Delta^{-1} \operatorname{div} u\]

The algorithm and example is take form the Phd thesis of Sabine Zaglmayr, p 145-150.

[1]:

import netgen.gui
from netgen.csg import *
from ngsolve import *

geom = CSGeometry()

cube1 = OrthoBrick(Pnt(-1,-1,-1), Pnt(1,1,1))
cube2 = OrthoBrick(Pnt(0,0,0), Pnt(2,2,2))
fichera = cube1-cube2
geom.Add (fichera)

# all edges defined by intersection of faces of
# cube2 and cube2 are marked for geometric edge refinement
geom.SingularEdge (cube2,cube2, 1)
geom.SingularPoint (cube2,cube2,cube2, 1)

mesh = Mesh(geom.GenerateMesh(maxh=0.5))
# two levels of hp-refinement
mesh.RefineHP(2, factor=0.2)
Draw (mesh)

[2]:

SetHeapSize(100*1000*1000)

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

a = BilinearForm(fes)
a += curl(u)*curl(v)*dx

m = BilinearForm(fes)
m += u*v*dx

apre = BilinearForm(fes)
apre += curl(u)*curl(v)*dx + u*v*dx
pre = Preconditioner(apre, "direct", inverse="sparsecholesky")

ndof = 32673

[3]:

with TaskManager():
    a.Assemble()
    m.Assemble()
    apre.Assemble()

    # build gradient matrix as sparse matrix (and corresponding scalar FESpace)
    gradmat, fesh1 = fes.CreateGradient()


    gradmattrans = gradmat.CreateTranspose() # transpose sparse matrix
    math1 = gradmattrans @ m.mat @ gradmat   # multiply matrices
    math1[0,0] += 1     # fix the 1-dim kernel
    invh1 = math1.Inverse(inverse="sparsecholesky")

    # build the Poisson projector with operator Algebra:
    proj = IdentityMatrix() - gradmat @ invh1 @ gradmattrans @ m.mat

    projpre = proj @ pre.mat

    evals, evecs = solvers.PINVIT(a.mat, m.mat, pre=projpre, num=12, maxit=20)

0 : [5.2016061125730335, 13.844259699441874, 17.005500584065551, 41.077594939418574, 45.081798966601532, 47.382379782770201, 55.57467853082067, 62.475697943371777, 64.591384900869528, 76.238450931850224, 79.626031367516831, 87.36061415556091]
1 : [3.2329398525264859, 6.0115389859946076, 6.0821761032449952, 11.978444313092586, 12.743245648604692, 14.61974605084036, 15.04122645086689, 15.751093688104062, 16.877934776950816, 19.947554329925993, 20.91671969361348, 27.16464982310524]
2 : [3.2206012255953302, 5.8852227257652849, 5.893091165670671, 10.864074925657167, 10.924514232090976, 12.478386683643222, 12.744800421829067, 13.318826601607629, 13.752489103021867, 14.468313566808503, 15.768248102933295, 21.985401691705434]
3 : [3.2204291185342138, 5.8807632547189925, 5.8815291591609391, 10.715038064773918, 10.737247816082862, 11.374966165617556, 12.406951796331434, 12.633641401736034, 13.464918883624032, 13.772913121190555, 14.134600556650746, 17.649282091658645]
4 : [3.2204266632333698, 5.8805734468934094, 5.8806162385781713, 10.696040140843809, 10.699249661014573, 10.84064030444468, 12.345413482616964, 12.430794585840941, 13.41439584600259, 13.576534937172838, 13.631798882477261, 15.723200905731977]
5 : [3.2204266284227998, 5.8805635854183489, 5.8805676974810819, 10.688929195612292, 10.694985449601942, 10.721014851900128, 12.327332147558534, 12.35956393708855, 13.389870872516276, 13.484966683212381, 13.504206393198396, 14.916836438020006]
6 : [3.2204266279135592, 5.8805624844709152, 5.8805660880883579, 10.687000163161317, 10.694503683687934, 10.699329514789305, 12.3215053269988, 12.333582933126785, 13.374452240516966, 13.44596553776535, 13.463053207706938, 14.561483655386315]
7 : [3.2204266279055918, 5.8805624124323117, 5.8805660313823838, 10.686490243997694, 10.694413942180597, 10.695348158731484, 12.31956519269389, 12.324166808390505, 13.366093268309157, 13.432633559811142, 13.442551962258472, 14.393713453859231]
8 : [3.2204266279054652, 5.8805624084633026, 5.8805660289385875, 10.686354488390773, 10.694369095799093, 10.694594086016126, 12.318889277090719, 12.32082247956447, 13.361968612858792, 13.427413018675686, 13.43248250917396, 14.310303293161137]
9 : [3.2204266279054616, 5.8805624082446526, 5.880566028825096, 10.686321218566048, 10.694315443033322, 10.694479510226241, 12.318626128185922, 12.319668706554094, 13.360002506099079, 13.425179877435587, 13.427668607320362, 14.266911996194994]
10 : [3.2204266279054661, 5.8805624082325849, 5.8805660288195174, 10.686313333573642, 10.694294267665299, 10.694463585703858, 12.318510965807711, 12.319286069075423, 13.359090138839644, 13.424186852717138, 13.425405185284284, 14.243767042862672]
11 : [3.2204266279054563, 5.8805624082319108, 5.8805660288192279, 10.686311479492751, 10.694288970892615, 10.694460478732861, 12.318462537992772, 12.319159601204097, 13.358673011764854, 13.423736071052394, 13.424350370463102, 14.231203104515847]
12 : [3.2204266279054599, 5.880562408231869, 5.8805660288192056, 10.686311043912443, 10.694287727048952, 10.694459788523194, 12.318443961331953, 12.319116711947887, 13.3584839239748, 13.423528469840837, 13.423862713389894, 14.224323429102764]
13 : [3.2204266279054528, 5.8805624082318531, 5.8805660288192367, 10.686310941500215, 10.694287436678152, 10.694459629564323, 12.318437189624097, 12.319101863497888, 13.358398423813732, 13.423430877036362, 13.423639469616253, 14.220533432096184]
14 : [3.2204266279054612, 5.8805624082317935, 5.8805660288192252, 10.686310917392385, 10.69428736870494, 10.694459592500394, 12.318434769095466, 12.319096670525905, 13.358359790717847, 13.423383971594559, 13.423538384586708, 14.218439020691296]
15 : [3.2204266279054563, 5.8805624082318699, 5.8805660288191044, 10.686310911713269, 10.69428735274775, 10.694459583813451, 12.318433909987046, 12.319094845906406, 13.358342319589045, 13.423361275518957, 13.42349284043919, 14.217278805889576]
16 : [3.2204266279054661, 5.8805624082318406, 5.8805660288191888, 10.68631091037309, 10.694287348989878, 10.694459581770802, 12.318433605712038, 12.319094203681226, 13.358334414202316, 13.423350419906832, 13.423472227099884, 14.216635518626958]
17 : [3.2204266279054599, 5.8805624082319614, 5.8805660288191905, 10.686310910056145, 10.694287348102669, 10.694459581289181, 12.31843349785901, 12.319093977204384, 13.358330832777716, 13.423345306176202, 13.423462815435522, 14.216278365188053]
18 : [3.2204266279054647, 5.8805624082317989, 5.8805660288190644, 10.686310909981552, 10.694287347891951, 10.694459581175437, 12.318433459617436, 12.319093897286328, 13.358329209396119, 13.423342922538481, 13.423458494844661, 14.216079973749826]
19 : [3.2204266279054616, 5.8805624082321017, 5.8805660288191932, 10.68631090996281, 10.694287347841883, 10.694459581148159, 12.318433445805207, 12.319093868942112, 13.358328472812442, 13.423341817536638, 13.423456471372022, 14.215969684583319]

[4]:

print ("Eigenvalues")
for lam in evals:
    print (lam)

Eigenvalues
3.22042662791
5.88056240823
5.88056602882
10.68631091
10.6942873478
10.6944595811
12.3184334458
12.3190938689
13.3583284728
13.4233418175
13.4234564714
14.2159696846

[5]:

Draw (mesh)
gfu = GridFunction(fes, multidim=len(evecs))
for i in range(len(evecs)):
    gfu.vecs[i].data = evecs[i]
Draw (Norm(gfu), mesh, "u", sd=4)

[ ]: