Neumann Laplace Direct Method

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, unknown Dirichlet data

Neumann Laplace Direct 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 unique solution \(u\in H^1(\Omega)\) can be written in the following form (representation forumula)

\[ 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) }}\,.\]

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

Define conforming finite element spaces for \(H^{\frac12}(\Gamma)\) and \(H^{-\frac12}(\Gamma)\) respectively:

fesL2 = SurfaceL2(mesh, order=2, dual_mapping=False) 
u,v = fesL2.TnT()
fesH1 = H1(mesh, order=2, definedon=mesh.Boundaries(".*"))
uH1,vH1 = fesH1.TnT()
print ("ndofL2 = ", fesL2.ndof, "ndof H1 = ", fesH1.ndof)

ndofL2 =  4332 ndof H1 =  1580

Project the given Neumann data \(u_1\) in finite element space and have a look at the boundary data:

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)
j = GridFunction(fesL2)
j.Interpolate(graduexa*n, definedon=mesh.Boundaries(".*"))

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

Boundary Integral Equation

Let’s carefully apply the Neumann trace \(\gamma_1\) to the representation formula and solve the resulting boundary integral equation for the Dirichlet trace \(u_0\) of \(u\) 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 &-& \displaystyle \int\limits_\Gamma v \cdot \gamma_1 \left(\mathrm{SL}(j)\right) \, \mathrm d s \\[2ex] & \Rightarrow & \mathrm{D}\,\mathrm{m} &=& \dfrac12 \, \mathrm{M} \, \mathrm j &-& \mathrm{K}' \, \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} = \left( \frac12 \,\mathrm{M} - \mathrm{K}' \right) \, \mathrm{j} \,, \]

where the linear operators are as follows

\(\mathrm{S}\) is a stabilisation matrix that removes the one dimensional kernel from hypersinbular operator.
\(\mathrm{D}\) is the hypersingular operator and the stabilisation \((\mathrm D + \mathrm{S})\) is regular and symmetric.
\(\mathrm{M}\) is a mass matrix.
\(\mathrm{K}'\) is the adjoint double layer operator.

Details

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

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

Recall that for any \(v_i\in H^{-\frac12}(\Gamma)\) and \(w_j \in H^{\frac12}(\Gamma)\) the double layer potential operators \(\mathrm K\) has been given by

\[ \begin{array}{rcl cl} 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}\]

And this means that the adjoint double layer potential operator equals its transposed, \(\mathrm K^\intercal\).

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

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

The stabilisation matrix \(S\) is needed to cope with the kernel of the hypersingular operator. The following stabilisation matrix is used here:

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

where \(v_i, v_j\) being \(H^{\frac12}(\Gamma)\)-conforming shape functions.

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 is declarative and close to their actual implementations (explicit formula are given below).

with TaskManager(): 
    D = LaplaceSL(curl(uH1)*ds)*curl(vH1)*ds
    K = LaplaceDL(uH1*ds)*v*ds
    M = BilinearForm( uH1 * v * ds(bonus_intorder=3)).Assemble()

Left and right hand side of linear system are computed and the solution is the Dirichlet data \(u_0\):

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

CG iteration 1, residual = 0.6734143938795208

CG iteration 2, residual = 0.3570586989565186

CG iteration 3, residual = 0.20905841633641048

CG iteration 4, residual = 0.17782730725965978

CG iteration 5, residual = 0.1240304646603016

CG iteration 6, residual = 0.09705555483495114

CG iteration 7, residual = 0.059686364640765925

CG iteration 8, residual = 0.048358396284000675

CG iteration 9, residual = 0.028887223089879264

CG iteration 10, residual = 0.02087591545709992

CG iteration 11, residual = 0.013576293553324575

CG iteration 12, residual = 0.008369259002401461

CG iteration 13, residual = 0.005519552288092554

CG iteration 14, residual = 0.003273842526139349

CG iteration 15, residual = 0.0021009042394241127

CG iteration 16, residual = 0.001230772030359319

CG iteration 17, residual = 0.0007221780720367607

CG iteration 18, residual = 0.00042551775688713883

CG iteration 19, residual = 0.0002320103253771422

CG iteration 20, residual = 0.00013264368940199342

CG iteration 21, residual = 7.176651710682202e-05

CG iteration 22, residual = 3.9139160079288725e-05

CG iteration 23, residual = 2.0404247412155772e-05

CG iteration 24, residual = 1.042677650806519e-05

CG iteration 25, residual = 5.204655489842847e-06

CG iteration 26, residual = 2.5870475716945516e-06

CG iteration 27, residual = 1.2572803365940249e-06

CG iteration 28, residual = 6.09225631215286e-07

CG iteration 29, residual = 2.8568758167594216e-07

CG iteration 30, residual = 1.2785128960165276e-07

CG iteration 31, residual = 5.716856867904242e-08

CG iteration 32, residual = 2.827875076787529e-08

CG iteration 33, residual = 1.8645284671804463e-08

CG iteration 34, residual = 1.8353008222638013e-08

CG iteration 35, residual = 2.1312773512214396e-08

CG iteration 36, residual = 1.684434996512971e-08

CG iteration 37, residual = 8.623457687244244e-09

CG iteration 38, residual = 3.8171683132670186e-09

Let’s have a look at the numerical and exact Dirichlet data and compare them quantitatively:

Draw (m, mesh, draw_vol=False, order=3);
mexa = GridFunction(fesH1)
mexa.Interpolate (uexa, definedon=mesh.Boundaries(".*"))
Draw (mexa, mesh, draw_vol=False, order=3);

print ("L2 error of surface gradients =", sqrt (Integrate ( (grad(mexa)-grad(m))**2, mesh, BND)))

L2 error of surface gradients = 0.004568783275440448

Evaluation 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)

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) # exact up to 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).