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 : [7.822091762752782, 24.10234221989844, 25.410210861561673, 31.591610086976, 41.04286265819465, 54.05306393208564, 55.57754470599475, 62.48922689789102, 73.5973762954086, 84.2113346521632, 88.22650435274521, 101.57701296671485]
1 : [3.2501484717633202, 6.125216794638108, 6.182311947956383, 11.524201703408531, 12.239391467664309, 13.362041719366227, 14.80507377723048, 15.24416731939922, 18.536522772033788, 18.83417233281189, 20.65939144569115, 23.476093201078196]
2 : [3.220774479247371, 5.88838121977756, 5.890395451030877, 10.773167306077193, 10.984625590013435, 11.204618750302684, 12.638158094680959, 13.123265344516328, 13.800547756008013, 14.995199238850853, 15.571911840846706, 17.48216978194901]
3 : [3.22039616265098, 5.880863278506782, 5.881043600703523, 10.703264199589414, 10.75365813932032, 10.83927913591117, 12.436787456620841, 12.591167821143655, 13.46744598828303, 13.973218319716791, 14.279890742556029, 16.93037982176147]
4 : [3.220391108732277, 5.880574208109167, 5.8805948348531905, 10.695015549372977, 10.703628085291545, 10.74602559651452, 12.365315433063007, 12.429770002810017, 13.41258343570958, 13.6290657402356, 14.026615076515677, 16.6644795627324]
5 : [3.2203910320881666, 5.880559916842052, 5.880570294367758, 10.692557890170661, 10.695326232548588, 10.712028436952988, 12.334863812929656, 12.379428178899833, 13.399760629323579, 13.507187416164152, 13.909165088376836, 16.3910781910024]
6 : [3.2203910307938863, 5.880558531346238, 5.880569271711317, 10.690148741453712, 10.694658496397235, 10.69952889464496, 12.324262124286859, 12.355591119237719, 13.394789054615616, 13.458331319354349, 13.826480883808635, 16.001910514120944]
7 : [3.220391030769952, 5.880558425991063, 5.8805692134166225, 10.688201871206966, 10.694580575766366, 10.695760990661684, 12.32078268723305, 12.341027972960482, 13.391306508177061, 13.438346613659048, 13.74937816928353, 15.513779089452942]
8 : [3.22039103076949, 5.880558418312408, 5.880569209656792, 10.687056211427002, 10.69456559727026, 10.694758763422122, 12.319633665711592, 12.331359871823649, 13.38827850940228, 13.430127743791495, 13.666503234404585, 15.038678640413853]
9 : [3.2203910307694783, 5.880558417786905, 5.880569209414908, 10.686597154749318, 10.694466636428984, 10.69457030327245, 12.319248022950363, 12.325288286793977, 13.385179207444821, 13.426655484756694, 13.58227792755703, 14.688877890216428]
10 : [3.220391030769484, 5.880558417753148, 5.88056920939987, 10.686447460300359, 10.694386377698597, 10.694567747898832, 12.319117074778093, 12.321905868095778, 13.381103350555872, 13.425137658139734, 13.511736466847163, 14.47586632489388]
11 : [3.2203910307694814, 5.880558417751133, 5.880569209398983, 10.686405049915098, 10.694364165698099, 10.694567368329603, 12.319071997361842, 12.320286361644039, 13.375465543103271, 13.424457985230685, 13.46630298850759, 14.358019821360292]
12 : [3.2203910307694823, 5.880558417751016, 5.880569209398922, 10.686394009872949, 10.694358423739923, 10.694567284811102, 12.319055986446434, 12.31959934973428, 13.369259852382497, 13.42414883510403, 13.442918931520389, 14.294004922649751]
13 : [3.2203910307694827, 5.880558417751002, 5.880569209398904, 10.686391282600455, 10.694357004646372, 10.694567265287642, 12.319049668196527, 12.319331987259377, 13.364427360687706, 13.4240058895651, 13.43233972089806, 14.259097727180253]
14 : [3.2203910307694796, 5.880558417750574, 5.880569209398955, 10.68639062698609, 10.69435666262042, 10.694567260680818, 12.319046756108182, 12.319233286508437, 13.361535979762573, 13.423936594268222, 13.427676199039626, 14.239940426412867]
15 : [3.2203910307694836, 5.880558417750923, 5.880569209398622, 10.686390471544922, 10.694356581247135, 10.694567259593004, 12.319045391861541, 12.319197855476183, 13.360036337379759, 13.423896453287828, 13.42559461247667, 14.229353486832718]
16 : [3.220391030769484, 5.880558417750749, 5.880569209398984, 10.686390434937259, 10.694356562012922, 10.694567259336056, 12.319044823596093, 12.319185273034654, 13.359311159056384, 13.423859997712796, 13.42466608665159, 14.223487935637632]
17 : [3.2203910307694805, 5.8805584177509065, 5.8805692093987965, 10.686390426331512, 10.694356557475741, 10.694567259275157, 12.31904460740204, 12.31918080999599, 13.358970892818785, 13.423810563084567, 13.424271838956209, 14.220219724640984]
18 : [3.220391030769482, 5.880558417751453, 5.880569209399031, 10.686390424322056, 10.694356556409977, 10.694567259261047, 12.319044529854956, 12.319179233842886, 13.35881463759433, 13.423754226006643, 13.424124234623621, 14.218407635461812]
19 : [3.220391030769485, 5.880558417751461, 5.880569209399395, 10.686390423854304, 10.69435655616041, 10.694567259257834, 12.319044502503612, 12.319178677070601, 13.358743434874382, 13.423712912962907, 13.424071532563206, 14.217402365024705]

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

Eigenvalues
3.220391030769485
5.880558417751461
5.880569209399395
10.686390423854304
10.69435655616041
10.694567259257834
12.319044502503612
12.319178677070601
13.358743434874382
13.423712912962907
13.424071532563206
14.217402365024705

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

BaseWebGuiScene