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.0432128525447082, 16.290102028230702, 29.432225284896955, 35.379873034570451, 40.358223178167286, 45.741924651842353, 58.703129887375368, 59.72233057611318, 63.373845569111296, 70.873535291543845, 78.424233275852274, 80.345590147027366]
1 : [3.2291370887121329, 6.1123367214509816, 6.3718141815916312, 11.627509350734101, 12.817582323519774, 13.981301724576236, 14.930392869074522, 16.700602473267214, 17.01713741545835, 19.745209506293197, 20.781588814730316, 26.042683455029469]
2 : [3.2205577711940352, 5.8913937153760649, 5.9041577487890615, 10.831766999719628, 11.054647546673111, 11.325162251078471, 12.656274560147311, 13.899420007688857, 13.931124443898135, 14.280906757559986, 16.663224912177267, 19.34231515472321]
3 : [3.2204289430632822, 5.8811467850223895, 5.8816366018434634, 10.718021677559333, 10.784266081373479, 10.84150510996967, 12.399662618584079, 12.926929635589914, 13.548090883373108, 13.650843639247748, 14.739557500251641, 16.256576878165774]
4 : [3.2204266710855221, 5.8805972698696918, 5.8806049672372085, 10.698413567490558, 10.716837893557157, 10.722783197950442, 12.341505299688208, 12.525109986587751, 13.453566756136855, 13.518938750878217, 13.759315356695387, 15.399254846768651]
5 : [3.2204266287076071, 5.8805639984497331, 5.8805677920642472, 10.694310248335119, 10.69544224030442, 10.699719902564507, 12.324875417126437, 12.386695717636862, 13.417433556411753, 13.461366785854844, 13.526054444690454, 14.921782478142909]
6 : [3.22042662792007, 5.8805624837219304, 5.8805661159618001, 10.688241192431526, 10.694473538100617, 10.695672522090444, 12.320176750215664, 12.342773462665962, 13.398457994036203, 13.437262465411949, 13.459220386508768, 14.623903272276884]
7 : [3.220426627905721, 5.8805624119911259, 5.8805660333595036, 10.686763737638206, 10.694340188009701, 10.69473234993071, 12.318929283806723, 12.327685795468783, 13.383683458922889, 13.428368755302902, 13.436104917092003, 14.447339852864337]
8 : [3.2204266279054821, 5.8805624084261279, 5.8805660290628321, 10.686417463624469, 10.694303208721241, 10.694520058188814, 12.318598116422562, 12.322199295242848, 13.372231465733254, 13.425194373726661, 13.428369447986205, 14.345888181704943]
9 : [3.2204266279054581, 5.8805624082421479, 5.880566028832515, 10.686336050879467, 10.694291497659131, 10.694473377818845, 12.318500353896816, 12.320198786271622, 13.365311038234186, 13.424092522959953, 13.425515826134903, 14.288599319196917]
10 : [3.2204266279054683, 5.8805624082324197, 5.8805660288199526, 10.68631685021202, 10.694288355309782, 10.694462799785828, 12.318463040959829, 12.319482014161297, 13.361676784653451, 13.423695289581238, 13.42432979116403, 14.25645138843125]
11 : [3.2204266279054647, 5.8805624082318877, 5.8805660288192501, 10.686312315790751, 10.694287587438119, 10.694460337901956, 12.318445727808827, 12.319230041291648, 13.359895361347508, 13.423539469322382, 13.4238087792526, 14.238481415810934]
12 : [3.2204266279054679, 5.8805624082318397, 5.8805660288192341, 10.686311243088092, 10.694287404576116, 10.694459759592499, 12.318438136365469, 12.319141848171745, 13.359052043406228, 13.423466099431256, 13.423581559252476, 14.228452712199021]
13 : [3.220426627905467, 5.8805624082312873, 5.8805660288191239, 10.686310989004262, 10.694287361266886, 10.69445962327859, 12.31843515452262, 12.319110831935754, 13.358660221114876, 13.423413806866119, 13.423497700851154, 14.222863829621376]
14 : [3.2204266279054616, 5.8805624082318566, 5.8805660288191408, 10.686310928766344, 10.694287351012132, 10.6944595911182, 12.318434053759123, 12.31909987482091, 13.358479883581804, 13.42337752555137, 13.42347162439189, 14.219749956836505]
15 : [3.2204266279054719, 5.8805624082318388, 5.8805660288191985, 10.686310914440497, 10.69428734858278, 10.694459583504916, 12.318433657675216, 12.319095991699951, 13.358397292345167, 13.423358356638053, 13.42346206391854, 14.218015524478954]
16 : [3.220426627905463, 5.880562408231877, 5.880566028819227, 10.686310911027224, 10.694287348006691, 10.694459581701313, 12.318433516469383, 12.319094613799262, 13.35835955934688, 13.423349083962382, 13.423458053672096, 14.217049422847369]
17 : [3.2204266279054754, 5.8805624082318619, 5.8805660288191461, 10.686310910213093, 10.694287347869929, 10.694459581273394, 12.318433466245718, 12.319094124111368, 13.358342330735409, 13.423344707196083, 13.423456276250976, 14.216510910507036]
18 : [3.2204266279054714, 5.8805624082317633, 5.8805660288188237, 10.686310910019058, 10.694287347838157, 10.694459581173311, 12.318433448399777, 12.319093949931306, 13.358334465993192, 13.423342657802895, 13.423455469883987, 14.216210598904906]
19 : [3.2204266279054656, 5.8805624082318442, 5.8805660288192767, 10.686310909972233, 10.694287347829674, 10.694459581147523, 12.318433442033969, 12.319093887872354, 13.358330867017589, 13.423341701785992, 13.423455098317731, 14.216042545337421]

[4]:

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

Eigenvalues
3.22042662791
5.88056240823
5.88056602882
10.68631091
10.6942873478
10.6944595811
12.318433442
12.3190938879
13.358330867
13.4233417018
13.4234550983
14.2160425453

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

[ ]: