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.6520448932278127, 13.603473931023476, 28.424604336856159, 36.164256606872442, 39.983570140291391, 47.557399823699711, 51.101363972858984, 58.21991430924988, 61.806380700938981, 74.915452203938401, 77.689320104490179, 89.651605770076458]
1 : [3.2308991900002004, 6.0070964652732259, 6.6404620349529973, 11.635422133737872, 12.244096013560952, 13.470544851351304, 14.886228771113883, 16.421943155406677, 16.936032215103367, 20.049293010238049, 22.737710788298727, 25.63881817006656]
2 : [3.2206129302462774, 5.8856716277088976, 5.921180255129471, 10.840765908856666, 10.953345842797328, 11.118637567642335, 12.642556295343994, 13.515883284666787, 14.258268117467102, 14.63136218897357, 16.848675890104637, 17.931579779580712]
3 : [3.2204292217075974, 5.8807720691897787, 5.8824708920215967, 10.721839837304337, 10.740880525988178, 10.763011317234628, 12.392276758443431, 12.59691756536404, 13.614860955128293, 13.738534641731095, 14.765758631309948, 16.688558066047886]
4 : [3.2204266669478505, 5.8805721559287765, 5.8806692372052716, 10.695623750211448, 10.702357209242862, 10.706815640933529, 12.337276031313595, 12.396372675042278, 13.464107930289925, 13.499238174029013, 14.239619935026091, 16.264140984588131]
5 : [3.2204266285342849, 5.8805628932977951, 5.8805726191375154, 10.689288202709642, 10.696052371498824, 10.696467809220458, 12.323547152636589, 12.349678718153477, 13.402818178935531, 13.439521503378687, 14.081619099469194, 15.739881611261136]
6 : [3.2204266279160385, 5.8805624339437621, 5.8805664718547703, 10.687234140333892, 10.694725576272788, 10.694816405259187, 12.320079071903084, 12.333782024393805, 13.378892173850495, 13.428225650288832, 13.977240685764951, 15.141561365175312]
7 : [3.2204266279056406, 5.8805624096518461, 5.8805660578667149, 10.686581333103129, 10.694392782068359, 10.694540474498217, 12.3191691648576, 12.326270202621572, 13.370702479943207, 13.425046185229197, 13.850235074421876, 14.675448344540669]
8 : [3.2204266279054621, 5.8805624083107908, 5.8805660306012424, 10.686384911083307, 10.694312768161851, 10.694479325884002, 12.318902742491659, 12.322240513067918, 13.367241818629752, 13.424026140940265, 13.694790754636649, 14.425094790597436]
9 : [3.2204266279054625, 5.8805624082362611, 5.8805660289238197, 10.686330001347574, 10.694293436242321, 10.694464431497142, 12.318786201077636, 12.320254891337632, 13.365109196396439, 13.423661490026319, 13.569146607035563, 14.316945668323971]
10 : [3.2204266279054639, 5.8805624082321239, 5.8805660288251564, 10.686315626062743, 10.69428879130475, 10.694460758468969, 12.318666575303704, 12.319447069011755, 13.363212355758543, 13.423519177937978, 13.494443606228769, 14.26795233761869]
11 : [3.2204266279054616, 5.8805624082318584, 5.8805660288195565, 10.686312048079619, 10.694287687960275, 10.694459864223163, 12.318545674298623, 12.319198204698814, 13.361441258629339, 13.423459972266524, 13.456554172243367, 14.243766591854007]
12 : [3.2204266279054625, 5.8805624082318504, 5.8805660288192119, 10.686311181131755, 10.694287427713048, 10.69445964879629, 12.318477652196821, 12.319127493129274, 13.360067311065029, 13.423434230918621, 13.43861615324465, 14.231081570250646]
13 : [3.2204266279054554, 5.8805624082318628, 5.8805660288192509, 10.686310974211315, 10.694287366568691, 10.694459597270619, 12.318449617787751, 12.319105299068603, 13.359210100787379, 13.423422590390183, 13.430361592692361, 14.224234187076773]
14 : [3.2204266279054612, 5.8805624082318406, 5.8805660288192172, 10.686310925151338, 10.694287352224775, 10.694459584983907, 12.318439204387488, 12.31909782813373, 13.358750644512345, 13.423416912514565, 13.426582792228356, 14.220481611180793]
15 : [3.2204266279054647, 5.880562408232521, 5.8805660288195316, 10.686310913547779, 10.694287348859641, 10.694459582056316, 12.318435475487677, 12.319095245190315, 13.358524694265983, 13.423413406783164, 13.424853543355802, 14.218411590818079]
16 : [3.2204266279054647, 5.8805624082317838, 5.8805660288194987, 10.686310910806215, 10.694287348069931, 10.694459581358513, 12.318434156020277, 12.319094342687576, 13.358418196399484, 13.42340976284575, 13.42406199117274, 14.217264612220685]
17 : [3.2204266279054585, 5.8805624082320342, 5.8805660288195627, 10.686310910157566, 10.694287347884519, 10.694459581192353, 12.318433690882854, 12.319094025683871, 13.358369021639977, 13.423403492125894, 13.423702851537875, 14.216627866291683]
18 : [3.2204266279054603, 5.8805624082318726, 5.8805660288192243, 10.686310910004988, 10.694287347840598, 10.69445958115235, 12.318433527274026, 12.319093914208084, 13.358346541450169, 13.423391167858982, 13.423547244069265, 14.216273676457279]
19 : [3.2204266279054643, 5.8805624082316426, 5.8805660288194135, 10.686310909968713, 10.69428734783064, 10.694459581143571, 12.3184334697533, 12.31909387499101, 13.358336307281524, 13.423373217451811, 13.423487987454623, 14.216076759895182]

[4]:

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

Eigenvalues
3.22042662791
5.88056240823
5.88056602882
10.68631091
10.6942873478
10.6944595811
12.3184334698
12.319093875
13.3583363073
13.4233732175
13.4234879875
14.2160767599

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

[ ]: