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)
u,v = fesH1.TnT()
a = BilinearForm(grad(u)*grad(v)*dx).Assemble()

gfudir = GridFunction(fesH1)
gfudir.Set ( mesh.BoundaryCF( { "topface" : 1, "botface" : -1 }), BND)

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

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"))
f2 = LinearForm(fesL2).Assemble()  # 0-vector
print ("L2-ndof = ", fesL2.ndof)

L2-ndof =  4035

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

[5]:

with TaskManager():
    V = SingleLayerPotentialOperator(fesL2, intorder=2*order)
    K = DoubleLayerPotentialOperator(fesH1, fesL2, trial_definedon=mesh.Boundaries("outer"),  test_definedon=mesh.Boundaries("outer"), intorder=2*order)
    D = HypersingularOperator(fesH1, definedon=mesh.Boundaries("outer"), intorder=2*order)
    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]])
rhs = BlockVector([res, 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:

[7]:

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

GMRes iteration 1, residual = 47.94329927058657
GMRes iteration 2, residual = 11.373090191398727
GMRes iteration 3, residual = 6.709311492494004
GMRes iteration 4, residual = 6.495804108604679
GMRes iteration 5, residual = 1.4328881666094593
GMRes iteration 6, residual = 0.8253014999039325
GMRes iteration 7, residual = 0.6024935444753226
GMRes iteration 8, residual = 0.3865476425807203
GMRes iteration 9, residual = 0.2753852464823687
GMRes iteration 10, residual = 0.25092009338551885
GMRes iteration 11, residual = 0.10003651813637562
GMRes iteration 12, residual = 0.09550113470010761
GMRes iteration 13, residual = 0.0425738660489936
GMRes iteration 14, residual = 0.03054274970551427
GMRes iteration 15, residual = 0.023933863792177674
GMRes iteration 16, residual = 0.023933859774625275
GMRes iteration 17, residual = 0.018760693724983762
GMRes iteration 18, residual = 0.016374848841108795
GMRes iteration 19, residual = 0.014108850687884695
GMRes iteration 20, residual = 0.009924660928686227
GMRes iteration 21, residual = 0.009920782010763407
GMRes iteration 22, residual = 0.007811734811770193
GMRes iteration 23, residual = 0.007766703224049966
GMRes iteration 24, residual = 0.006225361576706692
GMRes iteration 25, residual = 0.00620356903965731
GMRes iteration 26, residual = 0.003999570702980933
GMRes iteration 27, residual = 0.003884270768787932
GMRes iteration 28, residual = 0.0035496185564046195
GMRes iteration 29, residual = 0.0035459215913412182
GMRes iteration 30, residual = 0.003085763331710947
GMRes iteration 31, residual = 0.0030514890186890725
GMRes iteration 32, residual = 0.00271164456506776
GMRes iteration 33, residual = 0.002699189964614515
GMRes iteration 34, residual = 0.0023762747343303854
GMRes iteration 35, residual = 0.0021121224806499965
GMRes iteration 36, residual = 0.002007320868617572
GMRes iteration 37, residual = 0.001868629955303915
GMRes iteration 38, residual = 0.001848700470171609
GMRes iteration 39, residual = 0.0016043458161723713
GMRes iteration 40, residual = 0.0015902550156973335
GMRes iteration 41, residual = 0.0014588203818939125
GMRes iteration 42, residual = 0.001422038825446279
GMRes iteration 43, residual = 0.0012612978312431388
GMRes iteration 44, residual = 0.0011942142541429028
GMRes iteration 45, residual = 0.0010979317896149868
GMRes iteration 46, residual = 0.0010972029934491836
GMRes iteration 47, residual = 0.0009820600127087656
GMRes iteration 48, residual = 0.000865206956819346
GMRes iteration 49, residual = 0.0008166912482141482
GMRes iteration 50, residual = 0.0007499606389623837
GMRes iteration 51, residual = 0.0007493421049440593
GMRes iteration 52, residual = 0.0006249378461679438
GMRes iteration 53, residual = 0.0005734273613161808
GMRes iteration 54, residual = 0.0005321357107001589
GMRes iteration 55, residual = 0.0005132383405642408
GMRes iteration 56, residual = 0.00045951443359698374
GMRes iteration 57, residual = 0.00034788838986366655
GMRes iteration 58, residual = 0.0003431734861262124
GMRes iteration 59, residual = 0.00030745814162116427
GMRes iteration 60, residual = 0.0002995045179221391
GMRes iteration 61, residual = 0.0002327669351329046
GMRes iteration 62, residual = 0.00022878482842846092
GMRes iteration 63, residual = 0.00020593184418826095
GMRes iteration 64, residual = 0.00019438960881182753
GMRes iteration 65, residual = 0.00018032109467483107
GMRes iteration 66, residual = 0.00015370517570923407
GMRes iteration 67, residual = 0.00015320432036099453
GMRes iteration 68, residual = 0.00013170137383002658
GMRes iteration 69, residual = 0.00013039964238399578
GMRes iteration 70, residual = 0.00010479061032327055
GMRes iteration 71, residual = 9.437471995365411e-05
GMRes iteration 72, residual = 8.938739605848492e-05
GMRes iteration 73, residual = 7.940775356196922e-05
GMRes iteration 74, residual = 7.849790161869342e-05
GMRes iteration 75, residual = 6.686224944197573e-05
GMRes iteration 76, residual = 6.685674132566134e-05
GMRes iteration 77, residual = 5.5331369934307867e-05
GMRes iteration 78, residual = 5.466990317746458e-05
GMRes iteration 79, residual = 5.1504939394206276e-05
GMRes iteration 80, residual = 4.4656536251617656e-05
GMRes iteration 81, residual = 4.465441965301862e-05
GMRes iteration 82, residual = 3.8003990530320745e-05
GMRes iteration 83, residual = 3.776120437904957e-05
GMRes iteration 84, residual = 3.305286976432132e-05
GMRes iteration 85, residual = 3.196972018549111e-05
GMRes iteration 86, residual = 2.6396300449021422e-05
GMRes iteration 87, residual = 2.325098853631266e-05
GMRes iteration 88, residual = 2.2059807899664418e-05
GMRes iteration 89, residual = 1.797764113929488e-05
GMRes iteration 90, residual = 1.797632728037468e-05
GMRes iteration 91, residual = 1.4843928478239508e-05
GMRes iteration 92, residual = 1.4315784636439202e-05
GMRes iteration 93, residual = 1.3187362290631918e-05
GMRes iteration 94, residual = 1.1268316443781056e-05
GMRes iteration 95, residual = 1.1238480162838079e-05
GMRes iteration 96, residual = 9.33994439335063e-06
GMRes iteration 97, residual = 9.251556115845772e-06
GMRes iteration 98, residual = 8.234437248431841e-06
GMRes iteration 99, residual = 8.170491739726539e-06
GMRes iteration 100, residual = 6.881178120204716e-06
GMRes iteration 101, residual = 5.953015547417765e-06
GMRes iteration 102, residual = 5.86533907925411e-06
GMRes iteration 103, residual = 4.950597543583038e-06
GMRes iteration 104, residual = 4.937785643897347e-06
GMRes iteration 105, residual = 4.006101326381957e-06
GMRes iteration 106, residual = 3.961999015004813e-06
GMRes iteration 107, residual = 3.5548902447477054e-06
GMRes iteration 108, residual = 3.5138667556152845e-06
GMRes iteration 109, residual = 2.902224323574225e-06
GMRes iteration 110, residual = 2.5822033002841877e-06
GMRes iteration 111, residual = 2.4371752288286588e-06
GMRes iteration 112, residual = 2.099868481832127e-06
GMRes iteration 113, residual = 2.0534867523113927e-06
GMRes iteration 114, residual = 1.8206528427935217e-06
GMRes iteration 115, residual = 1.782074260722328e-06
GMRes iteration 116, residual = 1.4276309258569718e-06
GMRes iteration 117, residual = 1.4275188668549416e-06
GMRes iteration 118, residual = 1.1013563372257773e-06
GMRes iteration 119, residual = 1.0574315945188975e-06
GMRes iteration 120, residual = 9.70722967098118e-07

[8]:

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

The Neumann data:

[9]:

gfv = GridFunction(fesL2)
gfv.vec[:] = sol_sym[1]
Draw (gfv);

References:

[ ]: