{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "FEM-BEM Coupling\n",
    "==============\n",
    "\n",
    "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.\n",
    "\n",
    "In this demo we simulate a plate capacitor on an unbounded domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from netgen.occ import *\n",
    "from ngsolve.solvers import GMRes\n",
    "from ngsolve.webgui import Draw\n",
    "from ngsolve.bem import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "largebox = Box ((-2,-2,-2), (2,2,2) )\n",
    "eltop = Box ( (-1,-1,0.5), (1,1,1) )\n",
    "elbot = Box ( (-1,-1,-1), (1,1,-0.5))\n",
    "\n",
    "largebox.faces.name = \"outer\" # coupling boundary\n",
    "eltop.faces.name = \"topface\" # Dirichlet boundary\n",
    "elbot.faces.name = \"botface\" # Dirichlet boundary\n",
    "eltop.edges.hpref = 1\n",
    "elbot.edges.hpref = 1\n",
    "\n",
    "shell = largebox-eltop-elbot # FEM domain \n",
    "shell.solids.name = \"air\"\n",
    "\n",
    "mesh = shell.GenerateMesh(maxh=0.8)\n",
    "mesh.RefineHP(2)\n",
    "ea = { \"euler_angles\" : (-67, 0, 110) } \n",
    "Draw (mesh, clipping={\"x\":1, \"y\":0, \"z\":0, \"dist\" : 1.1}, **ea);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "On the **exterior** domain $\\Omega^c$, the solution can be expressed by the representation formula:\n",
    "\n",
    "$$ 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\\,, $$ \n",
    "\n",
    "where $\\gamma_0 u = u$ and $\\gamma_1 u = \\frac{\\partial u}{\\partial n}$ are Dirichlet and Neumann traces.\n",
    "These traces are related by the Calderon projector\n",
    "\n",
    "$$  \\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) $$.\n",
    "\n",
    "The $V$, $K$ are the single layer and double layer potential operators, and $D$ is the hypersingular operator.\n",
    "\n",
    "On the FEM domain we have the variational formulation\n",
    "\n",
    "$$\n",
    "\\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})\n",
    "$$\n",
    "\n",
    "We use Calderon's represenataion formula for the Neumann trace:\n",
    "\n",
    "$$\n",
    "\\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})\n",
    "$$\n",
    "\n",
    "\n",
    "To get a closed system, we use also the first equation of the Calderon equations. \n",
    "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:\n",
    "  \n",
    "  $$ \\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) \\,. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "Generate the finite element space for $H^1(\\Omega)$ and set the given Dirichlet boundary conditions on the surfaces of the plates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "order = 4\n",
    "fesH1 = H1(mesh, order=order, dirichlet=\"topface|botface\") \n",
    "print (\"H1-ndof = \", fesH1.ndof)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "fesL2 = SurfaceL2(mesh, order=order-1, dual_mapping=True, definedon=mesh.Boundaries(\"outer\")) \n",
    "print (\"L2-ndof = \", fesL2.ndof)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = fesH1 * fesL2\n",
    "u,dudn = fes.TrialFunction()\n",
    "v,dvdn = fes.TestFunction()\n",
    "\n",
    "a = BilinearForm(grad(u)*grad(v)*dx, check_unused=False).Assemble() \n",
    "\n",
    "gfudir = GridFunction(fes)\n",
    "gfudir.components[0].Set ( mesh.BoundaryCF( { \"topface\" : 1, \"botface\" : -1 }), BND)\n",
    "\n",
    "f = LinearForm(fes).Assemble()\n",
    "res = (f.vec - a.mat * gfudir.vec).Evaluate()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9",
   "metadata": {},
   "source": [
    "Generate the the single layer potential $V$, double layer potential $K$ and hypersingular operator $D$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = specialcf.normal(3)\n",
    "with TaskManager():\n",
    "    V = LaplaceSL(dudn*ds(\"outer\"))*dvdn*ds(\"outer\")\n",
    "    K = LaplaceDL(u*ds(\"outer\"))*dvdn*ds(\"outer\")\n",
    "    D = LaplaceSL(Cross(grad(u).Trace(),n)*ds(\"outer\"))*Cross(grad(v).Trace(),n)*ds(\"outer\")\n",
    "    M = BilinearForm(u*dvdn*ds(\"outer\"), check_unused=False).Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "Setup the coupled system matrix and the right hand side:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "sym = a.mat+D.mat - (0.5*M.mat+K.mat).T - (0.5*M.mat+K.mat) - V.mat\n",
    "rhs = res\n",
    "\n",
    "bfpre = BilinearForm(grad(u)*grad(v)*dx+1e-10*u*v*dx  + dudn*dvdn*ds(\"outer\") ).Assemble()\n",
    "pre = bfpre.mat.Inverse(freedofs=fes.FreeDofs(), inverse=\"sparsecholesky\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "Compute the solution of the coupled system:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager():\n",
    "    sol_sym = GMRes(A=sym, b=rhs, pre=pre, tol=1e-6, maxsteps=200, printrates=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes)\n",
    "gfu.vec[:] = sol_sym + gfudir.vec \n",
    "Draw(gfu.components[0], clipping={\"x\" : 1, \"y\":0, \"z\":0, \"dist\":0.0, \"function\" : True }, **ea, order=2); "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16",
   "metadata": {},
   "source": [
    "The Neumann data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (gfu.components[1], **ea);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18",
   "metadata": {},
   "source": [
    "**References:**\n",
    "\n",
    "- M. Costabel: [Principles of boundary element methods](https://pdf.sciencedirectassets.com/272705/1-s2.0-S0167797700X00068/1-s2.0-0167797787900141/main.pdf?X-Amz-Security-Token=IQoJb3JpZ2luX2VjEDQaCXVzLWVhc3QtMSJHMEUCIAUquWqfo6KKltdsWFjiJQlTK19Rp2RIAYPO8KSvFqbSAiEAgu56AezOrDSP1R6gQMPL%2B9KWRLYW5A5M0T1w4Y2RC00qvAUI7f%2F%2F%2F%2F%2F%2F%2F%2F%2F%2FARAFGgwwNTkwMDM1NDY4NjUiDIe%2Fn0sC87CTySnmmCqQBbf5We1Xu2u5VRqHY1MoJNp2fbRbYv5TXJXBw%2BAlELXUvtjL6kQQLfL%2FuHF%2FI2KQCy5ttpjLcg3P8QcjkhwZyjbvX6gpXmT91b7zizI8xy6EkVfrjxP0iagJO6EFGTw3WYLqotOxY4zWmKU7ud98ROPIDdKyYYWq8tKNmXUSXYAaZkU57bWq3Xc5GnETEQxOZOkTVZpvVkUm27HYuU4fKfzOxgbmo5Y3XPXYhjSBYCy%2BE4m3AwCTqbDbYYTj9fqTwBQeuOd%2FLuVtQdF9srJp3c%2FbZYtAzHHkgm1omg9nfF2PxS5u4x%2Ft7%2FDXKrFqHYiLgu3kqRQbnOy0at7pZ2YzGlLxf9cFGeot6M3%2BJtr%2BIugDj4aVzLgFnAQGDqMowFySYwSsjKNdH1QSmRzc0TMd35dI8OLEQR19omtvw%2BRZfwp4AQ%2BUIwhddL5JZagbvrPv8fRrzOQKau%2Bvq%2F9rjFTRF99fy8B6ZfxC8eB6hdgeUocrLJquPGca3u4XxwCyJM3LWRc2GIeOPiERKnso18n%2BDOQ7HcG61o2bttI3GqcHSvy2ZzE6J%2BLHBNu7fZKu65XcbofyLroXOilqJ7vUZcbcVSf44qofCfQcWCgnBnFGddsLjFHrrgLMKLmYW%2BT99%2BWRoUt8H5Hdaqf7xg41J0lTbNblNYaZQJuSo99uTRXZSzSchafApXOI1b84F%2BVo6GpI6tNdBlMNfbHHzuC3aPiSHt3cYWKNlrrR7u07tBhVpj7jKWo7aK3eUPmXN7qp6E6FO%2BLrUFD8jMGYwJoP18TkMh3y76leC3vooe9rfjymtOMC1C5vNng0KXLokiIuOnr0uQO5WFe8dxX1S3jW4cYwoI0GIEh1OpF%2FbgtiqH7qqCwqMIym3q0GOrEBwJoQ7RRJUF72%2FN%2FUhMl%2FEPq4CcQ0XjILMj4%2BLAlNFQ11%2Bw0dnvH5amm1JeWyxhsKB2Eytv4nvfDAezONlF3bpyyZ1wPec3USfEsuMDUzTi4JO1ifhIpBSv1x2C4IX9icEWH4f69RJHwUU6efONnUVB7b7Wy7%2B3sqMvTQAw1%2B8kmI2WqKr4iS2nHMNuPLE7qLIhyS%2BtFTclSnkidG5c4i%2F7AIUGF2DdYzAvfy2m9%2B0UP3&X-Amz-Algorithm=AWS4-HMAC-SHA256&X-Amz-Date=20240129T131454Z&X-Amz-SignedHeaders=host&X-Amz-Expires=300&X-Amz-Credential=ASIAQ3PHCVTYYOZJF7UN%2F20240129%2Fus-east-1%2Fs3%2Faws4_request&X-Amz-Signature=07de069e78def3fb3421bcdf984438efc15fa14475d662c9d6041f31c4fc9829&hash=071e43108fd6a95b78f41a85f877dc41fc07d17dec107a325cee6d820f08e663&host=68042c943591013ac2b2430a89b270f6af2c76d8dfd086a07176afe7c76c2c61&pii=0167797787900141&tid=spdf-a71b5532-b9d1-49eb-a596-5ba85ce194ee&sid=c44f2b54882a6847631b32a82fd7dc26326egxrqb&type=client&tsoh=d3d3LnNjaWVuY2VkaXJlY3QuY29t&ua=00025a5b0a530c0a07&rr=84d1bda7df9834b0&cc=de)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "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.13.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}