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.805905055547128, 24.468823844823696, 38.77401575634833, 42.35287378622451, 43.49334859549769, 47.10794239307522, 56.29890637091626, 62.66480012605601, 75.8203936960512, 84.82026483280802, 97.8056284961777, 102.90297344974414]
1 : [3.2417728818003475, 6.313486273375349, 6.690873833909696, 11.622157175694639, 12.696460506816, 12.808301460883845, 14.00062579488828, 15.495738214352345, 16.61977710354413, 19.81508266815309, 25.436871658038854, 27.41882812939409]
2 : [3.2207829548517104, 5.9023322572037005, 5.922499099016846, 10.785126210465748, 11.151236587704604, 11.334867763568266, 12.63999924637448, 13.044739404336976, 13.913064131884502, 15.338626834906673, 16.687809467575008, 19.79541845534263]
3 : [3.2203962875978274, 5.881904392752324, 5.882946288243623, 10.713653224483233, 10.779985246654451, 10.877888091231242, 12.393403744203109, 12.735451537616466, 13.573511364680103, 13.988245376916767, 14.75468183972399, 16.183227618863278]
4 : [3.220391243379386, 5.880629097636502, 5.8807192679686295, 10.699118393245348, 10.705296136353988, 10.742050645441234, 12.336771542901735, 12.56164581006459, 13.469599961511356, 13.657530152695127, 14.042780914894967, 14.888856113198354]
5 : [3.2203911736236845, 5.880562847565589, 5.880576851382139, 10.690713905319484, 10.695604350420826, 10.705725004320021, 12.323666703086475, 12.43756966947298, 13.431899646201163, 13.53020146565297, 13.685509242998265, 14.464067870903682]
6 : [3.22039117262611, 5.880558617392542, 5.880569293774316, 10.687480127938954, 10.694747947562334, 10.696923319495044, 12.320368058419582, 12.367851822126722, 13.414066730799181, 13.470072042796401, 13.519660130664946, 14.322982512706826]
7 : [3.22039117261119, 5.8805583774520915, 5.880568896124856, 10.686651836217804, 10.694534336847713, 10.694995500287202, 12.319455842180655, 12.33734367861679, 13.396282426313077, 13.441512726191004, 13.456871779293603, 14.267582636860032]
8 : [3.220391172610966, 5.8805583651864195, 5.880568874510879, 10.686450487637956, 10.694428454702528, 10.694629635116888, 12.319183450443717, 12.32564251813755, 13.378478726054047, 13.430442030030731, 13.437031365247227, 14.242514888149717]
9 : [3.220391172610956, 5.880558364557233, 5.880568873345878, 10.686404048484563, 10.694374205903365, 10.694578467308661, 12.319097086041312, 12.321423722111781, 13.367973135563101, 13.426618705571466, 13.429547775397031, 14.230130272204056]
10 : [3.220391172610957, 5.88055836452457, 5.880568873283339, 10.686393515367946, 10.694359635961574, 10.694569234895768, 12.319067727168983, 12.319944102828925, 13.362897818148442, 13.425131823538628, 13.426328499972568, 14.223700897403845]
11 : [3.220391172610956, 5.8805583645228605, 5.880568873279989, 10.686391132542505, 10.694356230378695, 10.6945672301205, 12.319055636329914, 12.319432813376828, 13.360572920161788, 13.424492765384958, 13.424900683186005, 14.220272243521721]
12 : [3.2203911726109538, 5.880558364522768, 5.880568873279795, 10.686390592815094, 10.694355450177921, 10.694566774299226, 12.319048518067, 12.319259128632774, 13.359526591759938, 13.424148386832172, 13.424318121724209, 14.218417308488903]
13 : [3.220391172610961, 5.880558364522759, 5.880568873279825, 10.686390470303643, 10.694355271671084, 10.694566669393888, 12.319044061645336, 12.319201191035384, 13.359058630019428, 13.423900276787485, 13.424146722043824, 14.217405980504447]
14 : [3.2203911726109538, 5.880558364522787, 5.880568873279829, 10.686390442437276, 10.694355230749045, 10.694566645146947, 12.319041880273165, 12.319181835520718, 13.358849789716556, 13.423773673568133, 13.424082669581866, 14.216852277820864]
15 : [3.2203911726109578, 5.8805583645228126, 5.8805688732798265, 10.686390436087752, 10.694355221344901, 10.694566639529228, 12.319041004338837, 12.31917527236045, 13.358756624387704, 13.423714701296497, 13.424054767599314, 14.216548277598678]
16 : [3.2203911726109546, 5.8805583645227, 5.880568873279774, 10.686390434638653, 10.694355219179329, 10.694566638225437, 12.31904068312039, 12.319173022553723, 13.358715071019919, 13.423687579229194, 13.424042250636429, 14.216381086055636]
17 : [3.220391172610955, 5.880558364522669, 5.880568873279808, 10.686390434307187, 10.694355218679602, 10.694566637922353, 12.31904056921327, 12.319172247049194, 13.358696526880028, 13.423675145492583, 13.424036580184268, 14.21628904246726]
18 : [3.2203911726109578, 5.880558364523925, 5.880568873280034, 10.68639043423121, 10.69435521856367, 10.694566637850857, 12.319040528706322, 12.319171978840185, 13.358688239303412, 13.423669189854108, 13.424033967820241, 14.216236021046]
19 : [3.2203911726109546, 5.880558364522759, 5.880568873279573, 10.68639043421362, 10.69435521853757, 10.694566637835115, 12.319040514938724, 12.319171884310782, 13.358684394785925, 13.42366658409081, 13.42403279461013, 14.216207842783662]

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

Eigenvalues
3.2203911726109546
5.880558364522759
5.880568873279573
10.68639043421362
10.69435521853757
10.694566637835115
12.319040514938724
12.319171884310782
13.358684394785925
13.42366658409081
13.42403279461013
14.216207842783662

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

BaseWebGuiScene