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

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

number of source points:  27830

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

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[8], line 2
      1 kappa = pi
----> 2 mp = SingularMLMultiPoleCF(Vec3D(0,0,0), 2, 20, kappa)
      4 ir = IntegrationRule(TRIG,5)
      5 pnts = mesh.MapToAllElements(ir, BND)

TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. ngsolve.bem.SingularMLMultiPoleCF(center: ngsolve.bla.Vec3D, r: float, kappa: float)

Invoked with: <ngsolve.bla.Vec3D object at 0x10f2da430>, 2, 20, 3.141592653589793

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