{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 3.5.1 Operator applications for Discontinuous Galerkin discretizations - linear transport"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In the previous units we already saw DG discretizations. However, we only solved stationary problems in these cases so that DG discretizations have only been treated in terms of solving linear systems (where HDG has some advantages over classical DG discretizations).\n",
    "\n",
    "Here, however, we treat DG discretizations in the case where operator applications are more relevant than the solution of linear systems, as e.g. in explicit time integration schemes. \n",
    "Note that for continuous Galerkin methods, in [unit 3.2](../unit-3.2-navierstokes/navierstokes.ipynb) we already saw simple ways on how to use operator applications within an IMEX scheme. We will develop this further and take a look at efficient ways to do it for DG methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In this unit we start with a model PDE problem: The unsteady scalar linear transport problem\n",
    "\n",
    "Find $u: [0,T] \\to V_D := \\{ u \\in L^2(\\Omega), b \\cdot \\nabla u \\in L^2(\\Omega), u|_{\\Gamma_{in}} = u_D\\}$, s.t.\n",
    "\\begin{equation}\n",
    "\\int_{\\Omega} \\partial_t u v +  b \\cdot \\nabla u v = \\int_{\\Omega} f v \\qquad \\forall v \\in V_0 = \\{ u \\in L^2(\\Omega), b \\cdot \\nabla u \\in L^2(\\Omega), u|_{\\Gamma_{in}} = 0\\}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We consider the unit square $(0,1)^2$ and a given advection velocity $b$. \n",
    "Accordingly the inflow boundary is $\\Gamma_{in} = \\{ x \\cdot y = 0\\}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# import libraries, define geometry and generate mesh\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import *\n",
    "from netgen.occ import *\n",
    "shape = Rectangle(1,1).Face()\n",
    "shape.edges.Min(X).name=\"left\"\n",
    "shape.edges.Max(X).name=\"right\"\n",
    "shape.edges.Min(Y).name=\"bottom\"\n",
    "shape.edges.Max(Y).name=\"top\"\n",
    "geom = OCCGeometry(shape, dim=2)\n",
    "mesh = Mesh(geom.GenerateMesh(maxh=0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We consider higher order discretizations and set $k=5$ here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "order = 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Warning**\n",
    "\n",
    "* We will make some timing measurements. \n",
    "* To avoid the impact of initialization overheads you may want to increase $k$ and decrease the mesh size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Convection and boundary values:\n",
    "\n",
    "* To make cases with `CoefficientFunctions` we use the `IfPos`-`CoefficientFunction`. `IfPos(a,b,c)`: `a` decides on the evaluation. If `a` is positive `b` is evaluated, otherwise `c`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "b = CF((1+sin(4*pi*y),2))\n",
    "ubnd = 0.1 * IfPos(x-0.125,IfPos(0.625-x,1+cos(8*pi*x),0),0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Draw(b,mesh,vectors = {\"grid_size\" : 32}, min = 2, max = 3, autoscale=False,\n",
    "     height = \"30vh\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We consider an Upwind DG discretization (in space):\n",
    "\n",
    "Find $u: [0,T] \\to V_h := \\bigoplus_{T\\in\\mathcal{T}_h} \\mathcal{P}^k(T)$ so that\n",
    "\n",
    "$$\n",
    "  \\sum_{T} \\int_T \\partial_t u v - b \\cdot \\nabla v u + \\int_{\\partial T} b_n \\hat{u} v = \\int_{\\Omega} f v, \\quad \\forall v \\in V_h.\n",
    "$$\n",
    "\n",
    "Here $\\hat{u}$ is the Upwind flux, i.e. $\\hat{u} = u$ on the outflow boundary part $\\partial T_{out} = \\{ x\\in \\partial T \\mid b(x) \\cdot n_T(x) \\geq 0 \\}$ of $T$ and $\\hat{u} = u^{other}$ else, with $u^{other}$ the value from the neighboring element."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "In the following we **only treat the explicit time integration case here** (take a look at [unit 2.8](../unit-2.8-DG/DG.ipynb) for solving linear systems with DG formulations)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Explicit time stepping with a DG formulation\n",
    "\n",
    "Upwind-in-space combined with an explicit Euler scheme in time yields:\n",
    "$$\n",
    "\\sum_{T} \\int_T u^{n+1} v = \\sum_{T} \\int_T u^{n} v + \\Delta t \\sum_{T} \\left\\{ \\int_T  b \\cdot \\nabla v u \n",
    "+ \\int_{\\partial T} b_n \\hat{u} v \\right\\} + \\Delta t \\int_{\\Omega} f v, \\quad \\forall v \\in V_h,\n",
    "$$\n",
    "\n",
    "\n",
    "$$\n",
    "  M u^{n+1} = M u^{n} - \\Delta t C(u^n) + \\Delta t f\n",
    "$$\n",
    "\n",
    "In our first example we set $u_0 = f = 0$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# roughly estimated time step:\n",
    "dt = 1e-3 / (order+1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Computing convection applications $C u^n$\n",
    "\n",
    "We recall: \n",
    "\n",
    "* We can define the bilinear form **without setting up a matrix** with storage. (`nonassemble=True`)\n",
    "* A `BilinearForm` is allowed to be **nonlinear in the 1st argument** (for operator applications)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "VT = L2(mesh,order=order)\n",
    "u,v = VT.TnT()\n",
    "c = BilinearForm(VT, nonassemble=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* To access the **neighbor** function we can use `u.Other()`, cf. [unit 2.8](../unit-2.8-DG/DG.ipynb).\n",
    "* To incorporate **boundary conditions** ($\\hat{u}$ on inflow boundaries), we can use the argument `bnd` of `.Other()` (not possible in [unit 2.8](../unit-2.8-DG/DG.ipynb) )\n",
    "* If there is no neighbor element (boundary!) the `CoefficientFunction` `bnd` is evaluated. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "n = specialcf.normal(mesh.dim)\n",
    "upw_flux = b*n * IfPos(b*n, u, u.Other(bnd=ubnd))\n",
    "dS = dx(element_boundary=True)\n",
    "c += - b * grad(v) * u * dx + upw_flux * v * dS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Bilinearform for operator applications\n",
    "Due to the `nonassemble=True` in the constructor of the `BilinearForm` we can use `c.mat * x` to realize the operator application without assembling a matrix first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(VT) \n",
    "#Draw(gfu,mesh, settings = { \"subdivision\" : 10 },\n",
    "#     min=-0.1,max=2.1,autoscale=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Operator application (equivalent to assemble and mult but faster)\n",
    "res = gfu.vec.CreateVector()\n",
    "res.data = c.mat * gfu.vec "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    },
    "tags": [
     "raises-exception"
    ]
   },
   "outputs": [],
   "source": [
    "# this does not work:\n",
    "c2 = BilinearForm(VT)\n",
    "res.data = c2.mat * gfu.vec "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Solving mass matrix problems (see also [unit-2.11](../unit-2.11-matrixfree/matrixfree.ipynb))\n",
    "\n",
    "We also have to deal with the mass matrix:\n",
    "\n",
    "* Need to invert the mass matrix\n",
    "* For DG methods the mass matrix is block diagonal, often even diagonal.\n",
    "* FESpace offers (if available):\n",
    "  * `Mass(rho)`: mass matrix as an operator  (not a sparse matrix)\n",
    "  * `Mass(rho).Inverse()`: inverse mass matrix as an operator (not a sparse matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "invm = VT.Mass(1).Inverse() # not a sparse matrix\n",
    "print(invm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The (simple) time loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def SolveAndVisualize(conv_op, t=0, tend=0.6, dt=dt, dt_sample=0.06, rhs_vec = None):\n",
    "    gfu.vec[:] = 0; cnt = 0;\n",
    "    sample_rate = int(floor(dt_sample*tend/dt))\n",
    "    gfu_t = GridFunction(gfu.space,multidim=0)\n",
    "    gfu_t.AddMultiDimComponent(gfu.vec)\n",
    "    if rhs_vec:\n",
    "        inv_rhs = invm * rhs_vec\n",
    "    while t < tend-0.5*dt:\n",
    "        if rhs_vec:\n",
    "            res = invm @ conv_op * gfu.vec + inv_rhs\n",
    "        else:\n",
    "            res = invm @ conv_op * gfu.vec\n",
    "        gfu.vec.data -= dt * res\n",
    "        t += dt; cnt += 1; print(\"\\r t=\",t,end=\"\")\n",
    "        if (cnt+1) % sample_rate == 0:\n",
    "            gfu_t.AddMultiDimComponent(gfu.vec)\n",
    "    return gfu_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%%time \n",
    "gfu_t = SolveAndVisualize(c.mat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Draw(gfu_t,mesh,interpolate_multidim=True,animate=True,\n",
    "     min=0, max=0.2, autoscale=False, deformation=True, height = \"60vh\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Skeleton formulation (cf. [unit 2.8](../unit-2.8-DG/DG.ipynb))\n",
    "\n",
    "So far we considered the DG formulation with integrals on the boundary of each element. Instead one could formulate the problem in terms of facet integrals where every facet appears only once. \n",
    "\n",
    "Then, the corresponding formulation is:\n",
    "\n",
    "Find $u: [0,T] \\to V_h := \\bigoplus_{T\\in\\mathcal{T}_h} \\mathcal{P}^k(T)$ so that\n",
    "$$\n",
    "  \\sum_{T} \\int_T \\partial_t u v - b \\cdot \\nabla v u + \\sum_{F\\in\\mathcal{F}^{inner}} \\int_{F} b_n \\hat{u} (v - v^{neighbor}) + \\\\\n",
    "  \\sum_{F\\in\\mathcal{F}^{bound}} \\int_{F} b_n \\hat{u} v = \\int_{\\Omega} f v , \\quad \\forall v \\in V_h.\n",
    "$$\n",
    "\n",
    "Here $\\mathcal{F}^{inner}$ is the set of interior facets and $\\mathcal{F}^{bound}$ is the set of boundary facets. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### integral types:\n",
    "The facet integrals are divided into inner facets and boundary facets.\n",
    "To obtain these integrals we combine the `DifferentialSymbol` `dx` or `ds` with `skeleton=True`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "dskel_inner  = dx(skeleton=True)\n",
    "dskel_bound  = ds(skeleton=True)\n",
    "# or:\n",
    "dskel_inflow = ds(skeleton=True, definedon=mesh.Boundaries(\"left|bottom\"))\n",
    "dskel_outflow = ds(skeleton=True, definedon=mesh.Boundaries(\"right|top\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Skeleton formulation of $C$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "c = BilinearForm(VT,nonassemble=True)\n",
    "c += -b * grad(v) * u * dx\n",
    "c += upw_flux * (v-v.Other()) * dskel_inner\n",
    "\n",
    "c += b*n * IfPos(b*n,u,ubnd) * v * dskel_bound\n",
    "#alternatively (if you know in/out beforehand):\n",
    "#c += b*n * ubnd * v * dskel_inflow\n",
    "#c += b*n * u * v * dskel_outflow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Next run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "gfu_t = SolveAndVisualize(c.mat)\n",
    "#Draw(gfu_t,mesh,interpolate_multidim=True,animate=True,\n",
    "#     min=0, max=0.2, autoscale=False, deformation=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## TraceOperator (Purpose: avoid neighbor access `Other()`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The trace operator maps from an `FESpace` to a compatible `FacetFESpace` by evaluating the traces. \n",
    "\n",
    "The **stored value** on a facet is the **sum of the traces** of the aligned elements. \n",
    "\n",
    "No need to access the neighbor (`.Other()`) later on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "VF = FacetFESpace(mesh, order=order)\n",
    "V = VT*VF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "trace = VT.TraceOperator(VF, False)\n",
    "#info output about the trace operator\n",
    "print(\"The space VT has {} dofs.\".format(VT.ndof))\n",
    "print(\"The trace space VF has {} dofs.\".format(VF.ndof))\n",
    "print(\"trace represents an {} x {} operator\".format(trace.height, trace.width))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "example application of the trace operator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "gfuF = GridFunction(VF)\n",
    "gfuF.vec.data = trace * gfu.vec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To make use of the trace operator, we split the previous (`element_boundary`) formulation into three parts:\n",
    "$$\n",
    "  \\sum_{T} \\int_T - b \\cdot \\nabla v u + \\sum_{T} \\int_{\\partial T} b_n \\hat{u}^* v\n",
    "  + \\sum_{F\\in\\mathcal{F}^{inflow}} \\int_{F} b_n u^{inflow} v \\quad \\text{ with } \\quad \\hat{u}^* = f(u^{left}+u^{right},u)\n",
    "$$\n",
    "\n",
    "* part 1: element local volume integrals (no interaction with neighbor elements) (on `VT`)\n",
    "* part 2: couplings on the facets (on `VT` $\\times$ `V`)\n",
    "* part 3: (rhs + ) boundary terms (on `VT`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Part 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "uT, vT = VT.TnT()\n",
    "c1 = BilinearForm(space=VT, nonassemble=True)                     # part 1\n",
    "c1 += -uT * b * grad(vT) * dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Part 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "uT,uF = V.TrialFunction() \n",
    "c2 = BilinearForm(trialspace=V, testspace=VT, nonassemble=True)   # part 2\n",
    "# here uf-u = u_me + u_other - u_me is the neighbor value\n",
    "# on bnd: uf = u_me and hence uf-u = 0\n",
    "c2 += b*n * IfPos(b*n, uT, uF-uT) * vT * dx(element_boundary=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Right hand side (Part 3):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "rhs = LinearForm(VT)                                              # part 3\n",
    "rhs += b*n * IfPos(b*n, 0, ubnd) * vT * dskel_bound # dskel_inflow\n",
    "rhs.Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Note that `c1` and `c2` act on different spaces. To make them compatible, we combine:\n",
    " * `TraceOperator`: VT $\\to$ VF\n",
    " * `Embedding`s: VT $\\to$ VT $\\times$ VF and VF $\\to$ VT $\\times$ VF (cf. [unit 2.10](../unit-2.10-dualbasis/dualbasis.ipynb))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The Embeddings are (with $n_F = \\operatorname{dim}(V_F)$, $n_T = \\operatorname{dim}(V_T)$, $n_V=n_F+n_T$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "$$\n",
    "\\texttt{embV}: V_T \\to V_T \\times V_F, \\qquad \n",
    "\\texttt{embV} = \n",
    "\\underbrace{\\left( \\begin{array}{c} I \\\\ 0 \\end{array} \\right)}_{\\texttt{embT} \\in \\mathbb{R}^{n_V \\times n_T}}\n",
    "+ \n",
    "\\underbrace{\\left( \\begin{array}{c} 0 \\\\ I \\end{array} \\right)}_{\\texttt{embF} \\in \\mathbb{R}^{n_V \\times n_F}} \\underbrace{\\operatorname{tr}}_{\\texttt{trace}\\in \\mathbb{R}^{n_F \\times n_T}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "embT = Embedding(V.ndof, V.Range(0))\n",
    "embF = Embedding(V.ndof, V.Range(1))\n",
    "embV = embT + embF @ trace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The operators can now be combined easily. Recall:\n",
    " * `c1`: `VT` $\\to$ `VT`' $\\quad$ $\\bullet$ `c2`: `VT` $\\times$ `VF` $\\to$ `VT`' $\\quad$ $\\bullet$ $\\Rightarrow$ `c`: `VT` $\\to$ `VT`'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "c = c1.mat + c2.mat @ embV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The time loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "gfu_t = SolveAndVisualize(c, rhs_vec=rhs.vec)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "#Draw(gfu_t,mesh,interpolate_multidim=True,animate=True,\n",
    "#     min=0, max=0.2, autoscale=False, deformation=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Speed up with `geom_free` integrals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "With $u(x) = \\hat{u}(\\hat{x})$, $v(x) = \\hat{v}(\\hat{x})$ (directly mapped) and $\\mathbf{w}(x) = \\frac{1}{|\\operatorname{det}(F)|} F \\cdot \\hat{\\mathbf{w}}(\\hat{x})$ (Piola mapped) there holds\n",
    "$$\n",
    " \\int_T \\mathbf{w} \\cdot \\nabla u \\ v  = \\int_{\\hat{T}} \\hat{\\mathbf{w}} \\cdot \\hat{\\nabla} \\hat{u} \\ \\hat{v}\n",
    "$$\n",
    "\n",
    "and similarly for the facet integrals. \n",
    "\n",
    "Integrals are independent of the \"physical element\" (up to ordering) similar to integrals in [unit-2.11](../unit-2.11-matrixfree/matrixfree.ipynb).\n",
    "\n",
    "$\\leadsto$ allows to prepare intermediate computations for a large class of elements at once."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Piola wind:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfwind = GridFunction(HDiv(mesh, order=order))\n",
    "gfwind.Set(b)\n",
    "b=gfwind; # <- overwrite old name \"b\" "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Same as before, but with  `geom_free=True` flag:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# part 1\n",
    "uT, vT = VT.TnT()\n",
    "c1 = BilinearForm(space=VT, nonassemble=True, geom_free=True)\n",
    "c1 += -uT * b * grad(vT) * dx\n",
    "# part 2\n",
    "uT,uF = V.TrialFunction() \n",
    "c2 = BilinearForm(trialspace=V, testspace=VT, nonassemble=True, geom_free=True)\n",
    "# here uf-u = u_me + u_other - u_me is the neighbor value\n",
    "c2 += b*n * IfPos(b*n, uT, uF-uT) * vT * dx(element_boundary=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "c = c1.mat + c2.mat @ embV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "gfu_t = SolveAndVisualize(c, rhs_vec=rhs.vec)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "#Draw(gfu_t,mesh,interpolate_multidim=True,animate=True,\n",
    "#     min=0, max=0.2, autoscale=False, deformation=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "And now with the TaskManager:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "with TaskManager():\n",
    "    gfu_t = SolveAndVisualize(c, rhs_vec=rhs.vec)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.10"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}