This page was generated from unit-11.3-bem-Helmholtz/BrakhageWerner.ipynb.

11.3.1 Helmholtz solver using Brakhage-Werner formulation

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 = (K + i \kappa V) \phi,\]

where \(\phi\) solve the boundary integral equation

\[\big( \tfrac{1}{2} + K + i \kappa V \big) \phi = u_{in} \qquad \text{on} \, \Gamma\]

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

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=1).Curve(order)
Draw (mesh);

fes_sphere = Compress(SurfaceL2(mesh, order=order, complex=True, definedon=mesh.Boundaries("sphere")))
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 =  4410 ndof_screen = 8010

kappa = 5

with TaskManager():
    V = HelmholtzSingleLayerPotentialOperator(fes_sphere, fes_sphere, kappa=kappa, intorder=10)
    K = HelmholtzDoubleLayerPotentialOperator(fes_sphere, fes_sphere, kappa=kappa, intorder=10)
    C = HelmholtzCombinedFieldOperator(fes_sphere, fes_sphere, kappa=kappa, intorder=10)
    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=10)

GMRes iteration 1, residual = 30.56280973696015
GMRes iteration 2, residual = 10.458547384840394
GMRes iteration 3, residual = 4.089611474126906
GMRes iteration 4, residual = 1.697464806944205
GMRes iteration 5, residual = 0.7290658300930473
GMRes iteration 6, residual = 0.34196426804845403
GMRes iteration 7, residual = 0.14934711993691918
GMRes iteration 8, residual = 0.06730848144077056
GMRes iteration 9, residual = 0.03243615905756445
GMRes iteration 10, residual = 0.015951069339732876
WARNING: GMRes did not converge to TOL

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

prostprocessing on screen

uscat = GridFunction(fes_screen)
with TaskManager():
    uscat.Set(1j*kappa*V.GetPotential(gfu)-K.GetPotential(gfu) , definedon=mesh.Boundaries("screen"))

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

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

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

Scattered field only: