{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Magnetostatics\n",
    "\n",
    "Computing the magnetic field induced by a wire coil\n",
    "\n",
    "* first, we solve an electric conductivity problem in the wire\n",
    "* then, we use the computed current to solve (a reduced) Maxwell equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.occ import *\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "## model of the coil:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "cyl = Cylinder((0,0,0), Z, r=0.01, h=0.03).faces[0]\n",
    "heli = Edge(Segment((0,0), (12*pi, 0.03)), cyl)\n",
    "ps = heli.start\n",
    "vs = heli.start_tangent\n",
    "pe = heli.end\n",
    "ve = heli.end_tangent\n",
    "\n",
    "e1 = Segment((0,0,-0.03), (0,0,-0.01))\n",
    "c1 = BezierCurve( [(0,0,-0.01), (0,0,0), ps-vs, ps])\n",
    "e2 = Segment((0,0,0.04), (0,0,0.06))\n",
    "c2 = BezierCurve( [pe, pe+ve, (0,0,0.03), (0,0,0.04)])\n",
    "spiral = Wire([e1, c1, heli, c2, e2])\n",
    "circ = Face(Wire([Circle((0,0,-0.03), Z, 0.001)]))\n",
    "coil = Pipe(spiral, circ)\n",
    "\n",
    "coil.faces.maxh=0.2\n",
    "coil.faces.name=\"coilbnd\"\n",
    "coil.faces.Max(Z).name=\"in\"\n",
    "coil.faces.Min(Z).name=\"out\"\n",
    "coil.faces.col=(184/256, 115/256, 51/256)\n",
    "coil.mat(\"coil\")\n",
    "crosssection = coil.faces.Max(Z).mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "ea = { \"euler_angles\" : (-130, -73, 0), \"radius\" : 0.025 }\n",
    "Draw (coil, **ea);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "box = Box((-0.04,-0.04,-0.03), (0.04,0.04,0.06))\n",
    "box.faces.name = \"outer\"\n",
    "air = box-coil\n",
    "air.mat(\"air\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "## mesh-generation of coil and air-box:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "geo = OCCGeometry(Glue([coil,air]))\n",
    "with TaskManager():\n",
    "    mesh = Mesh(geo.GenerateMesh(meshsize.coarse, maxh=0.01)).Curve(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "clipping = { \"clipping\" : { \"y\":1, \"z\":0, \"dist\":0.012} }\n",
    "Draw (mesh, **clipping, **ea);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9",
   "metadata": {},
   "source": [
    "checking mesh data materials and boundaries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh.ne, mesh.nv, mesh.GetMaterials(), mesh.GetBoundaries()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "## Solve a potential problem to determine current density in wire:\n",
    "on the domain $\\Omega_{\\text{coil}}$, solve for potential $\\Phi$ and current density $j$:\n",
    "\n",
    "\\begin{align*}\n",
    "j & = \\sigma \\nabla \\Phi \\\\\n",
    "\\operatorname{div} j & = 0\n",
    "\\end{align*}\n",
    "with electric conductivity $\\sigma$.\n",
    "\n",
    "\n",
    "port boundary conditions: \n",
    "\\begin{align*}\n",
    "\\Phi & = 0  \\qquad \\qquad \\text{on } \\Gamma_{\\text{out}},  \\\\\n",
    "j_n & = \\frac{1}{|\\Gamma_{in}|} \\quad \\qquad \\text{on } \\Gamma_{\\text{in}},\n",
    "\\end{align*}\n",
    "\n",
    "and $j_n=0$ else"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "fespot = H1(mesh, order=3, definedon=mesh.Materials(\"coil\"), dirichlet=\"out\")\n",
    "phi,psi = fespot.TnT()\n",
    "sigma = 58.7e6\n",
    "with TaskManager():\n",
    "    bfa = BilinearForm(sigma*grad(phi)*grad(psi)*dx).Assemble()\n",
    "    inv = bfa.mat.Inverse(freedofs=fespot.FreeDofs(), inverse=\"sparsecholesky\")\n",
    "    lff = LinearForm(1/crosssection*psi*ds(\"in\")).Assemble()\n",
    "    gfphi = GridFunction(fespot)\n",
    "    gfphi.vec.data = inv * lff.vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (gfphi, draw_vol=False, **clipping, **ea);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14",
   "metadata": {},
   "source": [
    "## Solve magnetostatic problem:\n",
    "\n",
    "current source is current from potential equation: find $u \\in H(\\operatorname{curl})$:\n",
    "\n",
    "$$\n",
    "\\int \\mu^{-1} \\operatorname{curl} u \\cdot \\operatorname{curl} v \\, dx =\n",
    "\\int j \\cdot v \\, dx\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = HCurl(mesh, order=2, nograds=True)\n",
    "print (\"HCurl dofs:\", fes.ndof)\n",
    "u,v = fes.TnT()\n",
    "mu = 4*pi*1e-7\n",
    "a = BilinearForm(1/mu*curl(u)*curl(v)*dx+1e-6/mu*u*v*dx)\n",
    "pre = preconditioners.BDDC(a)\n",
    "f = LinearForm(sigma*grad(gfphi)*v*dx(\"coil\"))\n",
    "with TaskManager():\n",
    "    a.Assemble()\n",
    "    f.Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "inv = solvers.CGSolver(a.mat, pre, plotrates=True)\n",
    "gfu = GridFunction(fes)\n",
    "with TaskManager():\n",
    "    gfu.vec.data = inv * f.vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "s = 0.04*2\n",
    "N = 15 \n",
    "p = [(-s+2*s*i/N,-s+2*s*j/N,-s+2*s*k/N) for i in range(1,N) for j in range(1,N) for k in range(1,N)]\n",
    "# \n",
    "fieldlines = curl(gfu)._BuildFieldLines(mesh, p, num_fieldlines=N**3//5, randomized=True, length=0.3)\n",
    "from ngsolve.webgui import FieldLines\n",
    "# fieldlines = FieldLines(curl(gfu), mesh.Materials(\".*\"), length=0.2, num_lines=100)\n",
    "\n",
    "# fieldlines = FieldLines(curl(gfu), mesh=mesh, start_points=p, length=0.2, num_lines=100)\n",
    "\n",
    "Draw(curl(gfu), mesh,  \"X\", draw_vol=False, draw_surf=True, objects=[fieldlines], \\\n",
    "     min=0, max=1e-4, autoscale=False, settings={\"Objects\": {\"Surface\": False}},\n",
    "    **ea, **clipping);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}