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: homogeneous Neumann bvp, hypersingular operator \(D\)

Neumann Laplace Indirect Method

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&=& j, \quad &\Gamma = \partial \Omega\,.\end{array} \right. \end{split}\]

The solution \(u\in H^1(\Omega)\) is given by

\[ (1) \quad \quad \quad 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}s_y}_{\displaystyle{ \mathrm{LaplaceDL}(m) }}\,,\]

where \(\mathrm m\) solves the boundary integral equation

\[ \left( \mathrm D + \mathrm S\right) \,\mathrm m = - \mathrm M \, \mathrm j\quad \text{on}\; \Gamma\,.\]

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

sp = Sphere( (0,0,0), 1)
mesh = Mesh( OCCGeometry(sp).GenerateMesh(maxh=0.3)).Curve(2)

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

ndof H1 =  630

Consider as Neumann data \(j = n\cdot \nabla \displaystyle \frac{1}{\|x-x_0\|}\), thus

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

n = specialcf.normal(3)
j = gradu_exa*n

Boundary Integral Equation

Let’s carefully apply the Neumann trace to (1) and derive a boundary integral equation for \(m\)

\[ \forall \, v\in H^{\frac12}(\Gamma): \quad \displaystyle \int\limits_\Gamma v \cdot \gamma_1 \left(\mathrm{DL}(m) \right) \, \mathrm d s = -\displaystyle \int\limits_\Gamma j \cdot v \,\mathrm d s \]

The operator on the right is called the hypersingular operator. To remove the 1d kernel of \(\mathrm D\) one needs a stabilization \(\mathrm S\), for instance

\[ S \in \mathbb R^{n\times n}, \quad S_{ij} = \int\limits_{\Gamma} v_i(x) \,\mathrm{d} s \cdot \int\limits_{\Gamma} u_j(x) \,\mathrm{d} s\]

Given this stabilization, the operator \((\mathrm D + \mathrm S)\) is regular and symmetric and the resulting linear system is uniquely solvable:

\[ \left( \mathrm D + \mathrm S\right) \,\mathrm m = - \mathrm M \,\mathrm j \]

Note that for a trial function \(u_j\) and a test function \(v_i\), the \(ij\)-th entry of \(\mathrm D\) reads

\[ D_{ij} = \frac{1}{4\pi} \int\limits_\Gamma\int\limits_\Gamma \frac{ \langle \mathrm{\boldsymbol {curl}}_\Gamma \,u_j(x), \mathrm{\boldsymbol{curl}}_\Gamma\, v_i(y) \rangle}{\|x-y\|} \, \mathrm{d} s_{y} \, \mathrm{d} s_x \]

Thus, we obtain the hypersingular operator \(\mathrm D\) in NGSBEM based on \(\mathrm{LaplaceSL}\) with modified test and trial functions:

D = LaplaceSL(curl(u)*ds)*curl(v)*ds

We compute \(\mathrm D + \mathrm S\), the right hand side and solve the resulting system of linear equations for \(\mathrm m\):

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

rhs = LinearForm(-j*vH1.Trace()*ds(bonus_intorder=3)).Assemble()

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

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

Evaluation of the solution

Use the represenation formula \((1)\) and evaluate the solution \(u\) wherever you want, for instance on a screen:

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

DLPotential = LaplaceDL( uH1*ds)
repformula = DLPotential(m)
fes_screen = H1(vismesh, order=5)
gf_screen = GridFunction(fes_screen)
with TaskManager():
    gf_screen.Set(repformula, definedon=vismesh.Boundaries("screen"))
Draw(repformula,vismesh)

BaseWebGuiScene

Draw(repformula - u_exa, vismesh) # this must be a constant 

BaseWebGuiScene

References

For details on the explicit representation of the hypersingular operator have a look into Numerische Näherungsverfahren für elliptische Randwertprobleme, p.127, p.259 (1st edition).