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.805905055837108, 24.468823847828265, 38.774015756493476, 42.35287378417515, 43.49334858885579, 47.10794241119643, 56.298906362897306, 62.6648001242203, 75.82039370073063, 84.8202648303027, 97.80562849293072, 102.9029734488957]
1 : [3.2417728818202227, 6.313486273643343, 6.690873834649835, 11.622157175424576, 12.696460506448272, 12.808301462681843, 14.000625794561847, 15.495738214552292, 16.619777102745132, 19.815082665085946, 25.436871664146352, 27.4188281296188]
2 : [3.22078295485238, 5.902332257216972, 5.922499099106113, 10.785126210418685, 11.151236588237781, 11.334867763447237, 12.639999246295455, 13.044739404306274, 13.91306413140521, 15.338626837384037, 16.687809468719717, 19.79541845354934]
3 : [3.220396287597829, 5.881904392753576, 5.88294628824511, 10.713653224473132, 10.779985246800937, 10.877888090957594, 12.393403744196613, 12.735451537504183, 13.573511364447887, 13.988245377617416, 14.754681838730999, 16.183227616455255]
4 : [3.220391243379386, 5.8806290976364615, 5.880719267968672, 10.699118393242195, 10.70529613638668, 10.742050645350153, 12.336771542895734, 12.561645809998971, 13.469599961451571, 13.65753015262637, 14.04278091494763, 14.888856112572414]
5 : [3.2203911736236877, 5.880562847565588, 5.8805768513821075, 10.690713905329526, 10.695604350419467, 10.70572500429955, 12.323666703083108, 12.437569669482215, 13.431899646211019, 13.530201465494445, 13.685509243326464, 14.46406787065888]
6 : [3.220391172626108, 5.880558617392532, 5.8805692937743155, 10.687480127942184, 10.694747947561899, 10.696923319490493, 12.320368058418167, 12.36785182214898, 13.41406673084819, 13.470072042696426, 13.51966013084148, 14.322982512597909]
7 : [3.220391172611192, 5.880558377452057, 5.880568896124848, 10.68665183621869, 10.694534336847529, 10.694995500286334, 12.319455842180112, 12.33734367863194, 13.396282426385456, 13.441512726150451, 13.45687177932954, 14.267582636812596]
8 : [3.2203911726109635, 5.880558365186419, 5.880568874510877, 10.686450487638151, 10.694428454702381, 10.69462963511683, 12.319183450443546, 12.325642518142955, 13.37847872609985, 13.430442030012394, 13.43703136524899, 14.242514888142358]
9 : [3.220391172610959, 5.880558364557214, 5.880568873345908, 10.686404048484366, 10.694374205903294, 10.694578467308627, 12.319097086040655, 12.321423722113297, 13.367973135621906, 13.426618705578422, 13.429547775392873, 14.230130272165704]
10 : [3.2203911726109573, 5.880558364524582, 5.8805688732833445, 10.68639351536791, 10.694359635961309, 10.694569234895695, 12.319067727169017, 12.31994410282556, 13.362897818195279, 13.425131823529574, 13.426328499895915, 14.223700897159695]
11 : [3.2203911726109578, 5.880558364522865, 5.8805688732799695, 10.686391132542687, 10.694356230377563, 10.69456723012046, 12.319055636336005, 12.319432813648195, 13.360572922331345, 13.42449276541656, 13.424900682687586, 14.220272240411301]
12 : [3.2203911726109564, 5.880558364522774, 5.880568873279811, 10.686390592815352, 10.694355450175841, 10.694566774297858, 12.31904851806154, 12.319259128475064, 13.35952659221752, 13.424148384664337, 13.42431811991552, 14.218417288813578]
13 : [3.220391172610956, 5.880558364522728, 5.880568873279767, 10.686390470298052, 10.694355271672688, 10.694566669392875, 12.31904406144969, 12.31920118743312, 13.359058586239604, 13.423900284901556, 13.42414671689998, 14.217405988689185]
14 : [3.220391172610956, 5.880558364522754, 5.880568873279798, 10.686390442436624, 10.6943552307503, 10.694566645146313, 12.319041880035915, 12.319181834502153, 13.358849766074865, 13.423773678304524, 13.424082658957833, 14.216852283652837]
15 : [3.220391172610958, 5.880558364522771, 5.8805688732797945, 10.686390436088686, 10.694355221345676, 10.694566639529022, 12.319041004200605, 12.31917527271754, 13.358756617935617, 13.423714703113108, 13.424054756599878, 14.216548288976268]
16 : [3.2203911726109533, 5.880558364522711, 5.880568873279868, 10.686390434639163, 10.694355219179242, 10.694566638225504, 12.319040683047785, 12.319173022748805, 13.358715069610565, 13.423687577200612, 13.4240422443787, 14.216381088828268]
17 : [3.2203911726109573, 5.880558364522778, 5.880568873279815, 10.686390434307057, 10.694355218673264, 10.69456663792152, 12.319040569212866, 12.319172247087177, 13.358696373408492, 13.423674709743809, 13.424036538893391, 14.216286784338408]
18 : [3.2203911726109586, 5.880558364522926, 5.880568873279996, 10.686390434231479, 10.69435521856259, 10.694566637851592, 12.31904052906862, 12.319171977603483, 13.35868806072468, 13.423669155819402, 13.424033954386909, 14.216236409099205]
19 : [3.220391172610959, 5.880558364522759, 5.880568873279848, 10.686390434214422, 10.694355218537114, 10.694566637835871, 12.319040515178672, 12.319171885027034, 13.35868438846117, 13.423666667170876, 13.424032785419149, 14.216209036488252]

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

Eigenvalues
3.220391172610959
5.880558364522759
5.880568873279848
10.686390434214422
10.694355218537114
10.694566637835871
12.319040515178672
12.319171885027034
13.35868438846117
13.423666667170876
13.424032785419149
14.216209036488252

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

BaseWebGuiScene