Laplace solver using direct formulation (Neumann)#
keys: Laplace hypersingular operator, electrostatics, homogeneous Neumann bvp, Dirichlet data
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)\) 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 the Dirichlet trace \(u_0\) of \(u\) 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.2)).Curve(2)
Define conforming finite element spaces for \(H^{\frac12}(\Gamma)\) and \(H^{-\frac12}(\Gamma)\) respectively:
fesL2 = SurfaceL2(mesh, order=2, dual_mapping=False) # TODO: check interpolation with option dual_mapping=True
u,v = fesL2.TnT()
fesH1 = H1(mesh, order=2, definedon=mesh.Boundaries(".*"))
uH1,vH1 = fesH1.TnT()
print ("ndofL2 = ", fesL2.ndof, "ndof H1 = ", fesH1.ndof)
ndofL2 = 4332 ndof H1 = 1600
Project the given Neumann data \(u_1\) in finite element space and have a look at the boundary data:
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)
u1 = GridFunction(fesL2)
u1.Interpolate(graduexa*n, definedon=mesh.Boundaries(".*"))
Draw (u1, mesh, draw_vol=False, order=3);
The discretisation of the variational formulation leads to a system of linear equations, ie
where the linear operators are as follows
\(\mathrm{D}\) is the hypersingular operator and the stabilisation \((\mathrm D + \mathrm{S})\) is regular and symmetric.
\(\mathrm{M}\) is a mass matrix.
\(\mathrm{K}'\) is the adjoint double layer operator.
The stabilisation matrix \(S\) is needed to cope with the kernel of the hypersingular operator. The following stabilisation matrix is used here:
where \(v_i, v_j\) being \(H^{\frac12}(\Gamma)\)-conforming shape functions.
vH1m1 = LinearForm(vH1*1*ds(bonus_intorder=3)).Assemble()
S = (BaseMatrix(Matrix(vH1m1.vec.Reshape(1))))@(BaseMatrix(Matrix(vH1m1.vec.Reshape(fesH1.ndof))))
Now we assemble the stabilised system matrix and the right hand side and solve for the Dirichlet data \(u_0\) with an iterative solver:
u0 = GridFunction(fesH1)
pre = BilinearForm(uH1*vH1*ds).Assemble().mat.Inverse(freedofs=fesH1.FreeDofs())
with TaskManager(): # pajetrace=1000*1000*1000):
D = HypersingularOperator(fesH1, intorder=12, leafsize=40, eta=3., eps=1e-4, method="aca", testhmatrix=False)
M = BilinearForm( v.Trace() * uH1 * ds(bonus_intorder=3)).Assemble()
Kt = DoubleLayerPotentialOperator(fesL2, fesH1, intorder=12, leafsize=40, eta=3., eps=1e-4, method="aca")
rhs = ( (0.5 * M.mat.T - Kt.mat) * u1.vec ).Evaluate()
CG(mat=D.mat+S, pre=pre, rhs=rhs, sol=u0.vec, tol=1e-8, maxsteps=200, initialize=False, printrates=True)
CG iteration 1, residual = 0.6734119976551735
CG iteration 2, residual = 0.3570708359436435
CG iteration 3, residual = 0.20888851335800046
CG iteration 4, residual = 0.1706301171838047
CG iteration 5, residual = 0.1292250234421183
CG iteration 6, residual = 0.10837761567750001
CG iteration 7, residual = 0.06658980509446069
CG iteration 8, residual = 0.05390070177634602
CG iteration 9, residual = 0.03616839125884154
CG iteration 10, residual = 0.023747331088605498
CG iteration 11, residual = 0.017135975646323597
CG iteration 12, residual = 0.010529013246749049
CG iteration 13, residual = 0.006843245189196467
CG iteration 14, residual = 0.004377148034700796
CG iteration 15, residual = 0.002792183423350529
CG iteration 16, residual = 0.0017274928252797245
CG iteration 17, residual = 0.001004732706787392
CG iteration 18, residual = 0.0005839403720913693
CG iteration 19, residual = 0.00034101898173912544
CG iteration 20, residual = 0.00019007063675770494
CG iteration 21, residual = 0.00010675896994169914
CG iteration 22, residual = 5.8059772565655696e-05
CG iteration 23, residual = 3.149811493549551e-05
CG iteration 24, residual = 1.7107557048998617e-05
CG iteration 25, residual = 9.037396099049056e-06
CG iteration 26, residual = 4.626967507663484e-06
CG iteration 27, residual = 2.4760648866817508e-06
CG iteration 28, residual = 1.33189265203938e-06
CG iteration 29, residual = 8.230311535063995e-07
CG iteration 30, residual = 6.745210316117556e-07
CG iteration 31, residual = 7.082765813255227e-07
CG iteration 32, residual = 7.260965308677649e-07
CG iteration 33, residual = 5.508925208061692e-07
CG iteration 34, residual = 2.9975673881111267e-07
CG iteration 35, residual = 1.4416911434107242e-07
CG iteration 36, residual = 6.545008867426668e-08
CG iteration 37, residual = 3.0354147814096156e-08
CG iteration 38, residual = 1.5207593599889786e-08
CG iteration 39, residual = 9.546096972222183e-09
CG iteration 40, residual = 8.580269058852484e-09
CG iteration 41, residual = 8.208573615402671e-09
CG iteration 42, residual = 5.6662151786155515e-09
Let’s have a look at the numerical and exact Dirichlet data and compare them quantitatively:
Draw (u0, mesh, draw_vol=False, order=3);
u0exa = GridFunction(fesH1)
u0exa.Interpolate (uexa, definedon=mesh.Boundaries(".*"))
Draw (u0exa, mesh, draw_vol=False, order=3);
print ("L2 error of surface gradients =", sqrt (Integrate ( (grad(u0exa)-grad(u0))**2, mesh, BND)))
L2 error of surface gradients = 0.004797701047764322