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

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.332178541306903, 21.177401462045808, 32.23493344523636, 46.45087047148528, 53.1378348687711, 58.958156450460514, 66.0871129265807, 73.18814222454257, 79.0162401541592, 81.85779772395354, 91.61841211527418, 101.902197650331]
1 : [3.2325472762261125, 6.103050180973414, 6.439573213653074, 12.092149276158048, 12.532612800605296, 13.888431528010354, 16.03272256699882, 17.63414922060643, 18.459313391365644, 19.719396558506983, 21.124484907288984, 25.77899973728244]
2 : [3.220505065350464, 5.891052707622656, 5.8961466352462075, 10.871528019702598, 10.96311339498989, 11.673503423003439, 12.526990777593593, 13.164302198369011, 14.185081702652605, 14.833305076887436, 16.2449502643787, 19.355865109561503]
3 : [3.2203812140997217, 5.881032701617127, 5.881254810190064, 10.721120691153937, 10.758825390612676, 10.90090566481315, 12.362355767712655, 12.528764785041249, 13.562854634596198, 13.735768026357679, 14.874696019159899, 17.550608487953422]
4 : [3.2203797652217454, 5.880559130744092, 5.880591761513983, 10.695978277480569, 10.712951182752915, 10.740605624513874, 12.331706183324822, 12.378233559175145, 13.416204670351496, 13.536267831498684, 14.192583462402016, 16.754081629427358]
5 : [3.22037974603072, 5.880536365844034, 5.88055643744775, 10.69042051905806, 10.699866783741738, 10.70377327122033, 12.32371077152091, 12.336816667942314, 13.380139888206315, 13.477377407170456, 13.851152359363937, 16.191600645754484]
6 : [3.220379745754425, 5.880534903292373, 5.880554828824124, 10.687679044192448, 10.69563599149584, 10.696637346343802, 12.320868378456243, 12.32437068041234, 13.369554771945365, 13.450713369114066, 13.676724979073533, 15.669273391882014]
7 : [3.2203797457502388, 5.88053482295406, 5.880554745817777, 10.686619688393447, 10.694606856210244, 10.695041086904059, 12.319570124256424, 12.320464251514567, 13.365006355532502, 13.436820581118795, 13.578643815202309, 15.194006459799322]
8 : [3.2203797457501726, 5.8805348186268676, 5.880554741437239, 10.68634733280694, 10.694356179748574, 10.694596829217632, 12.318887881632564, 12.319232158037419, 13.362446928061342, 13.429735898457215, 13.517142062668833, 14.825088042275615]
9 : [3.220379745750172, 5.8805348183931425, 5.880554741204548, 10.686283674083771, 10.69427760619781, 10.694487683864288, 12.318518903349549, 12.318861424307551, 13.360825797426912, 13.426268051750679, 13.476769795722682, 14.574845664444904]
10 : [3.2203797457501717, 5.88053481838052, 5.880554741192143, 10.686269015868774, 10.694254043449478, 10.694462467137996, 12.318353721502612, 12.318739919581779, 13.359717022299153, 13.424613572814145, 13.451753853204616, 14.421389527962118]
11 : [3.2203797457501704, 5.880534818379861, 5.880554741191465, 10.68626563741885, 10.694247769084132, 10.694456523457152, 12.318286556542589, 12.318698195742144, 13.358963248250003, 13.423834853774673, 13.437553088064544, 14.331708696832829]
12 : [3.2203797457501717, 5.8805348183798225, 5.880554741191457, 10.686264856149633, 10.694246202506951, 10.694455121118853, 12.318260698646636, 12.318683807854907, 13.358498626832679, 13.423469684459524, 13.43012832268607, 14.280676764437011]
13 : [3.220379745750169, 5.88053481837982, 5.880554741191433, 10.68626467486263, 10.694245822447762, 10.694454791827702, 12.318251076385096, 12.318678837028987, 13.358242421587317, 13.423295536402504, 13.426452052486072, 14.251929459870555]
14 : [3.2203797457501735, 5.880534818379796, 5.88055474119145, 10.68626463278017, 10.694245731736546, 10.694454714884765, 12.318247577500362, 12.31867712027997, 13.35811264464625, 13.423205215707362, 13.424699123619128, 14.235840855751]
15 : [3.220379745750173, 5.880534818379828, 5.880554741191505, 10.686264622971521, 10.694245710253263, 10.694454696960554, 12.31824632037606, 12.318676523812591, 13.358050162281739, 13.423143709962336, 13.42389473244303, 14.226864249233977]
16 : [3.2203797457501753, 5.880534818379657, 5.8805547411913555, 10.686264620679742, 10.694245705166736, 10.69445469278582, 12.318245870011609, 12.3186763159807, 13.358020768648133, 13.423084656388147, 13.423549954723864, 14.22183935762953]
17 : [3.2203797457501744, 5.8805348183798385, 5.880554741191274, 10.68626462014473, 10.6942457039789, 10.694454691818079, 12.318245710999353, 12.318676243468278, 13.358007296136153, 13.42303160632588, 13.42341731116358, 14.219058851416527]
18 : [3.2203797457501753, 5.880534818380067, 5.880554741191894, 10.686264620019578, 10.694245703699156, 10.694454691592579, 12.318245654296078, 12.318676218133264, 13.358001076310817, 13.422996376048248, 13.42336631040825, 14.217494449921118]
19 : [3.2203797457501735, 5.8805348183801325, 5.8805547411913635, 10.686264619990064, 10.694245703630164, 10.69445469153803, 12.318245632693595, 12.318676209025092, 13.357997889705219, 13.422975045614056, 13.423342843612472, 14.216478329851551]

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

Eigenvalues
3.2203797457501735
5.8805348183801325
5.8805547411913635
10.686264619990064
10.694245703630164
10.69445469153803
12.318245632693595
12.318676209025092
13.357997889705219
13.422975045614056
13.423342843612472
14.216478329851551

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

BaseWebGuiScene