{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# 11.3.1 Helmholtz solver using Brakhage-Werner formulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "Combined field integral equations combine single and double layer integral operators, one simple option is the Brakhage-Werner formulation.\n",
    "\n",
    "The solution is represented as\n",
    "\n",
    "$$\n",
    "u = (i \\kappa S - D) \\phi,\n",
    "$$\n",
    "\n",
    "where $\\phi$ solve the boundary integral equation\n",
    "$$\n",
    "    \\big( \\tfrac{1}{2} + K + i \\kappa V \\big) \\phi = u_{in} \\qquad \\text{on} \\, \\Gamma\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.occ import *\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from ngsolve.bem import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "kappa=10\n",
    "order=4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "screen = WorkPlane(Axes( (0,0,0), Z, X)).RectangleC(15,15).Face()\n",
    "\n",
    "sp = Fuse(Sphere( (0,0,0), pi).faces)\n",
    "screen.faces.name=\"screen\"\n",
    "sp.faces.name=\"sphere\"\n",
    "shape = Compound([screen,sp])\n",
    "\n",
    "mesh = shape.GenerateMesh(maxh=5/kappa).Curve(order)\n",
    "Draw (mesh);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes_sphere = Compress(SurfaceL2(mesh, order=order, complex=True, definedon=mesh.Boundaries(\"sphere\")))\n",
    "u,v = fes_sphere.TnT()\n",
    "fes_screen = Compress(SurfaceL2(mesh, order=order, dual_mapping=True, complex=True, definedon=mesh.Boundaries(\"screen\")))\n",
    "print (\"ndof_sphere = \", fes_sphere.ndof, \"ndof_screen =\", fes_screen.ndof)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager():\n",
    "    # V = HelmholtzSingleLayerPotentialOperator(fes_sphere, fes_sphere, kappa=kappa, intorder=10)\n",
    "    # K = HelmholtzDoubleLayerPotentialOperator(fes_sphere, fes_sphere, kappa=kappa, intorder=10)\n",
    "    # C = HelmholtzCombinedFieldOperator(fes_sphere, fes_sphere, kappa=kappa, intorder=10)\n",
    "    C = HelmholtzCF(u*ds(\"sphere\"), kappa)*v*ds\n",
    "    u,v  = fes_sphere.TnT()\n",
    "    Id = BilinearForm(u*v*ds).Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "lhs = 0.5 * Id.mat + C.mat\n",
    "source = exp(1j * kappa * x) \n",
    "rhs = LinearForm(-source*v*ds).Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes_sphere)\n",
    "pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()\n",
    "with TaskManager():\n",
    "    gfu.vec[:] = solvers.GMRes(A=lhs, b=rhs.vec, pre=pre, maxsteps=40, tol=1e-8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (gfu, order=5, min=-1, max=1);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "## prostprocessing on screen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "uscat = GridFunction(fes_screen)\n",
    "with TaskManager(pajetrace=10**8):\n",
    "    # uscat.Set(1j*kappa*V.GetPotential(gfu)-K.GetPotential(gfu) , definedon=mesh.Boundaries(\"screen\"))\n",
    "    # uscat.Set(HelmholtzCF(u*ds(\"sphere\"), kappa)(gfu)  , definedon=mesh.Boundaries(\"screen\"))\n",
    "    uscat.Set(1j*kappa*HelmholtzSL(u*ds(\"sphere\"),kappa)(gfu) - HelmholtzDL(u*ds(\"sphere\"), kappa)(gfu), \\\n",
    "              definedon=mesh.Boundaries(\"screen\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "print (\"Scattered field\")\n",
    "Draw (uscat, mesh, min=-1,max=1, animate_complex=True, order=4);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "uin = mesh.BoundaryCF( {\"screen\": source }, default=0)\n",
    "print (\"Total field\")\n",
    "Draw (uin-uscat, mesh, min=-1,max=1, animate_complex=True, order=4);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "Scattering from sphere with $D = 50 \\lambda$. About 5 min on Macbook Apple M4 Pro\n",
    "\n",
    "<img src=\"wave50.png\" alt=\"Alternative text\" width=\"800\" align=\"center\"/>\n"
   ]
  },
  {
   "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.13.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}