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.8165477377212902, 15.832532822754249, 19.724434185908592, 31.862460675017218, 46.949073422634001, 54.783939130059927, 58.21508575108566, 63.279545432070535, 69.058163714624584, 76.210153392858857, 82.417877278559871, 90.063537391909549]
1 : [3.2427682576720169, 6.0029300391342266, 6.0865952088207287, 11.207765516536222, 13.32111403561804, 13.733193578816728, 15.563214520637256, 17.157654891685027, 17.340015652666839, 19.561440530863074, 21.845462119294957, 24.086163765980924]
2 : [3.22089316805741, 5.8838776201287866, 5.8940825043965495, 10.72148236640651, 11.135842870019614, 11.4925032720954, 12.9173575311001, 13.578856195194668, 13.793098036722562, 14.192522583136922, 15.995123352785896, 18.985736014496236]
3 : [3.2204306293048348, 5.880680529729692, 5.8814287131445582, 10.694991815794255, 10.754542446125321, 10.864102627854269, 12.412828519809214, 13.043271870832678, 13.447224266035704, 13.540363781085842, 14.362576263856086, 16.923978035152846]
4 : [3.2204266844245826, 5.8805696049217069, 5.8806325024704078, 10.691566219603798, 10.700542000429991, 10.752461934559459, 12.341572922687032, 12.769187495779384, 13.394680540139916, 13.440638084125048, 13.916196976802924, 16.301856085601766]
5 : [3.2204266289244745, 5.8805644826384222, 5.8805696042405202, 10.688885231600352, 10.694999357806333, 10.715506788656391, 12.325592151236565, 12.597648121235769, 13.377023185743145, 13.426727006810374, 13.702227444509312, 15.760583143942462]
6 : [3.2204266279249629, 5.8805626743191146, 5.8805661709386277, 10.68761384902353, 10.694479839464012, 10.701372703113089, 12.321111077561703, 12.477027716653359, 13.371122066399874, 13.424138960411726, 13.584809679414205, 15.262123105804342]
7 : [3.2204266279058382, 5.8805624266931575, 5.8805660377345674, 10.686877779005343, 10.694396792461395, 10.696442713553207, 12.319709676351861, 12.39875292483789, 13.368035095192122, 13.42356410870382, 13.51555030423766, 14.870785334249982]
8 : [3.2204266279054656, 5.8805624093957318, 5.8805660293753581, 10.686497744997455, 10.694372173998035, 10.694950657323849, 12.319247429433247, 12.354835967341867, 13.36543166132749, 13.423420669407379, 13.47378408049075, 14.605242658514083]
9 : [3.2204266279054594, 5.880562408301949, 5.8805660288525941, 10.686362821272917, 10.694340182357935, 10.694557894566843, 12.319088877557254, 12.333579276874334, 13.363010408875017, 13.423379657888594, 13.449590635107517, 14.440992197225818]
10 : [3.2204266279054625, 5.8805624082359591, 5.8805660288211659, 10.68632421771836, 10.694306046261385, 10.694479083602296, 12.319029012713973, 12.324402589843203, 13.361065105438669, 13.423366071555801, 13.436475639563014, 14.344167286652096]
11 : [3.2204266279054674, 5.8805624082321133, 5.8805660288193335, 10.68631419903218, 10.694292340900923, 10.69446400728361, 12.318990961475411, 12.32075438426536, 13.359778497200507, 13.423360876340857, 13.429767812606119, 14.288439366780548]
12 : [3.2204266279054634, 5.8805624082318664, 5.8805660288192332, 10.686311708786706, 10.694288583392389, 10.694460634121215, 12.318898853688376, 12.319443555274427, 13.359048317571133, 13.423358534904866, 13.426461348712982, 14.256737591123605]
13 : [3.2204266279054674, 5.8805624082318948, 5.8805660288189561, 10.686311102415798, 10.694287646929888, 10.694459833872799, 12.318669413646749, 12.319153862188024, 13.358672442827679, 13.423357092703917, 13.424868348750806, 14.238827141717056]
14 : [3.2204266279054581, 5.8805624082318273, 5.8805660288192287, 10.686310956133777, 10.694287419687655, 10.694459641796097, 12.318523915605946, 12.319110142447252, 13.358489319891177, 13.423355625526238, 13.424111766476621, 14.228742010851672]
15 : [3.2204266279054634, 5.8805624082312322, 5.8805660288182873, 10.686310921018661, 10.694287365049131, 10.694459595688061, 12.318466483800384, 12.319099197292001, 13.358402722592098, 13.423353544743696, 13.42375657131427, 14.223076264339218]
16 : [3.2204266279054585, 5.8805624082317518, 5.8805660288191879, 10.686310912605409, 10.694287351950997, 10.694459584627253, 12.318445359404777, 12.319095713459699, 13.358362385049675, 13.423350627109345, 13.423591925370106, 14.219895985448444]
17 : [3.220426627905463, 5.8805624082317136, 5.8805660288194188, 10.686310910591486, 10.694287348814244, 10.694459581975861, 12.318437722279276, 12.31909451400001, 13.358343746166948, 13.423347352262134, 13.423516867503412, 14.218111767010836]
18 : [3.2204266279054625, 5.8805624082318353, 5.8805660288191852, 10.686310910109608, 10.69428734806351, 10.694459581340032, 12.318434975993759, 12.319094089812348, 13.358335162003321, 13.423344610210668, 13.42348303153433, 14.217111067629016]
19 : [3.2204266279054696, 5.8805624082312287, 5.880566028818861, 10.686310909994198, 10.694287347884771, 10.694459581187338, 12.318433989981886, 12.319093938311969, 13.35833121522322, 13.423342823305271, 13.423467749359295, 14.216549498191903]

[4]:

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

Eigenvalues
3.22042662791
5.88056240823
5.88056602882
10.68631091
10.6942873479
10.6944595812
12.31843399
12.3190939383
13.3583312152
13.4233428233
13.4234677494
14.2165494982

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

[ ]: