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 = 48.00652698516188
GMRES iteration 2, residual = 10.586564331902602
GMRES iteration 3, residual = 2.612405525139787
GMRES iteration 4, residual = 1.9198941715526345
GMRES iteration 5, residual = 0.41493202317732314
GMRES iteration 6, residual = 0.3902292839584079
GMRES iteration 7, residual = 0.1770295577280713
GMRES iteration 8, residual = 0.1378342707396139
GMRES iteration 9, residual = 0.0884645234132339
GMRES iteration 10, residual = 0.04757992631011475
GMRES iteration 11, residual = 0.04657352674349492
GMRES iteration 12, residual = 0.04487535406464283
GMRES iteration 13, residual = 0.02076476094045796
GMRES iteration 14, residual = 0.013027137422914815
GMRES iteration 15, residual = 0.012421138404594391
GMRES iteration 16, residual = 0.007390298784676543
GMRES iteration 17, residual = 0.007339759964231885
GMRES iteration 18, residual = 0.0061267891105417915
GMRES iteration 19, residual = 0.005585725266425443
GMRES iteration 20, residual = 0.003674537829928322
GMRES iteration 21, residual = 0.0035980449793216137
GMRES iteration 22, residual = 0.0029980985578670705
GMRES iteration 23, residual = 0.002950750367839531
GMRES iteration 24, residual = 0.001980865660431898
GMRES iteration 25, residual = 0.0019423894211251388
GMRES iteration 26, residual = 0.0015599392053133872
GMRES iteration 27, residual = 0.001425749283989699
GMRES iteration 28, residual = 0.001222769637347616
GMRES iteration 29, residual = 0.001086599054837066
GMRES iteration 30, residual = 0.0009854492735543404
GMRES iteration 31, residual = 0.0007496407661299212
GMRES iteration 32, residual = 0.0007494574711856704
GMRES iteration 33, residual = 0.0006133897311235304
GMRES iteration 34, residual = 0.0006123194841762512
GMRES iteration 35, residual = 0.0004931840353189407
GMRES iteration 36, residual = 0.00047807776045590786
GMRES iteration 37, residual = 0.00037069728538418917
GMRES iteration 38, residual = 0.0003700263273660692
GMRES iteration 39, residual = 0.0003126723075786722
GMRES iteration 40, residual = 0.00030934658809024653
GMRES iteration 41, residual = 0.00021775746677276268
GMRES iteration 42, residual = 0.00021366426314632373
GMRES iteration 43, residual = 0.00019135860116463194
GMRES iteration 44, residual = 0.0001776515816624419
GMRES iteration 45, residual = 0.00016128312563950083
GMRES iteration 46, residual = 0.0001382543910923052
GMRES iteration 47, residual = 0.0001376459385938241
GMRES iteration 48, residual = 0.00011138080559744412
GMRES iteration 49, residual = 0.00011127951320764979
GMRES iteration 50, residual = 8.604445391025813e-05
GMRES iteration 51, residual = 8.577737353120347e-05
GMRES iteration 52, residual = 7.217488975576836e-05
GMRES iteration 53, residual = 7.17959966705383e-05
GMRES iteration 54, residual = 5.6586074723698356e-05
GMRES iteration 55, residual = 4.7587272256104905e-05
GMRES iteration 56, residual = 4.6393643945233394e-05
GMRES iteration 57, residual = 3.795535140395034e-05
GMRES iteration 58, residual = 3.7839508681811295e-05
GMRES iteration 59, residual = 2.798586576888185e-05
GMRES iteration 60, residual = 2.754310058581305e-05
GMRES iteration 61, residual = 2.4190486685946766e-05
GMRES iteration 62, residual = 2.4122346216461425e-05
GMRES iteration 63, residual = 2.023690482482504e-05
GMRES iteration 64, residual = 1.9745654385343596e-05
GMRES iteration 65, residual = 1.5929189935082727e-05
GMRES iteration 66, residual = 1.440781736824174e-05
GMRES iteration 67, residual = 1.4098978668006347e-05
GMRES iteration 68, residual = 1.0407643281523522e-05
GMRES iteration 69, residual = 1.032250657081662e-05
GMRES iteration 70, residual = 8.203033601029585e-06
GMRES iteration 71, residual = 7.818669279694059e-06
GMRES iteration 72, residual = 7.03129804310105e-06
GMRES iteration 73, residual = 5.902721900840366e-06
GMRES iteration 74, residual = 5.555659809801821e-06
GMRES iteration 75, residual = 4.1637235132965374e-06
GMRES iteration 76, residual = 4.132576998831749e-06
GMRES iteration 77, residual = 3.6086844635612973e-06
GMRES iteration 78, residual = 3.3306377211991754e-06
GMRES iteration 79, residual = 2.991487291659059e-06
GMRES iteration 80, residual = 2.5567127784714704e-06
GMRES iteration 81, residual = 2.457647116213061e-06
GMRES iteration 82, residual = 2.135231050182665e-06
GMRES iteration 83, residual = 2.080520311873733e-06
GMRES iteration 84, residual = 1.6117645506916891e-06
GMRES iteration 85, residual = 1.5864934135043225e-06
GMRES iteration 86, residual = 1.3479705925963649e-06
GMRES iteration 87, residual = 1.2291379745924514e-06
GMRES iteration 88, residual = 1.1835873681584836e-06
GMRES iteration 89, residual = 9.07034493577789e-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:

[ ]: