{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# 2.14 Coupling of FEM with plane waves\n",
    "\n",
    "In this tutorial we define global basis functions via Python. As an example, we choose out-going plane waves in a wave-guide."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from netgen.occ import *\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "r1 = Rectangle(1,1).Face()\n",
    "r1.edges.Min(Y).name = \"dir\"\n",
    "r1.edges.Max(Y).name = \"dir\"\n",
    "r1.edges.Min(X).name = \"dir\"\n",
    "r1.edges.Max(X).name = \"coupling\"\n",
    "r1.faces.name = \"fem\"\n",
    "circ = Circle((0.4,0.7),0.02).Face()\n",
    "circ.faces.name = \"femsource\"\n",
    "\n",
    "circ2 = Circle((0.4,0.3),0.05).Face()\n",
    "circ2.faces.name = \"femobstacle\"\n",
    "\n",
    "\n",
    "r2 = MoveTo(1,0).Rectangle(3,1).Face()\n",
    "r2.edges.Min(Y).name = \"dir\"\n",
    "r2.edges.Max(Y).name = \"dir\"\n",
    "r2.faces.name = \"waves\"\n",
    "shape = Glue([r1-circ-circ2, circ,circ2,r2])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (shape);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "geo = OCCGeometry(shape, dim=2)\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.2)).Curve(3)\n",
    "Draw (mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "### outgoing plane waves:\n",
    "$$\n",
    "\\sin (k_{y,j} y) e^{i k_{x,j} x}\n",
    "$$\n",
    "with\n",
    "\\begin{eqnarray*}\n",
    "k_{y,j} & = & j \\pi \\\\\n",
    "k_{x,j} & = & \\sqrt{\\omega^2 - k_{y,j}^2}\n",
    "\\end{eqnarray*}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "Define them as multidimensional `CoefficientFunction`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "omega = 3.5*pi * (1+0j)\n",
    "\n",
    "ky = [j*pi for j in range(1,5)]\n",
    "kx = [sqrt(omega**2-kyi**2) for kyi in ky]\n",
    "k = zip(ky,kx)\n",
    "shapes = CF(tuple(sin(kyi*y)*exp(1j*kxi*x) for kyi,kxi in k))\n",
    "dshapesx = shapes.Diff(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (shapes[2], mesh, animate_complex=True );  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve.comp import GlobalSpace\n",
    "\n",
    "fesfem = H1(mesh, order=4, definedon=\"fem.*\", dirichlet=\"dir\", complex=True)\n",
    "profile=sin(pi*y)\n",
    "feswaves = GlobalSpace (mesh, definedon=\"waves\", basis = shapes) \n",
    "feswaves.AddOperator(\"dn\", BND, dshapesx) \n",
    "fes = fesfem*feswaves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes.ndof, fesfem.ndof, feswaves.ndof"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12",
   "metadata": {},
   "source": [
    "Bilinear-form:\n",
    "$$\n",
    "\\int_{\\Omega_{FEM}} \\nabla u \\nabla v - \\omega^2 \\varepsilon u v  dx +\n",
    "\\int_{\\gamma_{coupling}} - u \\partial_n v_w - v \\partial_n u_w + u_w \\partial_n v_w\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "u,v = fes.TnT()\n",
    "(uf,uw), (vf,vw) = fes.TnT()\n",
    "uwdx = uw.Operator(\"dn\")\n",
    "vwdx = vw.Operator(\"dn\")\n",
    "\n",
    "epsilon = mesh.MaterialCF( { \"obstacle\" : 10 }, 1)\n",
    "\n",
    "a = BilinearForm(fes)\n",
    "a += (grad(uf))*grad(vf)*dx(\"fem.*\") - omega**2*epsilon*uf*vf*dx(\"fem.*\")\n",
    "a += (-uf*vwdx-vf*uwdx+uw.Trace()*vwdx)*ds(\"coupling\")\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(1*vf*dx(\"femsource\")).Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes)\n",
    "gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec\n",
    "\n",
    "gfu2 = mesh.MaterialCF( { \"fem.*\" : gfu.components[0], \"waves\" : gfu.components[1] })\n",
    "Draw (gfu2, mesh, animate_complex=True );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}