{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2b9804e9-50a1-4d32-9fa9-b0ff5668355a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.occ import *\n",
    "from netgen.meshing import BoundaryLayerParameters\n",
    "import numpy as np\n",
    "from netgen.webgui import Draw as DrawGeo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "bbd8172e-b0ce-45c9-8f28-9c8fc13de073",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of boundary layers\n",
    "number_of_layers = 1\n",
    "# Order of elements\n",
    "Order=2\n",
    "# Choice of solver type\n",
    "#IterativeSolver=\"Direct\" # works ok with or without prismatic layers.\n",
    "IterativeSolver=\"Iterative\" # only appears to work without prismatic layers.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9d6a7d96-8002-427f-9a8c-9738d599e912",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d702723c9a204ded9381ae1a4ab8ae82",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='50vh', width='100%'), value={'gui_settings': {}, 'mesh_dim': 3, 'mesh_cente…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "BaseWebGuiScene"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Setup the mesh and materials (sphere of unit radius in a large box)\n",
    "\n",
    "material_name = ['object']\n",
    "mur = {\"object\": 4, \"air\":1 }\n",
    "sigma = {\"object\": 1e6}\n",
    "inorout = {\"object\": 1, \"air\":0 }\n",
    "alpha = 0.01\n",
    "\n",
    "# Info for boundary layers on sphere\n",
    "max_target_frequency = 1e8\n",
    "boundary_layer_material = material_name[0]\n",
    "\n",
    "\n",
    "# sphere\n",
    "sphere = Sphere(Pnt(0,0,0), r=1)\n",
    "pos_sphere = sphere - Box(Pnt(0,100,100), Pnt(-100,-100,-100))\n",
    "neg_sphere = sphere - Box(Pnt(0,100,100), Pnt(100,-100,-100))\n",
    "sphere = pos_sphere + neg_sphere\n",
    "sphere.bc('default')\n",
    "sphere.mat(material_name[0])\n",
    "sphere.maxh = 0.2\n",
    "\n",
    "# box\n",
    "dim=1e3 # \n",
    "box = Box(Pnt(-dim, -dim, -dim), Pnt(dim,dim,dim))\n",
    "box.mat('air')\n",
    "box.bc('outer')\n",
    "box.maxh=dim\n",
    "box=box-sphere\n",
    "\n",
    "# Glue joins two OCC objects together without interior elemements\n",
    "joined_object = Glue([sphere, box])\n",
    "\n",
    "mu0 = 4 * np.pi * 1e-7\n",
    "tau = (2/(max_target_frequency * sigma[\"object\"] * mu0 * mur[\"object\"]))**0.5 / alpha\n",
    "if number_of_layers > 0:\n",
    "    layer_thicknesses = [(2**n)*tau for n in range(number_of_layers)]\n",
    "\n",
    "    B = BoundaryLayerParameters(boundary=\".*\", thickness=layer_thicknesses, new_material=boundary_layer_material,\n",
    "        domain=boundary_layer_material, outside=False, disable_curving=False )\n",
    "    nmesh = OCCGeometry(joined_object).GenerateMesh(boundary_layers=[B])\n",
    "else:\n",
    "    nmesh = OCCGeometry(joined_object).GenerateMesh()\n",
    "\n",
    "from ngsolve import *\n",
    "mesh = Mesh(nmesh)\n",
    "DrawGeo(nmesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ee8bc69f-604b-4f5a-acff-2c4531667b69",
   "metadata": {},
   "outputs": [],
   "source": [
    "dom_nrs_metal = [0 if mat == \"air\" else 1 for mat in mesh.GetMaterials()]\n",
    "fes = HCurl(mesh, order=Order, dirichlet=\"outer\", gradientdomains=dom_nrs_metal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "af460feb-5690-4e21-b359-ae537be51ac5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "assembled matrices\n"
     ]
    }
   ],
   "source": [
    "# Coefficient functions\n",
    "mu_coef = [mur[mat] for mat in mesh.GetMaterials()]\n",
    "mu = CoefficientFunction(mu_coef)\n",
    "inout_coef = [inorout[mat] for mat in mesh.GetMaterials()]\n",
    "inout = CoefficientFunction(inout_coef)\n",
    "\n",
    "\n",
    "u = fes.TrialFunction()\n",
    "v = fes.TestFunction()\n",
    "\n",
    "# Weak form\n",
    "a = BilinearForm(fes)#,symmetric=True)\n",
    "a += SymbolicBFI((1/mu)*curl(u)*curl(v))\n",
    "\n",
    "mreg = BilinearForm(fes)#,symmetric=True)\n",
    "epsi=1e-8\n",
    "mreg += SymbolicBFI(inout*u*v) \n",
    "mreg += SymbolicBFI(epsi*u*v*(1-inout))\n",
    "\n",
    "apre = BilinearForm(fes)#,symmetric=True)\n",
    "apre += SymbolicBFI((1/mu)*curl(u)*curl(v))\n",
    "apre += SymbolicBFI(u*v*inout)\n",
    "apre += SymbolicBFI(epsi*u*v*(1-inout))\n",
    "\n",
    "pre = Preconditioner(apre, \"bddc\")\n",
    "\n",
    "with TaskManager():\n",
    "        a.Assemble()\n",
    "        mreg.Assemble()\n",
    "        apre.Assemble()\n",
    "        pre.Update()\n",
    "        \n",
    "        # build gradient matrix as sparse matrix (and corresponding scalar FESpace)\n",
    "        gradmat, fesh1 = fes.CreateGradient()\n",
    "\n",
    "        gradmattrans = gradmat.CreateTranspose() # transpose sparse matrix\n",
    "        math1 = gradmattrans @ mreg.mat @ gradmat   # multiply matrices\n",
    "\n",
    "        print(\"assembled matrices\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "22c02033-93aa-48bc-ab78-4b3a8a9c34cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "if IterativeSolver==\"Direct\":\n",
    "    invh1 = math1.Inverse(inverse=\"sparsecholesky\", freedofs=fesh1.FreeDofs()) # Note use of free DOFs only\n",
    "else:\n",
    "    math1smooth = math1.CreateSmoother(freedofs=fesh1.FreeDofs())\n",
    "\n",
    "    cgmath1 = CGSolver(\n",
    "                    mat=math1,\n",
    "                    pre=math1smooth,\n",
    "                    maxsteps=300,\n",
    "                    precision=1e-8\n",
    "                    )\n",
    "\n",
    "    class H1Inverse(BaseMatrix):\n",
    "\n",
    "                    def Mult(self, x, y):\n",
    "                        y.data = cgmath1 * x\n",
    "\n",
    "                    def CreateColVector(self):\n",
    "                        return math1.CreateColVector()\n",
    "\n",
    "                    def CreateRowVector(self):\n",
    "                        return math1.CreateRowVector()\n",
    "\n",
    "                    def Height(self):\n",
    "                        return math1.height\n",
    "\n",
    "                    def Width(self):\n",
    "                        return math1.width\n",
    "\n",
    "    invh1 = H1Inverse()\n",
    "# build the Poisson projector with operator Algebra:\n",
    "proj = IdentityMatrix() - gradmat @ invh1 @ gradmattrans @ mreg.mat\n",
    "projpre = proj @ pre.mat\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d8c1c92-8921-49b0-9871-d1b4aaf4244d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 : [3.6124009024222824e-05, 0.07174821886269699, 0.08877170763081393, 0.09245105361139802, 0.1025221394492977, 0.1081084472564487, 0.11875016749607516, 0.1272307406837432, 0.1311511392302223, 0.14115324328622145, 0.14746967958808027, 0.14994922553951784, 0.15318015122842366, 0.16480746756023382, 0.1685706790236387, 0.18665742193269094, 0.1949132238412333, 0.2071935834880672, 0.218802328627858, 0.22096541161394218]\n",
      "1 : [1.5270158353052643e-05, 0.030772322476561153, 0.03735080823560294, 0.03915102575220435, 0.04322149304919363, 0.04599031543777567, 0.05091753635878178, 0.05431916874125278, 0.056868284874484414, 0.05971947644060693, 0.06299717337494921, 0.06466217764925049, 0.06674169954115959, 0.06955729842758826, 0.07210080402905208, 0.08143858996686806, 0.08426156156516235, 0.08838509358800213, 0.09501363968091422, 0.09838668317387761]\n",
      "2 : [7.652991055616832e-06, 0.015310842078069415, 0.018356705016815217, 0.019170543883189808, 0.021473651045241147, 0.022731784530217583, 0.02506240271297434, 0.027039810188115868, 0.028911590489956594, 0.029422244431214984, 0.031185474860997438, 0.032472762192416654, 0.03360671472099675, 0.03551761216197608, 0.03608045735601893, 0.04089820457365759, 0.04241791479222807, 0.04394237887751302, 0.04726912709442861, 0.050886764244926906]\n"
     ]
    }
   ],
   "source": [
    "evals, evecs = solvers.PINVIT(a.mat, mreg.mat, pre=projpre, num=20, maxit=20)\n",
    "print(evals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1647c942-5e1a-4703-aa4d-9827c9104293",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9f05563-242f-4fcd-9ea4-83dcacaef9a4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c60663af-87d5-4a82-bd62-03e11b4b8892",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a87413e7-7897-4c64-94de-fb31135376fe",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24a89e0e-a82f-4635-a76b-522fe48edf43",
   "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
}