Laplace Demo 2

Laplace Demo 2#

keys: homogeneous Neumann bvp, double layer potential ansatz, hypersingular operator, electrostatics

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
from ngbem import *
from ngsolve import Projector, Preconditioner
from ngsolve.krylovspace import CG

Loading ngbem library

Neumann Boundary Value Problem		Double Layer Potential		Variational Formulation
\( \left\{ \begin{array}{rcl l} -\Delta u &=& 0, \quad &\Omega \\ \gamma_1 u&=& u_1, \quad &\Gamma \end{array} \right. \)	\(\quad \Rightarrow \quad\)	\( u(x) = \mathrm{DL}(m) \)	\(\quad \Rightarrow \quad\)	\(\left\langle v, \gamma_1 \left(\mathrm{DL}(m)\right) \right\rangle_{-\frac12} = \left\langle u_1, v\right\rangle_{-\frac12} \)

				\(\mathrm{D} \, \mathrm{m} = \mathrm{M} \, \mathrm{u}_1 \)

NG-BEM Python interface	symbol	FE trial space	FE test space
`SingleLayerPotentialOperator`	\(\mathrm V \)	`SurfaceL2`	`SurfaceL2`
`DoubleLayerPotentialOperator`	\(\mathrm K \)	\(\gamma_0\) `H1`	`SurfaceL2`
`HypersingularOperator`	\(\mathrm D\)	\(\gamma_0\) `H1`	\(\gamma_0\) `H1`
`DoubleLayerPotentialOperator`	\(\mathrm K'\)	`SurfaceL2`	\(\gamma_0\) `H1`

Mesh

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

Trial and Test Spaces

fesH1 = H1(mesh, order=1, definedon=mesh.Boundaries(".*"))
uH1,vH1 = fesH1.TnT()

Right Hand Side \(\;\mathrm{M}\mathrm{1}_0\)

uexa = 1/ sqrt( (x-1)**2 + (y-1)**2 + (z-1)**2 )
graduexa = CF( (uexa.Diff(x), uexa.Diff(y), uexa.Diff(z)) )

n = specialcf.normal(3)
u1exa = graduexa*n
Mu1 = LinearForm(u1exa*vH1.Trace()*ds(bonus_intorder=3)).Assemble()

System Matrix \(\, \mathrm{D}\)

D=HypersingularOperator(fesH1, intorder=12, leafsize=40, eta=3., eps=1e-11, method="aca")

Stabilization Matrix \(\, \mathrm{S}\)

vH1m1 = LinearForm(vH1*1*ds(bonus_intorder=3)).Assemble()
S = (BaseMatrix(Matrix(vH1m1.vec.Reshape(1))))@(BaseMatrix(Matrix(vH1m1.vec.Reshape(fesH1.ndof))))

Solve \(\; \left(\mathrm{D} + \mathrm{S}\right) \mathrm{m} = \mathrm{M}\mathrm{u}_1\)

m = GridFunction(fesH1)
pre = BilinearForm(uH1*vH1*ds).Assemble().mat.Inverse(freedofs=fesH1.FreeDofs()) 
with TaskManager(): 
    CG(mat = D.mat+S, pre=pre, rhs = Mu1.vec, sol=m.vec, tol=1e-8, maxsteps=200, initialize=False, printrates=False)
Draw (m, mesh, draw_vol=False, order=3);

Note: Details for instance in Numerische Näherungsverfahren für elliptische Randwertprobleme, p.127, p.259 (1st edition).