This page was generated from unit-1.7-helmholtz/helmholtz.ipynb.

1.7 Complex-valued waves

In NGSolve finite element spaces can be built and linear systems can be solved over the complex field. This tutorial shows how to compute the solution of the Helmholtz equation with impedance boundary conditions in complex arithmetic. The boundary value problem is to find u satisfying

Δuω2u=f in Ω

together with the impedance (outgoing) boundary condition

uniωu=0 on Ω

where i= 1j is the imaginary unit.

[1]:
import netgen.gui
from ngsolve import *
from netgen.geom2d import SplineGeometry
[2]:
# Geometry
geo = SplineGeometry()
geo.AddCircle((0.5, 0.5), 0.8,  bc="outer")
geo.AddRectangle((0.7, 0.3), (0.75, 0.7),
                 leftdomain=0, rightdomain=1, bc="scat")
mesh = Mesh(geo.GenerateMesh(maxh=0.05))

Declare a complex finite element space

[3]:
fes = H1(mesh, order=5, complex=True)
u, v = fes.TnT()
[4]:
# Wavenumber & source
omega = 100
pulse = 1e3*exp(-(100**2)*((x-0.5)*(x-0.5) + (y-0.5)*(y-0.5)))
Draw(pulse, mesh, 'pulse')

Forming the system

The weak form for uH1:

Ω[uˉvω2uˉv]dxiωΩuˉvds=Ωfˉv

for all v in H1.

[5]:
# Forms
a = BilinearForm(fes)
a += grad(u)*grad(v)*dx - omega**2*u*v*dx
a += -omega*1j*u*v * ds("outer")
a.Assemble()

f = LinearForm(fes)
f += pulse * v * dx
f.Assemble()

Solve

[6]:
gfu = GridFunction(fes, name="u")
gfu.vec.data = a.mat.Inverse() * f.vec
Draw(gfu)

Explore the GUI’s menu options in Visual tab:

  • Increase subdivions

  • Real and imaginary parts

  • View absolute value

  • Turn off Autoscale

  • Turn on Deformation

  • Turn on Periodic Animation