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}\]

The solution \(u\in H^1(\Omega)\) is given by

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

where \(\mathrm j \) solves the boundary integral equation

\[ \mathrm V\,\mathrm j = \mathrm M\,\mathrm m\quad \text{on} \; \Gamma\,.\]

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

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

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

Consider as Dirichlet data \(m = \displaystyle \frac{1}{\|x-x_0\|}\), thus

m = 1/ sqrt( (x-1)**2 + (y-1)**2 + (z-1)**2 )

Boundary Integral Equation

We carefully apply the Dirichlet trace to (1) and derive a boundary integral equation for the unknown density \(j\):

\[\forall \, v\in H^{-\frac12}(\Gamma): \quad \displaystyle \int_\Gamma \gamma_0 \left(\mathrm{LaplaceSL}(j)\right) \cdot v \, \mathrm d s = \displaystyle \int_\Gamma m \cdot v \, \mathrm d s \]

The operator on the left is known as single layer potential operator - in NGSBEM, you assemble it by

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

Compute the right hand side \(\mathrm M\,\mathrm m\) and \(\mathrm V\)

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=False)
Draw (j);

Evaluation Representation Formula

Use the represenation formula \((1)\) and evaluate the solution \(u\) wherever you want, for instance on a screen:

type(V)

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)

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

Details: The Single Layer Potential Operator

Given a trial function \(u_j\) and a test function \(v_i\), the \(ij\)-th entry of the single layer potential operator reads

\[ V_{ij} = \frac{1}{4\pi} \int_\Gamma \int_\Gamma \frac{1}{\|x-y\|} \, u_j(x) \, v_i(y) \, \mathrm d s_y \, \mathrm d s_x \]

References