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 : [4.5759499295781785, 16.264858622192783, 26.937513783806338, 36.75114402469773, 40.581172204540088, 47.226886175458873, 56.957013330302502, 62.182464742527024, 68.375677763635437, 76.431615822600506, 82.41300770962485, 88.528019564570045]
1 : [3.227491248658374, 6.066264701841714, 6.7683198761409296, 11.816114228198817, 12.432180199510585, 13.531775911122178, 14.820054441713573, 16.527241483852062, 17.331514694762305, 22.561634957064115, 24.655404397888216, 26.164883939159775]
2 : [3.2204909118354021, 5.8902135929414694, 5.9380903301744921, 10.880820087059471, 10.963632185302837, 11.29198683985036, 12.717991412219201, 13.509957385632003, 14.023513388695392, 16.074164864698993, 17.042644294163495, 20.25099858503949]
3 : [3.220427315978553, 5.8811072560469384, 5.8824379451621178, 10.726822926523237, 10.735309194697379, 10.807881422306398, 12.408600753881297, 12.964696352145848, 13.657867244469799, 13.773201795958668, 15.196434408038595, 16.11995602024318]
4 : [3.2204266367229102, 5.8805911882084745, 5.8806549593451081, 10.694534028809356, 10.702495342897665, 10.721850127989393, 12.341734264042358, 12.661674234366993, 13.426276692651554, 13.568079422333486, 14.36172902721499, 14.806470630416028]
5 : [3.2204266280348257, 5.880564360330105, 5.8805703248286942, 10.688563417922333, 10.69576846853184, 10.701383924229578, 12.32486716509686, 12.469419112124495, 13.375401685847052, 13.48940615877801, 13.853682328317056, 14.456284248418209]
6 : [3.220426627907468, 5.8805625365616168, 5.880566234564724, 10.686893516708162, 10.694576662406595, 10.696129481330585, 12.320590946049119, 12.376368641251299, 13.36312469934923, 13.453409336634456, 13.614901061620159, 14.333148324739991]
7 : [3.2204266279054918, 5.8805624154178648, 5.8805660396052613, 10.686451008537128, 10.694349543485801, 10.694852395462263, 12.319462278154917, 12.339468100269052, 13.359840076409, 13.437115250239771, 13.509501215692675, 14.277391231927696]
8 : [3.2204266279054581, 5.88056240862716, 5.8805660293995663, 10.686343621276341, 10.69430159842695, 10.694550727759102, 12.3191429938282, 12.325907120296758, 13.358856507422084, 13.429682351310017, 13.462428082459237, 14.249049090571393]
9 : [3.220426627905459, 5.8805624082536134, 5.8805660288508568, 10.686318500967531, 10.694290664526106, 10.694480701707159, 12.319033451252967, 12.321116972974457, 13.358528158998098, 13.42629093366129, 13.441176915621032, 14.233985427470953]
10 : [3.2204266279054696, 5.8805624082330272, 5.8805660288209838, 10.686312670910532, 10.694288121189468, 10.694464484703476, 12.318928913232989, 12.319509938691482, 13.358408169196814, 13.424741977914115, 13.431513525948482, 14.225805166355089]
11 : [3.2204266279054559, 5.8805624082319223, 5.88056602881933, 10.686311319161202, 10.694287528966331, 10.694460722013543, 12.318684467694956, 12.319159499559555, 13.358361351637104, 13.424035170096836, 13.427098430141212, 14.221326038470698]
12 : [3.2204266279054501, 5.8805624082318637, 5.8805660288192314, 10.686311005238474, 10.69428739048516, 10.694459846940125, 12.318526402499588, 12.319110983778989, 13.358342220646517, 13.423712582536369, 13.425075244587182, 14.218863201278197]
13 : [3.2204266279054621, 5.8805624082318362, 5.8805660288192145, 10.686310932201076, 10.694287357915652, 10.694459643144095, 12.318466310889846, 12.319099333410509, 13.35833414179381, 13.423564714089624, 13.424146889303161, 14.217507179501753]
14 : [3.2204266279054603, 5.8805624082317038, 5.8805660288188832, 10.686310915157973, 10.694287350220993, 10.694459595629878, 12.318444949553715, 12.319095715749299, 13.358330642869547, 13.423494549066405, 13.4237227126655, 14.216758296425681]
15 : [3.2204266279054545, 5.8805624082318051, 5.8805660288190928, 10.686310911177186, 10.6942873483962, 10.694459584531268, 12.318437468891242, 12.319094500295732, 13.358329104748485, 13.423452318661189, 13.423537604375941, 14.216345056530107]
16 : [3.2204266279054536, 5.8805624082318531, 5.8805660288190662, 10.686310910244131, 10.694287347962751, 10.694459581934259, 12.318434849018056, 12.319094079828115, 13.358328420765449, 13.423406370633582, 13.423478857467803, 14.216116264692339]
17 : [3.2204266279054603, 5.880562408231973, 5.8805660288191124, 10.686310910024972, 10.694287347859461, 10.694459581326091, 12.318433932246817, 12.319093932998662, 13.358328114194652, 13.423372906326625, 13.423464068643263, 14.215989498688716]
18 : [3.220426627905455, 5.8805624082319543, 5.8805660288187962, 10.686310909973718, 10.69428734783496, 10.694459581183782, 12.318433611424592, 12.31909388153189, 13.35832797619398, 13.423355932583197, 13.42345878980125, 14.215919239062512]
19 : [3.2204266279054585, 5.8805624082313379, 5.8805660288192021, 10.686310909961433, 10.694287347828579, 10.694459581149978, 12.318433499225343, 12.319093863480312, 13.358327913853236, 13.423347882658211, 13.423456577612439, 14.215880316941197]

[4]:

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

Eigenvalues
3.22042662791
5.88056240823
5.88056602882
10.68631091
10.6942873478
10.6944595811
12.3184334992
12.3190938635
13.3583279139
13.4233478827
13.4234565776
14.2158803169

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

[ ]: