NGSolve24

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(advanced:Hcurl)=\n",
    "# Hcurl and related problems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NGSolve contains a variety of finite elements for apporximating different function spaces:\n",
    "\n",
    "- `H1`: Lagrange elements approximating functions in $H^1$ \n",
    "- `HCurl`: piecewise Nedelec elements approximating functions in $H(curl)$\n",
    "- `HDiv`: piecewise Brezzi-Douglas-Marini elements approximating functions in $H(div)$ \n",
    "\n",
    "One can loop over the basis function and observe that they are build in a hierarchical way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.occ import *\n",
    "\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.5,quad_dominated=True) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "order = 3\n",
    "fes = H1(mesh, order=order)\n",
    "gfu = GridFunction(fes )\n",
    "gfut = GridFunction(fes , multidim=0)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "for i in range(len(gfu.vec)):\n",
    "\n",
    "    gfu.vec[:] = 0\n",
    "    gfu.vec[i] = 1\n",
    "\n",
    "    gfut.AddMultiDimComponent(gfu.vec)\n",
    "\n",
    "settings = {\"camera\": {\"transformations\": [\n",
    "        {\"type\": \"rotateX\", \"angle\": -45}]}}\n",
    "\n",
    "Draw(gfut,deformation=True, animate=True ,autoscale = True , min=-0.03, max=0.03, order = 3, settings=settings);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "order = 2\n",
    "fes = HCurl(mesh, order=order)\n",
    "gfu = GridFunction(fes)\n",
    "gfut = GridFunction(fes, multidim=0)\n",
    "\n",
    "\n",
    "for i in range(len(gfu.vec)):\n",
    "\n",
    "    gfu.vec[:] = 0\n",
    "    gfu.vec[i] = 1\n",
    "\n",
    "    gfut.AddMultiDimComponent(gfu.vec)\n",
    "\n",
    "Draw(gfut, deformation=False, animate=True, autoscale = True , order = 3,  vectors = { \"grid_size\":30}, settings=settings);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "order = 2\n",
    "fes = HDiv(mesh, order=order)\n",
    "gfu = GridFunction(fes)\n",
    "gfut = GridFunction(fes, multidim=0)\n",
    "\n",
    "\n",
    "for i in range(len(gfu.vec)):\n",
    "\n",
    "    gfu.vec[:] = 0\n",
    "    gfu.vec[i] = 1\n",
    "\n",
    "    gfut.AddMultiDimComponent(gfu.vec)\n",
    "\n",
    "Draw(gfut, deformation=False, animate=True, autoscale = True , order = 3,  vectors = { \"grid_size\": 30}, settings=settings);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The eigenvalue problem on the torus:\n",
    "\n",
    "For simply connected domains the kernel of the curl is given by the image of the gradient. \n",
    "\n",
    "In case of a torus the kernel is slightly larger than the image of the gradient.\n",
    "\n",
    "**problem**: find $u \\in H(curl)$ such that\n",
    "\\begin{align*}  \n",
    "\\text{curl} \\,u &= 0 \\\\\n",
    "u &\\neq \\nabla \\phi \\quad \\forall \\phi \\in H^1\n",
    "\\end{align*}\n",
    "\n",
    "To do so we solve the eigenvalue problem:\\\n",
    "Find $\\lambda \\in \\mathbb{R}$ and $u \\in H(curl)$ such that\n",
    "\\begin{align*}\n",
    "\\int_{\\Omega} \\text{curl} u \\cdot \\text{curl} v \\,dx = \\lambda \\int_{\\Omega} u \\cdot v \\,dx\\\\\n",
    "u \\cdot \\nabla \\phi =0  \\quad \\forall \\phi \\in H^1\n",
    "\\end{align*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start creating the (double) torus "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "circ = Circle( Pnt(2,0), 1.5).Face()\n",
    "axis = Axis(Pnt(0,0,0), Dir(0,1,0))#, Dir(1,0,0))\n",
    "angle = 360\n",
    "torus = circ.Revolve( axis, angle)\n",
    "torus2 = torus.Move( (4,0,0), )\n",
    "\n",
    "geo = torus.Rotate(Axis((0, 0, 0), Y), 180) + torus2\n",
    "\n",
    "\n",
    "mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=3)).Curve(4)\n",
    "\n",
    "\n",
    "settings = {\"camera\": {\"transformations\": [ {\"type\": \"rotateX\", \"angle\": 45}]}}\n",
    "Draw(mesh, settings=settings);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I will not go too much into detail here, to have a better understanding of the problem I recommend to read the Phd thesis of Sabine Zaglmayr, p 145-150."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SetHeapSize(100*1000*1000)\n",
    "\n",
    "fes = HCurl(mesh, order=2)\n",
    "\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\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The idea is to use the PINVIT (preconditioned inverse iteration) method and let the preconditioner take care of the kernel of the curl.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager():\n",
    "    a.Assemble()\n",
    "    m.Assemble()\n",
    "    apre.Assemble() # when we assemble the BilinearForm , it assembles the matrix as well\n",
    "\n",
    "    ### Advanced part of the code ###\n",
    "    \n",
    "    gradmat, fesh1 = fes.CreateGradient() # create gradient from H1 to HCurl\n",
    "    gradmattrans = gradmat.CreateTranspose() \n",
    "    math1 = gradmattrans @ m.mat @ gradmat   # create projection onto grad H1\n",
    "    math1[0, 0] += 1     # fix the 1-dim kernel\n",
    "    invh1 = math1.Inverse(inverse=\"sparsecholesky\") \n",
    "    # build the Poisson projector with operator Algebra:\n",
    "    proj = IdentityMatrix() - gradmat @ invh1 @ gradmattrans @ m.mat \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": [
    "\n",
    "gfu = GridFunction(fes, multidim=len(evecs))\n",
    "for i in range(len(evecs)):\n",
    "    gfu.vecs[i].data = evecs[i]\n",
    "Draw(gfu, mesh, order=2,clipping = {\"y\":-1, \"z\":0} ,vectors = { \"grid_size\": 30},settings=settings);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ngs24",
   "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}