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{array}{r}\{\begin{array}{rlcl}\mathrm{\Delta }u& =& 0,& \mathrm{\Omega }\subset {\mathbb{R}}^3,\\ {\gamma }_{1}u& =& u_1,& \mathrm{\Gamma }=\partial \mathrm{\Omega }.\end{array}\end{array}$$

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

$$u(x)=\underset{{\displaystyle \mathrm{DL}(m)}}{\underbrace{\underset{\mathrm{\Gamma }}{\int }{\displaystyle \frac{1}{4\pi }\frac{⟨n(y),x-y⟩}{\Vert x-y{\Vert }^3}m(y)\mathrm{d}{\sigma }_y}}}.$$

Let’s carefully apply the Neumann trace ${\gamma }_1$ to the representation formula and solve the resulting boundary integral equation for $m\in H^{\frac{1}{2}}(\mathrm{\Gamma })$ by discretisation of the following variational formulation:

$$\mathrm{\forall }v\in H^{\frac{1}{2}}(\mathrm{\Gamma }):{⟨v,{\gamma }_1(\mathrm{DL}(m))⟩}_{-\frac{1}{2}}={⟨u_1,v⟩}_{-\frac{1}{2}}.$$

Define the geometry $\mathrm{\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^{\frac{1}{2}}(\mathrm{\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

$$(\mathrm{D}+\mathrm{S})\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=HypersingularOperator(fesH1, intorder=12)
    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.0926249988826615
CG iteration 2, residual = 0.5729023155859267
CG iteration 3, residual = 0.34312522978550714
CG iteration 4, residual = 0.23038616941319026
CG iteration 5, residual = 0.17839889103315
CG iteration 6, residual = 0.11050474751791985
CG iteration 7, residual = 0.07041146814586027
CG iteration 8, residual = 0.04245777835218712
CG iteration 9, residual = 0.02694897679878841
CG iteration 10, residual = 0.015528300861547082
CG iteration 11, residual = 0.009047035564036587
CG iteration 12, residual = 0.004930817009662449
CG iteration 13, residual = 0.0026692120928693457
CG iteration 14, residual = 0.00132895570675758
CG iteration 15, residual = 0.0006492170614669755
CG iteration 16, residual = 0.00030590671473326313
CG iteration 17, residual = 0.00014049422573658983
CG iteration 18, residual = 6.30555600001843e-05
CG iteration 19, residual = 2.7095898951731424e-05
CG iteration 20, residual = 1.1675735602123292e-05
CG iteration 21, residual = 4.6096567946550486e-06
CG iteration 22, residual = 1.7636141668020139e-06
CG iteration 23, residual = 6.502535922153314e-07
CG iteration 24, residual = 2.2994849726730342e-07
CG iteration 25, residual = 8.588956399610995e-08
CG iteration 26, residual = 4.088879988656464e-08
CG iteration 27, residual = 3.604982487471187e-08
CG iteration 28, residual = 3.898176981005858e-08
CG iteration 29, residual = 2.0294282982351927e-08
CG iteration 30, residual = 8.991301715217601e-09

Note: the hypersingular operator is implemented as follows

$$D=⟨v,{\gamma }_1\mathrm{DL}(m){⟩}_{-\frac{1}{2}}=\frac{1}{4\pi }\underset{\mathrm{\Gamma }}{\int }\underset{\mathrm{\Gamma }}{\int }\frac{⟨{\mathbf{curl}}_{\mathrm{\Gamma }}m(x),{\mathbf{curl}}_{\mathrm{\Gamma }}v(y)⟩}{\Vert x-y\Vert }\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).