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 = 3
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 =  39139

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 =  2691

Generate the the single layer potential \(V\), double layer potential \(K\) and hypersingular operator \(D\):

[5]:

with TaskManager():
    V = SingleLayerPotentialOperator(fesL2, intorder=7)
    K = DoubleLayerPotentialOperator(fesH1, fesL2, trial_definedon=mesh.Boundaries("outer"),  test_definedon=mesh.Boundaries("outer"), intorder=7)
    D = HypersingularOperator(fesH1, definedon=mesh.Boundaries("outer"), intorder=7)
    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=40, printrates=True)

GMRes iteration 1, residual = 41.627412667014106
GMRes iteration 2, residual = 9.928873774991159
GMRes iteration 3, residual = 6.445427133153754
GMRes iteration 4, residual = 5.914629904071614
GMRes iteration 5, residual = 1.2199751335720221
GMRes iteration 6, residual = 0.797034594354077
GMRes iteration 7, residual = 0.4837866120144668
GMRes iteration 8, residual = 0.48322004881077285
GMRes iteration 9, residual = 0.2589276665593373
GMRes iteration 10, residual = 0.1971800820002907
GMRes iteration 11, residual = 0.08600297806784656
GMRes iteration 12, residual = 0.08599473481796722
GMRes iteration 13, residual = 0.040373843130155664
GMRes iteration 14, residual = 0.03863958814958185
GMRes iteration 15, residual = 0.021515510054387652
GMRes iteration 16, residual = 0.021457586214434072
GMRes iteration 17, residual = 0.01836923771255684
GMRes iteration 18, residual = 0.01805907256483053
GMRes iteration 19, residual = 0.011040139299988138
GMRes iteration 20, residual = 0.008745800929747747
GMRes iteration 21, residual = 0.008616040866316153
GMRes iteration 22, residual = 0.00737662796069357
GMRes iteration 23, residual = 0.007335467214129255
GMRes iteration 24, residual = 0.005893377495134481
GMRes iteration 25, residual = 0.0058661296334348795
GMRes iteration 26, residual = 0.003979099613660768
GMRes iteration 27, residual = 0.003935342958145997
GMRes iteration 28, residual = 0.0035819844055376576
GMRes iteration 29, residual = 0.003569513336137529
GMRes iteration 30, residual = 0.003272112326119759
GMRes iteration 31, residual = 0.003166009297680206
GMRes iteration 32, residual = 0.002745730614940562
GMRes iteration 33, residual = 0.002728140005140794
GMRes iteration 34, residual = 0.0024702926629182626
GMRes iteration 35, residual = 0.0024331340575076887
GMRes iteration 36, residual = 0.0018419652444669422
GMRes iteration 37, residual = 0.0017560242826115862
GMRes iteration 38, residual = 0.0016967484120416526
GMRes iteration 39, residual = 0.0015380509230964482
GMRes iteration 40, residual = 0.0015293962561646304
WARNING: GMRes did not converge to TOL

[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:

[ ]: