2D HCurl in 3D for surface PDE

More
2 years 3 months ago - 2 years 3 months ago #4489 by PKonig
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




 
Last edit: 2 years 3 months ago by PKonig.
More
2 years 3 months ago #4490 by joachim
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
 
More
2 years 2 months ago #4491 by PKonig
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 " (u*f*dx) * (v*g*dx) " or is there a more specific way to prescribe which integration runs over which function?

Greetings,
Philipp
More
2 years 2 months ago - 2 years 2 months ago #4493 by mneunteufel
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( u*f ), i.e. mu=Int( u*f ), mu is just a number. Then the equivalent saddle point problem 
Code:
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
Last edit: 2 years 2 months ago by mneunteufel.
More
1 year 9 months ago - 1 year 9 months ago #4661 by PKonig
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

 
Last edit: 1 year 9 months ago by PKonig.
More
1 year 9 months ago #4662 by christopher
Hi,
you can create a NumberSpace(mesh)**N, test and trialfunctions will then be Nx1 matrices:
Code:
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
Time to create page: 0.113 seconds