This page was generated from unit-11.2-bem-Laplace/Laplace_NtD_indirect.ipynb.

11.2.3 Neumann Laplace Indirect Method

keys: homogeneous Neumann bvp, hypersingular operator

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

Consider the Neumann boundary value problem

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

The solution \(u\in H^1(\Omega)\) of the above bvp can be written in the following form (representation forumula)

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

Let’s carefully apply the Neumann trace \(\gamma_1\) to the representation formula and solve the resulting boundary integral equation for \(m \in H^{\frac12}(\Gamma)\) by discretisation of the following variational formulation:

\[\forall \, v\in H^{\frac12}(\Gamma): \left\langle v, \gamma_1 \left(\mathrm{DL}(m)\right) \right\rangle_{-\frac12} = \left\langle u_1, v\right\rangle_{-\frac12} \,.\]

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.3)).Curve(2)

Define the finite element space for \(H^{\frac12}(\Gamma)\) for given geometry :

fesH1 = H1(mesh, order=2, definedon=mesh.Boundaries(".*"))
uH1,vH1 = fesH1.TnT()
print ("ndof H1 = ", fesH1.ndof)

ndof H1 =  630

Define the right hand side with given Neumann data \(u_1\):

uexa = 1/ sqrt( (x-1)**2 + (y-1)**2 + (z-1)**2 )
graduexa = CF( (uexa.Diff(x), uexa.Diff(y), uexa.Diff(z)) )

n = specialcf.normal(3)
u1exa = graduexa*n
rhs = LinearForm(u1exa*vH1.Trace()*ds(bonus_intorder=3)).Assemble()

The discretisation of the variational formulation leads to a system of linear equations, ie

\[\left( \mathrm{D} + \mathrm{S}\right) \, \mathrm{m}= \mathrm{rhs}\,,\]

where \(\mathrm{D}\) is the hypersingular operator and the stabilisation \((\mathrm D + \mathrm{S})\) is regular and symmetric.

vH1m1 = LinearForm(vH1*1*ds(bonus_intorder=3)).Assemble()
S = (BaseMatrix(Matrix(vH1m1.vec.Reshape(1))))@(BaseMatrix(Matrix(vH1m1.vec.Reshape(fesH1.ndof))))

m = GridFunction(fesH1)
pre = BilinearForm(uH1*vH1*ds).Assemble().mat.Inverse(freedofs=fesH1.FreeDofs())
with TaskManager():
    D=LaplaceSL(curl(uH1)*ds)*curl(vH1)*ds
    CG(mat = D.mat+S, pre=pre, rhs = rhs.vec, sol=m.vec, tol=1e-8, maxsteps=200, initialize=False, printrates=True)

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

CG iteration 1, residual = 1.092624998882661
CG iteration 2, residual = 0.5729010293389226
CG iteration 3, residual = 0.3431253018753176
CG iteration 4, residual = 0.23039240259684254
CG iteration 5, residual = 0.17840078129319534
CG iteration 6, residual = 0.11050148131747657
CG iteration 7, residual = 0.07041182879523007
CG iteration 8, residual = 0.042459174063529524
CG iteration 9, residual = 0.026947807306561482
CG iteration 10, residual = 0.015528094493476547
CG iteration 11, residual = 0.00904690105963536
CG iteration 12, residual = 0.00493067157055675
CG iteration 13, residual = 0.002669245037366956
CG iteration 14, residual = 0.0013289096044564173
CG iteration 15, residual = 0.0006491869037722831
CG iteration 16, residual = 0.0003059003574703774
CG iteration 17, residual = 0.0001404871255974491
CG iteration 18, residual = 6.305215386724583e-05
CG iteration 19, residual = 2.7094055093138668e-05
CG iteration 20, residual = 1.1674905025452704e-05
CG iteration 21, residual = 4.609501618646041e-06
CG iteration 22, residual = 1.7634976055472485e-06
CG iteration 23, residual = 6.502119936568242e-07
CG iteration 24, residual = 2.2998515713272648e-07
CG iteration 25, residual = 8.601673743623391e-08
CG iteration 26, residual = 4.1166656003751387e-08
CG iteration 27, residual = 3.654195281848217e-08
CG iteration 28, residual = 3.9446866076846766e-08
CG iteration 29, residual = 2.0417995295223827e-08
CG iteration 30, residual = 9.03407969080151e-09

Note: the hypersingular operator is implemented as follows

\[D = \langle v, \gamma_1 \mathrm{DL}(m) \rangle_{-\frac12} = \frac{1}{4\pi} \int\limits_\Gamma\int\limits_\Gamma \frac{ \langle \mathrm{\boldsymbol {curl}}_\Gamma \,m(x), \mathrm{\boldsymbol{curl}}_\Gamma\, v(y) \rangle}{\|x-y\|} \, \mathrm{d} \sigma_{y} \, \mathrm{d} \sigma_x\]

Details for instance in Numerische Näherungsverfahren für elliptische Randwertprobleme, p.127, p.259 (1st edition).