{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.4.1 Maxwell eigenvalue problem\n",
    "We solve the Maxwell eigenvalue problem\n",
    "\n",
    "$$\n",
    "\\int \\operatorname{curl} u \\, \\operatorname{curl} v\n",
    "= \\lambda \\int u v \n",
    "$$\n",
    "\n",
    "for $u, v \\; \\bot \\; \\nabla H^1$\n",
    "using a PINVIT solver from the ngsolve solvers module.\n",
    "\n",
    "The orthogonality to gradient fields is important to eliminate the huge number of zero eigenvalues. The orthogonal sub-space is implemented using a Poisson projection:\n",
    "\n",
    "$$\n",
    "P u = u - \\nabla \\Delta^{-1} \\operatorname{div} u\n",
    "$$\n",
    "\n",
    "The algorithm and example is take form the Phd thesis of Sabine Zaglmayr, p 145-150."
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.occ import *\n",
    "\n",
    "cube1 = Box( (-1,-1,-1), (1,1,1) )\n",
    "\n",
    "cube2 = Box( (0,0,0), (2,2,2) )\n",
    "cube2.edges.hpref=1    # mark edges for geometric refinement\n",
    "\n",
    "fichera = cube1-cube2\n",
    "\n",
    "Draw (fichera);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(OCCGeometry(fichera).GenerateMesh(maxh=0.4))\n",
    "mesh.RefineHP(levels=2, factor=0.2)\n",
    "Draw (mesh);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SetHeapSize(100*1000*1000)\n",
    "\n",
    "fes = HCurl(mesh, order=3)\n",
    "print (\"ndof =\", fes.ndof)\n",
    "u,v = fes.TnT()\n",
    "\n",
    "a = BilinearForm(curl(u)*curl(v)*dx)\n",
    "m = BilinearForm(u*v*dx)\n",
    "\n",
    "apre = BilinearForm(curl(u)*curl(v)*dx + u*v*dx)\n",
    "pre = Preconditioner(apre, \"direct\", inverse=\"sparsecholesky\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager():\n",
    "    a.Assemble()\n",
    "    m.Assemble()\n",
    "    apre.Assemble()\n",
    "\n",
    "    # build gradient matrix as sparse matrix (and corresponding scalar FESpace)\n",
    "    gradmat, fesh1 = fes.CreateGradient()\n",
    "    \n",
    "    \n",
    "    gradmattrans = gradmat.CreateTranspose() # transpose sparse matrix\n",
    "    math1 = gradmattrans @ m.mat @ gradmat   # multiply matrices \n",
    "    math1[0,0] += 1     # fix the 1-dim kernel\n",
    "    invh1 = math1.Inverse(inverse=\"sparsecholesky\")\n",
    "\n",
    "    # build the Poisson projector with operator Algebra:\n",
    "    proj = IdentityMatrix() - gradmat @ invh1 @ gradmattrans @ m.mat\n",
    "\n",
    "    projpre = proj @ pre.mat\n",
    "\n",
    "    evals, evecs = solvers.PINVIT(a.mat, m.mat, pre=projpre, num=12, maxit=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (\"Eigenvalues\")\n",
    "for lam in evals: \n",
    "    print (lam)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes, multidim=len(evecs))\n",
    "for i in range(len(evecs)):\n",
    "    gfu.vecs[i].data = evecs[i]\n",
    "Draw (Norm(gfu), mesh, order=4, min=0, max=2);"
   ]
  },
  {
   "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.13.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}