This page was generated from wta/fembem.ipynb.

FEM-BEM Coupling

The ngbem boundary element addon project initiated by Lucy Weggler (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 representation 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 = 53.25532390083373
GMRES iteration 2, residual = 11.204020448047672
GMRES iteration 3, residual = 2.624791062567705
GMRES iteration 4, residual = 1.9538056659157523
GMRES iteration 5, residual = 0.4224026151805053
GMRES iteration 6, residual = 0.39752546220394414
GMRES iteration 7, residual = 0.18178805778523824
GMRES iteration 8, residual = 0.14276786045224654
GMRES iteration 9, residual = 0.09379134926514958
GMRES iteration 10, residual = 0.049119570807404055
GMRES iteration 11, residual = 0.04817369825668193
GMRES iteration 12, residual = 0.046262474314091945
GMRES iteration 13, residual = 0.022592068478819193
GMRES iteration 14, residual = 0.013683108949373012
GMRES iteration 15, residual = 0.013135234018673995
GMRES iteration 16, residual = 0.007865742072748623
GMRES iteration 17, residual = 0.00784509683343906
GMRES iteration 18, residual = 0.006510891537841356
GMRES iteration 19, residual = 0.006010697481694803
GMRES iteration 20, residual = 0.003998779340045902
GMRES iteration 21, residual = 0.00394271137501964
GMRES iteration 22, residual = 0.0032959576927091452
GMRES iteration 23, residual = 0.0032676763957131733
GMRES iteration 24, residual = 0.0022444259552192093
GMRES iteration 25, residual = 0.0021980519789551126
GMRES iteration 26, residual = 0.0017818624071490915
GMRES iteration 27, residual = 0.001620143396078669
GMRES iteration 28, residual = 0.0014081996060172654
GMRES iteration 29, residual = 0.0012630114127908916
GMRES iteration 30, residual = 0.0011495135581383
GMRES iteration 31, residual = 0.0008726235776975727
GMRES iteration 32, residual = 0.0008724249548266068
GMRES iteration 33, residual = 0.000707549203856537
GMRES iteration 34, residual = 0.0007065256887381155
GMRES iteration 35, residual = 0.0005653262678264924
GMRES iteration 36, residual = 0.0005492501720333716
GMRES iteration 37, residual = 0.0004279646754533122
GMRES iteration 38, residual = 0.00042713841312304297
GMRES iteration 39, residual = 0.0003573354498777612
GMRES iteration 40, residual = 0.0003539012111881696
GMRES iteration 41, residual = 0.0002495237364927417
GMRES iteration 42, residual = 0.0002441042759853239
GMRES iteration 43, residual = 0.00021991310535018048
GMRES iteration 44, residual = 0.00020009350267417077
GMRES iteration 45, residual = 0.00018200103240535308
GMRES iteration 46, residual = 0.0001549541512790952
GMRES iteration 47, residual = 0.00015424369712231634
GMRES iteration 48, residual = 0.00012348601905165028
GMRES iteration 49, residual = 0.00012315028815622502
GMRES iteration 50, residual = 9.585974626740955e-05
GMRES iteration 51, residual = 9.550860400240718e-05
GMRES iteration 52, residual = 8.056738241764894e-05
GMRES iteration 53, residual = 8.003746981091063e-05
GMRES iteration 54, residual = 6.259698615535514e-05
GMRES iteration 55, residual = 5.3034259643549934e-05
GMRES iteration 56, residual = 5.179762066346828e-05
GMRES iteration 57, residual = 4.215166724041719e-05
GMRES iteration 58, residual = 4.206263910145165e-05
GMRES iteration 59, residual = 3.1136411749257e-05
GMRES iteration 60, residual = 3.0632022719325635e-05
GMRES iteration 61, residual = 2.6975329564522053e-05
GMRES iteration 62, residual = 2.6891259832926744e-05
GMRES iteration 63, residual = 2.251191808825054e-05
GMRES iteration 64, residual = 2.1856284813043015e-05
GMRES iteration 65, residual = 1.7843265227585427e-05
GMRES iteration 66, residual = 1.5987665802973065e-05
GMRES iteration 67, residual = 1.569849494733274e-05
GMRES iteration 68, residual = 1.1403072073395672e-05
GMRES iteration 69, residual = 1.129215033507851e-05
GMRES iteration 70, residual = 9.159660935551851e-06
GMRES iteration 71, residual = 8.673524519853627e-06
GMRES iteration 72, residual = 7.84727082910624e-06
GMRES iteration 73, residual = 6.4503854242587836e-06
GMRES iteration 74, residual = 6.112145774037074e-06
GMRES iteration 75, residual = 4.612589729819061e-06
GMRES iteration 76, residual = 4.567463357689726e-06
GMRES iteration 77, residual = 3.946178346343645e-06
GMRES iteration 78, residual = 3.616699239839648e-06
GMRES iteration 79, residual = 3.2577032235409573e-06
GMRES iteration 80, residual = 2.7764524469472387e-06
GMRES iteration 81, residual = 2.685457667467237e-06
GMRES iteration 82, residual = 2.3249636563810243e-06
GMRES iteration 83, residual = 2.265725891313645e-06
GMRES iteration 84, residual = 1.7588461609564698e-06
GMRES iteration 85, residual = 1.7201161106809023e-06
GMRES iteration 86, residual = 1.469101830421169e-06
GMRES iteration 87, residual = 1.3378421140976064e-06
GMRES iteration 88, residual = 1.2844909461142551e-06
GMRES iteration 89, residual = 9.86517971895477e-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: