{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maxwell Equations\n",
    "==="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.csg import *\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Geometric model using Netgen constructive solid geometry:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MakeGeometry():\n",
    "    geometry = CSGeometry()\n",
    "    box = OrthoBrick(Pnt(-1,-1,-1),Pnt(2,1,2)).bc(\"outer\")\n",
    "\n",
    "    core = OrthoBrick(Pnt(0,-0.05,0),Pnt(0.8,0.05,1))- \\\n",
    "           OrthoBrick(Pnt(0.1,-1,0.1),Pnt(0.7,1,0.9))- \\\n",
    "           OrthoBrick(Pnt(0.5,-1,0.4),Pnt(1,1,0.6)).maxh(0.2).mat(\"core\")\n",
    "    \n",
    "    coil = (Cylinder(Pnt(0.05,0,0), Pnt(0.05,0,1), 0.3) - \\\n",
    "            Cylinder(Pnt(0.05,0,0), Pnt(0.05,0,1), 0.15)) * \\\n",
    "            OrthoBrick (Pnt(-1,-1,0.3),Pnt(1,1,0.7)).maxh(0.2).mat(\"coil\")\n",
    "    \n",
    "    geometry.Add ((box-core-coil).mat(\"air\"), transparent=True)\n",
    "    geometry.Add (core, col=(0.5,0.5,0))\n",
    "    geometry.Add (coil, col=(0,1,0))\n",
    "    return geometry\n",
    "\n",
    "geo = MakeGeometry()\n",
    "# Draw (geo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(geo.GenerateMesh(maxh=0.5))\n",
    "mesh.Curve(5)\n",
    "Draw (mesh, clipping = { \"pnt\" : (0,0,0), \"vec\" : (0,1,0) })"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Magnetostatics:\n",
    "\n",
    "\n",
    "find $u \\in H(curl)$ such that\n",
    "$$\n",
    "\\int \\mu^{-1} \\operatorname{curl} u \\operatorname{curl} v = \\int j v\n",
    "$$\n",
    "for all $v \\in H(curl)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = HCurl(mesh, order=4, dirichlet=\"outer\", nograds = True)\n",
    "print (\"ndof =\", fes.ndof)\n",
    "u,v = fes.TnT()\n",
    "\n",
    "mur = { \"core\" : 1000, \"coil\" : 1, \"air\" : 1 }\n",
    "mu0 = 1.257e-6\n",
    "nu_coef = [ 1/(mu0*mur[mat]) for mat in mesh.GetMaterials() ]\n",
    "\n",
    "nu = CoefficientFunction(nu_coef)\n",
    "a = BilinearForm(fes, symmetric=True)\n",
    "a += nu*curl(u)*curl(v)*dx + 1e-6*nu*u*v*dx\n",
    "\n",
    "c = Preconditioner(a, type=\"bddc\")\n",
    "\n",
    "f = LinearForm(fes)\n",
    "f += CoefficientFunction((y,0.05-x,0)) * v * dx(\"coil\")\n",
    "\n",
    "u = GridFunction(fes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Distribute work by task-parallelization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager():\n",
    "    a.Assemble()\n",
    "    f.Assemble()\n",
    "    solver = CGSolver(mat=a.mat, pre=c.mat)\n",
    "    u.vec.data = solver * f.vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (curl(u), mesh, \"B-field\", draw_surf=False, \\\n",
    "      clipping = { \"pnt\" : (0,0,0), \"vec\" : (0,1,0), \"function\" : False },\n",
    "      vectors = { \"grid_size\" : 100 })"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}