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)

WebGLScene

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.

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.1761858690412437     
CG iteration 2, residual = 0.5482159540767939     
CG iteration 3, residual = 0.13178153060209566     
CG iteration 4, residual = 0.032853845892815245     
CG iteration 5, residual = 0.010919299980605078     
CG iteration 6, residual = 0.008294644409316855     
CG iteration 7, residual = 0.0038205039348878714     
CG iteration 8, residual = 0.002005738060233477     
CG iteration 9, residual = 0.0014503915027362026     
CG iteration 10, residual = 0.0011780601751354618     
CG iteration 11, residual = 0.0012971112821272856     
CG iteration 12, residual = 0.0007335456334032436     
CG iteration 13, residual = 0.0004604698434819254     
CG iteration 14, residual = 0.00034107937467262923     
CG iteration 15, residual = 0.0002950251360273075     
CG iteration 16, residual = 0.00028268445350020175     
CG iteration 17, residual = 0.00018488500826333723     
CG iteration 18, residual = 0.00012925883074687333     
CG iteration 19, residual = 8.502022723963184e-05     
CG iteration 20, residual = 6.54213818473314e-05     
CG iteration 21, residual = 5.9689393990408824e-05     
CG iteration 22, residual = 3.955464167448273e-05     
CG iteration 23, residual = 3.107788319811991e-05     
CG iteration 24, residual = 2.3558518902509442e-05     
CG iteration 25, residual = 1.5739201790557147e-05     
CG iteration 26, residual = 1.5138826377322623e-05     
CG iteration 27, residual = 9.634735016944813e-06     
CG iteration 28, residual = 1.0103579032642797e-05     
CG iteration 29, residual = 7.923944359993876e-06     
CG iteration 30, residual = 5.711770934231516e-06     
CG iteration 31, residual = 4.446253370139276e-06     
CG iteration 32, residual = 3.997513969568856e-06     
CG iteration 33, residual = 4.5825248981651105e-06     
CG iteration 34, residual = 2.6948330100374843e-06     
CG iteration 35, residual = 1.7794915589421282e-06     
CG iteration 36, residual = 1.4933102547839452e-06     
CG iteration 37, residual = 1.6699350307758e-06     

CG iteration 38, residual = 9.792247045429012e-07     
CG iteration 39, residual = 7.600226251795564e-07     
CG iteration 40, residual = 5.47922615772629e-07     
CG iteration 41, residual = 7.072751837760491e-07     
CG iteration 42, residual = 7.71089086308538e-07     
CG iteration 43, residual = 4.4900292851635705e-07     
CG iteration 44, residual = 2.610829087319441e-07     
CG iteration 45, residual = 2.0549474994708963e-07     
CG iteration 46, residual = 1.615166417989887e-07     
CG iteration 47, residual = 1.28786738521156e-07     
CG iteration 48, residual = 9.049652319796454e-08     
CG iteration 49, residual = 6.438935290375447e-08     
CG iteration 50, residual = 6.531489321167287e-08     
CG iteration 51, residual = 4.8142521703750766e-08     
CG iteration 52, residual = 3.9785333648790946e-08     
CG iteration 53, residual = 3.066050317790361e-08     
CG iteration 54, residual = 7.497694214380977e-08     
CG iteration 55, residual = 2.1701559500915442e-08     

Draw (j, order=3);