Forum Message



We have moved the forum to . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.


The forum is in read only mode.

2D HCurl in 3D for surface PDE

1 year 3 months ago - 1 year 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


Last edit: 1 year 3 months ago by PKonig.
1 year 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
1 year 3 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?

1 year 3 months ago - 1 year 3 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 
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.

Last edit: 1 year 3 months ago by mneunteufel.
9 months 3 weeks ago - 9 months 3 weeks 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?


Last edit: 9 months 3 weeks ago by PKonig.
9 months 3 weeks ago #4662 by christopher
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)

Time to create page: 0.168 seconds