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.

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.csg 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)

BaseWebGuiScene

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

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.805905053742296, 24.46882384880346, 38.77401575622128, 42.35287378023363, 43.49334859587597, 47.10794240589327, 56.2989063582672, 62.66480012031737, 75.82039370526336, 84.82026483343141, 97.80562848348872, 102.90297345062555]
1 : [3.241772881796691, 6.313486273615755, 6.690873835096278, 11.622157174805526, 12.696460506767005, 12.808301461508957, 14.000625794710533, 15.495738213673516, 16.619777103943232, 19.815082667349895, 25.43687165842567, 27.418828132942252]
2 : [3.220782954851952, 5.902332257216597, 5.9224990991050595, 10.785126210388837, 11.151236587759497, 11.33486776366769, 12.639999246559741, 13.044739404675274, 13.913064131513064, 15.338626836488064, 16.687809465641962, 19.79541845545134]
3 : [3.220396287597826, 5.881904392753476, 5.882946288245191, 10.71365322447372, 10.779985246679761, 10.877888091124243, 12.393403744248042, 12.735451537671628, 13.57351136442006, 13.988245377580853, 14.754681838111052, 16.183227616363133]
4 : [3.220391243379389, 5.880629097636515, 5.880719267968642, 10.699118393243943, 10.705296136361993, 10.742050645389547, 12.336771542902566, 12.561645810020794, 13.469599961388072, 13.657530153002158, 14.042780914131988, 14.8888561119473]
5 : [3.220391173623684, 5.8805628475655665, 5.880576851382124, 10.690713905322411, 10.695604350419618, 10.705725004305053, 12.323666703084221, 12.437569669424303, 13.431899646156639, 13.530201465766675, 13.685509242681695, 14.464067870358429]
6 : [3.2203911726261047, 5.880558617392533, 5.880569293774314, 10.687480127939438, 10.694747947562007, 10.69692331949155, 12.320368058418424, 12.367851822101247, 13.414066730784898, 13.470072042821343, 13.519660130550278, 14.322982512452434]
7 : [3.220391172611192, 5.880558377452078, 5.880568896124846, 10.68665183621796, 10.6945343368476, 10.694995500286533, 12.319455842180245, 12.337343678608772, 13.39628242630039, 13.441512726175915, 13.456871779278758, 14.267582636735336]
8 : [3.220391172610964, 5.880558365186417, 5.880568874510886, 10.686450487638064, 10.69442845470242, 10.694629635116838, 12.31918345044373, 12.325642518134856, 13.378478726043937, 13.430442030018323, 13.437031365254247, 14.242514888092073]
9 : [3.2203911726109578, 5.880558364557214, 5.880568873345901, 10.686404048484569, 10.694374205903408, 10.694578467308702, 12.319097086040903, 12.32142372211076, 13.367973135553889, 13.426618705556573, 13.429547775419808, 14.230130272225521]
10 : [3.22039117261096, 5.880558364524592, 5.880568873283353, 10.68639351536815, 10.69435963596133, 10.694569234895724, 12.319067727169555, 12.319944102832968, 13.362897818173586, 13.425131823496923, 13.426328499873453, 14.223700897165504]
11 : [3.220391172610957, 5.880558364522866, 5.880568873279996, 10.686391132542381, 10.694356230378819, 10.694567230120459, 12.319055636335186, 12.319432813193613, 13.360572918230275, 13.424492765095108, 13.424900682996558, 14.220272243321165]
12 : [3.2203911726109538, 5.88055836452276, 5.880568873279794, 10.68639059281571, 10.694355450179193, 10.694566774299355, 12.319048518065939, 12.319259128630156, 13.35952659044539, 13.424148387446845, 13.424318121778876, 14.21841731305947]
13 : [3.220391172610961, 5.880558364522761, 5.8805688732797305, 10.686390470304723, 10.694355271673093, 10.694566669394092, 12.319044061624071, 12.319201191180031, 13.359058637130893, 13.423900288392474, 13.424146721942623, 14.21740600117799]
14 : [3.2203911726109533, 5.88055836452275, 5.8805688732797625, 10.686390442439183, 10.694355230747771, 10.694566645146672, 12.319041880261166, 12.319181836143786, 13.358849797641518, 13.423773663152893, 13.424082668257743, 14.216852247642741]
15 : [3.2203911726109564, 5.880558364522773, 5.880568873279774, 10.686390436088491, 10.69435522134342, 10.694566639528471, 12.319041004279555, 12.319175272699288, 13.358756633953229, 13.423714676674196, 13.424054766284494, 14.216548163113615]
16 : [3.2203911726109546, 5.880558364522931, 5.880568873279845, 10.686390434639046, 10.694355219178993, 10.694566638225307, 12.31904068310516, 12.31917302266845, 13.358715072143228, 13.423687569470218, 13.42404224785341, 14.216381039643826]
17 : [3.2203911726109555, 5.880558364522976, 5.880568873279786, 10.686390434307466, 10.694355218678735, 10.69456663792146, 12.319040569003263, 12.319172247447158, 13.358696531480051, 13.423675020874965, 13.424036572984779, 14.216288167735867]
18 : [3.220391172610957, 5.880558364522868, 5.880568873279536, 10.68639043423188, 10.694355218564061, 10.694566637851771, 12.319040529225875, 12.319171979481704, 13.35868825706329, 13.423669367368927, 13.424034000866124, 14.216237669229502]
19 : [3.2203911726109538, 5.880558364522408, 5.880568873279774, 10.686390434214333, 10.69435521853809, 10.694566637835674, 12.319040515324058, 12.319171886407961, 13.358684546190746, 13.423666807471676, 13.424032825902136, 14.216210070726403]

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

Eigenvalues
3.2203911726109538
5.880558364522408
5.880568873279774
10.686390434214333
10.69435521853809
10.694566637835674
12.319040515324058
12.319171886407961
13.358684546190746
13.423666807471676
13.424032825902136
14.216210070726403

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

BaseWebGuiScene