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.94329927058771
GMRes iteration 2, residual = 11.373090191396354
GMRes iteration 3, residual = 6.709311492493795
GMRes iteration 4, residual = 6.495804108604396
GMRes iteration 5, residual = 1.4328881666094997
GMRes iteration 6, residual = 0.8253014999039162
GMRes iteration 7, residual = 0.6024935444752104
GMRes iteration 8, residual = 0.3865476425803826
GMRes iteration 9, residual = 0.2753852464821465
GMRes iteration 10, residual = 0.25092009338512716
GMRes iteration 11, residual = 0.10003651813621979
GMRes iteration 12, residual = 0.09550113470005434
GMRes iteration 13, residual = 0.04257386604903423
GMRes iteration 14, residual = 0.0305427497054748
GMRes iteration 15, residual = 0.02393386379216942
GMRes iteration 16, residual = 0.023933859774617032
GMRes iteration 17, residual = 0.018760693725363833
GMRes iteration 18, residual = 0.016374848841618876
GMRes iteration 19, residual = 0.01410885068807812
GMRes iteration 20, residual = 0.0099246609286574
GMRes iteration 21, residual = 0.009920782010736753
GMRes iteration 22, residual = 0.007811734811808763
GMRes iteration 23, residual = 0.007766703224089924
GMRes iteration 24, residual = 0.006225361576649208
GMRes iteration 25, residual = 0.006203569039581713
GMRes iteration 26, residual = 0.00399957070284462
GMRes iteration 27, residual = 0.0038842707687045527
GMRes iteration 28, residual = 0.003549618556317713
GMRes iteration 29, residual = 0.0035459215912725535
GMRes iteration 30, residual = 0.003085763331928212
GMRes iteration 31, residual = 0.0030514890188998007
GMRes iteration 32, residual = 0.0027116445649832577
GMRes iteration 33, residual = 0.0026991899645193254
GMRes iteration 34, residual = 0.0023762747340604455
GMRes iteration 35, residual = 0.002112122480390149
GMRes iteration 36, residual = 0.0020073208683879405
GMRes iteration 37, residual = 0.0018686299552511012
GMRes iteration 38, residual = 0.0018487004700786958
GMRes iteration 39, residual = 0.0016043458156370662
GMRes iteration 40, residual = 0.0015902550152434216
GMRes iteration 41, residual = 0.0014588203814829145
GMRes iteration 42, residual = 0.0014220388253017714
GMRes iteration 43, residual = 0.0012612978295011693
GMRes iteration 44, residual = 0.001194214252807001
GMRes iteration 45, residual = 0.0010979317891140867
GMRes iteration 46, residual = 0.0010972029930724997
GMRes iteration 47, residual = 0.0009820600119964907
GMRes iteration 48, residual = 0.0008652069571552921
GMRes iteration 49, residual = 0.0008166912476058391
GMRes iteration 50, residual = 0.0007499606377514711
GMRes iteration 51, residual = 0.0007493421035949162
GMRes iteration 52, residual = 0.0006249378456919045
GMRes iteration 53, residual = 0.0005734273616001876
GMRes iteration 54, residual = 0.0005321357100808667
GMRes iteration 55, residual = 0.000513238339893283
GMRes iteration 56, residual = 0.00045951443258111443
GMRes iteration 57, residual = 0.0003478883946373303
GMRes iteration 58, residual = 0.00034317349191061814
GMRes iteration 59, residual = 0.00030745814309624566
GMRes iteration 60, residual = 0.0002995045223094984
GMRes iteration 61, residual = 0.00023276694185413575
GMRes iteration 62, residual = 0.00022878483326488805
GMRes iteration 63, residual = 0.00020593186817601447
GMRes iteration 64, residual = 0.00019438961786531628
GMRes iteration 65, residual = 0.00018032112485665436
GMRes iteration 66, residual = 0.00015370518081481155
GMRes iteration 67, residual = 0.00015320432926314403
GMRes iteration 68, residual = 0.0001317013940294286
GMRes iteration 69, residual = 0.00013039965197685625
GMRes iteration 70, residual = 0.00010479061027489897
GMRes iteration 71, residual = 9.43747128836951e-05
GMRes iteration 72, residual = 8.93874027612684e-05
GMRes iteration 73, residual = 7.940779943533882e-05
GMRes iteration 74, residual = 7.849797729740186e-05
GMRes iteration 75, residual = 6.686231676812335e-05
GMRes iteration 76, residual = 6.685681014989007e-05
GMRes iteration 77, residual = 5.5331475560636595e-05
GMRes iteration 78, residual = 5.466997404202136e-05
GMRes iteration 79, residual = 5.150510946899554e-05
GMRes iteration 80, residual = 4.465655545163615e-05
GMRes iteration 81, residual = 4.465444291213635e-05
GMRes iteration 82, residual = 3.8004161896481976e-05
GMRes iteration 83, residual = 3.7761290647536305e-05
GMRes iteration 84, residual = 3.305306016592912e-05
GMRes iteration 85, residual = 3.1969715618314015e-05
GMRes iteration 86, residual = 2.6397369206323075e-05
GMRes iteration 87, residual = 2.3251174239256122e-05
GMRes iteration 88, residual = 2.20618401175677e-05
GMRes iteration 89, residual = 1.798025033917766e-05
GMRes iteration 90, residual = 1.7979111669387345e-05
GMRes iteration 91, residual = 1.4856923347400344e-05
GMRes iteration 92, residual = 1.4320284796237752e-05
GMRes iteration 93, residual = 1.3213917200747158e-05
GMRes iteration 94, residual = 1.1275605550096797e-05
GMRes iteration 95, residual = 1.1250190603968575e-05
GMRes iteration 96, residual = 9.348302181632187e-06
GMRes iteration 97, residual = 9.254620912494051e-06
GMRes iteration 98, residual = 8.239317118735207e-06
GMRes iteration 99, residual = 8.172715738322892e-06
GMRes iteration 100, residual = 6.886199115523116e-06
GMRes iteration 101, residual = 5.95347391544934e-06
GMRes iteration 102, residual = 5.8662640935286775e-06
GMRes iteration 103, residual = 4.950720215546106e-06
GMRes iteration 104, residual = 4.937968408893128e-06
GMRes iteration 105, residual = 4.006219158465873e-06
GMRes iteration 106, residual = 3.962064181378853e-06
GMRes iteration 107, residual = 3.554925686818733e-06
GMRes iteration 108, residual = 3.5138892284124515e-06
GMRes iteration 109, residual = 2.90223510259755e-06
GMRes iteration 110, residual = 2.582212282205764e-06
GMRes iteration 111, residual = 2.437187597952471e-06
GMRes iteration 112, residual = 2.0998742591514735e-06
GMRes iteration 113, residual = 2.053495769988097e-06
GMRes iteration 114, residual = 1.8206540381292895e-06
GMRes iteration 115, residual = 1.7820761322652612e-06
GMRes iteration 116, residual = 1.4276323763109805e-06
GMRes iteration 117, residual = 1.4275203679806913e-06
GMRes iteration 118, residual = 1.1013583774237323e-06
GMRes iteration 119, residual = 1.0574331882077212e-06
GMRes iteration 120, residual = 9.70723653983673e-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:

[ ]: