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

[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.453937861661431, 17.534205875619545, 26.47841523864205, 36.99726308458254, 43.140613399374296, 48.8562068989762, 53.71648896823856, 57.92768236425888, 60.96995589602054, 70.96380617064078, 72.79980569289992, 89.44499107008514]
1 : [3.2408195801556943, 6.094280243101815, 6.308235976921912, 11.791934500052575, 12.16040670501988, 13.501878067111276, 14.687839697143549, 15.71462505811975, 19.70486773734114, 20.800684711959597, 23.58971454628699, 26.142260349144767]
2 : [3.2208390367454918, 5.887760464736464, 5.902942003516158, 10.783529527015027, 11.07089473944574, 11.55399485862232, 12.789893363916732, 13.371076508251155, 14.048599846357698, 15.291151807312684, 17.682595132120746, 19.538388368874024]
3 : [3.2204925262457076, 5.880930085972806, 5.881955696930549, 10.70312618189063, 10.788134259624975, 11.023212712300715, 12.528711517999565, 12.722530541355333, 13.478090252471066, 13.881861864813587, 16.220822004307983, 17.205174259352532]
4 : [3.2204879474529124, 5.88064837934718, 5.880701127616578, 10.690672085763822, 10.722432031456966, 10.79417449769524, 12.443016903395156, 12.48843370604635, 13.423055788549554, 13.532917883885975, 14.872867943629752, 16.657552630939083]
5 : [3.2204878785551463, 5.880606720975348, 5.880637366312938, 10.68860866500417, 10.704412620299728, 10.720206378225448, 12.381271431355666, 12.407151623697777, 13.39870010548514, 13.442688940272069, 14.313526163836046, 16.20524741950666]
6 : [3.2204878773517227, 5.880599563329128, 5.880636529935391, 10.687774990099252, 10.697974317554873, 10.70072422944421, 12.342313194651855, 12.375426628710096, 13.376718847767435, 13.428049094540116, 14.06155434064751, 15.649533546805527]
7 : [3.2204878773276047, 5.880598989053484, 5.880636489135828, 10.687123389172115, 10.695582142776487, 10.696270139232391, 12.327604223467665, 12.353304736885063, 13.366794674530219, 13.425464872592828, 13.89931970915686, 15.091944242217432]
8 : [3.2204878773271153, 5.8805989461456205, 5.8806364870152725, 10.686832965995457, 10.694916939368138, 10.695057813161341, 12.322431478933934, 12.33782305075367, 13.362775841145169, 13.424844106533833, 13.752799329819428, 14.682765362925137]
9 : [3.2204878773270917, 5.88059894322521, 5.880636486903075, 10.68674264365595, 10.694687495970705, 10.694762219409911, 12.320653850143753, 12.328355525602467, 13.361053015961254, 13.424675262278873, 13.622571571665805, 14.454491693720316]
10 : [3.2204878773271073, 5.880598943041277, 5.8806364868970515, 10.686718351323828, 10.694581571644068, 10.69472680006504, 12.320043343491234, 12.323476172265792, 13.360220857879638, 13.424625014194673, 13.53077841006967, 14.34071661181562]
11 : [3.220487877327103, 5.880598943030316, 5.8806364868967105, 10.686712291799909, 10.69455356310315, 10.694718766772978, 12.319829218969344, 12.321317987412824, 13.359751491929652, 13.424608194721792, 13.477744073813657, 14.283483367891225]
12 : [3.2204878773271037, 5.88059894302971, 5.880636486896693, 10.686710835467517, 10.694546553165829, 10.694716937489703, 12.31974555079519, 12.32046367589619, 13.3594617741347, 13.42460081974672, 13.450306282411631, 14.253547858385428]
13 : [3.220487877327097, 5.880598943029653, 5.880636486896706, 10.686710492113749, 10.694544852875481, 10.694716520163238, 12.319704725904664, 12.32015325439601, 13.359289698768576, 13.424594891910855, 13.436902558305833, 14.237454145207293]
14 : [3.220487877327095, 5.880598943029687, 5.880636486896694, 10.686710411817534, 10.69454444672336, 10.6947164248155, 12.319683745992911, 12.320046445281145, 13.359196720341107, 13.424586170076177, 13.430541291186703, 14.228666562216944]
15 : [3.2204878773271015, 5.88059894302949, 5.880636486896668, 10.686710393094865, 10.694544350400745, 10.694716402997848, 12.319674508042944, 12.320009864103572, 13.359150360553798, 13.424570605031997, 13.427573722984555, 14.223832454572735]
16 : [3.2204878773271, 5.8805989430298355, 5.880636486895987, 10.686710388733433, 10.694544327616894, 10.694716397995219, 12.319670927023287, 12.319997123782723, 13.359128368767815, 13.424544727212764, 13.426211986055145, 14.221162065856138]
17 : [3.2204878773271033, 5.880598943029558, 5.8806364868970675, 10.686710387717042, 10.694544322228799, 10.694716396847404, 12.319669617464642, 12.319992632994623, 13.359118194290891, 13.42451048887475, 13.42560233075675, 14.219682082620402]
18 : [3.220487877327098, 5.880598943029606, 5.880636486896766, 10.686710387480348, 10.694544320957668, 10.694716396583711, 12.319669150004556, 12.31999104436804, 13.359113556855586, 13.424477941647208, 13.425337095249954, 14.218862292325133]
19 : [3.2204878773270984, 5.880598943029855, 5.880636486897068, 10.686710387424714, 10.694544320656826, 10.694716396523322, 12.319668984199684, 12.319990480706833, 13.35911145025091, 13.424455198742907, 13.425222052705513, 14.218407643879125]

[4]:

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

Eigenvalues
3.2204878773270984
5.880598943029855
5.880636486897068
10.686710387424714
10.694544320656826
10.694716396523322
12.319668984199684
12.319990480706833
13.35911145025091
13.424455198742907
13.425222052705513
14.218407643879125

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

[ ]: