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

[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 : [5.114703744846457, 18.256182852639533, 20.803231126807944, 34.585724229281013, 45.615956752418775, 48.882878831184911, 57.025617427141576, 57.636816553781635, 67.993864085956332, 70.349406070757482, 76.967348182948356, 87.491534241939249]
1 : [3.2273318752288271, 6.0297061798378806, 6.2011564110847406, 11.898112950570022, 13.074565361480143, 13.432909841476452, 14.982384954097659, 16.272206786745034, 19.08117777805327, 20.129548053208424, 22.932420713015969, 26.119945152177511]
2 : [3.2205317371723443, 5.885665131179084, 5.8983209838819075, 10.854983992231464, 11.205169447172407, 11.366526478598804, 12.84434303181466, 13.533640174924052, 14.112595150514233, 14.553843904481205, 15.865045900166695, 20.175650284204107]
3 : [3.2204554386390511, 5.8809043048795475, 5.8814362644245053, 10.731094603646101, 10.788537559816897, 10.905999301901511, 12.512094313352263, 12.716506033102039, 13.539612457781612, 13.766867681053762, 14.566948447925624, 16.301043909699711]
4 : [3.2204543052205281, 5.8806236289567613, 5.8806322345759945, 10.702609708568627, 10.711785610680833, 10.736310638565625, 12.386982807989094, 12.453040963222518, 13.434822464989386, 13.536670404635878, 14.167348884540438, 14.366497334281458]
5 : [3.2204542873288808, 5.8805856629754878, 5.8806183378429582, 10.693907230715633, 10.696407576854936, 10.698600665303962, 12.339250533141559, 12.361013349375478, 13.400209668876467, 13.453193522151079, 13.641622085493934, 14.266030977057129]
6 : [3.2204542870438657, 5.8805838636548762, 5.880617695199053, 10.687638947313783, 10.694895821264545, 10.695241260892812, 12.324585115360948, 12.332100841480104, 13.377432846546576, 13.431601023795467, 13.494045607361237, 14.240175547959021]
7 : [3.220454287039237, 5.8805837759973549, 5.8806176648437267, 10.686657623319983, 10.694568524980816, 10.694626712531177, 12.32079178317859, 12.323450720620537, 13.366332277000994, 13.42667214038015, 13.449723615689363, 14.229075752304045]
8 : [3.220454287039161, 5.8805837716033107, 5.8806176633444389, 10.686509999052111, 10.694459949518134, 10.694531945698719, 12.31988595135836, 12.320728230514925, 13.361977328353204, 13.425277969971033, 13.434413787971419, 14.223633574909401]
9 : [3.2204542870391619, 5.8805837713773261, 5.8806176632673699, 10.686486688268444, 10.694432145592966, 10.694514557840735, 12.319652857877008, 12.319852523273278, 13.360231923867627, 13.42480563045952, 13.428618271531057, 14.220794301030576]
10 : [3.220454287039165, 5.8805837713654752, 5.8806176632632923, 10.686482753171557, 10.6944257815654, 10.694510718481284, 12.319497253876373, 12.319655405284848, 13.359495189529664, 13.42461118952032, 13.426286084979719, 14.219265083319495]
11 : [3.2204542870391628, 5.8805837713648392, 5.8806176632630978, 10.686482040731869, 10.694424326845988, 10.694509844214133, 12.319408134774978, 12.319628871150957, 13.359173715533062, 13.424496831995391, 13.425328953617919, 14.218428200527264]
12 : [3.2204542870391686, 5.8805837713648081, 5.880617663263072, 10.686481902935904, 10.694423993224627, 10.694509642685375, 12.319376080436564, 12.319622261678537, 13.359030936673673, 13.424398923510207, 13.424956347825766, 14.217966088442866]
13 : [3.2204542870391681, 5.8805837713648232, 5.8806176632630587, 10.686481874752163, 10.694423916444434, 10.69450959588133, 12.31936497639329, 12.31962026488041, 13.358966875287241, 13.424322889557429, 13.424821191541247, 14.217709719146791]
14 : [3.2204542870391641, 5.8805837713648073, 5.8806176632630791, 10.686481868727775, 10.694423898737181, 10.694509584957897, 12.319361147587925, 12.319619618854995, 13.358937969375903, 13.424277934616976, 13.424770066715048, 14.217567077628555]
15 : [3.2204542870391721, 5.8805837713648561, 5.8806176632630649, 10.686481867398681, 10.694423894641265, 10.694509582399771, 12.319359825634399, 12.319619402760553, 13.358924873970752, 13.424254975922329, 13.424748951825864, 14.217487575285119]
16 : [3.2204542870391655, 5.8805837713648019, 5.8806176632630587, 10.686481867098625, 10.694423893691923, 10.694509581799261, 12.319359368955407, 12.319619329319154, 13.358918929104016, 13.424243901313465, 13.424739701904596, 14.21744321141534]
17 : [3.2204542870391668, 5.8805837713654228, 5.8806176632628633, 10.686481867030597, 10.694423893470997, 10.694509581658849, 12.319359210259645, 12.319619303938341, 13.358916225136621, 13.424238680519029, 13.42473552410935, 14.217418446279854]
18 : [3.2204542870391655, 5.8805837713648677, 5.8806176632628127, 10.686481867013443, 10.694423893419495, 10.694509581623706, 12.319359154828499, 12.319619295115903, 13.358914957696552, 13.424236215562656, 13.424733606812426, 14.217404606492057]
19 : [3.2204542870391681, 5.8805837713642957, 5.8806176632630045, 10.686481867009128, 10.694423893406952, 10.69450958161535, 12.319359135700145, 12.319619292032446, 13.358914400975229, 13.424235082859225, 13.424732720912901, 14.217396868471825]

[4]:

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

Eigenvalues
3.22045428704
5.88058377136
5.88061766326
10.686481867
10.6944238934
10.6945095816
12.3193591357
12.319619292
13.358914401
13.4242350829
13.4247327209
14.2173968685

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

[ ]: