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.0174048264160067, 14.491458837780387, 23.741060791268662, 29.533323897006461, 38.554226820851994, 40.976392471164644, 47.491446871150899, 50.959929509327161, 63.051129048543459, 70.064528670446066, 74.597435971991146, 87.377662046831787]
1 : [3.2258782253759497, 5.9938707077826416, 6.1621633209944386, 11.545222204825411, 11.785859136793917, 12.660891204973053, 13.992594612326588, 15.748078127425337, 16.563164810680369, 16.973832155595126, 22.278490544472632, 24.466141181052798]
2 : [3.2204993947998468, 5.8838501127757352, 5.890950800844859, 10.812716514096962, 10.960085803141626, 11.271481448634862, 12.594894928784299, 13.003918595964294, 13.857825744312162, 14.498613205915879, 16.291840803486668, 18.569333177903143]
3 : [3.2204275095640287, 5.8806919776905042, 5.8810343929158337, 10.715372823070869, 10.751330776580199, 10.863199949695789, 12.410170000811927, 12.469447552655627, 13.571244106973156, 14.142017001173459, 14.466943577475064, 15.716992808847493]
4 : [3.2204266393494412, 5.8805678392349714, 5.8805892028185704, 10.694244361531933, 10.702445447213103, 10.738679492156182, 12.348628777045212, 12.365635514302083, 13.477943913837667, 13.783322143975813, 13.919756595449423, 14.509725819157621]
5 : [3.2204266280615577, 5.8805626822656194, 5.8805671822844845, 10.688096328721702, 10.695605278921821, 10.704726326096214, 12.326454631702159, 12.335252100524333, 13.44384252904929, 13.56953315251001, 13.592666377281027, 14.289694432990618]
6 : [3.2204266279077216, 5.8805624218007404, 5.8805660888302853, 10.686697254572218, 10.694646944645937, 10.696651320496631, 12.320315662998468, 12.324781229541509, 13.427864215445684, 13.4490461846105, 13.48872487626902, 14.242967211108606]
7 : [3.2204266279055056, 5.8805624089008681, 5.8805660320330553, 10.686396242482717, 10.694488932609969, 10.694823627730555, 12.318873389497341, 12.321073375565911, 13.391639151651015, 13.428274316590523, 13.45178231855286, 14.227508616291558]
8 : [3.2204266279054701, 5.8805624082659662, 5.880566028993953, 10.686329973773972, 10.694403795374559, 10.694470106387742, 12.318544837417983, 12.319776856823662, 13.372243748872974, 13.425429447140813, 13.435964765903359, 14.221330344126331]
9 : [3.220426627905455, 5.8805624082336445, 5.8805660288287962, 10.686315230430198, 10.694314525051462, 10.694461154228756, 12.318464675162431, 12.319328985348214, 13.364294843572223, 13.424287247979722, 13.429034570864289, 14.218582808755588]
10 : [3.2204266279054679, 5.8805624082319508, 5.8805660288197439, 10.686311900222346, 10.694293590357836, 10.694459895119024, 12.318442775238948, 12.319175221444127, 13.360940094187358, 13.423782867401128, 13.425956188041974, 14.217264751193824]
11 : [3.2204266279054656, 5.8805624082318948, 5.8805660288192572, 10.686311138715912, 10.694288783346451, 10.694459647929941, 12.318436337220668, 12.319122154985575, 13.359487496541682, 13.42355380426371, 13.424575932867494, 14.216598096490396]
12 : [3.2204266279054541, 5.8805624082319499, 5.8805660288198265, 10.686310963096656, 10.69428767962588, 10.694459595896372, 12.318434370990486, 12.31910374701982, 13.358847353152791, 13.42344755598927, 13.423955508391272, 14.216248471952907]
13 : [3.2204266279054576, 5.8805624082319898, 5.8805660288191994, 10.686310922347969, 10.69428742474255, 10.69445958448139, 12.318433747309308, 12.31909732204433, 13.358562020750576, 13.423396325601095, 13.423676598754144, 14.216060648149755]
14 : [3.2204266279054585, 5.8805624082317722, 5.8805660288191817, 10.686310912853729, 10.694287365682964, 10.694459581908626, 12.318433543048796, 12.319095071301378, 13.35843382487284, 13.423369905230199, 13.423552481210805, 14.215958242258845]
15 : [3.2204266279054647, 5.8805624082316266, 5.8805660288189729, 10.686310910635861, 10.694287351976891, 10.694459581318744, 12.318433474483388, 12.319094281484645, 13.358375946004386, 13.423355627333482, 13.423498008203774, 14.215901910685295]
16 : [3.2204266279054594, 5.8805624082317873, 5.8805660288190582, 10.686310910116754, 10.694287348794209, 10.69445958118197, 12.318433451059626, 12.319094004339076, 13.358349734218514, 13.423348080082043, 13.423474144451953, 14.215870770249715]
17 : [3.2204266279054607, 5.8805624082329873, 5.8805660288187047, 10.686310909994887, 10.694287348052709, 10.694459581149649, 12.318433442930049, 12.319093906738013, 13.358337824845979, 13.423344291848922, 13.423463542573147, 14.215853459246786]
18 : [3.2204266279054563, 5.8805624082315795, 5.8805660288191817, 10.686310909966892, 10.69428734788009, 10.694459581142453, 12.318433440086411, 12.319093872349091, 13.358332404382217, 13.423342465938417, 13.423458778240461, 14.215843799110218]
19 : [3.2204266279054625, 5.880562408232195, 5.880566028818845, 10.686310909960605, 10.694287347840206, 10.694459581141667, 12.318433439091208, 12.319093860249886, 13.358329934134627, 13.423341609038854, 13.423456611148744, 14.215838438007433]

[4]:

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

Eigenvalues
3.22042662791
5.88056240823
5.88056602882
10.68631091
10.6942873478
10.6944595811
12.3184334391
12.3190938602
13.3583299341
13.423341609
13.4234566111
14.215838438

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

[ ]: