{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# 11.1.4 Biot Savart\n",
    "\n",
    "Let $G(x,y) = \\frac{1}{4 \\pi} \\frac{\\exp (i k |x-y|)}{|x-y|}$ be Green's function for the Helmholtz equation. \n",
    "\n",
    "For a given current path $j$ along a Curve $C$, the magnetic field in vacuum on the full space ${\\mathbb R}^3$ is\n",
    "\n",
    "$$\n",
    "H(x) = \\int_C  j(y) \\times \\nabla_y G(x,y) dl_y\n",
    "$$"
   ]
  },
  {
   "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\n",
    "\n",
    "vismesh = Mesh(OCCGeometry(Box((-5,-5,-5), (5,5,5))).GenerateMesh(maxh=1))\n",
    "for l in range(0):\n",
    "    vismesh.Refine()\n",
    "Draw (vismesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve.bem import BiotSavartCF, BiotSavartSingularMLCF, BiotSavartRegularMLCF"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "We evaluate the single layer integral using numerical integration on the surface mesh. Thus, we get a sum of many Green's functions, which is compressed using a multilevel-multipole. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "kappa = 0.01*pi\n",
    "mp = BiotSavartSingularMLCF((0,0,0), r=5, kappa=kappa)\n",
    "mp.expansion.AddCurrent( (1,0,-1), (1,0,1), 1, num=100)\n",
    "mp.expansion.AddCurrent( (1,0, 1), (-1,0,1), 1, num=100)\n",
    "mp.expansion.AddCurrent( (-1,0, 1), (-1,0,-1), 1, num=10)\n",
    "mp.expansion.AddCurrent( (-1,0, -1), (1,0,-1), 1, num=100)\n",
    "\n",
    "regmp = mp.CreateRegularExpansion((0,0,0),r=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "clipping = { \"function\" : False,  \"pnt\" : (0,0,0), \"vec\" : (0,0,-1) }\n",
    "Draw (regmp.real, vismesh, min=0, max=1, order=2,  vectors={\"grid_size\" : 40, \"offset\" : 0 }, clipping=clipping);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "visplane = WorkPlane(Axes( (0,0.2,0), Z, X)).RectangleC(5,5).Face()\n",
    "vismesh2 = Mesh(OCCGeometry(visplane).GenerateMesh(maxh=0.5))\n",
    "Draw (regmp.real[0], vismesh2, min=-0.1, max=0.1, order=10);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "kappa = 0.01*pi\n",
    "mp = BiotSavartSingularMLCF((0,0,0), r=5, kappa=kappa)\n",
    "\n",
    "coil = Cylinder((0,-1,0), Y, r=1, h=1, mantle=\"outer\") - Cylinder((0,-1,0), Y, r=0.5, h=1)\n",
    "coilmesh = Mesh(OCCGeometry(coil).GenerateMesh(maxh=0.3)).Curve(3)\n",
    "Draw (coilmesh)\n",
    "\n",
    "current = CF((z,0,-x))\n",
    "current /= Norm(current)\n",
    "mp.expansion.AddCurrentDensity(current, coilmesh.Materials(\".*\"))\n",
    "# mp.expansion.AddCurrentDensity(current, coilmesh.Boundaries(\"outer\"))\n",
    "\n",
    "regmp = mp.CreateRegularExpansion((0,0,0),r=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "clipping = { \"function\" : False,  \"pnt\" : (0,0,0), \"vec\" : (0,0,-1) }\n",
    "Draw (regmp.real, vismesh, min=0, max=0.2, order=2,  vectors={\"grid_size\" : 40, \"offset\" : 0 }, clipping=clipping);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (regmp.real[0], vismesh2, min=-0.1, max=0.1, order=10);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "the spikes *inside* the coil domain stem from numerical integration of the current source. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve.webgui import FieldLines\n",
    "N=14\n",
    "fl = FieldLines(mp.real, vismesh.Materials('.*'), num_lines=N**3/20, length=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "clipping = { \"function\" : False,  \"pnt\" : (0,0,0), \"vec\" : (0,0,-1) }\n",
    "Draw (regmp.real, vismesh, min=0, max=0.5, order=2,  vectors={\"grid_size\" : 40, \"offset\" : 0 }, clipping=clipping);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "settings = {\"objects\": { \"Surface\":False, \"Wireframe\":False}}\n",
    "clipping = {\"function\" : False, \"pnt\" : (0,0,0), \"vec\" : (0,0,-1) }\n",
    "vectors = {\"grid_size\" : 100, \"offset\": 0.2}\n",
    "Draw (mp.real, vismesh, \"X\", objects=[fl], min=0, max=0.5, autoscale=False, settings=settings, vectors=vectors, clipping=clipping);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "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
}