{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.10 Dual basis functions\n",
    "We use dual basis functions to define interpolation operators, define transfer operators between different finite element spaces, and auxiliary space preconditioners."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Canonical interpolation\n",
    "The canonical finite element interpolation operator is defined by specifying the degrees of freedom. For low order methods these are typically nodal values, while for high order methods these are most often moments. For example, the interpolation of a function $u$ onto the $p^{th}$ order triangle is given by: find $u_{hp} \\in V_{hp}$ such that\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "u_{hp} (V) & = & u(V) \\quad \\forall \\text{ vertices } V \\\\\n",
    "\\int_E u_{hp} q & = & \\int_E u q \\quad \\forall q \\in P^{p-2}(E) \\; \\forall \\text{ edges } E \\\\\n",
    "\\int_T u_{hp} q & = & \\int_T u q \\quad \\forall q \\in P^{p-3}(T) \\; \\forall \\text{ triangles } T\n",
    "\\end{eqnarray*}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "\n",
    "import matplotlib.pylab as plt\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The NGSolve 'Set' function does local projection, and simple averaging. In particular, this does not respect point values in mesh vertices. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=3, low_order_space=False)  \n",
    "\n",
    "func = x*x*x*x\n",
    "gfu = GridFunction(fes)\n",
    "gfu.Set(func)\n",
    "Draw (gfu)\n",
    "print (gfu.vec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most NGSolve finite element spaces provide now a \"dual\" operator, which delivers the moments (i.e. the dual space basis functions) instead of function values. The integrals over faces, edges and also vertices are defined by co-dimension 1 (=BND), co-dimension 2 (=BBND) or co-dimension 3 (=BBBND) integrals over the volume elements.\n",
    "We define a variational problem for canonical interpolation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u,v = fes.TnT()\n",
    "vdual = v.Operator(\"dual\")\n",
    "\n",
    "a = BilinearForm(fes)\n",
    "a += u*vdual*dx + u*vdual*dx(element_vb=BND) + \\\n",
    "    u*vdual*dx(element_vb=BBND)\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(fes)\n",
    "f += func*vdual*dx + func*vdual*dx(element_vb=BND) + \\\n",
    "    func*vdual*dx(element_vb=BBND)\n",
    "f.Assemble()\n",
    "\n",
    "# interpolation in vertices preserves values 0 and 1\n",
    "gfu.vec.data = a.mat.Inverse() * f.vec\n",
    "print (gfu.vec)\n",
    "Draw (gfu);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The vertex degrees of freedom vanish for edge and element basis functions, and the edge degrees of freedom vanish for element basis functions, but not vice-versa. Thus, the obtained matrix A is block-triangular:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.sparse as sp\n",
    "A = sp.csr_matrix(a.mat.CSR())\n",
    "plt.rcParams['figure.figsize'] = (4,4)\n",
    "plt.spy(A)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use proper block Gauss-Seidel smoothing for solving with that block triangular matrix by blocking the dofs for the individual vertices, edges and elements. Since the NGSolve Gauss-Seidel smoother reorders the order of smoothing blocks for parallelization, we have to take care to first compute vertex values, then edge values, and finally element values by running three different Gauss-Seidel sweeps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vblocks = [fes.GetDofNrs(vertex) for vertex in mesh.vertices]\n",
    "eblocks = [fes.GetDofNrs(edge) for edge in mesh.edges]\n",
    "fblocks = [fes.GetDofNrs(face) for face in mesh.faces]\n",
    "\n",
    "print (vblocks)\n",
    "print (eblocks)\n",
    "print (fblocks)\n",
    "\n",
    "vinv = a.mat.CreateBlockSmoother(vblocks)\n",
    "einv = a.mat.CreateBlockSmoother(eblocks)\n",
    "finv = a.mat.CreateBlockSmoother(fblocks)\n",
    "\n",
    "vinv.Smooth(gfu.vec, f.vec)\n",
    "einv.Smooth(gfu.vec, f.vec)\n",
    "finv.Smooth(gfu.vec, f.vec)\n",
    "print (gfu.vec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Embedding Finite Element Spaces\n",
    "This interpolation can be used to transform functions from one finite element space $V_{src}$ to another one $V_{dst}$. We use the dual space of the destination space:\n",
    "\n",
    "$$\n",
    "\\int_{node} u_{dst} v_{dual} = \\int_{node} u_{src} v_{dual} \\qquad\n",
    "\\forall \\, v_{dual} \\; \\forall \\, \\text{nodes}\n",
    "$$\n",
    "\n",
    "The left hand side leads to a non-symmetric square matrix, the right hand side to a rectangular matrix.\n",
    "\n",
    "As an example we implement the transformation from an vector valued $H^1$ space into $H(\\operatorname{div})$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.geom2d import unit_square\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))\n",
    "\n",
    "fesh1 = VectorH1(mesh, order=2)\n",
    "feshdiv = HDiv(mesh, order=2)\n",
    "\n",
    "gfuh1 = GridFunction(fesh1)\n",
    "gfuh1.Set ( (x*x,y*y) )\n",
    "\n",
    "gfuhdiv = GridFunction(feshdiv, name=\"uhdiv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build the matrices, and use a direct solver:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amixed = BilinearForm(trialspace=fesh1, testspace=feshdiv)\n",
    "ahdiv = BilinearForm(feshdiv)\n",
    "\n",
    "u,v = feshdiv.TnT()\n",
    "vdual = v.Operator(\"dual\")\n",
    "uh1 = fesh1.TrialFunction()\n",
    "\n",
    "dS = dx(element_boundary=True)\n",
    "ahdiv += u*vdual * dx + u*vdual * dS\n",
    "ahdiv.Assemble()\n",
    "\n",
    "amixed += uh1*vdual*dx + uh1*vdual*dS\n",
    "amixed.Assemble()\n",
    "\n",
    "# build transformation operator:\n",
    "transform = ahdiv.mat.Inverse() @ amixed.mat\n",
    "gfuhdiv.vec.data = transform * gfuh1.vec\n",
    "\n",
    "Draw (gfuh1)\n",
    "Draw (gfuhdiv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We implement a linear operator performing the fast conversion by Gauss-Seidel smoothing: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MyBlockInverse(BaseMatrix):\n",
    "    def __init__ (self, mat, eblocks, fblocks):\n",
    "        super(MyBlockInverse, self).__init__()\n",
    "        self.mat = mat\n",
    "        self.einv = mat.CreateBlockSmoother(eblocks)\n",
    "        self.finv = mat.CreateBlockSmoother(fblocks)\n",
    "        self.res = self.mat.CreateColVector()\n",
    "\n",
    "    def CreateRowVector(self):\n",
    "        return self.mat.CreateColVector()\n",
    "    def CreateColVector(self):\n",
    "        return self.mat.CreateRowVector()\n",
    "        \n",
    "    def Mult(self, x, y):\n",
    "        # y[:] = 0\n",
    "        # self.einv.Smooth(y,x)    #   y = y +  A_E^-1  (x - A y)\n",
    "        # self.finv.Smooth(y,x)    #   y = y +  A_E^-1  (x - A y)\n",
    "        \n",
    "        # the same, but we see how to transpose that\n",
    "        y.data = self.einv * x\n",
    "        self.res.data = x - self.mat * y\n",
    "        y.data += finv * self.res\n",
    "\n",
    "    def MultTrans(self, x, y):\n",
    "        y.data = self.finv.T * x\n",
    "        self.res.data = x - self.mat.T * y\n",
    "        y.data += einv.T * self.res\n",
    "\n",
    "\n",
    "eblocks = [feshdiv.GetDofNrs(edge) for edge in mesh.edges]\n",
    "fblocks = [feshdiv.GetDofNrs(face) for face in mesh.faces]\n",
    "\n",
    "transform = MyBlockInverse(ahdiv.mat, eblocks, fblocks) @ amixed.mat\n",
    "gfuhdiv.vec.data = transform * gfuh1.vec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Auxiliary Space Preconditioning\n",
    "Nepomnyaschikh 91, Hiptmair-Xu 07, ....\n",
    "\n",
    "Assume we have a complicated problem with some complicated discretization, and we have good preconditioners for a nodal finite element discretization for the Laplace operator. By auxiliary space preconditioning we can reuse the simple preconditioners for the complicated problems. It is simple, and works well in many cases.\n",
    "\n",
    "As a simple example, we precondition a DG discretization by an $H^1$ conforming method. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The DG discretization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "order=4\n",
    "fesDG = L2(mesh, order=order, dgjumps=True)\n",
    "u,v = fesDG.TnT()\n",
    "aDG = BilinearForm(fesDG)\n",
    "jump_u = u-u.Other()\n",
    "jump_v = v-v.Other()\n",
    "n = specialcf.normal(2)\n",
    "mean_dudn = 0.5*n * (grad(u)+grad(u.Other()))\n",
    "mean_dvdn = 0.5*n * (grad(v)+grad(v.Other()))\n",
    "alpha = 4\n",
    "h = specialcf.mesh_size\n",
    "aDG = BilinearForm(fesDG)\n",
    "aDG += grad(u)*grad(v) * dx\n",
    "aDG += alpha*order**2/h*jump_u*jump_v * dx(skeleton=True)\n",
    "aDG += alpha*order**2/h*u*v * ds(skeleton=True)\n",
    "aDG += (-mean_dudn*jump_v -mean_dvdn*jump_u) * dx(skeleton=True)\n",
    "aDG += (-n*grad(u)*v-n*grad(v)*u)* ds(skeleton=True)\n",
    "aDG.Assemble()\n",
    "\n",
    "fDG = LinearForm(fesDG)\n",
    "fDG += 1*v * dx\n",
    "fDG.Assemble()\n",
    "gfuDG = GridFunction(fesDG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The auxiliary $H^1$ discretization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fesH1 = H1(mesh, order=2, dirichlet=\".*\")\n",
    "u,v = fesH1.TnT()\n",
    "aH1 = BilinearForm(fesH1)\n",
    "aH1 += grad(u)*grad(v)*dx\n",
    "preH1 = Preconditioner(aH1, \"bddc\")\n",
    "aH1.Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "transform = fesH1.ConvertL2Operator(fesDG)\n",
    "pre = transform @ preH1.mat @ transform.T + aDG.mat.CreateSmoother()\n",
    "\n",
    "solvers.CG(mat=aDG.mat, rhs=fDG.vec, sol=gfuDG.vec, pre=pre, \\\n",
    "           maxsteps=200)\n",
    "\n",
    "Draw (gfuDG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.10.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}