{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# 2.1.7 Multigrid for hybrid methods\n",
    "\n",
    "Mixed methods for second order problems can often be reduced to the mesh facet, so called hybrid mixed methods. Simiar, hybrid DG methods introduce new variables on the facets, such that the bulk of element variables can be condensed out.\n",
    "\n",
    "We show how to setup a multigrid preconditioner for hybrid methods. Interesting applications are nearly incompressible materials, or Stokes, discretized by $H(\\operatorname{div})$-conforming HDG or hybrid mixed methods.\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": [
    "The  hybrid DG method:\n",
    "\n",
    "<img src=\"facetelement.png\" alt=\"Alternative text\" width=\"100\" align=\"center\"/>\n",
    "\n",
    "On $V_h = V_T \\times V_F = P^k({\\mathcal T}) \\times P^k ({\\mathcal F})$ we define the bilinear form\n",
    "\n",
    "$$\n",
    "\\DeclareMathOperator{\\Div}{div}\n",
    "  \\sum_T \\int_T \\nabla u \\nabla v\n",
    "- \\sum_T \\int_{\\partial T} \\tfrac{\\partial u}{\\partial n} (v-\\widehat v)\n",
    "- \\sum_T \\int_{\\partial T} \\tfrac{\\partial v}{\\partial n} (u-\\widehat u)\n",
    "+ \\frac{\\alpha p^2}{h} \\sum_F \\int_F (u-\\widehat u)(v-\\widehat v)\n",
    "$$\n",
    "\n",
    "Element variables can be condensed out, which leads to a system reduced to the Skeleton. \n",
    "\n",
    "When splitting a large triangle $T_H$ into small trianles, the functions on $\\partial T_H$ have a canonical representation on the facets of the fine triangles. However, facet variables on internal edges of $T_H$ are not defined by embedding. The `HarmonicProlongation` provides the energy optimal extension to the internal edges. To define energy optimal we need the energy defined by a bilinear form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "ngmesh = unit_square.GenerateMesh(maxh=2)\n",
    "mesh = Mesh(ngmesh)\n",
    "\n",
    "order = 3\n",
    "fes = L2(mesh, order=order) * FacetFESpace(mesh, order=order, hoprolongation=True, dirichlet=\".*\")\n",
    "(u,uhat), (v,vhat) = fes.TnT()\n",
    "\n",
    "n = specialcf.normal(2)\n",
    "h = specialcf.mesh_size\n",
    "dS = dx(element_vb=BND)\n",
    "\n",
    "HDGform = u*v*dx+ grad(u)*grad(v)*dx - n*grad(u)*(v-vhat)*dS - n*grad(v)*(u-uhat)*dS + 5*(order+1)**2/h*(u-uhat)*(v-vhat)*dS\n",
    "bfa = BilinearForm(HDGform, condense=True).Assemble()\n",
    "fes.SetHarmonicProlongation(bfa, inverse=\"sparsecholesky\")\n",
    "pre = preconditioners.MultiGrid(bfa, blocktype=[\"vertexpatch\"], cycle=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager():\n",
    "    for l in range(7):\n",
    "        mesh.Refine()\n",
    "        bfa.Assemble()\n",
    "        # pre.Update()\n",
    "        lam = EigenValues_Preconditioner(bfa.mat, pre)\n",
    "        print (\"l =\", l, \"ndof =\", fes.ndof, \"lam_min/lam_max = \", lam[0], lam[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = LinearForm (x*v*dx).Assemble()\n",
    "gfu = GridFunction(fes)\n",
    "gfu.vec[:]=0\n",
    "with TaskManager():\n",
    "    Solve (bfa*gfu==f, pre, lin_solver=solvers.CGSolver, printrates=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (gfu.components[0]);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "## Hybrid-mixed methods:\n",
    "\n",
    "Find $\\sigma, u, \\widehat u \\in \\Sigma_h \\times V_h \\times F_h$:\n",
    "\n",
    "$$\n",
    "\\DeclareMathOperator{\\Div}{div}\n",
    "\\begin{array}{ccccccll}\n",
    "\\int a \\sigma \\tau & + & \\sum_T \\int_T \\Div \\tau \\, u & + & \\sum_F \\int_F [\\tau_n] \\widehat u & = & 0 & \\forall \\, \\tau \\in \\Sigma \\\\\n",
    "\\int \\Div \\sigma \\, v &&&&& = & \\int f v & \\forall \\, v \\in V_h \\\\\n",
    "\\int [ \\sigma_n ] \\, \\widehat v &&&&& = & \\int_{\\Gamma_n} g \\widehat v & \\forall \\, \\widehat v \\in F_h\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "where $\\Sigma_h$ is an discontinuous $H(div)$ finite element space, $V_h$ a sub-space of $L_2$, and $F_h$ consists of polynomials on every facet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "ngmesh = unit_square.GenerateMesh(maxh=0.2)\n",
    "mesh = Mesh(ngmesh)\n",
    "\n",
    "order = 2\n",
    "\n",
    "fesSigma = PrivateSpace(Discontinuous(HDiv(mesh, order=order, RT=True)))\n",
    "fesL2 = L2(mesh, order=order)\n",
    "fesFacet = FacetFESpace(mesh, order=order, hoprolongation=True, dirichlet=\".*\")\n",
    "fes = fesSigma*fesL2*fesFacet\n",
    "\n",
    "(sigma, u,uhat), (tau, v,vhat) = fes.TnT()\n",
    "n = specialcf.normal(2)\n",
    "dS = dx(element_vb=BND)\n",
    "mixedform = -sigma*tau*dx\n",
    "mixedform += div(sigma)*v*dx - sigma*n*vhat*dS \n",
    "mixedform += div(tau)*u*dx - tau*n*uhat*dS \n",
    "\n",
    "bfa = BilinearForm(mixedform, condense=True).Assemble()\n",
    "fes.SetHarmonicProlongation(bfa, inverse=\"sparsecholesky\")\n",
    "pre = preconditioners.MultiGrid(bfa, blocktype=[\"vertexpatch\"], cycle=1)\n",
    "\n",
    "for l in range(5):\n",
    "    mesh.Refine()\n",
    "    bfa.Assemble()\n",
    "    lam = EigenValues_Preconditioner(bfa.mat, pre)\n",
    "    print (\"l =\", l, \"ndof =\", fes.ndof, \"lam_min/lam_max = \", lam[0], lam[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = LinearForm(x*v*dx).Assemble()\n",
    "gfu = GridFunction(fes)\n",
    "Solve(bfa*gfu==f, pre, solvers.CGSolver, printrates=True)\n",
    "Draw (gfu.components[1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "## Nearly incompressible materials, H(div)-conforming HDG\n",
    "\n",
    "[Lehrenfeld+Schöberl, 2016]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "ngmesh = unit_square.GenerateMesh(maxh=0.3)\n",
    "mesh = Mesh(ngmesh)\n",
    "\n",
    "order = 3\n",
    "\n",
    "fesT = HDiv(mesh, order=order, hoprolongation=True, dirichlet=\".*\")\n",
    "fesF = TangentialFacetFESpace(mesh, order=order, hoprolongation=True, highest_order_dc=True, dirichlet=\".*\")\n",
    "fes = fesT*fesF\n",
    "\n",
    "(u,uhat), (v,vhat) = fes.TnT()\n",
    "n = specialcf.normal(2)\n",
    "def tang(v): return v-(v*n)*n\n",
    "h = specialcf.mesh_size\n",
    "dS = dx(element_vb=BND)\n",
    "\n",
    "HDGform = InnerProduct(Grad(u),Grad(v))*dx - (Grad(u)*n)*tang(v-vhat)*dS - (Grad(v)*n)*tang(u-uhat)*dS \\\n",
    "    + 1*(order+1)**2/h*tang(u-uhat)*tang(v-vhat)*dS\n",
    "\n",
    "bfa = BilinearForm(HDGform + 1e3*div(u)*div(v)*dx, condense=True).Assemble()\n",
    "fes.SetHarmonicProlongation(bfa)\n",
    "pre = preconditioners.MultiGrid(bfa, smoother=\"block\", smoothingsteps=1, blocktype=[\"vertexpatch\"], cycle=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager():\n",
    "    for l in range(4):\n",
    "        mesh.Refine()\n",
    "        bfa.Assemble()\n",
    "        lam = EigenValues_Preconditioner(bfa.mat, pre)\n",
    "        print (\"l =\", l, \"ndof =\", fes.ndof, lam[0], lam[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager():\n",
    "    f = LinearForm ((0.5-y)*v[0]*dx).Assemble()\n",
    "gfu = GridFunction(fes)\n",
    "gfu.vec[:]=0\n",
    "with TaskManager(pajetrace=10**8):\n",
    "    solvers.BVP(bfa, f, gfu, pre, print=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (gfu.components[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "## Nearly incompressible materials / Stokes in 3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "ngmesh = unit_cube.GenerateMesh(maxh=2)\n",
    "mesh = Mesh(ngmesh)\n",
    "\n",
    "order = 2\n",
    "\n",
    "fesT = HDiv(mesh, order=order, hoprolongation=True, dirichlet=\".*\")\n",
    "fesF = TangentialFacetFESpace(mesh, order=order, hoprolongation=True, highest_order_dc=True, dirichlet=\".*\")\n",
    "fes = fesT*fesF\n",
    "\n",
    "(u,uhat), (v,vhat) = fes.TnT()\n",
    "n = specialcf.normal(3)\n",
    "def tang(v): return v-(v*n)*n\n",
    "h = specialcf.mesh_size\n",
    "dS = dx(element_vb=BND)\n",
    "\n",
    "HDGform = 0.001*u*v*dx+InnerProduct(Grad(u),Grad(v))*dx - (Grad(u)*n)*tang(v-vhat)*dS - (Grad(v)*n)*tang(u-uhat)*dS \\\n",
    "    + 5*(order+1)**2/h*tang(u-uhat)*tang(v-vhat)*dS\n",
    "\n",
    "bfa = BilinearForm(HDGform + 1e3*div(u)*div(v)*dx, condense=True).Assemble()\n",
    "fes.SetHarmonicProlongation(bfa)\n",
    "pre = preconditioners.MultiGrid(bfa, smoother=\"block\", smoothingsteps=3, blocktype=[\"edgepatch\"], cycle=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager():\n",
    "    for l in range(3):\n",
    "        mesh.Refine()\n",
    "        bfa.Assemble()\n",
    "        # pre.Update()\n",
    "        lam = EigenValues_Preconditioner(bfa.mat, pre)\n",
    "        print (\"l =\", l, \"ndof =\", fes.ndof, \"lam min/max = \", lam[0], lam[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18",
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager():\n",
    "    f = LinearForm ((0.5-y)*v[0]*dx).Assemble()\n",
    "gfu = GridFunction(fes)\n",
    "gfu.vec[:]=0\n",
    "\n",
    "with TaskManager():\n",
    "    Solve(bfa*gfu==f, pre, solvers.CGSolver, printrates=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "metadata": {},
   "outputs": [],
   "source": [
    "clipping = { \"function\" : True,  \"pnt\" : (0.5,0.5,0.51), \"vec\" : (0,0,-1) }\n",
    "Draw (gfu.components[0], order=2, clipping=clipping);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20",
   "metadata": {},
   "source": [
    "## Flow channel in 3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from netgen.occ import *\n",
    "from ngsolve.webgui import Draw\n",
    "from ngsolve.krylovspace import CGSolver\n",
    "\n",
    "box = Box((0,0,0), (2.5,0.41,0.41))\n",
    "box.faces.name=\"wall\"\n",
    "box.faces.Min(X).name=\"inlet\"\n",
    "box.faces.Max(X).name=\"outlet\"\n",
    "cyl = Cylinder((0.5,0.2,0), Z, h=0.41,r=0.05)\n",
    "cyl.faces.name=\"cyl\"\n",
    "shape = box-cyl\n",
    "\n",
    "mesh = shape.GenerateMesh(maxh=0.2).Curve(3)\n",
    "Draw (mesh);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22",
   "metadata": {},
   "outputs": [],
   "source": [
    "order = 2\n",
    "\n",
    "fesT = HDiv(mesh, order=order, hoprolongation=True, dirichlet=\"wall|inlet|cyl\")\n",
    "fesF = TangentialFacetFESpace(mesh, order=order, hoprolongation=True, highest_order_dc=True, dirichlet=\".*\")\n",
    "fes = fesT*fesF\n",
    "\n",
    "(u,uhat), (v,vhat) = fes.TnT()\n",
    "n = specialcf.normal(3)\n",
    "def tang(v): return v-(v*n)*n\n",
    "h = specialcf.mesh_size\n",
    "dS = dx(element_vb=BND)\n",
    "\n",
    "HDGform = 0.001*u*v*dx+InnerProduct(Grad(u),Grad(v))*dx - (Grad(u)*n)*tang(v-vhat)*dS - (Grad(v)*n)*tang(u-uhat)*dS \\\n",
    "    + 5*(order+1)**2/h*tang(u-uhat)*tang(v-vhat)*dS\n",
    "\n",
    "with TaskManager():\n",
    "    bfa = BilinearForm(HDGform + 1e3*div(u)*div(v)*dx, condense=True).Assemble()\n",
    "    fes.SetHarmonicProlongation(bfa)\n",
    "    pre = preconditioners.MultiGrid(bfa, inverse=\"sparsecholesky\", smoother=\"block\", smoothingsteps=1, blocktype=[\"edgepatch\"], cycle=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23",
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager():\n",
    "    for l in range(1):\n",
    "        mesh.Refine()\n",
    "        bfa.Assemble()\n",
    "        lam = EigenValues_Preconditioner(bfa.mat, pre)\n",
    "        print (\"l =\", l, \"ndof =\", fes.ndof, \"lam min/max = \", lam[0], lam[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24",
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes)\n",
    "\n",
    "uin = (1.5*4*y*(0.41-y)/(0.41*0.41)*z*(0.41-z)/0.41**2,0, 0)\n",
    "gfu.components[0].Set(uin, definedon=mesh.Boundaries(\"inlet\"))\n",
    "\n",
    "inv = CGSolver(bfa.mat, pre.mat, printrates=True, tol=1e-5)\n",
    "\n",
    "with TaskManager():\n",
    "    gfu.vec.data -= inv@bfa.mat * gfu.vec\n",
    "gfu.vec.data += bfa.harmonic_extension * gfu.vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25",
   "metadata": {},
   "outputs": [],
   "source": [
    "clipping = { \"function\" : True,  \"pnt\" : (1,0.2,0.2 ), \"vec\" : (0,0,-1.0) }\n",
    "Draw (gfu.components[0], mesh, order=2, clipping=clipping);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26",
   "metadata": {},
   "source": [
    "Pressure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (div(gfu.components[0]), mesh, order=2, clipping=clipping, draw_surf=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30",
   "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
}