Laplace FEM-BEM Coupling

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 *

keys: exterior bvp, transmission problem, FEM-BEM coupling, capacitor, electrostatics

Laplace FEM-BEM Coupling#

We consider a plate capacitor problem and solve it with a FEM-BEM coupling.

\(\begin{array}{r rcl l} \textnormal{FEM domain}: & - \Delta u &=& 0, \quad &\Omega \,, \\[1ex] \textnormal{BEM domain}: & - \Delta u &=& 0, \quad&\Omega^c\,,\\[1ex] \textnormal{Dirichlet condition}: & \gamma_0 u &=& F_\Gamma & \Gamma_0\,, \\[1ex] \textnormal{coupling condition}: & [ \gamma_0 u] &=& 0 & \Gamma\,, \\[1ex] \textnormal{coupling condition}: & [ \gamma_1 u ] &=& 0 & \Gamma\,, \\[1ex] \textnormal{radiation condition}: & \lim\limits_{|x| \to \infty} u(x) &=& \mathcal O\left( \displaystyle{ \frac{1}{|x|} }\right)\,.\end{array}\)

\(\quad\quad\quad\)

 

largebox = Box ((-2,-2,-2), (2,2,2) )
eltop = Box ( (-1,-1,0.5), (1,1,1) )
elbot = Box ( (-1,-1,-1), (1,1,-0.5))

largebox.faces.name = "outer" # coupling boundary
eltop.faces.name = "topface" # Dirichlet boundary
elbot.faces.name = "botface" # Dirichlet boundary
#eltop.edges.hpref = 1
#elbot.edges.hpref = 1
eltop.faces.maxh=0.3
elbot.faces.maxh=0.3

shell = largebox-eltop-elbot # FEM domain 
shell.solids.name = "air"

mesh = shell.GenerateMesh(maxh=1)
# mesh.RefineHP(2)
ea = { "euler_angles" : (-67, 0, 110) } 
Draw (mesh, clipping={"x":1, "y":0, "z":0, "dist" : 1.1}, **ea);

The coupling needs finite elements in the product space \(H^1(\Omega)\times H^{-\frac12}(\Gamma)\) where \(H^1(\Omega)\) discretise the solution inside the FEM domain \(\Omega\) and \(H^{-\frac12}(\Gamma)\) is needed to resolve the unknown Neumann trace of \(u\) on the coupling boundary \(\Gamma\).

order = 2
fesH1 = H1(mesh, order=order, dirichlet="topface|botface") # H^1 conforming elements with H^(1/2) conforming elements on boundary
fesL2 = SurfaceL2(mesh, order=order-1, dual_mapping=True, definedon=mesh.Boundaries("outer")) # H^(-1/2) conforming elements
fes = fesH1 * fesL2
fes.ndof

8030

u,dudn = fes.TrialFunction()
v,dvdn = fes.TestFunction()

Set the given Dirichlet boundary conditions on the surfaces of the plates:

gfFGamma = GridFunction(fes)
gfFGamma.components[0].Set ( mesh.BoundaryCF( { "topface" : 1, "botface" : -1 }), BND)

FEM-BEM-Coupling

In the exterior domain \(\Omega^+\) the following representation formula for the solution \(u\) holds:

\[ x \in \Omega^c: \quad u(x) = - \int\limits_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{1}{\| x-y\|} } \, \gamma_1 u (y)\, \mathrm{d}s_y + \int\limits_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{ n(y) \cdot (x-y) }{\| x-y\|^3} } \, \gamma_0 u (y)\, \mathrm{d}s_y\,. \]

Within the weak formulation of the interior boundary value problem the transmission boundary conditions are incorporated leading to the FEM-BEM formulation (details below):

\[\begin{split} \begin{array}{rcl} \displaystyle { \int\limits_{\Omega}} \nabla u(x) \cdot \nabla v(x) \, \mathrm{d} x - \displaystyle{ \int\limits_{\Gamma} } \left(\left( \frac12 - K^\intercal\right) \,\gamma_1 u - D \, \gamma_0 u\right) \cdot \gamma_0 v(x) \, \mathrm{d}s &=& \displaystyle{ \int\limits_{\Omega} }f(x) \cdot v(x) \, \mathrm{d} x \\[2ex] \left( \begin{array}{ccc } A_{\Omega\Omega} & A_{\Omega\Gamma} & 0 \\ A_{\Gamma\Omega} & A_{\Gamma\Gamma } + D & -\frac12 M^\intercal + K^\intercal \\ 0 & -\frac12 M + K & -V \end{array}\right) \left( \begin{array}{c} u \\ \gamma_0 u \\ \gamma_1 u \end{array}\right) &=& \left( \begin{array}{c} F_{\Omega} \\ F_{\Gamma}\\ 0 \end{array}\right) \end{array} \end{split}\]
A = BilinearForm(grad(u)*grad(v)*dx, check_unused=False).Assemble() # FEM in Omega
F = LinearForm(fes).Assemble()

with TaskManager():
    V = LaplaceSL(dudn*ds("outer"))*dvdn*ds("outer", bonus_intorder=12) # single layer
    K = LaplaceDL(u*ds("outer"))*dvdn*ds("outer", bonus_intorder=12) # double layer
    D = LaplaceSL(curl(u)*ds("outer"))*curl(v)*ds("outer", bonus_intorder=12) # hypersingular 
    M = BilinearForm(u*dvdn*ds("outer"), check_unused=False).Assemble() #mass matrix

Note all matrix contriubutions are build with respect to the number of dofs in the product space where only the relevant entries. This means for instance that

  • \(A\) is of size \(\verb-fes.ndof x fes.ndof-\) and only \(\verb-H1fes x H1fes-\) dofs are filled.

  • \(V\) is of size \(\verb-fes.ndof x fes.ndof-\) and only \(\verb-(H1fes bndry dofs) x (H1fes bndry dofs)-\) dofs are filled.

  • \(K\) is of size \(\verb-fes.ndof x fes.ndof-\) and only \(\verb-fesL2 x (H1fes bndry dofs)-\) dofs are filled.

Due to the product space structure the variational formulation translates such as given by its formula:

lhs = A.mat+D.mat - (0.5*M.mat+K.mat).T - (0.5*M.mat+K.mat) - V.mat
rhs = (F.vec - A.mat * gfFGamma.vec).Evaluate()

Compute the solution of the coupled system (symmetric version) and look at the solution:

bfpre = BilinearForm(grad(u)*grad(v)*dx+1e-10*u*v*dx  + dudn*dvdn*ds("outer") ).Assemble()
pre = bfpre.mat.Inverse(freedofs=fes.FreeDofs(), inverse="sparsecholesky")
with TaskManager():
    sol_sym = CG(mat=lhs, rhs=rhs, pre=pre, tol=1e-8, maxsteps=400, printrates=False)

gfu = GridFunction(fes)
gfu.vec[:] = sol_sym + gfFGamma.vec 
Draw(gfu.components[1], clipping={"x" : 1, "y":0, "z":0, "dist":0.0, "function" : True }, **ea, order=2);

Details

In the exterior domain \(\Omega^+\) the following representation formula for the solution \(u\) holds:

\[ x \in \Omega^c: \quad u(x) = - \int\limits_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{1}{\| x-y\|} } \, \gamma_1 u (y)\, \mathrm{d}s_y + \int\limits_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{ n(y) \cdot (x-y) }{\| x-y\|^3} } \, \gamma_0 u (y)\, \mathrm{d}s_y\,, \]

and the unique traces of \(u\) are related by the Calderon projector

\[\begin{split} \left( \begin{array}{c} \gamma_0 u \\ \gamma_1 u \end{array}\right) = \left( \begin{array}{cc} -V & \frac12 + K \\ \frac12 - K^\intercal & -D \end{array} \right) \left( \begin{array}{c} \gamma_1 u \\ \gamma_0 u \end{array}\right) \end{split}\]

The weak formulation of the interior boundary value problem reads

\[\begin{split} \begin{array}{r rcl} \forall v \in H^1(\Omega): \quad \quad & \displaystyle { \int\limits_{\Omega}} \nabla u(x) \cdot \nabla v(x) \, \mathrm{d} x - \displaystyle{ \int\limits_{\Gamma} } \gamma_1 u(x) \cdot \gamma_0 v(x) \, \mathrm{d}s &=& \displaystyle{ \int\limits_{\Omega} }f(x) \cdot v(x) \, \mathrm{d} x \\[1ex] \textnormal{substitution of } \gamma_1 u: & \displaystyle { \int\limits_{\Omega}} \nabla u(x) \cdot \nabla v(x) \, \mathrm{d} x - \displaystyle{ \int\limits_{\Gamma} } \left(\left( \frac12 - K^\intercal\right) \,\gamma_1 u - D \, \gamma_0 u\right) \cdot \gamma_0 v(x) \, \mathrm{d}s &=& \displaystyle{ \int\limits_{\Omega} }f(x) \cdot v(x) \, \mathrm{d} x \end{array} \end{split}\]

The coupling condition \([\gamma_1 u]=0\) is now implicitly build in and we used the second equation of the Calderon projector. Adding the first equation of the Calderon projector is known as symmetric coupling formulation. (Note: without the substitution we obtain the Nedelec coupling. The Nedelec coupling is not symmetric.)

References: