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.805905055203783, 24.468823841970732, 38.774015756435645, 42.352873784881105, 43.49334859448238, 47.10794240966373, 56.298906364777146, 62.664800128410356, 75.8203937056083, 84.82026482544933, 97.80562848462382, 102.90297344777693]
1 : [3.241772881797568, 6.313486273498762, 6.690873833742419, 11.622157174940877, 12.696460506784812, 12.808301463484682, 14.000625796081314, 15.495738214689338, 16.619777102284058, 19.815082661280982, 25.43687166265005, 27.418828133620142]
2 : [3.2207829548517752, 5.902332257213175, 5.922499099018044, 10.785126210366629, 11.151236588204208, 11.334867763994218, 12.639999246280544, 13.044739404579705, 13.913064131277714, 15.338626837096593, 16.68780946538464, 19.7954184622484]
3 : [3.2203962875978354, 5.881904392752992, 5.882946288245242, 10.71365322445965, 10.779985246769186, 10.877888091476132, 12.39340374417508, 12.735451537782561, 13.573511364472754, 13.988245377679048, 14.754681838310301, 16.183227622677094]
4 : [3.2203912433793866, 5.88062909763651, 5.88071926796872, 10.69911839323937, 10.70529613637703, 10.742050645528177, 12.336771542895208, 12.561645810153662, 13.469599961463693, 13.657530152902392, 14.042780914160804, 14.888856114728002]
5 : [3.2203911736236854, 5.880562847565554, 5.880576851382117, 10.690713905325135, 10.695604350419277, 10.705725004341469, 12.323666703084468, 12.437569669509548, 13.431899646181362, 13.530201465749707, 13.685509242666335, 14.46406787146713]
6 : [3.22039117262611, 5.880558617392552, 5.880569293774326, 10.687480127940404, 10.694747947562021, 10.69692331949966, 12.320368058418962, 12.367851822139183, 13.414066730777, 13.470072042825478, 13.5196601305553, 14.32298251292359]
7 : [3.220391172611192, 5.880558377452086, 5.880568896124853, 10.686651836218134, 10.69453433684773, 10.694995500288167, 12.319455842180371, 12.33734367862189, 13.396282426277747, 13.441512726182982, 13.456871779298014, 14.267582636958215]
8 : [3.2203911726109617, 5.8805583651864275, 5.880568874510894, 10.686450487638051, 10.694428454702622, 10.694629635117021, 12.319183450443594, 12.325642518140311, 13.378478726028373, 13.430442030024148, 13.437031365264485, 14.242514888197999]
9 : [3.220391172610959, 5.880558364557212, 5.880568873345895, 10.686404048484476, 10.694374205903303, 10.694578467308638, 12.319097086041129, 12.321423722115037, 13.367973135581789, 13.426618705576447, 13.42954777538619, 14.230130272161722]
10 : [3.220391172610957, 5.880558364524578, 5.880568873283342, 10.686393515367723, 10.694359635961357, 10.694569234895646, 12.319067727170069, 12.319944102797919, 13.362897818113606, 13.425131823588268, 13.42632849985173, 14.223700897140532]
11 : [3.220391172610957, 5.880558364522876, 5.880568873279988, 10.686391132542333, 10.694356230378258, 10.69456723012043, 12.319055636342483, 12.319432813397045, 13.360572920131032, 13.424492765477133, 13.424900682669179, 14.220272241312102]
12 : [3.2203911726109564, 5.880558364522763, 5.880568873279808, 10.686390592815343, 10.694355450176182, 10.694566774299547, 12.319048518090822, 12.319259129023488, 13.359526594398979, 13.42414838353558, 13.424318123023111, 14.218417301573549]
13 : [3.220391172610959, 5.880558364522735, 5.880568873279766, 10.686390470302985, 10.694355271663472, 10.69456666938956, 12.319044061344613, 12.319201188093869, 13.359058630997191, 13.423900271803511, 13.424146720107105, 14.217405932793552]
14 : [3.220391172610955, 5.880558364522807, 5.880568873279846, 10.68639044243814, 10.694355230743273, 10.694566645144327, 12.319041879951687, 12.319181833457138, 13.35884978513635, 13.423773649431803, 13.42408266643354, 14.216852149997601]
15 : [3.22039117261096, 5.880558364522716, 5.880568873279969, 10.6863904360882, 10.694355221342063, 10.694566639528627, 12.319041004258343, 12.319175272073515, 13.358756633166223, 13.42371467105448, 13.424054766592871, 14.216548127739914]
16 : [3.220391172610958, 5.880558364522374, 5.880568873279596, 10.686390434638355, 10.694355219178755, 10.694566638225316, 12.319040683068197, 12.319173021635569, 13.358715060316314, 13.423687567712001, 13.424042252458253, 14.21638098117391]
17 : [3.220391172610958, 5.880558364522792, 5.8805688732798265, 10.686390434307592, 10.694355218680462, 10.694566637922522, 12.31904056873523, 12.319172245980942, 13.3586964626965, 13.423675149685163, 13.424036526836819, 14.21628904029817]
18 : [3.220391172610954, 5.880558364522984, 5.88056887327977, 10.68639043423135, 10.694355218563098, 10.69456663785109, 12.319040528371763, 12.319171973830882, 13.358688080692666, 13.423669332665293, 13.424033940754557, 14.216236707338798]
19 : [3.220391172610954, 5.880558364522706, 5.880568873279493, 10.686390434214186, 10.694355218536641, 10.694566637835248, 12.319040514751405, 12.319171882489615, 13.358684417479598, 13.423666637222345, 13.424032783705929, 14.21620754177143]

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

Eigenvalues
3.220391172610954
5.880558364522706
5.880568873279493
10.686390434214186
10.694355218536641
10.694566637835248
12.319040514751405
12.319171882489615
13.358684417479598
13.423666637222345
13.424032783705929
14.21620754177143

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

BaseWebGuiScene