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.060341809139501, 21.04354065768085, 31.768231163600465, 35.374502601981916, 40.69516741556049, 46.96742055753494, 55.20089111867111, 61.32835921363698, 65.90151406775699, 70.81647192444339, 80.76123487336712, 102.59834160318533]
1 : [3.240538349139269, 6.038525555231951, 6.354490291544329, 11.488956257138021, 12.283123917999937, 13.562732470432477, 14.112984624438212, 16.494120748981533, 17.033156414135078, 17.77775798112224, 20.026918998295585, 31.13338711288589]
2 : [3.2207032734461962, 5.884894893601808, 5.892302293838372, 10.791799414199906, 11.037988826382934, 11.270780188503462, 12.564737778480985, 12.900313273639672, 13.689393541936063, 14.490171445393642, 15.215250338016945, 18.85672177304041]
3 : [3.2204299687452242, 5.8807788254255335, 5.880999556178517, 10.711248028441991, 10.769258085861559, 10.797159478585563, 12.36740136105479, 12.484566427964413, 13.429073800392597, 13.780478101254893, 14.419074870330524, 15.225214567965573]
4 : [3.220426667189892, 5.880577119121213, 5.880580765301611, 10.696406302511662, 10.705874929960249, 10.71340573386001, 12.329562767103685, 12.380456406446166, 13.381725924087522, 13.540867714392318, 13.920339214292, 14.398826631214765]
5 : [3.2204266283847227, 5.880563177990238, 5.880566675829788, 10.690867225873186, 10.695917217223748, 10.696543510097873, 12.321544938895373, 12.342116108457608, 13.367180719362278, 13.458399924951033, 13.624345660700012, 14.275721665797338]
6 : [3.2204266279117264, 5.880562443973561, 5.880566064033315, 10.687451448891487, 10.69463611709516, 10.694707590112111, 12.319635362442511, 12.327262690624533, 13.361293654039322, 13.433834954903846, 13.508519803444134, 14.24007717386731]
7 : [3.220426627905542, 5.880562409971955, 5.880566030753426, 10.686573280842618, 10.694360469341412, 10.694501806526647, 12.319100325933634, 12.321755509372565, 13.359235457171584, 13.426693213328074, 13.46086437429108, 14.226805711118818]
8 : [3.2204266279054528, 5.880562408320157, 5.8805660289262045, 10.686371029241933, 10.694303402903389, 10.694467705084106, 12.318875420810047, 12.319828585982188, 13.358603541962069, 13.424499337514535, 13.440227843090083, 14.221201718379596]
9 : [3.220426627905453, 5.880562408236479, 5.8805660288251635, 10.686324772488815, 10.694290992635938, 10.694461250635472, 12.318670192218782, 12.319272595375777, 13.358414619476369, 13.423772804965393, 13.431057380741256, 14.218597012763588]
10 : [3.220426627905452, 5.880562408232115, 5.880566028819544, 10.68631412594132, 10.694288190869486, 10.694459941071475, 12.318528796034048, 12.31914525767156, 13.358356752842674, 13.423513259964603, 13.426921747675348, 14.217302768903174]
11 : [3.220426627905454, 5.880562408231877, 5.8805660288192545, 10.686311659938905, 10.694287544966558, 10.694459661370859, 12.318468435548118, 12.319110740854208, 13.358338131811589, 13.423413547627955, 13.425041769986304, 14.216629518836669]
12 : [3.2204266279054536, 5.880562408231869, 5.880566028819231, 10.68631108548934, 10.694287394245828, 10.694459599433685, 12.318445965893655, 12.319099686524234, 13.358331754612507, 13.423372851178137, 13.424182808771793, 14.216269464881538]
13 : [3.2204266279054545, 5.880562408231856, 5.880566028819227, 10.686310951156855, 10.694287358803196, 10.69445958537488, 12.318437890132017, 12.319095904864907, 13.358329421741793, 13.423355415998245, 13.4237892056596, 14.21607352187977]
14 : [3.2204266279054488, 5.880562408231848, 5.880566028819226, 10.686310919644617, 10.694287350429736, 10.694459582130248, 12.31843501709465, 12.319094579954816, 13.358328514666235, 13.42334764559268, 13.423608514963822, 14.215965820882865]
15 : [3.220426627905452, 5.880562408231873, 5.8805660288192225, 10.686310912238973, 10.694287348445505, 10.694459581373216, 12.318433998123401, 12.319094111533554, 13.358328143747313, 13.423344063671797, 13.423525500925681, 14.215906274298353]
16 : [3.2204266279054505, 5.880562408231853, 5.880566028819129, 10.686310910495514, 10.694287347974258, 10.694459581195133, 12.318433636961464, 12.319093945356936, 13.358327986241449, 13.423342375388403, 13.423487341839744, 14.215873234086802]
17 : [3.2204266279054523, 5.880562408231916, 5.8805660288192705, 10.686310910084746, 10.694287347862202, 10.694459581153275, 12.318433508986514, 12.319093886302568, 13.358327917573083, 13.423341576037537, 13.423469794969451, 14.21585485826823]
18 : [3.2204266279054528, 5.880562408231356, 5.8805660288192465, 10.686310909987702, 10.694287347835669, 10.694459581143137, 12.318433463515849, 12.319093865217955, 13.358327887089395, 13.423341198505655, 13.423461685679166, 14.215844538247516]
19 : [3.2204266279054505, 5.88056240823188, 5.880566028819199, 10.686310909964861, 10.69428734782924, 10.694459581140778, 12.318433447399046, 12.319093857734371, 13.358327873430067, 13.423341021162697, 13.42345795226201, 14.215838742259113]

[4]:

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

Eigenvalues
3.2204266279054505
5.88056240823188
5.880566028819199
10.686310909964861
10.69428734782924
10.694459581140778
12.318433447399046
12.319093857734371
13.358327873430067
13.423341021162697
13.42345795226201
14.215838742259113

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

[ ]: