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 : [5.060341809013678, 21.04354065871958, 31.768231164494292, 35.37450260197868, 40.6951674156894, 46.967420557291014, 55.20089111932704, 61.328359213115185, 65.90151406845708, 70.81647192349716, 80.76123487340121, 102.59834160333568]
1 : [3.240538349140024, 6.0385255552333446, 6.354490291585591, 11.48895625716993, 12.283123917976257, 13.562732470580052, 14.1129846243706, 16.49412074905113, 17.033156414101743, 17.777757980785378, 20.026918998526348, 31.13338711283536]
2 : [3.22070327344623, 5.88489489360147, 5.892302293840525, 10.79179941420173, 11.03798882637868, 11.270780188502542, 12.564737778474312, 12.900313273695607, 13.689393541951386, 14.490171445244428, 15.215250338069241, 18.856721772212264]
3 : [3.2204299687452296, 5.880778825425507, 5.880999556178559, 10.711248028441783, 10.769258085858338, 10.797159478580403, 12.367401361054194, 12.484566427974615, 13.429073800394255, 13.780478101184354, 14.419074870271924, 15.225214567485358]
4 : [3.2204266671898933, 5.8805771191212015, 5.880580765301576, 10.696406302511557, 10.705874929959085, 10.71340573385759, 12.329562767103276, 12.380456406445383, 13.381725924088256, 13.540867714362326, 13.920339214120174, 14.398826631171918]
5 : [3.2204266283847267, 5.880563177990255, 5.880566675829804, 10.690867225872529, 10.695917217223082, 10.696543510098248, 12.321544938895148, 12.34211610845574, 13.367180719362477, 13.458399924939666, 13.624345660624886, 14.275721665789414]
6 : [3.220426627911728, 5.880562443973583, 5.880566064033329, 10.687451448891322, 10.694636117095245, 10.694707590112063, 12.31963536244241, 12.327262690623568, 13.361293654039358, 13.433834954899682, 13.508519803413614, 14.240077173863584]
7 : [3.220426627905537, 5.880562409971973, 5.880566030753435, 10.68657328084259, 10.694360469341413, 10.694501806526624, 12.319100325933602, 12.321755509372206, 13.359235457171643, 13.42669321332663, 13.460864374278199, 14.226805711117064]
8 : [3.2204266279054568, 5.880562408320143, 5.880566028926201, 10.686371029241943, 10.694303402903387, 10.694467705084103, 12.318875420809944, 12.319828585982139, 13.358603541962172, 13.424499337513728, 13.440227843086129, 14.221201718377431]
9 : [3.2204266279054528, 5.880562408236483, 5.880566028825153, 10.686324772488828, 10.694290992635956, 10.694461250635454, 12.318670192218752, 12.319272595375821, 13.3584146194765, 13.423772804965441, 13.431057380741578, 14.218597012760949]
10 : [3.2204266279054528, 5.880562408232111, 5.880566028819528, 10.686314125941287, 10.694288190869464, 10.69445994107148, 12.318528796034357, 12.319145257671439, 13.35835675284353, 13.423513259967326, 13.426921747677236, 14.217302768950342]
11 : [3.2204266279054554, 5.88056240823188, 5.8805660288192385, 10.686311659938939, 10.694287544966532, 10.694459661370882, 12.31846843554778, 12.319110740853722, 13.358338131811944, 13.423413547629895, 13.425041769979234, 14.216629518864941]
12 : [3.2204266279054536, 5.880562408231855, 5.8805660288192225, 10.686311085489258, 10.69428739424579, 10.694459599433664, 12.318445965903463, 12.319099686522158, 13.358331754622478, 13.423372851188073, 13.424182808510421, 14.216269464439813]
13 : [3.2204266279054528, 5.880562408231871, 5.88056602881921, 10.68631095115662, 10.69428735880316, 10.694459585374855, 12.318437890126006, 12.31909590485858, 13.358329421744033, 13.423355415980557, 13.4237892043755, 14.21607352153974]
14 : [3.2204266279054448, 5.880562408231864, 5.880566028819202, 10.68631091964588, 10.694287350429898, 10.694459582130264, 12.318435017172947, 12.319094579929912, 13.358328514664601, 13.423347646275383, 13.423608513979808, 14.21596582933164]
15 : [3.2204266279054488, 5.880562408231867, 5.8805660288191275, 10.686310912238778, 10.694287348445414, 10.694459581373124, 12.31843399802235, 12.319094111478112, 13.358328143734251, 13.423344063970035, 13.423525487507817, 14.215906278524717]
16 : [3.220426627905456, 5.8805624082318655, 5.880566028819165, 10.686310910495441, 10.694287347974237, 10.69445958119521, 12.318433636889118, 12.319093945322502, 13.358327986231162, 13.42334237549995, 13.42348733076236, 14.215873236010905]
17 : [3.220426627905452, 5.880562408231898, 5.880566028819201, 10.686310910084863, 10.694287347862362, 10.694459581153103, 12.31843350893693, 12.319093886284033, 13.358327917573668, 13.423341575849312, 13.423469787008905, 14.215854856829278]
18 : [3.2204266279054496, 5.880562408232052, 5.880566028819309, 10.68631090998707, 10.694287347835358, 10.694459581143203, 12.318433463003844, 12.319093865022044, 13.358327887002732, 13.423341194456077, 13.423461637863314, 14.215844599070758]
19 : [3.2204266279054536, 5.880562408232316, 5.880566028819049, 10.686310909964527, 10.694287347829125, 10.694459581140789, 12.318433447051335, 12.319093857601642, 13.358327872707632, 13.42334101225094, 13.423457860933748, 14.215838367081432]

[4]:

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

Eigenvalues
3.2204266279054536
5.880562408232316
5.880566028819049
10.686310909964527
10.694287347829125
10.694459581140789
12.318433447051335
12.319093857601642
13.358327872707632
13.42334101225094
13.423457860933748
14.215838367081432

[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)

[ ]: