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

11.2.4. Neumann Laplace Direct Method¶

keys: homogeneous Neumann bvp, hypersingular operator, unknown Dirichlet data

[1]:

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 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\|} } \, u_1(y)\, \mathrm{d}\sigma_y}_{\displaystyle{ \mathrm{SL}(u_1) }} + \underbrace{ \int\limits_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{\langle n(y), x-y \rangle}{\| x-y\|^3} } \, u_1( y)\, \mathrm{d}\sigma_y}_{\displaystyle{ \mathrm{DL}(u_1) }}\,.\]

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:

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

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

[2]:

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

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

[3]:

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 =  1446

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

[4]:

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

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

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

\[\left( \mathrm{D} + \mathrm{S}\right) \, \mathrm{u}_0 = \left( \frac12 \,\mathrm{M} - \mathrm{K}' \right) \, \mathrm{u}_1\,,\]

where the linear operators are as follows

\(\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.

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} \sigma \cdot \int\limits_{\Gamma} v_j(x) \,\mathrm{d} \sigma\,,\]

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

[5]:

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

Now we assemble the stabilised system matrix and the right hand side and solve for the Dirichlet data \(u_0\) with an iterative solver:

[6]:

u0 = GridFunction(fesH1)
pre = BilinearForm(uH1*vH1*ds).Assemble().mat.Inverse(freedofs=fesH1.FreeDofs())
with TaskManager():
    D = HypersingularOperator(fesH1, intorder=12)

    M = BilinearForm( v.Trace() * uH1 * ds(bonus_intorder=3)).Assemble()
    Kt = DoubleLayerPotentialOperator(fesL2, fesH1, intorder=12)
    rhs = ( (0.5 * M.mat.T - Kt.mat) * u1.vec ).Evaluate()

    CG(mat=D.mat+S, pre=pre, rhs=rhs, sol=u0.vec, tol=1e-8, maxsteps=200, printrates=True)

CG iteration 1, residual = 0.6734144317722747
CG iteration 2, residual = 0.3570539379523582
CG iteration 3, residual = 0.20848860722515933
CG iteration 4, residual = 0.16312663219143442
CG iteration 5, residual = 0.128573795755395
CG iteration 6, residual = 0.09481253218480892
CG iteration 7, residual = 0.06093958329330518
CG iteration 8, residual = 0.04820646021252343
CG iteration 9, residual = 0.02978674071219027
CG iteration 10, residual = 0.0212113018573045
CG iteration 11, residual = 0.013534492769773133
CG iteration 12, residual = 0.008289487261873922
CG iteration 13, residual = 0.005443966188702585
CG iteration 14, residual = 0.0032119264262769906
CG iteration 15, residual = 0.0020385485830059793
CG iteration 16, residual = 0.0011894769220963002
CG iteration 17, residual = 0.0006848870433749331
CG iteration 18, residual = 0.0004044330140602837
CG iteration 19, residual = 0.00022073436851857695
CG iteration 20, residual = 0.00012681943819321353
CG iteration 21, residual = 6.907297599162601e-05
CG iteration 22, residual = 3.736019764514392e-05
CG iteration 23, residual = 1.9761100557887124e-05
CG iteration 24, residual = 9.964135469280613e-06
CG iteration 25, residual = 4.933925387272878e-06
CG iteration 26, residual = 2.478719966857495e-06
CG iteration 27, residual = 1.190054871278054e-06
CG iteration 28, residual = 5.710376983671422e-07
CG iteration 29, residual = 2.6226614954729127e-07
CG iteration 30, residual = 1.20893486112704e-07
CG iteration 31, residual = 5.4874600490775255e-08
CG iteration 32, residual = 2.7248227480447033e-08
CG iteration 33, residual = 1.8187818269437897e-08
CG iteration 34, residual = 1.8157784457816937e-08
CG iteration 35, residual = 2.095186858848211e-08
CG iteration 36, residual = 1.5817717347673397e-08
CG iteration 37, residual = 7.829425924931581e-09
CG iteration 38, residual = 3.4516503033609765e-09

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

[7]:

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

[8]:

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

L2 error of surface gradients = 0.00444984782070144

[ ]:

[ ]: