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.

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.1761858690412437     
CG iteration 2, residual = 0.5482159540767941     
CG iteration 3, residual = 0.13178153060209583     
CG iteration 4, residual = 0.03285384589281531     
CG iteration 5, residual = 0.010919299980605096     
CG iteration 6, residual = 0.008294644409316855     
CG iteration 7, residual = 0.0038205039348878536     
CG iteration 8, residual = 0.002005738060233457     
CG iteration 9, residual = 0.0014503915027361967     
CG iteration 10, residual = 0.001178060175136705     
CG iteration 11, residual = 0.0012971112822167217     
CG iteration 12, residual = 0.0007335456395832783     
CG iteration 13, residual = 0.0004604739796476631     
CG iteration 14, residual = 0.0003425674772581703     
CG iteration 15, residual = 0.0004837529973170583     
CG iteration 16, residual = 0.0002249479719049113     
CG iteration 17, residual = 0.0001847536668739514     
CG iteration 18, residual = 0.000129258219882603     
CG iteration 19, residual = 8.50202267970241e-05     
CG iteration 20, residual = 6.542138184673892e-05     
CG iteration 21, residual = 5.9689393990406825e-05     
CG iteration 22, residual = 3.955464167448066e-05     
CG iteration 23, residual = 3.107788319816929e-05     
CG iteration 24, residual = 2.3558518914678287e-05     
CG iteration 25, residual = 1.5739209381761213e-05     
CG iteration 26, residual = 1.514039145843612e-05     
CG iteration 27, residual = 9.891936931625937e-06     
CG iteration 28, residual = 2.2675974122171828e-05     
CG iteration 29, residual = 6.413103427416558e-06     
CG iteration 30, residual = 5.71096223269957e-06     
CG iteration 31, residual = 4.463813204978003e-06     
CG iteration 32, residual = 4.601365244511545e-06     
CG iteration 33, residual = 4.253034293643534e-06     
CG iteration 34, residual = 2.6203526354129033e-06     
CG iteration 35, residual = 1.7807569381512877e-06     
CG iteration 36, residual = 1.5870023112260822e-06     
CG iteration 37, residual = 1.5733208601148372e-06     
CG iteration 38, residual = 9.762262508581441e-07     
CG iteration 39, residual = 7.604349441199408e-07     
CG iteration 40, residual = 6.005184670341102e-07     
CG iteration 41, residual = 9.652336345557704e-07     
CG iteration 42, residual = 7.147744181810041e-07     
CG iteration 43, residual = 4.0323567047631645e-07     

CG iteration 44, residual = 2.609060380095094e-07     
CG iteration 45, residual = 2.0549282548095206e-07     
CG iteration 46, residual = 1.6151661802657004e-07     
CG iteration 47, residual = 1.287867377572548e-07     
CG iteration 48, residual = 9.049652319010544e-08     
CG iteration 49, residual = 6.438935289736774e-08     
CG iteration 50, residual = 6.531489427179293e-08     
CG iteration 51, residual = 4.814289959998808e-08     
CG iteration 52, residual = 3.99998493386514e-08     
CG iteration 53, residual = 7.380766240725458e-08     
CG iteration 54, residual = 3.129785828108193e-08     
CG iteration 55, residual = 2.148017855908136e-08     

Draw (j, order=3);