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 the PINVIT solver from the ngsolve.solvers module.

The orthogonality to gradient fields is important to eliminate the huge number of zero eigenvalues. The projection onto the orthogonal sub-space is implemented by a Poisson projection:

\[P u = u - \nabla \Delta^{-1} \operatorname{div} u\]

The algorithm and the example are taken from the PhD thesis of Sabine Zaglmayr, pp 145-150.

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw

The domain is the Fichera cube. The eigenfunctions are singular at the re-entrant edges and at the re-entrant corner, where we want to apply geometric mesh refinement:

cube1 = Box( (-1,-1,-1), (1,1,1) )

cube2 = Box( (0,0,0), (2,2,2) )
cube2.edges.hpref=1    # mark edges for geometric refinement

fichera = cube1-cube2

Draw (fichera);

mesh = fichera.GenerateMesh(maxh=0.4)
mesh.RefineHP(levels=2, factor=0.2)
Draw (mesh);

We need the stiffness matrix and the mass matrix. The preconditioner is a sparse direct solver for the shifted problem:

fes = HCurl(mesh, order=3)
print ("ndof =", fes.ndof)
u,v = fes.TnT()

a = BilinearForm(curl(u)*curl(v)*dx)
m = BilinearForm(u*v*dx)

apre = BilinearForm(curl(u)*curl(v)*dx + u*v*dx)
pre = preconditioners.Direct(apre, inverse="sparsecholesky")

ndof = 34224

The discrete gradient is a sparse matrix mapping the coefficients of an \(H^1\) function to the coefficients of its gradient in \(H(curl)\). With it we assemble the Laplace matrix for the projection, and combine the projected preconditioner:

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 : [7.198343685792936, 23.643818856996614, 37.677223942105, 47.30298161790627, 50.45402361629184, 58.3136463174334, 63.106920012550574, 75.4799063857739, 83.58942658975575, 89.65970055585778, 102.28276798585539, 106.61233628740752]
1 : [3.236988801908727, 6.0310646004754895, 6.230525952659127, 11.905510217270626, 12.742935648541437, 14.719906189181435, 15.317887288930281, 15.745556255997673, 17.40994660031028, 18.535145529720936, 20.568469656570414, 24.861302538240754]
2 : [3.220853534475484, 5.885170944827662, 5.890562277053273, 10.804832765381429, 10.936019392859583, 11.485916981892776, 12.917511392678758, 13.63559329009197, 13.679306609279143, 14.357128236536354, 15.49295740147936, 18.608597042395314]
3 : [3.2207135653307635, 5.8808504315057135, 5.881200867149412, 10.704891109589173, 10.752772045346843, 10.841509165634374, 12.426692625998285, 13.280972456645625, 13.42717283318252, 13.7294204750973, 14.76723399640874, 16.01649590233829]
4 : [3.220711485551334, 5.88063941097229, 5.880652998766276, 10.691550277211363, 10.709054543686268, 10.734575548721661, 12.339959980520968, 13.001755854711027, 13.38480140515803, 13.548157893741843, 14.213536204184571, 14.584576206164957]
5 : [3.2207114497565232, 5.880607505621905, 5.880636974822869, 10.689230799393199, 10.697799784800644, 10.706582823719549, 12.322970775405167, 12.665956456383162, 13.37349528101254, 13.478406373209289, 13.668498257273159, 14.3860766806634]
6 : [3.220711449130025, 5.880605298568287, 5.88063632342976, 10.688138674027394, 10.695504381528098, 10.698061349875491, 12.319436965354157, 12.453220178884916, 13.368123891527711, 13.447878944793505, 13.510702102219305, 14.305329721944227]
7 : [3.2207114491192868, 5.880605175661456, 5.880636287257649, 10.687565256731299, 10.695043493384851, 10.695988155514891, 12.318654994419447, 12.36678809208487, 13.364016153179236, 13.434044075399386, 13.458128211816787, 14.264146998633489]
8 : [3.220711449119101, 5.880605168984393, 5.880636285246236, 10.687384281028965, 10.694946658559699, 10.695512899589716, 12.318470175138692, 12.335331932720486, 13.361099585247874, 13.427761551064402, 13.437757254251313, 14.242228675527938]
9 : [3.2207114491190976, 5.880605168622143, 5.880636285134871, 10.687338695825767, 10.694925439041484, 10.695403002039143, 12.318423052716549, 12.324253022544346, 13.359398131921507, 13.424895914194746, 13.42924001086856, 14.230343334399327]
10 : [3.2207114491191016, 5.880605168602372, 5.880636285128733, 10.687327886052648, 10.694920711938046, 10.69537739082836, 12.318409671500383, 12.320379311076426, 13.358527136245376, 13.42358447143147, 13.425512200251267, 14.223836264192016]
11 : [3.2207114491190976, 5.8806051686012575, 5.880636285128385, 10.687325356235558, 10.694919650413835, 10.69537140094254, 12.318404100183567, 12.319028890334632, 13.358109821180967, 13.42298135275006, 13.423838900817143, 14.22025391282808]
12 : [3.2207114491190985, 5.8806051686012, 5.880636285128358, 10.68732476472399, 10.694919410082179, 10.695369996165546, 12.318395200325636, 12.318564110590648, 13.357916059461536, 13.422701664650996, 13.423079652061942, 14.218275999736278]
13 : [3.2207114491190953, 5.8806051686012015, 5.880636285128352, 10.687324626207998, 10.694919355331395, 10.695369666005135, 12.31835544772538, 12.318438210203052, 13.35782732656678, 13.422566925101943, 13.422737625633925, 14.217181701062328]
14 : [3.220711449119098, 5.880605168601196, 5.880636285128361, 10.687324593708723, 10.694919342790502, 10.695369588259767, 12.318312784282059, 12.31842284249354, 13.357786932005125, 13.422492700830237, 13.42259243584295, 14.216575639633582]
15 : [3.2207114491190985, 5.880605168601226, 5.880636285128532, 10.68732458607144, 10.694919339909662, 10.695369569923756, 12.318295316431623, 12.318419963855824, 13.3577685879592, 13.42244697434878, 13.422537013643243, 14.216239635769426]
16 : [3.2207114491190985, 5.880605168601349, 5.880636285128535, 10.687324584274116, 10.694919339245583, 10.695369565593724, 12.318289003703592, 12.31841913184365, 13.357760263105215, 13.422421900366704, 13.42251541229861, 14.216053254795625]
17 : [3.220711449119094, 5.880605168601175, 5.880636285128348, 10.687324583850113, 10.694919339091845, 10.69536956456858, 12.318286767546832, 12.318418857391967, 13.357756478792659, 13.42240948823994, 13.422506215050733, 14.215949604833273]
18 : [3.220711449119098, 5.880605168600963, 5.880636285128132, 10.687324583750417, 10.694919339055279, 10.695369564326567, 12.31828597665604, 12.318418762792962, 13.357754737962827, 13.422403586845919, 13.422502010689186, 14.215890299833196]
19 : [3.2207114491190962, 5.880605168601237, 5.88063628512839, 10.68732458372688, 10.6949193390481, 10.695369564269862, 12.318285699576187, 12.318418729762397, 13.357753949981872, 13.422400810364211, 13.422500176895154, 14.215858318466738]

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

Eigenvalues
3.2207114491190962
5.880605168601237
5.88063628512839
10.68732458372688
10.6949193390481
10.695369564269862
12.318285699576187
12.318418729762397
13.357753949981872
13.422400810364211
13.422500176895154
14.215858318466738

We copy the eigenvectors into a multi-dimensional GridFunction for visualization:

gfu = GridFunction(fes, multidim=len(evecs))
for i in range(len(evecs)):
    gfu.vecs[i].data = evecs[i]
Draw (Norm(gfu), mesh, order=4, min=0, max=2);