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 : [5.060341808875264, 21.04354065785334, 31.768231163583664, 35.37450260035117, 40.69516741472638, 46.967420554979526, 55.200891114731846, 61.32835921594658, 65.90151406259643, 70.81647192357606, 80.76123487613884, 102.59834160498478]
1 : [3.240538349136266, 6.038525555223775, 6.3544902915307295, 11.488956257132685, 12.28312391800568, 13.56273247027122, 14.112984624373397, 16.49412074978619, 17.033156415044534, 17.777757979931504, 20.026918999272485, 31.13338711589483]
2 : [3.2207032734462344, 5.8848948936025245, 5.892302293839044, 10.791799414207489, 11.03798882641334, 11.270780188779945, 12.564737778480506, 12.900313273660599, 13.689393541920317, 14.490171445767619, 15.2152503383877, 18.85672177440141]
3 : [3.2204299687452345, 5.880778825425665, 5.880999556178424, 10.711248028443423, 10.769258085903783, 10.797159478625934, 12.367401361056162, 12.484566427963824, 13.429073800381065, 13.78047810146376, 14.419074870363309, 15.225214567806198]
4 : [3.220426667189925, 5.880577119121237, 5.8805807653016275, 10.696406302512106, 10.705874929968408, 10.713405733866175, 12.329562767104331, 12.380456406442173, 13.381725924085481, 13.540867714459576, 13.920339214126718, 14.398826631214966]
5 : [3.2204266283847427, 5.880563177990217, 5.880566675829729, 10.690867225873284, 10.695917217223135, 10.696543510101026, 12.321544938895588, 12.342116108454363, 13.367180719364242, 13.458399924957934, 13.62434566062323, 14.275721665799768]
6 : [3.2204266279117126, 5.880562443973469, 5.88056606403338, 10.68745144889165, 10.694636117095731, 10.694707590111987, 12.319635362442437, 12.327262690622835, 13.36129365404012, 13.43383495489981, 13.508519803414925, 14.240077173866771]
7 : [3.2204266279055194, 5.880562409971976, 5.880566030753558, 10.68657328084282, 10.694360469341701, 10.694501806526796, 12.319100325933825, 12.321755509371759, 13.3592354571716, 13.426693213324198, 13.46086437427702, 14.226805711117388]
8 : [3.2204266279054576, 5.8805624083201256, 5.880566028926249, 10.686371029242054, 10.694303402903492, 10.69446770508401, 12.318875420809938, 12.319828585981892, 13.358603541962005, 13.424499337512785, 13.440227843085497, 14.221201718381405]
9 : [3.220426627905429, 5.880562408236524, 5.880566028825162, 10.686324772488776, 10.694290992635999, 10.694461250635484, 12.318670192218818, 12.31927259537573, 13.358414619475692, 13.423772804963082, 13.4310573807433, 14.218597012765011]
10 : [3.2204266279054914, 5.880562408232126, 5.880566028819612, 10.686314125941342, 10.694288190869413, 10.694459941071484, 12.318528796034276, 12.319145257671359, 13.358356752844285, 13.423513259962824, 13.426921747674681, 14.21730276888645]
11 : [3.220426627905472, 5.880562408231919, 5.880566028819237, 10.68631165993891, 10.69428754496636, 10.694459661370812, 12.318468435548047, 12.319110740853894, 13.358338131812515, 13.423413547627552, 13.425041770005295, 14.216629518832796]
12 : [3.220426627905507, 5.880562408231844, 5.88056602881928, 10.686311085489466, 10.694287394245825, 10.694459599433657, 12.318445965907555, 12.31909968652162, 13.358331754611969, 13.42337285117249, 13.424182808452228, 14.216269464819888]
13 : [3.2204266279054705, 5.880562408231875, 5.880566028819295, 10.686310951156932, 10.694287358803304, 10.694459585374977, 12.318437890139634, 12.319095904860783, 13.358329421732957, 13.423355416053326, 13.423789205322636, 14.216073521776828]
14 : [3.2204266279054616, 5.8805624082319055, 5.88056602881925, 10.686310919645706, 10.694287350429859, 10.69445958213039, 12.318435017059285, 12.319094579902327, 13.358328514639956, 13.42334764632726, 13.423608498729688, 14.215965826838612]
15 : [3.22042662790545, 5.880562408231778, 5.880566028819256, 10.686310912239227, 10.694287348445457, 10.694459581373293, 12.31843399808524, 12.319094111478366, 13.35832814364686, 13.423344063939362, 13.423525489232631, 14.215906272248871]
16 : [3.220426627905453, 5.880562408231899, 5.880566028819275, 10.686310910495676, 10.694287347974376, 10.694459581195124, 12.318433636937694, 12.319093945327893, 13.35832798620366, 13.423342375557056, 13.423487331733485, 14.215873226451754]
17 : [3.2204266279054545, 5.880562408232065, 5.8805660288192225, 10.686310910084593, 10.694287347861964, 10.694459581153076, 12.318433508865628, 12.319093886267872, 13.358327917561441, 13.423341575723812, 13.423469769582352, 14.215854850983886]
18 : [3.2204266279054505, 5.880562408231923, 5.880566028819496, 10.686310909987583, 10.694287347835612, 10.694459581143228, 12.31843346347152, 12.319093865285758, 13.358327887110198, 13.423341198461442, 13.42346168971226, 14.215844618844512]
19 : [3.2204266279054465, 5.8805624082318575, 5.880566028819369, 10.686310909964906, 10.694287347829011, 10.694459581140812, 12.318433447255392, 12.319093857769015, 13.358327873382896, 13.423341019198416, 13.423457929896784, 14.215838795183585]

[4]:

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

Eigenvalues
3.2204266279054465
5.8805624082318575
5.880566028819369
10.686310909964906
10.694287347829011
10.694459581140812
12.318433447255392
12.319093857769015
13.358327873382896
13.423341019198416
13.423457929896784
14.215838795183585

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

[ ]: