Dirichlet Laplace Indirect Method

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

keys: electrostatics, homogeneous Dirichlet bvp, single layer potential \(V\)

Dirichlet Laplace Indirect Method#

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&=& m, \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}s_y }_{\displaystyle{ \mathrm{SL}(j) } }\]

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

sp = Sphere( (0,0,0), 1)
mesh = Mesh( OCCGeometry(sp).GenerateMesh(maxh=0.3)).Curve(3)

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

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

Boundary Integral Equation

We carefully apply the Dirichlet trace on the ansatz for \(u\) and get an integral equation for the unknown density \(j\in H^{-\frac12}(\Gamma)\). The integral equation may be read as follows: determine \(j\) such that

\[\begin{split}\begin{array}{r c ccc} \forall \, v\in H^{-\frac12}(\Gamma): && \displaystyle \int\limits_\Gamma \gamma_0 \left(\mathrm{SL}(j)\right) \cdot v \, \mathrm d s &=& \displaystyle \int\limits_\Gamma m \cdot v \, \mathrm d s \\[2ex] & \Rightarrow & \mathrm{V}\,\mathrm{j} &=& \mathrm{M} \, \mathrm m \,. \end{array} \end{split}\]

Details

For any \(u_j, v_i \in H^{-\frac12}(\Gamma)\), the single layer potential operator is implemented as follows

\[ V_{ij} = \langle v_i, \gamma_0 \mathrm{SL}(u_j) \rangle_{-\frac12} = \frac{1}{4\pi} \int\limits_\Gamma\int\limits_\Gamma \frac{1}{\|x-y\|} \, u_j(x) \, v_i(y) \, \mathrm{d} s_{y} \, \mathrm{d} s_x \]

You may assemble \(V\) it by

V = LaplaceSL(u*ds)*v*ds

Compute the given Dirichlet data \(m\), the right hand side and the matrix \(\mathrm V\) (single layer potential operator).

m = 1/ sqrt( (x-1)**2 + (y-1)**2 + (z-1)**2 )
rhs = LinearForm (m*v.Trace()*ds(bonus_intorder=3)).Assemble()
with TaskManager(pajetrace=1000*1000*1000):
    V = LaplaceSL(u*ds)*v*ds

Solve the linear system:

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

CG iteration 1, residual = 2.185632470081333     
CG iteration 2, residual = 0.5511118070234002     

CG iteration 3, residual = 0.1311310259046071     
CG iteration 4, residual = 0.031078990634675403     

CG iteration 5, residual = 0.007370161304676199     
CG iteration 6, residual = 0.0017872719313664826     

CG iteration 7, residual = 0.0005164640897886519     
CG iteration 8, residual = 0.0002492780566936905     

CG iteration 9, residual = 0.00027505280100685716     
CG iteration 10, residual = 0.0001804296883779487     

CG iteration 11, residual = 0.00013054631855996055     
CG iteration 12, residual = 3.7227312416526054e-05     

CG iteration 13, residual = 3.312346253569015e-05     
CG iteration 14, residual = 1.610364811063763e-05     

CG iteration 15, residual = 9.915860637612635e-06     
CG iteration 16, residual = 1.0658417152027922e-05     

CG iteration 17, residual = 7.466789251299928e-06     
CG iteration 18, residual = 7.402407545278815e-06     

CG iteration 19, residual = 5.4726515089507445e-06     
CG iteration 20, residual = 2.3569221270401377e-06     

CG iteration 21, residual = 1.647571884000155e-06     
CG iteration 22, residual = 1.713550652759644e-06     

CG iteration 23, residual = 1.3120434766960407e-06     
CG iteration 24, residual = 6.358196555359483e-07     

CG iteration 25, residual = 1.1080271990215965e-06     
CG iteration 26, residual = 2.6228375487977723e-07     

CG iteration 27, residual = 2.3839343618324335e-07     
CG iteration 28, residual = 2.521414429331238e-07     

CG iteration 29, residual = 1.3798656192496498e-07     
CG iteration 30, residual = 1.3285753304935006e-07     

CG iteration 31, residual = 7.050800576568768e-08     
CG iteration 32, residual = 1.8106090955115502e-07     

CG iteration 33, residual = 5.2934875936904487e-08     
CG iteration 34, residual = 3.767213790499984e-08     

CG iteration 35, residual = 3.122043118455898e-08     
CG iteration 36, residual = 2.9167854465250435e-08     

CG iteration 37, residual = 1.3858160148081166e-08     

type(V)
#for t in Timers():
#    if "ngbem" in t["name"]:
#        print (t)

ngsolve.bem.IntegralOperator

screen = WorkPlane(Axes( (0,0,0), Z, X)).RectangleC(0.5,0.5).Face()
screen.faces.name="screen"
vismesh = screen.GenerateMesh(maxh=0.05)
#Draw (vismesh);

SLPotential = LaplaceSL( u*ds)
repformula = SLPotential(j)
fes_screen = H1(vismesh, order=5)
gf_screen = GridFunction(fes_screen)
with TaskManager():
    gf_screen.Set(repformula, definedon=vismesh.Boundaries("screen"))
Draw(repformula,vismesh)

BaseWebGuiScene

References