20210301 - Preconditioning High-Order Finite Elements with P-Multigrid and BDDC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# <center>Preconditioning High-Order Finite Elements with P-Multigrid and BDDC</center>\n",
    "\n",
    "### <center>Jeremy L Thompson</center>\n",
    "\n",
    "### <center>Jed Brown</center>\n",
    "\n",
    "Local Fourier Analysis Toolkit is a Julia package for analysis and tuning of multigrid and smoothing techniques for finite element operators with arbitrary, user provided weak forms of partial differential equations\n",
    "\n",
    "<br>\n",
    "<center>\n",
    "    <a href=\"\n",
    "         https://www.github.com/jeremylt/LFAToolkit.jl\">github.com/jeremylt/LFAToolkit.jl</a> \n",
    "</center>\n",
    "<style>\n",
    "//change background settings {}\n",
    "div.slide-background {\n",
    "\tborder-top: 30px midnightblue solid;\n",
    "\tborder-bottom: 30px midnightblue solid;\n",
    "}\n",
    "</style>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <center>Overview</center>\n",
    "\n",
    "* Introduction\n",
    "* Mathematical Background\n",
    " - Operator Representation\n",
    " - Local Fourier Analysis\n",
    "* PDE Examples\n",
    " - Poisson Problem\n",
    " - Neo-Hookean Hyperelasticity\n",
    "* Smoother Examples\n",
    " - Jacobi\n",
    " - Chebyshev\n",
    "* Multigrid Examples\n",
    " - p-multigrid\n",
    " - h-multigrid\n",
    "* Development Roadmap\n",
    "* References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <center>Introduction</center>\n",
    "\n",
    "<center>\n",
    "    <a href=\"\n",
    "         https://www.github.com/CEED/libCEED/\">github.com/CEED/libCEED/</a> \n",
    "</center>\n",
    "\n",
    "libCEED Neo-Hookean Hyperelasticity miniapp\n",
    "\n",
    "* Neo-Hookean hyperelasticity at finite strain\n",
    "\n",
    "* matrix-free operator evaluation\n",
    "\n",
    "* fully unstructured hexahedral meshes\n",
    "\n",
    "* runtime selection of basis polynomial order\n",
    "\n",
    "* p-multigrid preconditioner for Jacobian solve\n",
    "\n",
    "* algebraic multigrid coarse grid solve\n",
    "\n",
    "* Chebyshev smoother based on operator diagonal\n",
    "\n",
    "Development of BDDC smoother requires parameter tuning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Introduction</center>\n",
    "\n",
    "Local Fourier Analysis (LFA) of multigrid and preconditioning methods\n",
    "* models error reduction across frequencies\n",
    "* predicts convergence of methods\n",
    "* facilitates parameter tuning\n",
    "\n",
    "The technique was first introduced by Brandt in 1977<sup>[2]</sup> for h-multigrid with finite differences\n",
    "\n",
    "Technique was later adapted for finite elements and additional preconditioners"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "LFAToolkit.jl does LFA for\n",
    "* user defined weak forms of PDE\n",
    "* with user defined or arbitrary order polynomial bases\n",
    "* for h and p multigrid with various smoothers\n",
    "* and facilitates parameter optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <center>Mathematical Background</center>\n",
    "\n",
    "Mathematical representation of finite element operators and Local Fourier Analysis in LFAToolkit.jl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Mathematical Background</center>\n",
    "### <center>Operator Representation</center>\n",
    "\n",
    "LFAToolkit uses a finite element operator representation designed to facilitate matrix-free implementations <sup>[3]</sup>\n",
    "\n",
    "This representation supports arbitrary second order PDEs\n",
    "\n",
    "<br>\n",
    "\n",
    "$$\n",
    "\\int_{\\Omega} v \\cdot f_0 \\left( u, \\nabla u \\right) + \\nabla v : f_1 \\left( u, \\nabla u \\right) = \\int_{\\Omega} f v, \\forall v \\in V\n",
    "$$\n",
    "\n",
    "for some suitable $V \\subseteq H_0^1 \\left( \\Omega \\right)$\n",
    "\n",
    "where\n",
    "\n",
    "* $\\cdot$ represents contraction over fields\n",
    "* $:$ represents contraction over fields and spatial dimensions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Mathematical Background</center>\n",
    "### <center>Operator Representation</center>\n",
    "\n",
    "This weak form can be discretized by selecting a finite element basis\n",
    "\n",
    "<br>\n",
    "\n",
    "$$\n",
    "\\color{Maroon}{\\mathbf{A}} u = \\sum_e \\mathcal{E}^T_e \\left[\n",
    "  \\left( \\color{RoyalBlue}{\\mathbf{B}}_0^e \\right)^T \\color{RoyalBlue}{\\mathbf{W}}^e \\Lambda \\left( \\color{YellowGreen}{\\mathbf{f_0}} \\left( u^e, \\nabla u^e \\right) \\right)\n",
    "  + \\sum_{i = 1}^d \\left( \\color{RoyalBlue}{\\mathbf{B}}_i^e \\right)^T \\color{RoyalBlue}{\\mathbf{W}}^e \\Lambda \\left( \\color{YellowGreen}{\\mathbf{f_1}} \\left( u^e, \\nabla u^e \\right) \\right)\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "<br>\n",
    "\n",
    "<p style=\"margin-left:10%;\">$u^e = \\color{RoyalBlue}{\\mathbf{B}}_0^e \\mathcal{E}^e u$ and $\\nabla u^e = \\lbrace \\color{RoyalBlue}{\\mathbf{B}}_i^e \\mathcal{E}^e u \\rbrace_{i = 1}^d$</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Mathematical Background</center>\n",
    "### <center>Operator Representation</center>\n",
    "\n",
    "This operator can be represented as\n",
    "\n",
    "<br>\n",
    "\n",
    "$$\n",
    "\\color{Maroon}{\\mathbf{A}} = \\mathbf{G}^T \\color{RoyalBlue}{\\mathbf{B}}^T \\color{YellowGreen}{\\mathbf{D}} \\color{RoyalBlue}{\\mathbf{B}} \\mathbf{G}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "* $\\color{Maroon}{\\mathbf{A}}$ - full finite element operator\n",
    "* $\\mathbf{G}$ - element restriction or assembly operator\n",
    "* $\\color{RoyalBlue}{\\mathbf{B}}$ - element basis operator (providing $\\color{RoyalBlue}{\\mathbf{B}}_i$)\n",
    "* $\\color{YellowGreen}{\\mathbf{D}}$ - diagonal pointwise weak form operator (from $\\color{YellowGreen}{\\mathbf{f_0}}$, $\\color{YellowGreen}{\\mathbf{f_1}}$, and element geometry)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Mathematical Background</center>\n",
    "### <center>Operator Representation</center>\n",
    "\n",
    "This representation facilitates matrix-free implementation, such as in libCEED<sup>[1]</sup>\n",
    "<center>\n",
    "<a href=\"https://github.com/CEED/libCEED\">github.com/CEED/libCEED</a>\n",
    "<img src=\"https://raw.githubusercontent.com/CEED/libCEED/main/doc/img/libCEED.png\">\n",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Mathematical Background</center>\n",
    "### <center>Local Fourier Analysis</center>\n",
    "\n",
    "LFA investigates the local properties of the descretized operator\n",
    "\n",
    "Consider a scalar Toeplitz operator $L_h$ on the infinite 1D grid $G_h$\n",
    "\n",
    "<br>\n",
    "\n",
    "$$\n",
    "L_h \\mathrel{\\hat{=}} \\left[ s_\\kappa \\right]_h \\left( \\kappa \\in V \\right)\\\\\n",
    "L_h w_h \\left( x \\right) = \\sum_{\\kappa \\in V} s_\\kappa w_h \\left( x + \\kappa h \\right)\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "* $V \\subset \\mathcal{Z}$ is an index set\n",
    "* $s_\\kappa \\in \\mathcal{R}$ are constant coefficients\n",
    "* $w_h \\left( x \\right)$ is a $l^2$ function on $G_h$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Mathematical Background</center>\n",
    "### <center>Local Fourier Analysis</center>\n",
    "\n",
    "If for all $\\varphi \\left( \\theta, x \\right)$\n",
    "\n",
    "<br>\n",
    "\n",
    "$$\n",
    "L_h \\varphi \\left( \\theta, x \\right) = \\tilde{L}_h \\left( \\theta \\right) \\varphi \\left( \\theta, x \\right)\n",
    "$$\n",
    "\n",
    "then $\\tilde{L}_h \\left( \\theta \\right) = \\sum_{\\kappa \\in V} s_\\kappa e^{\\imath \\theta \\kappa}$ is the **symbol** of $L_h$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Mathematical Background</center>\n",
    "### <center>Local Fourier Analysis</center>\n",
    "\n",
    "For a system of equations, such as the $p \\times p$ linear system of equations for a scalar problem on a single $p$ order finite element\n",
    "\n",
    "<br>\n",
    "\n",
    "$$\n",
    "\\tilde{L}_h =\n",
    "\\begin{bmatrix}\n",
    "    \\tilde{L}_h^{1, 1}  &&  \\cdots  &&  \\tilde{L}_h^{1, p}  \\\\\n",
    "    \\vdots              &&  \\vdots  &&  \\vdots              \\\\\n",
    "    \\tilde{L}_h^{p, 1}  &&  \\cdots  &&  \\tilde{L}_h^{p, p}  \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Multi-component equations are represented in a similar fashion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Mathematical Background</center>\n",
    "### <center>Operator Representation + Local Fourier Analysis</center>\n",
    "\n",
    "We combine these two formulations with a Fourier mode localization operator in 1D\n",
    "\n",
    "<br>\n",
    "\n",
    "$$\n",
    "\\tilde{\\color{Maroon}{\\mathbf{A}}} \\left( \\mathbf{\\theta} \\right) = \\mathbf{Q}^T \\left( \\color{Maroon}{\\mathbf{A}}_e \\odot \\left[ e^{\\imath \\left( x_i - x_j \\right) \\theta / h} \\right] \\right) \\mathbf{Q}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "* $\\color{Maroon}{\\mathbf{A}}_e = \\color{RoyalBlue}{\\mathbf{B}}^T \\color{YellowGreen}{\\mathbf{D}} \\color{RoyalBlue}{\\mathbf{B}}$\n",
    "* $\\mathbf{Q} = \n",
    "\\begin{bmatrix}\n",
    "    I    \\\\\n",
    "    e_0  \\\\\n",
    "\\end{bmatrix} =\n",
    "\\begin{bmatrix}\n",
    "    1       &&  0       &&  \\cdots  &&  0       \\\\\n",
    "    0       &&  1       &&  \\cdots  &&  0       \\\\\n",
    "    \\vdots  &&  \\vdots  &&  \\vdots  &&  \\vdots  \\\\\n",
    "    0       &&  0       &&  \\cdots  &&  1       \\\\\n",
    "    1       &&  0       &&  \\cdots  &&  0       \\\\\n",
    "\\end{bmatrix}\n",
    "$\n",
    "\n",
    "Higher dimensions use with a Kronecker product of the 1D localization operator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <center>PDE Examples</center>\n",
    "\n",
    "PDEs can be represented in LFAToolkit.jl by providing bases, a mesh, and a pointwise representation of the weak form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>PDE Example</center>\n",
    "### <center>Poisson Problem</center>\n",
    "\n",
    "As the classic example, let's look at the scalar Poisson problem\n",
    "\n",
    "<br>\n",
    "\n",
    "$$\n",
    "\\int \\nabla v \\nabla u = \\int f v\n",
    "$$\n",
    "\n",
    "<br>\n",
    "\n",
    "$$\n",
    "\\color{Maroon}{\\mathbf{A}}_e = \\color{RoyalBlue}{\\mathbf{B}}^T \\color{YellowGreen}{\\mathbf{D}} \\color{RoyalBlue}{\\mathbf{B}}\n",
    "$$\n",
    "\n",
    "$\\color{RoyalBlue}{\\mathbf{B}}$ represents the finite element basis, such as a H1 Lagrange basis on Gauss-Lobatto points\n",
    "\n",
    "$\\color{YellowGreen}{\\mathbf{D}}$ represents the pointwise application of the weak form"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hideOutput": true,
    "hidePrompt": true,
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# setup\n",
    "using LFAToolkit\n",
    "\n",
    "# mesh\n",
    "dim = 1\n",
    "mesh = Mesh1D(1.0)\n",
    "# basis\n",
    "p = 3\n",
    "ncomp = 1\n",
    "basis = TensorH1LagrangeBasis(p, p+1, ncomp, dim)\n",
    "# weak form\n",
    "function diffusionweakform(du::Array{Float64}, w::Array{Float64})\n",
    "    return dv = du*w[1]\n",
    "end\n",
    "\n",
    "# operator\n",
    "inputs = [\n",
    "    OperatorField(basis, [EvaluationMode.gradient]),\n",
    "    OperatorField(basis, [EvaluationMode.quadratureweights]),\n",
    "]\n",
    "outputs = [OperatorField(basis, [EvaluationMode.gradient])]\n",
    "A = Operator(diffusionweakform, mesh, inputs, outputs)\n",
    "\n",
    "# compute symbols\n",
    "θ = [π]\n",
    "symbol = computesymbols(A, θ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     9,
     11,
     26,
     54
    ],
    "hideCode": true,
    "hidePrompt": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "using LFAToolkit\n",
    "using LinearAlgebra\n",
    "using Plots\n",
    "\n",
    "# setup\n",
    "p = 3\n",
    "dimension = 1\n",
    "\n",
    "mesh = []\n",
    "if dimension == 1\n",
    "   mesh = Mesh1D(1.0)\n",
    "elseif dimension == 2\n",
    "   mesh = Mesh2D(1.0, 1.0)\n",
    "end\n",
    "\n",
    "# operator\n",
    "diffusion = GalleryOperator(\"diffusion\", p, p, mesh)\n",
    "\n",
    "# compute full operator symbols\n",
    "numberruns = 250\n",
    "maxeigenvalue = 0\n",
    "θ_min = -π/2\n",
    "θ_max = 3π/2\n",
    "\n",
    "# compute and plot smoothing factor\n",
    "# -- 1D --\n",
    "if dimension == 1\n",
    "    # compute\n",
    "    maxeigenvalues = zeros(numberruns)\n",
    "\n",
    "    for i in 1:numberruns\n",
    "        θ = [θ_min + (θ_max - θ_min)*i/numberruns]\n",
    "        A = computesymbols(diffusion, θ)\n",
    "        eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "        maxeigenvalues[i] = max(eigenvalues...)\n",
    "        maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i])\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    plot(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        maxeigenvalues,\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        ylabel=\"spectral radius\",\n",
    "        linewidth=3,\n",
    "        legend=:none,\n",
    "        title=\"Diffusion Operator Symbol\"\n",
    "    )\n",
    "    ylims!(0.0, max(maxeigenvalues...) * 1.1)\n",
    "# -- 2D --\n",
    "elseif dimension == 2\n",
    "    # compute\n",
    "    maxeigenvalues = zeros(numberruns, numberruns)\n",
    "\n",
    "    for i in 1:numberruns, j in 1:numberruns\n",
    "        θ = [\n",
    "            θ_min + (θ_max - θ_min)*i/numberruns,\n",
    "            θ_min + (θ_max - θ_min)*j/numberruns\n",
    "        ]\n",
    "        A = computesymbols(diffusion, θ)\n",
    "        eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "        maxeigenvalues[i, j] = max(eigenvalues...)\n",
    "        maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i, j])\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    heatmap(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        xrange,\n",
    "        ylabel=\"θ/π\",\n",
    "        maxeigenvalues,\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        title=\"Diffusion Operator Symbol\",\n",
    "        transpose=true,\n",
    "        aspect_ratio=:equal\n",
    "    )\n",
    "    xlims!(θ_min/π, θ_max/π)\n",
    "    ylims!(θ_min/π, θ_max/π)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     9,
     11,
     26,
     54
    ],
    "hideCode": true,
    "hidePrompt": true,
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "using LFAToolkit\n",
    "using LinearAlgebra\n",
    "using Plots\n",
    "\n",
    "# setup\n",
    "p = 3\n",
    "dimension = 2\n",
    "\n",
    "mesh = []\n",
    "if dimension == 1\n",
    "   mesh = Mesh1D(1.0)\n",
    "elseif dimension == 2\n",
    "   mesh = Mesh2D(1.0, 1.0)\n",
    "end\n",
    "\n",
    "# operator\n",
    "diffusion = GalleryOperator(\"diffusion\", p, p, mesh)\n",
    "\n",
    "# compute full operator symbols\n",
    "numberruns = 250\n",
    "maxeigenvalue = 0\n",
    "θ_min = -π/2\n",
    "θ_max = 3π/2\n",
    "\n",
    "# compute and plot smoothing factor\n",
    "# -- 1D --\n",
    "if dimension == 1\n",
    "    # compute\n",
    "    maxeigenvalues = zeros(numberruns)\n",
    "\n",
    "    for i in 1:numberruns\n",
    "        θ = [θ_min + (θ_max - θ_min)*i/numberruns]\n",
    "        A = computesymbols(diffusion, θ)\n",
    "        eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "        maxeigenvalues[i] = max(eigenvalues...)\n",
    "        maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i])\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    plot(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        maxeigenvalues,\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        ylabel=\"spectral radius\",\n",
    "        linewidth=3,\n",
    "        legend=:none,\n",
    "        title=\"Diffusion Operator Symbol\"\n",
    "    )\n",
    "    ylims!(0.0, max(maxeigenvalues...) * 1.1)\n",
    "# -- 2D --\n",
    "elseif dimension == 2\n",
    "    # compute\n",
    "    maxeigenvalues = zeros(numberruns, numberruns)\n",
    "\n",
    "    for i in 1:numberruns, j in 1:numberruns\n",
    "        θ = [\n",
    "            θ_min + (θ_max - θ_min)*i/numberruns,\n",
    "            θ_min + (θ_max - θ_min)*j/numberruns\n",
    "        ]\n",
    "        A = computesymbols(diffusion, θ)\n",
    "        eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "        maxeigenvalues[i, j] = max(eigenvalues...)\n",
    "        maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i, j])\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    heatmap(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        xrange,\n",
    "        ylabel=\"θ/π\",\n",
    "        maxeigenvalues,\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        title=\"Diffusion Operator Symbol\",\n",
    "        transpose=true,\n",
    "        aspect_ratio=:equal\n",
    "    )\n",
    "    xlims!(θ_min/π, θ_max/π)\n",
    "    ylims!(θ_min/π, θ_max/π)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>PDE Example</center>\n",
    "### <center>Neo-Hookean Hyperelasticity</center>\n",
    "\n",
    "The strong form of the static balance of linear-momentum at finite strain is given by\n",
    "\n",
    "<br>\n",
    "\n",
    "$$\n",
    "-\\nabla_x \\cdot \\mathbf{P} - \\rho_0 \\mathbf{g} = \\mathbf{0}\n",
    "$$\n",
    "\n",
    "where $-\\nabla_k$ is the gradient with respect to the reference configuration, $\\mathbf{P}$ is the first Piola-Kirchhoff stress tensor, $\\rho_0$ is the reference mass density, and $\\mathbf{g}$ is the forcing function\n",
    "\n",
    "The first Piola-Kirchhoff stress tensor is given by\n",
    "\n",
    "<br>\n",
    "\n",
    "$$\n",
    "\\mathbf{P} = \\mathbf{F} \\mathbf{S}\n",
    "$$\n",
    "\n",
    "where $\\mathbf{F}$ is the deformation gradient and $\\mathbf{S}$ is the second Piola-Kirchhoff stress tensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "hideCode": false,
    "hideOutput": true,
    "hidePrompt": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# setup\n",
    "using LFAToolkit\n",
    "using LinearAlgebra\n",
    "\n",
    "# fem parameters\n",
    "dim = 3\n",
    "mesh = Mesh3D(1.0, 1.0, 1.0)\n",
    "p = 3\n",
    "ncomp = 3\n",
    "basis = TensorH1LagrangeBasis(p, p+1, ncomp, dim)\n",
    "\n",
    "# constants\n",
    "e = 1E6                     # Young's modulus\n",
    "ν = 0.3                     # Poisson's ratio\n",
    "K = e/(3*(1 - 2*ν))         # bulk modulus\n",
    "λ = e*ν/((1 + ν)*(1 - 2*ν)) # Lamé parameters\n",
    "μ = e/(2*(1 + ν))\n",
    "\n",
    "# state vector\n",
    "gradu = [1; 2; 3]*ones(1, 3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "hideCode": false,
    "hideOutput": true,
    "hidePrompt": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Neo-Hookean weak form\n",
    "function neohookeanweakform(deltadu::Array{Float64}, w::Array{Float64})\n",
    "    # dP = dF S + F dS\n",
    "\n",
    "    F = gradu + I\n",
    "    J = det(F)\n",
    "    # Green-Lagrange strain tensor\n",
    "    E = (gradu*gradu' + gradu'*gradu)/2\n",
    "    # right Cauchy-Green tensor\n",
    "    C = 2*E + I\n",
    "    C_inv = C^-1\n",
    "    # second Piola-Kirchhoff\n",
    "    S = λ*log(J)*C_inv + 2*μ*C_inv*E\n",
    "\n",
    "    deltadu = deltadu'\n",
    "    dF = deltadu + I\n",
    "    deltaE = (deltadu*deltadu' + deltadu'*deltadu)/2\n",
    "    dS = λ*sum(C_inv.*deltaE)*C_inv + 2*(μ - λ*log(J))*C_inv*deltaE*C_inv\n",
    "    dP = dF*S + F*dS\n",
    "\n",
    "    return [dP']\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "hideCode": false,
    "hideOutput": true,
    "hidePrompt": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# linearized Neo-Hookean operator\n",
    "inputs = [\n",
    "    OperatorField(\n",
    "        basis,\n",
    "        [EvaluationMode.gradient],\n",
    "        \"gradent of deformation\"\n",
    "        ),\n",
    "    OperatorField(\n",
    "        basis,\n",
    "        [EvaluationMode.quadratureweights],\n",
    "        ),\n",
    "]\n",
    "outputs = [\n",
    "    OperatorField(\n",
    "        basis,\n",
    "        [EvaluationMode.gradient],\n",
    "        \"test function gradient of deformation\",\n",
    "    ),\n",
    "]\n",
    "A = Operator(neohookeanweakform, mesh, inputs, outputs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <center>Smoother Examples</center>\n",
    "\n",
    "The error propagation operator for smoothers is given by\n",
    "\n",
    "<br>\n",
    "\n",
    "$$\n",
    "\\mathbf{S} = \\mathbf{I} - \\mathbf{M}^{-1} \\color{Maroon}{\\mathbf{A}}\n",
    "$$\n",
    "\n",
    "where $\\mathbf{M}^{-1}$ is the smoothing operator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Smoother Example</center>\n",
    "### <center>Jacobi</center>\n",
    "\n",
    "Let's consider simple Jacobi smoothing\n",
    "\n",
    "<br>\n",
    "\n",
    "$$\n",
    "\\mathbf{M}^{-1} = \\omega \\space \\text{diag} \\left( \\color{Maroon}{\\mathbf{A}} \\right)^{-1}\n",
    "$$\n",
    "\n",
    "The error propagation **symbol** is given by\n",
    "\n",
    "<br>\n",
    "\n",
    "$$\n",
    "\\tilde{\\mathbf{S}} \\left( \\omega, \\mathbf{\\theta} \\right) = \\mathbf{I} - \\omega \\space \\text{diag} \\left( \n",
    "\\tilde{\\color{Maroon}{\\mathbf{A}}} \\right) \\tilde{\\color{Maroon}{\\mathbf{A}}} \\left( \\mathbf{\\theta} \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hideOutput": true,
    "hidePrompt": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# setup\n",
    "using LFAToolkit\n",
    "\n",
    "# mesh\n",
    "dim = 2\n",
    "mesh = Mesh2D(1.0, 1.0)\n",
    "\n",
    "# operator\n",
    "A = GalleryOperator(\"diffusion\", p, p+1, mesh)\n",
    "S = Jacobi(A)\n",
    "\n",
    "# compute symbols\n",
    "ω = [1.0] # smoothing weight\n",
    "θ = [π, π]\n",
    "symbol = computesymbols(S, ω, θ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     9,
     11,
     29,
     61
    ],
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "using LFAToolkit\n",
    "using LinearAlgebra\n",
    "using Plots\n",
    "\n",
    "# setup\n",
    "p = 3\n",
    "dimension = 1\n",
    "\n",
    "mesh = []\n",
    "if dimension == 1\n",
    "   mesh = Mesh1D(1.0)\n",
    "elseif dimension == 2\n",
    "   mesh = Mesh2D(1.0, 1.0)\n",
    "end\n",
    "\n",
    "# operator\n",
    "diffusion = GalleryOperator(\"diffusion\", p, p, mesh)\n",
    "\n",
    "# Jacobi smoother\n",
    "jacobi = Jacobi(diffusion)\n",
    "\n",
    "# full operator symbols\n",
    "numberruns = 250\n",
    "maxeigenvalue = 0\n",
    "θ_min = -π/2\n",
    "θ_max = 3π/2\n",
    "\n",
    "# compute and plot smoothing factor\n",
    "# -- 1D --\n",
    "if dimension == 1\n",
    "    # setup\n",
    "    ω = [0.77]\n",
    "    maxeigenvalues = zeros(numberruns)\n",
    "\n",
    "    # compute\n",
    "    for i in 1:numberruns\n",
    "        θ = [θ_min + (θ_max - θ_min)*i/numberruns]\n",
    "        if abs(θ[1] % 2π) >  π/256\n",
    "            A = computesymbols(jacobi, ω, θ)\n",
    "            eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "            maxeigenvalues[i] = max(eigenvalues...)\n",
    "            maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    plot(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        maxeigenvalues,\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        ylabel=\"spectral radius\",\n",
    "        linewidth=3,\n",
    "        legend=:none,\n",
    "        title=\"Jacobi Operator Symbol\"\n",
    "    )\n",
    "    ylims!(0.0, max(maxeigenvalues...) * 1.1)\n",
    "# -- 2D --\n",
    "elseif dimension == 2\n",
    "    # setup\n",
    "    ω = [1.0]\n",
    "    maxeigenvalues = zeros(numberruns, numberruns)\n",
    "\n",
    "    # compute\n",
    "    for i in 1:numberruns, j in 1:numberruns\n",
    "        θ = [\n",
    "            θ_min + (θ_max - θ_min)*i/numberruns,\n",
    "            θ_min + (θ_max - θ_min)*j/numberruns\n",
    "        ]\n",
    "        if sqrt(abs(θ[1] % 2π)^2 + abs(θ[2] % 2π)^2) > π/128\n",
    "            A = computesymbols(jacobi, ω, θ)\n",
    "            eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "            maxeigenvalues[i, j] = max(eigenvalues...)\n",
    "            maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i, j])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    heatmap(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        xrange,\n",
    "        ylabel=\"θ/π\",\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        maxeigenvalues,\n",
    "        title=\"Jacobi Operator Symbol\",\n",
    "        transpose=true,\n",
    "        aspect_ratio=:equal\n",
    "    )\n",
    "    xlims!(θ_min/π, θ_max/π)\n",
    "    ylims!(θ_min/π, θ_max/π)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     9,
     11,
     29,
     61
    ],
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "using LFAToolkit\n",
    "using LinearAlgebra\n",
    "using Plots\n",
    "\n",
    "# setup\n",
    "p = 3\n",
    "dimension = 2\n",
    "\n",
    "mesh = []\n",
    "if dimension == 1\n",
    "   mesh = Mesh1D(1.0)\n",
    "elseif dimension == 2\n",
    "   mesh = Mesh2D(1.0, 1.0)\n",
    "end\n",
    "\n",
    "# operator\n",
    "diffusion = GalleryOperator(\"diffusion\", p, p, mesh)\n",
    "\n",
    "# Jacobi smoother\n",
    "jacobi = Jacobi(diffusion)\n",
    "\n",
    "# full operator symbols\n",
    "numberruns = 250\n",
    "maxeigenvalue = 0\n",
    "θ_min = -π/2\n",
    "θ_max = 3π/2\n",
    "\n",
    "# compute and plot smoothing factor\n",
    "# -- 1D --\n",
    "if dimension == 1\n",
    "    # setup\n",
    "    ω = [1.0]\n",
    "    maxeigenvalues = zeros(numberruns)\n",
    "\n",
    "    # compute\n",
    "    for i in 1:numberruns\n",
    "        θ = [θ_min + (θ_max - θ_min)*i/numberruns]\n",
    "        if abs(θ[1] % 2π) >  π/128\n",
    "            A = computesymbols(jacobi, ω, θ)\n",
    "            eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "            maxeigenvalues[i] = max(eigenvalues...)\n",
    "            maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    plot(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        maxeigenvalues,\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        ylabel=\"spectral radius\",\n",
    "        linewidth=3,\n",
    "        legend=:none,\n",
    "        title=\"Jacobi Operator Symbol\"\n",
    "    )\n",
    "    ylims!(0.0, max(maxeigenvalues...) * 1.1)\n",
    "# -- 2D --\n",
    "elseif dimension == 2\n",
    "    # setup\n",
    "    ω = [1.0]\n",
    "    maxeigenvalues = zeros(numberruns, numberruns)\n",
    "\n",
    "    # compute\n",
    "    for i in 1:numberruns, j in 1:numberruns\n",
    "        θ = [\n",
    "            θ_min + (θ_max - θ_min)*i/numberruns,\n",
    "            θ_min + (θ_max - θ_min)*j/numberruns\n",
    "        ]\n",
    "        if sqrt(abs(θ[1] % 2π)^2 + abs(θ[2] % 2π)^2) > π/128\n",
    "            A = computesymbols(jacobi, ω, θ)\n",
    "            eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "            maxeigenvalues[i, j] = max(eigenvalues...)\n",
    "            maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i, j])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    heatmap(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        xrange,\n",
    "        ylabel=\"θ/π\",\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        maxeigenvalues,\n",
    "        title=\"Jacobi Operator Symbol\",\n",
    "        transpose=true,\n",
    "        aspect_ratio=:equal\n",
    "    )\n",
    "    xlims!(θ_min/π, θ_max/π)\n",
    "    ylims!(θ_min/π, θ_max/π)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Smoother Example</center>\n",
    "### <center>Chebyshev</center>\n",
    "\n",
    "Let's consider Chebyshev polynomial smoothing\n",
    "\n",
    "<br>\n",
    "\n",
    "$$\n",
    "\\mathbf{M}^{-1} = \\rho_q \\left( \\text{diag} \\left( \\color{Maroon}{\\mathbf{A}} \\right)^{-1} \\color{Maroon}{\\mathbf{A}} \\right)\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$\\rho$ is derived from the Chebyshev polynomials<sup>[5]</sup> and $q$ is the polynomial order of the smoother\n",
    "\n",
    "The error propagation **symbol** is given by\n",
    "\n",
    "<br>\n",
    "\n",
    "$$\n",
    "\\tilde{\\mathbf{S}} \\left( q, \\mathbf{\\theta} \\right) = \\mathbf{I} - \\rho_q \\left( \\text{diag} \\left( \\tilde{\\color{Maroon}{\\mathbf{A}}} \\right)^{-1} \\tilde{\\color{Maroon}{\\mathbf{A}}} \\left( \\mathbf{\\theta} \\right) \\right) \\tilde{\\color{Maroon}{\\mathbf{A}}} \\left( \\mathbf{\\theta} \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hideOutput": true,
    "hidePrompt": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# setup\n",
    "using LFAToolkit\n",
    "\n",
    "# mesh\n",
    "dim = 2\n",
    "mesh = Mesh2D(1.0, 1.0)\n",
    "\n",
    "# operator\n",
    "A = GalleryOperator(\"diffusion\", p, p+1, mesh)\n",
    "S = Chebyshev(A)\n",
    "\n",
    "# compute symbols\n",
    "ρ = [3] # Chebyshev order\n",
    "θ = [π, π]\n",
    "symbol = computesymbols(S, ρ, θ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     9,
     11,
     29,
     61
    ],
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "using LFAToolkit\n",
    "using LinearAlgebra\n",
    "using Plots\n",
    "\n",
    "# setup\n",
    "p = 3\n",
    "dimension = 1\n",
    "\n",
    "mesh = []\n",
    "if dimension == 1\n",
    "   mesh = Mesh1D(1.0)\n",
    "elseif dimension == 2\n",
    "   mesh = Mesh2D(1.0, 1.0)\n",
    "end\n",
    "\n",
    "# operator\n",
    "diffusion = GalleryOperator(\"diffusion\", p, p, mesh)\n",
    "\n",
    "# Chebyshev smoother\n",
    "chebyshev = Chebyshev(diffusion)\n",
    "\n",
    "# full operator symbols\n",
    "numberruns = 250\n",
    "maxeigenvalue = 0\n",
    "θ_min = -π/2\n",
    "θ_max = 3π/2\n",
    "\n",
    "# compute and plot smoothing factor\n",
    "# -- 1D --\n",
    "if dimension == 1\n",
    "    # setup\n",
    "    ω = [3]\n",
    "    maxeigenvalues = zeros(numberruns)\n",
    "\n",
    "    # compute\n",
    "    for i in 1:numberruns\n",
    "        θ = [θ_min + (θ_max - θ_min)*i/numberruns]\n",
    "        if abs(θ[1] % 2π) >  π/256\n",
    "            A = computesymbols(chebyshev, ω, θ)\n",
    "            eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "            maxeigenvalues[i] = max(eigenvalues...)\n",
    "            maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    plot(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        maxeigenvalues,\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        ylabel=\"spectral radius\",\n",
    "        linewidth=3,\n",
    "        legend=:none,\n",
    "        title=\"Chebyshev Operator Symbol\"\n",
    "    )\n",
    "    ylims!(0.0, max(maxeigenvalues...) * 1.1)\n",
    "# -- 2D --\n",
    "elseif dimension == 2\n",
    "    # setup\n",
    "    ω = [3]\n",
    "    maxeigenvalues = zeros(numberruns, numberruns)\n",
    "\n",
    "    # compute\n",
    "    for i in 1:numberruns, j in 1:numberruns\n",
    "        θ = [\n",
    "            θ_min + (θ_max - θ_min)*i/numberruns,\n",
    "            θ_min + (θ_max - θ_min)*j/numberruns\n",
    "        ]\n",
    "        if sqrt(abs(θ[1] % 2π)^2 + abs(θ[2] % 2π)^2) > π/128\n",
    "            A = computesymbols(chebyshev, ω, θ)\n",
    "            eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "            maxeigenvalues[i, j] = max(eigenvalues...)\n",
    "            maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i, j])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    heatmap(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        xrange,\n",
    "        ylabel=\"θ/π\",\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        maxeigenvalues,\n",
    "        title=\"Chebyshev Operator Symbol\",\n",
    "        transpose=true,\n",
    "        aspect_ratio=:equal\n",
    "    )\n",
    "    xlims!(θ_min/π, θ_max/π)\n",
    "    ylims!(θ_min/π, θ_max/π)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     9,
     11,
     29,
     61
    ],
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "using LFAToolkit\n",
    "using LinearAlgebra\n",
    "using Plots\n",
    "\n",
    "# setup\n",
    "p = 3\n",
    "dimension = 2\n",
    "\n",
    "mesh = []\n",
    "if dimension == 1\n",
    "   mesh = Mesh1D(1.0)\n",
    "elseif dimension == 2\n",
    "   mesh = Mesh2D(1.0, 1.0)\n",
    "end\n",
    "\n",
    "# operator\n",
    "diffusion = GalleryOperator(\"diffusion\", p, p, mesh)\n",
    "\n",
    "# Chebyshev smoother\n",
    "chebyshev = Chebyshev(diffusion)\n",
    "\n",
    "# full operator symbols\n",
    "numberruns = 250\n",
    "maxeigenvalue = 0\n",
    "θ_min = -π/2\n",
    "θ_max = 3π/2\n",
    "\n",
    "# compute and plot smoothing factor\n",
    "# -- 1D --\n",
    "if dimension == 1\n",
    "    # setup\n",
    "    ω = [3]\n",
    "    maxeigenvalues = zeros(numberruns)\n",
    "\n",
    "    # compute\n",
    "    for i in 1:numberruns\n",
    "        θ = [θ_min + (θ_max - θ_min)*i/numberruns]\n",
    "        if abs(θ[1] % 2π) >  π/128\n",
    "            A = computesymbols(chebyshev, ω, θ)\n",
    "            eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "            maxeigenvalues[i] = max(eigenvalues...)\n",
    "            maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    plot(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        maxeigenvalues,\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        ylabel=\"spectral radius\",\n",
    "        linewidth=3,\n",
    "        legend=:none,\n",
    "        title=\"Chebyshev Operator Symbol\"\n",
    "    )\n",
    "    ylims!(0.0, max(maxeigenvalues...) * 1.1)\n",
    "# -- 2D --\n",
    "elseif dimension == 2\n",
    "    # setup\n",
    "    ω = [3]\n",
    "    maxeigenvalues = zeros(numberruns, numberruns)\n",
    "\n",
    "    # compute\n",
    "    for i in 1:numberruns, j in 1:numberruns\n",
    "        θ = [\n",
    "            θ_min + (θ_max - θ_min)*i/numberruns,\n",
    "            θ_min + (θ_max - θ_min)*j/numberruns\n",
    "        ]\n",
    "        if sqrt(abs(θ[1] % 2π)^2 + abs(θ[2] % 2π)^2) > π/128\n",
    "            A = computesymbols(chebyshev, ω, θ)\n",
    "            eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "            maxeigenvalues[i, j] = max(eigenvalues...)\n",
    "            maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i, j])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    heatmap(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        xrange,\n",
    "        ylabel=\"θ/π\",\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        maxeigenvalues,\n",
    "        title=\"Chebyshev Operator Symbol\",\n",
    "        transpose=true,\n",
    "        aspect_ratio=:equal\n",
    "    )\n",
    "    xlims!(θ_min/π, θ_max/π)\n",
    "    ylims!(θ_min/π, θ_max/π)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <center>Multigrid Examples</center>\n",
    "\n",
    "The error propagation of a multigrid method is given by\n",
    "\n",
    "<br>\n",
    "\n",
    "$$\n",
    "\\mathbf{M}_{TMG} = \\mathbf{S}_f \\left( \\mathbf{I} - \\mathbf{P}_{ctof} \\color{Maroon}{\\mathbf{A}}_c^{-1} \\mathbf{R}_{ftoc} \\color{Maroon}{\\mathbf{A}}_f \\right) \\mathbf{S}_f\n",
    "$$\n",
    "\n",
    "where the grid transfer opertors $\\mathbf{P}_{ctof}$ and $\\mathbf{R}_{ftoc}$ determine the type of multigrid method\n",
    "\n",
    "Error propagation symbol is given by the symbols of the components"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Multigrid Example</center>\n",
    "### <center>P-Multigrid</center>\n",
    "\n",
    "Prolongation given by interpolating low order basis functions to high order basis nodes, restriction by the transpose\n",
    "\n",
    "<br>\n",
    "\n",
    "$$\n",
    "\\mathbf{P}_{ctof} = \\mathbf{P}_f^T \\color{YellowGreen}{\\mathbf{D}}_{scale} \\color{RoyalBlue}{\\mathbf{B}}_{ctof} \\mathbf{P}_c, \\ \\ \\mathbf{R}_{ftoc} = \\mathbf{P}_{ctof}^T\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     4,
     15,
     27,
     42,
     53
    ],
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "using LFAToolkit\n",
    "using Plots\n",
    "\n",
    "# Lagrange plotting function\n",
    "function lagrange(x, i, p)\n",
    "    points = LFAToolkit.lobattoquadrature(p, false)\n",
    "    xi = points[i]\n",
    "    return prod([if (xi == xj) 1 else (x - xj)/(xi - xj) end for xj in points])\n",
    "end\n",
    "\n",
    "# coarse element\n",
    "x = -1:0.05:1\n",
    "pcoarse = 3\n",
    "coarse = [[lagrange(xx, i, pcoarse) for xx in x] for i in 1:pcoarse]\n",
    "fine = [[lagrange(xx, i, 4) + 1.25 for xx in x] for i in 1:4]\n",
    "plot(\n",
    "    x,\n",
    "    coarse,\n",
    "    yticks = false,\n",
    "    yshowaxis = false,\n",
    "    grid = false,\n",
    "    legend = false,\n",
    "    color = :blue,\n",
    "    size = (600, 300),\n",
    ")\n",
    "xi = LFAToolkit.lobattoquadrature(pcoarse, false)\n",
    "yi = zeros(pcoarse)\n",
    "scatter!(\n",
    "    xi,\n",
    "    yi,\n",
    "    yticks = false,\n",
    "    yshowaxis = false,\n",
    "    grid = false,\n",
    "    legend = false,\n",
    "    markersize = 5,\n",
    "    markercolor = :black,\n",
    ")\n",
    "\n",
    "# fine element\n",
    "shift = 1.25\n",
    "pfine = 5\n",
    "fine = [[lagrange(xx, i, pfine) + shift for xx in x] for i in 1:pfine]\n",
    "plot!(\n",
    "    x,\n",
    "    fine,\n",
    "    yticks = false,\n",
    "    yshowaxis = false,\n",
    "    grid = false,\n",
    "    legend = false,\n",
    "    color = :green\n",
    ")\n",
    "xi = LFAToolkit.lobattoquadrature(pfine, false)\n",
    "yi = ones(pfine)*shift\n",
    "scatter!(\n",
    "    xi,\n",
    "    yi,\n",
    "    yticks = false,\n",
    "    yshowaxis = false,\n",
    "    grid = false,\n",
    "    legend = false,\n",
    "    markersize = 5,\n",
    "    markercolor = :black,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hideOutput": true,
    "hidePrompt": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# setup\n",
    "using LFAToolkit\n",
    "\n",
    "# mesh\n",
    "dim = 2\n",
    "mesh = Mesh2D(1.0, 1.0)\n",
    "\n",
    "# setup\n",
    "finep = 3\n",
    "coarsep = 2\n",
    "ncomp = 1\n",
    "\n",
    "# diffusion operators\n",
    "finediffusion = GalleryOperator(\"diffusion\", finep, finep, mesh)\n",
    "coarsediffusion = GalleryOperator(\"diffusion\", coarsep, finep, mesh)\n",
    "\n",
    "# prolongation basis\n",
    "ctofbasis = TensorH1LagrangeBasis(coarsep, finep, ncomp, dim, lagrangequadrature = true)\n",
    "\n",
    "# Chebyshev smoother\n",
    "S = Chebyshev(finediffusion)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hideOutput": true,
    "hidePrompt": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# p-multigrid preconditioner\n",
    "M = PMultigrid(finediffusion, coarsediffusion, S, [ctofbasis])\n",
    "\n",
    "# compute symbols\n",
    "ρ = [3]\n",
    "v = [1, 1]\n",
    "θ = [π, π]\n",
    "symbol = computesymbols(M, ρ, v, θ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     11,
     13,
     37,
     70
    ],
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "using LFAToolkit\n",
    "using LinearAlgebra\n",
    "using Plots\n",
    "\n",
    "# setup\n",
    "finep = 3\n",
    "coarsep = 2\n",
    "numbercomponents = 1\n",
    "dimension = 1\n",
    "\n",
    "mesh = []\n",
    "if dimension == 1\n",
    "   mesh = Mesh1D(1.0)\n",
    "elseif dimension == 2\n",
    "   mesh = Mesh2D(1.0, 1.0)\n",
    "end\n",
    "\n",
    "ctofbasis = TensorH1LagrangePProlongationBasis(coarsep, finep, numbercomponents, dimension)\n",
    "\n",
    "# diffusion operators\n",
    "finediffusion = GalleryOperator(\"diffusion\", finep, finep, mesh)\n",
    "coarsediffusion = GalleryOperator(\"diffusion\", coarsep, finep, mesh)\n",
    "\n",
    "# Jacobi smoother\n",
    "chebyshev = Chebyshev(finediffusion)\n",
    "\n",
    "# p-multigrid preconditioner\n",
    "multigrid = PMultigrid(finediffusion, coarsediffusion, chebyshev, [ctofbasis])\n",
    "\n",
    "# full operator symbols\n",
    "numberruns = 250\n",
    "maxeigenvalue = 0\n",
    "θ_min = -π/2\n",
    "θ_max = 3π/2\n",
    "\n",
    "# compute and plot smoothing factor\n",
    "# -- 1D --\n",
    "if dimension == 1\n",
    "    # setup\n",
    "    ω = [3]\n",
    "    v = [1, 1]\n",
    "    maxeigenvalues = zeros(numberruns)\n",
    "\n",
    "    # compute\n",
    "    for i in 1:numberruns\n",
    "        θ = [θ_min + (θ_max - θ_min)*i/numberruns]\n",
    "        if abs(θ[1]  % 2π) >  π/256\n",
    "            A = computesymbols(multigrid, ω, v, θ)\n",
    "            eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "            maxeigenvalues[i] = max(eigenvalues...)\n",
    "            maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    plot(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        maxeigenvalues,\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        ylabel=\"spectral radius\",\n",
    "        linewidth=3,\n",
    "        legend=:none,\n",
    "        title=\"P-Multigrid Error Symbol\"\n",
    "    )\n",
    "    ylims!(0.0, max(maxeigenvalues...) * 1.1)\n",
    "# -- 2D --\n",
    "elseif dimension == 2\n",
    "    # setup\n",
    "    ω = [3]\n",
    "    v = [1, 1]\n",
    "    maxeigenvalues = zeros(numberruns, numberruns)\n",
    "\n",
    "    # compute\n",
    "    for i in 1:numberruns, j in 1:numberruns\n",
    "        θ = [\n",
    "            θ_min + (θ_max - θ_min)*i/numberruns,\n",
    "            θ_min + (θ_max - θ_min)*j/numberruns\n",
    "        ]\n",
    "        if sqrt(abs(θ[1] % 2π)^2 + abs(θ[2] % 2π)^2) > π/128\n",
    "            A = computesymbols(multigrid, ω, v, θ)\n",
    "            eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "            maxeigenvalues[i, j] = max(eigenvalues...)\n",
    "            maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i, j])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    heatmap(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        xrange,\n",
    "        ylabel=\"θ/π\",\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        maxeigenvalues,\n",
    "        title=\"P-Multigrid Error Symbol\",\n",
    "        transpose=true,\n",
    "        aspect_ratio=:equal\n",
    "    )\n",
    "    xlims!(θ_min/π, θ_max/π)\n",
    "    ylims!(θ_min/π, θ_max/π)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     11,
     13,
     37,
     70
    ],
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "using LFAToolkit\n",
    "using LinearAlgebra\n",
    "using Plots\n",
    "\n",
    "# setup\n",
    "finep = 3\n",
    "coarsep = 2\n",
    "numbercomponents = 1\n",
    "dimension = 2\n",
    "\n",
    "mesh = []\n",
    "if dimension == 1\n",
    "   mesh = Mesh1D(1.0)\n",
    "elseif dimension == 2\n",
    "   mesh = Mesh2D(1.0, 1.0)\n",
    "end\n",
    "\n",
    "ctofbasis = TensorH1LagrangePProlongationBasis(coarsep, finep, numbercomponents, dimension)\n",
    "\n",
    "# diffusion operators\n",
    "finediffusion = GalleryOperator(\"diffusion\", finep, finep, mesh)\n",
    "coarsediffusion = GalleryOperator(\"diffusion\", coarsep, finep, mesh)\n",
    "\n",
    "# Jacobi smoother\n",
    "chebyshev = Chebyshev(finediffusion)\n",
    "\n",
    "# p-multigrid preconditioner\n",
    "multigrid = PMultigrid(finediffusion, coarsediffusion, chebyshev, [ctofbasis])\n",
    "\n",
    "# full operator symbols\n",
    "numberruns = 250\n",
    "maxeigenvalue = 0\n",
    "θ_min = -π/2\n",
    "θ_max = 3π/2\n",
    "\n",
    "# compute and plot smoothing factor\n",
    "# -- 1D --\n",
    "if dimension == 1\n",
    "    # setup\n",
    "    ω = [3]\n",
    "    v = [1, 1]\n",
    "    maxeigenvalues = zeros(numberruns)\n",
    "\n",
    "    # compute\n",
    "    for i in 1:numberruns\n",
    "        θ = [θ_min + (θ_max - θ_min)*i/numberruns]\n",
    "        if abs(θ[1]  % 2π) >  π/256\n",
    "            A = computesymbols(multigrid, ω, v, θ)\n",
    "            eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "            maxeigenvalues[i] = max(eigenvalues...)\n",
    "            maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    plot(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        maxeigenvalues,\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        ylabel=\"spectral radius\",\n",
    "        linewidth=3,\n",
    "        legend=:none,\n",
    "        title=\"P-Multigrid Error Symbol\"\n",
    "    )\n",
    "    ylims!(0.0, max(maxeigenvalues...) * 1.1)\n",
    "# -- 2D --\n",
    "elseif dimension == 2\n",
    "    # setup\n",
    "    ω = [3]\n",
    "    v = [1, 1]\n",
    "    maxeigenvalues = zeros(numberruns, numberruns)\n",
    "\n",
    "    # compute\n",
    "    for i in 1:numberruns, j in 1:numberruns\n",
    "        θ = [\n",
    "            θ_min + (θ_max - θ_min)*i/numberruns,\n",
    "            θ_min + (θ_max - θ_min)*j/numberruns\n",
    "        ]\n",
    "        if sqrt(abs(θ[1] % 2π)^2 + abs(θ[2] % 2π)^2) > π/128\n",
    "            A = computesymbols(multigrid, ω, v, θ)\n",
    "            eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "            maxeigenvalues[i, j] = max(eigenvalues...)\n",
    "            maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i, j])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    heatmap(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        xrange,\n",
    "        ylabel=\"θ/π\",\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        maxeigenvalues,\n",
    "        title=\"P-Multigrid Error Symbol\",\n",
    "        transpose=true,\n",
    "        aspect_ratio=:equal\n",
    "    )\n",
    "    xlims!(θ_min/π, θ_max/π)\n",
    "    ylims!(θ_min/π, θ_max/π)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Multigrid Example</center>\n",
    "### <center>H-Multigrid</center>\n",
    "\n",
    "H-multigrid prolongation is interpolation from a single element to a *macro-element* with multiple sub *micro-elements*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     4,
     27,
     43
    ],
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "using LFAToolkit\n",
    "using Plots\n",
    "\n",
    "# Lagrange plotting function\n",
    "function lagrange(x, i, p)\n",
    "    points = LFAToolkit.lobattoquadrature(p, false)\n",
    "    xi = points[i]\n",
    "    return prod([if (xi == xj) 1 else (x - xj)/(xi - xj) end for xj in points])\n",
    "end\n",
    "\n",
    "# coarse element\n",
    "x = -1:0.05:1\n",
    "pcoarse = 2\n",
    "coarse = [[lagrange(xx, i, pcoarse) for xx in x] for i in 1:pcoarse]\n",
    "fine = [[lagrange(xx, i, 4) + 1.25 for xx in x] for i in 1:4]\n",
    "plot(\n",
    "    x,\n",
    "    coarse,\n",
    "    yticks = false,\n",
    "    yshowaxis = false,\n",
    "    grid = false,\n",
    "    legend = false,\n",
    "    color = :blue,\n",
    "    size = (600, 300),\n",
    ")\n",
    "xi = LFAToolkit.lobattoquadrature(pcoarse, false)\n",
    "yi = zeros(pcoarse)\n",
    "scatter!(\n",
    "    xi,\n",
    "    yi,\n",
    "    yticks = false,\n",
    "    yshowaxis = false,\n",
    "    grid = false,\n",
    "    legend = false,\n",
    "    markersize = 5,\n",
    "    markercolor = :black,\n",
    ")\n",
    "\n",
    "# left fine element\n",
    "shift = 1.25\n",
    "pfine = 2\n",
    "x = x ./ 2 .- 0.5\n",
    "fine = [[yi .+ shift] for yi in coarse]\n",
    "plot!(\n",
    "    x,\n",
    "    fine,\n",
    "    yticks = false,\n",
    "    yshowaxis = false,\n",
    "    grid = false,\n",
    "    legend = false,\n",
    "    color = :green\n",
    ")\n",
    "xi = LFAToolkit.lobattoquadrature(pfine, false) ./ 2 .- 0.5\n",
    "yi = ones(pfine)*shift\n",
    "scatter!(\n",
    "    xi,\n",
    "    yi,\n",
    "    yticks = false,\n",
    "    yshowaxis = false,\n",
    "    grid = false,\n",
    "    legend = false,\n",
    "    markersize = 5,\n",
    "    markercolor = :black,\n",
    ")\n",
    "\n",
    "# right fine element\n",
    "shift = 1.25\n",
    "x = x .+ 1\n",
    "fine = [[yi .+ shift] for yi in coarse]\n",
    "plot!(\n",
    "    x,\n",
    "    fine,\n",
    "    yticks = false,\n",
    "    yshowaxis = false,\n",
    "    grid = false,\n",
    "    legend = false,\n",
    "    color = :green\n",
    ")\n",
    "xi = LFAToolkit.lobattoquadrature(pfine, false) ./ 2 .+ 0.5\n",
    "yi = ones(pfine)*shift\n",
    "scatter!(\n",
    "    xi,\n",
    "    yi,\n",
    "    yticks = false,\n",
    "    yshowaxis = false,\n",
    "    grid = false,\n",
    "    legend = false,\n",
    "    markersize = 5,\n",
    "    markercolor = :black,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "hideCode": false,
    "hideOutput": true,
    "hidePrompt": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# setup\n",
    "using LFAToolkit\n",
    "\n",
    "# mesh\n",
    "dim = 2\n",
    "mesh = Mesh2D(1.0, 1.0)\n",
    "\n",
    "# bases \n",
    "ncomp = 1\n",
    "p = 2\n",
    "numelemfine1d = 2\n",
    "\n",
    "# weak form\n",
    "function diffusionweakform(du::Array{Float64}, w::Array{Float64})\n",
    "    dv = du*w[1]\n",
    "    return [dv]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     17
    ],
    "hideCode": false,
    "hideOutput": true,
    "hidePrompt": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# operators\n",
    "# -- fine\n",
    "basis = TensorH1LagrangeMacroBasis(p, p+1, ncomp, dim, numelemfine1d);\n",
    "inputs = [\n",
    "    OperatorField(basis, [EvaluationMode.gradient]),\n",
    "    OperatorField(basis, [EvaluationMode.quadratureweights]),\n",
    "];\n",
    "outputs = [OperatorField(basis, [EvaluationMode.gradient])];\n",
    "finediffusion = Operator(diffusionweakform, mesh, inputs, outputs);\n",
    "# -- coarse\n",
    "basis = TensorH1LagrangeBasis(p, p+1, ncomp, dim);\n",
    "inputs = [\n",
    "    OperatorField(basis, [EvaluationMode.gradient]),\n",
    "    OperatorField(basis, [EvaluationMode.quadratureweights]),\n",
    "];\n",
    "outputs = [OperatorField(basis, [EvaluationMode.gradient])];\n",
    "coarsediffusion = Operator(diffusionweakform, mesh, inputs, outputs);\n",
    "\n",
    "# prolongation basis\n",
    "ctofbasis = TensorH1LagrangeHProlongationBasis(p, ncomp, dim, numelemfine1d);\n",
    "\n",
    "# Chebyshev smoother\n",
    "S = Chebyshev(finediffusion)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [],
    "hideCode": false,
    "hideOutput": true,
    "hidePrompt": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# h-multigrid preconditioner\n",
    "M = HMultigrid(finediffusion, coarsediffusion, S, [ctofbasis])\n",
    "\n",
    "# compute symbols\n",
    "ρ = [3]\n",
    "v = [1, 1]\n",
    "θ = [π, π]\n",
    "symbol = computesymbols(M, ρ, v, θ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     10,
     12,
     20,
     26,
     34,
     55,
     88
    ],
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "using LFAToolkit\n",
    "using LinearAlgebra\n",
    "using Plots\n",
    "\n",
    "#setup\n",
    "p = 2\n",
    "numbercomponents = 1\n",
    "dimension = 1\n",
    "\n",
    "mesh = []\n",
    "if dimension == 1\n",
    "   mesh = Mesh1D(1.0)\n",
    "elseif dimension == 2\n",
    "   mesh = Mesh2D(1.0, 1.0)\n",
    "end\n",
    "\n",
    "numberfineelements1d = 2\n",
    "ctofbasis = TensorH1LagrangeHProlongationBasis(p, numbercomponents, dimension, numberfineelements1d);\n",
    "\n",
    "# operators\n",
    "function diffusionweakform(du::Array{Float64}, w::Array{Float64})\n",
    "    dv = du*w[1]\n",
    "    return [dv]\n",
    "end\n",
    "# -- fine\n",
    "basis = TensorH1LagrangeMacroBasis(p, p + 1, numbercomponents, dimension, numberfineelements1d);\n",
    "inputs = [\n",
    "    OperatorField(basis, [EvaluationMode.gradient]),\n",
    "    OperatorField(basis, [EvaluationMode.quadratureweights]),\n",
    "];\n",
    "outputs = [OperatorField(basis, [EvaluationMode.gradient])];\n",
    "finediffusion = Operator(diffusionweakform, mesh, inputs, outputs);\n",
    "# -- coarse\n",
    "basis = TensorH1LagrangeBasis(p, p + 1, numbercomponents, dimension);\n",
    "inputs = [\n",
    "    OperatorField(basis, [EvaluationMode.gradient]),\n",
    "    OperatorField(basis, [EvaluationMode.quadratureweights]),\n",
    "];\n",
    "outputs = [OperatorField(basis, [EvaluationMode.gradient])];\n",
    "coarsediffusion = Operator(diffusionweakform, mesh, inputs, outputs);\n",
    "\n",
    "# Chebyshev smoother\n",
    "chebyshev = Chebyshev(finediffusion)\n",
    "\n",
    "# h-multigrid preconditioner\n",
    "multigrid = HMultigrid(finediffusion, coarsediffusion, chebyshev, [ctofbasis])\n",
    "\n",
    "# full operator symbols\n",
    "numberruns = 250\n",
    "maxeigenvalue = 0\n",
    "θ_min = -1π/2\n",
    "θ_max = 3π/2\n",
    "\n",
    "# compute and plot smoothing factor\n",
    "# -- 1D --\n",
    "if dimension == 1\n",
    "    # setup\n",
    "    ω = [3]\n",
    "    v = [1, 1]\n",
    "    maxeigenvalues = zeros(numberruns)\n",
    "\n",
    "    # compute\n",
    "    for i in 1:numberruns\n",
    "        θ = [θ_min + (θ_max - θ_min)*i/numberruns]\n",
    "        if abs(θ[1]  % 2π) >  π/256\n",
    "            A = computesymbols(multigrid, ω, v, θ)\n",
    "            eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "            maxeigenvalues[i] = max(eigenvalues...)\n",
    "            maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    plot(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        maxeigenvalues,\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        ylabel=\"spectral radius\",\n",
    "        linewidth=3,\n",
    "        legend=:none,\n",
    "        title=\"H-Multigrid Error Symbol\"\n",
    "    )\n",
    "    ylims!(0.0, max(maxeigenvalues...) * 1.1)\n",
    "# -- 2D --\n",
    "elseif dimension == 2\n",
    "    # setup\n",
    "    ω = [3]\n",
    "    v = [1, 1]\n",
    "    maxeigenvalues = zeros(numberruns, numberruns)\n",
    "\n",
    "    # compute\n",
    "    for i in 1:numberruns, j in 1:numberruns\n",
    "        θ = [\n",
    "            θ_min + (θ_max - θ_min)*i/numberruns,\n",
    "            θ_min + (θ_max - θ_min)*j/numberruns\n",
    "        ]\n",
    "        if sqrt(abs(θ[1] % 2π)^2 + abs(θ[2] % 2π)^2) > π/128\n",
    "            A = computesymbols(multigrid, ω, v, θ)\n",
    "            eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "            maxeigenvalues[i, j] = max(eigenvalues...)\n",
    "            maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i, j])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    heatmap(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        xrange,\n",
    "        ylabel=\"θ/π\",\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        maxeigenvalues,\n",
    "        title=\"H-Multigrid Error Symbol\",\n",
    "        transpose=true,\n",
    "        aspect_ratio=:equal\n",
    "    )\n",
    "    xlims!(θ_min/π, θ_max/π)\n",
    "    ylims!(θ_min/π, θ_max/π)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     10,
     12,
     20,
     26,
     34,
     55,
     88
    ],
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "using LFAToolkit\n",
    "using LinearAlgebra\n",
    "using Plots\n",
    "\n",
    "#setup\n",
    "p = 2\n",
    "numbercomponents = 1\n",
    "dimension = 2\n",
    "\n",
    "mesh = []\n",
    "if dimension == 1\n",
    "   mesh = Mesh1D(1.0)\n",
    "elseif dimension == 2\n",
    "   mesh = Mesh2D(1.0, 1.0)\n",
    "end\n",
    "\n",
    "numberfineelements1d = 2\n",
    "ctofbasis = TensorH1LagrangeHProlongationBasis(p, numbercomponents, dimension, numberfineelements1d);\n",
    "\n",
    "# operators\n",
    "function diffusionweakform(du::Array{Float64}, w::Array{Float64})\n",
    "    dv = du*w[1]\n",
    "    return [dv]\n",
    "end\n",
    "# -- fine\n",
    "basis = TensorH1LagrangeMacroBasis(p, p + 1, numbercomponents, dimension, numberfineelements1d);\n",
    "inputs = [\n",
    "    OperatorField(basis, [EvaluationMode.gradient]),\n",
    "    OperatorField(basis, [EvaluationMode.quadratureweights]),\n",
    "];\n",
    "outputs = [OperatorField(basis, [EvaluationMode.gradient])];\n",
    "finediffusion = Operator(diffusionweakform, mesh, inputs, outputs);\n",
    "# -- coarse\n",
    "basis = TensorH1LagrangeBasis(p, p + 1, numbercomponents, dimension);\n",
    "inputs = [\n",
    "    OperatorField(basis, [EvaluationMode.gradient]),\n",
    "    OperatorField(basis, [EvaluationMode.quadratureweights]),\n",
    "];\n",
    "outputs = [OperatorField(basis, [EvaluationMode.gradient])];\n",
    "coarsediffusion = Operator(diffusionweakform, mesh, inputs, outputs);\n",
    "\n",
    "# Chebyshev smoother\n",
    "chebyshev = Chebyshev(finediffusion)\n",
    "\n",
    "# h-multigrid preconditioner\n",
    "multigrid = HMultigrid(finediffusion, coarsediffusion, chebyshev, [ctofbasis])\n",
    "\n",
    "# full operator symbols\n",
    "numberruns = 250\n",
    "maxeigenvalue = 0\n",
    "θ_min = -1π/2\n",
    "θ_max = 3π/2\n",
    "\n",
    "# compute and plot smoothing factor\n",
    "# -- 1D --\n",
    "if dimension == 1\n",
    "    # setup\n",
    "    ω = [3]\n",
    "    v = [1, 1]\n",
    "    maxeigenvalues = zeros(numberruns)\n",
    "\n",
    "    # compute\n",
    "    for i in 1:numberruns\n",
    "        θ = [θ_min + (θ_max - θ_min)*i/numberruns]\n",
    "        if abs(θ[1]  % 2π) >  π/256\n",
    "            A = computesymbols(multigrid, ω, v, θ)\n",
    "            eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "            maxeigenvalues[i] = max(eigenvalues...)\n",
    "            maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    plot(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        maxeigenvalues,\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        ylabel=\"spectral radius\",\n",
    "        linewidth=3,\n",
    "        legend=:none,\n",
    "        title=\"H-Multigrid Error Symbol\"\n",
    "    )\n",
    "    ylims!(0.0, max(maxeigenvalues...) * 1.1)\n",
    "# -- 2D --\n",
    "elseif dimension == 2\n",
    "    # setup\n",
    "    ω = [3]\n",
    "    v = [1, 1]\n",
    "    maxeigenvalues = zeros(numberruns, numberruns)\n",
    "\n",
    "    # compute\n",
    "    for i in 1:numberruns, j in 1:numberruns\n",
    "        θ = [\n",
    "            θ_min + (θ_max - θ_min)*i/numberruns,\n",
    "            θ_min + (θ_max - θ_min)*j/numberruns\n",
    "        ]\n",
    "        if sqrt(abs(θ[1] % 2π)^2 + abs(θ[2] % 2π)^2) > π/128\n",
    "            A = computesymbols(multigrid, ω, v, θ)\n",
    "            eigenvalues = [abs(val) for val in eigvals(A)]\n",
    "            maxeigenvalues[i, j] = max(eigenvalues...)\n",
    "            maxeigenvalue = max(maxeigenvalue, maxeigenvalues[i, j])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # plot\n",
    "    #println(\"max eigenvalue: \", maxeigenvalue)\n",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    heatmap(\n",
    "        xrange,\n",
    "        xlabel=\"θ/π\",\n",
    "        xtickfont=font(12, \"Courier\"),\n",
    "        xrange,\n",
    "        ylabel=\"θ/π\",\n",
    "        ytickfont=font(12, \"Courier\"),\n",
    "        maxeigenvalues,\n",
    "        title=\"H-Multigrid Error Symbol\",\n",
    "        transpose=true,\n",
    "        aspect_ratio=:equal\n",
    "    )\n",
    "    xlims!(θ_min/π, θ_max/π)\n",
    "    ylims!(θ_min/π, θ_max/π)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <center>Development Roadmap</center>\n",
    "\n",
    "<center>\n",
    "<a href=\"https://www.github.com/jeremylt\">github.com/jeremylt</a>\n",
    "<br>\n",
    "<a href=\"https://www.github.com/jeremylt/LFAToolkit.jl\">github.com/jeremylt/LFAToolkit.jl</a>\n",
    "</center>\n",
    "\n",
    "* v0.2 released with preconditioners in this talk\n",
    "\n",
    "* BDDC in current development<sup>[4]</sup>\n",
    "\n",
    "* Parameter optimization examples\n",
    "\n",
    "* Improved interoperability with libCEED.jl\n",
    "\n",
    "* Contributors and friendly users welcome"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <center>Questions?</center>\n",
    "\n",
    "<center>\n",
    "<a href=\"https://www.github.com/jeremylt\">github.com/jeremylt</a>\n",
    "<br>\n",
    "<a href=\"https://www.github.com/jeremylt/LFAToolkit.jl\">github.com/jeremylt/LFAToolkit.jl</a>\n",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <center>References</center>\n",
    "\n",
    "[1] A. Abdelfattah, V. Barra, N. Beams, J. Brown, J.S. Camier, V. Dobrev, Y. Dudouit,L. Ghaffari, T. Kolev, D. Medina, T. Rathnayake, J. L. Thompson,  and S. Tomov, *libCEED User Manual*, Zenodo (2020), https://doi.org/10.5281/zenodo.4302737.\n",
    "\n",
    "[2] A. Brandt, *Multi-level adaptive solutions to boundary-value problems*, Math. Comp., 31(138) (1977), pp. 33-390.\n",
    "\n",
    "[3] J. Brown, *Efficient nonlinear solvers for nodal high-order finite elements in 3D*, Journal of Scientific Computing, 45 (2010), pp. 48-63.\n",
    "\n",
    "[4] J. Brown, Y. He, S. MacLachlan, *Local Fourier Analysis of Balancing Domain Decomposition by Constraints Algorithms*, SIAM Journal on Scientific Computing, 41 (2019), pp. S346-S369.\n",
    "\n",
    "[5] M. Gutknecht and S. Röllin,  *The Chebyshev iteration revisited*, Parallel Computing, 28 (2002), pp. 263-283."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "hide_code_all_hidden": false,
  "kernelspec": {
   "display_name": "Julia 1.5.2",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}