Laplace Plate Capacitor

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

\( \left\{ \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} \right. \quad\quad\quad\)

The exterior solution is represented by the representation formula

\[ u(x)=-Vj(x)+Km(x), \qquad Vj(x)=\int_\Gamma G(x,y)j(y)\,dy, \qquad Km(x)=\int_\Gamma \frac{\partial G(x,y)}{\partial n_y}m(y)\,dy, \]

where the Dirichlet trace \(m=\gamma_0 u\) is prescribed and the Neumann trace \(j\) is unknown. Applying the trace operator we need to solve

\[ Vj=\left(K-\frac12 I\right)m. \]

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 *

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
b1.edges.hpref = 1
b2.edges.hpref = 1

shell = largebox-b1-b2 # Omega^c 
shape = Compound([b1,b2])
mesh = shape.GenerateMesh(maxh=1)

order = 3
fesH1 = H1(mesh, order=order, definedon=mesh.Boundaries(".*")) 
uH1,vH1 = fesH1.TnT()
fesL2 = SurfaceL2(mesh, order=order-1, dual_mapping=True) 
u,v = fesL2.TnT()

m = GridFunction(fesH1)
m[BND("top")] = 1
m[BND("bot")] = -1
Draw (m, mesh, draw_vol=False, order=3, euler_angles=(90,0,0) )

WebGLScene

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

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

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

Evaluation of the Solution

screen = WorkPlane(Axes((0,0,0), X, Z)).RectangleC(4, 4).Face() - Box((-1,-1,0.5), (1,1,1)) - Box((-1,-1,-1), (1,1,-0.5))
mesh_screen = Mesh(OCCGeometry(screen).GenerateMesh(maxh=0.15)).Curve(1)
fes_screen = H1(mesh_screen, order=3)
gf_screen = GridFunction(fes_screen)
with TaskManager():
    gf_screen.Set(-SLPotential(j) + DLPotential(m), definedon=mesh_screen.Boundaries(".*"), dual=False)
ea = {"euler_angles": (0,90,180)} 
Draw(gf_screen, **ea)

WebGLScene