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.4737646365557158, 13.901614220037544, 19.418284337333688, 40.312373155888956, 42.682978262889769, 49.319953784648597, 50.55510339828956, 59.288294824938575, 65.347975874583796, 70.191970727590359, 72.420654167832538, 82.621104988754794]
1 : [3.2273253727708502, 6.018110246937642, 6.022669742969418, 11.980318350958981, 12.188096543786223, 14.074351962614179, 14.587174500896859, 15.732585505456575, 17.083676251294115, 18.020932300016575, 19.264555715090303, 23.468158692373695]
2 : [3.2205092844252006, 5.8849584602555325, 5.887228876908372, 10.837136982959077, 10.976412997979363, 11.575780990315444, 12.734613760421519, 13.188941350289332, 13.638197849152213, 14.034740027339007, 15.493541765256307, 16.996735682053952]
3 : [3.2204276980942392, 5.8807527684718686, 5.8809527838131874, 10.715350940174469, 10.750828474780411, 10.970785364864618, 12.379164098796593, 12.527380515005282, 13.463019573834234, 13.600166611420109, 14.524988804077415, 14.874017882250698]
4 : [3.2204266416950338, 5.8805744720311006, 5.8805841285372011, 10.695358599748731, 10.702465138999445, 10.7651382135981, 12.328817765334282, 12.381488819791425, 13.431534919145935, 13.489398883770349, 13.744621995073148, 14.564355071057898]
5 : [3.2204266280881226, 5.8805635850427684, 5.8805664306159882, 10.688895661935462, 10.695517861377187, 10.710829046432446, 12.320657938583979, 12.33881023599699, 13.421550356477082, 13.442092781970279, 13.513025120878249, 14.403173425359036]
6 : [3.2204266279079934, 5.880562471988573, 5.8805660486687641, 10.686955868283201, 10.694589601516048, 10.698110292936944, 12.319184465037859, 12.32533113503063, 13.398642309220639, 13.425775185174796, 13.45029489555964, 14.31723512703689]
7 : [3.2204266279054887, 5.8805624117029076, 5.8805660298376985, 10.686461298502204, 10.694410038650647, 10.695241580585392, 12.318864749590475, 12.320975672858728, 13.37689707982382, 13.423980248862115, 13.433988644121486, 14.271274889444966]
8 : [3.2204266279054585, 5.8805624084220671, 5.8805660288729431, 10.686344795944352, 10.694349815256233, 10.694607210535235, 12.31871319055004, 12.319598811844807, 13.366600077095798, 13.423590270004278, 13.428011467260378, 14.246378545006278]
9 : [3.2204266279054576, 5.8805624082423167, 5.880566028822118, 10.686318558173728, 10.694308760534364, 10.694487301541598, 12.318570106766384, 12.319227534491619, 13.362020222064391, 13.423485361181921, 13.425463388546511, 14.232754329422907]
10 : [3.220426627905471, 5.8805624082324472, 5.8805660288194073, 10.686312653172077, 10.694292844804618, 10.694465615651227, 12.318486797144677, 12.319134517702386, 13.359983860972596, 13.423454590074668, 13.424322030464008, 14.225236815512083]
11 : [3.2204266279054656, 5.880562408231885, 5.880566028819227, 10.68631131095148, 10.694288663327466, 10.694460978519347, 12.318452816488168, 12.319107366957178, 13.359073786688395, 13.423443498219759, 13.423802740864947, 14.221070088630734]
12 : [3.2204266279054607, 5.880562408231838, 5.8805660288192083, 10.686311002823055, 10.694287658321981, 10.694459910074672, 12.318440332545272, 12.319098508811193, 13.358665022754257, 13.423434174850231, 13.423569237258224, 14.218752983420396]
13 : [3.2204266279054581, 5.8805624082318992, 5.8805660288195076, 10.686310931566334, 10.69428742097797, 10.694459658966828, 12.318435877513584, 12.319095479780632, 13.358480663817344, 13.423410484369898, 13.4234812305211, 14.217463527446974]
14 : [3.2204266279054576, 5.8805624082318353, 5.880566028819187, 10.68631091500982, 10.694287365135219, 10.694459599626379, 12.318434301237835, 12.319094426121731, 13.35839739839283, 13.423378118200532, 13.423462021349639, 14.216743905785762]
15 : [3.220426627905451, 5.880562408231957, 5.8805660288191888, 10.686310911141835, 10.694287351926043, 10.694459585535267, 12.318433743752506, 12.319094055709119, 13.358359546333217, 13.423358793816007, 13.423457482087533, 14.216342036883942]
16 : [3.2204266279054536, 5.8805624082318282, 5.8805660288190236, 10.686310910235662, 10.694287348798516, 10.694459582185791, 12.318433546579694, 12.31909392506914, 13.358342309335006, 13.423349292326341, 13.423455924181738, 14.216117293708018]
17 : [3.2204266279054594, 5.8805624082317562, 5.8805660288192065, 10.686310910023087, 10.694287348057669, 10.694459581388861, 12.318433476786311, 12.319093878913851, 13.35833445446197, 13.423344799087127, 13.423455290953671, 14.215991579356935]
18 : [3.220426627905459, 5.8805624082319312, 5.8805660288192785, 10.686310909973127, 10.694287347882044, 10.694459581199345, 12.318433452098072, 12.319093862609929, 13.358330872700231, 13.423342697650542, 13.423455014261709, 14.215921229685735]
19 : [3.2204266279054532, 5.8805624082318078, 5.8805660288191195, 10.68631090996138, 10.694287347840026, 10.694459581153986, 12.31843344335889, 12.319093856842542, 13.358329237934724, 13.423341718967091, 13.42345488940553, 14.215881859082225]

[4]:

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

Eigenvalues
3.22042662791
5.88056240823
5.88056602882
10.68631091
10.6942873478
10.6944595812
12.3184334434
12.3190938568
13.3583292379
13.423341719
13.4234548894
14.2158818591

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

[ ]: