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.62741266701414
GMRes iteration 2, residual = 9.928874716575775
GMRes iteration 3, residual = 6.4454282417011814
GMRes iteration 4, residual = 5.914631696346634
GMRes iteration 5, residual = 1.2199739868981212
GMRes iteration 6, residual = 0.7970348672124856
GMRes iteration 7, residual = 0.48377678613842096
GMRes iteration 8, residual = 0.48321174328456074
GMRes iteration 9, residual = 0.25892815233153654
GMRes iteration 10, residual = 0.1971795785062445
GMRes iteration 11, residual = 0.08600238802575388
GMRes iteration 12, residual = 0.08599415993940367
GMRes iteration 13, residual = 0.040332855501947644
GMRes iteration 14, residual = 0.03862202800581738
GMRes iteration 15, residual = 0.021514919943539122
GMRes iteration 16, residual = 0.02145705127878397
GMRes iteration 17, residual = 0.018366954600472846
GMRes iteration 18, residual = 0.018056189886762657
GMRes iteration 19, residual = 0.011039232225375091
GMRes iteration 20, residual = 0.008745646307376278
GMRes iteration 21, residual = 0.008615807362841045
GMRes iteration 22, residual = 0.007377140799106825
GMRes iteration 23, residual = 0.007335584062179396
GMRes iteration 24, residual = 0.005893254134233849
GMRes iteration 25, residual = 0.005865901577006014
GMRes iteration 26, residual = 0.003979408650225095
GMRes iteration 27, residual = 0.003936063788575678
GMRes iteration 28, residual = 0.0035816489154286875
GMRes iteration 29, residual = 0.0035692669307070984
GMRes iteration 30, residual = 0.003271792164223266
GMRes iteration 31, residual = 0.0031658948348597072
GMRes iteration 32, residual = 0.0027463935586110955
GMRes iteration 33, residual = 0.0027284330650108827
GMRes iteration 34, residual = 0.002470624113494183
GMRes iteration 35, residual = 0.0024332604233046934
GMRes iteration 36, residual = 0.0018420002628144168
GMRes iteration 37, residual = 0.001756321463898502
GMRes iteration 38, residual = 0.0016966546910784435
GMRes iteration 39, residual = 0.0015380292334515394
GMRes iteration 40, residual = 0.0015293456552090768
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:

[ ]: