{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# 2.1.6. Multigrid for high order finite element spaces\n",
    "\n",
    "Officially released with NGSolve-2504, most high order spaces provide high order polynomial preserving prolongation operators. \n",
    "Current restrictions:\n",
    "* require simplicial meshes\n",
    "* refinement by bisection\n",
    "* require uniform polynomial order\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from ngsolve.la import EigenValues_Preconditioner"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "Many FESpaces provide now a high-order accurate prolongation operator. It has to be enabled by the flag `hoprolongation=True`, maybe becoming default in future."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))\n",
    "gfu = GridFunction(H1(mesh, order=10, hoprolongation=True))\n",
    "gfu.Set (sin(50*x*y))\n",
    "Draw (gfu, order=10);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "for l in range(2):\n",
    "    mesh.Refine()\n",
    "Draw (gfu, order=5);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "## Multigrid preconditioners for high order spaces\n",
    "\n",
    "If the high order prolongation is enabled, the multigrid preconditioner uses the high order discretization on the mesh hierarchy. If not, the coarse grid spaces use the lowest order spaces in the background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))\n",
    "\n",
    "fes = H1(mesh, order=5, hoprolongation=True, dirichlet=\".*\")\n",
    "u,v = fes.TnT()\n",
    "\n",
    "a = BilinearForm(grad(u)*grad(v)*dx).Assemble()\n",
    "pre = preconditioners.MultiGrid(a, blocktype=[\"vertexpatch\"])\n",
    "# pre = preconditioners.MultiGrid(a, smoother=\"block\", blocktype=[\"vertexedge\",\"face\"])\n",
    "# pre = preconditioners.MultiGrid(a, smoother=\"block\", blocktype=[\"vertpatch\"])\n",
    "\n",
    "for l in range(4):\n",
    "    mesh.Refine()\n",
    "    a.Assemble()\n",
    "    lam = EigenValues_Preconditioner(a.mat, pre)\n",
    "    print (mesh.levels, fes.ndof, lam[0], lam[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = LinearForm(x*v*dx).Assemble()\n",
    "gfu = GridFunction(fes)\n",
    "\n",
    "Solve(a * gfu == f, dirichlet=x*(1-x), lin_solver=solvers.CGSolver, pre=pre, printrates=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "ea = { \"euler_angles\" : (-70, 0,-55) }\n",
    "Draw (gfu, deformation=True, **ea);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9",
   "metadata": {},
   "source": [
    "## High order with static condensation\n",
    "\n",
    "For high order methods, static condensation may save a lot of computation. However, canonical prolongation for the skeleton variables only does not preserve high order polynomials. Here the `HarmonicProlongation` comes into play: It prolongates functions on the boundaries of the coarse grid elements, and then solves local Dirichlet problems for the dofs on the fine-grid skeleton inside the coarse grid elements.\n",
    "\n",
    "The Dirichlet problem is solved for the problem-specific bilinear form, which has to be provided when enabling the `HarmonicProlongation`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))\n",
    "\n",
    "fes = H1(mesh, order=5, hoprolongation=True, dirichlet=\".*\")\n",
    "u,v = fes.TnT()\n",
    "\n",
    "a = BilinearForm(grad(u)*grad(v)*dx, condense=True).Assemble()\n",
    "fes.SetHarmonicProlongation(a)\n",
    "pre = preconditioners.MultiGrid(a, blocktype=\"vertexpatch\")\n",
    "\n",
    "for l in range(5):\n",
    "    mesh.Refine()\n",
    "    a.Assemble()\n",
    "    lam = EigenValues_Preconditioner(a.mat, pre)\n",
    "    print (mesh.levels, fes.ndof, lam[0], lam[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = LinearForm(x*v*dx).Assemble()\n",
    "gfu = GridFunction(fes)\n",
    "\n",
    "Solve(a * gfu == f, dirichlet=x*(1-x), lin_solver=solvers.CGSolver, pre=pre, printrates=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (gfu, deformation=True, **ea);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "This example shows the result of a `HarmonicProlongation`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(unit_square.GenerateMesh(maxh=2))\n",
    "\n",
    "fes = H1(mesh, order=8, hoprolongation=True)\n",
    "u,v = fes.TnT()\n",
    "\n",
    "a = BilinearForm(grad(u)*grad(v)*dx, condense=True).Assemble()\n",
    "fes.SetHarmonicProlongation(a)\n",
    "pre = preconditioners.MultiGrid(a, smoother=\"block\", blocktype=\"vertexpatch\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes)\n",
    "gfu.Set(sin(10*x))\n",
    "Draw(gfu, order=5)\n",
    "mesh.Refine()\n",
    "Draw(gfu, order=5)\n",
    "gfu.vec.data += a.harmonic_extension * gfu.vec\n",
    "Draw (gfu, order=5);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16",
   "metadata": {},
   "source": [
    "##  Nearly incompressible elasticity and Stokes\n",
    "\n",
    "The Scott-Vogelius element:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))\n",
    "\n",
    "fes = VectorH1(mesh, order=4, hoprolongation=True, dirichlet=\".*\")\n",
    "u,v = fes.TnT()\n",
    "\n",
    "eq = InnerProduct(Grad(u),Grad(v)) * dx + 1e4 * div(u)*div(v)* dx\n",
    "a = BilinearForm(eq).Assemble()\n",
    "\n",
    "# a = BilinearForm( (Grad(u)|Grad(v)) * dx).Assemble()\n",
    "\n",
    "pre = preconditioners.MultiGrid(a, blocktype=[\"vertexpatch\"])\n",
    "\n",
    "with TaskManager():\n",
    "    for l in range(5):\n",
    "        mesh.Refine()\n",
    "        a.Assemble()\n",
    "        lam = EigenValues_Preconditioner(a.mat, pre)\n",
    "        print (mesh.levels, fes.ndof, lam[0], lam[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18",
   "metadata": {},
   "source": [
    "## A robust method with reduced integration\n",
    "\n",
    "We define the bilinear form as\n",
    "\n",
    "$$\n",
    "A(u,v) = \\int_\\Omega \\nabla u : \\nabla v \\, dx + \\frac{1}{\\varepsilon} \\int_\\Omega P_{L_2}^Q \\operatorname{div} u \\operatorname{div} v \\, dx,\n",
    "$$\n",
    "\n",
    "where $P_{L_2}^Q$ is an element-wise $L_2$-projector into a lower order space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))\n",
    "\n",
    "fesp = L2(mesh, order=0)\n",
    "fes = VectorH1(mesh, order=2, hoprolongation=True, dirichlet=\".*\")\n",
    "u,v = fes.TnT()\n",
    "\n",
    "eq = InnerProduct(Grad(u),Grad(v)) * dx + 1e4 * Interpolate(div(u),fesp)*div(v)* dx\n",
    "a = BilinearForm(eq).Assemble()\n",
    "\n",
    "# a = BilinearForm( (Grad(u)|Grad(v)) * dx).Assemble()\n",
    "fes.SetHarmonicProlongation(a)\n",
    "pre = preconditioners.MultiGrid(a, blocktype=[\"vertexpatch\"])\n",
    "\n",
    "for l in range(4):\n",
    "    mesh.Refine()\n",
    "    a.Assemble()\n",
    "    lam = EigenValues_Preconditioner(a.mat, pre)\n",
    "    print (mesh.levels, fes.ndof, lam[0], lam[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = LinearForm((x-0.5)*v[1]*dx).Assemble()\n",
    "gfu = GridFunction(fes)\n",
    "Solve (a*gfu==f, pre, solvers.CGSolver, printrates=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (gfu);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22",
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.occ import *\n",
    "shape = MoveTo(0,0).LineTo(1,0,\"in\").LineTo(1,1).LineTo(2,1).LineTo(3,0).LineTo(4,1).LineTo(5,1) \\\n",
    "    .LineTo(5,2,\"out\").LineTo(4,2).LineTo(3,1).LineTo(2,2).LineTo(1,2).Rotate(180).Arc(1,90).Close().Face()\n",
    "mesh = shape.GenerateMesh(dim=2, maxh=0.25).Curve(3)\n",
    "Draw (mesh)\n",
    "print (mesh.GetBoundaries())\n",
    "\n",
    "fesp = L2(mesh, order=0)\n",
    "fes = VectorH1(mesh, order=2, hoprolongation=True, dirichlet=\"in|default\")\n",
    "u,v = fes.TnT()\n",
    "\n",
    "eq = InnerProduct(Grad(u),Grad(v)) * dx + 1e4 * Interpolate(div(u),fesp)*div(v)* dx\n",
    "a = BilinearForm(eq).Assemble()\n",
    "\n",
    "# a = BilinearForm( (Grad(u)|Grad(v)) * dx).Assemble()\n",
    "fes.SetHarmonicProlongation(a)\n",
    "pre = preconditioners.MultiGrid(a, blocktype=[\"vertexpatch\"])\n",
    "\n",
    "for l in range(4):\n",
    "    mesh.Refine()\n",
    "    a.Assemble()\n",
    "    print (mesh.levels, fes.ndof)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23",
   "metadata": {},
   "source": [
    "from [Dissertation J. Schöberl](https://www.tuwien.at/index.php?eID=dumpFile&t=f&f=256729&token=a532668f99d52b812999d002e22655734632a80e),   page 116."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = LinearForm(fes)\n",
    "gfu = GridFunction(fes)\n",
    "Solve (a*gfu==f, pre, solvers.CGSolver, dirichlet=CF((0,x*(1-x))) | mesh.Boundaries(\"in\"), printrates=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (gfu);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}