2D HCurl in 3D for surface PDE

NGSolve
1

Hello everone,

I am trying to simulate some nano-photonics device/optical guide. On the input surface (i.e. the boundary on which the light enters the device) the boundary conditions (optical modes) are the solution of a surface eigenvalue problem given by:

\nabla_{\parallel} \times (\nabla_{\parallel} \times E_{\parallel}) - \beta^2 \, \nabla_{\bot} E_{z} + (\beta^2 - k^2) E_{\parallel} = 0 \beta^2 (-\delta_{\parallel} E_{z} + div_{\parallel} E_{\parallel} - k^2 E_{z}) = 0

with f_{\parallel} being only the X and Y component of the respective vector/operator, \beta^2 the eigenvalue and E = (E_{\parallel}, E_z) the eigenvector. Here all components only depend on x and y but not on z.

I would like to directly compute the solution of the EVP on the respective boundary of the 3-dimensional mesh and then use E_{\parallel} for the boundary values of the main problem.

I looked up tutorial 6.1.2 for the definition of surface PDEs and I’m using the Arnoldi eigensolver accoridng to 1.7.1.

My main question is: Can I acces the 2-dimensional HCurl FE space (only dependend on x and y) for the E_{\parallel} component on a 3-dimensional mesh? Or is there an elegant way to force E_z to be 0?

Currently I’m using a compound FE space consisting of HCurl * H1 (since E_z has more regularity than E_{\parallel}).

This approach was presented by N. Lebbe page 6.

Thanks for any input

Greetings
Philipp

2

Hi Philipp,

you can define the restriction of an HCurl-space to a part of the boundary.
dofs are assigned to all edges of the mesh, but only dofs at the port-boundary are marked as freedofs.
You can build boundary-mass and boundary curl-stiffness matrices using the boundary trace operator.

portbnd = geo.Boundaries("bndname")
fes = HCurl(mesh, definedon=portbnd)
print (fes.FreeDofs())

u,v = fes.TnT()
mass = BilinearForm(u.Trace()*v.Trace()*ds(portbnd)).Assemble()
stiff = BilinearForm(curl(u).Trace()*curl(v).Trace()*ds(portbnd)).Assemble()
  • Joachim
3

Hello Joachim,

tanks for your response. I just noticed a small mistake in my code - fixing it already gave me the correct result.

Maybe one more question: Is it possible (if so how) to define a bilinear form which looks like the product of two integrals, i. e.

a(Test, Trial) = Int( Test * f )*dx * Int( Trial * g)*dx

for given f and g?

Is it as simple as " (ufdx) * (vgdx) " or is there a more specific way to prescribe which integration runs over which function?

Greetings,
Philipp

4

Hi Philipp,
a(u, v) = Int( u * f )dx * Int( v * g)dx
would lead to a completely dense matrix. A possibility is to introduce another unknown mu representing the integral Int( uf ), i.e. mu=Int( uf ), mu is just a number. Then the equivalent saddle point problem

fes1 = H1(mesh, order=1) fes2 = NumberSpace(mesh) fes = fes1*fes2 (u,mu), (v,lam) = fes.TnT() a( (u,mu), (v,lam) ) = mu*g*v*dx + (f*u-mu)*lam*dx

should work.

Best,
Michael

5

Hello again,

new problem, same topic: I’d like to extend this problem to multiple (given) functions f and g, i.e. list of functions (f_i, g_i)_i. At first I defined a separate NumberSpace for each f_i, yet this quickly became very cumbersome. Therefore, since the number of f_i will be computed as part of the problem, I’d like to make it a bit more adaptive.
Initially I thought I could simply make use of the “dim” flag during the definition of the NumberSpace and then take the product space of HCurl and NumberSpace. However, as soon as I set dim greater than 1, I get the error “Product space needs same dimension for all components”.

Is there an elegant way to work around this?

Greetings,
Philipp

6

Hi,
you can create a NumberSpace(mesh)**N, test and trialfunctions will then be Nx1 matrices:


from netgen.occ import unit_square
from ngsolve import *


mesh = Mesh(unit_square.GenerateMesh(maxh=1))

fes = H1(mesh) * NumberSpace(mesh)**5
(u, uN), (v,vN) = fes.TnT()


a = BilinearForm(fes)
a += u * v * dx
a += uN.trans * vN * dx

a.Assemble()
print(a.mat)

Best
Christopher