This page was generated from unit-2.14-globalspaces/planewavecoupling.ipynb.

2.14 Coupling of FEM with plane waves

In this tutorial we define global basis functions via Python. As an example, we choose out-going plane waves in a wave-guide.

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw

r1 = Rectangle(1,1).Face()
r1.edges.Min(Y).name = "dir"
r1.edges.Max(Y).name = "dir"
r1.edges.Min(X).name = "dir"
r1.edges.Max(X).name = "coupling"
r1.faces.name = "fem"
circ = Circle((0.4,0.7),0.02).Face()
circ.faces.name = "femsource"

circ2 = Circle((0.4,0.3),0.05).Face()
circ2.faces.name = "femobstacle"


r2 = MoveTo(1,0).Rectangle(3,1).Face()
r2.edges.Min(Y).name = "dir"
r2.edges.Max(Y).name = "dir"
r2.faces.name = "waves"
shape = Glue([r1-circ-circ2, circ,circ2,r2])

The left side is discretized by finite elements. On the right side we use (propagating and evanescent) plane waves. On the interface both fields are coupled by a ultra-weak variational formulation.

Draw (shape);

geo = OCCGeometry(shape, dim=2)
mesh = Mesh(geo.GenerateMesh(maxh=0.2)).Curve(3)
Draw (mesh);

outgoing plane waves:

\[\sin (k_{y,j} y) e^{i k_{x,j} x}\]

with \begin{eqnarray*} k_{y,j} & = & j \pi \\ k_{x,j} & = & \sqrt{\omega^2 - k_{y,j}^2} \end{eqnarray*}

Define them as multidimensional CoefficientFunction:

omega = 3.5*pi * (1+0j)

ky = [j*pi for j in range(1,5)]
kx = [sqrt(omega**2-kyi**2) for kyi in ky]
k = zip(ky,kx)
shapes = CF(tuple(sin(kyi*y)*exp(1j*kxi*x) for kyi,kxi in k))
dshapesx = shapes.Diff(x)

Draw (shapes[2], mesh, animate_complex=True );

from ngsolve.comp import GlobalSpace

fesfem = H1(mesh, order=4, definedon="fem.*", dirichlet="dir", complex=True)
profile=sin(pi*y)
feswaves = GlobalSpace (mesh, definedon="waves", basis = shapes)
feswaves.AddOperator("dn", BND, dshapesx)
fes = fesfem*feswaves

fes.ndof, fesfem.ndof, feswaves.ndof

(3685, 3681, 4)

Bilinear-form:

\[\int_{\Omega_{FEM}} \nabla u \nabla v - \omega^2 \varepsilon u v dx + \int_{\gamma_{coupling}} - u \partial_n v_w - v \partial_n u_w + u_w \partial_n v_w\]

u,v = fes.TnT()
(uf,uw), (vf,vw) = fes.TnT()
uwdx = uw.Operator("dn")
vwdx = vw.Operator("dn")

epsilon = mesh.MaterialCF( { "obstacle" : 10 }, 1)

a = BilinearForm(fes)
a += (grad(uf))*grad(vf)*dx("fem.*") - omega**2*epsilon*uf*vf*dx("fem.*")
a += (-uf*vwdx-vf*uwdx+uw.Trace()*vwdx)*ds("coupling")
a.Assemble()

f = LinearForm(1*vf*dx("femsource")).Assemble()

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec

gfu2 = mesh.MaterialCF( { "fem.*" : gfu.components[0], "waves" : gfu.components[1] })
Draw (gfu2, mesh, animate_complex=True );