- Thank you received: 0
2D HCurl in 3D for surface PDE
2 years 4 months ago - 2 years 4 months ago #4489
by PKonig
2D HCurl in 3D for surface PDE was created 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
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 4 months ago by PKonig.
2 years 4 months ago #4490
by joachim
Replied by joachim on topic 2D HCurl in 3D for surface PDE
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
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
2 years 4 months ago #4491
by PKonig
Replied by PKonig on topic 2D HCurl in 3D for surface PDE
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
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
- mneunteufel
- Offline
- Premium Member
Less
More
- Thank you received: 59
2 years 3 months ago - 2 years 3 months ago #4493
by mneunteufel
Replied by mneunteufel on topic 2D HCurl in 3D for surface PDE
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
should work.
Best,
Michael
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 3 months ago by mneunteufel.
1 year 10 months ago - 1 year 10 months ago #4661
by PKonig
Replied by PKonig on topic 2D HCurl in 3D for surface PDE
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
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 10 months ago by PKonig.
- christopher
- Offline
- Administrator
Less
More
- Thank you received: 101
1 year 10 months ago #4662
by christopher
Replied by christopher on topic 2D HCurl in 3D for surface PDE
Hi,
you can create a NumberSpace(mesh)**N, test and trialfunctions will then be Nx1 matrices:
Best
Christopher
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.110 seconds