BEM Laplace Capacitor

from ngsolve import *
from netgen.occ import *
import netgen.meshing as meshing
from ngsolve.krylovspace import CG, GMRes
from ngsolve.webgui import Draw
from ngsolve.bem import *

BEM Laplace Capacitor#

We consider a plate capacitor problem and solve it with a BEM. The setting and notation looks like this:


$$\begin{array}{rcl l} - \Delta u &=& 0, \quad &\mathrm{in} \; \Omega^c \,, \\[1ex] \gamma_0 u &=& m\,, & \mathrm{on}\;\Gamma\,, \\[1ex] \lim\limits_{\|x\| \to \infty} u(x) &=& \mathcal O\left( \displaystyle{ \frac{1}{\|x\|} }\right)\,, & \|x\|\to \infty \,. \end{array} $$	$\quad\quad\quad$

largebox = Box ((-2,-2,-2), (2,2,2) )
b1 = Box ( (-1,-1,0.5), (1,1,1) )
b2 = Box ( (-1,-1,-1), (1,1,-0.5))

largebox.faces.name = "outer" 
b1.faces.name = "top" # part of Gamma
b2.faces.name = "bot" # part of Gamma
shell = largebox-b1-b2 # Omega^c 
shape = Compound([b1,b2])
mesh = Mesh(OCCGeometry(shape).GenerateMesh(maxh=1))
Draw (mesh);

order = 3
fesH1 = H1(mesh, order=order, definedon=mesh.Boundaries(".*")) # trace space, i.e. H^(1/2)(Gamma) conforming elements 
uH1,vH1 = fesH1.TnT()
fesL2 = SurfaceL2(mesh, order=order-1, dual_mapping=True) # trace space, i.e. H^(-1/2)(Gamma) conforming elements 
u,v = fesL2.TnT()
print ("L2-ndof = ", fesL2.ndof, "LH1-ndof = ", fesH1.ndof)

L2-ndof =  576 LH1-ndof =  437

utop = GridFunction(fesH1)
utop.Interpolate(1, definedon=mesh.Boundaries("top"))
ubot = GridFunction(fesH1)
ubot.Interpolate(-1, definedon=mesh.Boundaries("bot"))
m = utop.vec + ubot.vec
Draw (utop + ubot, mesh, draw_vol=False, order=3)

BaseWebGuiScene

Compute boundary integral operators and the mass matrix for right hand side. Solve the linear system of equations:

j = GridFunction(fesL2)
with TaskManager():
    SLPotential = LaplaceSL(u*ds(bonus_intorder=3))
    DLPotential = LaplaceDL(uH1*ds(bonus_intorder=3))
    V = SLPotential*v*ds(bonus_intorder=3)
    K = DLPotential*v*ds(bonus_intorder=3)
    M = BilinearForm(uH1*v*ds(bonus_intorder=3)).Assemble()

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

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

Evaluate the solution on a vertical screen

screen = WorkPlane(Axes((0,0,0), X, Z)).RectangleC(4, 4).Face() - Box((-1.1,-1.1,0.4), (1.1,1.1,1.1)) - Box((-1.1,-1.1,-1.1), (1.1,1.1,-0.4))
mesh_screen = Mesh(OCCGeometry(screen).GenerateMesh(maxh=0.25)).Curve(1)
fes_screen = H1(mesh_screen, order=3)
gf_screen = GridFunction(fes_screen)
print ("ndofscreen=", fes_screen.ndof)
with TaskManager():
    gf_screen.Set(-SLPotential(j)+DLPotential(utop)+DLPotential(ubot), definedon=mesh_screen.Boundaries(".*"), dual=False)
Draw (gf_screen)

ndofscreen= 2210

BaseWebGuiScene

Testing the BEM solution with FEM

Evaluate the solution $u$ on a box that wraps the place capacitor. We use this data as boundary data for the FEM computation:

mesh_shell = Mesh(OCCGeometry(shell).GenerateMesh(maxh=0.5)).Curve(1)
fes_shell = H1(mesh_shell, order=3)
gf_shell = GridFunction(fes_shell)
with TaskManager():
    gf_shell.Set(-SLPotential(j)+DLPotential(utop)+DLPotential(ubot), definedon=mesh_shell.Boundaries("outer"), dual=False)
Draw(gf_shell)

BaseWebGuiScene

Solve the FEM variational formulation with given Dirichlet conditions on the plates and on the ficticious box:

fesH1d = H1(mesh_shell, order=order, dirichlet="outer|top|bot")
ud,vd = fesH1d.TnT()
ad = BilinearForm(grad(ud)*grad(vd)*dx).Assemble()
fd = LinearForm(fesH1d).Assemble()
gfud = GridFunction(fesH1d)

utop = GridFunction(fesH1d)
utop.Interpolate(1, definedon=mesh_shell.Boundaries("top"))
ubot = GridFunction(fesH1d)
ubot.Interpolate(-1, definedon=mesh_shell.Boundaries("bot"))

Have a look at the FEM solution on a clipping plane through the box:

r = fd.vec - ad.mat * (gf_shell.vec + utop.vec + ubot.vec)
gfud.vec.data = gf_shell.vec + utop.vec + ubot.vec
gfud.vec.data += ad.mat.Inverse(freedofs=fesH1d.FreeDofs()) * r
Draw(gfud, clipping={"y":1, "z":0, "dist":0.0, "function" : True }); # turn the clippling plane on in the control panel