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.805905055748913, 24.468823849161403, 38.77401575576791, 42.352873789188436, 43.49334858993122, 47.10794239899552, 56.29890637043648, 62.66480012430691, 75.820393699844, 84.82026483114946, 97.80562847595765, 102.90297344574111]
1 : [3.2417728817978544, 6.3134862735176, 6.690873833652178, 11.622157175845732, 12.696460506609334, 12.80830146157237, 14.000625795191999, 15.495738212696908, 16.619777103433734, 19.815082667028964, 25.436871654149943, 27.418828132131324]
2 : [3.2207829548516704, 5.902332257207094, 5.922499099010635, 10.785126210447848, 11.151236587818497, 11.334867763841903, 12.639999246457842, 13.044739404132173, 13.913064132073593, 15.338626835734877, 16.687809464104333, 19.795418461213035]
3 : [3.2203962875978323, 5.881904392752279, 5.882946288245299, 10.7136532244744, 10.779985246701726, 10.877888091578049, 12.393403744199603, 12.735451537483685, 13.573511364752843, 13.988245377767287, 14.754681838826183, 16.183227624184454]
4 : [3.220391243379392, 5.880629097636493, 5.88071926796877, 10.699118393243028, 10.705296136371787, 10.742050645591291, 12.336771542895928, 12.561645810016826, 13.469599961537398, 13.65753015357598, 14.042780914431615, 14.888856115663371]
5 : [3.2203911736236805, 5.8805628475656055, 5.880576851382126, 10.69071390532471, 10.69560435041984, 10.705725004364446, 12.323666703084568, 12.437569669479215, 13.431899646223295, 13.530201466193954, 13.685509242912904, 14.464067871807206]
6 : [3.2203911726261083, 5.8805586173925475, 5.88056929377431, 10.687480127940303, 10.694747947562297, 10.696923319505977, 12.320368058419113, 12.367851822139, 13.414066730820105, 13.470072043026716, 13.519660130705073, 14.322982513055205]
7 : [3.2203911726111927, 5.880558377452089, 5.880568896124844, 10.686651836218159, 10.694534336848001, 10.694995500289476, 12.31945584218045, 12.337343678624846, 13.396282426324257, 13.441512726234588, 13.456871779388731, 14.26758263701748]
8 : [3.2203911726109604, 5.880558365186413, 5.880568874510891, 10.686450487638055, 10.694428454702962, 10.69462963511718, 12.319183450443589, 12.32564251814109, 13.378478726058322, 13.430442030038733, 13.437031365315123, 14.242514888239079]
9 : [3.220391172610955, 5.880558364557221, 5.880568873345898, 10.686404048484599, 10.694374205903605, 10.694578467308723, 12.319097086041246, 12.321423722113245, 13.36797313555924, 13.426618705575455, 13.429547775446753, 14.230130272302448]
10 : [3.2203911726109604, 5.880558364524588, 5.880568873283353, 10.686393515367776, 10.694359635961497, 10.694569234895768, 12.319067727168948, 12.319944102815583, 13.362897818121736, 13.425131823546149, 13.426328499963796, 14.223700897351442]
11 : [3.2203911726109546, 5.88055836452286, 5.880568873279998, 10.686391132542319, 10.694356230378967, 10.69456723012046, 12.3190556363306, 12.319432813227476, 13.3605729188584, 13.424492765296359, 13.424900683214865, 14.220272244022006]
12 : [3.220391172610956, 5.880558364522785, 5.88056887327982, 10.686390592814202, 10.69435545017608, 10.694566774298817, 12.31904851806399, 12.319259128309408, 13.359526590080527, 13.424148384883882, 13.424318121319043, 14.218417300696297]
13 : [3.2203911726109578, 5.880558364522772, 5.880568873279801, 10.68639047030473, 10.69435527167177, 10.69456666939346, 12.31904406169159, 12.319201191705105, 13.359058651348128, 13.423900286192728, 13.424146725274825, 14.217405979354426]
14 : [3.220391172610956, 5.88055836452276, 5.880568873279769, 10.686390442438432, 10.694355230748885, 10.6945666451465, 12.319041880273534, 12.319181835873165, 13.358849805580235, 13.423773675242508, 13.424082672077764, 14.216852259092256]
15 : [3.2203911726109573, 5.8805583645227815, 5.880568873279887, 10.686390436088212, 10.694355221344352, 10.694566639529034, 12.319041004358377, 12.31917527273193, 13.358756638692523, 13.423714696403655, 13.424054769871496, 14.21654823873062]
16 : [3.220391172610955, 5.8805583645225, 5.880568873279934, 10.686390434638827, 10.694355219179059, 10.694566638225549, 12.319040683125609, 12.319173022656306, 13.35871507881084, 13.423687575176277, 13.424042252065536, 14.21638104932492]
17 : [3.2203911726109573, 5.880558364522712, 5.880568873279866, 10.686390434307548, 10.694355218679771, 10.694566637922586, 12.319040569207797, 12.319172247053595, 13.358696529105194, 13.423675145426806, 13.424036579459278, 14.216289032042347]
18 : [3.220391172610956, 5.880558364523185, 5.880568873279789, 10.68639043423133, 10.694355218564615, 10.694566637850237, 12.319040526573678, 12.319171968173084, 13.358687669695344, 13.423669438289293, 13.424033689276486, 14.216237454317593]
19 : [3.2203911726109586, 5.880558364522438, 5.880568873279899, 10.686390434213552, 10.694355218537652, 10.694566637834702, 12.319040513400934, 12.31917187815353, 13.358683883752096, 13.42366652570552, 13.424032535311522, 14.216206770984144]

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

Eigenvalues
3.2203911726109586
5.880558364522438
5.880568873279899
10.686390434213552
10.694355218537652
10.694566637834702
12.319040513400934
12.31917187815353
13.358683883752096
13.42366652570552
13.424032535311522
14.216206770984144

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

BaseWebGuiScene