lecture-5.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 5: Finite element method (FEM)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## FEM ideology:\n",
    "\n",
    "* Discretize the domain: triangles (tetrahedra, quadrilaterals, ...)\n",
    "<img width=\"33%\" src=\"fig/FEM_mesh_example.jpg\">\n",
    "* Discretize the functions: linear on each triangle (i.e., they are piecewise linear, and in general, piecewise polynomial)\n",
    "* Discretize the equations: Galerkin, collocation\n",
    "* Solve: same methods as for FDM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# FEM in 1D:\n",
    "\n",
    "* Start with a domain $\\Omega = [0,1]$\n",
    "* Mesh (triangulation) ${\\mathcal T}$ consisting of\n",
    "    - segments $[x_{i}, x_{i+1}]$ (where $x_0,...,x_N$ are nodes)\n",
    "<img src=\"fig/FEM_mesh_1D.png\">\n",
    "\n",
    "* Idea: the material consists of small triangles (segments, tetrahedra) that act like basic building blocks of a material\n",
    "\n",
    "* Functions that are piecewise linear (more precisely, piecewise affine) with respect to $\\mathcal T$:\n",
    "$$\n",
    "    u_h(x) = \\sum_{i=1}^N c_i v_i(x)\n",
    "$$\n",
    "<img src=\"fig/FEM_basis_1D.png\">\n",
    "\n",
    "* In 2D it looks like this:\n",
    "<img width=\"50%\" src=\"fig/FEM_basis_2D.png\">\n",
    "\n",
    "* Such space of function is denoted as ${\\mathcal P}^1(\\mathcal T)$\n",
    "\n",
    "* Often we'll use\n",
    "$$\n",
    "\\mathcal P^1_0(\\mathcal T)\n",
    ":=\n",
    "\\{u\\in\\mathcal P^1(\\mathcal T) : u|_{\\Omega} = 0 \\}\n",
    "$$\n",
    "\n",
    "* Discretize the equations..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Discretization of Equations (Galerkin's method):\n",
    "\n",
    "First we need the variational (weak) form\n",
    "\n",
    "#### Variational form\n",
    "\n",
    "* Start with the equation (now I like to have minus in front of $\\Delta$):\n",
    "\\\\[\n",
    "\\begin{align*}\n",
    "-\\Delta u &= f\n",
    "\\\\\n",
    "u|_{\\Gamma} &= 0\n",
    "\\end{align*}\n",
    "\\\\]\n",
    "\n",
    "* Multiply the equation by a function $v$ (called the test function) such that $v|_{\\Gamma}=0$ and integrate (Green's theorem):\n",
    "$$\n",
    "-\\int_{\\Omega} (\\Delta u) v\n",
    "=- \\int_{\\Gamma} \\nabla u v \\cdot n + \\int_{\\Omega} \\nabla u \\cdot \\nabla v\n",
    "=\\int_{\\Omega} f v\n",
    "\\qquad\\forall v\n",
    "$$\n",
    "\n",
    "* because $v=0$ on the boundary, we have\n",
    "$$\n",
    "\\int_{\\Omega} \\nabla u \\cdot \\nabla v\n",
    "=\\int_{\\Omega} f v\n",
    "\\qquad\\forall v\n",
    "$$\n",
    "* This is called the **variational form**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Alternative way for Variational Form\n",
    "\n",
    "* Recall from Lecture 2 that the energy of interaction of a two-dimensional network of springs was something like (in 1D):\n",
    "\\\\[\n",
    "E^{\\rm 1d}_h(u) = \\sum_x k (u(x+h) - u(x))^2/2 - f(x) u(x),\n",
    "\\\\]\n",
    "where $f(x)$ is the external force acting on the mass $x$.\n",
    "\n",
    "* if we let $k=1/h$, scale $f$ accordingly, and let $h\\to 0$ then we have that\n",
    "$$\n",
    "E(u) = \\int_\\Omega |\\nabla u|^2/2 - \\int_\\Omega f u\n",
    "$$\n",
    "(this formula works in any dimension)\n",
    "\n",
    "* Now, we can ask a question:\n",
    "\\\\[\n",
    "\\text{find $u$ such that $E(u)$ is minimal}\n",
    "\\\\]\n",
    "\n",
    "* To find the equation for $u$, one can take a functional derivative (equivalent to the principle of virtual displacements of mechanics):\n",
    "$$\n",
    "\\langle \\partial E(u), v \\rangle\n",
    ":= \\lim_{\\epsilon\\to 0}\\frac{E(u+\\epsilon v)-E(u)}{\\epsilon}\n",
    "= \\int_\\Omega \\nabla u \\cdot \\nabla v - \\int_\\Omega f v\n",
    "= 0,\n",
    "\\qquad\\forall v\n",
    "$$\n",
    "so we derived the variational form without the strong form\n",
    "\n",
    "* $\\langle \\partial E(u), v \\rangle$ is called a Gateaux derivative (see also the Frechet derivative)\n",
    "\n",
    "* In general, the variational form is obtained directly from the energy minimum principle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### What if we have Neumann or mixed BCs\n",
    "\n",
    "(BC = boundary condition)\n",
    "\n",
    "* Suppose we have a problem\n",
    "$$\n",
    "\\begin{align*}\n",
    "-\\Delta u &= f\\qquad\\text{on $\\Omega$}\n",
    "\\\\\n",
    "u &= 0\\qquad\\text{on $\\Gamma_1$}\n",
    "\\\\\n",
    "u_n &= 0\\qquad\\text{on $\\Gamma_2$}\n",
    "\\end{align*}\n",
    "$$\n",
    "where $u_n$ is the normal derivative and $\\Gamma_1 \\cup\\Gamma_2 = \\partial\\Omega$.\n",
    "\n",
    "* The trick is then to introduce the space $X = \\{\\text{function }u : u|_{\\Gamma_1}=0 \\}$. Then we have, for $v\\in X$\n",
    "$$\n",
    "-\\int_{\\Omega} \\Delta u v\n",
    "= - \\int_{\\Gamma_1} u_n v - \\int_{\\Gamma_2} u_n v + \\int_{\\Omega} \\nabla u\\cdot \\nabla v\n",
    "$$\n",
    "\n",
    "* Magically, the first integral =0 because $v=0$ on $\\Gamma_1$ and the second integral =0 because $u_n=0$ on $\\Gamma_2$.\n",
    "\n",
    "* **Same variational formulation**, only spaces are different"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Galerkin Method, cont'd\n",
    "\n",
    "* Denote $A(u,v) = \\int_{\\Omega} \\nabla u\\cdot\\nabla v$ and $F(v) := \\int_\\Omega f v$.\n",
    "\n",
    "* Discrete equations:\n",
    "\\\\[\n",
    "A(u_h, v_h) = F(v_h)\n",
    "\\qquad \\forall v_h\\in \\mathcal P^1_0(\\mathcal T)\n",
    "\\\\]\n",
    "(That's it!)\n",
    "\n",
    "* But to implement, we do\n",
    "$$\n",
    "A(u_h, v_\\ell) = F(v_\\ell),\n",
    "\\qquad \\ell=1,...,N-1\n",
    "$$\n",
    "(i.e., it is enough to test only with basis functions)\n",
    "\n",
    "* Then we substitute $u_h = \\sum_{k=1}^{N-1} c_k v_k$:\n",
    "$$\n",
    "\\sum_{k=1}^{N-1} c_k A(v_k, v_\\ell) = F(v_\\ell) \\qquad \\ell=1,\\ldots,N-1,\n",
    "$$\n",
    "\n",
    "* Thus, $A(v_k, v_\\ell)$ are the elements of the matrix (often called the **stiffness matrix**), and $F(v_\\ell)$ are the components of the right-hand side\n",
    "\n",
    "* It remains to solve the linear system..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### An example of a practical problem\n",
    "\n",
    "* We want to compute the resistance of pumping the oil (very practical for Russia! :) )\n",
    "\n",
    "<img src=\"fig/borehole.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### The corresponding model problem\n",
    "\n",
    "* Cross-section of a tube:<br>\n",
    "<img src=\"fig/borehole_simple.png\">\n",
    "\n",
    "* Here $P=1$ is the (dimensionless) pressure drop on the tube\n",
    "\n",
    "* $\\int_\\Omega u$ will then be the total fluid flux through the tube"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## FEM in 2D\n",
    "\n",
    "* The mesh consists of triangles (tetrahedra, quadrilaterals, ...)\n",
    "<img width=\"33%\" src=\"fig/FEM_mesh_example.jpg\">\n",
    "\n",
    "    - The mesh is defined by, typically, 3 arrays:\n",
    "\n",
    "    <tt>nodes:</tt> <br>\n",
    "    $x_1$, $y_1$ <br>\n",
    "    $x_2$, $y_2$ <br>\n",
    "    ... <br>\n",
    "    $x_N$, $y_N$ <br>\n",
    "    (i.e., the first node has coords $(x_1, y_1)$, etc.)\n",
    "\n",
    "    <tt>triangles:</tt> <br>\n",
    "    $n_1$, $n_2$, $n_3$ <br>\n",
    "    $m_1$, $m_2$, $m_3$ <br>\n",
    "    ... <br>\n",
    "    (i.e., the first triangle has nodes $n_1$, $n_2$, $n_3$ as vertices, etc.)\n",
    "\n",
    "    <tt>faces:</tt> <br>\n",
    "    $n_1$, $n_2$\n",
    "    ... <br>\n",
    "    (i.e., the first face **on the boundary** has nodes $n_1$, $n_2$, etc. )\n",
    "\n",
    "* Basis functions look like this:\n",
    "\n",
    "<img src=\"fig/FEM_basis_2D.png\">    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Assembling the stiffness matrix\n",
    "\n",
    "Classical way:\n",
    "\n",
    "* Notice that $\\nabla v_k$ are piecewise constant:\n",
    "$$\n",
    "A_{k,\\ell}\n",
    "= \\int_{\\Omega} \\nabla v_k \\cdot \\nabla v_\\ell\n",
    "= \\sum_{T\\in\\mathcal T} |T| (\\nabla v_k|_T) \\cdot (\\nabla v_\\ell|_T),\n",
    "$$\n",
    "where $|T|$ is the volume of $T$\n",
    "\n",
    "* The algorithm would the look something like this:\n",
    "```\n",
    "for k = 1 to n\n",
    "  for l = 1 to n\n",
    "    for all T\n",
    "      if(k and l are nodes of T)\n",
    "      // otherwise the integral is zero\n",
    "        A(k,l) += |T|\n",
    "           *(\\nabla v_k|_T)\n",
    "           *(\\nabla v_l|_T)\n",
    "```\n",
    "\n",
    "* The loops in the algorithm are often reversed:\n",
    "```\n",
    "for all T\n",
    "    for k, vertices of T\n",
    "        for l, vertices of T\n",
    "            A(k,l) += <<as before>>\n",
    "```\n",
    "\n",
    "* Advantage: can now loop only over 3 vertices of each triangle, not waisting time\n",
    "* Now the we only need to consider the geometry of a single triangle\n",
    "\n",
    "<img src=\"fig/FEM_triangle.png\">\n",
    "\n",
    "* We have\n",
    "$$\n",
    "2 |T| = \\det \\left(\n",
    "\\begin{array}{cc}\n",
    "x _2 - x _1 & y _2 - y _1 \\\\\n",
    "x _3 - x _1 & y _3 - y _1 \\\\\n",
    "\\end{array}\n",
    "\\right)\n",
    "$$\n",
    "Let $T$'s nodes being $n_1$, $n_2$, $n_3$. Denote $\\eta_i := v_{n_i}$.\n",
    "Then it can be shown\n",
    "$$\n",
    "\\nabla \\eta _j = \\frac{1}{2|T|}\n",
    "\\left(\n",
    "\\begin{array}{cc}\n",
    "y _{j+1} - y _{j+2} \\\\\n",
    "x _{j+2} - x _{j+1} \\\\\n",
    "\\end{array}\n",
    "\\right)\n",
    "$$\n",
    "Here we mean $(x_4,y_4) = (x_1,y_1)$, $(x_5,y_5) = (x_2,y_2)$\n",
    "\n",
    "* One can do calculations and derive that the matrix $M_{j,k} = |T| (\\nabla \\eta_j)\\cdot(\\nabla \\eta_k)$ can be evaluated as\n",
    "$$\n",
    "M = \\frac{|T|}{2} G G^T\n",
    "\\qquad\\text{where}\\qquad\n",
    "G =\n",
    "\\left(\n",
    "\\begin{array}{ccc}\n",
    "1 & 1 & 1 \\\\\n",
    "x _1 & x _2 & x _3 \\\\\n",
    "y _1 & y _2 & y _3\n",
    "\\end{array}\n",
    "\\right)^{-1}\n",
    "\\left(\n",
    "\\begin{array}{ccc}\n",
    "0 & 0 \\\\\n",
    "1 & 0 \\\\\n",
    "0 & 1\n",
    "\\end{array}\n",
    "\\right)\n",
    "$$\n",
    "\n",
    "* the pseudo-code would then look like this:\n",
    "```\n",
    "for all T\n",
    "    calculate M\n",
    "    for k=1..3\n",
    "        for l=1..3\n",
    "            A(nodes(k),nodes(l)) += M(k,l)\n",
    "```\n",
    "\n",
    "* More details in the <a href=\"http://www2.mathematik.hu-berlin.de/~cc/cc_homepage/download/1999-AJ_CC_FS-50_Lines_of_Matlab.pdf\">Remarks around 50 lines of Matlab: short finite element\n",
    "implementation</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Assembling the right-hand side (forces)\n",
    "\n",
    "* To evaluate forces, and avoid the exact integration in a triangle\n",
    "$$\n",
    "\\int_T f \\eta_i\n",
    "$$\n",
    "one can approximate\n",
    "$$\n",
    "\\int_T f \\eta_i \\approx f(x_S, y_S) \\int_\\Omega \\eta_i,\n",
    "$$\n",
    "where $(x_S,y_S)$ is the barycenter (i.e., center of mass) of the triangle.\n",
    "The code (in 2D) would the look like\n",
    "```\n",
    "for all T\n",
    "  calculate f(xS, yS)\n",
    "  for k=1..3\n",
    "    f(nodes(k)) += 1/3 * area(T) * f(xS, yS)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### How about boundary conditions?\n",
    "\n",
    "* What do you think?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### How about boundary conditions?\n",
    "\n",
    "* The code should actually look like this:\n",
    "\n",
    "```\n",
    "A = zero matrix for free nodes\n",
    "for k = free nodes\n",
    "    for l = free nodes\n",
    "        for all T\n",
    "            <<SAME>>\n",
    "```\n",
    "This fills only the needed rows & columns of the matrix\n",
    "\n",
    "* The loops in the algorithm are often reversed:\n",
    "```\n",
    "A = zero matrix for all nodes\n",
    "for all T\n",
    "    <<SAME>>\n",
    "remove rows and columns from A corr. to non-free nodes\n",
    "```\n",
    "\n",
    "* The **free nodes** are those that are on $\\Gamma_2$ in\n",
    "$$\n",
    "\\begin{align*}\n",
    "-\\Delta u &= 0\\qquad\\text{on $\\Omega$}\n",
    "\\\\\n",
    "u &= 0\\qquad\\text{on $\\Gamma_1$}\n",
    "\\\\\n",
    "u_n &= 0\\qquad\\text{on $\\Gamma_2$}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "* Alternatively, one can **replace** the corresponding row and column by\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "\\cdot & 0 & \\cdot & \\cdot \\\\\n",
    "0 & 1 & 0 & 0\\\\\n",
    "\\cdot & 0 & \\cdot & \\cdot \\\\\n",
    "\\cdot & 0 & \\cdot & \\cdot \\\\\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "* **Exercise**: think why"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Algorithm\n",
    "\n",
    "* Let's think together what the algorithm would look like (and draft on a whiteboard)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## More advanced stuff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Error analysis\n",
    "\n",
    "* If $u$ is the exact solution and $u_h$ is the approximate solution then one can prove that\n",
    "$$\n",
    "\\|\\nabla u_h - \\nabla u\\|_{L^2}\n",
    "\\leq\n",
    "h \\|\\nabla^2 u\\|_{L^2},\n",
    "$$\n",
    "where $h$ is the **maximal** size of the triangle and $\\nabla^2 u$ is the <a href=\"http://en.wikipedia.org/wiki/Hessian_matrix\">Hessian matrix</a>.\n",
    "\n",
    "* Comes from a truly beautiful argument:\n",
    "    - One can show that $u_h$ is the minimizer of\n",
    "    $$\n",
    "    \\int \\nabla (v_h-u) \\cdot \\nabla (v_h-u) = \\|\\nabla v_h - \\nabla u\\|_{L^2}^2\n",
    "    $$\n",
    "    over all $v_h\\in \\mathcal P_1(\\mathcal T)$ (let's prove it).\n",
    "    - It remains to invoke a geometry-related argument that given $u$\n",
    "    there exists $v_h$ such that $\\|\\nabla v_h - \\nabla u\\| \\leq h \\|\\nabla^2 u\\|$.\n",
    "    This is achieved with, e.g., the Cl&eacute;ment interpotation (warning: google for that on your own risk!)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Mesh refinement / Adaptivity\n",
    "\n",
    "<img src=\"fig/borehole_degen.png\">\n",
    "\n",
    "<img src=\"fig/Lshape_refin.png\">\n",
    "\n",
    "* Often done through the a posteriori error estimate, e.g.,\n",
    "$$\n",
    "\\|\\nabla u_h - \\nabla u\\|^2 \\leq C \\sum_{{\\rm edge}} h_{{\\rm edge}} \\Big[\\frac{\\nabla u_h}{h_{{\\rm edge}}}\\Big]^2_{{\\rm edge}}\n",
    "$$\n",
    "where $[\\bullet]_{\\rm edge}$ is the jump across the edge (recall that $\\nabla u_h$ is piecewise constant)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Mesh quality\n",
    "\n",
    "* Mesh should not have very sharp angles (otherwise the approximation will be bad)!\n",
    "\n",
    "<img src=\"fig/fig_mesh_quality.gif\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Collocation\n",
    "\n",
    "* If we have non-constant coefficients, we use the quadrature rules\n",
    "$$\n",
    "\\int_T k(x,y) \\nabla u\\cdot\\nabla v\n",
    "\\approx \\sum_{i} w_i k(x_i,y_i) \\nabla u(x_i,y_i)\\cdot\\nabla v(x_i,y_i)\n",
    "$$\n",
    "\n",
    "* must be careful to take the quadrature rule of sufficiently high order"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Quad(rilateral) elements:\n",
    "\n",
    "<img src=\"fig/FEM_quad_basis.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### High-order elements\n",
    "\n",
    "* $\\mathcal P^p(\\mathcal T)$: polynomials of power $p$\n",
    "* $hp$-refinement\n",
    "\n",
    "A good source for these is <BR>\n",
    "C. Schwab, \"p- and hp- Finite Element Methods: Theory and Applications to Solid and Fluid Mechanics\", Clarendon Press, 1998."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### More stuff:\n",
    "\n",
    "* Discontinous Galerkin, $\\mathcal P^0$ elements\n",
    "* mixed FEM, generalized FEM, extended FEM, what-not FEM..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Adv/Disadv of FEM\n",
    "\n",
    "* Good for simple physics, solid mechanics\n",
    "* Good energy (and other) conservation properties\n",
    "\n",
    "But...\n",
    "\n",
    "* Not easy for fluid flow (must work hard for stability)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##### Questions?\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}