This page was generated from unit-11.1-fmm/layerpotentials.ipynb.

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

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