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

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 : [7.822724748852332, 14.737952312200752, 24.325572835958095, 43.873650237143664, 49.56750803943967, 58.372016641528425, 67.64803658027547, 68.97661297835458, 73.78261210478703, 79.83141581479963, 87.65300356490404, 105.64272847610339]
1 : [3.2513016736898424, 5.99631021749792, 6.188064306160371, 11.976761101169236, 12.935054636676897, 13.836081393549488, 14.78119062167469, 15.893304235733844, 16.877821696613054, 19.32978953546393, 22.0172983366227, 26.727681506451773]
2 : [3.2208288910199587, 5.88411930117543, 5.894154922348247, 10.822715656065892, 11.029262447429367, 11.214705799794377, 12.651256404638694, 13.274832975066145, 14.048275535183338, 15.429389635680987, 15.811202524885443, 19.678189448377587]
3 : [3.2204108525000987, 5.880694109654364, 5.881359849343331, 10.718849709810458, 10.74697927675992, 10.787327571735648, 12.431204349902634, 12.82071457166385, 13.606266334205529, 13.92036430506761, 14.825216978754137, 17.35118008575142]
4 : [3.2204040819007167, 5.880575104167088, 5.880639759636684, 10.696542628488569, 10.699546677410149, 10.714875508440484, 12.360986931854322, 12.61630529980066, 13.468257518714728, 13.560721761031786, 14.419577601897698, 16.06002542702825]
5 : [3.220403962403185, 5.880567302619671, 5.880595433322504, 10.688331320118474, 10.695423577215374, 10.699419894574133, 12.334182372223044, 12.494368049164887, 13.409960562569719, 13.471804862878683, 14.186353356517987, 15.090993869915152]
6 : [3.2204039602541683, 5.880566393330402, 5.880592847433293, 10.686806844324494, 10.694618413495265, 10.695665285143146, 12.32439682539465, 12.413229702981633, 13.386824843718063, 13.443465517555499, 13.944078259081886, 14.544355752455902]
7 : [3.2204039602156924, 5.880566329769076, 5.880592684844842, 10.68645018127875, 10.694444044200978, 10.694756182957052, 12.32097420294457, 12.363170790091647, 13.377142793454585, 13.432389415336178, 13.701733239092674, 14.35182231363885]
8 : [3.2204039602150196, 5.880566325906187, 5.880592675055868, 10.686360638126342, 10.694384961954553, 10.694558943349046, 12.319808916629787, 12.33744371947085, 13.371826563401813, 13.427682104366838, 13.554112148367839, 14.282161383858401]
9 : [3.2204039602150103, 5.880566325684012, 5.880592674493732, 10.6863389337744, 10.694363383084008, 10.694520734643781, 12.31941422782044, 12.326300093681215, 13.367625799321898, 13.425606609297551, 13.481894579421812, 14.250743534695717]
10 : [3.22040396021501, 5.880566325671631, 5.8805926744623145, 10.686333818088835, 10.69435767897066, 10.694512599578347, 12.31927942059609, 12.321955380096401, 13.364048824573725, 13.42467523985514, 13.449395794131055, 14.234894167033984]
11 : [3.2204039602150103, 5.88056632567097, 5.88059267446058, 10.686332628070083, 10.694356337787012, 10.69451074491658, 12.319229929301423, 12.32035708645573, 13.3615264406589, 13.424252904118, 13.435235315318696, 14.226491993376355]
12 : [3.2204039602150063, 5.880566325670899, 5.880592674460485, 10.686332352080674, 10.694356029593314, 10.694510313808605, 12.319207035876397, 12.319789248185968, 13.36007151862034, 13.424059303411848, 13.429043968933314, 14.221940143018204]
13 : [3.2204039602150085, 5.880566325670902, 5.880592674460491, 10.686332288055333, 10.694355958800775, 10.694510213149828, 12.319194005092799, 12.319593028873317, 13.359330689705587, 13.423967809955688, 13.4262965129094, 14.219448025314566]
14 : [3.2204039602150094, 5.880566325670918, 5.880592674460414, 10.686332273164197, 10.69435594246097, 10.694510189592684, 12.319187368356747, 12.319525844403145, 13.358976260504965, 13.423919876939744, 13.425065185713676, 14.218075849834324]
15 : [3.220403960215009, 5.880566325670937, 5.880592674460452, 10.686332269699601, 10.694355938675654, 10.694510184080164, 12.319184588153666, 12.31950267541184, 13.358811994485725, 13.423887922117125, 13.424517268455206, 14.217318796828371]
16 : [3.220403960215005, 5.880566325671022, 5.880592674460492, 10.686332268893207, 10.694355937796587, 10.694510182788074, 12.319183548000634, 12.319494606354935, 13.358737029441627, 13.423861639570893, 13.424279953835825, 14.216901349478363]
17 : [3.220403960215007, 5.880566325671091, 5.880592674460489, 10.686332268704557, 10.694355937591139, 10.69451018248471, 12.319183175458514, 12.319491770769858, 13.35870289407413, 13.423841412944999, 13.424179718855044, 14.216670061711804]
18 : [3.22040396021501, 5.880566325671062, 5.880592674460572, 10.686332268660726, 10.694355937543454, 10.694510182413211, 12.319183044601967, 12.319490775314593, 13.358687455691042, 13.423828776482756, 13.424137396114368, 14.216539458871186]
19 : [3.2204039602150085, 5.880566325670964, 5.8805926744606385, 10.686332268650215, 10.694355937532052, 10.694510182396609, 12.319182998810094, 12.319490425497243, 13.358680482525028, 13.423822025635875, 13.424119032078682, 14.216468789610422]

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

Eigenvalues
3.2204039602150085
5.880566325670964
5.8805926744606385
10.686332268650215
10.694355937532052
10.694510182396609
12.319182998810094
12.319490425497243
13.358680482525028
13.423822025635875
13.424119032078682
14.216468789610422

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

BaseWebGuiScene