2020-03-07-fem-heat-equation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finite element method for 1D heat equation\n",
    "\n",
    "A Jupyter Notebook with source code is located at [this](https://github.com/krashkov/Blog) GitHub repository.\n",
    "\n",
    "**Contents**\n",
    "\n",
    "* Classical and weak solutions of PDEs\n",
    "* The basics of FEM\n",
    "    * Galerkin approximations\n",
    "    * Basis functions\n",
    "    * Final system of equations\n",
    "* Numerical example: 1D heat transfer problem\n",
    "    * Analytical solution\n",
    "    * Numerical solution with FEM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.integrate\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Classical and weak solutions of PDEs\n",
    "\n",
    "In this section an intuitive explanation of classical and weak solutions of elliptic second-order PDEs is given. However, for simplicity we will consider an ODE instead on an open bounded domain $\\Omega$ with a mixed condition on the boundary $\\partial \\Omega$:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "& - \\frac{d}{dx} \\left( k(x) \\frac{du(x)}{dx} \\right) + c(x) \\frac{du(x)}{dx}\n",
    "  + b(x) u(x) = f(x),\\ x \\in (0,1) \\\\\n",
    "& \\alpha_0 \\frac{du(0)}{dx} + \\beta_0 u(0) = \\gamma_0\\\\\n",
    "& \\alpha_1 \\frac{du(1)}{dx} + \\beta_1 u(1) = \\gamma_1\\\\\n",
    "\\end{align}\n",
    "\\tag{1}\n",
    "$$\n",
    "\n",
    "Roughly speaking, minimal smoothness solution of a PDE which satisfies both equations and boundary condition is called a *classical solution*. In other words, if a PDE contains $n$-th partial derivative, then this derivative must exist. A classical solution for a general linear strictly elliptic PDE is a function $u$ such that $u \\in C^{2, \\alpha} (\\bar{\\Omega})$, i.e. it is twice differentiable with $\\alpha$*-Holder continuous* second-order derivatives [1].\n",
    "\n",
    "However, many PDEs do not admit of sufficiently smooth solutions. Alternative description which allows to relax the assumptions is called a *weak formulation*. To derive a weak formulation of (1) let us test this equation against any sufficiently smooth function $\\phi(x)$:\n",
    "\n",
    "$$\n",
    "- \\int_0^1 \\frac{d}{dx} \\left( k \\frac{du}{dx} \\right) \\phi dx\n",
    "+ \\int_0^1 c \\frac{du}{dx} \\phi dx\n",
    "+ \\int_0^1 b u \\phi dx\n",
    "= \\int_0^1 f \\phi dx\n",
    "$$\n",
    "\n",
    "Itegrating by parts the first term in the left-hand side of the equation gives the following relation:\n",
    "\n",
    "$$\n",
    "- ku'\\phi|_0^1\n",
    "+ \\int_0^1 (ku'\\phi' + c u' \\phi + b u \\phi ) dx\n",
    "= \\int_0^1 f \\phi dx\n",
    "$$\n",
    "\n",
    "Note that now $u$ de facto is not required to have a second derivative at every single point. But in what class of functions to look for to find a suitable choice for $u$ and $\\phi$ so that it will make sense? It turns out, that a *Sobolev space* $H^1$ would be a reasonable answer. For more information about space of solutions and test functions of second-order elleptic PDEs, please, refer to [1]. Thus a weak formulation of the problem (1) is as follows.\n",
    "\n",
    "**Weak formulation**. Find $u \\in H^1$ such that\n",
    "\n",
    "$$\n",
    "- ku'\\phi|_0^1\n",
    "+ \\int_0^1 (ku'\\phi' + c u' \\phi + b u \\phi ) dx\n",
    "= \\int_0^1 f \\phi dx,\\ \n",
    "\\forall \\phi \\in H^1\n",
    "\\tag{2}\n",
    "$$\n",
    "\n",
    "A function $u \\in H^1$ that satisfies (2) is called a *weak solution*.\n",
    "\n",
    "In addition, in the new formulation of the problem it is also not required of $k(x)$ to be continuous since the existence of a second derivative is not implied. In order to obtain compatibility condition, let us assume that $k(x)$ is discontinuous at some $x=b$. The following holds:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "0\n",
    "&= - ku'\\phi|_0^b + \\int_0^b (ku'\\phi' + c u' \\phi + b u \\phi - f \\phi) dx\n",
    "   - ku'\\phi|_b^1 + \\int_b^0 (ku'\\phi' + c u' \\phi + b u \\phi - f \\phi) dx\\\\\n",
    "&= - ku'\\phi|_{x=b-0} + ku'\\phi|_{x=b+0}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Since $\\phi$ is arbitrary:\n",
    "\n",
    "$$\n",
    "ku'|_{x=b-0} = ku'|_{x=b+0}\n",
    "\\tag{3}\n",
    "$$\n",
    "\n",
    "Discontinuity of the second kind of $f$ leads to the same result (3), but discontinuity of the third kind result in a jump of a *flux* $W = -ku'$ at this point [2]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 The basics of FEM\n",
    "\n",
    "Finite element method is a numerical approach used to find an approximate solution of a weak formulation of PDEs or ODEs. Thus, it is possible to deal with poorly behaved functions which are quite common in practical applications.\n",
    "\n",
    "However, a weak formulation of the problem does not seem to be any simpler since a space of possible solutions $u$ is supposed to be infinite-dimensional. So the question is: how to approximate a weak formulation in such a way that it will be well-posed?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Galerkin approximations\n",
    "\n",
    "The idea behind approximation of a weak formulation is to replace an infinite-dimensional space of solution and test functions with a finite-dimensional subspace. Let us denote this subspace as $S^h \\subset H^1$, $\\dim S^h = n$, and select a convinient basis $\\{\\psi_0(x), \\cdots, \\psi_{n-1}(x)\\}$ in $S^h$. Decomposition on the basis functions is as follows:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\phi_h(x) &= \\sum_{i=0}^{n-1} \\zeta_i \\psi_i(x)\\\\\n",
    "u_h(x)    &= \\sum_{i=0}^{n-1} u_i \\psi_i(x)\n",
    "\\end{align}\n",
    "\\tag{4}\n",
    "$$\n",
    "\n",
    "So the goal is to find a set of coefficients $(u_0, \\cdots, u_{n-1})$ since they fully define a finite element approximation $u_h \\in S^h$. Substitution of $u$ and $\\phi$ in Eq. (2) with expressions in Eq. (4) gives:\n",
    "\n",
    "$$\n",
    "\\sum_{i=0}^{n-1} \\zeta_i \\left[\n",
    "  - k(1) u'(1) \\psi_i(1)\n",
    "  + k(0) u'(0) \\psi_i(0)\n",
    "  + \\sum_{j=0}^{n-1} u_j \\int_0^1 \\left( k \\psi_j' \\psi_i' + c \\psi_j' \\psi_i + b  \\psi_j \\psi_i \\right)dx\n",
    "  - \\int_0^1 f \\psi_i dx \n",
    "\\right] = 0\\\\\n",
    "\\sum_{i=0}^{n-1} \\zeta_i \\left[\n",
    "    k(1) \\frac{\\gamma_1 - \\beta_1 u(1)}{\\alpha_1} \\psi_i(1)\n",
    "  - k(0) \\frac{\\gamma_0 - \\beta_0 u(0)}{\\alpha_0} \\psi_i(0)\n",
    "  + \\sum_{j=0}^{n-1} K_{ij} u_j - F_i\n",
    "\\right] = 0\n",
    "$$\n",
    "\n",
    "where the following notation is utilised $K_{ij} = \\int_0^1 \\left( k \\psi_j' \\psi_i' + c \\psi_j' \\psi_i + b  \\psi_j \\psi_i \\right)dx$ and $F_i = \\int_0^1 f \\psi_i dx$. $K$ is called *stiffness matrix* and $F$ is a *load vector*. Since $\\zeta_i$ are arbitrary the following equation holds:\n",
    "\n",
    "$$\n",
    "  k(1) \\frac{\\gamma_1 - \\beta_1 u(1)}{\\alpha_1} \\psi_i(1)\n",
    "- k(0) \\frac{\\gamma_0 - \\beta_0 u(0)}{\\alpha_0} \\psi_i(0)\n",
    "+ \\sum_{j=0}^{n-1} K_{ij} u_j\n",
    "= F_i,\\ \\forall i \\in \\{0, \\cdots, n-1\\}\n",
    "\\tag{5}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Basis functions\n",
    "\n",
    "In fact, there are different ways to define basis functions. Here we will consider a concrete example and use it later for a numerical problem in section 3. Let us introduce a mesh $\\mathcal{T}_h$:\n",
    "\n",
    "$$\n",
    "\\mathcal{T}_h = \\{x_i = ih,\\ h = \\frac{1}{n - 1},\\ i \\in \\{0,\\cdots, n-1\\}\\}\n",
    "$$\n",
    "\n",
    "Each element of length $h_k$ is called a *finite element*. Basis functions are defined over finite elements in such a way that $u_h(x_i) = u_i,\\ \\forall i$. Hence:\n",
    "\n",
    "$$\n",
    "u(0) = u_0,\\ u(1) = u_{n-1}\n",
    "\\tag{6}\n",
    "$$\n",
    "\n",
    "Now let us introduce a set of basis functions on $[0,1]$ also known as \"hats\":\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "&\\Phi(t) =\n",
    "  \\begin{cases}\n",
    "    1 - |t|,\\ &|t| \\leq 1 \\\\ 0,\\ &|t| > 1\n",
    "  \\end{cases}\\\\\n",
    "&\\psi_i(x) = \n",
    "  \\Phi\\left(\\frac{x - x_i}{h}\\right),\\ \\forall i \\in \\{0, \\cdots, n-1\\}\n",
    "\\end{align}\n",
    "\\tag{7}\n",
    "$$\n",
    "\n",
    "It is possible to verify that:\n",
    "\n",
    "$$\n",
    "\\psi_j(x_i) = \\delta_{ij}\\ \\Rightarrow\\ u_h(x_i) = u_i\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Final system of equations\n",
    "\n",
    "Taking into consideration all above mentioned, we get a system of $(n-1)$ linear equations with $(n-1)$ unknowns.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "i = 0:          &\\ - k_0 \\frac{\\beta_0 u_0}{\\alpha_0} + \\sum_{j=0}^{n-1} K_{0j} u_j \n",
    "                   = F_0 + k_0\\frac{\\gamma_0}{\\alpha_0}\\\\\n",
    "i\\neq\\{0, n-1\\}:&\\ \\sum_{j=0}^{n-1} K_{ij} u_j = F_i\\\\\n",
    "i = n-1:        &\\ k_{n-1} \\frac{\\beta_1 u_{n-1}}{\\alpha_1} + \\sum_{j=0}^{n-1} K_{n-1,j} u_j\n",
    "                   = F_{n-1} - k_{n-1}\\frac{\\gamma_1}{\\alpha_1}\\\\\n",
    "\\end{align}\n",
    "\\tag{8}\n",
    "$$\n",
    "\n",
    "Each basis function (7) intersects only with $\\psi_{i-1}, \\psi_{i+1}$. Thus (8) is in fact tridiagonal system of linear equations:\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "K_{00} - \\frac{k_0\\beta_0}{\\alpha_0} & K_{01} & 0      & \\cdots & 0\\\\\n",
    "K_{10} & K_{11} & K_{12} & \\cdots & 0\\\\\n",
    "0      & K_{21} & K_{22} & \\cdots & 0\\\\\n",
    "\\cdots & \\cdots & \\cdots & \\cdots & \\cdots\\\\\n",
    "0      & 0      & 0      & \\cdots & K_{n-1,n-1} + \\frac{k_{n-1}\\beta_1}{\\alpha_1}\\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "u_0\\\\\n",
    "u_1\\\\\n",
    "u_2\\\\\n",
    "\\cdots\\\\\n",
    "u_{n-1}\n",
    "\\end{bmatrix}\n",
    "= \\begin{bmatrix}\n",
    "F_0 + k_0\\frac{\\gamma_0}{\\alpha_0}\\\\\n",
    "F_1\\\\\n",
    "F_2\\\\\n",
    "\\cdots\\\\\n",
    "F_{n-1} - k_{n-1}\\frac{\\gamma_1}{\\alpha_1}\n",
    "\\end{bmatrix}\n",
    "\\tag{9}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 Numerical example: 1D heat transfer problem\n",
    "\n",
    "Consider the following heat equation:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "& -\\frac{d}{dx}\\left(k\\frac{du}{dx}\\right) = f,\\ x \\in (0,1)\\\\\n",
    "& u'(0) = 0,\\ u(1) = \\gamma_1\n",
    "\\end{align}\n",
    "\\tag{10}\n",
    "$$\n",
    "\n",
    "where $u$ represents temperature and $f$ is power distribution of heat sources inside the material. Also suppose that:\n",
    "\n",
    "$$\n",
    "f(x) =\n",
    "\\begin{cases}\n",
    "f_1, &x < 0.5 \\\\ f_2, &x \\geq 0.5\n",
    "\\end{cases}\n",
    "\\qquad\n",
    "k(x) =\n",
    "\\begin{cases}\n",
    "k_1, &x < 0.5 \\\\ k_2, &x \\geq 0.5\n",
    "\\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    if (x < 0.5): return 2.\n",
    "    else:         return 0.\n",
    "    \n",
    "def k(x):\n",
    "    if (x < 0.5): return 0.01\n",
    "    else:         return 0.04\n",
    "    \n",
    "gamma1 = 273.\n",
    "\n",
    "x0 = 0.\n",
    "x1 = 1.\n",
    "n  = 5\n",
    "\n",
    "mesh = np.linspace(x0, x1, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Analytical solution\n",
    "\n",
    "On each domain $[0, 0.5)$ and $[0.5, 1]$ Eq. (10) is simply an ODE with constant coefficients. \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "x < 0.5:&   \\quad u(x) = c_1 + c_2 x - \\frac{f_1 x^2}{2 k_1}\\\\\n",
    "x \\geq 0.5:&\\quad u(x) = c_3 + c_4 x - \\frac{f_2 x^2}{2 k_2}\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Since $u$ is de facto temperature, it should be continuous at $x = 0.5$, so apart from boundary conditions the following equations should be used to find unknown constants:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "u(x)|_{x=0.5-0} &= u(x)|_{x=0.5+0}\\\\\n",
    "k_1 u'(x)|_{x=0.5-0} &= k_2 u'(x)|_{x=0.5+0}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "One may obtain that:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "c_1 &= \\frac{f_1k_1+0.5f_2k_1+0.5f_1k_2+4.k_1k_2\\gamma_1}{4k_1k_2},\\ &c_2 &= 0,\\\\\n",
    "c_3 &= \\frac{f_1+2k_2\\gamma_1}{2k_2},\\ &c_4 &= -\\frac{f_1-f_2}{2k_2}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exact_solver(x, f1, f2, k1, k2, gamma1):\n",
    "    c1 =   (f1*k1 + 0.5*f2*k1 + 0.5*f1*k2 + 4.*k1*k2*gamma1)/(4.*k1*k2)\n",
    "    c2 =   0.\n",
    "    c3 =   (f1 + 2.*k2*gamma1)/(2.*k2)\n",
    "    c4 = - (f1 - f2)/(2.*k2)\n",
    "    \n",
    "    if (x < 0.5): return c1 + c2*x - f1*x*x/(2.*k1)\n",
    "    else:         return c3 + c4*x - f2*x*x/(2.*k2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "exact_solution = np.array([\n",
    "    exact_solver(x, f(0), f(1), k(0), k(1), gamma1) for x in np.linspace(x0, x1, 10**3)\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Numerical solution with FEM\n",
    "\n",
    "As a basis we will use \"hat\" functions defined by (7)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def psi(x, i, mesh):\n",
    "    h = mesh[1] - mesh[0]\n",
    "    t = (x - mesh[i])/h\n",
    "    if (abs(t) <= 1): return (1. - abs(t))\n",
    "    else:             return 0.\n",
    "\n",
    "def psi_prime(x, i, mesh):\n",
    "    h = mesh[1] - mesh[0]\n",
    "    t = (x - mesh[i])/h\n",
    "    if (abs(t) <= 1.):\n",
    "        if (t <= 0.): return  1./h\n",
    "        else:         return -1./h\n",
    "    else: return 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to simplify (9) even further using boundary conditions from (10). Since $u_{n-1} = \\gamma_1$ there are only $n-2$ equations:\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "K_{00} & K_{01} & 0      & \\cdots & 0\\\\\n",
    "K_{10} & K_{11} & K_{12} & \\cdots & 0\\\\\n",
    "0      & K_{21} & K_{22} & \\cdots & 0\\\\\n",
    "\\cdots & \\cdots & \\cdots & \\cdots & \\cdots\\\\\n",
    "0      & 0      & 0      & \\cdots & K_{n-2,n-2}\\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "u_0\\\\\n",
    "u_1\\\\\n",
    "u_2\\\\\n",
    "\\cdots\\\\\n",
    "u_{n-2}\n",
    "\\end{bmatrix}\n",
    "= \\begin{bmatrix}\n",
    "F_0\\\\\n",
    "F_1\\\\\n",
    "F_2\\\\\n",
    "\\cdots\\\\\n",
    "F_{n-2} - K_{n-2,n-1}\\gamma_1\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Certainly, this system should be solved with the tridiagonal matrix algorithm. However, in principle, different possibilities for $f$, $k$ as well as basis functions are possible, so here a general implementation is presents, i.e. each element of the stiffness matrix and load vector is constructed by integrating corresponding functions over finite elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fem_solver(gamma1, f, k, x0, x1, n, mesh, psi, psi_prime, exact_solver):\n",
    "    \n",
    "    # Construct K\n",
    "    stiffness_matrix = np.zeros(shape=(n-1,n-1))\n",
    "    for i in range(n-1):\n",
    "        for j in range(n-1):\n",
    "            stiffness_matrix[i][j] = scipy.integrate.quad(\n",
    "                lambda x: k(x)*psi_prime(x, i, mesh)*psi_prime(x, j, mesh),\n",
    "                x0, x1, epsabs = 1e-4, limit=200, points = [mesh[i-1], mesh[i], mesh[i+1]]\n",
    "            )[0]\n",
    "        \n",
    "    # Construct F\n",
    "    load_vector = np.zeros(shape=(n-1))\n",
    "    for i in range(n-1):\n",
    "        load_vector[i] = scipy.integrate.quad(\n",
    "            lambda x: f(x)*psi(x, i, mesh),\n",
    "            x0, x1, epsabs = 1e-4, limit=200, points = [mesh[i-1], mesh[i], mesh[i+1]]\n",
    "        )[0]\n",
    "    load_vector[n-2] = load_vector[n-2] - gamma1 * scipy.integrate.quad(\n",
    "        lambda x: k(x)*psi_prime(x, n-1, mesh)*psi_prime(x, n-2, mesh),\n",
    "        x0, x1, epsabs = 1e-4, limit=200, points = [mesh[n-3], mesh[n-2], mesh[n-1]]\n",
    "    )[0]\n",
    "        \n",
    "    # Solve linear system\n",
    "    res = np.linalg.solve(stiffness_matrix, load_vector)\n",
    "    res = np.append(res, gamma1)\n",
    "    \n",
    "    # Find exact solution for a given mesh\n",
    "    exact_solution = np.array([exact_solver(x, f(0), f(1), k(0), k(1), gamma1) for x in mesh])\n",
    "    \n",
    "    return {'res': res, 'tol': np.linalg.norm(res - exact_solution)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fem_solution = fem_solver(gamma1, f, k, x0, x1, n, mesh, psi, psi_prime, exact_solver)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us compare exact and numerical solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEJCAYAAACHRBAhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd3hUxdfA8e9JSAiEkERKQEMISJEqSChKEQWkSO/Si6I0wYYdu9gR9PfSpIgivSNFRJSiIFWkKL0XEUgIhBCSnPePXRCRsgubbMr5PM99cnduOzeEnMzMvTOiqhhjjDE34uPtAIwxxqQPljCMMca4xBKGMcYYl1jCMMYY4xJLGMYYY1ySxdsBpJTcuXNrZGSkt8Mwxph0Zd26dX+rap6rbcuwCSMyMpK1a9d6OwxjjElXRGTftbZZk5QxxhiXWMIwxhjjEksYxhhjXGIJwxhjjEsybKe3McYAJCcnc/DgQc6ePevtUNKMwMBAwsPD8fFxr85gCcMYk6H9/fffiAjFixd3+xdkRpScnMyhQ4f4+++/yZs3r1vH2nfPGJOhRUdHExYWZsnCycfHh7CwMGJiYtw+1moYV/hr71988tInRBWMIn9ofnyz+pIlaxb8Av3Injv7v5asQVm9Ha4x5gaSkpLw8/Pzdhhpip+fH4mJiW4f55WEISIBwDIgqzOGaar6moj0AfoDdwJ5VPVv5/4CDAEaAHFAF1VdnxKxrfltDQGTAtiiW9jCluvuGxASQEihEEILhRJSKITbitxGWNkw8pbJa8nEmDTE8SvEXHSz3w9v1TDOAw+q6hkR8QNWiMgCYCUwD/jxiv3rA0WdS2VgmPOrxz3c5GEqxVZi1K+jGLlqJH9F/0WRoCJ0Ld6V+nnrI6eFuL/jiDseR/S+aKL3RHN863G2f7udpPNJl84TWjiUsLvDuKPSHURUi+D2qNvJEmAVOmNM+uWV32DqmObvjPOjn3NRVd0AV81+TYDxzuNWiUiIiORX1SMpEV+ewDy89MBLPFfjOWZsm8GQ1UN4euvTvLbzNbqV70bfJn0pf1v5f99TshJzIIZjm45x7LdjHNt0jKMbj/LHzD8A8PX35faKtxNRPYI7H7qTiKoR+Pr7pkT4xhiTIsRbU7SKiC+wDigC/E9Vn79s214g6rImqXnAe6q6wvl5CfC8qq694pw9gB4AERERFfbtu+aQKG5bc2gNQ1YPYcqWKSQmJ/JwsYfpV7kftQrVum71Lu7vOA78fID9K/azf8V+Dq85THJiMv5B/hSuXZgi9YtQ7OFiBN0e5LFYjTH/2LZtGyVKlPB2GGnOtb4vIrJOVaOudozXEsalAERCgJlAX1Xd7Czby78TxrfAoCsSxgBVXXet80ZFRWlKDD54JPYIw9cOZ/i64fx19i9K5SnFk5WfpEPZDmT3y37D48/HnmfPkj3sWLCDnQt2cvrAaRAoWL0gpdqUokSLEuQIy+HxuI3JrCxhXN3NJAyvP2emqtE4+izqXWe3g0CByz6HA4dTMKxryh+UnzceeIN9/fcxrsk4/H39eXze44R/Es7zi59nf8z+6x6fNSgrdzW9i0YjGtF/X396/t6T+1+7n7PHzzK/93w+uf0Txtcaz29f/caFuAupdFfGGG+IjIzko48+omzZsgQHB9OmTRvi4+MZN24c1apV+9e+IsLOnTsB6NKlC7169aJ+/frkyJGDqlWrcvToUfr3709oaCh33XUXGzZs8Hi83npKKg9wQVWjRSQbUBt4/zqHzAH6iMgkHJ3dMSnVf+GqgCwBdC7XmU53d2LF/hUM/XUoH/3yER//8jHNSjSjX+V+VC1Q9brNVSJC3tJ5yVs6L/cPvJ/jW46zefJmNk/czKxOs1jQdwFl2pfhnkfvIX/5/Kl4d8ZkXP0X9mfj0Y0peo1y+crxab1PXdp3ypQpLFy4kICAAKpWrcq4ceMICAhw6bhFixZRqlQpGjRowL333ssbb7zBxx9/zGuvvcbTTz/N0qVLb/VW/sVbNYz8wFIR2QSsARar6jwReVJEDuKoQWwSkS+c+88HdgM7gVFAL28EfTUiQvWC1Znaaip7+u3hmXufYcnuJVQfW52oUVGM/2085xPPu3SevKXz8uBbD9J3e186L+1MsYbF2DhmIyPvGcmoiqPYNGETSQlJNzyXMSb9ePLJJ7n99tu57bbbaNSoERs3upbMmjVrRoUKFQgICKBZs2YEBATQqVMnfH19adOmTcapYajqJqD8VcqHAkOvUq5A71QI7ZZEBEfwfp33GXj/QL7e9DVDfx1K51mdeW7xczxR4Ql6VuxJvhz5bnge8REia0YSWTOSc5+d4/cJv7Pmf2uY2WEm3w/4nop9KhL1eBTZbsuWCndlTMbi6l/+qSVfvn9+J2TPnp3Dh11rbQ8LC7u0ni1btv98PnPmzNUOuyVe78PIiAL9A3k86nE299zMdx2+o+LtFXlz2ZtEDI6g48yOrD3semd8ttBsVOpTiV5betHu23bkKZmHH176gcEFBrOw/0JiD8em4J0YY7whMDCQuLi4S5+PHj3qxWj+YQkjBYkIde6sw7x289jeZzs9o3oy+4/ZVBxVkapjqjJlyxQuJLnWsS0+QtEGRem4uCNPbHqCkq1K8uvnvzKk8BDm953P6YOnU/hujDGp5e6772bLli1s3LiR+Ph4Xn/9dW+HBFjCSDVFcxVlSP0hHHz6IJ/W/ZRjZ47RZlobCg8tzKDlgzgRd8Llc4WVCaPpuKb03d6Xsh3Ksm74OobeOZRve3/LmaOer4YaY1JXsWLFGDhwILVr16Zo0aL/eWLKW7z+HkZKSan3MDwlKTmJ+TvmM2T1EJbsWUJAlgA6lOnAk5WfpExYGbfOFb03muWDlrNxzEZ8s/py33P3cd8z9+Gfwz+Fojcm/bD3MK4uXb6HkVn5+vjSqHgjvu/0PZt7bqZT2U5M+H0CZYeXpdb4Wsz+YzZJya49ERUSGUKjEY3otbUXResX5afXf2JokaGsHb6W5MTkFL4TY0xmYQkjDSiVtxQjGo3gwFMHeK/We+w4sYOmk5tS7PNiDP5lMDHxro1bn6toLlpNbUX3X7qTq2guvu35LSPKj2DfMs8NkWKMybwsYaQhubLn4vlqz7O7326mtJzC7UG38/R3TxM+OJy+8/uy/cR2l84TXiWcLsu60HpGa87Hnmfc/eOY0WEGsUfsiSpjzM2zhJEGZfHJQqtSrVjedTlrH1tL8xLNGbl+JMU/L06DCQ1YtHMRN+p7EhFKNCtB7629qfFqDbZO3crnxT/nl09+IemCvfxnjHGfJYw0rsLtFfiy6Zfs77+f1+9/nfVH1lNvQj1K/l9Jhq0ZxtmE609s75fdjwfefIBeW3oRUS2C7575ji8qf8HRjWnjuW5jTPphCSOdCMsRxms1X2Nf/32MbzqeQL9Aes3vRfjgcJ797ln2Ru+97vG3FbmNdt+2o9W0VsQeimVUxVEsHbiUxPPuT9NojMmcLGGkM1mzZKXj3R1Z89gaVnZbyUN3PsSnqz7lzqF30nxyc37a+9M1m6tEhJItStJray9KP1KaZW8tY2SFkRz69VAq34UxJj2yhJFOiQj3FbiPyS0ns6ffHgbcN4Cf9v1EzS9rUn5EecZuGEt8YvxVj82eKzvNxjej3bftOB9zntH3juaHV36wvg1jzHVZwsgACgQXYFDtQRx46gCjGo0iSZPoNqcbBQYX4JUfXuFw7NUHMyvaoCg9N/fk7k53s/yd5YytNpaTO0+mcvTGmPTCEkYGkt0vO4/e8yibntjEkk5LuK/Afby7/F0KflqQdtPbsfrg6v8cExAcQJOxTWg5pSUndpxgeLnhbBiz4YZPYRlj0qbIyEi+//77FDm3JYwMSER4sNCDzG47mx19d9CnYh++3fEtVUZXocoXVZj4+8T/DHpYqlUpem7qyR2V7mBO9zlMbTWV+OirN2kZYzInSxgZ3J233cngeoM5+NRBhtYbyslzJ2k3ox2RQyJ5e9nbHD97/NK+OcNz0un7TtR+vzZ/zv6TkRVGcmSDVyc2NCbDO3z4MC1atCBPnjwUKlSIoUMdUwI1aNCAZ5555tJ+bdq0oVu3bpc+jxo1ihIlShAUFETJkiVZv349HTt2ZP/+/TRq1IgcOXLwwQcfeDRWG3wwk0nWZBbuXMiQ1UP4btd3ZPXNSrsy7ehXuR9357v70n4HfjnAtNbTOHv8LA0+b0D57uWvO92sMWnVlYPsLey/MMXfQ8pXLh/1Pq13w/2Sk5OpWLEiTZo04YUXXuDgwYPUrl2bYcOGcffdd1O2bFkmTZrEkSNHePXVV/ntt98ICgpi6tSp9O/fn1mzZhEVFcWuXbvw8/OjYMGCREZG8sUXX1C7du3rXjvdDD4oIgEi8quI/CYiW0TkDWd5IRFZLSI7RGSyiPg7y7uIyHER2ehcHvVG3BmBj/jQoGgDFnVYxNZeW+lariuTt0ym3Ihy1BxXkxnbZpCUnESBewvQY30PCtYoyNzH5jK762wuxLk2d4cxxjVr1qzh+PHjDBw4EH9/fwoXLsxjjz3GpEmTyJcvH8OHD6dz587069eP8ePHExQUBMAXX3zBgAEDqFixIiJCkSJFKFiwYIrH65UpWoHzwIOqekZE/IAVIrIAeBoYrKqTRGQ40B0Y5jxmsqr28VK8GVKJPCUY1nAY79Z6l9EbRvP5r5/TYkoLCgYXpE+lPnQv3532C9qz7K1l/PTmTxxZf4S2s9sSWijU26Ebc9Nc+cs/tezbt4/Dhw8TEhJyqSwpKYnq1asD0LBhQ/r06UPx4sX/NSfGgQMHuPPOO1M9Xq/UMNTh4kw/fs5FgQeBac7yL4GmXggv0wnNFsqz9z3Lzid3Mr31dAqGFOS5xc8RPjicPgv7kK9PPtovaM/pg6cZVXEUe3/c6+2QjckQChQoQKFChYiOjr60xMbGMn/+fABefvllSpQowZEjR5g4ceK/jtu1a9dVz5mSTcde6/QWEV8R2Qj8BSwGdgHRqnpxrIqDwB2XHdJCRDaJyDQRKXCNc/YQkbUisvb48eNX28VcRxafLDQv0ZyfuvzE+h7raV2qNaM3jKbE/0rQ+3hvik4uSmDeQL6q8xVrhq3xdrjGpHuVKlUiZ86cvP/++5w7d46kpCQ2b97MmjVrWLZsGWPHjmX8+PGMHz+evn37cuiQY1SGRx99lI8++oh169ahquzcuZN9+xzTGISFhbF79+6UCVhVvboAIcBSoDqw87LyAsDvzvVcQFbn+hPADzc6b4UKFdTcumNnjulbP72l+T/Kr7yOlvqglL5z3zv6Oq/r3CfmamJCordDNOa6tm7d6u0QruvQoUPatm1bDQsL05CQEK1cubLOnj1bCxYsqBMnTry034ABA7ROnTqanJysqqrDhg3TYsWKaWBgoJYqVUrXr1+vqqqzZs3SAgUKaHBwsH744YfXvO61vi/AWr3G79U08ZSUiLwGxAHPA/lUNVFE7gVeV9W6V+zrC5xU1eDrndOekvKshKQEpm2dxpDVQ1hzYA0NfmxAxWUVCasWRuc5nckWms3bIRpzVTZF69Wlp6ek8ohIiHM9G1Ab2IajptHSuVtnYLZzn/yXHd7Yua9JRf6+/rQr047Vj67m58d+JujJIGY3n83hnw/zRpk32LbZ/kmMyei81YeRH1gqIpuANcBiVZ2Ho4bxtIjsxNEMNdq5/5POx29/A54EunghZuNUJbwKE1tMZO7YueiHipwQxlcdz741NhWsMRlZmmiSSgnWJJV6Js2dxK8dfiXoQhAdZ3akSN0i3g7JmEusSerq0k2TlMlY2jZqi98IP44FH2PCwxPYMGaDt0MyxqQASxjGI95p/Q7bB25nd+HdzOk+hxXvr/B2SMZcklFbUm7WzX4/LGEYj8jik4UJHSawuPti9kbtZckLS/j+xe/tP6rxuoCAAE6cOGE/i06qyokTJwgICHD7WG8NDWIyoPxB+ZnYZiK1z9WmZ2hPVr63kvjoeB7+38OIjw1caLwjPDycgwcPYi/z/iMgIIDw8HC3j7OEYTzq/sj7eafOOzzP8wzKN4h1w9eRcDqBJuOa4Ovn6+3wTCbk5+dHoUKFvB1GhmBNUsbjnrvvOZrc1YRXi75KkeeL8Ps3vzOlxRQS4xNvfLAxJs2yhGE8TkQY13QcBYMLMiDvAKoPrs72eduZ1HSSJQ1j0jFLGCZFhASEMK31NE6eO8m7t71Lw5EN2fXdLiY1mcSFczavhjHpkSUMk2LK5SvH/zX4P5bsWcLMwjNpPLoxuxZb0jAmvbKEYVJU1/Jd6V6+O+8sf4dDVQ/RZEwTdn+/25KGMemQJQyT4j6r/xnl8pWj48yOBDcJ/idpNJ5k074ak464nDBEJExEOorIRyLyhfNrRxHJl5IBmvQvm182preejqrScmpL7upwF03GNmH3kt2Op6fOW0e4MenBDROGiJQQkWnAVqAjjulUjzq/dgS2OGfBK5mikZp0rXBoYcY3G8/6I+vpt6Af5TqXo9GoRuxcuJMZ7WeQnJjs7RCNMTfgyot744APgfaqev7KjSLiDzTBMRT5vR6NzmQojYs35vmqz/P+yvepGlGVTt07kRCbwKKnFjGn+xyajG1ib4Qbk4bdMGGoauWL6yLiq6pJV2xPAKY6F2Ou6+0H32b1odU8Me8JyucrT5X+VTgfe54fB/6If5A/9T+rn6KT2Btjbp67nd4LRSRHikRiMoUsPlmY2GIiIQEhtJjSgpj4GGq8UoN7n72XNf9bw5KXlng7RGPMNbibMDYCK0Xk9osFIlJDRJZ7NiyTkeXLkY/JLSez+9Ruus3pBkCdD+pQ4YkKrHxvJSves6HRjUmL3EoYqvocMAxH0mgrIouBKcAsd84jIgEi8quI/OacevUNZ3khEVktIjtEZLKzfwQRyer8vNO5PdKd65m0p3rB6rxX+z1mbJvB4FWDEREe/t/DlGlXhiUvLmHjuI3eDtEYc4WbeQ/jZ+A0MAE4ABRS1Y/dPMd54EFVvRsoB9QTkSrA+8BgVS0KnAK6O/fvDpxS1SLAYOd+Jp175t5naHZXMwYsHsCK/SsQH6HJ2CYUrl2YOY/OYceCHd4O0RhzGbcShojMAH4EpgPNgbrAg+5eVB3OOD/6ORd1nmuas/xLoKlzvYnzM87ttcR6RtM9EWFsk7EUCi1E66mtOXbmGL7+vrSe3pqwsmFMbTmVQ2sOeTtMY4yTuzWMHUARVX1TVWcD9YDPRaS3uxcWEV8R2Qj8BSwGdgHRqnrxLa6DwB3O9Ttw1GZwbo8Bcrl7TZP2BAcEM731dKLjo3lk+iMkJieSNWdW2s9vT2BYIN88/A0ndpzwdpjGGNzvw3heVU9e9vl3oBrwmLsXVtUkVS0HhAOVgBJX28359Wq1if/MtygiPURkrYistdm10o+yYWUZ9vAwlu5dyqs/vApAjnw56LCwAyh8Xfdrzhw7c4OzGGNS2i2PJaWqh4Dqt3B8NI5mripAiIhcfDckHDjsXD8IFABwbg8GTnIFVR2pqlGqGpUnT56bDcl4QedynXnsnsd4b+V7zPlzDgC5iuXikXmPcPbYWb5p8A0JZxK8HKUxmZsrQ4M8KSJZb7Bbgog86epFRSSPiIQ417MBtYFtwFKgpXO3zsBs5/oc52ec239Qm9E9wxlafyj35L+HTjM7sfvUbgDCK4fTamorjm48yvR200lOsiFEjPEWV2oY+YCdIjJCRNqJSAURKeb8+oiIjMDRt5HXjevmB5aKyCZgDbBYVecBzwNPi8hOHH0Uo537jwZyOcufBl5w41omnQjIEsC0VtMQEVpOaUl8YjwARRsUpd6Qemyfu53FAxZ7OUpjMi9x5Q91EckNdAHqA2WAEByPvW4C5gPjVTVN9UxGRUXp2rVrvR2GuQnzts+j0cRGPFr+UUY1HnWpfH7f+az5fA0NRzSkQo8KXozQmIxLRNapatTVtrky+CCq+jfwkXMxJkU1LNaQF6u9yKAVg6gaUZUu5boAUG9wPU7tOsW3vb4ltHAohWsX9m6gxmQyNoGSSZPefOBNHoh8gJ7f9uS3o78B4JPFh5aTWpKnRB6mtJzC8W32JJwxqcndF/fOiMgqERkpIr1FpLqIBKdUcCbzujhIYWhAKC2ntiQmPgaArDmz8si8R8iSNQsTG04k7u84L0dqTObhbg2jOPAmsAfHo7TTgJMiss/TgRkTliOMKa2msOfUHrrM7sLF/raQgiG0ndOW2MOxTGkxhaQLSTc4kzHGE9x9ce+Qqs5X1UGq2hbHG9hDAHfHkjLGJdUiqvFBnQ+Y9ccsPv7lnx+z8MrhNB7dmH3L9rHoqUVejNCYzOOW+jCcw3Q8AzTyTDjG/NdTVZ6iRYkWvPD9Cyzbt+xSeZl2Zbj3Gcc8GutHr/dihMZkDu72YXwiIp1FpPzFoced5yjo+dCMcRARxjQZw5233UmbaW04eubopW2136tN4TqFmd9rPgdXHfRilMZkfO7WMOKAVsBc4IyIbAF24nhL25gUkzNrTqa1mkZMfAxtprUhMdkxRuXFJ6dyhudkcvPJxB6J9XKkxmRc7vZhvKKqDVU1HLgd6AV0xTHUuTEpqkxYGUY0HMGyfct4ecnLl8qz3ZaNNrPacD7mPFNaTCHxfOJ1zmKMuVk33Yehqn+r6k+q+qOq2mMqJlV0vLsjj1d4nA9+/oDZf8y+VB5WJoymXzbl4C8Hmd9nPjbUmDGeZy/umXTn03qfUiF/BTrP6syuk7sulZdsWZLqL1dnwxcb2DBmgxcjNCZjsoRh0p2ALAFMaz0NH/GhxZQWnLtw7tK2mm/UdHSC957PkQ1HvBilMRmPJQyTLkWGRPJVs6/47dhv9Jnf51K5j68PzSc0J3vu7ExtOZX46HgvRmlMxnLLCUNExohINxHx9URAxrjq4WIP83L1lxmzcQxjNoy5VB6YJ5BWU1oRsz+GWV1mWX+GMR7iiRqGAO1wDHVuTKp6o+Yb1CpUi97ze7Px6MZL5QXuK0CdD+vw5+w/+fmjn70YoTEZhyemaO2qqrWBch6Ixxi3+Pr48k2Lb8iVLRctprQgOj760rbK/SpTsmVJlry4hH3LbLgzY26Vx/owVPWCp85ljDvyBuZlSqsp7I/ZT5dZ/wxSKCI0Ht2Y0MKhTGszjTNHz3g5UmPSN7cThojUEZHRIjLX+TlKRB70fGjGuO6+AvfxUZ2PmP3nbD78+cNL5VlzZqX19NbEx8Qzo/0MmxPcmFvg7lhSfYFhOObwruEsPge87eZ5CojIUhHZJiJbRKSfs/xuEflFRH4XkbkiktNZHiki50Rko3MZ7s71TObwZOUnaV2qNS8ueZEf9/54qTysTBj1P6vPnh/2sPL9ld4L0Jh0zt0aRn+gtqq+B1z8U+0PHPNkuCMReEZVSwBVgN4iUhL4AnhBVcsAM4HnLjtml6qWcy5PuHk9kwmICF80+oKitxWl7bS2HIn95z2M8t3KU6pNKZYOXMqBnw94MUpj0i93E0YQcPF/28VnFf2ABHdOoqpHVHW9cz0Wx+CFd+BIPBfHr14MtHAzPpPJBWUNYnrr6cQmxNJmWhsuJDm61kSEhiMaEhwRzPR20+39DGNugrsJYznwwhVlTwJLbzYAEYkEygOrgc1AY+emVkCBy3YtJCIbROQnEal+jXP1EJG1IrL2+HGb7zmzKpW3FCMbjmT5/uW8tOSlS+UBwQG0mNiC2EOxzH1srr2fYYybbqZJqpmI7AWCRORPHL/Yn76Zi4tIDmA60F9VTwPdcDRPrcNRm7lYczkCRKhqeee1vrnYv3E5VR2pqlGqGpUnT56bCclkEO3LtqdnVE8++uUjZm6beak8vHI4D77zIFunbWX9KJt0yRh3ZHF1R+eb3NuB24CyQASO5qlfVdXtR09ExA9HspigqjMAVPUP4CHn9mLAw87y88B55/o6EdkFFAPWuntdk3kMrjuYtYfX0mV2F0rnLU3RXEUBuO/Z+9izZA8L+y2kQNUC5C2V18uRGpM+uFzDcA5hvh0IVdXVqjpVVVfdZLIQYDSwTVU/uaw8r/OrD/AKMNz5Oc/FoUdEpDBQFNjt7nVN5pI1S1amtppKFp8stJzakrgLcQCIj9D0y6ZkzZmVaW2mcSHOXiEyxhXuNklNAOY5p2mtJSIPXlzcPE9VoCPw4GWPyjYAHhGR7TievDoMjHXuXwPYJCK/AdOAJ1T1pJvXNJlQwZCCfN3sa34/9ju95/e+1G+RI18Omn3VjONbjrN4wGIvR2lM+iDudPyJyJ5rbFJVLeyZkDwjKipK1661FivjMHDpQN5a9hajGo3i0XsevVS+6OlFrBq8ivYL21OkbhEvRmhM2iAi61Q16qrbMuqTIpYwzOWSkpOoP6E+y/Yt4+fuP3NP/nsASIxPZGTUSM6dPEfP33uSPVd2L0dqjHddL2HYfBgmU/D18WVC8wnkCcxDyyktOXXuFABZArLQ/OvmxP0dx7zH59mjtsZch7tDg7x5rSWlAjTGU/IE5mFqq6kcPH2QzrM6k+x8XiNfuXw88NYDbJu+jU1f2Sj9xlyLuzWMAlcsFYFngTs9HJcxKaJKeBU+fuhj5m6fy/sr3r9Uft+z9xFRPYL5feYTvTf6OmcwJvNyK2E45764fKkPNMcxNpQx6UKfSn1oW7otryx9haV7HIMU+Pj60Gx8MwBmdpppo9oacxWe6MP4DmjqgfMYkypEhFGNRlEsVzHaTm/LodOHAAiJDKH+Z/XZv3w/v3z8i5ejNCbtcbcPo/AVS2kcQ5vb8J8mXcnhn4PpradzNuHsvwYpvLvT3ZRoXoIfXvmBY5uOeTlKY9IWd2sYO3HMhbHTuawCqgOdPRyXMSmuZJ6SjGo0ipUHVvLC944xNS+OapstNBuzuswi6UKSl6M0Ju1wtw/DR1V9nV99VDWHqlZX1XUpFaAxKemRMo/Qu2JvPln1CdO3Tgcge+7sPDzsYY5uOMqK91Z4OUJj0g53m6SevUb5TY1Wa0xa8PFDH1Ppjkp0nd2V7Se2A1CieQlKty3NsreWWdOUMU7uNkkNvEb5K7caiDHecnGQQn9ff1pO+WeQwvqf1bemKWMu41LCuGyAQV8ReeDyQQdF5FEgNmXDNPM0UJIAABtCSURBVCZlRQRHMKH5BDb/tZme3/ZEVa1pypgruDofxmjn1wBgzGXlChwD+noyKGO8oW6Rugy8fyBv/PQGVQtUpUeFHv9qmrqryV2ElQ3zdpjGeI1LNQxVLaSqhXBMdlTosqWwqt6rqnNSOE5jUsWrNV6l7p116bugL+sOO57lsKYpYxzcfUqqk4iEiUgjEekqIt0uLikVoDGpydfHl6+bf01YYBgtp7bk5LmT1jRljJO7T0k1BXYBbwIjcDRFjcAxGZIxGULu7LmZ1noah04fouPMjiRr8r+apv7a8pe3QzTGK9x9SuptoKuqlgfOOr/2AOw9DJOhVLqjEoPrDmb+jvkMWj4IgHpD65E1Z1bm9ZiHJtsw6CbzcTdhRKjq1CvKvgQ6uXMSESkgIktFZJuIbBGRfs7yu0XkFxH5XUTmikjOy455UUR2isifIlLXzbiNcVuvir14pPQjDPxxIEt2LyEwTyB1P6nLgZ8PsHa4Tc5lMh93E8ZfInLxMZG9InIvjqHNfd08TyLwjKqWAKoAvUWkJPAF8IKqlgFmAs8BOLe1BUoB9YD/ExF3r2mMW0SEkY1GUjxXcR6Z/giHTh+ibMeyFK5TmO9f+J7TB097O0RjUpW7CWMUUM25PhhYCvwG/J87J1HVI6q63rkeC2wD7gCKA8ucuy0GWjjXmwCTVPW8qu7BMY5VJTdjN8ZtFwcpjLsQR+tprUlMTqTh8IYkJybzba9vbYY+k6m4mzA+VNXpAKo6HigGVFDVV282ABGJBMoDq4HNQGPnplY4JmkCRzK5fETcg86yK8/VQ0TWisja48eP32xIxvxLiTwlGN14ND8f+JkBiwcQWjiUB958gO1zt7N12lZvh2dMqnE5YTibgM6KSNaLZaq6X1W33ezFRSQHMB3or6qngW44mqfWAUFAwsVdr3L4f/60U9WRqhqlqlF58uS52bCM+Y82pdvQt1JfPl39KVO3TKVK/yrkvyc/C/ou4Nypc94Oz5hU4XLCUNUkYDuQyxMXFhE/HMligqrOcF7jD1V9SFUrABNxPMILjhpFgcsODwcOeyIOY1z10UMfUSW8Ct3mdGNH9A4afdGIuL/jWPzcYm+HZkyqcLdJagIwT0Q6i0ity8eUcuckIiI4hhvZpqqfXFae1/nVB8eAhsOdm+YAbUUkq4gUAooCv7oZuzG3xN/XnyktpxCQJYAWU1qQs1RO7n36XjaM3sDeH/d6OzxjUpy402knInuusUlVtbAb56kGLAd+By5OnvwSjkTQ2/l5BvCiOgMUkZdxNFkl4mjCWnC9a0RFRenatfboo/G8xbsWU/frurQv257RD41meNnhiK/Qc1NPsgS4OjybMWmTiKxT1airbsuoT3lYwjAp6a2f3mLgjwMZ9vAwHop+iK/qfMX9r99Pzddqejs0Y27J9RKGu01SiEgdERktInOdnyu42yRlTHr3co2XqV+kPv0W9uNEiROUbluaFYNWcHLnSW+HZkyKcXcsqb7AMBzzetdwFsfjGDLEmEzDR3z4qtlX5MuRj5ZTWxL1dhS+/r7M7zPf3s0wGZa7NYz+QG1VfY9/+h7+wPHCnTGZSq7suZjWahpHzxyl56qe1HyzJrsW7WLb9Jt+0tyYNM3dhBHEPy/QXfwzyo9/3pcwJlOpeEdFPq37KQt2LmBxmcXkK5ePhf0Xcj72vLdDM8bj3E0Yy4AXrih7EscQIcZkSk9EPUH7Mu15fcXr5Hk5D7GHY/nx9R+9HZYxHuduwugLNBORvUCQiPyJYwiPpz0dmDHphYgwouEISuYpyeN7Hqd45+KsHrKaY5uOeTs0YzzK3Rn3jgAVgTZAO6AzUFlVj6ZAbMakG4H+gUxvPZ34xHj+r9z/ERAawLc9v7V5M0yG4vZjteqwWlWnquoqVU2+8VHGZHzFcxdnTOMxLI9ezvGOxznw8wE2jtvo7bCM8Rh3H6v1F5E3RWSHiJx1fn1LRAJSKkBj0pNWpVrRr3I/Pg76mGz3ZGPxgMXEnYjzdljGeIS7NYxhwIM4OrorOr/ej5vzYRiTkX1Q5wOqRFRhWLVhxEfH88PLP3g7JGM8wt2E0RRoqKoLVHWrczynps7FGINzkMJWUzgTfoY/7v+DdSPXcXitDa5s0j93E8ZRIPsVZdmAI54Jx5iMITxnOJNaTGJm5ZkkBic6ZuezDnCTzrmbML4CForIYyJSX0R6APOB8Tc71LkxGVWtwrV4pd4rzH5wNofXHGb96PXeDsmYW+Kp4c0v59ZQ5ynFRqs1aUGyJtP4m8bkfTkvRU4Xof+O/mTPdWUl3Zi0w2Oj1apqIRcWrycLY9IKH/FhfPPxbGi7gfMx55n//Hxvh2TMTXP7PQxjjHtuy3Ybo54cxZoqa9g8ZjMHfj1w44OMSYPcfQ8jWEReFZEZIvLd5UtKBWhMRhB1exSN3mnEmcAzjOk8xjrATbrkbg1jKlAT+AGYfMVijLmOJ+5/grjucfAHvN3+bQ4ftkdtTfribqf3aSCXql64pYuKFADGA/lwzKsxUlWHiEg5YDgQgGPu7l6q+quI1ARmAxc73Weo6pvXu4Z1epu06GzCWV64/wVyr8pNgl8CF2pfoPnLzalStYq3QzMG8OwUrSuAErceEonAM6paAqgC9BaRksAHwBuqWg4Y6Px80XJVLedcrpssjEmrAv0D+eyXz3jg+wdIqJFA1sVZWVRtEf3L9GfCyAkkJSV5O0RjrimLm/t3AeaLyGrgX2M3u/NL3Dnq7RHneqyIbAPuwDEpU07nbsGA1dlNhlSjVg1q1KrB/r37+fKdL/Gf4s/Ox3cy4JUBhHUMo/uL3cmVO5e3wzTmX9xtkhoFNAaWA+cu26Sq2ummAhCJxDExU2kcSWMRIDhqP/ep6j5nk9R04CCOJPKsqm65yrl6AD0AIiIiKuzbt+9mQjIm1cXHxzN+6Hi2jdhGyO4Qzmc9j9ZX2rzShnIVynk7PJOJXK9Jyt2EEQsUc9YQPBFYDuAn4B1VnSEiQ4GfVHW6iLQGeqhqbRHJCSSr6hkRaQAMUdWi1zu39WGY9Oq7ed+x6INFBK4MxCfZh5h7YqjWvxrN2jfDx8eehDcpy5MJ4zeglqr+7YGg/IB5wCJV/cRZFgOEqKqKiAAxqprzKsfuBaKuF4clDJPe7di+g6/f+pr4mfFkP5udmNtjiOgaQdfnuhIcHOzt8EwG5cmE8SzQHPiM//ZhuDyGszMZfAmcVNX+l5VvA3qq6o8iUgv4QFUriEg+4JgzkVQCpgEF9TrBW8IwGcWZM2cY98k4dn2xi5ADIZzLfo4sjbLQ/tX2lCjliWdQjPmHJxPGtcaScmv8KBGphqMf5Hccj9UCvAScBobg6IyPx/FY7ToR6QP0xPF01TngaVX9+XrXsIRhMprk5GTmTZvHjx//SNCaIBCIrRRLrWdrUb9ZfWuuMh7hsYSRnljCMBnZ5k2bmfjWRJLnJRMQH0B0RDRFehShS/8uBAYGejs8k455NGGISB2gLZBXVRuJSAUg2J0mqdRgCcNkBqdOnWLsB2M5OO4gwUeDicsRR7YW2ej0aicK32njgBr3eezFPRHpi2Oa1h1ADWdxPPD2LUVojLkpoaGhPD3oaT469BElxpYgoUgC8qUwrvg4nqn5DEsWLvF2iCYDcbcPYxeOp6T2isgpVQ0VEV/gL1VNU28ZWQ3DZFbrfl3H1Hem4rvQF/8Ef04VOUXpnqXp1KcT/v7+3g7PpHGe7PT+C8ivqkkiclJVbxORAGCPqub3ULweYQnDZHZ/HfuLsYPGcvzr4wSdCOJMyBmC2wTT5eUuhBcI93Z4Jo3y5FhSy4AXrih7Elh6M4EZY1JO3rC8PP/p8ww6MojC/1eYxNsTSRqRxPA7h/Nc3edY+dNKb4do0hl3axj5gblAbhzDeOzG8ShsI1U9miIR3iSrYRjzXyt/XMnsQbPx/8Efv0Q/TpU8RYU+FWj7aFv8/Py8HZ5JA265SUpEHlHVic51ASoCBYEDwK+qmny9473BEoYx13bwwEHGvTOOmMkx5IjOQWyuWPJ0yEPXl7qSN29eb4dnvMgTCeP01YboSMssYRhzY/Hx8Xz9v6/ZMnwLITtDSPBPIKluEq1eaUWFShW8HZ7xAk8kjFhVDfJ4ZCnIEoYx7lmyYAkL3l9A9hXZ8U3yJbpcNPf2u5eWnVraW+SZiCcSRhzwMI5hx6/KXtwzJmPYvWs3498az7np58h+Jjsx+WII7xJO1wFdCQ0N9XZ4JoV5ImEkAfu4dsJwayyp1GAJw5hbc/bsWcYNHsfOUTsJ2R9CfEA8Po18aPdKO0qVLeXt8EwKsT4MY8xNS05OZsHMBSz5aAlBvwaBQmzFWGo+U5OGLRtac1UGYwnDGOMR2zZvY8LbE0icm0i2uGxEF4imcPfCdH2mKzly5PB2eMYDrNPbGONRMTExjP1wLPvH7Cf4SDBxgXEENAugw6sdKFrsupNhmjTOhjc3xqSI5ORkZk6YyYpPVxC8Pphkn2TO3neWhwY8RN1Gdb0dnrkJljCMMSlu49qNTH5nMrJAyHo+K6cKn6LE4yXo/GRnAgICvB2ecZElDGNMqjnx9wnGDBrD0a+OkvN4Ts7mPEtQ6yA6v9yZiMgIb4dnbsCTgw96hIgUEJGlIrJNRLaISD9neTkRWSUiG0VkrXP+bsRhqIjsFJFNInKPN+I2xtxYrty5eO7j5/jgyAcUGVGEhIgEkr9IZlTRUTxX+zmW/bDM2yGam+SVGoZzEMP8qrpeRIKAdUBT4FNgsKouEJEGwABVrelc7ws0ACoDQ1S18vWuYTUMY9KOVStWMfOdmfgt8cPvgh+n7jpFuV7laP9Eexv0MI1JczUMVT2iquud67HANhyj3ypw8fHdYOCwc70JMF4dVgEhzqRjjEkHqlSrwvsL3qf3nt749fLD75gf+57cx4v5XmRQ30EcO3rM2yEaF3i9D0NEInHMs1EaR9JYhOONch/gPlXdJyLzgPdUdYXzmCXA86q69opz9QB6AERERFTYt29fat2GMcYNCQkJfP1/X7Np2CZCt4eS4J9AYp1EWrzSgkpVKnk7vEwtzdUwLhKRHMB0oL+qngZ6Ak+pagHgKWD0xV2vcvh/Mp2qjlTVKFWNypMnT0qFbYy5Rf7+/nTr341P//yU+xfdz/lq5/Ff6M+CexfQ/+7+TPxiIklJSd4O01zBazUMEfED5gGLVPUTZ1kMEKKq6px3I0ZVc4rICODHy+bk+BOoqapHrnV+68MwJn3Zu2cv498ez5mpZwiMDSQmbwy3d7qdbi9247bbbvN2eJlGmqthOJPBaGDbxWThdBi437n+ILDDuT4H6OR8WqoKjkRyzWRhjEl/IgtFMnD0QN449gZhg8LQ7MrZj87y4R0f8lKLl/ht42/eDjHT89ZTUtWA5cDvwMXZ+l7CMd3rECALEA/0UtV1zgTzOVAPiAO6Xtl/cSWrYRiT/i2cs5DFHywm8JdAfJJ9iImKocZTNWjStokNephC7MU9Y0y6tv2P7Xz91tckzEpwDHp4RzSR3SLp+mxXcuZMV+OipnmWMIwxGcLp06cZ9/E49ozeQ8ihEM5lP4dfUz86vNKB4iWKezu8DMEShjEmQ0lOTmbOlDn89MlPBK9xDHp4psoZ6jxXh7qN61pz1S2whGGMybA2bdzE5Lcmkzw/mYD4AKIjoynWoxhd+nchW7Zs3g4v3bGEYYzJ8E6ePMnY98Zy+MvD5PzLMehhYMtAOr/SmchCkd4OL92whGGMyTSSkpKY+uVUVg9ZTcimEBJ9EzlX4xwNX2hIzYdqeju8NM8ShjEmU/r1l1+Z/s50sizOgn+CP6eKnaJsz7J06NUBf39/b4eXJlnCMMZkaseOHmPsu2P5e8LfBJ0M4kzoGUIeCaHrS125/Y7bvR1emmIJwxhjgAsXLvDNiG/Y8L8NhP4RSoJfAhdqXaD5y82pUq2Kt8NLEyxhGGPMFZYvXc6cd+cQ8GMAWRKzcKr0KSr2rUjb7m3x9fX1dnheYwnDGGOu4cD+A4x7exyxk2MJPB3I6TynCesQRreXupE7d25vh5fqLGEYY8wNxMfHM/6z8Wwbvo2Q3SGcz3oerae0ebUN5SqU83Z4qcYShjHGuGHxvMUs/HAhgSsD8UnyIaZ8DFWfqkrz9s0z/FvkljCMMeYm7Nqxi/FvjSd+RjzZz2YnJn8MEV0j6DqgK8HBwd4OL0VYwjDGmFtw5swZxn0yjl1f7CLkgGPQwyyNstDu1XaULFXS2+F5lCUMY4zxgOTkZL6d/i1LP15K0K9BIBBbKZYHn32QBs0aZIjmKksYxhjjYVs2beGbt78hea5z0MOIaIo8VoQuT3UhMDDQ2+HdNEsYxhiTQk6dOsW4D8dxYOwBgo8GE5cjjmzNs9FpYCcK31nY2+G5Lc0lDBEpAIwH8uGYonWkqg4RkcnAxVlQQoBoVS0nIpHANuBP57ZVqvrE9a5hCcMYk5qSk5OZ/tV0fv70Z0I2hpDkm8TZamep/3x9atev7e3wXJYWE0Z+IL+qrheRIGAd0FRVt162z8dAjKq+6UwY81S1tKvXsIRhjPGWdb+uY9o70/BZ6OMY9LDIKUo/UZoOvTsQEBDg7fCuK80ljP8EITIb+FxVFzs/C7AfeFBVd1jCMMakR8ePH2fsu2M59vUxcv6dkzMhZwhuE0yXl7sQXiDc2+FdVZpOGM5ksAworaqnnWU1gE8uBu3cZwuwHTgNvKKqy69yrh5AD4CIiIgK+/btS/kbMMaYG7hw4QKTvpjEus/XEbo1lAtZLnD+gfM0eakJ1WpW83Z4/5JmE4aI5AB+At5R1RmXlQ8Ddqrqx87PWYEcqnpCRCoAs4BSFxPM1VgNwxiTFq38aSWz352N/w/++CX6carkKSr0qUDbR9vi5+fn7fDSZsIQET9gHrBIVT+5rDwLcAiooKoHr3Hsj8CzqnrNjGAJwxiTlh08cJAv3/2S6EnR5IjOQWyuWHK3z023l7qRNyyv1+JKcwnD2UfxJXBSVftfsa0e8KKq3n9ZWR7nvkkiUhhYDpRR1ZPXuoYlDGNMepCQkMBXn3/F78N+J3RnKAn+CSQ9lETLV1oSVfmqv7dTVFpMGNVw/NL/HcdjtQAvqep8ERmH47HZ4Zft3wJ4E0gEkoDXVHXu9a5hCcMYk94sWbiEBe8tIPuK7Pgm+RJdLpp7+91Ly04tU+0t8jSXMFKDJQxjTHq1e9duvnr7K+KmxZH9THZi8sVwR6c76PZCN0JDQ1P02pYwjDEmHTp79ixffvolO0buIGR/CPEB8fg09KHdq+0oVbZUilzTEoYxxqRjycnJLJy9kO8//J6g1UGgEFsxlprP1KRhy4Yeba6yhGGMMRnEti3b+Obtb7gw5wLZ4rIRXSCawt0K0/XZruTIkeOWz28JwxhjMpiYmBjGfjiW/WP3E3w4mLjAOAKaBdDh1Q4ULVb0ps9rCcMYYzKo5ORkZk2cxfLBywleF0yyTzKJTRN5d/q7N3W+6yWM9D/bhzHGZGI+Pj40b9+cwWsH03RtUy40vkBoZMo8SZUlRc5qjDEm1ZWrUI5yM8ul2PmthmGMMcYlljCMMca4xBKGMcYYl1jCMMYY4xJLGMYYY1xiCcMYY4xLLGEYY4xxiSUMY4wxLsmwQ4OIyHFg3y2cIjfwt4fCSQ8y2/2C3XNmYffsnoKqmudqGzJswrhVIrL2WuOpZESZ7X7B7jmzsHv2HGuSMsYY4xJLGMYYY1xiCePaRno7gFSW2e4X7J4zC7tnD7E+DGOMMS6xGoYxxhiXWMIwxhjjkkydMESknoj8KSI7ReSFq2zPKiKTndtXi0hk6kfpWS7c89MislVENonIEhEp6I04PelG93zZfi1FREUk3T+C6co9i0hr57/1FhH5JrVj9DQXfrYjRGSpiGxw/nw38EacniIiY0TkLxHZfI3tIiJDnd+PTSJyzy1fVFUz5QL4AruAwoA/8BtQ8op9egHDnettgcnejjsV7vkBILtzvWdmuGfnfkHAMmAVEOXtuFPh37kosAEIdX7O6+24U+GeRwI9neslgb3ejvsW77kGcA+w+RrbGwALAAGqAKtv9ZqZuYZRCdipqrtVNQGYBDS5Yp8mwJfO9WlALRGRVIzR0254z6q6VFXjnB9XAeGpHKOnufLvDPAW8AEQn5rBpRBX7vkx4H+qegpAVf9K5Rg9zZV7ViCncz0YOJyK8Xmcqi4DTl5nlybAeHVYBYSISP5buWZmThh3AAcu+3zQWXbVfVQ1EYgBcqVKdCnDlXu+XHccf6GkZze8ZxEpDxRQ1XmpGVgKcuXfuRhQTERWisgqEamXatGlDFfu+XWgg4gcBOYDfVMnNK9x9//7DWW5pXDSt6vVFK58xtiVfdITl+9HRDoAUcD9KRpRyrvuPYuIDzAY6JJaAaUCV/6ds+BolqqJoxa5XERKq2p0CseWUly550eAcar6sYjcC3zlvOfklA/PKzz++ysz1zAOAgUu+xzOf6uol/YRkSw4qrHXqwKmda7cMyJSG3gZaKyq51MptpRyo3sOAkoDP4rIXhxtvXPSece3qz/bs1X1gqruAf7EkUDSK1fuuTswBUBVfwECcAzSl1G59P/dHZk5YawBiopIIRHxx9GpPeeKfeYAnZ3rLYEf1NmblE7d8J6dzTMjcCSL9N6uDTe4Z1WNUdXcqhqpqpE4+m0aq+pa74TrEa78bM/C8YADIpIbRxPV7lSN0rNcuef9QC0AESmBI2EcT9UoU9ccoJPzaakqQIyqHrmVE2baJilVTRSRPsAiHE9YjFHVLSLyJrBWVecAo3FUW3fiqFm09V7Et87Fe/4QyAFMdfbv71fVxl4L+ha5eM8Ziov3vAh4SES2AknAc6p6wntR3xoX7/kZYJSIPIWjaaZLev4DUEQm4mhSzO3sl3kN8ANQ1eE4+mkaADuBOKDrLV8zHX+/jDHGpKLM3CRljDHGDZYwjDHGuMQShjHGGJdYwjDGGOMSSxjGGGNcYgnDGGOMSyxhGGOMcYklDGOMMS6xhGFMKhGRD0Rk5mWfP3ROUuXnzbiMcZW96W1MKhGRXDgm+amJY5DD3kA1VY3xZlzGuMoShjGpSEReB5rjGPm4mqoeuP4RxqQd1iRlTOraAJQBXrRkYdIbq2EYk0pEpAyO0VS/A8JVtbaXQzLGLVbDMCYViMgdwFzgCaAXUEZEano1KGPcZAnDmBQmIjlxzE3wiarOUdU4HPOOvOPdyIxxjzVJGWOMcYnVMIwxxrjEEoYxxhiXWMIwxhjjEksYxhhjXGIJwxhjjEssYRhjjHGJJQxjjDEusYRhjDHGJf8PMb4FMQSdECoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(6,4))\n",
    "plt.plot(mesh, fem_solution['res'], label='num', color='green')\n",
    "plt.plot(np.linspace(x0, x1, 10**3), exact_solution, label='exct', color = 'purple')\n",
    "plt.xlabel(r\"$x$\", fontsize=12)\n",
    "plt.ylabel(r\"Temperature, $u(x)$\", fontsize=12)\n",
    "plt.legend(loc=\"best\", fontsize=12)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*References*\n",
    "\n",
    "1. D. Gilbarg and N. Trudinger. Elliptic partial differential equations of second order, 1998\n",
    "2. Finite Element Method Course Notes by T. LaForce, Stanford University, 2006 [[link](https://web.stanford.edu/class/energy281/FiniteElementMethod.pdf)]\n",
    "3. A. Larcher and N. Degirmenci. Lecture Notes: The Finite Element Method, KTH Institute of Technology, 2013 [[link](https://www.kth.se/social/upload/5261b9c6f276543474835292/main.pdf)]\n",
    "4. P. Arbenz. The Poisson Equation, ETH University [[link](http://people.inf.ethz.ch/arbenz/FEM17/pdfs/0-19-852868-X.pdf)]\n",
    "5. M. Asadzadeh. An Introduction to the Finite Element Method (FEM) for Differential Equations, 2016 [[link](http://www.math.chalmers.se/~mohammad/teaching/PDEbok/draft_FEM_version6.pdf)]"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}