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.

[1]:

import netgen.gui
from netgen.csg import *
from ngsolve 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)

[2]:

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

[3]:

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 : [4.6898535365197942, 15.984971330974986, 26.690290636810534, 31.737251842694576, 44.586227579189718, 47.708983785928467, 55.742041044003749, 62.521110476936123, 67.599839944728046, 75.359563420170062, 78.919166039157844, 88.194869592938559]
1 : [3.2296338274377372, 6.0158608786252827, 6.2680815009821496, 11.689110066856381, 12.360804427284862, 12.932028451948286, 14.486646571980366, 15.609488331606634, 17.458210915017755, 20.385687810499025, 23.304539890035073, 25.854753252395202]
2 : [3.2205657903259, 5.8851174014450125, 5.9002857899351007, 10.8364946949124, 10.93496076683622, 11.199112460764038, 12.752307384529812, 13.569235221875894, 13.965885071668124, 15.226534500911491, 16.1355449622889, 18.493179199186276]
3 : [3.2204283484873146, 5.8807255597680754, 5.8820974930061993, 10.70684829490243, 10.750525096695156, 10.844877377326977, 12.530143551158126, 12.769660255225535, 13.558194487273939, 13.923612374412775, 14.107235924991169, 17.286621830307613]
4 : [3.2204266502609187, 5.8805699933920446, 5.8806918058498008, 10.690700691457787, 10.710410069468784, 10.730264868056372, 12.421311639952229, 12.476707167787854, 13.458299245952039, 13.579309296295488, 13.796715380271776, 16.89948722179772]
5 : [3.2204266282250504, 5.8805641038893652, 5.8805753056788141, 10.687760483703844, 10.698399745997259, 10.704498011244089, 12.350919967141966, 12.41385296086446, 13.431251564938508, 13.476635332302253, 13.673307422770986, 16.729813673368987]
6 : [3.2204266279106135, 5.8805628043581271, 5.8805665419491246, 10.68671658728241, 10.695296647026627, 10.697987862183689, 12.329015162317212, 12.379373153614246, 13.422683189013064, 13.439026404075735, 13.604118396306793, 16.588568133237697]
7 : [3.2204266279055513, 5.8805624448273814, 5.8805660665029116, 10.686408916780026, 10.694636221981469, 10.695708924167493, 12.322308770856997, 12.3567243271204, 13.415900772260496, 13.426908643368124, 13.555959703411583, 16.389924465965255]
8 : [3.220426627905459, 5.8805624111626642, 5.8805660317641024, 10.686333945388473, 10.694487773195668, 10.694830251274336, 12.320154544615766, 12.341686429448382, 13.409318747796499, 13.424435441808878, 13.518281045767571, 16.088284267351668]
9 : [3.220426627905463, 5.8805624084528993, 5.8805660290403265, 10.686316351677918, 10.694429262254946, 10.694515125482047, 12.319437105517286, 12.331836840426964, 13.404143238675511, 13.423773263527941, 13.487232388482536, 15.689538863138486]
10 : [3.2204266279054661, 5.8805624082474104, 5.8805660288347701, 10.686312205622077, 10.694343601524999, 10.694467262467477, 12.319191543797082, 12.325611770087189, 13.398462429678592, 13.42353640632582, 13.462017834809412, 15.263034709085515]
11 : [3.220426627905467, 5.8805624082328798, 5.8805660288202501, 10.686311219977741, 10.694304620449035, 10.694461393690876, 12.319103617372116, 12.321952864890667, 13.390480874068901, 13.423439579994177, 13.443872568534383, 14.894731027253442]
12 : [3.2204266279054545, 5.8805624082319508, 5.8805660288192971, 10.686310984139114, 10.694292149813254, 10.694460042191013, 12.319061363590437, 12.320025805380109, 13.380562242664016, 13.423397347572914, 13.433312993908721, 14.627478410147882]
13 : [3.2204266279054576, 5.8805624082316372, 5.8805660288197412, 10.686310927679598, 10.694288595190113, 10.69445969735647, 12.318946266323604, 12.319217325683464, 13.371589274741277, 13.42337820864282, 13.428099507538874, 14.455339512798584]
14 : [3.2204266279054572, 5.8805624082318415, 5.880566028819203, 10.686310914180849, 10.694287657741587, 10.694459609786977, 12.318660298741465, 12.319111716206221, 13.365417727665672, 13.423369110142549, 13.425637224909206, 14.352014754172494]
15 : [3.2204266279054519, 5.880562408231885, 5.8805660288192074, 10.686310910962705, 10.694287422778773, 10.694459588078271, 12.318518915404081, 12.319098978651647, 13.361859350308626, 13.423364199542133, 13.424473152533166, 14.292307488674345]
16 : [3.2204266279054585, 5.880562408231933, 5.8805660288193469, 10.686310910196749, 10.694287365695809, 10.694459582802258, 12.318464548668448, 12.319095559025822, 13.360014629844198, 13.423360579065138, 13.423923646437201, 14.258518762800632]
17 : [3.2204266279054572, 5.8805624082329242, 5.8805660288193939, 10.686310910014738, 10.694287352052752, 10.694459581535984, 12.318444586347592, 12.319094443931901, 13.359114215597977, 13.42335669509033, 13.423666864079326, 14.239577184027908]
18 : [3.2204266279054599, 5.8805624082315013, 5.8805660288195645, 10.686310909971562, 10.694287348823597, 10.694459581233877, 12.318437406338928, 12.319094060483383, 13.35868983464629, 13.423352151392551, 13.423549449327021, 14.22902527182818]
19 : [3.2204266279054687, 5.8805624082317598, 5.8805660288197412, 10.686310909960774, 10.694287348061399, 10.694459581162123, 12.31843484526453, 12.319093926355981, 13.358493267223789, 13.423347787043802, 13.423496965102609, 14.223155579578094]

[4]:

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

Eigenvalues
3.22042662791
5.88056240823
5.88056602882
10.68631091
10.6942873481
10.6944595812
12.3184348453
12.3190939264
13.3584932672
13.423347787
13.4234969651
14.2231555796

[5]:

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

[ ]: