Laplace solver using indirect formulation (Neumann)#
keys: Laplace hypersingular operator, electrostatics, homogeneous Neumann bvp
from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
from libbem import *
from ngsolve import Projector, Preconditioner
from ngsolve.krylovspace import CG
Loading ngbem library
Consider the Neumann boundary value problem
The unique solution \(u\in H^1(\Omega)\) of the above bvp can be written in the following form (representation forumula)
Let’s carefully apply the Neumann trace \(\gamma_1\) to the representation formula and solve the resulting boundary integral equation for \(m \in H^{\frac12}(\Gamma)\) by discretisation of the following variational formulation:
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(2)
Define the finite element space for \(H^{\frac12}(\Gamma)\) for given geometry :
fesH1 = H1(mesh, order=1, definedon=mesh.Boundaries(".*"))
uH1,vH1 = fesH1.TnT()
print ("ndof H1 = ", fesH1.ndof)
ndof H1 = 173
Define the right hand side with given Neumann data \(u_1\):
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
rhs = LinearForm(u1exa*vH1.Trace()*ds(bonus_intorder=3)).Assemble()
The discretisation of the variational formulation leads to a system of linear equations, ie
where \(\mathrm{D}\) is the hypersingular operator and the stabilisation \((\mathrm D + \mathrm{S})\) is regular and symmetric.
vH1m1 = LinearForm(vH1*1*ds(bonus_intorder=3)).Assemble()
S = (BaseMatrix(Matrix(vH1m1.vec.Reshape(1))))@(BaseMatrix(Matrix(vH1m1.vec.Reshape(fesH1.ndof))))
m = GridFunction(fesH1)
pre = BilinearForm(uH1*vH1*ds).Assemble().mat.Inverse(freedofs=fesH1.FreeDofs())
with TaskManager():
D=HypersingularOperator(fesH1, intorder=12, leafsize=40, eta=3., eps=1e-11,
method="aca", testhmatrix=False)
CG(mat = D.mat+S, pre=pre, rhs = rhs.vec, sol=m.vec, tol=1e-8, maxsteps=200, initialize=False, printrates=True)
Draw (m, mesh, draw_vol=False, order=3);
CG iteration 1, residual = 1.0916223099897648
CG iteration 2, residual = 0.5870733996299382
CG iteration 3, residual = 0.35509458993776527
CG iteration 4, residual = 0.20439552637076824
CG iteration 5, residual = 0.10974300697070706
CG iteration 6, residual = 0.05522792544595984
CG iteration 7, residual = 0.025779925625166934
CG iteration 8, residual = 0.01249642797823975
CG iteration 9, residual = 0.008069904459550295
CG iteration 10, residual = 0.0029851625858649334
CG iteration 11, residual = 0.0010472414030450602
CG iteration 12, residual = 0.0003275190006576015
CG iteration 13, residual = 8.626799554097373e-05
CG iteration 14, residual = 1.9643810613358547e-05
CG iteration 15, residual = 2.9768861436102513e-06
CG iteration 16, residual = 8.178046730000429e-07
CG iteration 17, residual = 5.752611230527317e-07
CG iteration 18, residual = 2.864452917631025e-07
CG iteration 19, residual = 1.974886049253747e-07
CG iteration 20, residual = 1.1820019284473308e-07
CG iteration 21, residual = 7.213478334825916e-08
CG iteration 22, residual = 3.967268044732498e-08
CG iteration 23, residual = 1.9860987977673407e-08
CG iteration 24, residual = 1.0113232350425156e-08
Note: the hypersingular operator is implemented as follows
Details for instance in Numerische Näherungsverfahren für elliptische Randwertprobleme, p.127, p.259 (1st edition).