This page was generated from wta/fembem.ipynb.

FEM-BEM Coupling

The ngbem boundary element addon project initiated by Lucy Weggeler (see https://weggler.github.io/ngbem/intro.html) is now partly integrated into core NGSolve. Find a short and sweet introduction to the boundary element method there.

In this demo we simulate a plate capacitor on an unbounded domain.

[1]:

from ngsolve import *
from netgen.occ import *
from ngsolve.solvers import GMRes
from ngsolve.webgui import Draw
from ngsolve.bem import *

[2]:

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

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

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

On the exterior domain \(\Omega^c\), the solution can be expressed by the representation formula:

\[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\,,\]

where \(\gamma_0 u = u\) and \(\gamma_1 u = \frac{\partial u}{\partial n}\) are Dirichlet and Neumann traces. These traces 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 \(V\), \(K\) are the single layer and double layer potential operators, and \(D\) is the hypersingular operator.

On the FEM domain we have the variational formulation

\[\int_{\Omega_\text{FEM}} \nabla u \nabla v \, dx - \int_\Gamma \gamma_1 u v \, ds = 0 \qquad \forall \, v \in H^1(\Omega_\text{FEM})\]

We use Calderon’s represenataion formula for the Neumann trace:

\[\int_{\Omega_\text{FEM}} \nabla u \nabla v \, dx - \int_\Gamma \left( \left( \tfrac{1}{2} - K^\intercal\right) \,\gamma_1 u - D \, \gamma_0 u\right) v = 0 \qquad \forall \, v \in H^1(\Omega_\text{FEM})\]

To get a closed system, we use also the first equation of the Calderon equations. To see the structure of the discretized system, the dofs are split into degrees of freedom inside \(\Omega\), and those on the boundary \(\Gamma\). The FEM matrix \(A\) is split accordingly. We see, the coupled system is symmetric, but indefinite:

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

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

[3]:

order = 4
fesH1 = H1(mesh, order=order, dirichlet="topface|botface")
print ("H1-ndof = ", fesH1.ndof)

H1-ndof =  90703

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"))
print ("L2-ndof = ", fesL2.ndof)

L2-ndof =  4035

[5]:

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

a = BilinearForm(grad(u)*grad(v)*dx, check_unused=False).Assemble()

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

f = LinearForm(fes).Assemble()
res = (f.vec - a.mat * gfudir.vec).Evaluate()

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

[6]:

n = specialcf.normal(3)
with TaskManager():
    V = LaplaceSL(dudn*ds("outer"))*dvdn*ds("outer")
    K = LaplaceDL(u*ds("outer"))*dvdn*ds("outer")
    D = LaplaceSL(Cross(grad(u).Trace(),n)*ds("outer"))*Cross(grad(v).Trace(),n)*ds("outer")
    M = BilinearForm(u*dvdn*ds("outer"), check_unused=False).Assemble()

Setup the coupled system matrix and the right hand side:

[7]:

sym = a.mat+D.mat - (0.5*M.mat+K.mat).T - (0.5*M.mat+K.mat) - V.mat
rhs = res

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

Compute the solution of the coupled system:

[8]:

with TaskManager():
    sol_sym = GMRes(A=sym, b=rhs, pre=pre, tol=1e-6, maxsteps=200, printrates=True)

GMRes iteration 1, residual = 47.943299270587204
GMRes iteration 2, residual = 10.539369746936108
GMRes iteration 3, residual = 2.6115871565308733
GMRes iteration 4, residual = 1.9195757191363212
GMRes iteration 5, residual = 0.4166916580572418
GMRes iteration 6, residual = 0.39297352574789196
GMRes iteration 7, residual = 0.17795761889271197
GMRes iteration 8, residual = 0.1381652601529096
GMRes iteration 9, residual = 0.08671805759934241
GMRes iteration 10, residual = 0.04706002593758194
GMRes iteration 11, residual = 0.04534479530149084
GMRes iteration 12, residual = 0.04442191913711668
GMRes iteration 13, residual = 0.01670150254309103
GMRes iteration 14, residual = 0.012683670900487982
GMRes iteration 15, residual = 0.011115768703880067
GMRes iteration 16, residual = 0.00705122598313693
GMRes iteration 17, residual = 0.0067145802604505435
GMRes iteration 18, residual = 0.005947766472067566
GMRes iteration 19, residual = 0.004793409918545511
GMRes iteration 20, residual = 0.0035486000830410396
GMRes iteration 21, residual = 0.0031764124769938173
GMRes iteration 22, residual = 0.0028262168204223966
GMRes iteration 23, residual = 0.0024657041484887894
GMRes iteration 24, residual = 0.0018632042929135734
GMRes iteration 25, residual = 0.0018471765308543317
GMRes iteration 26, residual = 0.0013478969077251636
GMRes iteration 27, residual = 0.0013444294385622885
GMRes iteration 28, residual = 0.0010720995980173698
GMRes iteration 29, residual = 0.0010545592009086762
GMRes iteration 30, residual = 0.0007620541686120474
GMRes iteration 31, residual = 0.0006873897894447038
GMRes iteration 32, residual = 0.0006388026389506508
GMRes iteration 33, residual = 0.0005522682691807808
GMRes iteration 34, residual = 0.00048353717750979314
GMRes iteration 35, residual = 0.00044792987834669097
GMRes iteration 36, residual = 0.0003647407064139316
GMRes iteration 37, residual = 0.0003467884081898194
GMRes iteration 38, residual = 0.0003021956208435088
GMRes iteration 39, residual = 0.0002885819324165296
GMRes iteration 40, residual = 0.00023060160549384884
GMRes iteration 41, residual = 0.00019995191132517228
GMRes iteration 42, residual = 0.00019138778671354067
GMRes iteration 43, residual = 0.00015891304880227188
GMRes iteration 44, residual = 0.00015494149385500544
GMRes iteration 45, residual = 0.0001238328487838517
GMRes iteration 46, residual = 0.0001238119182382527
GMRes iteration 47, residual = 9.622952550442424e-05
GMRes iteration 48, residual = 9.34971377512853e-05
GMRes iteration 49, residual = 7.689469193091164e-05
GMRes iteration 50, residual = 7.458496855696329e-05
GMRes iteration 51, residual = 6.382336305437319e-05
GMRes iteration 52, residual = 6.236932620184402e-05
GMRes iteration 53, residual = 4.9562812531511316e-05
GMRes iteration 54, residual = 4.124180340964926e-05
GMRes iteration 55, residual = 4.0183431261235996e-05
GMRes iteration 56, residual = 3.2049208552207483e-05
GMRes iteration 57, residual = 3.184712755060231e-05
GMRes iteration 58, residual = 2.4094093883075078e-05
GMRes iteration 59, residual = 2.4009479477308206e-05
GMRes iteration 60, residual = 2.11991217270018e-05
GMRes iteration 61, residual = 2.1157691656665414e-05
GMRes iteration 62, residual = 1.727389165560683e-05
GMRes iteration 63, residual = 1.7055554352732797e-05
GMRes iteration 64, residual = 1.3764608786215311e-05
GMRes iteration 65, residual = 1.2844403571977285e-05
GMRes iteration 66, residual = 1.1069384621956165e-05
GMRes iteration 67, residual = 8.716932637317827e-06
GMRes iteration 68, residual = 8.571988881485973e-06
GMRes iteration 69, residual = 6.935376742891925e-06
GMRes iteration 70, residual = 6.934584630748257e-06
GMRes iteration 71, residual = 5.211342934246865e-06
GMRes iteration 72, residual = 5.074443689421805e-06
GMRes iteration 73, residual = 4.230290151022309e-06
GMRes iteration 74, residual = 3.663195241258226e-06
GMRes iteration 75, residual = 3.5897824592906982e-06
GMRes iteration 76, residual = 2.815391008038292e-06
GMRes iteration 77, residual = 2.8153854351511085e-06
GMRes iteration 78, residual = 2.267174778751487e-06
GMRes iteration 79, residual = 2.2424486925513235e-06
GMRes iteration 80, residual = 1.9715398474103937e-06
GMRes iteration 81, residual = 1.9167706956946506e-06
GMRes iteration 82, residual = 1.7678126981659424e-06
GMRes iteration 83, residual = 1.368533596610436e-06
GMRes iteration 84, residual = 1.3588969690414093e-06
GMRes iteration 85, residual = 1.0769583639778505e-06
GMRes iteration 86, residual = 1.0765771466164101e-06
GMRes iteration 87, residual = 8.757140339837792e-07

[9]:

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

The Neumann data:

[10]:

Draw (gfu.components[1], **ea);

References:

[ ]: