Dirichlet Laplace Indirect Method using Quads

Dirichlet Laplace Indirect Method using Quads#

keys: homogeneous Dirichlet bvp, single layer potential, quads

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:

sp = Glue(Sphere( (0,0,0), 1).faces)
mesh = Mesh(OCCGeometry(sp).GenerateMesh(maxh=0.5, quad_dominated=True)).Curve(3)
Draw (mesh)

BaseWebGuiScene

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

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

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

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\):

j = GridFunction(fesL2)
pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()
SetNumThreads(1)
with TaskManager():
    V = LaplaceSL(u*ds(bonus_intorder=4), use_fmm=False)*v*ds(bonus_intorder=4)
    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.1761858690412383     
CG iteration 2, residual = 0.5482159540767929     
CG iteration 3, residual = 0.13178153060209558     
CG iteration 4, residual = 0.03285384589281548     
CG iteration 5, residual = 0.010919299980605566     
CG iteration 6, residual = 0.008294644409317674     
CG iteration 7, residual = 0.0038205039348886446     
CG iteration 8, residual = 0.002005738060234026     
CG iteration 9, residual = 0.0014503915027366291     
CG iteration 10, residual = 0.0011780601751356253     
CG iteration 11, residual = 0.0012971112821141527     
CG iteration 12, residual = 0.0007335456324858695     
CG iteration 13, residual = 0.0004604692295797129     
CG iteration 14, residual = 0.00034085779168566756     
CG iteration 15, residual = 0.00024045650242522147     
CG iteration 16, residual = 0.00038271190858650523     
CG iteration 17, residual = 0.0001852942812881629     
CG iteration 18, residual = 0.0001292607402280529     
CG iteration 19, residual = 8.502022862315073e-05     
CG iteration 20, residual = 6.542138184915657e-05     
CG iteration 21, residual = 5.968939399039654e-05     
CG iteration 22, residual = 3.9554641674473404e-05     
CG iteration 23, residual = 3.1077883198118536e-05     
CG iteration 24, residual = 2.355851890233216e-05     
CG iteration 25, residual = 1.5739201598472318e-05     
CG iteration 26, residual = 1.5138785279509823e-05     
CG iteration 27, residual = 9.627841820366277e-06     
CG iteration 28, residual = 6.706899614829621e-06     
CG iteration 29, residual = 1.6242765915279338e-05     
CG iteration 30, residual = 5.742786458938375e-06     
CG iteration 31, residual = 4.482398371010459e-06     
CG iteration 32, residual = 5.042535429520525e-06     
CG iteration 33, residual = 4.0059059891993395e-06     
CG iteration 34, residual = 2.608980291121406e-06     
CG iteration 35, residual = 1.7806042511497512e-06     
CG iteration 36, residual = 1.5817344065230353e-06     
CG iteration 37, residual = 1.5780094579936627e-06     
CG iteration 38, residual = 9.763341712760746e-07     
CG iteration 39, residual = 7.603860982448828e-07     
CG iteration 40, residual = 5.759425320405671e-07     
CG iteration 41, residual = 6.826640827190862e-07     
CG iteration 42, residual = 4.0321911440785017e-07     
CG iteration 43, residual = 8.535485980897427e-07     
CG iteration 44, residual = 2.676282380535318e-07     
CG iteration 45, residual = 2.0549931923519701e-07     
CG iteration 46, residual = 1.615166248385285e-07     
CG iteration 47, residual = 1.2878673778499882e-07     
CG iteration 48, residual = 9.049652319051307e-08     

CG iteration 49, residual = 6.438935290493145e-08     
CG iteration 50, residual = 6.531489323191935e-08     
CG iteration 51, residual = 4.814251981858665e-08     
CG iteration 52, residual = 3.97842161679993e-08     
CG iteration 53, residual = 3.007442320954777e-08     
CG iteration 54, residual = 2.1751259693572604e-08     

Draw (j, order=3);