Layer potentials

34. Layer potentials#

Let \(G(x,y) = \frac{1}{4 \pi} \frac{\exp (i k |x-y|)}{|x-y|}\) be Green’s function for the Helmholtz equation. For a surface \(\Gamma\), and a scalar function \(\rho\) defined on \(\Gamma\), we define the single layer potential as

\[ (V \rho) (x) = \int_{\Gamma} G(x,y) \rho(y) dy \]
from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw

# the boundary Gamma
face = WorkPlane(Axes((0,0,0), -Y, Z)).RectangleC(1,1).Face()
mesh = Mesh(OCCGeometry(face).GenerateMesh(maxh=0.1))
Draw (mesh);

# the visualization mesh
visplane = WorkPlane(Axes((0,0,0), Z, X)).RectangleC(5,5).Face()
vismesh = Mesh(OCCGeometry(visplane).GenerateMesh(maxh=0.1))
Draw (vismesh);

from ngsolve.bem import SingularMLMultiPoleCF, RegularMLMultiPoleCF
from ngsolve.bla import Vec3D

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Cell In[2], line 1
----> 1 from ngsolve.bem import SingularMLMultiPoleCF, RegularMLMultiPoleCF
      2 from ngsolve.bla import Vec3D

ImportError: cannot import name 'SingularMLMultiPoleCF' from 'ngsolve.bem' (/Applications/Netgen.app/Contents/Resources/lib/python3.14/site-packages/ngsolve/ngslib.so)

We evaluate the single layer integral using numerical integration on the surface mesh. Thus, we get a sum of many Green’s functions, which is compressed using a multilevel-multipole.

kappa = 3*pi
mp = SingularMLMultiPoleCF(Vec3D(0,0,0), r=3, kappa=kappa)

ir = IntegrationRule(TRIG,20)
pnts = mesh.MapToAllElements(ir, BND)
# vals = (20*x)(pnts).flatten()
vals = CF(1)(pnts).flatten()

# find the integration weights:  sqrt(F^T F)*weight_ref
F = specialcf.JacobianMatrix(3,2)
weightCF = sqrt(Det (F.trans*F))
weights = weightCF(pnts).flatten()
for j in range(len(ir)):
    weights[j::len(ir)] *= ir[j].weight
print ("number of source points: ", len(pnts))
for p,vi,wi in zip(pnts, vals, weights):
    mp.mlmp.AddCharge (Vec3D(x(p), y(p), z(p)), vi*wi)

mp.mlmp.Calc()

Draw (mp, vismesh, min=-0.02,max=0.02, animate_complex=True, order=2);

We see that the single layer potential is continuous across the surface \(\Gamma\), but has a kink at \(\Gamma\). This shows that the normal derivative is jumping (exactly by \(\rho\)).

regmp = RegularMLMultiPoleCF(mp, Vec3D(0,0,0.00),r=5)

Draw (regmp, vismesh, min=-0.02,max=0.02, animate_complex=True, order=2);

Draw (mp.real-regmp.real, vismesh, min=-1e-5,max=1e-5, animate_complex=False, order=2);

34.1. Double layer potential#

the double layer potential is

\[ (K \rho) (x) = \int_{\Gamma} n_y \cdot \frac{\partial G}{\partial y}(x,y) \rho(y) dy \]

the name comes from adding a charge density \(\tfrac{1}{2 \varepsilon} \rho\) at the offset surface \(\Gamma + \varepsilon n\), and a second charge density \(\tfrac{-1}{2 \varepsilon} \rho\) at \(\Gamma - \varepsilon n\), and passing to the limit, i.e.

\[ (K \rho) (x) = \lim_{\varepsilon \rightarrow 0} \int_{\Gamma} \frac{1}{2 \varepsilon } (G(x,y+\varepsilon n) - G(x,y-\varepsilon n) \big) \rho(y) dy, \]

This pair of charges defines a charge dipole in normal direction.

kappa = pi
mp = SingularMLMultiPoleCF(Vec3D(0,0,0), 2, 20, kappa)

ir = IntegrationRule(TRIG,5)
pnts = mesh.MapToAllElements(ir, BND)
vals = (20*x)(pnts).flatten()

F = specialcf.JacobianMatrix(3,2)
weightCF = sqrt(Det (F.trans*F))
weights = weightCF(pnts).flatten()
for j in range(len(ir)):
    weights[j::len(ir)] *= ir[j].weight

for p,vi,wi in zip(pnts, vals, weights):
    mp.mlmp.AddDipole (Vec3D(x(p), y(p), z(p)), Vec3D(0,1,0), vi*wi)
mp.mlmp.Calc()

Draw (mp.real, vismesh, min=-1,max=1, animate_complex=True, order=2)
Draw (mp, vismesh, min=-1,max=1, animate_complex=True, order=2);

Now we see that the double layer potential is discontinuous (with jump exactly equal to \(\rho\)), and the normal derivative is continuous.

These layer potentials are the foundation for the boundary element method, see ngbem

Contents