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.94329927058598
GMRES iteration 2, residual = 10.547738042890783
GMRES iteration 3, residual = 2.6122774443199845
GMRES iteration 4, residual = 1.9193598537426715
GMRES iteration 5, residual = 0.4151312135075365
GMRES iteration 6, residual = 0.3904448743752739
GMRES iteration 7, residual = 0.17719618613351396
GMRES iteration 8, residual = 0.13791749556355007
GMRES iteration 9, residual = 0.08865763341000958
GMRES iteration 10, residual = 0.04770629758401662
GMRES iteration 11, residual = 0.04662109067768837
GMRES iteration 12, residual = 0.044970971312833015
GMRES iteration 13, residual = 0.020803869025729033
GMRES iteration 14, residual = 0.013029445657634564
GMRES iteration 15, residual = 0.012428849092442934
GMRES iteration 16, residual = 0.007375363713895484
GMRES iteration 17, residual = 0.0073286714915262305
GMRES iteration 18, residual = 0.006116844119486991
GMRES iteration 19, residual = 0.0055747504431094596
GMRES iteration 20, residual = 0.003669977355663017
GMRES iteration 21, residual = 0.0035919757619537252
GMRES iteration 22, residual = 0.0029943239014043744
GMRES iteration 23, residual = 0.002946929400715947
GMRES iteration 24, residual = 0.001979267927169113
GMRES iteration 25, residual = 0.001940866488549286
GMRES iteration 26, residual = 0.0015575489596132738
GMRES iteration 27, residual = 0.0014246204450290769
GMRES iteration 28, residual = 0.001220442512079464
GMRES iteration 29, residual = 0.0010847044561632833
GMRES iteration 30, residual = 0.0009837865749438162
GMRES iteration 31, residual = 0.0007490400520314075
GMRES iteration 32, residual = 0.0007488559796138245
GMRES iteration 33, residual = 0.0006134924989505811
GMRES iteration 34, residual = 0.0006124109584833781
GMRES iteration 35, residual = 0.0004939642690218078
GMRES iteration 36, residual = 0.00047858686160329075
GMRES iteration 37, residual = 0.00037174187502318176
GMRES iteration 38, residual = 0.0003709704896931077
GMRES iteration 39, residual = 0.00031421102531648334
GMRES iteration 40, residual = 0.00031040251753934004
GMRES iteration 41, residual = 0.00021837357945550724
GMRES iteration 42, residual = 0.00021449524405554044
GMRES iteration 43, residual = 0.00019217128046502004
GMRES iteration 44, residual = 0.00017880247628143133
GMRES iteration 45, residual = 0.00016236092092361998
GMRES iteration 46, residual = 0.0001393585424121001
GMRES iteration 47, residual = 0.00013877785900205094
GMRES iteration 48, residual = 0.00011238293927326688
GMRES iteration 49, residual = 0.0001122642292852316
GMRES iteration 50, residual = 8.692521861922205e-05
GMRES iteration 51, residual = 8.664928872818024e-05
GMRES iteration 52, residual = 7.291903344109864e-05
GMRES iteration 53, residual = 7.253473755928916e-05
GMRES iteration 54, residual = 5.704763944971428e-05
GMRES iteration 55, residual = 4.7964544398149326e-05
GMRES iteration 56, residual = 4.6725280152449115e-05
GMRES iteration 57, residual = 3.8199258451370194e-05
GMRES iteration 58, residual = 3.807476463584527e-05
GMRES iteration 59, residual = 2.8152110655246814e-05
GMRES iteration 60, residual = 2.7724815845884198e-05
GMRES iteration 61, residual = 2.4337832357387283e-05
GMRES iteration 62, residual = 2.4272918487358856e-05
GMRES iteration 63, residual = 2.036836917492015e-05
GMRES iteration 64, residual = 1.987315615723596e-05
GMRES iteration 65, residual = 1.6038071691296974e-05
GMRES iteration 66, residual = 1.447703672022613e-05
GMRES iteration 67, residual = 1.4163716842714894e-05
GMRES iteration 68, residual = 1.0409809064802316e-05
GMRES iteration 69, residual = 1.0328804506665193e-05
GMRES iteration 70, residual = 8.190112716821977e-06
GMRES iteration 71, residual = 7.817834143420838e-06
GMRES iteration 72, residual = 7.018532976374305e-06
GMRES iteration 73, residual = 5.897819858177622e-06
GMRES iteration 74, residual = 5.545471214116591e-06
GMRES iteration 75, residual = 4.159607901710078e-06
GMRES iteration 76, residual = 4.128810904976767e-06
GMRES iteration 77, residual = 3.6120160262373897e-06
GMRES iteration 78, residual = 3.35157784155421e-06
GMRES iteration 79, residual = 2.9846775091539334e-06
GMRES iteration 80, residual = 2.5627033721733283e-06
GMRES iteration 81, residual = 2.4522367055776266e-06
GMRES iteration 82, residual = 2.1344262381718388e-06
GMRES iteration 83, residual = 2.0774974697898625e-06
GMRES iteration 84, residual = 1.611673536960832e-06
GMRES iteration 85, residual = 1.5868596963214546e-06
GMRES iteration 86, residual = 1.3487736097536368e-06
GMRES iteration 87, residual = 1.2296439913259345e-06
GMRES iteration 88, residual = 1.1847030105778293e-06
GMRES iteration 89, residual = 9.075403848051566e-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:

[ ]: