This page was generated from unit-11.2-bem-Laplace/Laplace_DtN_direct.ipynb.
11.2.2. Dirichlet Laplace Direct Method¶
keys: homogeneous Dirichlet bvp, double and single layer potential, unknown Neumann data
[1]:
from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
from ngsolve.bem import *
from ngsolve.solvers import CG
Consider the Dirichlet boundary value problem
Let us choose the following ansatz for the solution \(u\in H^1(\Omega)\) (direct ansatz)
and solve for the Neumann data \(u_1 \in H^{-\frac12}(\Gamma)\) by the boundary element method, i.e.,
Define the geometry \(\Omega \subset \mathbb R^3\) and create a mesh:
[2]:
sp = Sphere( (0,0,0), 1)
mesh = Mesh( OCCGeometry(sp).GenerateMesh(maxh=0.2)).Curve(4)
Create the finite element spaces for \(H^{-\frac12}(\Gamma)\) and \(H^{\frac12}(\Gamma)\) according to the given mesh:
[3]:
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
Compute the interpolation of exact Dirichlet data \(u_0\) in finite element space:
[4]:
uexa = 1/ sqrt( (x-1)**2 + (y-1)**2 + (z-1)**2 )
u0 = GridFunction(fesH1)
u0.Interpolate (uexa)
Draw (u0, mesh, draw_vol=False, order=3);
The discretisation of the above variational formulation leads to a system of linear equations, ie
where the linear operators are as follows
\(\mathrm{V}\) is the single layer operator. \(\mathrm V\) is regular and symmetric.
\(\mathrm{M}\) is a mass matrix.
\(\mathrm{K}\) is the double layer operator.
We approximate the linear operators as hmatrices and solve for the Neumann data \(u_1\) with an iterative solver:
[5]:
u1 = GridFunction(fesL2)
pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()
with TaskManager():
V = SingleLayerPotentialOperator(fesL2, intorder=12)
M = BilinearForm( uH1 * v.Trace() * ds(bonus_intorder=3)).Assemble()
K = DoubleLayerPotentialOperator(fesH1, fesL2, intorder=12)
rhs = ( (0.5 * M.mat + K.mat)*u0.vec).Evaluate()
CG(mat = V.mat, pre=pre, rhs = rhs, sol=u1.vec, tol=1e-8, maxsteps=50, initialize=False, printrates=True)
CG iteration 1, residual = 0.2688244350833476
CG iteration 2, residual = 0.06553517738161971
CG iteration 3, residual = 0.01553700664638541
CG iteration 4, residual = 0.003677550794065572
CG iteration 5, residual = 0.0009473297863994385
CG iteration 6, residual = 0.0013803247805949576
CG iteration 7, residual = 0.00019659306335135563
CG iteration 8, residual = 5.781718855344875e-05
CG iteration 9, residual = 4.643518368928038e-05
CG iteration 10, residual = 3.399771904942353e-05
CG iteration 11, residual = 1.7058052504809743e-05
CG iteration 12, residual = 9.988909760937879e-06
CG iteration 13, residual = 4.669062226494854e-06
CG iteration 14, residual = 2.487876899453433e-06
CG iteration 15, residual = 5.52060778427861e-06
CG iteration 16, residual = 2.0923515652974644e-06
CG iteration 17, residual = 1.9107087754028916e-06
CG iteration 18, residual = 1.787457210484154e-06
CG iteration 19, residual = 9.145440804862722e-07
CG iteration 20, residual = 6.252906187712046e-07
CG iteration 21, residual = 5.046629337958127e-07
CG iteration 22, residual = 3.8856765207698746e-07
CG iteration 23, residual = 3.7945880145017083e-07
CG iteration 24, residual = 8.134174222074597e-07
CG iteration 25, residual = 2.0911860530728172e-07
CG iteration 26, residual = 1.008877821047814e-07
CG iteration 27, residual = 8.446741259131523e-08
CG iteration 28, residual = 5.0345322110237105e-08
CG iteration 29, residual = 5.520758267486288e-08
CG iteration 30, residual = 2.6787677538226487e-08
CG iteration 31, residual = 5.112962426457841e-08
CG iteration 32, residual = 3.460545120812853e-08
CG iteration 33, residual = 2.5817833881928984e-08
CG iteration 34, residual = 2.5905744662953545e-08
CG iteration 35, residual = 1.3211024818147368e-08
CG iteration 36, residual = 7.520154524559476e-09
CG iteration 37, residual = 4.500489919174989e-09
CG iteration 38, residual = 4.85949059612833e-09
CG iteration 39, residual = 3.081855310512987e-09
CG iteration 40, residual = 3.5451235455473237e-09
CG iteration 41, residual = 5.558273694481586e-09
CG iteration 42, residual = 2.0671917166251457e-09
[6]:
Draw (u1, mesh, draw_vol=False, order=3);
Let’s have a look at the exact Neumann data and compute the error of the numerical solution:
[7]:
graduexa = CF( (uexa.Diff(x), uexa.Diff(y), uexa.Diff(z)) )
n = specialcf.normal(3)
u1exa = graduexa*n
Draw (u1exa, mesh, draw_vol=False, order=3);
print ("L2-error =", sqrt (Integrate ( (u1exa-u1)**2, mesh, BND)))
L2-error = 0.0002493869498611152
[ ]:
[ ]: