This page was generated from wta/fembem.ipynb.

FEM-BEM Coupling¶

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

[1]:

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 *

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

:math:`begin{array}{r rcl l} textnormal{FEM domain}: & -: Delta u &=& 0, quad &Omega ,, \[1ex] textnormal{BEM domain}: & - Delta u &=& 0, quad&Omega^c,,\[1ex] textnorm
al{Dirichlet condition}: & gamma_0 u &=& u_0 & Gamma_0,, \[1ex] textnormal{coupling condition}: & [ gamma_0 u] &=& 0: & Gamma,, \[1ex] textnormal{coupling condition}: & [ gamma_1 u ] &=& 0 & Gamma,, \[1ex] textnormal{radiation con

dition}: & limlimits_{|x| to infty} u(x) &=& mathcal Oleft( displaystyle{ frac{1}{|x|} }right),.end{array}`

\(\quad\quad\quad\)

[2]:

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))

b1.name = "top"
b2.name = "bot"
b1.faces.name = "topface" # Dirichlet boundary
b2.faces.name = "botface" # Dirichlet boundary
b1.edges.maxh=0.1
b2.edges.maxh=0.1
shell = largebox-b1-b2 # FEM domain Omega
shell.name = "air"
largebox.faces.name = "outer" # coupling boundary

shape = shell
mesh = Mesh(OCCGeometry(shape).GenerateMesh(maxh=0.5))

Draw (mesh, clipping={"x":1, "y":0, "z":0, "dist" : 1.1});

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}\sigma_y + \int\limits_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{\langle n(y) , x-y\rangle }{\| x-y\|^3} } \, \gamma_0 u (y)\, \mathrm{d}\sigma_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}} \langle \nabla u(x), \nabla v(x)\rangle \, \mathrm{d} x - \displaystyle{ \int\limits_{\Gamma} } \gamma_1 u(x) \cdot \gamma_0 v(x) \, \mathrm{d}\sigma &=& \displaystyle{ \int\limits_{\Omega} }f(x) \cdot v(x) \, \mathrm{d} x \\[1ex] \textnormal{substitution of } \gamma_1 u: & \displaystyle { \int\limits_{\Omega}} \langle \nabla u(x), \nabla v(x)\rangle \, \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}\sigma &=& \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. To understand better the structure of the discretized system the dofs are splitted in degrees of freedom inside \(\Omega\) and those on the boundary \(\Gamma\).

\[\begin{split}\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{split}\]

Note: without the substitution we obtain the Nedelec coupling. The Nedelec coupling is not symmetric.

Generate the finite element space for \(H^1(\Omega)\) and set the given Dirichlet boundary conditions on the surfaces of the plates:

[3]:

order = 3
fesH1 = H1(mesh, order=order, dirichlet="topface|botface")
print ("H1-ndof = ", fesH1.ndof)
u,v = fesH1.TnT()
a = BilinearForm(grad(u)*grad(v)*dx).Assemble() # system matrix of variational formulation in Omega
dirtopface = GridFunction(fesH1)
dirbotface = GridFunction(fesH1)
dirtopface.Set(1, definedon=mesh.Boundaries("topface")) # boundary condition on upper plate
dirbotface.Set(-1, definedon=mesh.Boundaries("botface")) # boundary condition on lower plate
r = LinearForm(fesH1).Assemble()
f = r.vec - a.mat * (dirtopface.vec + dirbotface.vec)

H1-ndof =  138139

The finite element space \(\verb-fesH1-\) provides \(H^{\frac12}(\Gamma)\) conforming element to discretize the Dirichlet trace on the coupling boundary \(\Gamma\). However we still need \(H^{-\frac12}(\Gamma)\) conforming elements to discretize the Neumann trace of \(u\) on the coupling boundary. Here it is:

[4]:

fesL2 = SurfaceL2(mesh, order=order-1, dual_mapping=True, definedon=mesh.Boundaries("outer")) # H^(-1/2) conforming elements
f2 = LinearForm(fesL2).Assemble()  # 0-vector
print ("L2-ndof = ", fesL2.ndof)

L2-ndof =  8410

Generate the the single layer potential \(V\), double layer potential \(K\) and hypersingular operator \(D\):

[5]:

with TaskManager():
    V = SingleLayerPotentialOperator(fesL2, intorder=7)
    K = DoubleLayerPotentialOperator(fesH1, fesL2, trial_definedon=mesh.Boundaries("outer"),  test_definedon=mesh.Boundaries("outer"), intorder=7)
    D = HypersingularOperator(fesH1, definedon=mesh.Boundaries("outer"), intorder=7)
    M = BilinearForm(fesH1.TrialFunction()*fesL2.TestFunction().Trace()*ds("outer")).Assemble()

Setup the coupled system matrix and the right hand side:

[6]:

sym = BlockMatrix ([[a.mat+D.mat, (-0.5*M.mat+K.mat).T], [(-0.5*M.mat+K.mat), -V.mat]])
nonsym = BlockMatrix ([[a.mat, -M.mat.T ], [(-0.5*M.mat+K.mat), -V.mat]])
rhs = BlockVector([f, f2.vec])

l2mass = BilinearForm( fesL2.TrialFunction().Trace()*fesL2.TestFunction().Trace()*ds).Assemble()
astab = BilinearForm((grad(u)*grad(v) + 1e-10 * u * v)*dx).Assemble()
pre = BlockMatrix ([[astab.mat.Inverse(freedofs=fesH1.FreeDofs(), inverse="sparsecholesky"), None], \
                    [None, l2mass.mat.Inverse(freedofs=fesL2.FreeDofs())] ])

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

[7]:

with TaskManager(): # pajetrace=10**9):
    # sol_sym = CG(mat=sym, rhs=rhs, pre=pre, tol=1e-8, maxsteps=400, printrates=True)
    sol_sym = GMRes(A=sym, b=rhs, pre=pre, tol=1e-8, maxsteps=20, printrates=True)

GMRes iteration 1, residual = 49.554175657311035
GMRes iteration 2, residual = 8.415409608259521
GMRes iteration 3, residual = 3.5689498547253393
GMRes iteration 4, residual = 3.215550680466106
GMRes iteration 5, residual = 0.7095237505462254
GMRes iteration 6, residual = 0.46106165933017274
GMRes iteration 7, residual = 0.3328209200177098
GMRes iteration 8, residual = 0.3283674902952477
GMRes iteration 9, residual = 0.1602269194553177
GMRes iteration 10, residual = 0.13361636275522146
GMRes iteration 11, residual = 0.049133164911870364
GMRes iteration 12, residual = 0.047868273698335004
GMRes iteration 13, residual = 0.027484972365377308
GMRes iteration 14, residual = 0.024831234240233573
GMRes iteration 15, residual = 0.008778597303005724
GMRes iteration 16, residual = 0.008772578394420086
GMRes iteration 17, residual = 0.007676942906678719
GMRes iteration 18, residual = 0.007676169562480614
GMRes iteration 19, residual = 0.005131342614623752
GMRes iteration 20, residual = 0.004078891687600916
WARNING: GMRes did not converge to TOL

[8]:

gfu = GridFunction(fesH1)
gfu.vec[:] = sol_sym[0] + dirtopface.vec + dirbotface.vec
Draw(gfu, clipping={"y":1, "z":0, "dist":0.0, "function" : True });

References:

Principles of boundary element methods (symmetric coupling)
On the Coupling of Boundary Integral and Finite Element Methods (non-symmetric coupling)

[ ]:

[ ]: