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)

 Start Findpoints
 main-solids: 1
 Found points 14
 Analyze spec points
 Find edges
 21 edges found
 Check intersecting edges
 CalcLocalH: 53 Points 0 Elements 0 Surface Elements
 Start Findpoints
 Analyze spec points
 Find edges
 21 edges found
 Check intersecting edges
 Start Findpoints
 Analyze spec points
 Find edges
 21 edges found
 Check intersecting edges
 Surface 1 / 9
 load internal triangle rules
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 34 elements, 63 points
 Surface 2 / 9
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 36 elements, 74 points
 Surface 3 / 9
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 20 elements, 77 points
 Surface 4 / 9
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 6 elements, 77 points
 Surface 5 / 9
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 6 elements, 77 points
 Surface 6 / 9
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 26 elements, 83 points
 Surface 7 / 9
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 22 elements, 87 points
 Surface 8 / 9
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 6 elements, 87 points
 Surface 9 / 9
 Surface meshing done
 0 illegal triangles
 Optimize Surface
 Edgeswapping, topological
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Combine improve
 Smoothing
 38 elements, 99 points
 SplitSeparateFaces
 Check subdomain 1 / 1

 Meshing subdomain 1 of 1
 Use internal rules
 Delaunay meshing
 number of points: 94
 blockfill local h
 number of points: 116
 Points: 116
 Elements: 599
 Tree data entries per element: 1.5025
 Tree nodes per element: 0.0283806
 SwapImprove
 SwapImprove
 SwapImprove
 SwapImprove
 0 degenerated elements removed
 Remove intersecting
 Remove outer
 116 points, 364 elements
 start tetmeshing
 Use internal rules
 Success !
 116 points, 364 elements
 Remove Illegal Elements
 Volume Optimization
 CombineImprove
 ImproveMesh
 SplitImprove
 ImproveMesh
 SwapImprove
 SwapImprove2
 ImproveMesh
 CombineImprove
 ImproveMesh
 SplitImprove
 ImproveMesh
 SwapImprove
 SwapImprove2
 ImproveMesh
 CombineImprove
 ImproveMesh
 SplitImprove
 ImproveMesh
 SwapImprove
 SwapImprove2
 ImproveMesh
 Update mesh topology
 Update clusters
 HP Refinement called, levels = 2
 Curve elements, order = 1
 Update mesh topology
 Update clusters

 undefined elements update classification: 0
 non-implemented in update classification: 0
  Start new hp-refinement: step 1
  in HP-REFINEMENT with fac1 0.2
 0 wrong prisms, 108 right prisms
 Update mesh topology
 Update clusters

 undefined elements update classification: 0
 non-implemented in update classification: 0
  Start new hp-refinement: step 2
  in HP-REFINEMENT with fac1 0.2
 0 wrong prisms, 191 right prisms
 Update mesh topology
 Update clusters

 PrepareElements no more to do for actual refinement 3
  HP-Refinement done with 2 refinement steps.

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

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)

assemble VOL element 507/507
assemble VOL element 507/507
assemble VOL element 507/507
Update Direct Solver Preconditioner0 : [7.112482703170652, 24.49913826316374, 34.081714586920306, 48.076314293830656, 49.554152514734206, 59.64161425548686, 70.7657725304168, 73.6828103277423, 83.38920909147443, 86.35756158441711, 100.14220886655038, 109.54562388708857]
1 : [3.277339714274206, 6.152771237150206, 6.950869243937089, 12.053801914592805, 12.708736618944144, 13.477940688197997, 15.293505396756348, 18.27682544912982, 20.060299477186646, 21.542971881776808, 24.628694721748868, 26.556918968139023]
2 : [3.221300820313879, 5.89169992041499, 5.934478171707845, 10.932674455719475, 11.04329310387672, 11.251732279520676, 12.847355533150637, 13.41228481768711, 14.280538541766198, 15.096436854003052, 18.177198038610193, 19.463982312757086]
3 : [3.220364279667241, 5.881034914886239, 5.883284541316069, 10.74271084460103, 10.772577499013021, 10.805814060098333, 12.433415601613909, 12.626333608583465, 13.553558380777787, 13.88718437816403, 14.906344142559435, 16.369459567426603]
4 : [3.2203515658745094, 5.88060027948238, 5.88070948287469, 10.70200282815241, 10.705640219444739, 10.720602910103525, 12.341702361755486, 12.432192846316974, 13.419859128485838, 13.560571612462182, 13.983055451829651, 15.128880724399579]
5 : [3.2203513836070115, 5.88057306792141, 5.880584761963929, 10.689208297756151, 10.696880267259775, 10.700407379685608, 12.324030806037033, 12.362905439448392, 13.378420775438661, 13.469764303888608, 13.650872527012968, 14.625852190860659]
6 : [3.2203513809705155, 5.880566770170538, 5.880583589320373, 10.686867334566598, 10.694975095620402, 10.69579383858378, 12.32047750611579, 12.335703681081197, 13.366363383473818, 13.440243787013191, 13.520422128784409, 14.416333951420372]
7 : [3.2203513809308606, 5.880566438492955, 5.880583547755376, 10.686442012678839, 10.694562929430623, 10.694732800716208, 12.319728749231388, 12.325343928945431, 13.362218671009055, 13.429993934517343, 13.466427907324295, 14.319108277303407]
8 : [3.2203513809302557, 5.8805664213755815, 5.880583545966608, 10.686360701588715, 10.694470431759452, 10.694492295113408, 12.319559845637247, 12.321552076351578, 13.360423060310346, 13.426329925519477, 13.443238237836365, 14.270719631053002]
9 : [3.220351380930242, 5.880566420481258, 5.880583545885012, 10.686344212782895, 10.69443020324334, 10.694456900275156, 12.31951139413323, 12.320202944356762, 13.3595933540199, 13.42494174140927, 13.432988256639932, 14.245533540424974]
10 : [3.220351380930246, 5.880566420434145, 5.880583545881084, 10.686340713942192, 10.694417952532412, 10.694451905280017, 12.319474931444388, 12.319748801518585, 13.35921388188637, 13.424381046590067, 13.428385448345374, 14.232141909619294]
11 : [3.2203513809302446, 5.880566420431613, 5.880583545880891, 10.686339948145825, 10.69441507597037, 10.694450828892714, 12.31942801279143, 12.319623780527504, 13.359042763157476, 13.424142667179204, 13.426298666875873, 14.224929431373733]
12 : [3.220351380930243, 5.880566420431497, 5.8805835458809, 10.686339777059407, 10.69441440635287, 10.694450589270764, 12.319397582769634, 12.319593962181761, 13.358966069728613, 13.424037518014064, 13.425348689717026, 14.221022053180887]
13 : [3.220351380930242, 5.88056642043144, 5.880583545880865, 10.686339738325252, 10.694414250334516, 10.69445053545169, 12.319384870687289, 12.319585552896752, 13.35893171322995, 13.42398986641364, 13.424915321496572, 14.21889640913101]
14 : [3.2203513809302433, 5.880566420431504, 5.880583545880838, 10.686339729479636, 10.694414213918167, 10.694450523324027, 12.319380171763806, 12.319582838595664, 13.358916295532806, 13.423967923949402, 13.424717428147954, 14.217738190600148]
15 : [3.2203513809302464, 5.880566420431519, 5.8805835458809295, 10.686339727447123, 10.694414205402524, 10.694450520580377, 12.319378492483237, 12.319581913181503, 13.358909360755549, 13.423957773201828, 13.424626874190206, 14.217105171630466]
16 : [3.2203513809302438, 5.880566420431835, 5.880583545881193, 10.686339726978703, 10.694414203409618, 10.694450519959691, 12.319377898805426, 12.319581591391572, 13.358906237817662, 13.423953089687725, 13.424585428160695, 14.21675966328625]
17 : [3.220351380930244, 5.880566420431147, 5.880583545880826, 10.686339726870235, 10.69441420294249, 10.6944505198187, 12.319377689365034, 12.319581478565155, 13.358904828184071, 13.423950933617101, 13.424566403621819, 14.216570309559769]
18 : [3.220351380930246, 5.880566420431275, 5.880583545881191, 10.68633972684527, 10.694414202826312, 10.694450519786388, 12.319377613182619, 12.319581438214547, 13.358904119426127, 13.423949876443686, 13.424556965538821, 14.216463604638054]
19 : [3.2203513809302455, 5.8805664204317765, 5.880583545880759, 10.686339726839277, 10.694414202804818, 10.694450519778577, 12.319377587958384, 12.319581424401626, 13.358903848449758, 13.423949430465504, 13.424552894680291, 14.216396565890442]

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

Eigenvalues
3.2203513809302455
5.8805664204317765
5.880583545880759
10.686339726839277
10.694414202804818
10.694450519778577
12.319377587958384
12.319581424401626
13.358903848449758
13.423949430465504
13.424552894680291
14.216396565890442

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

BaseWebGuiScene