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

1.7 Complex-valued waves

In NGSolve finite element spaces can be built over the complex field and the resulting complex matrix systems can be solved. 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

\[-\Delta u - \omega^2 u = f\qquad \text{ in } \Omega\]

together with the impedance (outgoing) boundary condition

\[\frac{\partial u }{ \partial n} - i \omega u = 0 \quad \text{ on } \partial \Omega\]

where \(i =\) 1j is the imaginary unit.

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *
from netgen.geom2d import SplineGeometry

air = Circle((0.5, 0.5), 0.8).Face()
air.edges.name = 'outer'
scatterer = MoveTo(0.7, 0.3).Rectangle(0.05, 0.4).Face()
scatterer.edges.name = 'scat'
geo = OCCGeometry(air - scatterer, dim=2)
mesh = Mesh(geo.GenerateMesh(maxh=0.05))
Draw(mesh);

Declare a complex-valued finite element space

fes = H1(mesh, order=5, complex=True)
u, v = fes.TnT()

omega = 100
pulse = 5e4*exp(-(40**2)*((x-0.5)*(x-0.5) + (y-0.5)*(y-0.5)))
Draw(pulse, mesh, order=3)

BaseWebGuiScene

Forming the system

The weak form for \(u \in H^1\):

\[\int_\Omega\big[ \nabla u \cdot \nabla \bar v - \omega^2 u \bar v \big] \, dx - i \,\omega\, \int_{\partial \Omega} u \bar v \, ds = \int_{\Omega} f \bar v\]

for all \(v\) in \(H^1\).

# 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(pulse * v * dx).Assemble();

Solve

gfu = GridFunction(fes, name="u")
gfu.vec.data = a.mat.Inverse() * f.vec
Draw(gfu, mesh, min=-1, max=1, order=3, animate_complex=True);

Open controls and explore webgui’s further viewing options such as viewing real and imaginary parts, viewing absolute value, and viewing with deformation turned on.