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.1785526131537924, 24.11186408771013, 33.07240989928315, 39.455820797146956, 51.4810996298624, 56.8619894703534, 61.39943941959258, 74.5559765872194, 77.76955879517354, 91.50320353772041, 108.5334853680693, 110.00589316827069]
1 : [3.241238722504738, 6.119382557470123, 6.718576304355177, 11.466554426641851, 12.59103195764537, 14.309430967609927, 15.812239349892394, 16.24753575204879, 19.472015291349777, 22.704271376187073, 24.709174983589012, 27.083385242256405]
2 : [3.220773299071192, 5.88770831876466, 5.918637537482139, 10.80870169660269, 11.116526242445453, 11.171096572334383, 12.791705119673457, 13.497026107840034, 14.487021344519947, 15.77002641658361, 17.20966690811096, 19.92829368347476]
3 : [3.2204310695045573, 5.880862583506379, 5.882184264757154, 10.716495730968328, 10.759533768274759, 10.848031628216956, 12.457571457434012, 12.897609685328264, 13.521913779050331, 14.48133939367378, 16.014011485865314, 17.893688308939016]
4 : [3.2204266879672874, 5.880578060180525, 5.880635110420944, 10.694928124421093, 10.706568299358775, 10.752717566019435, 12.36417653634213, 12.640895293884673, 13.398551183100603, 14.154446768376346, 14.896322296081424, 16.989696477422033]
5 : [3.2204266288008303, 5.880564085912902, 5.8805684462817265, 10.689482275371732, 10.697435871951063, 10.716123309511524, 12.336653880074689, 12.511511698464199, 13.372266270776683, 13.996946540581776, 14.170510106407763, 16.3391934807208]
6 : [3.2204266279199323, 5.880562508944361, 5.880566139150837, 10.687922991764859, 10.695160698349067, 10.701539087582653, 12.326329479076996, 12.42719732826717, 13.365025196544813, 13.717712748224868, 13.966373274626504, 15.661490418306013]
7 : [3.220426627905696, 5.880562413183012, 5.8805660351366855, 10.687046947903054, 10.694551517649936, 10.696535728351499, 12.321904559090605, 12.3756246332153, 13.362381663769785, 13.550220166145984, 13.862532098577955, 15.036409895371795]
8 : [3.220426627905456, 5.88056240847208, 5.880566029174295, 10.686562019954778, 10.694394931677401, 10.694991019107293, 12.320078222336948, 12.345668872095207, 13.360939280854137, 13.477221870383906, 13.734691762768023, 14.617864725260624]
9 : [3.22042662790545, 5.880562408243535, 5.880566028838662, 10.686380846888508, 10.694337986219036, 10.694573154821354, 12.319380548583613, 12.33014201644682, 13.35996135631104, 13.446296513489603, 13.607711853547237, 14.407920135507094]
10 : [3.2204266279054536, 5.880562408232409, 5.880566028820263, 10.686328279546869, 10.694304243513857, 10.694482481497953, 12.319125581989727, 12.323041571579413, 13.359312293620757, 13.433171577672832, 13.519007652343767, 14.311956935679842]
11 : [3.2204266279054523, 5.880562408231899, 5.880566028819259, 10.686314996573964, 10.694291699786515, 10.694464530094871, 12.319012001682212, 12.320202768092159, 13.358894503199684, 13.42757868482871, 13.469632589128386, 14.26600255522715]
12 : [3.220426627905451, 5.880562408231867, 5.880566028819216, 10.686311846549868, 10.6942883803947, 10.69446067757679, 12.318839584683577, 12.31927959385041, 13.358632583052016, 13.42517865537963, 13.445002657491408, 14.242648334691966]
13 : [3.2204266279054496, 5.88056240823186, 5.880566028819228, 10.68631112202731, 10.694287586480288, 10.694459824638777, 12.318601142587424, 12.319129892240623, 13.358480431108406, 13.424144632842722, 13.433329411028703, 14.230326540210804]
14 : [3.220426627905449, 5.880562408231882, 5.880566028819248, 10.686310957551136, 10.694287402302203, 10.694459634937367, 12.31849257134218, 12.319104026975198, 13.358400334573586, 13.423698225299422, 13.42792667793194, 14.223694695726362]
15 : [3.2204266279054483, 5.880562408231852, 5.880566028819131, 10.686310920578773, 10.694287360197611, 10.694459592940609, 12.318453653536046, 12.31909703446934, 13.358361273250129, 13.4235044757358, 13.42546357575223, 14.22010295078033]
16 : [3.220426627905452, 5.880562408231898, 5.880566028819097, 10.686310912337929, 10.694287350643291, 10.69445958374477, 12.318440322084308, 12.319094886725505, 13.358343037805362, 13.42342058723034, 13.424352654899842, 14.218158717203474]
17 : [3.2204266279054488, 5.88056240823154, 5.880566028819021, 10.686310910492821, 10.694287348470333, 10.694459581716364, 12.31843578145041, 12.319094192841003, 13.35833471007285, 13.423383039096155, 13.423852783809778, 14.217101968042874]
18 : [3.2204266279054514, 5.880562408232618, 5.8805660288193025, 10.686310910077983, 10.694287347973948, 10.694459581267392, 12.318434235296332, 12.319093965422486, 13.358330941294371, 13.42336473709362, 13.423629041277897, 14.216525454129284]
19 : [3.2204266279054488, 5.880562408232063, 5.880566028819317, 10.686310909976925, 10.694287347852386, 10.694459581159933, 12.318433635904364, 12.319093880739013, 13.358329091791552, 13.42335397663019, 13.42351304657985, 14.216143967037148]

[4]:

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

Eigenvalues
3.2204266279054488
5.880562408232063
5.880566028819317
10.686310909976925
10.694287347852386
10.694459581159933
12.318433635904364
12.319093880739013
13.358329091791552
13.42335397663019
13.42351304657985
14.216143967037148

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

[ ]: