20210402 - LocalFourierAnalysisOfPMultigridWithLFAToolkitjl.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# <center>Local Fourier Analysis of P-Multigrid<br>with LFAToolkit.jl</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",
    "* Agressive Coarsening\n",
    " - 1D Scalar Poisson\n",
    " - 2D Scalar Poisson\n",
    " - Neo-Hookean Hyperelasticity\n",
    "* Development Roadmap\n",
    "* References"
   ]
  },
  {
   "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>\n",
    "\n",
    "Note: improved libCEED.jl integration on development roadmap"
   ]
  },
  {
   "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": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Poisson Example</center>\n",
    "\n",
    "<center>1D quadratic elements</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     9,
     11,
     26,
     54
    ],
    "hideCode": true,
    "hidePrompt": true,
    "scrolled": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "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": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Poisson Example</center>\n",
    "\n",
    "<center>2D quadratic elements</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     9,
     11,
     26,
     54
    ],
    "hideCode": true,
    "hidePrompt": true,
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "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": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Jacobi Example</center>\n",
    "\n",
    "<center>1D quadratic elements, ω = 0.77</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     9,
     11,
     61
    ],
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "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": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Jacobi Example</center>\n",
    "\n",
    "<center>2D quadratic elements, ω = 1.0</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     9,
     11,
     29
    ],
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "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": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Chebyshev Example</center>\n",
    "\n",
    "<center>1D quadratic elements, ω = 3</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     9,
     11,
     29,
     61
    ],
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "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": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Chebyshev Example</center>\n",
    "\n",
    "<center>2D quadratic elements, ω = 3</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     9,
     11,
     29,
     61
    ],
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "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": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Multigrid Example</center>\n",
    "\n",
    "<center>1D quadratic to linear elements, 3rd order Chebyshev smoother</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     11,
     13,
     37,
     70
    ],
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "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": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Multigrid Example</center>\n",
    "\n",
    "<center>2D quadratic to linear elements, 3rd order Chebyshev smoother</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     11,
     13,
     37,
     70
    ],
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "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": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Multigrid Example</center>\n",
    "\n",
    "<center>1D linear elements, 3rd order Chebyshev smoother</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     10,
     12,
     20,
     26,
     34,
     55,
     88
    ],
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "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": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Multigrid Example</center>\n",
    "\n",
    "<center>2D linear elements, 3rd order Chebyshev smoother</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     10,
     12,
     20,
     26,
     34,
     55,
     88
    ],
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "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": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <center>Agressive Coarsening</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Agressive Coarsening</center>\n",
    "### <center>1D Scalar Diffusion - coarsen to linear</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "using LFAToolkit\n",
    "using LinearAlgebra\n",
    "using Plots\n",
    "\n",
    "# setup\n",
    "numberpairs = 4\n",
    "numbersmoothers = 3\n",
    "maxeigenvalues = zeros(numberpairs, numbersmoothers)\n",
    "\n",
    "for (finep, coarsep) in [(3, 2), (5, 2), (7, 2), (9, 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 = 100\n",
    "maxeigenvalue = 0\n",
    "θ_min = -π/2\n",
    "θ_max = 3π/2\n",
    "\n",
    "# compute effect of smoother\n",
    "for w in 3:5\n",
    "# -- 1D --\n",
    "if dimension == 1\n",
    "    # setup\n",
    "    ω = [w]\n",
    "    v = [1, 1]\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[floor(Int, (finep-1)/2), w-2] = max(maxeigenvalues[floor(Int, (finep-1)/2), w-2], eigenvalues...)\n",
    "        end\n",
    "    end\n",
    "# -- 2D --\n",
    "elseif dimension == 2\n",
    "    # setup\n",
    "    ω = [w]\n",
    "    v = [1, 1]\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[floor(Int, (finep-1)/2), w-2] = max(maxeigenvalues[floor(Int, (finep-1)/2), w-2], eigenvalues...)\n",
    "        end\n",
    "    end\n",
    "end\n",
    "end\n",
    "end\n",
    "\n",
    "# plot\n",
    "xrange = [2, 4, 6, 8]\n",
    "bar(\n",
    "    xrange,\n",
    "    xlabel=\"fine p\",\n",
    "    xtickfont=font(12, \"Courier\"),\n",
    "    maxeigenvalues[:, 1],\n",
    "    ytickfont=font(12, \"Courier\"),\n",
    "    ylabel=\"spectral radius\",\n",
    "    label=\"ω = 3\",\n",
    "    title=\"P-Multigrid with Chebyshev Smoothing Symbol\"\n",
    ")\n",
    "bar!(\n",
    "    xrange,\n",
    "    xlabel=\"fine p\",\n",
    "    xtickfont=font(12, \"Courier\"),\n",
    "    maxeigenvalues[:, 2],\n",
    "    ytickfont=font(12, \"Courier\"),\n",
    "    ylabel=\"spectral radius\",\n",
    "    label=\"ω = 4\",\n",
    "    title=\"P-Multigrid with Chebyshev Smoothing Symbol\"\n",
    ")\n",
    "bar!(\n",
    "    xrange,\n",
    "    xlabel=\"fine p\",\n",
    "    xtickfont=font(12, \"Courier\"),\n",
    "    maxeigenvalues[:, 3],\n",
    "    ytickfont=font(12, \"Courier\"),\n",
    "    ylabel=\"spectral radius\",\n",
    "    label=\"ω = 5\",\n",
    "    title=\"P-Multigrid with Chebyshev Smoothing Symbol\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Agressive Coarsening</center>\n",
    "### <center>1D Scalar Diffusion - Chebyshev Spectrum</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     10,
     12,
     61
    ],
    "hideCode": true,
    "hideOutput": false,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "using LFAToolkit\n",
    "using LinearAlgebra\n",
    "using Plots\n",
    "\n",
    "plot()\n",
    "\n",
    "# setup\n",
    "for finep in [3, 5, 7, 9]\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",
    "# diffusion operators\n",
    "finediffusion = GalleryOperator(\"diffusion\", finep, finep, mesh)\n",
    "\n",
    "# Jacobi smoother\n",
    "chebyshev = Chebyshev(finediffusion)\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",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    finepm1 = finep-1\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",
    "        label=\"p = \" * string(finepm1),\n",
    "        title=\"Chebyshev Smoother Error Symbol, ω = 3\"\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=\"P-Multigrid Error Symbol\",\n",
    "        transpose=true,\n",
    "        aspect_ratio=:equal\n",
    "    )\n",
    "    xlims!(θ_min/π, θ_max/π)\n",
    "    ylims!(θ_min/π, θ_max/π)\n",
    "end\n",
    "  \n",
    "end\n",
    "\n",
    "plot!(legend=:bottomleft)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Agressive Coarsening</center>\n",
    "### <center>1D Scalar Diffusion - Chebyshev Spectrum</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     12
    ],
    "hideCode": true,
    "hideOutput": false,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "using LFAToolkit\n",
    "using LinearAlgebra\n",
    "using Plots\n",
    "\n",
    "plot()\n",
    "\n",
    "# setup\n",
    "for finep in [3, 5, 7, 9]\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",
    "# diffusion operators\n",
    "finediffusion = GalleryOperator(\"diffusion\", finep, finep, mesh)\n",
    "\n",
    "# Jacobi smoother\n",
    "chebyshev = Chebyshev(finediffusion)\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",
    "    ω = [4]\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",
    "    xrange = θ_min/π:(θ_max - θ_min)/π/(numberruns-1):θ_max/π\n",
    "    finepm1 = finep-1\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",
    "        label=\"p = \" * string(finepm1),\n",
    "        title=\"Chebyshev Smoother Error Symbol, ω = 4\"\n",
    "    )\n",
    "    ylims!(0.0, max(maxeigenvalues...) * 1.1)\n",
    "# -- 2D --\n",
    "elseif dimension == 2\n",
    "    # setup\n",
    "    ω = [4]\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=\"P-Multigrid Error Symbol\",\n",
    "        transpose=true,\n",
    "        aspect_ratio=:equal\n",
    "    )\n",
    "    xlims!(θ_min/π, θ_max/π)\n",
    "    ylims!(θ_min/π, θ_max/π)\n",
    "end\n",
    "  \n",
    "end\n",
    "\n",
    "plot!(legend=:bottomleft)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Agressive Coarsening</center>\n",
    "### <center>2D Scalar Diffusion - coarsen to linear</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "using LFAToolkit\n",
    "using LinearAlgebra\n",
    "using Plots\n",
    "\n",
    "# setup\n",
    "numberpairs = 4\n",
    "numbersmoothers = 3\n",
    "maxeigenvalues = zeros(numberpairs, numbersmoothers)\n",
    "\n",
    "for (finep, coarsep) in [(3, 2), (5, 2), (7, 2), (9, 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 = 100\n",
    "maxeigenvalue = 0\n",
    "θ_min = -π/2\n",
    "θ_max = 3π/2\n",
    "\n",
    "# compute effect of smoother\n",
    "for w in 3:5\n",
    "# -- 1D --\n",
    "if dimension == 1\n",
    "    # setup\n",
    "    ω = [w]\n",
    "    v = [1, 1]\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[floor(Int, (finep-1)/2), w-2] = max(maxeigenvalues[floor(Int, (finep-1)/2), w-2], eigenvalues...)\n",
    "        end\n",
    "    end\n",
    "# -- 2D --\n",
    "elseif dimension == 2\n",
    "    # setup\n",
    "    ω = [w]\n",
    "    v = [1, 1]\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[floor(Int, (finep-1)/2), w-2] = max(maxeigenvalues[floor(Int, (finep-1)/2), w-2], eigenvalues...)\n",
    "        end\n",
    "    end\n",
    "end\n",
    "end\n",
    "end\n",
    "\n",
    "# plot\n",
    "xrange = [2, 4, 6, 8]\n",
    "bar(\n",
    "    xrange,\n",
    "    xlabel=\"fine p\",\n",
    "    xtickfont=font(12, \"Courier\"),\n",
    "    maxeigenvalues[:, 1],\n",
    "    ytickfont=font(12, \"Courier\"),\n",
    "    ylabel=\"spectral radius\",\n",
    "    label=\"ω = 3\",\n",
    "    title=\"P-Multigrid with Chebyshev Smoothing Symbol\"\n",
    ")\n",
    "bar!(\n",
    "    xrange,\n",
    "    xlabel=\"fine p\",\n",
    "    xtickfont=font(12, \"Courier\"),\n",
    "    maxeigenvalues[:, 2],\n",
    "    ytickfont=font(12, \"Courier\"),\n",
    "    ylabel=\"spectral radius\",\n",
    "    label=\"ω = 4\",\n",
    "    title=\"P-Multigrid with Chebyshev Smoothing Symbol\"\n",
    ")\n",
    "bar!(\n",
    "    xrange,\n",
    "    xlabel=\"fine p\",\n",
    "    xtickfont=font(12, \"Courier\"),\n",
    "    maxeigenvalues[:, 3],\n",
    "    ytickfont=font(12, \"Courier\"),\n",
    "    ylabel=\"spectral radius\",\n",
    "    label=\"ω = 5\",\n",
    "    title=\"P-Multigrid with Chebyshev Smoothing Symbol\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Agressive Coarsening</center>\n",
    "### <center>Neo-Hookean Hyperelasticity</center>"
   ]
  },
  {
   "attachments": {
    "Screenshot_2021-04-01%2020210402%20-%20LocalFourierAnalysisOfPMultigridwithLFAToolkitjl-Copy1%20-%20Jupyter%20Notebook.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAYAAAByNR6YAAAgAElEQVR4nO2df5wT9YH3h0eB5Yeb4A+gwpmAKygKu5RDaFGyilUEgUW85yxVdkHlPGi72yuCSK/ZO22FJbIBtCrXEuSsrFXIop6eWpulqMu2SiwtAgUMWGWL/JiFXdhFWD7PHzzfYZKdZJOQbCbz/bxfr/frBUl2ZjKZmbwzM5koIIQQQgghKUXJ9AQQQgghhFgNBhYhhBBCSIphYBHTsGfPHrjdbtTW1sb1+F/84hdwu91pmZYDBw7A7XYjEAjE9fjy8nKUl5enZVr0/PGPf4Tb7cbnn38e1+PLy8vh8/nSMi01NTWaRvh8PpSXlyMUCqV83OkcdryEQqG0zt90EAqFtNcsk/Mu3SSyfOzfvx9utxu///3v0z9hcRIMBrXXSVXVTE9OGOlc5t9//3243W58+eWXaRl+R8PAsiiBQACKorTRbrejqKgIa9asSXg4Tqcz5mOnTJmiPTaZFfDtt9+GoiiorKzUxl1eXh41coYNGwZFUXDmzJmEx9UeW7duhaIocUeTeN7p5tlnn4WiKPjggw8AnJ9H0d5IFEWBy+VK2fiDwSCKiooMl63CwsKw5crlckFRlLgjNRHSOex4EetGKudvuvB6vXA6nYavW1FRUdbFlojbaK9/IsvH+++/D0VR8POf/zy1E5kgqqqipKQEdru9zWvkdDpRUlKS0ekTpHOZr6iogKIo+Oijj9Iy/I6GgWVR9GHkcrk08/PztdvjWWEjQ83v9xs+LhQKhT0uFYHldruhKErUvVSzZs1K24r+17/+FS6XK+4Q7ajA2rhxI1wuF/785z8DOD+Por2RpHJj6PP5tOeZn5+P0tJSuN1uuN3usLgWeyAZWOagtLRUe22Ki4u116y0tFTbHmRyPiaDmPfRtg2JLB/btm2Dy+XCiy++mNqJTABVVVFQUABFUeBwOMLWreLiYjgcjg7ZvsQDAyt+zPGKkZSjD6NI9G+Uu3fvjms4YgUvLi42fFxlZaX2xttRgWUmOiqwIumowPL7/dpzFK9PJKqqYsqUKdq0MLAyj/jgY7PZEAwGoz4m2/ZgpTKwzIB++xntkGC0D7cdDQMrfhhYFiVWYAHnN0DthZD+TUTspTDaADidTjgcDu0NXwxXnEsQbaNRU1Ojbfj1gRUMBlFSUqLtadOf73P69GkA585HMtqAVldXo7y8HF6vV3vjEOeeNDc3h02Xfjr05xY0NjYiEAhg3759YcNWVRVr1qxBeXk5qqurtecVT2CJaTCaF9HOt1BVFTU1NaivrwcAfPnllwgEAmhoaEBNTY02j7xer+E5UfqNoaqq2rxZs2ZNQud2iMNL0eJKj5jHkW9ywWBQe12ivdkD5+aT1+vVptPozT9y2OL107/mgsjX2gijx4RCIW1+iddbP+zIwIq2bBgROVw9Bw8eRCAQwM6dO6P+fX19PQKBAP72t7/FfF4ijBN5Q4ycF2JajZ5TMBiE1+uN+joJxOtj9HwjCYVC2nyM9np6vV7DbYOYvkSWvYaGBgQCAezfv18bv/4cNTG+8vLymMttPNudaMS7PTaavliPEduNyG2Pfn5EDkfMf7FNiSTWMh8L/XyMtrwwsEhW0F5gRYZQe8NxuVzanq/IvwkGg1AURdutrX9Me58k9SurPrDE3xkpNhKR52Dpd7PrLSsr06ZLrNRi+MFgMOxvxLQYnYMVCATanB9ht9vbndcCMf8iI0X/95GfUsUn21deeQVA+DlY0eaPfjrEcwoGg4bTHusNQyDepB0OR7uP1SPmsd/vR2FhYZtpjNz7oKpq1PO7ysrKog472msuKC4ujrmsi708+udXVlYWdd6KZTly3Yicv06ns8389fv9hudC6R/7xRdfQFEUDBgwIOq8feCBB6AoCt59992Yr0EygaVfNyKnVSwz0V4ro21De89XT7T5XlRU1CaeYr027S17+uUj8hws/Z5z8QEmclr0JLLdiYb48BrPBxixvS0sLGx3eJs3bw57Tn6/P+brJsJVb+SpJPplPp7XNdZ6HXmOKwOLZAXp2IMFADabDQUFBWGPEW9goVAoZYEVDAa14RYXFyMQCGiKPViRgSU2pC6XC4FAAKFQCH6/HzabTdvIRwZWQUEB8vPz4ff7teEDbQMrFAppb6Bi+gKBAIqLi8PeWGMh3sinTJkSdruYZ4pyLlL1iA2lOJSrDywxfjFN+nmkn79OpxN2ux2lpaXa/eK8nPa+uACcP4cn2uHhaIh5bLfbw+axz+eDzWZr88Yj3qTE6xcIBOD3+7XDzvogixx2IBBAMBjUXm/9m1V7b0hi/ovHiygR0ywOoYn5Jt5AxLohXn/xrdNAIKC9bvpxisfbbDbt9QoGg9r47Xa7FhFiWf7www/bTG9zczN69OiBnj174tSpUzFfA/Hcjd7MoiHmrdPp1A756pe1wsJClJSUhL2mlZWV2nzX7zUR80a/DovlKTLwxXyw2Wzw+/3a6ylefxE2wWBQ++ARuW2IjLBYy55YT6IFlt1uh8Ph0P7W7/drp0rot5tiuRXzKhQKaeOJ3O5EQzwfp9MZ12FNMU+Mhmv0gUE8p8jtnRivWP71y6Z+3uufr36Zj9wWiueg39MplmUx3mAwGPY66NdrBhbJCmIFlv58mvb2YEQGlj6mBGIjBrTdM5ZsYOmHFe08C31g6d8UY82LyMCK9sk+MrDE8zaaFv0J3u1hdLJqQUEBHA4H8vPz28Sr3W6HzWbT/h/5LcJ4zsGKNt1iHrS38RePS/RcOPF3Rq+J+PQrglL8PzI+gXOfgG02G+x2e5vDjzabrc1hK334CGK9IYk3QTEcEQDtvdHplyujDyritRbTLMZjtM5FLutifvzLv/xLm8e++OKLUBQFs2bNijl9Av1J7iKQysrKoh7SEfPWKKjFczKa7+I5iL2tYp2J/NCgf6x4vVVVjbpNEq+//jUR8769c7DiWfaiBVasZUtsN8SwjLYjRtudaKiqGvYFJPFt72iH0yOfg57IDwz62xwOR5vnpF8+Ipd5Eej65xdrmyK2heI+MQ+MxiuGbbfb0djYCICBRbIE/cotzn0oKSkJ25Udzx6JyA2K/nAgcH5FF28wmQqs9g4DRb7B6g8hGBEZWGKeGZ1Xo99L0B5iOsX8EG8spaWl2oZOf+5HZHQkG1hG0y3+tr3IvtDAMnpNIjfc+sNSRohp3bFjR9jjjd5ggLavd7Q3JBHm+nXB6A3KCP2bhxH6DyPisUYBKVAURQvsxsZG5OTkoEePHm3O37ntttugKAo2bdoUc/r06PdG6LXb7fB6vWGPjfVaxPqgIZ7jypUrAUDby2EUF/qgAs6/PtG2SeI1Ea9fvIEVz7IXLbCiLVv6eBfzI9p2JFbYR6KqKtxutxaxegsKCtq8HuJDRySRHxj0zynW62YUo5HPV/y/vW2hWI7b2y6LIHvrrbcAMLBIlqAPrEjz8/PD3jzEyeJ6xQmfkYEFnPsUKw4tid2/YmXLVGC1N57IPW/tvaHrA0vsco+2ARLPI57AEm8k4jmJN3i/36/9W8w7sftev3FKJrCiBUB7fyu40MCK57WP/DAQqVjOxLDaG7bRczN6QxIbeP3j9MFcUFAQdW+P0bphNA36Q+eFhYVRn2PkMnTfffdBURRUVVVpt3355ZdQFAX9+vUzHGd7qKoKv9+P0tLSsDdy/fZAzNtYz8noDVPMj9LSUi2g9HtfIxHxIcIi1jIWOa/jDax4lr1Y52DFGnY84xHLV6Lf0hSH0YqLi7W9d3a7PWw44gOZ/rUw+sCgf05GIdjecixeJ0GsbYq4P975E7k8MbBIVhDrEGEkIlT0io2N0con3vj150EIMh1Y7e0BiQysaOgDq70NkHge8cxrEWtiWGIjqf82opifRhvnZAKrvQBoL7BiHeqJRTKB1Z4XEliRb0hG56oIgsFg2KFf8QYX+aWHRAMrHgXvvPMOFEXB+PHjtdt+9rOfQVEU/PSnPzUcZ6Loz4cSxBNYRvNdH1jxrDNiPPrAiranw+yBFe0bo5F7pZNBf/hQvwdULL/68/yMPjDon1Os1y3aaxW5PCSyLWRgEUuSSGB9+eWX2om8woaGhrDh6FcosWKLQwD6lSeRwIqMjXTuwRL3JxNYYjojz4/SE++8BsI/uYuTTvXTKfYOGr35ZyKwxF63eE6I15NMYCU67GiHZoyeW+QbUiKHAktLS9uclJtMYCW6F/CKK66Aoij46quvAAADBgyAoij47LPPEhpOLMS8F+fBpCKwxB4so0NYAnG6gpX3YEVud5JF7JmK3Abphx/rA0OqAyvWtoCBdR4GlkVJJLDiGU7kyic+KUWuzJErTOQ5R3oiv0J+IYHV3mOjfYswGpHnYIk3V6NPqonOa/GpVuwB1O8Z0s8/RWm7qz8TgaU/yTiRix0m8iYnojPe4YthR9urJt68I9/Y9G9IRueqxCLy/JJEAks8NlakG/Hoo49CURQsWbIEH374IRRFwU033ZTQMNpDLLuR38AzIt7AAmD4TVFBoudgiXUm8hysaK9/RwWWmC6jx6qqGvM8tESI9iUecbv+EjlGHxhSHVjR1pvI89v02zojxKF/cTkJBhbJCtIdWOJNI/KQXGRgxdpY6S+rALQNLBEg8QSW/hspkdOk/ypysoEV6+Re/TVe4kFsFMUbvH6jJ+a32DBHHjKJFljRDq2kIrCA8/OwvWtneb3euA/j6adNvMHGChD9Bl0MW39pA4GYh0af5PV7AqK9oUcLLhEF4k0ukcACzn8DL5FI3bFjBxRFwbXXXovZs2dDURT813/9V9x/L9bTaOivZSRIVWDF8y1CMf/b+xZh5N7yeKMg3YEltjtOp7NNROnPq2svsPx+f8zQF9uYaN/stNvtMT8wpCOwYm0L9b8nK+ZPtG8R2mw2fouQZBfpDqxoRL7hi93WinL+auPV1dUoLCzU9lpECyz9oUj9ycDiW1WR18HSn09SVlYWdnJ05PlMiQZWKBTSPpGXl5drV44uKirSbo93XuvfTIz+Rj+8yA1zZGDpY00/jwSpCizg/Bumopy7+OCaNWu0+VBeXt4mGBN5k9M/3ul0hg27uroaZWVlhucJieuyVVdXIxgMYs2aNVHjVKA/udto2lwuF0pKSlBdXY2amhoEAgFtmTXaixJvYOnXycgrkK9ZswaFhYWGb1ojR44MW16OHTtmOD4j9HvO9MttdXV12EU09eNNVWDp15mioiJt3GK8NpstbPkWw7bb7drrX11dHXadKT1i2GJd1//oeUcFFpDYdicaYpji0gyRy4WYX7Eu8REtwPSPSVVgiXUoclso7tPHlJgHYj0Vz0t/LS0BA4tkBWYJLCD8U7JQ/OZWrMASf6sPDv0ntMjAEuPXP15chDLyZNNEAws494krclrERi/ReR0Zl3rEBsnoW4uRgQUg7PygyOlIZWABCLvYYqQ2mw1utzvqT+VEEjltqqq2uWaTXqMYi3b5gVjnOonn3d7lFYzU741JNLCAc8tQrPlnFIUrV67UHvPd73436vMyIhQKafMq1mumJ1WBFev5OhyOdmMh2nzXjy9y2MnE/YUGlvgbo+2O0XUDjaisrIy6XMSaX0D4h9j2znVKVWC5XC7DbXN+fr7hHkijdcpo2WNgkazg9OnTUFU1od+bizUcsQu3PZqbm6GqapsrTKuqimAwiNraWu1Kz8C5jYP4vxhXS0uL4bDF8zl79iwA4Pjx4zEP5+gRGwkxrsbGxpjz5syZM1BVtc01iJqbm7Fjxw4EAueu2CzuFyeZxosYv9FvlIl5aDTPW1paoKqqdjV7PUaveazXTozHaFjtUV9fr82HYDBo+NzFc4w2/GjT1tzcrB3aipzP+vHrb6+vr0dtba32My6xiHWuSuQ01NbWora21nAa2ls3xPxtbW1tc199fb12BexgMBi2TkQilsVY60Y8iOcT6zUDYq8bsZYZMT9OnjzZ5j7x+oh5GQv9679jx464tmGRy1oiy56YbvH6iucY7fcDY82faOf8xbsdVlU1bNnfsWNHzGVDjDPWBwYgvtct2nIc+Xz12+zI7WF7z81o26lHbN/0H5qzGQYWsTxiD1Oiv6VHrIm4Ov6FfvggJBaxzgVMJe2dh0kyBwOLWIbKykqUlZVpvxof61g/kY9QKKSd/5PoNb0IiYbRdsfr9bZ7LuCFoqoqqqur+YHBxDCwiGWIdv6G0bF+Ihf685DE+X+EpIJY2510faiLHGci30wlHQcDi1gK/XF+cQ5DtHMpiDyI853EbxkSkko6erujv0wOPyyYFwYWIYQQQkiKYWARQgghhKQYBhYhhBBCSIphYBFCCCGEpBgGFiGEEEJIimFgEUIIIYSkGAYWIYQQQkiKYWBdAA0NDTh06FCmJ4MQQgghJoOBdQF4vV7tJzdaWlrw9ddfZ3iKSEdz8uRJy/wwKYmfpqYm7UfHiTzE+6P3xDqcPXsWTU1NSf0tA+sCYGARBpacMLDkhIElHwysDMHAIgwsOWFgyQkDSz4YWBmCgUUYWHLCwJITBpZ8MLAyBAOLMLDkhIElJwws+WBgZQgGFmFgyQkDS04YWPJhicDas2cPZs2ahT59+uCSSy7BzTffjDfffDPl4/F4PFAUxdCcnJyEhsXAIgwsOWFgyQkDSz6yPrA2bdoEm82G6dOno6amBtu2bUNFRQVsNhuWLVuW0nF5PB7MmTMHhw8fNjQRGFiEgSUnDCw5YWDJR9YH1m233YYpU6a0ub2mpgZdu3ZFKBRK2bg8Hg/mzZuXkmExsAgDS046IrB27tyJH/7whygrK6Mmce7cuRmfBnreH/7wh9i9e3da18OsD6x169ahtrbW8L4hQ4bghRdeSNm4GFgklTCw5KQjAuuW274D5aZZUP6vh1Jq5JgS3HbH+LSuh1kfWLEYO3YsFi9enLLhMbBIKmFgyUlHBNb4iXdB+b4fyqpTVAZ/vgs5V12Pi7r1pPHatRu6dOmCnj17tmuPHj3g8XgSXg8tG1gtLS2w2+0pPdnd4/HA5XJhwoQJuPLKK2Gz2eByuVBVVZXwsBhYhIElJwwsmnJLX8c/jr0NjY2NNA2+8sormDBhQsLroWUDa8GCBbj++uvR2tqasmGuWLECTqcTq1evxq5du7B3716sWrUKvXr1wsKFCxMaltfrxfe+9z0Eg0HU1dXhj3/8I4LBIILBIE6ePImzZ89Si3vixAmcPn0649NBO9bGxka0tramdRzjJzCwpLL0ddxYeEfK3utIOG+88QYmTpyY8HrY2tqKxsZG7f+JYNrAWrZsGXJzc7F9+/aUDvfs2bNoaWlpc3tdXR06d+6Mbdu2xT0sr9eLb3zjGxg6dChuuOEGDB06VPPPf/4zTpw4QS3u4cOHcfz48YxPB+1Yv/rqKzQ1NaV1HN8ZfycDSyYZWGnljTfewPjx4xNeD5uamvDVV19p/08EUwbWkiVL0LNnT2zevLlDxzt27FgsXbo07sfzECHhIUI54SFCysDKLsQerEQ5e9ZChwgXLVqE3NxcbNq0qcPHPXPmTMydOzfuxzOwCANLThhYlIGVXUgdWK2trXj44YfRu3dvbN26Vbt9586d2LNnT4dMg8vlQkVFRdyPZ2ARBpacMLAoAyu7kDawWlpaMG3aNAwePLhNTM2bNw/l5eUx/766uhoDBw7E8OHD8fHHHyc1DbW1tcjJycH+/fvj/hsGFmFgyQkDizKwsgspA6uxsRFjx47FmDFjcODAAbS0tIRZVlbWbmCNGzdO+z1BETxGHD16FIMGDcLChQuxYcMGbN++HXV1daioqEBubm5C518BDCzCwJIVBhZlYGUXUgbWJ598EvXHl4XtBdbGjRuRl5eHESNGhB1ejKS1tRXvvfceHnjgAeTn5yM3Nxf9+/dHUVER3njjjYSnnYFFGFhywsCiHRVYfr8ffa/sh/4OZ0YtGDES+/btS+syn06kDKxshoFFGFhywsCiHRVY991fDGXak1B+viuj9vzmBLzyyitpXebTCQMry2BgEQaWnDCwaEcF1v3FJVBK/ivj05c7elpGA8vn86GwsBCFhYUoKipCIBBI6O8ZWFkGA4swsOSEgUVTLgMrJj6fD4FAAIFAAJWVlVAUBfX19XH/PQMry2BgEQaWnDCwaMrN4sAS8aNHVdV2z5++EFwuF9566624H8/AyjIYWISBJScMLJpyszSwVFWF3W5HMBhsc5/D4YDf729z+6FDh+B2u6N69OjRqOuFqqqorq6Gw+HgHiwrw8AiDCw5YWDRlJulgeX3++FwOAyX4dLSUhQXF7e5PdnA8vl8cDgcUBQFXq83ofWJgZVlMLAIA0tOGFg05WZpYLndbrhcLu3/JSUl2t6syPtSRSgUgsPh4CFCK8PAIgwsOWFg0ZRrkcBSFEU7HytaYO3cuTPmtS8/++yzdtcPt9sNn88X7+rEwMo2GFiEgSUnDCyacrM0sAKBAOx2O1RVRTAYhMPhgNvtBgAUFhbG/HWVRAiFQmH/djqdCX2rkYGVZTCwCANLThhYNOVmaWABQHFxMZxOp3ZSu8PhQEFBAfLz88PC6EJwOBxwOp0oKCho92fxjGBgZRkMLMLAkhMGFk25WRxYAHD69GmcPn067P+pprm5Gc3NzUn9LQMry2BgEQaWnDCwaMqNElj/7i5HL+e16HvD6Iz6jUH52Lx5c1qX+XTCwMoyGFiEgSUnDCzaUYHV0tKC2trajPunP/0prct7umFgZRkMLMLAkhMGFu2owCKpgYGVZTCwCANLThhYtKMCa9OmTbjxm8Nx0+hRGfWuO+/AwYMH07rMpxMGVpbBwCIMLDlhYNGOCqz777sP3//HAVh/z8iMevPVV2b0x54vFAZWlsHAIgwsOWFg0Y4KrJLiGXjqO9fj89LbM+pdNwwwTWD97Gc/w8svv5zQ3zCwsgwGFmFgyQkDizKwMsNTTz2FK664AtOmTUvo7xhYWQYDizCw5ISBRRlY5wkGg20uKKqqKmpqalK6Tuzbtw8jRoyA2+3GPffck9DfMrCyDAYWYWDJCQOLMrDOY7fbtR941uNwOLTfJdRz4sQJBAKBqEa7mKjL5UIwGMTjjz/OwLI6DCzCwJITBhZlYJ0jEAjA4XAYLsOlpaUoLi5uc/v+/fvhcrmieuDAgTZ/8+yzz+KRRx4BAAaWDDCwCANLThhYlIF1DrfbDZfLpf2/tLRU25sVeV+yfPHFFygoKND+z8CSAAYWYWDJCQOLMrDOERlRiqJohwWjBdbOnTuhKEpUP/vss7DH//KXvzR8XN++feNenxhYWQYDizCw5ISBRRlY5/D7/bDb7VBVFaFQCDabDZWVlQCA4uJiw0OEFwr3YEkAA4swsOSEgUUZWOfJz89HQUEBnE4nfD4fHA4HCgsLYbPZ2ny7MBUwsCSAgUUYWHLCwKIMrHBCoRDq6+sBAM3NzQiFQlG/DZgJGFhZBgOLMLDkhIFFOyqwSn/4Q4wa2A9TR1ybUfMH9Mfbb7+d1mU+nTCwsgwGFmFgyQkDi3ZUYDU1NWHdunUZ96233krr8p5uGFhZBgOLMLDkhIFFOyqwSGpgYGUZDCzCwJITBhbtqMD69NNPMe2fpuHe6fdm1H+d+69Jh4YZYGBlGQwswsCSEwYW7ajAmv696eh7y1UY/K/DM+oV1/U1xY89JwsDK8tgYBEGlpwwsGhHBdaM4hkYNLsAN784KaM6xlyd0cDy+XwoLy8PMxEYWFkGA4swsOSEgUUZWB2L+J1Ct9utmQgMrCyDgUUYWHLCwKIMrPMEg8E2FxRVVRU1NTUpWx9EXCULAyvLYGARBpacMLAoA+s8TqdT+/1BPTabTfvhZz0nTpxAIBCIqtEFSl0uFwoLC1FeXp5UuDGwsgwGFmFgyQkDizKwzhEIBGCz2QyX4eLiYu09Us/+/fu1Q35GHjhwIOqw3G639lM8icDAyjIYWISBJScMLMrAOofb7YbL5dL+7/P5tMOFkfelClVV4XA4Err4KQMry2BgEQaWnDCwKAPrHJERpSiKdrgwWmAlswcrEpfLBZ/PF+/qxMDKNhhYhIElJwwsysA6h8/ng9PpBHBuz5LNZtPCp7i4GMXFxW3+5syZM1BVNaqtra0x1w1VVeF0OrF58+a41ycGVpbBwCIMLDlhYFEG1nny8/NRVFSEgoICuN1u2O12lJWVRT3JPRmcTidKSkpQUlICp9OJKVOmJPT3DKwsg4FFGFhywsCiDKzzNDc3Y8eOHaivrwdwbg/Tjh07oKpqytaH+vp6BINBBINBbTyJwMDKMhhYhIElJwws2lGBVVxSjF75vdFvwtUZ9dKrr4Df70/rMp9OGFhZBgOLMLDkhIFFOyqwjh07hqVLl2bcRE4oNyMMrCyDgUUYWHLCwKIdFVgkNTCwsgwGFmFgyQkDizKwsgupA2vPnj2YNWsW+vTpg0suuQQ333wz3nzzzbSNb+XKlRg0aBC6d++OgoKCpI4tM7AIA0tOGFiUgZVdSBtYmzZtgs1mw/Tp01FTU4Nt27ahoqICNpsNy5YtS/n41q5di6uuugqbNm2CqqpYv349evXqhdra2oSGw8AiDCw5YWBRBlZ2IW1g3XbbbYbXtKipqUHXrl3b/Er3hXD48GH07Nmzzd4xj8eDvLy8hIbFwCIMLDlhYFEGVnYhbWCtW7cu6t6jIUOG4IUXXkjZuKqqqjB06NA2t586dQrdunXD9u3b4x4WA4swsOSEgUUZWNmFtIEVi7Fjx2Lx4sUpG97s2bMxe/Zsw/tuuukmrFy5Mu5hMbAIA0tOGFiUgZVdMLAiaGlpgd1uT+nJ7uPHj4fH4zG8b9asWZg/f37cw2JgEQaWnDCwaDoCa+B1Q7F48WKaBktKSnDbbbclvB5aNrAWLFiA66+/vt0ffkyE0aNH4/nnnze8r6ysDA899FDcw2JgEQaWnDCwaMp9rhnKvcugjJ9H43V4Ea7OuwYLFixo1/nz5yd1PrclA2vZsihR/XgAACAASURBVGXIzc1N6JyoeBg1ahRWrVpleF9ZWRkefPDBuIfl9XoxYsQI/OAHP8CcOXMwd+5c/OAHP8APfvAD7Nu3Dy0tLdTiqqqKEydOZHw6aMd65MgRNDc3p3Uct985gYFFaSznvIo7J05K63rY3NyMI0eOaP9PBFMG1pIlS9CzZ09s3rw55cMeP348li5danjfrFmz8Mgjj8Q9LK/Xi8LCQng8HixevBhLliyBx+OBx+NBfX09vv76a2pxjx07hubm5oxPB+1YVVXFqVOn0jqOO+6cyMCiNJZzXsWEuyandT08deoUVFXV/p8IpgusRYsWITc3F5s2bUrL8B966KGoJ7mPGTMGy5cvj3tYPERIeIhQTniIkFITOOdVTJzc9hJPqcQShwhbW1vx8MMPo3fv3ti6dat2+86dO7Fnz56UjWfdunUYNmxYm9tPnTqFnJwcbNu2Le5hMbAIA0tOGFiUmkAGVvu0tLRg2rRpGDx4cJuYmjdvHsrLy2P+fXV1NQYOHIjhw4fj448/jvnYQ4cOoUePHnjnnXfCbq+oqMCAAQMS2mgysAgDS04YWJSaQAZWbBobGzF27FiMGTMGBw4caHOCWVlZWbuBNW7cOCiKAkVRtOCJxQsvvICBAweitrYWTU1N8Pv9sNlseP/99xOadgYWYWDJCQOLUhPIwIrNJ598osVRNNsLrI0bNyIvLw8jRowIO7wYi+XLlyMvLw/dunXDsGHD8OqrryY87QwswsCSEwYWpSaQgWVdGFiEgSUnDCxKTSADy7owsAgDS04YWJSaQAaWdWFgEQaWnDCwKDWBDCzrwsAiDCw5YWBRagIZWNaFgUUYWHLCwKLUBDKwrAsDizCw5ISBRakJZGBZFwYWYWDJCQOLUhPIwLIuDCzCwJITBhalJpCBZV0YWISBJScMLEpNIAPLujCwCANLThhYlJpABpZ1YWARBpacMLAoNYEMLOvCwCIMLDlhYFFqAhlY1oWBRRhYcsLAotQEMrCsCwOLMLDkhIFFqQlkYFkXBhZhYMkJA4tSE8jAsi4MLMLAkhMGFqUmkIFlXRhYhIElJwwsSk0gA8u6MLAIA0tOGFiUmkAGlnVhYBEGlpwwsCg1gQws68LAIgwsOWFgUWoCGVjWhYFFGFhywsCi1AQysKwLA4swsOSEgUWpCWRgWRcGFmFgyQkDi1ITyMCyLgwswsCSEwYWpSaQgWVdGFiEgSUnDCxKTSADy7owsAgDS04YWJSaQAaWdWFgEQaWnDCwKDWBDCzrkonAeumll7BmzRpqEletWoXVq1dnfDroef/yl7+kfT1kYFFqAq0aWD6fD36/X/u/qqooKiqCoigoKiqCqqrJDjpr6OjAWr16NeyD/hF9x91PTWKfW8+Z6emg5+xdOB3f6H9VWtdDgIFFqSm0amDZ7Xb4fD7t/8XFxbDZbHC73XA4HCgsLEx20FlDRweW1+tFzu3fz/xCTalZXXkU3br3SOt6CDCwKDWFVg0sRVEQDAa1/9vtdrjdbgBAMBiEoig4ffp0soPPChhYlJpMBhal8mjlwAqFQgCAQCDQJrgURbH8YUIGFqUmk4FFqTxaObDEOVilpaVwOBzafaqqMrDSND4GFqUxZGBRKo9WDazi4mLY7XYUFhZCURQtNADA7/dDUaz/BUUGFqUmk4FFqTxaNbBUVYXb7YbL5YLb7Q7bW+V2u1FcXJzsoLMGBhalJpOBRak8WjWwPvjgAwQCgZhaHQYWpSaTgUWpPFo1sPr06QNFUWJqdRhYlJpMBhal8mjVwDLag7Vhwwb86Ec/gtPpxK9+9atkB501MLAoNZkMLErl0aqBFYvq6mp885vfTMegTQUDi1KTycCiVB5lDKy///3vUBQFJ06cSMfgTQMDi1KTycCiVB5lDCxx4dHm5uZ0DN40MLAoNZkMLErl0aqBZXQO1ltvvYWf/vSn6N+/P2644YZkB501MLAoNZkMLErl0aqBFetbhNdddx1qa2uTHXTWwMCi1GQysCiVR6sGVkNDA1RVDbOhoQGnTp1KdpBZBwOLUpPJwKJUHq0SWMFgEG63G5s3bwYALF26FG63O6ZWh4FFqclkYFEqj1YJLJ/PB0VR8OSTTwLghUYBBhalppOBRak8WiWwOoJjx47B4/Ggf//+ePrpp9MyDo/HEzUIc3JyEhoWA4tSk8nAolQerRJYgUAA5eXlCRkvX3zxBebPnw+73Y7Jkydj+PDhaQ2sOXPm4PDhw4YmAgOLUpPJwKJUHq0SWG63u91DgskeInzqqacwY8YM/OUvfwEAFBUVpTWw5s2bl5JhMbAoNZkMLErl0SqBFUkwGITD4YDb7UYwGAQAhEIh+Hw+OBwO+P3+ZAfNwIoxPgYWpTFkYFEqj1YNLKfTicrKSsP7/H4/7HZ7soNmYMUYHwOL0hgysCiVRysGlqqqUBRF23NlOGBFgaqqSU1UugPL5XJhwoQJuPLKK2Gz2eByuVBVVZXwsBhYlJpMBhal8mjFwALOBZTP5zO8z+/3Q1EUnD59OqlhpzOwVqxYAafTidWrV2PXrl3Yu3cvVq1ahV69emHhwoUJDYuBRanJZGBRKo9WDazS0lLY7XZ4vV6EQiEA5/ZsrVmzBna7Hfn5+ckOOq2BdfbsWbS0tLS5va6uDp07d8a2bdviHpbX68Ull1wCu90Om80Gm80Gu90Ou92Ojz76CI2NjSl1yZIlDCxKY7nyKHK690j5uhfpwYMHcfz48bSO4zt33MnAojSWc17F+AkT07oeHj9+HAcPHtT+nwhJB5aqqpgyZYrhtwddLlfShweB9AZWLMaOHYulS5fG/Xiv14uHH34Yqqri73//O7766ivtZ4NaW1tTPn3cg0VpO3IPFqXyaNU9WILdu3fjrbfegs/nwyuvvILdu3df6CAzFlgzZ87E3Llz4348DxFSajIZWJTKo5UDSxwSvNALjUaSqcByuVyoqKiI+/EMLEpNJgOLUnm0amAFg0HY7fa0/BZhooFVXV2NgQMHYvjw4fj444+TGmdtbS1ycnKwf//+uP+GgUWpyWRgUSqPVg2sgoIC5OfnQ1VVuN1uuN1uAOd+UsfhcET9hqERp0+fRktLi+akSZNQWVkZdlusjdm4ceO0qBPBY8TRo0cxaNAgLFy4EBs2bMD27dtRV1eHiooK5ObmJnT+FcDAotR0MrAolUerBpaiKAgEAgAQFlhA4hcaLSsra/dnd95///2of79x40bk5eVhxIgR2Lp1a9THtba24r333sMDDzyA/Px85Obmon///igqKsIbb7wR9/QKGFiUmkwGFqXyaOXAEhcadbvdbfYcXciFRrMFBhalJpOBRak8WjmwxB6syspKFBQUaPcFg0EoipLwNSOyDQYWpSaTgUWpPFo1sFwul3ZYMBQKQVEUFBYWory8HHa7HQ6HI9lBZw0MLEpNJgOLUnm0amCFQqGw3yIMBAJwuVyw2WyYMmVKzN8ptAoMLEpNJgOLUnm0amARBhalppOBRak8WjGwvvzySyiKgnfffTepkVoFBhalJpOBRak8WjGwGhsboSgK3n777aRGahUYWJSaTAYWpfJoxcACgEmTJmHq1KnJ/rklYGBRajIZWJTKo1UDa/369bj00ksxbtw4PPfcc/jd736HQCAQptVhYFFqMhlYlMqjVQOrT58+7V593eowsCg1mQwsSuXRqoHV0NAAVVVjanUYWJSaTAYWpfJo1cAiDCxKTScDi1J5ZGBZFwYWpSaTgUWpPDKwrAsDi1KTycCiVB4ZWNaFgUWpyWRgUSqPDCzrwsCi1GQysCiVR6sEVn19fZvrXLWn1WFgUWoyGViUyqNVAsvn87V73SteB4uBRWlGZWBRKo9WCazGxkaEQqGEtDoMLEpNJgOLUnm0SmDFIhQKoaampo1Wh4FFqclkYFEqj1YNLFVVUVJSwkOEDCxKzSMDi1J5tGpglZaWwmazobKyEsXFxSguLobP50NxcTHsdjt8Pl+yg84aGFiUmkwGFqXyaNXAUhRFiyi32w23263d53a7UVBQkOygswYGFqUmk4FFqTxaObCCwSCAtoGlqioURUFjY2Oyg88KGFiUmkwGFqXyaOXAEte6crvdmDLl/JMUgaWqarKDzwoYWJSaTAYWpfJo1cByOBzaIcJAIABFUeD1elFTU4OSkhLYbLZkB501MLAoNZkMLErl0aqBVVlZGXZYsLi4OOwbhDzJPT3jY2BRGkMGFqXyaNXAMqK+vh6hUMjyhwYFDCxKTSYDi1J5tGJgHT16FC6XCx9++GFSI7UKDCxKTSYDi1J5tGJgHTp0CIqi4I9//GNSI7UKDCxKTSYDi1J5tGJgAcD111+Pp59+Otk/twQMLEpNJgOLUnm0amBt2bIFeXl5qK2tTXYQWQ8Di1KTycCiVB6tGlh9+vSJ+TuE/C3C9IyPgUVpDBlYlMqjVQOrqqoKPp8vplaHgUWpyWRgUSqPVg0swsCi1HQysCiVR6sGVp8+faJepqG2tpaHCNM0PgYWpTFkYFEqj1YOrA8++MDwvs2bNzOw0jQ+BhalMWRgUSqPVg6saHuw/uM//gP9+/dPdtBZAwOLUpPJwKJUHq0UWJWVle1+c1D4s5/9LKkJyiYYWJSaTAYWpfJopcCqra2F2+2G2+1Gjx49MGvWLO3/Qq/Xi08++SSpick2GFiUmkwGFqXyaKXA0lNVVYWDBw8m++eWgIFFqclkYFEqj1YNrBtvvBGPPvqo4X1erxdDhw5NdtBZAwOLUpPJwKJUHq0aWIqiRD3J/fPPP4eiKGhpaUl28FkBA4tSk8nAolQerRhYX3zxBRRFwZ49ewzvP3HiBBRFwfHjx5OaqGyBgUWpyWRgUSqPVgyspqYmKIqC1157zfD+P/zhD1AUBSdOnEhqorIFBhalJpOBRak8WjGwAKCwsBA33HADjhw5Enb7iRMncMstt+Caa65JdtBZAwOLUpPJwKJUHq0aWFu2bMGll16Kbt26YfLkyViwYAGmTp2K3NxcKIqC119/PeFhHjt2DB6PB/3798fTTz+d7KTFxcqVKzFo0CB0794dBQUF8Pv9CQ+DgUWpyWRgUSqPVg0sADhw4ACmTp2KSZMmaRcY7d27N9avX5/QcL744gvMnz8fdrsdkydPxvDhw9MaWGvXrsVVV12FTZs2QVVVrF+/Hr169UJtbW1Cw2FgUWoyGViUyqOVAytVPPXUU5gxYwb+8pe/AACKiorSFliHDx9Gz5498eabb4bd7vF4kJeXl9CwGFiUmkwGFqXyaOXACgQCKCoq0vZeCcRV3ZMlnYFVVVVleI2uU6dOoVu3bti+fXvcw2JgUWoyGViUyqNVA8vv90NRFLhcLhQXF4cFls/ng91uT3bQaQ2s2bNnY/bs2Yb33XTTTVi5cmXcw2JgUWoyGViUyqNVA8vpdGpxEQgEwgJLVVUoioLm5uakhp3OwBo/fjw8Ho/hfbNmzcL8+fPjHhYDi1KTycCiVB6tGliKoiAYDAJoG1jiflVVkxp2OgNr9OjReP755w3vKysrw0MPPRT3sLxeL6ZMmYLXXnsN69evx4YNG/Daa6/htddeQ0NDA86cOZNSly1bxsCiNJb/P7BSve5FeuzYMZw+fTqt47hjwkQGFqWxnPMqJk6anNb18PTp0zh27Jj2/0S4oMAKBAIA2gZWMBiEoihobGxMatjpDKxRo0Zh1apVhveVlZXhwQcfjHtYXq8X11xzDe68807cfvvtuOOOO3DnnXfizjvvxK5du3Dy5MmUunTpUgYWpbFceRQ53XukfN2L9NChQzhx4kRax3H7+DsZWJTGcs6ruHPiXWldD0+cOIFDhw5p/0+EpAPL4XCgqKgIQNvAKioqgsPhSHbQaT9EuHTpUsP7Zs2ahUceeSTuYfEQIaUmk4cIKZVHqx4iFCe5FxYWoqSkBIqiwOv1orCwEIqiJHXhTkE6A+uhhx6KepL7mDFjsHz58riHxcCi1GQysCiVR6sGFnBuz1V+fr52mQZFUeBwOC4oroD0Bta6deswbNiwNrefOnUKOTk52LZtW9zDYmBRajIZWJTKo5UDS4+qqkmfcxVJooFVXV2NgQMHYvjw4fj4449jPvbQoUPo0aMH3nnnnbDbKyoqMGDAgIQ2mgwsSk0mA4tSebRyYKmqivLychQWFmqWlJQgFAolNJzTp0+jpaVFc9KkSaisrAy7LdbGbNy4cdoeNBE8sXjhhRcwcOBA1NbWoqmpCX6/HzabDe+//35C083AotRkMrAolUerBlYwGITdbtcuNiq02Wxh3zCMh7KysrDDjEbGip+NGzciLy8PI0aMwNatW+Ma5/Lly5GXl4du3bph2LBhePXVV+OeXgEDi1KTycCiVB6tGlgFBQXIz89vs7dKVVW4XK4LupJ7tsDAotRkMrAolUerBlasvVShUCjsQqRWhYFFqclkYFEqj1YOrGgBJX4qh4GV+vExsCiNIQOLUnm0amAVFxejpKTE8D632438/PxkB501MLAoNZkWCqw7JkyE8tCLULxfUUqNfOAFTJg0Oa3rYUYCq7q6Gk6nE06nE+Xl5ZqFhYWw2+2orq5GTU2NZrK/S2hmGFiUmkwLBdY111yDizt3ReecHpRSAy/u3AWDBw1K63qYsUOEiZjItwqzBQYWpSbTQoHVt09vrJ48HJ+X3k4pNfCXdxXgG337pnU9NMWFRmWEgUWpyWRgUSqNUgVWKBTCmjVrLH9yu4CBRanJZGBRKo2WDayCgoKwq6aLH38W+ny+ZAedNTCwKDWZDCxKpdGygRV5XpW48GgwGITb7eaFRtM0PgYWpTFkYFEqjZYMLHGdK/HNQHFhUb/fH3Z/c3NzUhOVLTCwKDWZDCxKpdGSgQWE78Hy+XxhwSXut+KlGfQwsCg1mQwsSqXRsoHlcDhQUlKCYDCoHR4UiD1a3IOV+vExsCiNIQOLUmm0bGBFntSuPx+rsrISNpst2UFnDQwsSk0mA4tSabRsYAHnzrUKhUJtDgUGg0EpLtXAwKLUZDKwKJVGSweW7DCwKDWZDCxKpZGBZWEYWJSaTAYWpdLIwLIwDCxKTSYDi1JpZGBZGAYWpSaTgUWpNDKwLExHB9aSJUtwUf/rodw0i1Jq5LfuR+fOXdK6HgIMLErNIAPLwnR0YE2dOhXf+ofLsGTcEEqpgY8XXotOSvo3awwsSjMvA8vCZCKwvjfMkfGFmlKzunPOrfg/nTqldT0EGFiUmkEGloVhYFFqLhlYlMojA8vCMLAoNZcMLErlkYFlYRhYlJpLBhal8sjAsjAMLErNJQOLUnlkYFkYBhal5pKBRak8MrAsDAOLUnPJwKJUHhlYFoaBRam5ZGBRKo8MLAvDwKLUXDKwKJVHBpaFYWBRai4ZWJTKIwPLwjCwKDWXDCxK5ZGBZWEYWJSaSwYWpfLIwLIwDCxKzSUDi1J5ZGBZGAYWpeaSgUWpPDKwLAwDi1JzycCiVB4ZWBaGgUWpuWRgUSqPDCwLw8Ci1FwysCiVRwaWhWFgUWouGViUyiMDy8IwsCg1lwwsSuWRgWVhGFiUmksGFqXyyMCyMAwsSs0lA4tSeWRgWRgGFqXmkoFFqTwysCwMA4tSc8nAolQeGVgWhoFFqblkYFEqjwysBPjtb3+LUaNGoUePHhgwYACeeOKJlG/EPB4PFEUxNCcnJ6FhMbAoNZcMLErlkYEVJzt27IDdbsfatWtx9OhR1NXVYciQIXjqqadSOh6Px4M5c+bg8OHDhiYCA4tSc8nAolQeGVhxMnr0aMyfPz/strq6OuTk5GDv3r0pG4/H48G8efNSMiwGFqXmkoFFqTwysOKgvr4enTp1MtyDdOutt+KZZ55J2bgYWJRaVwYWpfLIwIqDl156CYMGDTK87yc/+QmmTZuWsnExsCi1rgwsSuWRgRUHixcvxl133WV439q1a3HjjTembFwejwculwsTJkzAlVdeCZvNBpfLhaqqqoSHxcCi1FwysCiVRwZWHDz66KP47ne/a3hfdXU1rrnmmpSNa8WKFXA6nVi9ejV27dqFvXv3YtWqVejVqxcWLlyY0LAYWJSaSwYWpfLIwIqDBQsWYPr06Yb3VVdXIy8vL2XjOnv2LFpaWtrcXldXh86dO2Pbtm1xD8vr9eLaa6/F3XffjSlTpqCoqAh333037r77buzevRvNzc0pdfLkyQwsSmMoAivV616khw8fxsmTJ9M6jj69r2BgURrDX95VgG/06ZPW9fDkyZM4fPiw9v9EMEVgLV68GBMnTjS8b+3atRg5cmSHTMfYsWOxdOnSuB/v9XoxceJEVFVV4b//+7/x61//GlVVVaiqqsKRI0dw+vTplFpUVMTAojSGIrBSve5F2tDQgK+//jqt4+jTm3uwKI2l2IOVzvXw66+/RkNDg/b/RDBFYP3617+OepL7okWLMHXq1A6ZjpkzZ2Lu3LlxP56HCCk1lzxESKk88hBhHBw4cACdOnXCkSNH2tx3yy23YMWKFR0yHS6XCxUVFXE/noFFqblkYFEqjwysOBk1ahQee+yxsNu2bNmCLl26YPfu3TH/trq6GgMHDsTw4cPx8ccfJzX+2tpa5OTkYP/+/XH/DQOLUnPJwKJUHhlYcfLpp5+iV69eePnll9HY2IiPPvoI1113HZYsWdLu344bN077PUERPEYcPXoUgwYNwsKFC7FhwwZs374ddXV1qKioQG5ubkLnXwEMLErNJgOLUnlkYCXAu+++i5EjR6J79+5wOBwoLy+PayO2ceNG5OXlYcSIEdi6dWvUx7W2tuK9997DAw88gPz8fOTm5qJ///4oKirCG2+8kfD0MrAoNZcMLErlkYFlYRhYlJpLBhal8sjAsjAMLErNJQOLUnlkYFkYBhal5pKBRak8MrAsDAOLUnPJwKJUHhlYFoaBRam5ZGBRKo8MLAvDwKLUXDKwKJVHBpaFYWBRai4ZWJTKIwPLwjCwKDWXDCxK5ZGBZWEYWJSaSwYWpfLIwLIwDCxKzSUDi1J5ZGBZGAYWpeaSgUWpPDKwLAwDi1JzycCiVB4ZWBaGgUWpuWRgUSqPDCwLw8Ci1FwysCiVRwaWhWFgUWouGViUyiMDy8IwsCg1lwwsSuWRgWVhMhFY3/6Hy1Bx2xBKqYFP3HItOinp36wxsCjNvAwsC9PRgTV8+HBcMsCGPoVXUUoN7H1zf3S6iHuwKJVBBpaF6ejAGjp0KPqNc+DmFydRSg389i/vhNKJe7AolUEGloVhYFFqLhlYlMojA8vCMLAoNZcMLErlkYFlYRhYlJpLBhal8sjAsjAMLErNpZUC64rLL8OUwX3xo1FXU0oNnDyoL3pfcUVa10MGVoZgYFFqLq0UWF26dsEV3+qHq+4eRCk18IrRV6JrTpe0rocMrAzBwKLUXFopsLrmdMWQH9+Y8XlKqVkd8qORyOnWNa3rIQMrQzCwKDWXDCxK5ZGBZWEYWJSaSwYWpfLIwLIwDCxKzSUDi1J5ZGBZGAYWpeaSgUWpPDKwLAwDi1JzycCiVB4ZWBaGgUWpuWRgUSqPDCwLw8Ci1FwysCiVRwaWhWFgUWouGViUyiMDy8IwsCg1lwwsSuWRgWVhGFiUmksGFqXyyMCyMAwsSs0lA4tSeWRgWRgGFqXmkoFFqTwysCwMA4tSc8nAolQeGVgWhoFFqblkYFEqjwwsC8PAotRcMrAolUcGloVhYFFqLhlYlMojA8vCMLAoNZcMLErlkYFlYRhYlJpLBhal8sjAsjAMLErNJQOLUnlkYFkYBhal5pKBRak8MrAsDAOLUnPJwKJUHhlYCfDb3/4Wo0aNQo8ePTBgwAA88cQTaduIrVy5EoMGDUL37t1RUFAAv9+f8DAYWJSaSwYWpfLIwIqTHTt2wG63Y+3atTh69Cjq6uowZMgQPPXUUykf19q1a3HVVVdh06ZNUFUV69evR69evVBbW5vQcBhYlJpLBhal8sjAipPRo0dj/vz5YbfV1dUhJycHe/fuTdl4Dh8+jJ49e+LNN98Mu93j8SAvLy+hYTGwKDWXDCxK5ZGBFQf19fXo1KkTDh8+3Oa+W2+9Fc8880zKxlVVVYWhQ4e2uf3UqVPo1q0btm/fHvewGFiUmksGFqXyyMCKg5deegmDBg0yvO8nP/kJpk2blrJxzZ49G7Nnzza876abbsLKlSvjHhYDi1JzycCiVB4ZWHGwePFi3HXXXYb3rV27FjfeeGPKxjV+/Hh4PB7D+2bNmtXmMGUsGFiUmksGFqXyyMCKg0cffRTf/e53De+rrq7GNddck7JxjR49Gs8//7zhfWVlZXjooYfiHpbX68WcOXPQ1NSEI0eOQFVVNDU1pW3jO3ToUHyj8B/w7V/eSSk1cNQz37FMYHXu2hnX/nBExucppWb12rnfRFcGVmwWLFiA6dOnG95XXV2d8MnnsRg1ahRWrVpleF9ZWRkefPDBuIfl9XrRo0cPQz/66CM0Njam1H//939Hp4s6QemkUEqNVBTkdM9J+boX6cGDB3H8+PG0jmPQ4MGZn5+UmllFwZAbrk/renj8+HEcPHhQ+38imCKwFi9ejIkTJxret3btWowcOTJl4xo/fjyWLl1qeN+sWbPwyCOPxD2sjj5ESMzHyZMncebMmUxPBulgOmIPFjEfib7Bkuwn6/dg/frXv456kvuiRYswderUlI3roYceinqS+5gxY7B8+fK4h8XAIgwsOWFgyQkDSz6yPrAOHDiATp064ciRI23uu+WWW7BixYqUjWvdunUYNmxYm9tPnTqFnJwcbNu2Le5hMbAIA0tOGFhywsCSj6wPLODcuVGPPfZY2G1btmxBly5dsHv37ph/W11djYEDB2L48OH4+OOPYz720KFD6NGjB955552w2ysqKjBgwICENpoMLMLAkhMGlpwwsOTDEoH16aefolevXnj55ZfR2NiIjz76CNdddx2WLFnS7t+OGzcOBat0zwAACghJREFUinLuhDcRPLF44YUXMHDgQNTW1qKpqQl+vx82mw3vv/9+QtOsD6y6urp2Q5BYj9/97neor6/P9GSQDua1115LeqNLspeqqqpMTwLpYJqamvD6668n9bemCSwAePfddzFy5Eh0794dDocD5eXlcX1K3LhxI/Ly8jBixAhs3bo1rnEtX74ceXl56NatG4YNG4ZXX3014enVB9bMmTPx3HPPJTwMkt3ceuutePvttzM9GaSDcTgc+OyzzzI9GaQDOXPmDC666KJMTwbpYHbv3o2BAwcm9bemCqxsg4FFGFhywsCSDwaWnDCwMgQDizCw5ISBJR8MLDlhYGWIDz74AHa7HTk5OejcuTO6dOmCnJwcKpEXX3wxX3cJvfjii9G1a9eMTwftOLt27YqLLroo49NBO/51v/jii7X//+lPf4q7ERhYF0hLSwuam5sppZRSanET+fYwA4sQQgghJMUwsAghhBBCUgwDixBCCCEkxTCwCCGEEEJSDAOLEEIIISTFMLAukDNnzuA3v/kNRo4ciXvuuSfTk0M6gD179mDWrFno06cPLrnkEtx888148803Mz1ZJM28/vrrmDp1Kvr27Yt+/fph5syZOHDgQKYni3Qgq1evhqIoePrppzM9KSSNeDwe7ef3Is3JyYl7OAysJGlqasLKlSsxYMAAFBQUYMqUKQwsCdi0aRNsNhumT5+OmpoabNu2DRUVFbDZbFi2bFmmJ4+kiTlz5qBfv354+umnsWvXLnz66acoLi5GXl4ejh8/nunJIx1AY2Mj+vXrh8GDBzOwLI7H48GcOXNw+PBhQ+OFgZUkW7ZsQWFhIf73f/8XwLmrujOwrM9tt92GKVOmtLm9pqYGXbt2RSgU6viJImln27ZtbULq7NmzuPHGG+H1ejM0VaQjWbRoEe677z7ce++9DCyL4/F4MG/evAseDgMrRTCw5GDdunWora01vG/IkCF44YUXOniKSCZZuHAhZsyYkenJIGlm//79uPTSS/H5558zsCSAgWUyGFhk7NixWLx4caYng3QgM2fORFlZWaYng6SZe++9F4899pj2bwaWtWFgmQwGlty0tLTAbrfzZHeJUFUV/fr1w3vvvZfpSSFp5MMPP0SfPn20Q8QMLOvj8XjgcrkwYcIEXHnllbDZbHC5XKiqqkpoOAysFMHAkpsFCxbg+uuvR2tra6YnhaSZpqYm1NXV4ZZbbsG//du/ZXpySBo5e/YsRo0aheeee067jYFlfVasWAGn04nVq1dj165d2Lt3L1atWoVevXph4cKFcQ+HgZUiGFjysmzZMuTm5mL79u2ZnhTSARQWFmpf2b799tvx97//PdOTRNLEiy++iCFDhuDMmTPabQws63P27Fm0tLS0ub2urg6dO3fGtm3b4hoOAytFMLDkZMmSJejZsyc2b96c6UkhHUhzczM++OADfOtb38KNN97IPZcW5OTJk+jfv3+bw/4MLLkZO3Ysli5dGtdjGVgpgoElH4sWLUJubi42bdqU6UkhGeLYsWPo3r076urqMj0pJMU8+eSTuOmmm6CqKhoaGjTvvvtuVFRUoKGhAadOncr0ZJIOZubMmZg7d25cj2VgpQgGljy0trbi4YcfRu/evbF161bt9p07d2LPnj0ZnDKSCfLz8+H3+zM9GSTFPPjgg1Gv5i3kniz5cLlcqKioiOuxDKwUwcCSg5aWFkybNg2DBw9uE1Pz5s1DeXl5hqaMZIJ9+/ahS5cu+OKLLzI9KaSD4CFCeamtrUVOTg72798f1+MZWEnS2tqKlpYWzaVLl2Lq1Klht+lPjCTZT2NjI8aOHYsxY8bgwIEDYa91S0sLysrKGFgW5Ouvv8Z3vvMdPP744/j973+PAwcO4MCBA3jppZcwYMAAPPHEE5meRJJmvv76a209/6d/+idUVlZyG29Rjh49ikGDBmHhwoXYsGEDtm/fjrq6OlRUVCA3Nzfu868ABlbSVFdXt7v7mBtea/HJJ5+0+5ozsKzJtm3b8OMf/xjf/OY3kZubi969e2PcuHH4n//5n0xPGukA7r33Xh4ilITW1la89957eOCBB5Cfn4/c3Fz0798fRUVFeOONNxIaFgOLEEIIISTFMLAIIYQQQlIMA4sQQgghJMUwsAghhBBCUgwDixBCCCEkxTCwCCGEEEJSDAOLECIV+/fvR01NDf72t78BALZs2aL9mxBCUgUDixAiBU1NTfj2t7+tXcPI4/EAAK6++mrt34QQkioYWIQQKViyZAlyc3Nx8ODBsNuXL1+ODz74IENTRQixKgwsQogU/PM//zOGDx+e6ckghEgCA4sQYmmOHz8Ol8uFyy+/HD179oTL5YLL5dJ+sHX69Ol4+eWXAQCTJ0/GL37xC7z55pvo27cvFEXB6NGj8Yc//CFsmM8++yzy8/PRrVs3XHnllSgpKcHJkydjTseTTz6JkpISfPrpp8jPz4eiKBg4cCCeeuqp9DxxQkhGYWARQixPIBDA2LFjMWjQIAQCAQQCAZw4cQJA+DlYl112GfLy8nDDDTfg+eefx9q1a+FwOGC327XHu91uKIqC++67D+vXr4fX60Xfvn3xrW99K+Y0zJgxA5dddhkuu+wy/PznP8fGjRsxdepUKIqCZ555Jr0zgBDS4TCwCCFScM8992DEiBFtbo8MLLvdjuPHj2v3/+lPf4KiKPjtb3+Lffv2QVEU/Od//mfYMHbv3g1FUfD2229HHf+MGTOgKAp+9atfhd1+yy234IorrriQp0YIMSEMLEKIFMQbWDNmzGjzmG7dusHn88Hn80FRFOzbt6/NY6677jp4vd6o458xYwa6devW5vbf/OY3UBQFX331VSJPhxBichhYhBApiCewLr/8cvz4xz9u85jLL78czz77LBYtWqRd5sHIxx57LOr4Z8yYgauvvrrN7bW1tVAUBX/9618v4NkRQswGA4sQIgWpCKyf/vSnUBQFR44cSXj8999/P/r06dPm9nfffReKomgn3RNCrAEDixAiBakIrA0bNkBRFDz33HMJj1+cgxV5ePHxxx+HoihoampKeJiEEPPCwCKESEEqAgsABg8ejEsuuQQrVqzQvlkInLtg6datW6OOXwTWuHHjtEs6bN26Fbm5uXjwwQcv6LkRQswHA4sQIgWpCqy9e/di6NChUBQFPXr00K6xpSgKtmzZEnX8M2bMwJAhQ/D9738fl1xyCYYPHw5FUTBp0qSwUCOEWAMGFiFECpqamsIuvyA4duwYWlpaAAANDQ2GFwxtaGjQHiMIhUKa+/btw6lTp2KOXwQWAGzfvh2bNm3C559/nuzTIYSYHAYWIYR0APfff78WWIQQ68PAIoSQDoCBRYhcMLAIIaQD+P3vfw+/35/pySCEdBAMLEIIIYSQFMPAIoQQQghJMQwsQgghhJAUw8AihBBCCEkx/w9cUnIilD2tMQAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "![Screenshot_2021-04-01%2020210402%20-%20LocalFourierAnalysisOfPMultigridwithLFAToolkitjl-Copy1%20-%20Jupyter%20Notebook.png](attachment:Screenshot_2021-04-01%2020210402%20-%20LocalFourierAnalysisOfPMultigridwithLFAToolkitjl-Copy1%20-%20Jupyter%20Notebook.png)"
   ]
  },
  {
   "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.1 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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <center>Backup</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": true,
    "hideOutput": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Agressive Coarsening</center>\n",
    "### <center>1D Scalar Diffusion - coarsen by half</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "using LFAToolkit\n",
    "using LinearAlgebra\n",
    "using Plots\n",
    "\n",
    "# setup\n",
    "numberpairs = 4\n",
    "numbersmoothers = 3\n",
    "maxeigenvalues = zeros(numberpairs, numbersmoothers)\n",
    "\n",
    "for (finep, coarsep) in [(3, 2), (5, 3), (7, 4), (9, 5)]\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 = 100\n",
    "maxeigenvalue = 0\n",
    "θ_min = -π/2\n",
    "θ_max = 3π/2\n",
    "\n",
    "# compute effect of smoother\n",
    "for w in 3:5\n",
    "# -- 1D --\n",
    "if dimension == 1\n",
    "    # setup\n",
    "    ω = [w]\n",
    "    v = [1, 1]\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[coarsep-1, w-2] = max(maxeigenvalues[coarsep-1, w-2], eigenvalues...)\n",
    "        end\n",
    "    end\n",
    "# -- 2D --\n",
    "elseif dimension == 2\n",
    "    # setup\n",
    "    ω = [w]\n",
    "    v = [1, 1]\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[coarsep-1, w-2] = max(maxeigenvalues[coarsep-1, w-2], eigenvalues...)\n",
    "        end\n",
    "    end\n",
    "end\n",
    "end\n",
    "end\n",
    "\n",
    "# plot\n",
    "xrange = [2, 4, 6, 8]\n",
    "bar(\n",
    "    xrange,\n",
    "    xlabel=\"fine p\",\n",
    "    xtickfont=font(12, \"Courier\"),\n",
    "    maxeigenvalues[:, 1],\n",
    "    ytickfont=font(12, \"Courier\"),\n",
    "    ylabel=\"spectral radius\",\n",
    "    label=\"ω = 3\",\n",
    "    title=\"P-Multigrid with Chebyshev Smoothing Symbol\"\n",
    ")\n",
    "bar!(\n",
    "    xrange,\n",
    "    xlabel=\"fine p\",\n",
    "    xtickfont=font(12, \"Courier\"),\n",
    "    maxeigenvalues[:, 2],\n",
    "    ytickfont=font(12, \"Courier\"),\n",
    "    ylabel=\"spectral radius\",\n",
    "    label=\"ω = 4\",\n",
    "    title=\"P-Multigrid with Chebyshev Smoothing Symbol\"\n",
    ")\n",
    "bar!(\n",
    "    xrange,\n",
    "    xlabel=\"fine p\",\n",
    "    xtickfont=font(12, \"Courier\"),\n",
    "    maxeigenvalues[:, 3],\n",
    "    ytickfont=font(12, \"Courier\"),\n",
    "    ylabel=\"spectral radius\",\n",
    "    label=\"ω = 5\",\n",
    "    title=\"P-Multigrid with Chebyshev Smoothing Symbol\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## <center>Agressive Coarsening</center>\n",
    "### <center>2D Scalar Diffusion - coarsen by half</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": true,
    "hidePrompt": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "using LFAToolkit\n",
    "using LinearAlgebra\n",
    "using Plots\n",
    "\n",
    "# setup\n",
    "numberpairs = 4\n",
    "numbersmoothers = 3\n",
    "maxeigenvalues = zeros(numberpairs, numbersmoothers)\n",
    "\n",
    "for (finep, coarsep) in [(3, 2), (5, 3), (7, 4), (9, 5)]\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 = 100\n",
    "maxeigenvalue = 0\n",
    "θ_min = -π/2\n",
    "θ_max = 3π/2\n",
    "\n",
    "# compute effect of smoother\n",
    "for w in 3:5\n",
    "# -- 1D --\n",
    "if dimension == 1\n",
    "    # setup\n",
    "    ω = [w]\n",
    "    v = [1, 1]\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[coarsep-1, w-2] = max(maxeigenvalues[coarsep-1, w-2], eigenvalues...)\n",
    "        end\n",
    "    end\n",
    "# -- 2D --\n",
    "elseif dimension == 2\n",
    "    # setup\n",
    "    ω = [w]\n",
    "    v = [1, 1]\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[coarsep-1, w-2] = max(maxeigenvalues[coarsep-1, w-2], eigenvalues...)\n",
    "        end\n",
    "    end\n",
    "end\n",
    "end\n",
    "end\n",
    "\n",
    "# plot\n",
    "xrange = [2, 4, 6, 8]\n",
    "bar(\n",
    "    xrange,\n",
    "    xlabel=\"fine p\",\n",
    "    xtickfont=font(12, \"Courier\"),\n",
    "    maxeigenvalues[:, 1],\n",
    "    ytickfont=font(12, \"Courier\"),\n",
    "    ylabel=\"spectral radius\",\n",
    "    label=\"ω = 3\",\n",
    "    title=\"P-Multigrid with Chebyshev Smoothing Symbol\"\n",
    ")\n",
    "bar!(\n",
    "    xrange,\n",
    "    xlabel=\"fine p\",\n",
    "    xtickfont=font(12, \"Courier\"),\n",
    "    maxeigenvalues[:, 2],\n",
    "    ytickfont=font(12, \"Courier\"),\n",
    "    ylabel=\"spectral radius\",\n",
    "    label=\"ω = 4\",\n",
    "    title=\"P-Multigrid with Chebyshev Smoothing Symbol\"\n",
    ")\n",
    "bar!(\n",
    "    xrange,\n",
    "    xlabel=\"fine p\",\n",
    "    xtickfont=font(12, \"Courier\"),\n",
    "    maxeigenvalues[:, 3],\n",
    "    ytickfont=font(12, \"Courier\"),\n",
    "    ylabel=\"spectral radius\",\n",
    "    label=\"ω = 5\",\n",
    "    title=\"P-Multigrid with Chebyshev Smoothing Symbol\"\n",
    ")"
   ]
  }
 ],
 "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
}