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.

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.csg 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)

BaseWebGuiScene

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 = 32580

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.822091761907151, 24.102342217391183, 25.410210864089226, 31.59161008380949, 41.04286265795193, 54.05306393303934, 55.57754470187537, 62.489226899201384, 73.59737629868052, 84.21133465211875, 88.22650434877126, 101.57701298119346]
1 : [3.2501484717532985, 6.125216794672185, 6.182311947938549, 11.524201703379122, 12.23939146786013, 13.362041719406646, 14.80507377761941, 15.244167320032584, 18.536522771343467, 18.834172333565363, 20.65939144739967, 23.476093198088545]
2 : [3.22077447924718, 5.888381219782118, 5.89039545102236, 10.77316730607426, 10.984625590010687, 11.204618750276888, 12.638158094582213, 13.123265344500265, 13.800547755949482, 14.995199239120742, 15.571911838182874, 17.48216978235594]
3 : [3.220396162650975, 5.880863278506517, 5.8810436007038955, 10.703264199588986, 10.753658139277151, 10.83927913592552, 12.436787456498916, 12.591167821094105, 13.46744598829433, 13.973218319095759, 14.279890742738763, 16.930379822339038]
4 : [3.220391108732281, 5.880574208109178, 5.88059483485326, 10.695015549373402, 10.703628085272637, 10.74602559654693, 12.365315432961292, 12.42977000285817, 13.412583435713325, 13.62906574002197, 14.026615076737274, 16.66447956394288]
5 : [3.2203910320881675, 5.880559916842044, 5.880570294367756, 10.692557890164743, 10.695326232544247, 10.712028436979733, 12.334863812887155, 12.379428178940824, 13.399760629320207, 13.50718741606393, 13.909165088652987, 16.391078193349195]
6 : [3.2203910307938877, 5.880558531346219, 5.880569271711292, 10.690148741448663, 10.694658496396652, 10.699528894662288, 12.324262124271995, 12.355591119271866, 13.394789054604127, 13.45833131930213, 13.826480884169046, 16.00191051784656]
7 : [3.220391030769951, 5.880558425991055, 5.880569213416626, 10.688201871206639, 10.694580575766166, 10.69576099066912, 12.320782687228178, 12.34102797299317, 13.391306508159236, 13.438346613633467, 13.749378169767372, 15.513779093917922]
8 : [3.2203910307694934, 5.880558418312427, 5.880569209656789, 10.687056211430223, 10.69456559727016, 10.694758763426604, 12.319633665709452, 12.33135987188569, 13.388278509397068, 13.43012774377998, 13.666503235491154, 15.038678648114129]
9 : [3.2203910307694827, 5.880558417786918, 5.880569209414909, 10.686597154751967, 10.69446663643104, 10.69457030327261, 12.31924802294976, 12.3252882868672, 13.385179207473072, 13.426655484749658, 13.582277929379064, 14.688877898805584]
10 : [3.220391030769483, 5.880558417753172, 5.880569209399849, 10.686447459752136, 10.69438637736644, 10.694567747893926, 12.319117075880325, 12.321905845732527, 13.381103306583219, 13.425137659715302, 13.511735784499745, 14.47586404988724]
11 : [3.2203910307694765, 5.880558417751135, 5.88056920939898, 10.686405049580955, 10.694364165500167, 10.6945673683272, 12.319071998106335, 12.320286339121056, 13.3754653984621, 13.424457986936353, 13.466302119181274, 14.358017112009902]
12 : [3.2203910307694836, 5.8805584177510015, 5.880569209398923, 10.686394009952576, 10.694358423774517, 10.694567284811125, 12.31905598419932, 12.319599354997687, 13.369259836256912, 13.42414883466161, 13.44291891290163, 14.294003411351786]
13 : [3.220391030769481, 5.880558417750987, 5.880569209398913, 10.686391282688955, 10.694357004701157, 10.694567265288011, 12.319049665899053, 12.31933199952156, 13.364427422931762, 13.424005882025426, 13.432339996335699, 14.25909675965832]
14 : [3.2203910307694827, 5.88055841775096, 5.880569209398855, 10.686390627019076, 10.694356662637725, 10.694567260680607, 12.319046755143788, 12.319233293267049, 13.36153606817073, 13.423936582305881, 13.427676399466224, 14.239940081604882]
15 : [3.220391030769482, 5.88055841775147, 5.880569209399466, 10.686390471533734, 10.694356581244616, 10.694567259593102, 12.319045391247059, 12.31919785230786, 13.360036023999733, 13.42389646054184, 13.425594284368275, 14.229350710146019]
16 : [3.220391030769482, 5.880558417751095, 5.880569209399022, 10.68639043492443, 10.694356562006861, 10.694567259336205, 12.319044822900246, 12.319185266678153, 13.359310538010389, 13.423859971649673, 13.424665362108325, 14.223481074636084]
17 : [3.220391030769481, 5.880558417751015, 5.880569209398479, 10.686390426328362, 10.694356557473302, 10.694567259275406, 12.31904460769587, 12.319180807488495, 13.358970824789479, 13.423810502639068, 13.424271746952975, 14.220219849525547]
18 : [3.220391030769478, 5.88055841775053, 5.880569209399283, 10.686390424281889, 10.694356556398388, 10.69456725926115, 12.319044527722275, 12.319179197278057, 13.358809190395212, 13.423752951461172, 13.424119931834682, 14.218349557953244]
19 : [3.220391030769482, 5.880558417750857, 5.880569209399024, 10.686390423834565, 10.694356556156233, 10.694567259257983, 12.319044498568552, 12.319178648581387, 13.358735837966691, 13.423710411838206, 13.424066331786355, 14.2172719159992]

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

Eigenvalues
3.220391030769482
5.880558417750857
5.880569209399024
10.686390423834565
10.694356556156233
10.694567259257983
12.319044498568552
12.319178648581387
13.358735837966691
13.423710411838206
13.424066331786355
14.2172719159992

gfu = GridFunction(fes, multidim=len(evecs))
for i in range(len(evecs)):
    gfu.vecs[i].data = evecs[i]
Draw (Norm(gfu), mesh, "u")

BaseWebGuiScene