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)\) 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}s_y}_{\displaystyle{ \mathrm{DL}(m) }}\,.\]

NGSolve Solution

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

Boundary Integral Equation

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:

\[\begin{split}\begin{array}{r c ccc} \forall \, v\in H^{\frac12}(\Gamma): && \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 \\[2ex] & \Rightarrow & \mathrm{D}\,\mathrm{m} &=& -\mathrm{M} \, \mathrm j \,. \end{array} \end{split}\]

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

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

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

Details

For any \(u_j, v_i \in H^{\frac12}(\Gamma)\), hypersingular operator is implemented as follows (reference below)

\[ D_{ij} = \big\langle v_i, \gamma_1 \mathrm{DL}(u_j) \big\rangle = \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 \]

This explicit represenation of \(D\) is the singular layer potential operator with vector valued surface curl operator applied on trial and test functions. You may assemble \(D\) it by

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

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
rhs = LinearForm(-j*vH1.Trace()*ds(bonus_intorder=3)).Assemble()

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

The assembly of the hypersingular operator \(D\) is declarative and close to their actual implementations (explicit formula are given below).

D = LaplaceSL(curl(uH1)*ds) * curl(vH1)*ds

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=True)

CG iteration 1, residual = 1.0926249988826613     
CG iteration 2, residual = 0.5729010293377133     

CG iteration 3, residual = 0.34312530188543783     
CG iteration 4, residual = 0.23039240265445773     

CG iteration 5, residual = 0.17840078149360278     
CG iteration 6, residual = 0.11050148175913184     

CG iteration 7, residual = 0.0704118295186313     
CG iteration 8, residual = 0.042459175466083104     

CG iteration 9, residual = 0.026947808618547873     
CG iteration 10, residual = 0.01552809516341134     

CG iteration 11, residual = 0.009046901180573764     
CG iteration 12, residual = 0.004930671448195088     

CG iteration 13, residual = 0.0026692449307565798     
CG iteration 14, residual = 0.001328909553687845     

CG iteration 15, residual = 0.0006491868795976759     

CG iteration 16, residual = 0.0003059003457644725     
CG iteration 17, residual = 0.00014048712019643013     

CG iteration 18, residual = 6.305215120164812e-05     
CG iteration 19, residual = 2.7094053530708934e-05     

CG iteration 20, residual = 1.1674904095289665e-05     
CG iteration 21, residual = 4.609501314776319e-06     

CG iteration 22, residual = 1.7634974779953164e-06     
CG iteration 23, residual = 6.502119613097517e-07     

CG iteration 24, residual = 2.2998514250342298e-07     
CG iteration 25, residual = 8.601672886033162e-08     

CG iteration 26, residual = 4.116665990094554e-08     
CG iteration 27, residual = 3.654200605955962e-08     

CG iteration 28, residual = 3.9447025996258865e-08     
CG iteration 29, residual = 2.0418076149332225e-08     

CG iteration 30, residual = 9.034054982639904e-09     

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

Evaluation of the solution

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)
Draw(repformula - u_exa, vismesh) # this is 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).