from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
from ngsolve.bem import *

keys: Acoustic Scattering, Brakhage-Werner formulation, Helmholtz equation, Dirichlet data

Helmholtz Problem

We consider a Helmholtz problem with a plane wave as initial condition

\begin{align} \left{ \begin{array}{rcl l} \Delta u + \kappa^2 u &= 0, \quad &\Omega \subset \mathbb R^3,,\ u_{in} &= e^{i\kappa x}.\end{array} \right. \end{align}

Combined field integral equations combine single and double layer integral operators, one simple option is the Brakhage-Werner formulation.

The solution is represented as

\[ u = (i \kappa S - D) \phi, \]

where \(\phi\) solve the boundary integral equation $\( \big( \tfrac{1}{2} + K + i \kappa V \big) \phi = u_{in} \qquad \text{on} \, \Gamma = \partial \Omega \)$

The single layer and double layer operator are defined using the fundamental solution of the Helmholtz operator,

\[ G(x,y) = \frac{e^{i\kappa \|x-y\|}}{4\pi\,\|x-y\|}. \]

The layer potentials have the same form as in the Laplace case, except that the Laplace kernel is replaced by the Helmholtz kernel \(G\).

kappa=10
order=4

screen = WorkPlane(Axes( (0,0,0), Z, X)).RectangleC(15,15).Face()

sp = Fuse(Sphere( (0,0,0), pi).faces)
screen.faces.name="screen"
sp.faces.name="sphere"
shape = Compound([screen,sp])

mesh = shape.GenerateMesh(maxh=5/kappa).Curve(order)
Draw (mesh);

fes_sphere = Compress(SurfaceL2(mesh, 
                                order=order, 
                                complex=True, 
                                definedon=mesh.Boundaries("sphere")))
u,v = fes_sphere.TnT()
fes_screen = Compress(SurfaceL2(mesh, order=order, 
                                dual_mapping=True, 
                                complex=True, 
                                definedon=mesh.Boundaries("screen")))
print ("ndof_sphere = ", fes_sphere.ndof, "ndof_screen =", fes_screen.ndof)

ndof_sphere =  17010 ndof_screen = 31110

with TaskManager():
    C = HelmholtzCF(u*ds("sphere"), kappa)*v*ds
    u,v  = fes_sphere.TnT()
    Id = BilinearForm(u*v*ds).Assemble()

lhs = 0.5 * Id.mat + C.mat
source = exp(1j * kappa * x) 
rhs = LinearForm(source*v*ds).Assemble()

gfu = GridFunction(fes_sphere)
pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()
with TaskManager():
    gfu.vec[:] = solvers.GMRes(A=lhs, b=rhs.vec, pre=pre, maxsteps=40, tol=1e-8, 
                               printrates=False)

Draw (gfu, order=5, min=-1, max=1);

Evaluation of the Solution

uscat = GridFunction(fes_screen)
with TaskManager():
    uscat.Set(1j*kappa*HelmholtzSL(u*ds("sphere"),kappa)(gfu) - HelmholtzDL(u*ds("sphere"), kappa)(gfu), \
              definedon=mesh.Boundaries("screen"))

print ("Scattered field")
Draw (uscat, mesh, min=-1,max=1, animate_complex=True, order=4);

Scattered field

uin = mesh.BoundaryCF( {"screen": source }, default=0)
print ("Total field")
Draw (uin+uscat, mesh, min=-1,max=1, animate_complex=True, order=4);

Total field