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.943299270587374
GMRes iteration 2, residual = 11.37309019139713
GMRes iteration 3, residual = 6.709311492493946
GMRes iteration 4, residual = 6.495804108604595
GMRes iteration 5, residual = 1.432888166609496
GMRes iteration 6, residual = 0.8253014999039172
GMRes iteration 7, residual = 0.6024935444751981
GMRes iteration 8, residual = 0.3865476425730732
GMRes iteration 9, residual = 0.2753852464823924
GMRes iteration 10, residual = 0.25092009338554716
GMRes iteration 11, residual = 0.100036518136328
GMRes iteration 12, residual = 0.09550113470008298
GMRes iteration 13, residual = 0.04257386604916562
GMRes iteration 14, residual = 0.030542749705618356
GMRes iteration 15, residual = 0.02393386379211333
GMRes iteration 16, residual = 0.023933859774561094
GMRes iteration 17, residual = 0.018760693725179054
GMRes iteration 18, residual = 0.016374848841446285
GMRes iteration 19, residual = 0.014108850688120265
GMRes iteration 20, residual = 0.00992466092863145
GMRes iteration 21, residual = 0.009920782010712973
GMRes iteration 22, residual = 0.007811734811764174
GMRes iteration 23, residual = 0.007766703224033538
GMRes iteration 24, residual = 0.006225361576606542
GMRes iteration 25, residual = 0.006203569039529069
GMRes iteration 26, residual = 0.003999570702862017
GMRes iteration 27, residual = 0.003884270768720016
GMRes iteration 28, residual = 0.0035496185564396344
GMRes iteration 29, residual = 0.0035459215913820545
GMRes iteration 30, residual = 0.003085763331684475
GMRes iteration 31, residual = 0.0030514890186567234
GMRes iteration 32, residual = 0.002711644565026214
GMRes iteration 33, residual = 0.0026991899645646363
GMRes iteration 34, residual = 0.0023762747339659803
GMRes iteration 35, residual = 0.0021121224803445017
GMRes iteration 36, residual = 0.002007320868147828
GMRes iteration 37, residual = 0.0018686299551601618
GMRes iteration 38, residual = 0.0018487004698324388
GMRes iteration 39, residual = 0.001604345815215127
GMRes iteration 40, residual = 0.0015902550149328753
GMRes iteration 41, residual = 0.0014588203814043979
GMRes iteration 42, residual = 0.0014220388253412036
GMRes iteration 43, residual = 0.001261297830581665
GMRes iteration 44, residual = 0.0011942142537612825
GMRes iteration 45, residual = 0.0010979317892685042
GMRes iteration 46, residual = 0.0010972029931398198
GMRes iteration 47, residual = 0.0009820600128631657
GMRes iteration 48, residual = 0.0008652069575475605
GMRes iteration 49, residual = 0.0008166912490914975
GMRes iteration 50, residual = 0.0007499606390021204
GMRes iteration 51, residual = 0.0007493421050183485
GMRes iteration 52, residual = 0.0006249378469087029
GMRes iteration 53, residual = 0.0005734273618410702
GMRes iteration 54, residual = 0.0005321357108820569
GMRes iteration 55, residual = 0.0005132383404753959
GMRes iteration 56, residual = 0.0004595144340187492
GMRes iteration 57, residual = 0.00034788839256089403
GMRes iteration 58, residual = 0.0003431734893607763
GMRes iteration 59, residual = 0.0003074581428660961
GMRes iteration 60, residual = 0.00029950452193204313
GMRes iteration 61, residual = 0.00023276693956946916
GMRes iteration 62, residual = 0.00022878483165270132
GMRes iteration 63, residual = 0.00020593185931006632
GMRes iteration 64, residual = 0.00019438961560342598
GMRes iteration 65, residual = 0.00018032111217709668
GMRes iteration 66, residual = 0.0001537051705057386
GMRes iteration 67, residual = 0.00015320431617917035
GMRes iteration 68, residual = 0.0001317013908767591
GMRes iteration 69, residual = 0.00013039965122131516
GMRes iteration 70, residual = 0.00010479060511074911
GMRes iteration 71, residual = 9.437471210143418e-05
GMRes iteration 72, residual = 8.93874142320095e-05
GMRes iteration 73, residual = 7.940780243608517e-05
GMRes iteration 74, residual = 7.849798235290602e-05
GMRes iteration 75, residual = 6.6862328396764e-05
GMRes iteration 76, residual = 6.685682214559569e-05
GMRes iteration 77, residual = 5.533149625429129e-05
GMRes iteration 78, residual = 5.466998824303849e-05
GMRes iteration 79, residual = 5.150509066531941e-05
GMRes iteration 80, residual = 4.465656762469078e-05
GMRes iteration 81, residual = 4.4654452271914963e-05
GMRes iteration 82, residual = 3.80041105501308e-05
GMRes iteration 83, residual = 3.776130024162954e-05
GMRes iteration 84, residual = 3.305288766129751e-05
GMRes iteration 85, residual = 3.196971598478737e-05
GMRes iteration 86, residual = 2.6396299143075812e-05
GMRes iteration 87, residual = 2.3250898374692935e-05
GMRes iteration 88, residual = 2.2059411209165145e-05
GMRes iteration 89, residual = 1.797695656250345e-05
GMRes iteration 90, residual = 1.7975588993564153e-05
GMRes iteration 91, residual = 1.4840860220916414e-05
GMRes iteration 92, residual = 1.4314584442481067e-05
GMRes iteration 93, residual = 1.318092897187391e-05
GMRes iteration 94, residual = 1.1266687634340902e-05
GMRes iteration 95, residual = 1.1235890036178718e-05
GMRes iteration 96, residual = 9.338132710098977e-06
GMRes iteration 97, residual = 9.250845678955001e-06
GMRes iteration 98, residual = 8.233660502698412e-06
GMRes iteration 99, residual = 8.170181843316822e-06
GMRes iteration 100, residual = 6.880072878701697e-06
GMRes iteration 101, residual = 5.952924253691693e-06
GMRes iteration 102, residual = 5.865156779281997e-06
GMRes iteration 103, residual = 4.950606062857788e-06
GMRes iteration 104, residual = 4.937794613939789e-06
GMRes iteration 105, residual = 4.006097527725753e-06
GMRes iteration 106, residual = 3.96199898688428e-06
GMRes iteration 107, residual = 3.554876321657848e-06
GMRes iteration 108, residual = 3.51385408839433e-06
GMRes iteration 109, residual = 2.9022258758196686e-06
GMRes iteration 110, residual = 2.582204940006798e-06
GMRes iteration 111, residual = 2.4371797133424416e-06
GMRes iteration 112, residual = 2.099863494022419e-06
GMRes iteration 113, residual = 2.053482822254799e-06
GMRes iteration 114, residual = 1.8206528826649568e-06
GMRes iteration 115, residual = 1.7820743005014386e-06
GMRes iteration 116, residual = 1.4276317043152141e-06
GMRes iteration 117, residual = 1.4275196676966592e-06
GMRes iteration 118, residual = 1.1013570261973385e-06
GMRes iteration 119, residual = 1.0574313231595184e-06
GMRes iteration 120, residual = 9.707233199542e-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:

[ ]: