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/docu-ngsbem/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.94329927058916
GMRes iteration 2, residual = 10.547738043198812
GMRes iteration 3, residual = 3.4483285857125687
GMRes iteration 4, residual = 3.3925820325576543
GMRes iteration 5, residual = 0.5788213577921321
GMRes iteration 6, residual = 0.5480288625880633
GMRes iteration 7, residual = 0.30725442713268775
GMRes iteration 8, residual = 0.28418746451106636
GMRes iteration 9, residual = 0.1824812078399675
GMRes iteration 10, residual = 0.15076712396700756
GMRes iteration 11, residual = 0.1035901619331933
GMRes iteration 12, residual = 0.06769258395135722
GMRes iteration 13, residual = 0.04957676436234772
GMRes iteration 14, residual = 0.04077663433357896
GMRes iteration 15, residual = 0.029549247788983362
GMRes iteration 16, residual = 0.021772314576842842
GMRes iteration 17, residual = 0.00961715794731315
GMRes iteration 18, residual = 0.009543689762620112
GMRes iteration 19, residual = 0.007177978091955789
GMRes iteration 20, residual = 0.006301263592913071
GMRes iteration 21, residual = 0.005055306825730399
GMRes iteration 22, residual = 0.005054997659405165
GMRes iteration 23, residual = 0.0036880715288590583
GMRes iteration 24, residual = 0.003676864426682677
GMRes iteration 25, residual = 0.002741143442978121
GMRes iteration 26, residual = 0.002584536847490454
GMRes iteration 27, residual = 0.002120620977142287
GMRes iteration 28, residual = 0.0018677183248110327
GMRes iteration 29, residual = 0.0016414095576520044
GMRes iteration 30, residual = 0.0015190051603135007
GMRes iteration 31, residual = 0.0013581605813299841
GMRes iteration 32, residual = 0.0011735638395737832
GMRes iteration 33, residual = 0.0010714593948860393
GMRes iteration 34, residual = 0.000887863259104315
GMRes iteration 35, residual = 0.0008809312918839086
GMRes iteration 36, residual = 0.0007096729192859276
GMRes iteration 37, residual = 0.0007090872559797519
GMRes iteration 38, residual = 0.0005151523371321297
GMRes iteration 39, residual = 0.0005068078610138329
GMRes iteration 40, residual = 0.00042209685707496905
GMRes iteration 41, residual = 0.0004085284067215565
GMRes iteration 42, residual = 0.00034712205995494406
GMRes iteration 43, residual = 0.00030989617765736694
GMRes iteration 44, residual = 0.0002831358574481815
GMRes iteration 45, residual = 0.00026616094047990114
GMRes iteration 46, residual = 0.0002536862517620302
GMRes iteration 47, residual = 0.00022354523491784645
GMRes iteration 48, residual = 0.00021509103762914073
GMRes iteration 49, residual = 0.00018305148965133762
GMRes iteration 50, residual = 0.00018222568212668776
GMRes iteration 51, residual = 0.00015588341615155636
GMRes iteration 52, residual = 0.00015575325517926944
GMRes iteration 53, residual = 0.00012174293221802565
GMRes iteration 54, residual = 0.00011889739973499524
GMRes iteration 55, residual = 0.00010362884554101464
GMRes iteration 56, residual = 9.575242256303911e-05
GMRes iteration 57, residual = 8.649895735144994e-05
GMRes iteration 58, residual = 7.242181433823682e-05
GMRes iteration 59, residual = 7.157450665931478e-05
GMRes iteration 60, residual = 5.90449827073974e-05
GMRes iteration 61, residual = 5.9039211548073216e-05
GMRes iteration 62, residual = 4.5735890742383344e-05
GMRes iteration 63, residual = 4.379661521678233e-05
GMRes iteration 64, residual = 3.9167922409584904e-05
GMRes iteration 65, residual = 3.384628449411992e-05
GMRes iteration 66, residual = 3.2776905794666235e-05
GMRes iteration 67, residual = 2.6732864319467625e-05
GMRes iteration 68, residual = 2.6715947512677667e-05
GMRes iteration 69, residual = 2.3007595156296237e-05
GMRes iteration 70, residual = 2.240834057533737e-05
GMRes iteration 71, residual = 1.9611600590249218e-05
GMRes iteration 72, residual = 1.8027032908821283e-05
GMRes iteration 73, residual = 1.7047810675435233e-05
GMRes iteration 74, residual = 1.300136152907334e-05
GMRes iteration 75, residual = 1.3001069941380829e-05
GMRes iteration 76, residual = 1.1003843614700916e-05
GMRes iteration 77, residual = 1.0766569571240568e-05
GMRes iteration 78, residual = 9.105788260342922e-06
GMRes iteration 79, residual = 8.390870659815246e-06
GMRes iteration 80, residual = 7.625740727052213e-06
GMRes iteration 81, residual = 6.677095448192493e-06
GMRes iteration 82, residual = 6.589380454426059e-06
GMRes iteration 83, residual = 5.1437674752334716e-06
GMRes iteration 84, residual = 5.035758369183206e-06
GMRes iteration 85, residual = 4.278266760698381e-06
GMRes iteration 86, residual = 4.0151177110996215e-06
GMRes iteration 87, residual = 3.8039859813853494e-06
GMRes iteration 88, residual = 3.240823744376326e-06
GMRes iteration 89, residual = 3.1863873797667784e-06
GMRes iteration 90, residual = 2.4850582267890464e-06
GMRes iteration 91, residual = 2.4128179762504803e-06
GMRes iteration 92, residual = 2.2270510940001902e-06
GMRes iteration 93, residual = 1.9475212888090293e-06
GMRes iteration 94, residual = 1.8639868506496959e-06
GMRes iteration 95, residual = 1.5663649227706934e-06
GMRes iteration 96, residual = 1.5310935683152356e-06
GMRes iteration 97, residual = 1.2919090812007111e-06
GMRes iteration 98, residual = 1.2912904951295219e-06
GMRes iteration 99, residual = 9.998111603537845e-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:

[ ]: