Helmholtz equation with robin condition

Dear NGSolve developers,

I am a newbie to NGSolve, apologies if this sounds like a silly question.

I am trying to solve Helmholtz equation (\Delta u + \omega^2 u = 0 in \Omega with robin bc: \nabla u \cdot \mathbf{n} + i \omega u = g on \Gamma). I set my exact solution to be a plane propagating wave travelling in +ve x direction, i.e., u = \exp (i \omega x). This is imposed via robin bc, so the weak form reads:
find u \in H^1 (\Omega), s.t.
\int_\Omega \nabla u \cdot \nabla v - \omega^2 u v d\mathbf{x} + i\omega \int_{\Gamma} u v = \int_{\Gamma} g v
for all v in H^1(\Omega)

The NGSolve example code I am following is here. That example code has a non-zero source term and a zero right hand side for the robin bc. While I have opposite of that, i.e., a zero source term and a non-zero term in right hand side in robin bc.

Here is my python script:

[code]# following: 1.7 Complex-valued waves — NGS-Py 6.2.2302-87-ga5a5eff3b documentation
from ngsolve import *
from netgen.geom2d import unit_square

ngsglobals.msg_level = 1

generate a triangular mesh of mesh-size 0.2

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

set wavenumber

waveno = 1.0

H1-conforming finite element space

fes = H1(mesh, order=1, complex=True)

define trial- and test-functions

u = fes.TrialFunction()
v = fes.TestFunction()

Forms

a = BilinearForm(fes)
a += SymbolicBFI(grad(u)grad(v)-waveno**2uv)
a += SymbolicBFI(waveno
1juv, definedon=mesh.Boundaries(“outer”))
a.Assemble()

is this correct way to get normals

normals = specialcf.normal(mesh.dim)

rhs

f = LinearForm(fes)
f += SymbolicLFI((normals[0]+1.0)(exp(1jwaveno*x)1jwaveno) * v, definedon=mesh.Boundaries(“outer”))
f.Assemble()
#solve
gfu = GridFunction(fes, name=“u”)
gfu.vec.data = a.mat.Inverse() * f.vec
Draw(gfu)
[/code]

I must be doing something stupid when I assemble my rhs linear form as f.vec are all zeros. Perhaps the way I use normal vector is not correct or I don’t understand how to assemble the rhs when source term is not present. What am doing wrong here? Many thanks.

Hi, from a first glance I can see that you use the wrong name for the boundary conditions. The unit square has bcnames “top”, “bottom”, “right” and “left” and not “outer”. This may be the reason why your terms do not do anything. Try replacing “outer” with “top|left|right|bottom”.
You can get the available boundary names of a mesh with

print(mesh.GetBoundaries())

Best
Christopher