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

\[-\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.geom2d import SplineGeometry

# 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))

 Generate Mesh from spline geometry
 Boundary mesh done, np = 118
 CalcLocalH: 118 Points 0 Elements 0 Surface Elements
 Meshing domain 1 / 1
 load internal triangle rules
 Surface meshing done
 Edgeswapping, topological
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Update mesh topology
 Update clusters

Declare a complex finite element space

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

# 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')

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(fes)
f += pulse * v * dx
f.Assemble()

assemble VOL element 1788/1788
assemble BND element 118/118

<ngsolve.comp.LinearForm at 0x7f0ab7515f30>

assemble VOL element 1788/1788

Solve

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

call pardiso ... done

BaseWebGuiScene

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