from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
from ngsolve.bem import *
from ngsolve import Projector, Preconditioner
from ngsolve.krylovspace import CG

keys: electrostatics, homogeneous Dirichlet bvp, double layer potential \(K\) and single layer potential \(V\)

Dirichlet Laplace Direct Method

Consider the Dirichlet boundary value problem

\[\begin{split} \left\{ \begin{array}{rcl l} \Delta u &=& 0, \quad &\Omega \subset \mathbb R^3\,,\\ \gamma_0 u&=& m, \quad &\Gamma = \partial \Omega\,.\end{array} \right. \end{split}\]

Let us choose the following ansatz for the solution \(u\in H^1(\Omega)\) (direct ansatz)

\[ u(x) = \underbrace{ \int\limits_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{1}{\| x-y\|} } \, j(y)\, \mathrm{d}s_y}_{\displaystyle{ \mathrm{ SL}(j) }} - \underbrace{ \int\limits_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{\langle n(y) , x-y\rangle }{\| x-y\|^3} } \, m(y)\, \mathrm{d}s_y}_{\displaystyle{ \mathrm{DL}(m) }}\]

Define the geometry \(\Omega \subset \mathbb R^3\) and create a mesh:

sp = Sphere( (0,0,0), 1)
mesh = Mesh( OCCGeometry(sp).GenerateMesh(maxh=0.2)).Curve(4)

Create the finite element spaces for \(H^{-\frac12}(\Gamma)\) and \(H^{\frac12}(\Gamma)\) according to the given mesh:

fesL2 = SurfaceL2(mesh, order=3, dual_mapping=True)
u,v = fesL2.TnT()
fesH1 = H1(mesh, order=4)
uH1,vH1 = fesH1.TnT()
print ("ndofL2 = ", fesL2.ndof, "ndof H1 = ", fesH1.ndof)

ndofL2 =  7220 ndof H1 =  21241

Compute the interpolation of exact Dirichlet data \(u_0\) in finite element space:

uexa = 1/ sqrt( (x-1)**2 + (y-1)**2 + (z-1)**2 )
m = GridFunction(fesH1)
m.Interpolate (uexa)
Draw (m, mesh, draw_vol=False, order=3)

BaseWebGuiScene

Boundary Integral Equation

We carefully apply the Dirichlet trace on the ansatz for \(u\), we get a boundary integral equation for the Neumann data \(j \in H^{-\frac12}(\Gamma)\) which can be read as follows: determine \(j\) such that

\[\begin{split}\begin{array}{r c ccccc} \forall \, v\in H^{-\frac12}(\Gamma): && \displaystyle \int\limits_\Gamma \gamma_0 \left(\mathrm{SL}(j)\right) \cdot v \, \mathrm d s &=& \displaystyle \int\limits_\Gamma m \cdot v \, \mathrm d s &+& \displaystyle \int\limits_\Gamma \gamma_0 \left(\mathrm{DL}(m)\right)\cdot v \, \mathrm d s \\[2ex] & \Rightarrow & \mathrm{V}\,\mathrm{j} &=& \displaystyle \frac12\,\mathrm{M} \, \mathrm m &+& \mathrm K\, \mathrm m \,. \end{array} \end{split}\]

The discretisation of the above variational formulation is the boundary element method and it leads to a system of linear equations, ie

\[ \mathrm{V} \, \mathrm{j} = \left( \frac12 \,\mathrm{M} + \mathrm{K} \right) \, \mathrm{m}\,, \]

where the linear operators are as follows

• \(\mathrm{V}\) is the single layer operator. \(\mathrm V\) is regular and symmetric.

• \(\mathrm{M}\) is a mass matrix.

• \(\mathrm{K}\) is the double layer operator.

Details

For any \(v_i, u_j\in H^{-\frac12}(\Gamma)\) and \(w_j \in H^{\frac12}(\Gamma)\) the layer potential operators are implemented as follows

\[\begin{split} \begin{array}{rcl cl} V_{ij} &=& \langle v_i, \gamma_0 \mathrm{SL}(u_j) \rangle_{-\frac12} &=& \displaystyle \frac{1}{4\pi}\, \displaystyle \int\limits_\Gamma \int\limits_\Gamma \frac{1}{\|x-y\|} \, u_j(x) \, v_i(y) \, \mathrm{d} s_{y} \, \mathrm{d} s_x\,, \\[2ex] K_{ij} &=& \langle v_i, \gamma_0 \mathrm{DL}(w_j) \rangle_{-\frac12} &=& \displaystyle \frac{1}{4\pi} \int\limits_\Gamma\int\limits_\Gamma n_y\cdot \nabla_y \frac{1}{\|x-y\|} \, w_j(x) \, v_i(y) \, \mathrm{d} s_{y} \, \mathrm{d} s_x \end{array}\end{split}\]

You may assemble \(V\) and \(K\) it by

V = LaplaceSL(u*ds)*v*ds
K = LaplaceDL(w*ds)*v*ds

We generate the linear operators and solve for the Neumann data \(j\) using an iterative solver:

j = GridFunction(fesL2)
pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()
with TaskManager():
    V = LaplaceSL(u*ds)*v*ds
    K = LaplaceDL(uH1*ds)*v*ds
    M = BilinearForm(uH1*v*ds(bonus_intorder=3)).Assemble()
    
    rhs = ( (0.5 * M.mat + K.mat)*m.vec).Evaluate()
    CG(mat = V.mat, pre=pre, rhs = rhs, sol=j.vec, tol=1e-8, maxsteps=50, initialize=False, printrates=True)

Draw (j, mesh, draw_vol=False, order=3);

CG iteration 1, residual = 0.268824422385237     

CG iteration 2, residual = 0.06553516843504772     

CG iteration 3, residual = 0.015537002782918578     

CG iteration 4, residual = 0.0036775737524357664     

CG iteration 5, residual = 0.0009541450800419587     

CG iteration 6, residual = 0.0013667099821600933     

CG iteration 7, residual = 0.0001965670972102257     

CG iteration 8, residual = 5.780268069992089e-05     

CG iteration 9, residual = 4.635667988240774e-05     

CG iteration 10, residual = 3.400005858136626e-05     

CG iteration 11, residual = 1.7052395938330515e-05     

CG iteration 12, residual = 9.985434752221569e-06     

CG iteration 13, residual = 4.681261131406988e-06     

CG iteration 14, residual = 7.113358041663103e-06     

CG iteration 15, residual = 2.3904715149327493e-06     

CG iteration 16, residual = 2.0860619160176337e-06     

CG iteration 17, residual = 1.9082985895408567e-06     

CG iteration 18, residual = 1.779846278067858e-06     

CG iteration 19, residual = 9.154756115316215e-07     

CG iteration 20, residual = 6.251931615232174e-07     

CG iteration 21, residual = 5.05398203105611e-07     

CG iteration 22, residual = 5.253261069708412e-07     

CG iteration 23, residual = 4.981958341841999e-07     

CG iteration 24, residual = 3.8132848149698614e-07     

CG iteration 25, residual = 2.0579513254780794e-07     

CG iteration 26, residual = 1.0082995179585555e-07     

CG iteration 27, residual = 8.453288957781652e-08     

CG iteration 28, residual = 5.6118560372402515e-08     

CG iteration 29, residual = 9.925617533739608e-08     

CG iteration 30, residual = 8.751468220985332e-08     

CG iteration 31, residual = 4.602690850806891e-08     

CG iteration 32, residual = 2.7339424790280165e-08     

CG iteration 33, residual = 2.01747790126891e-08     

CG iteration 34, residual = 2.555118265967903e-08     

CG iteration 35, residual = 1.3194208556726591e-08     

CG iteration 36, residual = 7.544684405739223e-09     

CG iteration 37, residual = 2.2661931012568926e-08     

CG iteration 38, residual = 5.114872305783289e-09     

CG iteration 39, residual = 6.47815031131876e-09     

CG iteration 40, residual = 4.45093472117372e-09     

CG iteration 41, residual = 2.4100121389551915e-09     

Let’s have a look at the exact Neumann data and compute the error of the numerical solution:

graduexa = CF( (uexa.Diff(x), uexa.Diff(y), uexa.Diff(z)) )
n = specialcf.normal(3)
jexa = graduexa*n
Draw (jexa, mesh, draw_vol=False, order=3);
print ("L2-error =", sqrt (Integrate ( (jexa-j)**2, mesh, BND)))

L2-error = 0.00025059471684472587

Project repformula in H1

Evaluate the solution on a screen inside \(\Omega\)

screen = WorkPlane(Axes( (0,0,0), Z, X)).RectangleC(0.5,0.5).Face()
screen.faces.name="screen"
vismesh = screen.GenerateMesh(maxh=0.05)

LSPotential = LaplaceSL( u*ds )
LDPotential = LaplaceDL( uH1*ds )

repformula = LSPotential(j) - LDPotential(m)
fes_screen = H1(vismesh, order=5)
gf_screen = GridFunction(fes_screen)
with TaskManager():
    gf_screen.Set(repformula, definedon=vismesh.Boundaries("screen"))
Draw(gf_screen)
Draw(repformula - uexa, vismesh)

BaseWebGuiScene

type(LSPotential)

ngsolve.bem.PotentialOperator