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.060341809021796, 21.043540658215584, 31.768231163900285, 35.374502601547974, 40.695167415895824, 46.96742055775477, 55.20089111876471, 61.328359214249325, 65.90151406954449, 70.81647192119154, 80.761234873308, 102.59834160396409]
1 : [3.2405383491398103, 6.038525555225776, 6.354490291479194, 11.488956257204713, 12.283123918024172, 13.562732470549372, 14.11298462449869, 16.494120749237048, 17.033156413834387, 17.77775798071447, 20.02691899793314, 31.133387113387073]
2 : [3.2207032734461993, 5.884894893601546, 5.892302293837165, 10.791799414214466, 11.037988826424616, 11.270780188471035, 12.56473777845509, 12.900313273631479, 13.689393541938175, 14.490171445478229, 15.215250337880542, 18.856721772659387]
3 : [3.220429968745232, 5.880778825425536, 5.880999556178469, 10.711248028444707, 10.769258085869843, 10.79715947858373, 12.367401361048469, 12.48456642797067, 13.429073800396653, 13.780478101278064, 14.419074870278811, 15.225214567765319]
4 : [3.2204266671898947, 5.880577119121202, 5.880580765301571, 10.696406302512322, 10.705874929961814, 10.713405733860286, 12.329562767102471, 12.380456406447365, 13.38172592409087, 13.540867714391728, 13.920339214245033, 14.398826631174023]
5 : [3.220426628384726, 5.880563177990247, 5.880566675829788, 10.69086722587354, 10.695917217223922, 10.69654351009792, 12.321544938895114, 12.34211610845763, 13.367180719364153, 13.45839992494765, 13.62434566068161, 14.275721665785902]
6 : [3.2204266279117246, 5.880562443973586, 5.880566064033337, 10.68745144889159, 10.694636117095278, 10.69470759011206, 12.319635362442485, 12.327262690624488, 13.361293654040034, 13.433834954901798, 13.508519803437371, 14.240077173862915]
7 : [3.2204266279055456, 5.880562409971965, 5.880566030753433, 10.686573280842616, 10.694360469341438, 10.694501806526624, 12.31910032593364, 12.321755509372544, 13.359235457171863, 13.426693213326967, 13.46086437428854, 14.226805711116974]
8 : [3.2204266279054505, 5.8805624083201335, 5.880566028926194, 10.686371029241933, 10.694303402903394, 10.69446770508409, 12.31887542081007, 12.319828585982133, 13.358603541962093, 13.424499337514431, 13.440227843088502, 14.221201718380296]
9 : [3.220426627905451, 5.880562408236457, 5.880566028825161, 10.686324772488797, 10.694290992635956, 10.694461250635479, 12.318670192218729, 12.319272595375802, 13.358414619476665, 13.423772804966159, 13.431057380737705, 14.218597012765951]
10 : [3.220426627905453, 5.880562408232086, 5.880566028819557, 10.686314125941264, 10.694288190869456, 10.694459941071491, 12.318528796034268, 12.319145257671078, 13.358356752843944, 13.423513259967464, 13.426921747659636, 14.217302768947492]
11 : [3.220426627905458, 5.880562408231844, 5.880566028819228, 10.686311659938898, 10.694287544966528, 10.694459661370862, 12.318468435547532, 12.319110740853448, 13.358338131812369, 13.423413547630258, 13.425041769957337, 14.216629518880724]
12 : [3.2204266279054505, 5.880562408231852, 5.880566028819207, 10.686311085489292, 10.694287394245789, 10.694459599433618, 12.318445965898887, 12.319099686522078, 13.358331754619977, 13.423372851181659, 13.424182808603916, 14.216269464521238]
13 : [3.220426627905452, 5.8805624082318495, 5.880566028819218, 10.6863109511573, 10.694287358803251, 10.694459585374847, 12.31843789014494, 12.319095904865833, 13.358329421752941, 13.423355416079056, 13.423789206583328, 14.216073521589783]
14 : [3.2204266279054514, 5.88056240823184, 5.880566028819221, 10.686310919646095, 10.694287350429926, 10.694459582130264, 12.318435017176363, 12.319094579940588, 13.358328514623588, 13.423347646471322, 13.42360851474584, 14.21596583310655]
15 : [3.220426627905451, 5.88056240823189, 5.8805660288190635, 10.686310912239474, 10.694287348445494, 10.694459581373103, 12.318433998179218, 12.319094111519851, 13.35832814371205, 13.423344064228855, 13.42352550198477, 14.21590628129104]
16 : [3.2204266279054505, 5.880562408232025, 5.880566028819208, 10.686310910495205, 10.694287347973876, 10.694459581195128, 12.318433636955056, 12.319093945332542, 13.358327986219859, 13.42334237546033, 13.423487339995404, 14.215873233647622]
17 : [3.220426627905455, 5.880562408231221, 5.880566028819004, 10.686310910084613, 10.694287347862225, 10.694459581153875, 12.318433508900851, 12.319093886269638, 13.35832791755526, 13.42334157557371, 13.4234697801655, 14.215854853063288]
18 : [3.2204266279054576, 5.880562408231718, 5.880566028819565, 10.686310909986812, 10.69428734783524, 10.694459581143263, 12.318433462436898, 12.319093864459447, 13.358327886369459, 13.42334119250432, 13.423461581789951, 14.215844604639585]
19 : [3.220426627905456, 5.880562408231206, 5.880566028819229, 10.686310909963838, 10.694287347828881, 10.694459581140705, 12.318433445969372, 12.319093856825335, 13.35832787266954, 13.423341010834509, 13.423457718396559, 14.21583887677079]

[4]:

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

Eigenvalues
3.220426627905456
5.880562408231206
5.880566028819229
10.686310909963838
10.694287347828881
10.694459581140705
12.318433445969372
12.319093856825335
13.35832787266954
13.423341010834509
13.423457718396559
14.21583887677079

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

[ ]: