4.1. Hcurl and related problems

NGSolve contains a variety of finite elements for apporximating different function spaces:

• H1: Lagrange elements approximating functions in \(H^1\)

• HCurl: piecewise Nedelec elements approximating functions in \(H(curl)\)

• HDiv: piecewise Brezzi-Douglas-Marini elements approximating functions in \(H(div)\)

One can loop over the basis function and observe that they are build in a hierarchical way.

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

mesh = Mesh(unit_square.GenerateMesh(maxh=0.5,quad_dominated=True) )

order = 3
fes = H1(mesh, order=order)
gfu = GridFunction(fes )
gfut = GridFunction(fes , multidim=0)




for i in range(len(gfu.vec)):

    gfu.vec[:] = 0
    gfu.vec[i] = 1

    gfut.AddMultiDimComponent(gfu.vec)

settings = {"camera": {"transformations": [
        {"type": "rotateX", "angle": -45}]}}

Draw(gfut,deformation=True, animate=True ,autoscale = True , min=-0.03, max=0.03, order = 3, settings=settings);

order = 2
fes = HCurl(mesh, order=order)
gfu = GridFunction(fes)
gfut = GridFunction(fes, multidim=0)


for i in range(len(gfu.vec)):

    gfu.vec[:] = 0
    gfu.vec[i] = 1

    gfut.AddMultiDimComponent(gfu.vec)

Draw(gfut, deformation=False, animate=True, autoscale = True , order = 3,  vectors = { "grid_size":30}, settings=settings);

order = 2
fes = HDiv(mesh, order=order)
gfu = GridFunction(fes)
gfut = GridFunction(fes, multidim=0)


for i in range(len(gfu.vec)):

    gfu.vec[:] = 0
    gfu.vec[i] = 1

    gfut.AddMultiDimComponent(gfu.vec)

Draw(gfut, deformation=False, animate=True, autoscale = True , order = 3,  vectors = { "grid_size": 30}, settings=settings);

4.1.1. The eigenvalue problem on the torus:

For simply connected domains the kernel of the curl is given by the image of the gradient.

In case of a torus the kernel is slightly larger than the image of the gradient.

problem: find \(u \in H(curl)\) such that

\[\begin{align*} \text{curl} \,u &= 0 \\ u &\neq \nabla \phi \quad \forall \phi \in H^1 \end{align*}\]

To do so we solve the eigenvalue problem:
Find \(\lambda \in \mathbb{R}\) and \(u \in H(curl)\) such that

\[\begin{align*} \int_{\Omega} \text{curl} u \cdot \text{curl} v \,dx = \lambda \int_{\Omega} u \cdot v \,dx\\ u \cdot \nabla \phi =0 \quad \forall \phi \in H^1 \end{align*}\]

Start creating the (double) torus

circ = Circle( Pnt(2,0), 1.5).Face()
axis = Axis(Pnt(0,0,0), Dir(0,1,0))#, Dir(1,0,0))
angle = 360
torus = circ.Revolve( axis, angle)
torus2 = torus.Move( (4,0,0), )

geo = torus.Rotate(Axis((0, 0, 0), Y), 180) + torus2


mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=3)).Curve(4)


settings = {"camera": {"transformations": [ {"type": "rotateX", "angle": 45}]}}
Draw(mesh, settings=settings);

I will not go too much into detail here, to have a better understanding of the problem I recommend to read the Phd thesis of Sabine Zaglmayr, p 145-150.

# SetHeapSize(100*1000*1000)

fes = HCurl(mesh, order=2)

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 = Preconditioner(apre, "direct", inverse="sparsecholesky")

ndof = 34236

The idea is to use the PINVIT (preconditioned inverse iteration) method and let the preconditioner take care of the kernel of the curl.

with TaskManager():
    a.Assemble()
    m.Assemble()
    apre.Assemble() # when we assemble the BilinearForm , it assembles the matrix as well

    ### Advanced part of the code ###
    
    gradmat, fesh1 = fes.CreateGradient() # create gradient from H1 to HCurl
    gradmattrans = gradmat.CreateTranspose() 
    math1 = gradmattrans @ m.mat @ gradmat   # create projection onto grad H1
    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 : [0.11627436414698765, 0.20132032534854807, 2.9355567270685587, 3.5485560750946505, 4.1157422825294265, 4.487444737796731, 5.051517140116547, 5.480374752642389, 5.685824567166254, 6.5737473968640545, 7.231290902062009, 8.093067536090913]
1 : [0.001976658184892824, 0.002867832192860354, 1.3940862055353267, 1.6989243173789697, 1.7380215916491915, 1.848342430242387, 1.9699926174297602, 2.0925881047666612, 2.1484074889048626, 2.268114424359725, 2.82030667932415, 3.5806936071283286]

2 : [4.4134861095586106e-05, 8.595650706102788e-05, 1.1787001043817897, 1.5486302541360606, 1.6016621601680097, 1.635197857232276, 1.7229356830049554, 1.7333854046237276, 1.7628838377036253, 1.9474806304554806, 2.0092278933702135, 2.592519077397292]

3 : [1.155591801368637e-06, 2.6613274266178428e-06, 1.1257382463867385, 1.5160668871845888, 1.5393129218782737, 1.6066841179250204, 1.6355323478101835, 1.6401282824594212, 1.7117621769493931, 1.7780435377434, 1.8877160604113874, 2.1563335672541086]

4 : [3.246670949962612e-08, 8.719742240260279e-08, 1.1128056278744467, 1.5054753172597242, 1.511484586383966, 1.5937721508577936, 1.603772040420154, 1.6169700979350004, 1.665138050530807, 1.7152571084367232, 1.8646841118590656, 2.004501160247964]
5 : [9.469013804447221e-10, 3.0758044532929223e-09, 1.109732140378982, 1.4966820213566137, 1.5047241061758674, 1.5863816439699074, 1.5931123540910865, 1.6081028093305505, 1.6318688470982605, 1.7087794879736475, 1.855085234116962, 1.9442539976138666]

6 : [2.872961427292602e-11, 1.1632997475915567e-10, 1.1090008926198065, 1.4921546312362126, 1.503325023050903, 1.5828219423758458, 1.5891855338428218, 1.6007320431813685, 1.6222944864465376, 1.7067299978791786, 1.8504958539340899, 1.9153260136859174]

7 : [8.919841363371992e-13, 4.570286247564072e-12, 1.108821152823768, 1.4904413302633441, 1.502836696654507, 1.5813155258843403, 1.5873530480698326, 1.5963968961681763, 1.6196019779304074, 1.705823543098469, 1.8480529818488955, 1.8996869296474184]

8 : [2.8353522680401584e-14, 1.8485772942017741e-13, 1.1087754366015317, 1.4898176570974577, 1.5026489218384365, 1.5806989296861294, 1.5863834874230116, 1.5944032724106316, 1.6186324404250878, 1.7053688918318073, 1.8466516790178051, 1.8905906503249814]

9 : [1.0325956844241659e-15, 8.255414586861689e-15, 1.108763300046072, 1.4895932102206912, 1.5025726771678154, 1.5804457427689422, 1.5858794849615254, 1.5935308598472842, 1.618226927078801, 1.7051219299399911, 1.8457997664496306, 1.8850424267265051]
10 : [1.3928343015317575e-16, 2.9334937384032954e-15, 1.1087599437195974, 1.4895127763893232, 1.5025400485809348, 1.5803407973122694, 1.5856298162951814, 1.593146824583838, 1.6180459613869433, 1.704978399284936, 1.8452528184399652, 1.8815684070186671]
11 : [2.1487855781682066e-16, 7.441432770721513e-16, 1.1087589772141757, 1.4894839318851538, 1.5025252633634119, 1.5802968106339756, 1.5855094814899031, 1.5929751772193899, 1.6179627365703393, 1.7048894394823713, 1.8448776734791545, 1.8793543983307388]
12 : [-1.4118101015686767e-17, 2.183595438311754e-15, 1.1087586893839305, 1.489473562732161, 1.5025181668468652, 1.5802781163118491, 1.5854519690281812, 1.5928975590380166, 1.6179239353117687, 1.7048308791545987, 1.8446002501345171, 1.8779294258725625]
13 : [-9.25955573628027e-16, -1.0051212105443186e-16, 1.1087586012116493, 1.4894698194534126, 1.5025145730778047, 1.5802700131901366, 1.5854243144339952, 1.5928621441590323, 1.6179056968114827, 1.7047902233118886, 1.8443785993298565, 1.8770049127020612]

14 : [-1.0090112210878777e-16, 3.7026092041992674e-16, 1.1087585736309447, 1.4894684616869445, 1.5025126691186914, 1.5802664055850582, 1.585410806452539, 1.5928458711672537, 1.617897079794346, 1.7047607600113017, 1.8441892618330626, 1.8764018519526287]
15 : [-6.545296211549369e-15, -1.4462100321250894e-16, 1.1087585648645455, 1.4894679664463537, 1.502511623901795, 1.5802647413364728, 1.5854040476047806, 1.5928383402259163, 1.6178929890637015, 1.7047386961374444, 1.8440189807745735, 1.8760063750354394]
16 : [1.5421538959581364e-16, 1.778127694224692e-15, 1.1087585620459326, 1.4894677844890167, 1.502511035025646, 1.580263937961865, 1.5854005515776803, 1.5928348249718596, 1.617891033413385, 1.7047217877120702, 1.8438605401730843, 1.8757449793957748]
17 : [2.9361703190135946e-16, 2.6656751113160524e-15, 1.1087585611332313, 1.489467717486013, 1.5025106967401398, 1.580263530877109, 1.585398686595118, 1.5928331705256666, 1.6178900977935726, 1.7047086014526605, 1.8437093482614275, 1.8755722224803406]
18 : [4.0717188418338763e-16, 2.131636787747978e-14, 1.1087585608360038, 1.4894676926410053, 1.502510499807289, 1.580263313311323, 1.5853976539013208, 1.592832381801781, 1.6178896468406714, 1.7046981843315816, 1.843562505440353, 1.875457494498738]

19 : [1.8387075734669504e-16, 1.4211226319783884e-14, 1.1087585607370307, 1.4894676833500304, 1.5025103803044517, 1.580263186191435, 1.5853970392552676, 1.5928319970162923, 1.6178894263708627, 1.7046895374834352, 1.8434115231235144, 1.8753801835934505]

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

Eigenvalues
1.8387075734669504e-16
1.4211226319783884e-14
1.1087585607370307
1.4894676833500304
1.5025103803044517
1.580263186191435
1.5853970392552676
1.5928319970162923
1.6178894263708627
1.7046895374834352
1.8434115231235144
1.8753801835934505

gfu = GridFunction(fes, multidim=len(evecs))
for i in range(len(evecs)):
    gfu.vecs[i].data = evecs[i]
Draw(gfu, mesh, order=2,clipping = {"y":-1, "z":0} ,vectors = { "grid_size": 30},settings=settings);