{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Tutorial 6: Poisson equation (with DG)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "In this tutorial, we will learn\n",
    " - How to solve a simple PDE with a DG method\n",
    " - How to compute jumps and averages of quantities on the mesh skeleton\n",
    " - How to implement the method of manufactured solutions\n",
    " - How to integrate error norms\n",
    " - How to generate Cartesian meshes in arbitrary dimensions"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Problem statement\n",
    "\n",
    "The goal of this tutorial is to solve a PDE using a Discontinuous Galerkin (DG) formulation. For\n",
    "simplicity, we take the Poisson equation on the unit cube $\\Omega \\doteq\n",
    "(0,1)^3$ as the model problem, namely"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "$$\n",
    "\\left\\lbrace\n",
    "\\begin{aligned}\n",
    "-\\Delta u = f  \\ &\\text{in} \\ \\Omega,\\\\\n",
    "u = g \\ &\\text{on}\\ \\Gamma \\doteq \\partial\\Omega,\\\\\n",
    "\\end{aligned}\n",
    "\\right.\n",
    "$$\n",
    "where $f$ is the source term and $g$ is the prescribed Dirichlet boundary\n",
    "function. In this tutorial, we follow the method of manufactured solutions\n",
    "since we want to illustrate how to compute discretization errors. We take\n",
    "$u(x) = 3 x_1 + x_2^2 + 2 x_3^3 + x_1 x_2 x_3 $ as the exact solution of the\n",
    "problem, for which $f(x)= -2 - 12x_3 $ and $g(x) = u(x)$. The selected\n",
    "manufactured solution $u$ is a third order multi-variate polynomial, which\n",
    "can be represented exactly by the FE discretization that we are going to\n",
    "define below. In this scenario, the discretization error has to be close to\n",
    "the machine precision. We will use this result to validate the proposed\n",
    "implementation.\n",
    "\n",
    "## Numerical Scheme\n",
    "\n",
    "We consider a DG formulation to approximate the problem. In particular, we\n",
    "consider the symmetric interior penalty method (see, e.g. [1], for specific\n",
    "details). For this formulation, the approximation space is made of\n",
    "discontinuous piece-wise polynomials, namely\n",
    "\n",
    "$$\n",
    "V \\doteq \\{ v\\in L^2(\\Omega):\\ v|_{T}\\in Q_p(T) \\text{ for all } T\\in\\mathcal{T}  \\},\n",
    "$$\n",
    "where $\\mathcal{T}$ is the set of all cells $T$ of the FE mesh, and $Q_p(T)$\n",
    "is a polynomial space of degree $p$ defined on a generic cell $T$. For\n",
    "simplicity, we consider Cartesian meshes in this tutorial. In this case, the\n",
    "space $Q_p(T)$ is made of multi-variate polynomials up to degree $p$ in each\n",
    "spatial coordinate.\n",
    "\n",
    "In order to write the weak form of the problem, we need to introduce some\n",
    "notation. The sets of interior and boundary facets associated with the FE\n",
    "mesh $\\mathcal{T}$ are denoted here as $\\mathcal{F}_\\Lambda$ and\n",
    "$\\mathcal{F}_{\\Gamma}$ respectively. In addition, for a given\n",
    "function $v\\in V$ restricted to the interior facets $\\mathcal{F}_\\Lambda$, we\n",
    "introduce the well known jump and mean value operators,\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\lbrack\\!\\lbrack v\\ n \\rbrack\\!\\rbrack &\\doteq v^+\\ n^+ + v^- n^-,\\\\\n",
    "\\{\\! \\!\\{\\nabla v\\}\\! \\!\\} &\\doteq \\dfrac{ \\nabla v^+ + \\nabla v^-}{2},\n",
    "\\end{aligned}\n",
    "$$\n",
    "with $v^+$, and $v^-$ being the restrictions of $v\\in V$ to the cells $T^+$,\n",
    "$T^-$ that share a generic interior facet in $\\mathcal{F}_\\Lambda$, and $n^+$,\n",
    "and $n^-$ are the facet outward unit normals from either the perspective of\n",
    "$T^+$ and $T^-$ respectively.\n",
    "\n",
    "With this notation, the weak form associated with the interior penalty\n",
    "formulation of our problem reads: find $u\\in V$ such that $a(u,v) = l(v)$ for\n",
    "all $v\\in V$. The bilinear and linear forms $a(\\cdot,\\cdot)$ and $l(\\cdot)$\n",
    "have contributions associated with the bulk of $\\Omega$, the boundary facets\n",
    "$\\mathcal{F}_{\\Gamma}$, and the interior facets $\\mathcal{F}_\\Lambda$,\n",
    "namely\n",
    " $$\n",
    "\\begin{aligned}\n",
    "a(u,v) &= a_{\\Omega}(u,v) + a_{\\Gamma}(u,v) + a_{\\Lambda}(u,v),\\\\\n",
    "l(v) &= l_{\\Omega}(v) + l_{\\Gamma}(v).\n",
    "\\end{aligned}\n",
    "$$\n",
    "These contributions are defined as\n",
    "$$\n",
    "\\begin{aligned}\n",
    "a_{\\Omega}(u,v) &\\doteq \\sum_{T\\in\\mathcal{T}} \\int_{T} \\nabla v \\cdot \\nabla u \\ {\\rm d}T,\n",
    "\\\\\n",
    "l_{\\Omega}(v) &\\doteq \\int_{\\Omega} v\\ f \\ {\\rm d}\\Omega,\n",
    "\\end{aligned}\n",
    "$$\n",
    "for the volume,\n",
    "$$\n",
    "\\begin{aligned}\n",
    "a_{\\Gamma}(u,v)\n",
    "   \\doteq\n",
    "   & - \\sum_{F\\in\\mathcal{F}_{\\Gamma}} \\int_{F} v\\ (\\nabla u \\cdot n)  \\ {\\rm d}F\n",
    "\\\\ & - \\sum_{F\\in\\mathcal{F}_{\\Gamma}} \\int_{F} (\\nabla v \\cdot n)\\ u  \\ {\\rm d}F\n",
    "\\\\ & + \\sum_{F\\in\\mathcal{F}_{\\Gamma}} \\dfrac{\\gamma}{|F|} \\int_{F} v\\ u \\ {\\rm d}F,\n",
    "\\\\\n",
    "l_{\\Gamma}(v)\n",
    "\\doteq\n",
    "   & - \\sum_{F\\in\\mathcal{F}_{\\Gamma}} \\int_{F} (\\nabla v \\cdot n)\\ g  \\ {\\rm d}F\n",
    "\\\\ & + \\sum_{F\\in\\mathcal{F}_{\\Gamma}} \\dfrac{\\gamma}{|F|} \\int_{F} v\\ g \\ {\\rm d}F,\n",
    "\\end{aligned}\n",
    "$$\n",
    "for the boundary facets and,\n",
    "$$\n",
    "\\begin{aligned}\n",
    "a_{\\Lambda}(u,v)\n",
    "\\doteq\n",
    "   & - \\sum_{F\\in\\mathcal{F}_{\\Lambda}} \\int_{F} \\lbrack\\!\\lbrack v\\ n \\rbrack\\!\\rbrack\\cdot \\{\\! \\!\\{\\nabla u\\}\\! \\!\\} \\ {\\rm d}F\n",
    "\\\\ & - \\sum_{F\\in\\mathcal{F}_{\\Lambda}} \\int_{F} \\{\\! \\!\\{\\nabla v\\}\\! \\!\\}\\cdot \\lbrack\\!\\lbrack u\\ n \\rbrack\\!\\rbrack \\ {\\rm d}F\n",
    "\\\\ & + \\sum_{F\\in\\mathcal{F}_{\\Lambda}} \\dfrac{\\gamma}{|F|} \\int_{F} \\lbrack\\!\\lbrack v\\ n \\rbrack\\!\\rbrack\\cdot \\lbrack\\!\\lbrack u\\ n \\rbrack\\!\\rbrack \\ {\\rm d}F,\n",
    "\\end{aligned}\n",
    "$$\n",
    " for the interior facets. In previous expressions, $|F|$ denotes the diameter\n",
    " of the face $F$ (in our Cartesian grid, this is equivalent to the\n",
    " characteristic mesh size $h$), and $\\gamma$ is a stabilization parameter that\n",
    " should be chosen large enough such that the bilinear form $a(\\cdot,\\cdot)$ is\n",
    " stable and continuous. Here, we take $\\gamma = p\\ (p+1)$ as done in the\n",
    " numerical experiments in reference [2].\n",
    "\n",
    "## Manufactured solution\n",
    "\n",
    "We start by loading the Gridap library and defining the manufactured solution\n",
    "$u$ and the associated source term $f$ and Dirichlet function $g$."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Gridap\n",
    "u(x) = 3*x[1] + x[2]^2 + 2*x[3]^3 + x[1]*x[2]*x[3]\n",
    "f(x) = -2 - 12*x[3]\n",
    "g(x) = u(x)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We also need to define the gradient of $u$ since we will compute the $H^1$\n",
    "error norm later. In that case, the gradient is simply defined as"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "∇u(x) = VectorValue(3        + x[2]*x[3],\n",
    "                    2*x[2]   + x[1]*x[3],\n",
    "                    6*x[3]^2 + x[1]*x[2])"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "In addition, we need to tell the Gridap library that the gradient of the\n",
    "function `u` is available in the function `∇u` (at this moment `u` and `∇u`\n",
    "are two standard Julia functions without any connection between them). This\n",
    "is done by adding an extra method to the function `gradient` (aka `∇`)\n",
    "defined in Gridap:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "import Gridap: ∇\n",
    "∇(::typeof(u)) = ∇u"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    " Now, it is possible to recover function `∇u` from function `u` as `∇(u)`. You\n",
    " can check that the following expression evaluates to `true`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "∇(u) === ∇u"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Cartesian mesh generation\n",
    " In order to discretize the geometry of the unit cube, we use the Cartesian mesh generator available in Gridap."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "L = 1.0\n",
    "domain = (0.0, L, 0.0, L, 0.0, L)\n",
    "n = 4\n",
    "partition = (n,n,n)\n",
    "model = CartesianDiscreteModel(domain,partition)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The type `CartesianDiscreteModel` is a concrete type that inherits from\n",
    "`DiscreteModel`, which is specifically designed for building Cartesian\n",
    "meshes. The `CartesianDiscreteModel` constructor takes a tuple containing\n",
    "limits of the box we want to discretize plus a tuple with the number of cells\n",
    "to be generated in each direction (here $4\\times4\\times4$ cells). You can\n",
    "write the model in vtk format to visualize it (see next figure)."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    " writevtk(model,\"model\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "![](../assets/dg_discretization/model.png)\n",
    "\n",
    " Note that the `CaresianDiscreteModel` is implemented for arbitrary\n",
    " dimensions. For instance, the following lines build a\n",
    " `CartesianDiscreteModel` for the unit square $(0,1)^2$ with 4 cells per\n",
    " direction"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "domain2D = (0.0, L, 0.0, L)\n",
    "partition2D = (n,n)\n",
    "model2D = CartesianDiscreteModel(domain2D,partition2D)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "You could also generate a mesh for the unit tesseract $(0,1)^4$ (i.e., the\n",
    "unit cube in 4D). Look how the 2D and 3D models are built and just follow the\n",
    "sequence.\n",
    "\n",
    "## FE spaces\n",
    "\n",
    "On top of the discrete model, we create the discontinuous space $V$ as\n",
    "follows"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "order = 3\n",
    "V = TestFESpace(model,\n",
    "                ReferenceFE(lagrangian,Float64,order),\n",
    "                conformity=:L2)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We have select a Lagrangian, scalar-valued interpolation of order $3$ within\n",
    "the cells of the discrete model. Since the cells are hexahedra, the resulting\n",
    "Lagrangian shape functions are tri-cubic polynomials. In contrast to previous\n",
    "tutorials, where we have constructed $H^1$-conforming (i.e., continuous) FE\n",
    "spaces, here we construct a $L^2$-conforming (i.e., discontinuous) FE space.\n",
    "That is, we do not impose any type of continuity of the shape function on the\n",
    "cell boundaries, which leads to the discontinuous FE space $V$ of the DG\n",
    "formulation. Note also that we do not pass any information about the\n",
    "Dirichlet boundary to the `TestFESpace` constructor since the Dirichlet\n",
    "boundary conditions are not imposed strongly in this example.\n",
    "\n",
    "From the `V` object we have constructed in previous code snippet, we build\n",
    "the trial FE space as usual."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "U = TrialFESpace(V)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Note that we do not pass any Dirichlet function to the `TrialFESpace`\n",
    "constructor since we do not impose Dirichlet boundary conditions strongly\n",
    "here.\n",
    "\n",
    "## Numerical integration\n",
    "\n",
    "Once the FE spaces are ready, the next step is to set up the numerical\n",
    "integration. In this example, we need to integrate in three different\n",
    "domains: the volume covered by the cells $\\mathcal{T}$ (i.e., the\n",
    "computational domain $\\Omega$), the surface covered by the boundary facets\n",
    "$\\mathcal{F}_{\\Gamma}$ (i.e., the boundary $\\Gamma = \\partial \\Omega$), and\n",
    "the surface covered by the interior facets $\\mathcal{F}_{\\Lambda}$ (i.e. the\n",
    "so-called mesh skeleton). In order to integrate in $\\Omega$ and on its\n",
    "boundary $\\Gamma$, we use `Triangulation` and `BoundaryTriangulation` objects\n",
    "as already discussed in previous tutorials."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "Ω = Triangulation(model)\n",
    "Γ = BoundaryTriangulation(model)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here, we do not pass any boundary identifier to the `BoundaryTriangulation`\n",
    "constructor. In this case, an integration mesh for the entire boundary\n",
    "$\\Gamma$ is constructed by default (which is just what we need in this\n",
    "example).\n",
    "\n",
    "In order to generate an integration mesh for the interior facets\n",
    "$\\mathcal{F}_{\\Lambda}$, we use a new type of `Triangulation` referred to as\n",
    "`SkeletonTriangulation`. It can be constructed from a `DiscreteModel` object\n",
    "as follows:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "Λ = SkeletonTriangulation(model)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "As any other type of `Triangulation`, an `SkeletonTriangulation` can be\n",
    "written into a vtk file for its visualization (see next figure, where the\n",
    "interior facets $\\mathcal{F}_\\Lambda$ are clearly observed)."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "writevtk(Λ,\"strian\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "![](../assets/dg_discretization/skeleton_trian.png)\n",
    "\n",
    "Once we have constructed the triangulations needed in this example, we define\n",
    "the corresponding quadrature rules by passing the triangulations\n",
    "together with the desired degree to the `Measure` function."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "degree = 2*order\n",
    "\n",
    "dΩ = Measure(Ω,degree)\n",
    "dΓ = Measure(Γ,degree)\n",
    "dΛ = Measure(Λ,degree)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We still need a way to represent the unit outward normal vector to the\n",
    "boundary $\\Gamma$, and the unit normal vector on the interior faces\n",
    "$\\mathcal{F}_\\Lambda$. This is done with the `get_normal_vector` getter."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "n_Γ = get_normal_vector(Γ)\n",
    "n_Λ = get_normal_vector(Λ)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The `get_normal_vector` getter takes either a boundary or a skeleton\n",
    "triangulation and returns an object representing the normal vector to the\n",
    "corresponding surface. For boundary triangulations, the returned normal\n",
    "vector is the unit outwards one, whereas for skeleton triangulations the\n",
    "orientation of the returned normal is arbitrary. In the current\n",
    "implementation (Gridap v0.5.0), the unit normal is outwards to the cell with\n",
    "smaller id among the two cells that share an interior facet in\n",
    "$\\mathcal{F}_\\Lambda$.\n",
    "\n",
    "## Weak form\n",
    "\n",
    "With these ingredients we can define the different terms in the weak form.\n",
    "First, we start with the terms $a_\\Omega(\\cdot,\\cdot)$ , and\n",
    "$l_\\Omega(\\cdot)$ associated with integrals in the volume $\\Omega$. This is\n",
    "done as in the tutorial for the Poisson equation."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "a_Ω(u,v) = ∫( ∇(v)⊙∇(u) )dΩ\n",
    "l_Ω(v) = ∫( v*f )dΩ"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The terms $a_{\\Gamma}(\\cdot,\\cdot)$ and $l_{\\Gamma}(\\cdot)$ associated with\n",
    "integrals on the boundary $\\Gamma$ are defined using an analogous approach:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "h = L / n\n",
    "γ = order*(order+1)\n",
    "a_Γ(u,v) = ∫( - v*(∇(u)⋅n_Γ) - (∇(v)⋅n_Γ)*u + (γ/h)*v*u )dΓ\n",
    "l_Γ(v)   = ∫(                - (∇(v)⋅n_Γ)*g + (γ/h)*v*g )dΓ"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Note that in the definition of the functions `a_Γ` and `b_Γ`, we have used\n",
    "the object `n_Γ` representing the outward unit normal to the boundary\n",
    "$\\Gamma$. The code definition of `a_Γ` and `b_Γ` is indeed very close to\n",
    "the mathematical definition of the forms $a_{\\Gamma}(\\cdot,\\cdot)$ and\n",
    "$b_{\\Gamma}(\\cdot)$.\n",
    "\n",
    "Finally, we need to define the term $a_\\Lambda(\\cdot,\\cdot)$ integrated on\n",
    "the interior facets $\\mathcal{F}_\\Lambda$,"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "a_Λ(u,v) = ∫( - jump(v*n_Λ)⊙mean(∇(u))\n",
    "              - mean(∇(v))⊙jump(u*n_Λ)\n",
    "              + (γ/h)*jump(v*n_Λ)⊙jump(u*n_Λ) )dΛ"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Note that the arguments `v`, `u` of function `a_Λ` represent a test and trial\n",
    "function *restricted* to the interior facets $\\mathcal{F}_\\Lambda$. As\n",
    "mentioned before in the presentation of the DG formulation, the restriction\n",
    "of a function $v\\in V$ to the interior faces leads to two different values\n",
    "$v^+$ and $v^-$ . In order to compute jumps and averages of the quantities\n",
    "$v^+$ and $v^-$, we use the functions `jump` and `mean`, which represent the\n",
    "jump and mean value operators $\\lbrack\\!\\lbrack \\cdot \\rbrack\\!\\rbrack$ and\n",
    "$\\{\\! \\!\\{\\cdot\\}\\! \\!\\}$ respectively. Note also that we have used the\n",
    "object `n_Λ` representing the unit normal vector on the interior facets. As a\n",
    "result, the notation used to define function `a_Λ` is very close to the\n",
    "mathematical definition of the terms in the bilinear form\n",
    "$a_\\Lambda(\\cdot,\\cdot)$.\n",
    "\n",
    "Once the different terms of the weak form have been defined, we build and solve the FE problem."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "a(u,v) = a_Ω(u,v) + a_Γ(u,v) + a_Λ(u,v)\n",
    "l(v) = l_Ω(v) + l_Γ(v)\n",
    "\n",
    "op = AffineFEOperator(a, l, U, V)\n",
    "uh = solve(op)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Discretization error\n",
    "\n",
    "We end this tutorial by quantifying the discretization error associated with\n",
    "the computed numerical solution `uh`. In DG methods a simple error indicator\n",
    "is the jump of the computed (discontinuous) approximation on the interior\n",
    "faces. We compute and visualize the jump of these values as follows (see next\n",
    "figure):"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "writevtk(Λ,\"jumps\",cellfields=[\"jump_u\"=>jump(uh)])"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Note that the jump of the numerical solution is very small, close to the\n",
    "machine precision (as expected in this example with manufactured solution).\n",
    "![](../assets/dg_discretization/jump_u.png)\n",
    "\n",
    " A more rigorous way of quantifying the error is to measure it with a norm.\n",
    " Here, we use the $L^2$ and $H^1$ norms, namely\n",
    " $$\n",
    "\\begin{aligned}\n",
    " \\| w \\|_{L^2}^2 & \\doteq \\int_{\\Omega} w^2 \\ \\text{d}\\Omega, \\\\\n",
    " \\| w \\|_{H^1}^2 & \\doteq \\int_{\\Omega} w^2 + \\nabla w \\cdot \\nabla w \\ \\text{d}\\Omega.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The discretization error can be computed in this example as the difference of\n",
    "the manufactured and numerical solutions."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "e = u - uh"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We compute the error norms as follows. First, we implement the integrands of\n",
    "the norms we want to compute."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "l2(u) = sqrt(sum( ∫( u⊙u )*dΩ ))\n",
    "h1(u) = sqrt(sum( ∫( u⊙u + ∇(u)⊙∇(u) )*dΩ ))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Then, we compute the corresponding integrals with the `integrate` function."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "el2 = l2(e)\n",
    "eh1 = h1(e)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The `integrate` function returns a lazy object representing the contribution\n",
    "to the integral of each cell in the underlying triangulation. To end up with\n",
    "the desired error norms, one has to sum these contributions and take the\n",
    "square root. You can check that the computed error norms are close to machine\n",
    "precision (as one would expect)."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "tol = 1.e-10\n",
    "@assert el2 < tol\n",
    "@assert eh1 < tol"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## References\n",
    "\n",
    "[1] D. N. Arnold, F. Brezzi, B. Cockburn, and L. Donatella Marini. Unified analysis of discontinuous Galerkin methods for elliptic problems. *SIAM Journal on Numerical Analysis*, 39 (5):1749–1779, 2001. doi:[10.1137/S0036142901384162](http://dx.doi.org/10.1137/S0036142901384162).\n",
    "\n",
    "[2] B. Cockburn, G. Kanschat, and D. Schötzau. An equal-order DG method for the incompressible Navier-Stokes equations. *Journal of Scientific Computing*, 40(1-3):188–210, 2009. doi:[10.1007/s10915-008-9261-1](http://dx.doi.org/10.1007/s10915-008-9261-1)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "---\n",
    "\n",
    "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
   ],
   "metadata": {}
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.7"
  },
  "kernelspec": {
   "name": "julia-1.6",
   "display_name": "Julia 1.6.7",
   "language": "julia"
  }
 },
 "nbformat": 4
}