{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.4 Maxwell's Equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Peter Monk: \"Finite Elements for Maxwell's Equations\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Magnetostatic field generated by a permanent magnet\n",
    "\n",
    "magnetic flux $B$, magnetic field $H$, given magnetization $M$:\n",
    "\n",
    "$$\n",
    "\\DeclareMathOperator{\\Grad}{grad}\n",
    "\\DeclareMathOperator{\\Curl}{curl}\n",
    "\\DeclareMathOperator{\\Div}{div}\n",
    "B = \\mu (H + M), \\quad \\Div B = 0, \\quad \\Curl H = 0\n",
    "$$\n",
    "\n",
    "Introducing a vector-potential $A$ such that $B = \\Curl A$, and putting equations together we get\n",
    "\n",
    "$$\n",
    "\\Curl \\mu^{-1} \\Curl A = \\Curl M\n",
    "$$\n",
    "\n",
    "In weak form: Find $A \\in H(\\Curl)$ such that\n",
    "\n",
    "$$\n",
    "\\int \\mu^{-1} \\Curl A \\Curl v = \\int M \\Curl v \\qquad  \\forall \\, v \\in H(\\Curl)\n",
    "$$\n",
    "\n",
    "Usually, the permeability $\\mu$ is given as \n",
    "$\\mu = \\mu_r \\mu_0$, with $\\mu_0 = 4 \\pi 10^{-7}$ the permeability of vacuum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.occ import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Geometric model and meshing of a bar magnet:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# box = OrthoBrick(Pnt(-3,-3,-3),Pnt(3,3,3)).bc(\"outer\")\n",
    "# magnet = Cylinder(Pnt(-1,0,0),Pnt(1,0,0), 0.3) * OrthoBrick(Pnt(-1,-3,-3),Pnt(1,3,3))\n",
    "# air = box - magnet\n",
    "box = Box( (-3,-3,-3), (3,3,3))\n",
    "box.faces.name = \"outer\"\n",
    "\n",
    "magnet = Cylinder((-1,0,0),X, r=0.3, h=2)\n",
    "magnet.mat(\"magnet\")\n",
    "magnet.faces.col = (1,0,0)\n",
    "\n",
    "air = box-magnet\n",
    "air.mat(\"air\")\n",
    "shape = Glue([air,magnet])\n",
    "geo = OCCGeometry(shape)\n",
    "\n",
    "Draw (shape, clipping={ \"z\" : -1, \"function\":True})\n",
    "                       \n",
    "mesh = Mesh(geo.GenerateMesh(maxh=2, curvaturesafety=1))\n",
    "mesh.Curve(3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh.GetMaterials(), mesh.GetBoundaries()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define space, forms and preconditioner.\n",
    "\n",
    "* To obtain a regular system matrix, we regularize by adding a very small $L_2$ term.\n",
    "* We solve magnetostatics, so we can gauge by adding and arbitrary gradient field. A cheap possibility is to delete all basis-functions which are gradients (flag 'nograds')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = HCurl(mesh, order=3, dirichlet=\"outer\", nograds=True)\n",
    "print (\"ndof =\", fes.ndof)\n",
    "u,v = fes.TnT()\n",
    "\n",
    "from math import pi\n",
    "mu0 = 4*pi*1e-7\n",
    "mur = mesh.MaterialCF({\"magnet\" : 1000}, default=1)\n",
    "\n",
    "a = BilinearForm(fes)\n",
    "a += 1/(mu0*mur)*curl(u)*curl(v)*dx + 1e-8/(mu0*mur)*u*v*dx\n",
    "c = Preconditioner(a, \"bddc\")\n",
    "\n",
    "f = LinearForm(fes)\n",
    "mag = mesh.MaterialCF({\"magnet\" : (1,0,0)}, default=(0,0,0))\n",
    "f += mag*curl(v) * dx(\"magnet\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assemble system and setup preconditioner using task-parallelization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager():\n",
    "    a.Assemble()\n",
    "    f.Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, declare GridFunction and solve by preconditioned CG iteration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes)\n",
    "with TaskManager():\n",
    "    solvers.CG(sol=gfu.vec, rhs=f.vec, mat=a.mat, pre=c.mat, printrates=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (curl(gfu), mesh, \"B-field\", draw_surf=False, \\\n",
    "      clipping = { \"z\" : -1, \"function\":True}, \\\n",
    "      vectors = { \"grid_size\":50}, min=0, max=2e-5);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.10.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}