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

11.2.2. Dirichlet Laplace Direct Method¶

keys: homogeneous Dirichlet bvp, double and single layer potential, unknown Neumann data

[1]:

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

Consider the Dirichlet boundary value problem

\[\begin{split}\left\{ \begin{array}{rcl l} \Delta u &=& 0, \quad &\Omega \subset \mathbb R^3\,,\\ \gamma_0 u&=& u_0, \quad &\Gamma = \partial \Omega\,.\end{array} \right.\end{split}\]

Let us choose the following ansatz for the solution \(u\in H^1(\Omega)\) (direct ansatz)

\[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_0(y)\, \mathrm{d}\sigma_y}_{\displaystyle{ \mathrm{DL}(u_0) }}\]

and solve for the Neumann data \(u_1 \in H^{-\frac12}(\Gamma)\) by the boundary element method, i.e.,

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

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

[2]:

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

Create the finite element spaces for \(H^{-\frac12}(\Gamma)\) and \(H^{\frac12}(\Gamma)\) according to the given mesh:

[3]:

fesL2 = SurfaceL2(mesh, order=3, dual_mapping=True)
u,v = fesL2.TnT()
fesH1 = H1(mesh, order=4)
uH1,vH1 = fesH1.TnT()
print ("ndofL2 = ", fesL2.ndof, "ndof H1 = ", fesH1.ndof)

ndofL2 =  7220 ndof H1 =  21241

Compute the interpolation of exact Dirichlet data \(u_0\) in finite element space:

[4]:

uexa = 1/ sqrt( (x-1)**2 + (y-1)**2 + (z-1)**2 )
u0 = GridFunction(fesH1)
u0.Interpolate (uexa)
Draw (u0, mesh, draw_vol=False, order=3);

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

\[\mathrm{V} \, \mathrm{u}_1 = \left( \frac12 \,\mathrm{M} + \mathrm{K} \right) \, \mathrm{u}_0\,,\]

where the linear operators are as follows

\(\mathrm{V}\) is the single layer operator. \(\mathrm V\) is regular and symmetric.
\(\mathrm{M}\) is a mass matrix.
\(\mathrm{K}\) is the double layer operator.

We approximate the linear operators as hmatrices and solve for the Neumann data \(u_1\) with an iterative solver:

[5]:

u1 = GridFunction(fesL2)
pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()

with TaskManager():
    V = SingleLayerPotentialOperator(fesL2, intorder=12)

    M = BilinearForm( uH1 * v.Trace() * ds(bonus_intorder=3)).Assemble()
    K = DoubleLayerPotentialOperator(fesH1, fesL2, intorder=12)

    rhs = ( (0.5 * M.mat + K.mat)*u0.vec).Evaluate()
    CG(mat = V.mat, pre=pre, rhs = rhs, sol=u1.vec, tol=1e-8, maxsteps=50, initialize=False, printrates=True)

CG iteration 1, residual = 0.26882443508338744
CG iteration 2, residual = 0.06553517738164745
CG iteration 3, residual = 0.015537006646393927
CG iteration 4, residual = 0.0036775507945853254
CG iteration 5, residual = 0.0009473299290083554
CG iteration 6, residual = 0.001380324500133172
CG iteration 7, residual = 0.0001965930628898029
CG iteration 8, residual = 5.7817188553410047e-05
CG iteration 9, residual = 4.643518369001456e-05
CG iteration 10, residual = 3.399771904913404e-05
CG iteration 11, residual = 1.7058052504812094e-05
CG iteration 12, residual = 9.988909760938636e-06
CG iteration 13, residual = 4.669062224474859e-06
CG iteration 14, residual = 2.4878753502716016e-06
CG iteration 15, residual = 5.520624230728082e-06
CG iteration 16, residual = 2.092351591536956e-06
CG iteration 17, residual = 1.910708775491434e-06
CG iteration 18, residual = 1.787457210545918e-06
CG iteration 19, residual = 9.145440804791857e-07
CG iteration 20, residual = 6.252906187726185e-07
CG iteration 21, residual = 5.046629337801763e-07
CG iteration 22, residual = 3.8856764941159945e-07
CG iteration 23, residual = 3.7945833890639973e-07
CG iteration 24, residual = 8.134049483660121e-07
CG iteration 25, residual = 2.0911889490982828e-07
CG iteration 26, residual = 1.0088778212629894e-07
CG iteration 27, residual = 8.446741259154526e-08
CG iteration 28, residual = 5.03453221123945e-08
CG iteration 29, residual = 5.5207582766197555e-08
CG iteration 30, residual = 2.678767951488126e-08
CG iteration 31, residual = 5.1129608464688227e-08
CG iteration 32, residual = 3.460408363198751e-08
CG iteration 33, residual = 2.581839946916422e-08
CG iteration 34, residual = 2.5905747755004283e-08
CG iteration 35, residual = 1.3211024761808669e-08
CG iteration 36, residual = 7.520154524258728e-09
CG iteration 37, residual = 4.5004899216867595e-09
CG iteration 38, residual = 4.859490757214749e-09
CG iteration 39, residual = 3.0818619428720632e-09
CG iteration 40, residual = 3.5451068474293164e-09
CG iteration 41, residual = 5.557929992739343e-09
CG iteration 42, residual = 2.0672106541006666e-09

[6]:

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

Let’s have a look at the exact Neumann data and compute the error of the numerical solution:

[7]:

graduexa = CF( (uexa.Diff(x), uexa.Diff(y), uexa.Diff(z)) )
n = specialcf.normal(3)
u1exa = graduexa*n
Draw (u1exa, mesh, draw_vol=False, order=3);
print ("L2-error =", sqrt (Integrate ( (u1exa-u1)**2, mesh, BND)))

L2-error = 0.000249386949866592

[ ]:

[ ]: