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

11.2.1 Dirichlet Laplace Indirect MethodΒΆ

keys: homogeneous Dirichlet bvp, single layer potential, ACA

[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

We 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)\) (indirect ansatz)

\[u(x) = \underbrace{ \int\limits_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{1}{\| x-y\|} } \, j(y)\, \mathrm{d}\sigma_y }_{\displaystyle{ \mathrm{SL}(j) } }\]

and solve for the density \(j\in H^{-\frac12}(\Gamma)\) by the boundary element method, i.e. the numerical solution of the variational formulation

\[\forall \, v\in H^{-\frac12}(\Gamma): \quad \left\langle \gamma_0 \left(\mathrm{SL}(j)\right), v \right\rangle_{-\frac12} = \left\langle u_0, 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.3)).Curve(4)

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

[3]:

fesL2 = SurfaceL2(mesh, order=4, dual_mapping=True)
u,v = fesL2.TnT()

Define Dirichlet data \(u_0\) and compute the right hand side vector:

[4]:

u0 = 1/ sqrt( (x-1)**2 + (y-1)**2 + (z-1)**2 )
rhs = LinearForm (u0*v.Trace()*ds(bonus_intorder=3)).Assemble()

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

\[\mathrm{V} \, \mathrm{j} = \mathrm{rhs} \,,\]

where \(\mathrm{V}\) is the single layer potential operator. \(\mathrm V\) is regular and symmetric.

Demo 1: Assemble the single layer operator \(V\) as dense matrix and solve for unknwon density \(j\):

[5]:

j = GridFunction(fesL2)
pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()
with TaskManager():
    # V = SingleLayerPotentialOperator(fesL2, intorder=10)
    V = LaplaceSL(u*ds)*v*ds
    CG(mat = V.mat, pre=pre, rhs = rhs.vec, sol=j.vec, tol=1e-8, maxsteps=200, initialize=False, printrates=True)

CG iteration 1, residual = 2.185606954606291
CG iteration 2, residual = 0.5510990298605161
CG iteration 3, residual = 0.13112741187311436
CG iteration 4, residual = 0.03107892126080936
CG iteration 5, residual = 0.007372481975124272
CG iteration 6, residual = 0.0017929371359740678
CG iteration 7, residual = 0.0005259800974227861
CG iteration 8, residual = 0.00025807777881830575
CG iteration 9, residual = 0.00028323794974845196
CG iteration 10, residual = 0.00012036647881309768
CG iteration 11, residual = 0.00023507220246880676
CG iteration 12, residual = 4.4923691859943966e-05
CG iteration 13, residual = 4.801474719537957e-05
CG iteration 14, residual = 2.7825229462712193e-05
CG iteration 15, residual = 1.4627805195039993e-05
CG iteration 16, residual = 1.4781005107860767e-05
CG iteration 17, residual = 9.771555021457255e-06
CG iteration 18, residual = 6.175956781518063e-06
CG iteration 19, residual = 1.5387391127938983e-05
CG iteration 20, residual = 3.2801501501799858e-06
CG iteration 21, residual = 2.970067358732096e-06
CG iteration 22, residual = 3.493169500409302e-06
CG iteration 23, residual = 2.5937980739596076e-06
CG iteration 24, residual = 1.2604332726056218e-06
CG iteration 25, residual = 8.372899332359421e-07
CG iteration 26, residual = 1.2440384360425716e-06
CG iteration 27, residual = 5.872586630360098e-07
CG iteration 28, residual = 4.4146186963787075e-07
CG iteration 29, residual = 5.781842335656235e-07
CG iteration 30, residual = 3.526775204323348e-07
CG iteration 31, residual = 1.778225530958499e-07
CG iteration 32, residual = 1.5441947126512715e-07
CG iteration 33, residual = 2.1911326525546733e-07
CG iteration 34, residual = 1.2550257371311674e-07
CG iteration 35, residual = 8.723802303667268e-08
CG iteration 36, residual = 6.888878659916675e-08
CG iteration 37, residual = 1.0272145097732246e-07
CG iteration 38, residual = 3.4372701134325734e-08
CG iteration 39, residual = 4.0240015693504245e-08
CG iteration 40, residual = 5.7328898023203763e-08
CG iteration 41, residual = 1.8547250338573198e-08

[6]:

Draw (j, order=3);

Notes:

[ ]: