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.

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.csg 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)

BaseWebGuiScene

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

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.805905053557727, 24.46882385220502, 38.77401575718461, 42.352873780160536, 43.49334859296999, 47.1079424082101, 56.298906350370245, 62.664800112220526, 75.82039370743651, 84.820264829446, 97.80562848061757, 102.90297344604832]
1 : [3.241772881792549, 6.313486273589134, 6.690873834396933, 11.62215717462148, 12.696460506428936, 12.80830146193376, 14.000625794327643, 15.495738213832214, 16.619777102585733, 19.815082664233863, 25.43687165736177, 27.418828134239938]
2 : [3.2207829548518103, 5.902332257211285, 5.922499099050752, 10.785126210332146, 11.151236587873235, 11.334867763619272, 12.639999246386878, 13.044739404798252, 13.913064131050742, 15.33862683635111, 16.687809464049813, 19.795418460951044]
3 : [3.2203962875978274, 5.881904392752026, 5.882946288245713, 10.713653224451534, 10.77998524672856, 10.877888091283657, 12.393403744226944, 12.735451537671215, 13.573511364417032, 13.988245377347704, 14.754681837180192, 16.183227622229257]
4 : [3.2203912433793893, 5.880629097636433, 5.880719267968796, 10.699118393236883, 10.705296136375114, 10.742050645487575, 12.336771542900971, 12.561645810025066, 13.469599961451973, 13.65753015289985, 14.042780913735383, 14.888856115109574]
5 : [3.220391173623686, 5.8805628475655665, 5.880576851382128, 10.690713905326929, 10.69560435041753, 10.705725004337529, 12.323666703083884, 12.437569669463471, 13.431899646204155, 13.530201465751418, 13.68550924270335, 14.464067871656248]
6 : [3.220391172626103, 5.880558617392552, 5.880569293774309, 10.687480127941127, 10.694747947561515, 10.696923319499748, 12.320368058418257, 12.36785182213166, 13.414066730820906, 13.470072042800483, 13.519660130634412, 14.322982513006714]
7 : [3.2203911726111873, 5.8805583774520995, 5.880568896124865, 10.686651836218347, 10.694534336847541, 10.694995500288282, 12.319455842180062, 12.337343678621965, 13.396282426324634, 13.441512726147426, 13.456871779340359, 14.267582636995526]
8 : [3.220391172610964, 5.880558365186421, 5.880568874510881, 10.686450487638094, 10.694428454702628, 10.694629635117057, 12.319183450443424, 12.325642518139631, 13.378478726055414, 13.430442030010328, 13.437031365279593, 14.242514888220525]
9 : [3.2203911726109573, 5.880558364557229, 5.880568873345912, 10.686404048484595, 10.694374205903493, 10.694578467308686, 12.319097086041149, 12.32142372211176, 13.367973135562126, 13.426618705570887, 13.429547775422371, 14.230130272256945]
10 : [3.2203911726109573, 5.8805583645245925, 5.880568873283345, 10.686393515367744, 10.694359635961476, 10.69456923489574, 12.31906772716871, 12.319944102815635, 13.36289781812346, 13.425131823546728, 13.42632849995805, 14.223700897307303]
11 : [3.220391172610955, 5.8805583645228525, 5.880568873279972, 10.686391132542319, 10.694356230379032, 10.694567230120487, 12.3190556363335, 12.319432813267074, 13.360572919299806, 13.424492765607834, 13.424900683197745, 14.220272244031595]
12 : [3.220391172610956, 5.880558364522762, 5.880568873279804, 10.686390592814206, 10.694355450179522, 10.694566774299851, 12.319048518073986, 12.319259128543099, 13.35952659023366, 13.424148389262566, 13.424318122655793, 14.218417319279645]
13 : [3.220391172610956, 5.880558364522726, 5.88056887327979, 10.686390470303529, 10.694355271674338, 10.69456666939384, 12.319044061534912, 12.319201190193585, 13.359058625166288, 13.423900294667302, 13.424146719214088, 14.217406020980851]
14 : [3.2203911726109524, 5.880558364522732, 5.880568873279838, 10.686390442439302, 10.694355230750734, 10.694566645146967, 12.319041880213463, 12.3191818359678, 13.358849792885522, 13.423773682755916, 13.424082667315329, 14.216852314616933]
15 : [3.2203911726109538, 5.880558364522738, 5.880568873279882, 10.68639043608861, 10.694355221345054, 10.694566639529203, 12.319041004316285, 12.319175272754784, 13.358756630695014, 13.42371470122686, 13.42405476642059, 14.216548278114589]
16 : [3.220391172610954, 5.880558364522679, 5.880568873279799, 10.686390434638957, 10.694355219180249, 10.694566638225666, 12.319040683111593, 12.319173022600483, 13.358715070692165, 13.423687580888872, 13.424042249795393, 14.216381099701126]
17 : [3.2203911726109564, 5.8805583645229715, 5.880568873279794, 10.686390434307711, 10.694355218680107, 10.69456663792268, 12.319040569190499, 12.3191722469747, 13.358696522203735, 13.423675146209028, 13.424036578237269, 14.216289049491316]
18 : [3.2203911726109538, 5.880558364522677, 5.880568873279719, 10.686390434231512, 10.694355218564642, 10.694566637852194, 12.319040529054645, 12.319171976895378, 13.358688126283383, 13.423669432731485, 13.424033977150636, 14.216238207911566]
19 : [3.220391172610957, 5.880558364524128, 5.880568873280457, 10.68639043421188, 10.694355218538123, 10.694566637836846, 12.319040514224007, 12.319171871651434, 13.358683684126435, 13.423666810438696, 13.424032726635769, 14.216209246923961]

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

Eigenvalues
3.220391172610957
5.880558364524128
5.880568873280457
10.68639043421188
10.694355218538123
10.694566637836846
12.319040514224007
12.319171871651434
13.358683684126435
13.423666810438696
13.424032726635769
14.216209246923961

gfu = GridFunction(fes, multidim=len(evecs))
for i in range(len(evecs)):
    gfu.vecs[i].data = evecs[i]
Draw (Norm(gfu), mesh, "u")

BaseWebGuiScene