Dirichlet Laplace Direct 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, double layer potential \(K\) and single layer potential \(V\)

Dirichlet Laplace Direct Method#

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) }} - \underbrace{ \int\limits_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{\langle n(y) , x-y\rangle }{\| x-y\|^3} } \, m(y)\, \mathrm{d}s_y}_{\displaystyle{ \mathrm{LaplaceDL}(m) }}\,,\]

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

\[ \mathrm V \, \mathrm j= \left( \frac12 \mathrm M + \mathrm K\right) \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.2)).Curve(4)

fesL2 = SurfaceL2(mesh, order=3, dual_mapping=True)
u,v = fesL2.TnT()
fesH1 = H1(mesh, order=4)
uH1,vH1 = fesH1.TnT()
print ("ndofL2 = ", fesL2.ndof, "ndof H1 = ", fesH1.ndof)

ndofL2 =  7220 ndof H1 =  21241

Consider as Dirichlet data \(m = \displaystyle \frac{1}{\|x-x_0\|}\) and compute its interpolation in \(H^{\frac12}(\Gamma)\):

uexa = 1/ sqrt( (x-1)**2 + (y-1)**2 + (z-1)**2 )
m = GridFunction(fesH1)
m.Interpolate (uexa)
Draw (m, mesh, draw_vol=False, order=3)

BaseWebGuiScene

Boundary Integral Equation

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

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

The discretisation of the above variational formulation is the boundary element method and it leads to a system of linear equations, ie

\[ \mathrm{V} \, \mathrm{j} = \left( \frac12 \,\mathrm{M} + \mathrm{K} \right) \, \mathrm{m}\,, \]

where

\(\mathrm{V}\) is the single layer operator - it is regular and symmetric.
\(\mathrm{K}\) is the double layer operator.
\(\mathrm{M}\) is a mass matrix.

In NGSBEM, \(V\) and \(K\) are assembled by

V = LaplaceSL(u*ds)*v*ds
K = LaplaceDL(w*ds)*v*ds

Compute the linear operators and solve the linear system of equations for the Neumann data \(j\)

j = GridFunction(fesL2)
pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()
with TaskManager():
    V = LaplaceSL(u*ds)*v*ds
    K = LaplaceDL(uH1*ds)*v*ds
    M = BilinearForm(uH1*v*ds(bonus_intorder=3)).Assemble()
    
    rhs = ( (0.5 * M.mat + K.mat)*m.vec).Evaluate()
    CG(mat = V.mat, pre=pre, rhs = rhs, sol=j.vec, tol=1e-8, maxsteps=50, 
       initialize=False, printrates=False)

Draw (j, mesh, draw_vol=False, order=3);

Let’s have a look at the exact Neumann data and compute the error of the numerical solution:

graduexa = CF( (uexa.Diff(x), uexa.Diff(y), uexa.Diff(z)) )
n = specialcf.normal(3)
jexa = graduexa*n
#Draw (jexa, mesh, draw_vol=False, order=3);
print ("L2-error =", sqrt (Integrate ( (jexa-j)**2, mesh, BND)))

L2-error = 0.0002505947168449771

Evaluation of Representation formula

Use the represenation formula \((1)\) and evaluate the solution \(u\) on a screen inside \(\Omega\)

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

LSPotential = LaplaceSL( u*ds )
LDPotential = LaplaceDL( uH1*ds )

repformula = LSPotential(j) - LDPotential(m)
fes_screen = H1(vismesh, order=5)
gf_screen = GridFunction(fes_screen)
with TaskManager():
    gf_screen.Set(repformula, definedon=vismesh.Boundaries("screen"))
Draw(gf_screen)

BaseWebGuiScene

Have a look at the error in the solution \(u\) on the screen:

Draw(repformula - uexa, vismesh)

BaseWebGuiScene

type(LSPotential)

ngsolve.bem.PotentialOperator

Details: The Single and Double Layer Potential Operators

For any trial functions \(u_j, w_j\) and test function \(v_i\), the layer potential operators are implemented as follows

\[\begin{split} \begin{array}{rcl} V_{ij} &=& \displaystyle \frac{1}{4\pi}\, \displaystyle \int\limits_\Gamma \int\limits_\Gamma \frac{1}{\|x-y\|} \, u_j(x) \, v_i(y) \, \mathrm{d} s_{y} \, \mathrm{d} s_x\,, \\[2ex] K_{ij} &=& \displaystyle \frac{1}{4\pi} \int\limits_\Gamma\int\limits_\Gamma n_y\cdot \nabla_y \frac{1}{\|x-y\|} \, w_j(x) \, v_i(y) \, \mathrm{d} s_{y} \, \mathrm{d} s_x \end{array}\end{split}\]