FiniteElementApproach_1D.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "#-------------------------------------------------------\n",
    "def Lagrange_polynomial(x, i, points):\n",
    "    ''' Compute lagrangian interpolation polynomial \n",
    "    over domain x at specific points '''\n",
    "    p = 1\n",
    "    for k in range(len(points)):\n",
    "        if k != i:\n",
    "            p *= (x - points[k]) / (points[i] - points[k])\n",
    "    return p\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def Lagrange_polynomial_basis(x, points, N):\n",
    "    ''' Compute complete basis of lagrange polynomials \n",
    "    at specified points '''\n",
    "    return [Lagrange_polynomial(x, i, points) for i in range(N+1)]\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def Chebyshev_nodes(a, b, N):\n",
    "    ''' Computes Chebyshev nodes in a specific range'''\n",
    "    from math import cos, pi\n",
    "    return [0.5*(a+b) + 0.5*(b-a)*cos((2.*i+1.)/(2.*(N+1.))*pi) for i in range(N+1)]\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def Uniform_nodes(a, b, N):\n",
    "    '''Computes uniformly spaced nodes'''\n",
    "    h = (b-a) / N\n",
    "    return [a + i*h for i in range(N+1)]\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def init_elements(Omega, N, n, **kwargs):\n",
    "    ''' Function initializes a 2D array which defines N+1 global \n",
    "        1D finite elements, with n+1 nodes inside each element\n",
    "        in the domain Omega = (x0, x1) '''\n",
    "    node_spacing = kwargs.get('node_spacing', 'uniform')\n",
    "    if (node_spacing == 'uniform'):\n",
    "        nodes = Uniform_nodes(Omega[0], Omega[1], n + N*n)\n",
    "    elif (node_spacing == 'chebyshev'):\n",
    "        nodes = Chebyshev_nodes(Omega[0], Omega[1], n + N*n)\n",
    "    # Initialize elements\n",
    "    elements = [[i-j%(N+1) for i in range(j*(n+1), (j+1)*(n+1))] for j in range(N+1)]\n",
    "    return (nodes, elements)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finite Element Method in 1D using Least Squares Method\n",
    "\n",
    "##  Element matrix \n",
    "\n",
    "### Indice notation\n",
    "- $e$: Element indice\n",
    "- $r$: local node indice\n",
    "- $i$: global node indice\n",
    "- $i = q(e,r)$: Mapping of local node number $r$ in element $e$ to the global node number $i$. Mathematical notation for $i = elements[e][r]$\n",
    "\n",
    "### Element matrix definition\n",
    "$(d+1)\\times(d+1)$-Matrix, that stores the local contributions \n",
    "$$\n",
    "A_{i,j}^{(e)} \n",
    "= \\int_{\\Omega^{(e)}} \\phi_i \\phi_j \\text{d}x \n",
    "= \\int_{\\Omega^{(e)}} \\phi_{q(e,r)} \\phi_{q(e,s)} \\text{d}x\n",
    "$$\n",
    "\n",
    "Let $I_d = \\{0,...,d\\}$ be the valid indices of $r$ and $s$. We then introduce the notation\n",
    "$$\n",
    "\\tilde{A}^{(e)} = \\left\\{\\tilde{A}^{(e)}_{r,s}\\right\\}, \\qquad r,s \\in I_d\n",
    "$$\n",
    "For example, $d=2$ leads to\n",
    "$$\n",
    "\\tilde{A}^{(e)} = \n",
    "\\begin{bmatrix}\n",
    "\\tilde{A}^{(e)}_{0,0} & \\tilde{A}^{(e)}_{0,1} & \\tilde{A}^{(e)}_{0,2} \\\\\n",
    "\\tilde{A}^{(e)}_{1,0} & \\tilde{A}^{(e)}_{1,1} & \\tilde{A}^{(e)}_{1,2} \\\\\n",
    "\\tilde{A}^{(e)}_{2,0} & \\tilde{A}^{(e)}_{2,1} & \\tilde{A}^{(e)}_{2,2} \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "### Affine mapping\n",
    "Instead of computing the integrals in $\\tilde{A}^{(e)}$ for every element domain $\\Omega^{(e)} = [x_L, x_R]$, a reference element is introduced over the domain $[-1,1]$. \n",
    "\n",
    "A linear affine mapping from $X \\in [-1,1]$ to $x \\in \\Omega^{(e)}$ leads to\n",
    "$$\n",
    "x \n",
    "= \\frac{1}{2} (x_L + x_R) + \\frac{1}{2} (x_R - x_L) X \n",
    "= x_m + \\frac{1}{2} h X \n",
    "$$\n",
    "\n",
    "The integration over the reference element is a matter of changing the integration variable from $x$ to $X$. \n",
    "Introducing the basis function \n",
    "$$\n",
    "\\tilde{\\phi}_r(X) =  \\phi_{q(e,r)}\\left(x(X)\\right)\n",
    "$$\n",
    "associated with the local node $r$ in the reference element.\n",
    "The element matrix entry for the reference element reads then\n",
    "$$\n",
    "\\tilde{A}^{(e)}_{r,s} \n",
    "= \\int_{\\Omega^{(e)}} \\phi_{q(e,r)} \\phi_{q(e,s)} \\text{d}x\n",
    "= \\int_{-1}^{1} \\tilde{\\phi}_r(X) \\tilde{\\phi}_s(X) \\frac{\\text{d}x}{\\text{d}X} \\text{d}X\n",
    "= \\int_{-1}^{1} \\tilde{\\phi}_r(X) \\tilde{\\phi}_s(X) \\text{det}(J)  \\text{d}X\n",
    "$$\n",
    "Hereby the factor $\\frac{\\text{d}x}{\\text{d}X}$ is the determinant of the Jacobian matrix of the mapping between the coordinate systemxs. For 1D problems, this factor reduces to $\\text{det}(J) = h/2$.\n",
    "The corresponding entries for the element right hand side of the linear equation system reads\n",
    "$$\n",
    "\\tilde{b}_{r}^{(e)} \n",
    "= \\int_{\\Omega^{(e)}} f(x) \\phi_{q(e,r)}(x) \\text{d}x\n",
    "= \\int_{-1}^1 f\\left(x(X)\\right) \\tilde{\\phi}_r(X) \\text{det}(J)  \\text{d}X\n",
    "$$\n",
    "\n",
    "The $\\tilde{\\phi}_r(X)$ fuctions are the Lagrangian polynomials defined through the local nodes of the reference element.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Analytical soltion to element matrix and element rhs (only for d=1): \n",
      "h: 1\n",
      "A[i][i]: 0.3333333333333333\n",
      "A[i][j]: 0.16666666666666666\n",
      "b[0]: 0.08333333333333334\n",
      "b[1]: 0.08333333333333331\n",
      "\n",
      "Numerical soltion to element matrix and element rhs: \n",
      "[[0.333334 0.166666]\n",
      " [0.166666 0.333334]]\n",
      "[0.083333 0.083333]\n"
     ]
    }
   ],
   "source": [
    "# https://github.com/hplgit/INF5620/blob/master/src/fem/fe_approx1D.py\n",
    "import numpy as np\n",
    "from scipy import integrate\n",
    "\n",
    "#=======================================================\n",
    "# Implementation of the Finite Element Method in 1D\n",
    "# for the approximateion of a function f(x)\n",
    "#=======================================================\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def Lagrange_polynomial(x, i, points):\n",
    "    ''' Compute lagrangian interpolation polynomial \n",
    "    over domain x at specific points '''\n",
    "    p = 1\n",
    "    for k in range(len(points)):\n",
    "        if k != i:\n",
    "            p *= (x - points[k]) / (points[i] - points[k])\n",
    "    return p\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def phi_r(r, X, d):\n",
    "    ''' Defines a Lagrangian basis function for the  \n",
    "    local node r of an reference element of \n",
    "    degree d '''\n",
    "    nodes = np.linspace(-1., 1., d+1)\n",
    "    return Lagrange_polynomial(X, r, nodes)\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def basis(X, d=1):\n",
    "    ''' Create all basis functions for a reference \n",
    "    element of degree d\n",
    "    '''\n",
    "    phi = [phi_r(r, X, d) for r in range(d+1)]\n",
    "    return phi\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def element_matrix(X, phi, Omega_e):\n",
    "    ''' Creates the element matrix for an element over\n",
    "    its local domain Omega_e=[x0_e, x1_e] and its \n",
    "    basis functions stored in phi \n",
    "    '''\n",
    "    n = len(phi)\n",
    "    A_e = np.zeros((n,n))\n",
    "    h = Omega_e[1] - Omega_e[0]\n",
    "    detJ = 0.5*h\n",
    "    for r in range(n):\n",
    "        for s in range(n):\n",
    "            A_e[r, s] = np.trapz(phi[r]*phi[s]*detJ, X)\n",
    "            A_e[s, r] = A_e[r, s]\n",
    "    return A_e\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def element_rhs(X, f, phi, Omega_e):\n",
    "    n = len(phi)\n",
    "    b_e = np.zeros(n)\n",
    "    h = Omega_e[1] - Omega_e[0]\n",
    "    detJ = 0.5*h\n",
    "    # Compute mapping from f(x) to f(x(X))\n",
    "    x = np.linspace(Omega_e[0], Omega_e[1], len(X))\n",
    "    f_map = f(x)\n",
    "    for r in range(n):\n",
    "        b_e[r] = np.trapz(f_map*phi[r]*detJ, X)\n",
    "    return b_e\n",
    "\n",
    "\n",
    "\n",
    "# Test the fem-functions\n",
    "#-------------------------------------------------------\n",
    "d = 1\n",
    "Omega_e = (0, 1)\n",
    "X = np.linspace(-1, 1, 501)\n",
    "f = lambda x: x*(1-x)\n",
    "\n",
    "# Analytical solutions \n",
    "#-------------------------------------------------------\n",
    "print(\"\\nAnalytical soltion to element matrix and element rhs (only for d=1): \")\n",
    "h = Omega_e[1] - Omega_e[0]\n",
    "xm = 0.5 * (Omega_e[1] + Omega_e[0])\n",
    "print(\"h: \" + str(h))\n",
    "print(\"A[i][i]: \" + str(h/3))\n",
    "print(\"A[i][j]: \" + str(h/6))\n",
    "print(\"b[0]: \" + str(-(1./24.) * h**3 + (xm/6.)*h**2 - (1./12.)*h**2 - 0.5*xm**2*h + 0.5*h*xm))\n",
    "print(\"b[1]: \" + str(-(1./24.) * h**3 - (xm/6.)*h**2 + (1./12.)*h**2 - 0.5*xm**2*h + 0.5*h*xm))\n",
    "#-------------------------------------------------------\n",
    "\n",
    "# Numerical solution\n",
    "#-------------------------------------------------------\n",
    "print(\"\\nNumerical soltion to element matrix and element rhs: \")\n",
    "phi = basis(X, d)\n",
    "A_e = element_matrix(X, phi, Omega_e)\n",
    "b_e = element_rhs(X, f, phi, Omega_e)\n",
    "\n",
    "print(A_e)\n",
    "print(b_e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Global matrix / Finite Element Assembly Matrix\n",
    "The global matrix entries $A_{i,j}$ are the sum of the local element matrix entries \n",
    "$$\n",
    "A_{i,j} = \\int_{\\Omega} \\phi_i \\phi_j \\text{d}x = \\sum_{e=0}^{N} A_{i,j}^{(e)}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAEKCAYAAAALjMzdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xl4TOf7x/H3k8kmEQSxL0HVGmtiK2LfiSpqKVVKtdp+W1RpUVVdlZZSVUqqtNbaSositiKJNXYaSiyVpJYQ2SbP749J/VJbgkzOzOR+XddczJwz53wy1dxznvMsSmuNEEIIIeyfk9EBhBBCCJE1pKgLIYQQDkKKuhBCCOEgpKgLIYQQDkKKuhBCCOEgpKgLIYQQDsJqRV0pNUcpdVkpdeg+23srpQ6mPf5QSlW3VhYhhBAiJ7DmlXow0OYB208DgVrrasAHwLdWzCKEEEI4PGdrHVhrvVUp5fuA7X+ke7oLKGGtLEIIIUROYLWi/pAGAL9mZseCBQtqX19f66YRQgghbMiePXtitNY+Ge1neFFXSjXFUtQbPmCfQcAggFKlShEeHp5N6YQQQgjjKaX+ysx+hvZ+V0pVA2YDQVrr2Pvtp7X+Vmvtr7X29/HJ8IuKEEIIkSMZVtSVUqWAn4E+WusTRuUQQgghHIXVmt+VUj8BTYCCSqko4D3ABUBr/Q0wFigAfK2UAkjRWvtbK48QQgjh6KzZ+71nBttfBF601vmFEEI4luTkZKKiokhISDA6itW4u7tTokQJXFxcHun9hneUE0IIITIjKioKLy8vfH19SWvhdShaa2JjY4mKiqJMmTKPdAyZJlYIIYRdSEhIoECBAg5Z0AGUUhQoUOCxWiKkqAshhLAbjlrQ//W4P580vwshxANERUWxatUqrl69Sp48eahZsyYNGjRw+OIi7JMUdSGEuI9hw4YxefLk/7zm4eHBhZkzyZsvH3ToYFAyIe5Nmt+FEOIOWmsAXnnlFSZMmMCxY8dISEjgUmgov/n5kbdPH3RQEJNffpm4uDiD0wrx/6SoCyFEOsHBwTz77LMkJCRQrlw53n33XSqUK4fb119TuEkTGh06BJ9/zu7SpXnrm2+oX7s2MTExRscW2WzYsGFUrlyZgQMHEhgYiNlsvu++SUlJNG7cmJSUFKvnkqIuhBBp1q9fz8CBA7ly5QpOTmm/Hvftg3r1YOhQaNYMjhyBYcOot2EDv3l58eepU7Ru0YJr164ZG15km8jISHbs2MGRI0eoUaMGXbp0wWQy3Xd/V1dXmjdvzqJFi6yeTe6pCyEEll/UXbt2pXLlyixbtgzX5GR491344gsoWBAWL4auXeHfDnLlytFyxQp+btmSoIMH6dixIxs3bnzkSUPEQ3rjDdi/P2uPWaMGfPnlA3c5fvw4LVq0ICUlhZo1awKwfPny29ubNm3KO++8Q8uWLRk9ejTXr19n6tSpdO7cmVGjRtG7d++szXwHKepCiBxPa83gwYMBWL16NXl27oTBg+HMGRg0CD75BLy9735js2a0nTKF7197jRd37WLfvn3UqVMne8OLbFWhQgWef/55fH196du3L6VKlSL9cuDvv/8+Y8eO5fLly+zbt49Vq1YBULVqVcLCwqyeT4q6ECLH++uvv9i3bx8fv/MOpUaNgh9/hIoVYetWaNTowW8eMoSeBw7QdPZsipw+DVLUs0cGV9TWFBERQVBQEDExMeTLl+8/2xo3bozWmsmTJxMSEnK7Wd5kMuHq6kpcXBxeXl5WyyZFXQiR4/n6+nL8jz/I+9RTcO0ajBsHI0eCm1vGb1YKpk+nyLFj6H79WPjnn7QYOBBZJtpxHT58mCpVqpCYmHjX7G8RERFcvHiRggUL3lW8ExMTcXd3t2o26SgnhMjRDh48SEpKCvl/+QVTdDRs2QLvvZe5gv4vV1dYtozT3t70ffddRrz2mvUCC0PFxcXh4uKCh4cH3t7emM3m24X94sWL9O7dm5UrV+Lp6cm6detuvy82NhYfHx+r97mQoi6EyLGuXLlCo0aNeON//4O5cy1N5/XqPdrBChWi7Jo1DHd2JnjRIrZs2JC1YYVNOHToEFWrVr39vFWrVmzfvp34+Hi6dOnCpEmTqFSpEmPGjGHcuHG399u8eTPt2rWzej4p6kKIHGvq1Klcv36dgQ0bQkQEvPDC4x2wZk3GzJmDL/By9+4kJSZmRUxhQ+rXr8+SJUtuP3/11Vf5/vvv8fDwYOfOnbRs2RKw3FvfuXPn7f1+/PFHBg0aZPV8UtSFEDnStWvX+PLLLwkKCqL6jh3g7g49ejz2cT369GFq9+4cvXqV7/r0yYKkwpbVrFmTpk2bZjj5TOfOnalQoYLV80hRF0LkSN9//z1Xr17l3eHDLb3dn34a7ujJ/Kg6/PgjfUqUoOiyZbBxY5YcU9iu/v37Zzj5TN++fbMlixR1IUSOtGHDBgICAgi4cAGuXHn8pvd0lMnEvCNH6FypEnTrBpcuZdmxhXgQKepCiBxp5cqVlolB5s6FkiUtU8BmJS8v4r//ns+uXCF22rSsPbYQ9yFFXQiR46SmpuLk5EQRsxnWr4fnn4cHNJ8+qtPu7rwNTJ0+HVJTs/z4QtxJiroQIke5evUqvr6+LF26FObNsxTbfv2scq4qVarQ2d+fr65eJe6XX6xyDiHSk6IuhMhRFi9ezLlz5yjj62tpem/cGMqVs9r5Rk2ezBVg5pgxVjuHEP+Soi6EyFEWL17Mk08+Sa2EBDh5Mks7yN1LnUaNaF6yJJMOHiTh3DmrnksIKepCiBwjOjqazZs3061bN1RwMHh6WpZTtbJ3xo2jMvD3zJlWP5fI2aSoCyFyjOXLl5Oamkq3Dh1g0SLo3h1y57b6eZv178/G+vUpvWwZaG318wlj3Lp1i8DAwAwnomncuDEpKSlWySBFXQiRY/j7+/Puu+9S7fhxuHHD6k3v/zFwIBeOHePkwoXZd06RrebMmUOXLl0ynIimefPmLFq0yCoZrFbUlVJzlFKXlVKH7rNdKaWmKqVOKaUOKqVqWSuLEEIA1KpViwkTJqDmzoUnnoCGDbPt3OZnniFAKd56++1sO6eja9KkyV2Pr7/+GoD4+Ph7bg8ODgYgJibmrm2ZVa9ePc6cOQPA+fPn8ff3B2DBggUEBQXd3q9p06ZsSFvYZ/To0bz++usAdO7cmQULFjzmT39v1lxPPRiYBsy7z/a2QPm0R11gRtqfQgiR5Xbv3g1AnYIFUVu2wIQJlrXQs4kpTx5eqFmTj/buJXL/fsrWqJFt5xZZR2vN2bNnKV26NGBZutfPz4+kpCQiIyPx9fW9ve/777/P2LFjuXz5Mvv27bNMdgRUrVqVsLAwq+SzWlHXWm9VSvk+YJcgYJ7WWgO7lFL5lFJFtdYXrZVJCJFzffrpp+zZs4czfftaink2zcWd3ssffsgnbdsy8623+FSWZn1sISEh993m4eHxwO0FCxZ84Pb7OXXqFGXKlEGlfSH8t6jHxMSQ7461Axo3bozWmsmTJxMSEnK7Wd5kMuHq6kpcXBxeXl4PneFBjLynXhxIP74jKu01IYTIUklJSfz++++0ad0aNW8etGxpmRo2mxVv04YOefMyNyRElmW1UxEREfj5+d1+Hh4eTrVq1ciVKxcJCQl37Xvx4kXc3NzuKt6JiYm4u7tneT4ji/q92r3u2S1UKTVIKRWulAqPjo62ciwhhKPZuXMncXFxtC1eHM6ezd4Ocnd46bnniE1JYff33xuWQTy6f/75h1y5cgFw9OhR1qxZg5+fH97e3pjN5tuF/eLFi/Tu3ZuVK1fi6enJunXrbh8jNjYWHx8fXFxcsjyfkUU9Ckj/VbkEcOFeO2qtv9Va+2ut/X18fLIlnBDCcfz66684OzvT7PBhy/KqnTsblqXV+PH85e5Oo/BwwzKIR9e6dWs2btxI9+7dWbJkCQUKFKBw4cIAtGrViu3btxMfH0+XLl2YNGkSlSpVYsyYMYwbN+72MTZv3ky7du2skk9pK46ZTLun/ovWuuo9trUHXgXaYekgN1VrXSejY/r7++tw+Z9BCPEQ6tevj5vJRMiePZar9LQe0oZ54QVYuhR94QIqi++pOrKjR49SqVIlo2Pc1759+5g8eTI//PDDA/fr0qULH3/8MRUqVLjn9nv9nEqpPVpr/4wyWHNI20/ATqCCUipKKTVAKTVYKTU4bZe1QCRwCpgFvGKtLEKInC0kJIR57dtDQoKhTe//MvfvT7sbNxjTo4fRUUQWqlmzJk2bNs1w8pnOnTvft6A/LqteqVuDXKkLIR5JvXqWCWciIrJ1KNs9aU27PHk4lJTEmVu3cHKSecAyw9av1LOKTV6pCyGELRg/fjyfDB0Ku3dbrtKNLugAStHn6ac5l5TEljlzjE4jHIgUdSGEw9JaM2vWLPatXQsmEzz3nNGRbgv68EO8gB8mTzY6il2xt9blh/W4P58UdSGEwzp9+jRRUVEEXrwI7dtDWi9lW+BRsiRdy5RhydGjxMfGGh3HLri7uxMbG+uwhV1rTWxs7GONX7fmNLFCCGGoLVu2ABB4/bpNdJC70+Dhw6kyZAipy5fDiy8aHcfmlShRgqioKBx5vhJ3d3dKlCjxyO+XjnJCCIf1wgsvsPqnn7js7o5TdDRYYbKPx6I1PPkkFCsGaV9AhLgX6SgnhMjxPDw8CHJ3x6lRI9sr6ABKEd+3L/O3buXvHTuMTiMcgBR1IYTDmv7ee3x37Ro0amR0lPs626QJfYAfR482OopwAFLUhRAOKSUlBbZvtzzJxnXTH1bFRo0I8PZm3vbtkJRkdBxh56SoCyEc0quvvkq9IUPA3R38M7wVaag+zz7L/pQUDn31ldFRhJ2Toi6EcEg7d+4kX3w81K0Lrq5Gx3mgHu+9hzPwwzffGB1F2Dkp6kIIh3Pjxg0OHTpE3bg4m76f/i+fIkVoWbo0EZGRcPOm0XGEHZOiLoRwOOHh4aSmplJXa5u+n57ekq+/Zm1qKvz2m9FRhB2Toi6EcDi7d+8GoI5SUL++wWkyx7NVKyhQAL18udFRhB2Toi6EcDi1atViRMmSFKxZE/LkMTpO5jg7M6lMGar/9BNaesGLRyRFXQjhcFoGBvJpTIxd3E9Pr+BTTxGRmkrYzJlGRxF2Soq6EMKhxMXFcXz5clJv3bK7oh709tu4AEu++87oKMJOSVEXQjiU9evXU7FHD8LBbjrJ/Stf0aK0KlKExYcOoc1mo+MIOyRFXQjhUHbt2oWrkxPVy5WzqaVWM6t7UBBnzWZCg4ONjiLskBR1IYRDCQsLo4ZSuAUGGh3lkXQaOZLhSlEoLMzoKMIOSVEXQjgMrTX79+6lltlsd/fT/5XP15eJLVpQZvNmy9KsQjwEKepCCIdx+vRprsXFURPs7n56eimdOvH7iROclIloxEOSoi6EcBiFChVieaNGtPXxgXLljI7zyG62aEF74JsPPzQ6irAzUtSFEA4jd+7cdD5zhpJNmoBSRsd5ZHkrVqRVvnwsCwtDSxO8eAhS1IUQDmPFrFmEnztnt/fT03u6eXP+SkpivzTBi4cgRV0I4TBeGTGCqeAQRb3j8OE4ASumTjU6irAjVi3qSqk2SqnjSqlTSqmR99heSim1WSm1Tyl1UCnVzpp5hBCO6++//+bi1avUdHMDPz+j4zw2n3r1eMrDgy27dhkdRdgRqxV1pZQJmA60BSoDPZVSle/YbTSwWGtdE+gBfG2tPEIIx7Zv3z4AalarBiaTwWmyxqIXX2Tj9esQG2t0FGEnrHmlXgc4pbWO1FonAQuBoDv20cC/SyjlBS5YMY8QwoHt27EDgBqtWxucJOsU7dMHU2oqrF5tdBRhJ6xZ1IsD59I9j0p7Lb1xwHNKqShgLfCaFfMIIRzY/s2bKQPka9XK6ChZp3ZtvsiXj16jRxudRNgJaxb1e40nuXNsRk8gWGtdAmgH/KCUuiuTUmqQUipcKRUeHR1thahCCHs319+fdS4uEBBgdJSsoxTXn3yShefPc/n0aaPTCDtgzaIeBZRM97wEdzevDwAWA2itdwLuQME7D6S1/lZr7a+19vfx8bFSXCGEPfPYtYvydeqAu7vRUbJU54ED0cCqiRONjiLsgDWLehhQXilVRinliqUj3Ko79jkLNAdQSlXCUtTlUlwI8VAuRUYyMjSUE1WqGB0ly1V7/nl8nZxYserOX59C3M1qRV1rnQK8CqwDjmLp5X5YKTVeKdUpbbdhwECl1AHgJ6CflumThBAP6eDSpXyqNZfseGrY+1EuLjxdsSIbzp8n7p9/jI4jbJyzNQ+utV6LpQNc+tfGpvv7EeApa2YQQji+E1u2APBk0J0DbBxDt759iRk5kuvr1uHVs6fRcYQNkxnlhBB270REBF5OThR+8kmjo1hF/ddfZ56HB8W3bzc6irBxUtSFEPYtJYUTFy5QvkABlB0v4vJAuXKhW7fm8JIlJCUkGJ1G2DAp6kII+7ZvH3+bzTxZvrzRSaxq/RNPUDU6mpBvvzU6irBhUtSFEPZt+3b2AsHz5xudxKoav/EGnsDyuXONjiJsmBR1IYR927YNVbYsbmXKGJ3EqnIVK0abQoVYGRFBqtlsdBxho6SoCyHsl9Zs3bSJ3lpz8eJFo9NYXee2bbloNhO6ZInRUYSNkqIuhLBfx4+z+9o1fjx9mly5chmdxurav/UWzsCKGTOMjiJslFXHqQshhFVt28ZJwCd/fvLly2d0GqvzrlKF9RUrUjsuzugowkbJlboQwn5t384JFxeerFTJ6CTZpmm/fuTZtw/Onct4Z5HjSFEXQtivbds47uREeQcfzpaeDgriU+DHMWOMjiJskBR1IYR9On+ehNOn8cmfn+rVqxudJtuoihVZ7O7OtBUrjI4ibJAUdSGEfdq2DXfg4KpVvPHGG0anyVZBAQHsunaNv48fNzqKsDFS1IUQ9mn7dsidG2rUMDpJtuv80ktoYPVnnxkdRdgYKepCCPu0bRufFClCq3btyGkrNvv16IGvycTKtWsz3lnkKBkOaVNKuQMdgEZAMeAWcAhYo7U+bN14QghxD3FxEBHBzief5Pz58467kMt9KJOJbn5+HIuIQMfHozw8jI4kbMQDr9SVUuOAHUB9YDcwE1gMpACfKKU2KKWqWTukEEL8x969oDXHbtygUg4azpbep59+yiqzGbVxo9FRhA3J6Eo9TGs97j7bJiulCgGlsjaSEEJkIDSUJODPS5foVrGi0WkMoZo0gbx5ubVkCbk6djQ6jrARD7xS11qvgdtN8P+hlCqotb6stQ63VjghhLinsDBOFSuG2WymYg4t6ri68kHp0pSeP58UWWNdpMlsR7kwpVS9f58opZ4B/rBOJCGEyEBYGOaqVWnbtm2OGqN+p0otWhCtNTtkjXWRJrNFvRfwlVJqolJqATAQaGa9WEIIcR/R0XDmDH4tW7J27Vr8/PyMTmSY1sOG4QqsDA42OoqwEZkq6lrrCOBDYDDQFHhVax1lzWBCCHFPYWEApNaubXAQ43kVK0YLHx9WRkSgU1ONjiNsQKaKulLqO+ANoBrwArBaKTXEmsGEEOKeQkPByYn6I0bw/PPPG53GcEFt2hCZksKhn382OoqwAZltfj8ENNVan9ZarwPqAbWsF0sIIe4jLIzUihU5fPRojlhuNSNBI0YwBSgaGmp0FGEDMtv8/oVON2WT1vqa1nqA9WIJIcQ9aA2hoZytXJmbN29SpUoVoxMZrnDVqrzesCEF1683OoqwARlNPrNaKdVRKeVyj21llVLjlVL9rRdPCCHS+esviInhsI8PgBT1NNfbtOH7Awe4JFfrOV5GV+oDsUwPe0wpFaaUWquU2qSUisQyu9werfUcq6cUQgiw3E8HDjtb5s2Som4RVbs2/YDlEycaHUUY7IEzymmtLwEjgBFKKV+gKJa5309oreMzOrhSqg0wBTABs7XWn9xjn+7AOEADB7TWvR7uRxBC5BhhYeDqSo3WrRnq4iL31NNUat2aJ1xdWblxIy8bHUYYSllrdSOllAk4AbQEooAwoKfW+ki6fcpjmUu+mdb6ilKqkNb68oOO6+/vr8PDZRI7IXKkJk0gIQF27TI6ic0ZXr8+U3ftIubPP8lTtqzRcUQWU0rt0Vr7Z7RfZoe0xSmlrqc9EpRSZqXU9QzeVgc4pbWO1FonAQuBoDv2GQhM11pfAciooAshcjCzGcLDSfX35/Tp06TKuOz/6DxwIMnAr9IEn6Nltve7l9Y6T9rDHXgGmJbB24oD59I9j0p7Lb0ngSeVUjuUUrvSmuvvopQapJQKV0qFR0dHZyayEMLRHDsGN28S6etL2bJlCZZZ1P6jft+++Dg5Eb5hg9FRhIEyO079P7TWK8h4mth7LXB8Z1u/M1AeaAL0BGYrpe66Saa1/lZr7a+19vdJ6/UqhMhh/u0klysXIJ3k7mRydubYiy8y8fx5uHnT6DjCIJltfu+S7tFVKfUJdxfoO0UBJdM9LwFcuMc+K7XWyVrr08BxLEVeCCH+KywM8uTh8NWrAFSuXNngQLYnf48elj4H69YZHUUYJLNX6h3TPVoDcdx9f/xOYUB5pVQZpZQr0ANYdcc+K7DMJY9SqiCW5vjITGYSQuQkoaHg78+hw4cpVaoUXl5eRieyPY0a0c/VlVFjxxqdRBjkgUPa/qW1fuFhD6y1TlFKvQqswzKkbY7W+rBSajwQrrVelbatlVLqCGAG3tJaxz7suYQQDi4hAQ4ehKFDOfzrr1StWtXoRLbJ2ZkrhQuz6cgRPkpKQrm6Gp1IZLNMFfVHpbVeC6y947Wx6f6ugaFpDyHs06pVJIeEcLJqVS6WKIGHlxc+Pj6UK1cOpe7VtUQ8tAMHIDkZ6tTh/Xr18PT0NDqRzQrq1IlV06ezb84cag0ebHQckc2sWtSFcGgpKfD22/w1eTLlgeQ7No955x3Gf/ihEckcT9pyqwQE0LlkyQfvm8N1fOstnKZPZ/l330lRz4GkqAvxkFJSUpg9aRJnp0/no3PnKP3aa3zk40PR2FhKHDrErT/+4OytW9SZOhXOniXUz48v9+zhi6lTKVy4sNHx7VNoKBQpwslbt/h7+3bq1auHs7P8+roXn9KlaVywID/v388HWoO0FuUoj/R/hVIqCLiktd6dxXmEsGnHjh2jV1AQ+06cINDJiZTvv8e5b1+Gp98pMRF+/x2WLYOVKzk2fz7LgZBVq1g4ciSNx4wBp0caTZpzhYVBQABz5s5l0qRJ3Lhxw+hENu2lrl058M03JO/ejUu9ekbHEdnoUX+z1AVGK6V+zcowQtiyed9/T+3q1Tl74gRLihRh8969OPfte/eObm7Qvj3MmQOXLtF3wwZ2d+uGV3IyzcaNY96zz2Z/eHt27RocPw516nDgwAEqV66Mq3QAe6AeEybwscmEyy+/GB1FZLNHnXzmHa11R61126wOJIQtunDqFC8PGECdpCQOtGhB16NHUdWrZ/xGFxdo0YJqixcTHhND00KFeGHpUrYsWGD90I5izx7LOuoBARw4cIBq1aoZncj2FShAUsOGbJd/ZzlOZiefKaSUelopNUQp1V8pVUcpJe2HwuFpreHkSYo9/TTbzWbWjx9P8XXr4BFWB/PKl4+Vf/zBB25uNJgzx1KoRMbSOsnFlC3LhQsXqJ6ZL1OC6d7eNDpzhtObNhkdRWSjBxZmpVRTpdQ6YA3QFsvSq5WB0UCEUup9pVQe68cUIvulpqYyuE0bZlSrBhcvUnPdOlwe8364R7lyvDN5Mi6bNhEzcyZnz57NwsQOKjQUypXj4DnLUhJypZ45nYcNA+Dnzz83OInITg9celUpNRH4Smt9128epZQz0AEwaa2XWS/if8nSqyI7aK15vW5dpoWFMbpoUT7YuRNKl86ag5vNpDZoQK09ezBVrcqOXbtwd3fPmmM7olKloGFD4mfPZv/+/fj5+clscplUK3du3FNT+ePmTekFb+eyZOlVrfVb9yroadtStNYrsrOgC5FdJvbty7SwMIZWrsz4P//MuoIOYDLhNGsWH2jN3gMHeOedd7Lu2I7m0iU4dw4CAvDw8KBBgwZS0B9Cl2bN2HnrFhd+/93oKCKbZPae+g9KqbzpnvsqpTZaL5YQxlk6ezZvz5/Ps/nzM3HfPlTaqmBZqlo1Og4fzmBgypQp7NmzJ+vP4Qj+nXSmTh2++uortmzZYmweO/PMqFEArJY11nOMzN4c3A7sVkq1U0oNBNYDX1ovlhAG0ZpzU6bwlFLM/f13nKw5dOq99/i4VCkKKcXAF18kJSXFeueyV6GhYDKRXLUqw4cPZ82aNUYnsiuV6tdnd8OGDDx0CMxmo+OIbJCpoq61ngm8CKwExgONtdarrRlMCEPMn8+bhw4R8skn5KpZ07rn8vAg38yZTDWbKRYfT1xcnHXPZ4/CwqBKFY6dPUtSUpL0fH8EdV57DaeLFyEkxOgoIhtktvm9DzAH6AsEA2uVUvJ/l3AYqamp9H/2WdYNHgwNG+Kc1nPY6tq0oeuzz7L69Gm8L1/OnnPaC60tRb1OHfbu3QtA7dq1DQ5lf261aMHrLi4smTDB6CgiG2S2+f0ZoKHW+iet9ShgMJbiLoRDmDplCnMXL+Z4cjIEB4PJlG3nVl9+ifL05ETfvlz5559sO6/Ni4yEf/6BgAD27NmDp6cn5cuXNzqV3XH39uZXd3e+3bYNbt0yOo6wssw2v3fWWl9O9zwUy1SxQti9AwcO8PaIEXQCXvviCyhXLnsDFCnCsTffpEJoKEuGyirEt6XrJHf8+HFq1qyJKRu/bDkKpRQ9OnVik9nMpR9+MDqOsLKMxqmPBr7WWt/z8kEp1Qzw0Fpn2wTDMk5dZKVbt27hX706/5w6xcGGDfEJCTFksRVtNlPO05NKqamsuXABChbM9gw2Z+hQmDEDrl9HOztz9epVvL29jU5ll45ERFClWjWmVq3KaxERRscRjyBLxqkDEcBqpdRGpdREpdQIpdTYtCFuEUBHQFZqE3Zr/rx5HDl5ku/d3fGZP9+w1dOUyURQ9+5sTE7mxuuvG5LB5oSFQa1a4OKCUkoK+mOo7OeHX8GC/HToEMTEGB1HWFFGk8+s1Fo/heUe+mHABFwH5gN1tNZvaq2jrR9TCOtc+9nfAAAgAElEQVR4MSGBP4BW06ZZZi4zUFD//iQC63/6CTbm8GkgUlIsC7kEBLB27VpeeOEFrly5YnQquzZowAAqASkLFxodRVhRRs3vP2it+yil/qe1npKNue5Lmt9FVkhKSuLSjh2UatcOmjeH1asNn0YzJSWFQoUK0SklheDChSEiAnLq9LEHD0L16rBgAUPDw5kxYwZxcXE4Ozsbncx+aQ1+fpbFiLZvNzqNeEhZ1fxeWylVGuivlPJWSuVP/8iaqEJkv08++ojKLVpw1s0NZs0yvKADODs788svv/DFvHlw6hR8/LHRkYwTGmr5M63ne/Xq1aWgPy6l0L17c3jHDsvIAuGQMirq3wC/ARWBPXc85HJZ2KVDhw4xYcIEOqWmUmrGDCha1OhItzVo0ADvzp2hSxeYNg0SE42OZIywMMiXj9SyZdm3b5+MT88ic52dqQoc/VImBHVUGd1Tn6q1rgTM0VqX1VqXSfcom00ZhcgyZrOZAb16kddsZkqHDtCjh9GR7jJt2jS+KVHCMkZ7xQqj4xgjLAwCAjj155/ExcVRq1YtoxM5hLbPPYcCFv74o6U5XjiczI5Tf9naQYTIDtOnTiU0IoKpXl74zJljE83ud1qzZg2fr1mDLlUKZs82Ok72u3XLck+9Th1iYmKoUKECderUMTqVQyhatChNKlTgp9hYtCwi5JCMGb8jhEHOLFtGG6BHcDD4+Bgd556CgoL4888/OdqxI/z+O5w5Y3Sk7LV/v2XxkYAAGjRowLFjx/Dz8zM6lcPoOXgwJ4G9kycbHUVYgRR1kXNcuMDk/ftZ3bEjqksXo9PcV8eOHQFY6elpaUmYO9fgRNksXSc5kfWe6dsXZ6VYvGqVZeigcChWLepKqTZKqeNKqVNKqZEP2K+rUkorpTLsri/Eo9i7dy/hAwdCcjLOX3xhdJwHKl68OAEBAawICYFWrSxFPSctmxkWBsWLk1igACVKlGBuTvtSY2X58+fn9/ffZ9zNm7Bpk9FxRBazWlFXSpmA6UBboDLQUylV+R77eQGvIzPTCSsxm80M6NGDZ9auJfl//8v+ud0fQffu3cmdOzeJffvCuXOwYYPRkbJPaCgEBLB//37Onz9P3rx5jU7kcALfeotc+fLB/PlGRxFZzJpX6nWAU1rrSK11ErAQCLrHfh8AnwEJVswicrBvZsxg/8mTfJ4vHy5jxhgdJ1OGDRvGxo0bcXvmGcs88N99Z3Sk7HHlCpw8CXXqsHu35Xt+3bqydlSWc3fnmypVGLVwIdy8aXQakYWsWdSLA+fSPY9Ke+02pVRNoGR2Lggjcpbo6GhGv/02zYGuU6aAl5fRkTJFpfXKv5aQAH36wMqVEJ0DZmT+d7bIgABCQ0MpVqwYxYsXf/B7xCM5UrgwXyQnc/Wnn4yOIrKQNYv6vcYK3R4YqZRyAr4AhmV4IKUGKaXClVLh0TnhF5vIMqNHjOBGfDxfVa+Oeu45o+M8lB9++IGCBQtyvkMHSE7OGU2l27dbFtVJu1KXq3TreX7kSBKBRVNsYgZwkUWsWdSjgJLpnpcALqR77gVUBUKUUmeAesCqe3WW01p/q7X211r7+9joMCRhm3wjI3kbqPTtt4atwPaoAgICSElJYfmRI1C3rmXMuqNPGLJlC9SsSWru3HTu3Jlu3boZnchh1fL3p2rBgsw5dAguXzY6jsgi1vwtFwaUV0qVUUq5Aj2AVf9u1Fpf01oX1Fr7aq19gV1AJ621TD8rssbJk4zauZMJ/fqBHU5eUrFiRapUqcLSpUvhxRfhyBHY7cD9SRMSYNcuCAzEycmJiRMn0rNnT6NTOSylFC8OGEAocEDGrDsMqxV1rXUK8CqwDjgKLNZaH1ZKjVdKdbLWeYUACAkJYVnPnmhXV7teGKVr165s3bqVv5s2BU9Px+4wFxpqmes+MJBz586RmFPnvc9GfUaMoL2XF+Y1a4yOIrKIVdsjtdZrtdZPaq3Laa0/THttrNZ61T32bSJX6SIrpKSkMOT55xm5Zw/J77wDRYoYHemRde3aFa01yzdsgO7dYeFCuHHD6FjWsXWrZbKdRo147rnnaNasmdGJHF7+/Pn5ZcwYah06ZBl1IOyefd1kFCITZs+cyZGzZ/msSBFch2XYD9OmValShW+++YZ27drBgAGWgr54sdGxrGPLFvDzI9HDg9DQUOkkl1169iQK2D9pktFJRBZQ2s463vj7++vwcLmgF/d2/fp1nihenIo3brBl1SpU2pSrDkFrqFwZ8ueHHTuMTpO1kpMhXz4YMICdPXvSoEEDli1bRhcbns7XkVTNnRuv5GR23rwJsm69TVJK7dFaZzjrqlypC4fy8ejRRN+4weR69VAdOhgdJ0uYzWZ+/PFHNm7aZLla/+MPOHrU6FhZKzwc4uMhMJAdaV9YnnrqKYND5RwDn32WXUlJhH/2mdFRxGOSoi4cSs3DhxmhFP5z59rksqqPwsnJibFjxzJx4kTo29dyJeVoHea2bLH82bgx27dvp3z58hQuXNjYTDlIv88/J7dSTJ00yfGHTTo4KerCcezdS/fNm/n0jTegYkWj02QZpRTdunXj999/J1op6NgR5s2DpCSjo2WdrVsttxZ8fBg1ahSTZYhVtsrr7U2/wEAW/vMPl37+2eg44jFIURcOISw0lIldupBYoACMHWt0nCzXq1cvzGYzixcvtjTBR0fDLw4yu3JKimUmucaNActc7x0c5NaJPXlt6lQAtk6YYHAS8TikqAu7p7Xmzb59mfTXXySNHWvpcOVg/Pz8qFq1KgsWLIDWraF4ccdpgt+/H+LiIDCQXbt28csvv5Cammp0qhznST8/LowcSff9++HQIaPjiEckRV3YvWXz57Pj+HE+KFkSr1deMTqO1fTu3Zvr169zMzER+vWD336DqCijYz2+f++nBwYyZcoUBg0adHtBG5G9Cg4fDh4e3PrkE6OjiEckRV3YtcTEREa8+ip+QP+FC8FkMjqS1QwbNoyIiAg8PT2hf39ITYXgYKNjPb4tW6B8eXSRImzevJmmTZtKUTdKgQK8XLYszX78Ec6fNzqNeARS1IVdmzZqFKevX+fz9u0xNWhgdByrcnFxQSlFQkICukwZaNoU5syxFHd7lZoK27ZBYCBHjx7l77//lpnkDFa1a1d2ac3ukSONjiIegRR1Yb/MZhr89htveXjQ6ocfjE6TLUJCQihUqBD79u2zLPJy+jSEhBgd69FFRMDVq9C4MZs2bQKQom6wvkOHksfFha8WLYLr142OIx6SFHVhv2bNov7Ro3w2cyZ4exudJltUq1aNhIQES4e5p5+2dAqcPdvoWI8u3f300NBQSpcuTZkyZYzNlMN5eXnRv2tXFiUnEzVxotFxxEOSoi7s0tFt2xjyv//xT8OG0Lu30XGyTf78+WnXrh0//fQTZldXy8/+889w5YrR0R7Nli3g6wulShEcHHx7NjlhrP999BEa+OrLLy1T+Aq7IUVd2KW3evZkflIS5s8+c5iZ4zKrd+/eXLx4kZCQEHjhBctypUuWGB3r4WltmXQmMBCwzJxXvHhxg0MJAF9fX5aOHs3IGzcsKwMKuyFFXdidDZ99xprz5xndvDk+9esbHSfbdejQAS8vL0sTfK1aUKkSzJ9vdKyHd/QoxMRA48YEBwczaNAgkuWq0GZ0Hj8e7ypVYOJEmTrWjkhRF3bFHB/P0DFjKOPszOtLlxodxxC5cuVi5syZvPrqq5ZWiuees/QgP3PG6GgPJ9399Pnz5/PHH3/g4uJibCbx/5RiU7t2tIuIIHHNGqPTiEySoi7sypyePTmUlMRnb7+NmwPOHJdZPXv2pFatWpYnvXpZ/lywwLhAj2LLFihenBuFCrFt2zbatGljdCJxB3OTJvwKBL/1ltFRRCZJURf249QpWvz2G2MrV+aZDz4wOo3hwsPD+fDDDy0dzRo3tjTB20szqdaWoh4YyOaQEJKSkmjbtq3RqcQdWrRtS92SJfnk2DGSQ0ONjiMyQYq6sA9aw5AhlHFz4/0NG2TGMWDjxo2MHj2aY8eOWZrgjx2DPXuMjpU5J0/CpUsQGMivv/6Kp6cnDRs2NDqVuINSijETJ3IGmP/660bHEZkgRV3YhTPTpxO0fj2Rb74JxYoZHccm9OvXD2dnZ2bPng3duoGrq/10mNu61fJn48YUKFCAXr164ebmZmwmcU/tunenVqFCfLR7Nyl//ml0HJEBpe2luS6Nv7+/Dg8PNzqGyE7XrtGzSBFWJiZyPDKSkr6+RieyGc888wxbt24lKioKt969LR3mzp8HZ2ejoz1Ynz6wYQNcvJjjhiTaow0LFhDVty99X30V05QpRsfJkZRSe7TW/hntJ1fqwuZtHzCAhQkJDO/fXwr6HQYOHEhMTAwrV660NMFfvmwplrbs3/vpjRsTExuLvV1Y5EQte/fmhV69MH33nf1OdJRDSFEXNi1l926GLFtGydy5eVuuEO7SsmVLatWqxbVr16BdO8if3/ab4M+cgXPnIDCQbt260bp1a6MTiUxIHTqUKTdvMu+ll4yOIh5AirqwXWYzs7t35yDwxfTpliVHxX+YTCbCw8MZOHCg5Z569+6wfDnExRkd7f7SxqfH+PmxdetW6tata3AgkRmqRg1+9vbm7WXLiJerdZslRV3YrgkT6HX2LDP696dLnz5Gp7FZSilSU1M5ffq0pQn+1i1LYbdVW7dC/vysPnmS1NRUnn76aaMTiUxQSvHhuHFcSk1lcrduRscR9yEd5YRt2rSJ1ObNcerTB77/XjpTZeDFF19k7dq1nP3rL5wrVIAnnoD1642OdW/lykG1anQymzlw4ABnzpyRIYp2pEvJkqyPiuLU1q0UadTI6Dg5hk10lFNKtVFKHVdKnVJKjbzH9qFKqSNKqYNKqY1KqdLWzCPsxKVL7OjaFT9XV04MHSoFPRM6duzIxYsXWbV6teVqfeNGuHDB6Fh3i4qCyEhu1KvH+vXr6dy5sxR0O/Pp4sUkAu917w6pqUbHEXewWlFXSpmA6UBboDLQUylV+Y7d9gH+WutqwFLgM2vlEXbCbCalRw+GXL1KXP78FC9f3uhEdqFDhw6ULl2aKVOmWIp6air89JPRse6Wdj/dtWlTlixZwkvS6crulK9fn0m9etHz0iWYNcvoOOIO1rxSrwOc0lpHaq2TgIVAUPodtNabtdbxaU93ASWsmEfYg/HjmbplCwe0ZvJXX0nnuEwymUy89tprbN26lf3x8VCnjm32gt+6FfLmxbV2bTp27Ejlynd+zxf24PX582nSvDmMGGGZF0HYDGsW9eLAuXTPo9Jeu58BwK/32qCUGqSUCldKhUdHR2dhRGFTfv+dk+PH867JRMeOHXnmmWeMTmRXBgwYgKenJ8HBwZar9f374dAho2P915YtxNWrx9j33+fcuXMZ7y9sk1LcmjKFN27eZKF0dLQp1izq97pRds9eeUqp5wB/YOK9tmutv9Va+2ut/X18fLIworAZFy9C795M8fbGLXduZsyYIfdaH1K+fPnYtm0bn3/+OfToASaTbV2tX7oEx4+zPF8+PvjgAynqds61YkV2FC3K0LAwrtvSv7MczppFPQoome55CeCunjtKqRbAu0AnrXWiFfMIW5WSAj17wo0bfBkSwrZt2yhe/EGNOuJ+atasibOzM7pgQWjTxrIcq610Zkqb733+6dOUKVOG+vXrGxxIPA6TycT0RYu4BIx+6SW4etXoSALrFvUwoLxSqoxSyhXoAaxKv4NSqiYwE0tBv2zFLMKWvf8+F7ZsIeazz3CuVg0/Pz+jE9m15cuXU7lyZeKeecbS2zytc5rhtmzhUq5cbAwPp1evXtIS4wDqNGjAkGefZVp8PDv79TM6jsCKRV1rnQK8CqwDjgKLtdaHlVLjlVKd0nabCOQGliil9iulVt3ncMJRrV+PnjCBF4oVo+7kySQnJxudyO4VLVqUY8eOMevvv8HLy3aa4LduZWGpUqSmptK7d2+j04gs8tGsWZTw8uK1lSvRISFGx8nxZPIZYZwLF6BGDWa7uDDwwgWmTZvGkCFDjE7lEJo2bcrJkyeJbNYM15UrLfezc+UyLlBMDPj4MLZpUzYmJrJjxw7jsogs98emTRTq148ncuWCAwfA3d3oSA7HJiafEeK+UlKgRw+O37jB/65coVmzZrz88stGp3IYI0eO5Pz588zz8YHr12H1amMDbdsGwPgPPmDrv2upC4fRoFkznvjuO/SJE8S8847RcXI0KerCGO+9R+K2bfQsVIhcHh7MmzcPJyf555hVWrVqRUBAABOWLSO5WDHjm+C3bOGimxsEBGAymYzNIqyjZUveqFSJ+l98wc1du4xOk2PJb1GR/dasgY8+IqFPH0rVqMF3330nvd2zmFKKKVOmEBwcjEvv3vDrr5YmcCOkpPDPsmWUTU5moiyf69C6fPIJfwJvdOwIZrPRcXIkKeoie61bB127Qo0a5J05k+XLlxMUFJTx+8RDq1+/Pk2aNIE+fSy3OxYtMibIihXMi4oiITWVVq1aGZNBZIvATp0YFRTE7JgYFkhveENIURfZ55dfoFMnIkqXpqmnJ1GxsTKsycqSkpJ4fdYsZhQvblgTfMrkyUxzdqZ+vXpUr17dkAwi+7y/ZAmNvL15af58jm3caHScHEeKusgeK1ZAly5cqlSJjvHxHI+MlHur2cDFxYWjR48y+soVruzaBadOZW+A0FCW7dzJnykpvDViRPaeWxjC2cWFn379lRJKET1kCFy7ZnSkHEWKurC+JUugWzdu1qhBRycnomNjWb16NUWLFjU6mcNTSjFp0iSu3LrFBIDvvsveAF9+yQyTiQrly8ttlhykeN26HF6+nEaRkdC2LcTFGR0px5CiLqxrwQLo0YPEOnXo7u3N3gMHWLRoEbVr1zY6WY5RrVo1BgwYwFdOTpz68kvLmPXsEBUFS5awavBgFi9dKqMbchhTUBB64UI+3rWLMVWqwM2bRkfKEeT/MmE9339v6aQVGMi1+fM5fe4c33zzDR06dDA6WY7zwQcf4JYrFyMSE2H8+Gw5p/7qK8xmM3mGD6datWrZck5hY55+msimTZlw7hyzateG+PiM3yMei8woJ6xj1ix46SXimzbFtHQpbt7eJCYm4ubmZnSyHGvFihX4LV1KuUWL4MgRKF/eeie7eZNVhQszymRizYED+Pr6Wu9cwqYlJycTFBDA+gMHWFGjBh127pQZ5x6BzCgnjPP11zBoEJeaNSPw6lVefP11ACnoBuvcuTPlPv8c7eJC0qhRVj1X8pw5jLh5k1Rvb5mDIIdzcXFh8fbt1PD1pdv+/YQ0aQKJsiCntUhRF1nryy9hyBA2N2hAwLFjHD56lK5duxqdSqRJLVSIZ0qV4qVly2DPHiudJJVZH3zAceDTKVNwcXGxznmE3cidOze/hoZSvnhxzuzeDd27Q1KS0bEckhR1kTW0hk8+4fqbb/LmE0/QfOdOPDw92b59u/R6tiFOTk5U6dSJYGChlSYH+fvHHxkbHU1gpUp07NQp4zeIHMHHx4fwyEj6TZ8Oq1ZxrVs3y6RIIktJUReP76+/uNWiBYwaRdLTTzM3OpqXXnqJvXv3UqtWLaPTiTu899FH1Pf15aVDhzg9b16WH3/qmDHcAL7+6SeZXEj8h6urK7zyCltffRXfVatY27y5TCebxaSoi0eXmkr4iBEEli1L/ZAQ+PprCi5dyp9//smMGTPw9PQ0OqG4B2dnZ3787TdQiq6DB3MzK8cQHzzI+DNn2PrKK1SW2ePEfVQaO5ayxYsTtHUr3zVsKIU9C0lRF4/kxr59vFmyJHUnTuSEszPdhw0jZeBAcHKiQIECRscTGfCtUIEf33yTi7ducfqbb7LkmKdOneL8hx9i8vCgzgcfZMkxhWPy8fFh85EjNC9fnhd37WJkjRqYZYKaLCFD2sTDSUnh5Dvv0HziRM4BLzdtysc//0zefPmMTiYeltlMfI0aeNy6hT58GFxdH7m5PDY2lnoBAXidOcOel15CzZiRxWGFI0pOTua1+vWZuWcPP/j48FxwMLRrZ3QsmyRD2kSWu/bHH1C/Pr4TJ/JU8eLsWLmSrzdtkoJur0wmPD79FP78k3c7dWLw4MGYH6EZ9ObNm3Tp0oVzZ88yTWvUG29YIaxwRC4uLswIC+OXjz6iV/780L4914KCLLMRikciRV1k6K+TJ+lbvToVnnqK66dP47JoET+dO0cD6dls/9q2RTdujNP27Xz77bc8++yzXHuIBTguXbpEYGAg27dvJ9jTkwbt20OFClYMLByNUor2o0bhdPAgF95+m/KrVjGybFkSPv1Uesc/Ainq4r4iT5xgSNu2PFmhAksOHqRf5cqo8HDLGFPp1ewYlEJ9+ikT4uOZ1Lo1K1asoHr16mzdujVTbx8+fDhHjx5l5auv0uP6dZCrdPGoXF3JM3o0Qc8+y6fJyVQfOZLNFSrAzp1GJ7MvWmu7etSuXVsLK/vrL338lVe0E2gX0C96euqz331ndCphTU8/rXXu3PqPNWt0mTJldK5cufSVK1e01lonJyff3i0hIUGvXbtWnzhxQmutdWxsrA4PC9Paz8/ySE01JL5wLOvXrdPlihTRgO4GOmnAAK1jY42OZSggXGeiRsqVugAg4cYNFgwbxsRKlaBMGZ6cMYMplStzZs4cZl29Ssn+/Y2OKKzpo48gPp7669dz8OBBVq1aRb60vhI1atSgXLlylC5dGm9vb9q1a8fw4cMByJ8/P7WvXYOICMtVurTgiCzQslUrIiIjGffOO3hWqYJLcDBUqMA/06ZBaqrR8Wya9H7PwbTW7PvtN4LHjWN+WBhXtMbP2Zl9b72FadAgkEU4cpaBAy0r6x0/DmXKAJZ/I+PHj+fo0aO4u7vj7e1Ny5YtadasGe7/LsrRoQOEhcFff8lCHcI6Dh7k2PPPU2P/frrkysVL7dvTeOhQVL16OeaLZGZ7vxvenP6wD2l+fzxms1mbo6K0XrxYf1Spkga0K+ieRYrojePGaXNCgtERhVGiorR2d9f6uecy/57jx7UGrd97z2qxhNBa6wtRUfp/bdvqvM7OGtAVQU/Ol0/feP11rUNDHf7WD5lsfpcrdUeXkkLUxo1sXLCA37dt4/dz5wg2m2kNHCtQgD/q1qXzhAnkr1nT6KTCFowcCZ99Bvv2QWZmhBsyBGbPhrNnoXBh6+cTOV58fDyLg4P5ZvJk9kZGEmMykSclhZ1FipC7bVuqvvIKqnZth7uCz+yVulWLulKqDTAFMAGztdaf3LHdDZgH1AZigWe11mcedEwp6g+W+Pff3Ni6lQKHDnFp0yYCd+zgRNp/Yx8nJ5oXL87rTz9N/V69oFYtkBW0RHpXrkDZslC0KNSoAcnJlmFFKSn3/vuBA9CzJ8yda3RykQOdOXMG37x5YeVKGr75JjuuXqUIUNfDg7rVq9O4USOeCgqCypXBzufTMLyoK6VMwAmgJRAFhAE9tdZH0u3zClBNaz1YKdUDeFpr/eyDjpsVRd1/wgZibty97F/B3K6Ej275WMe2Kq3h8mXLvcuzZ+HsWZb8/jsRkZEc+ftvjsbFcdJs5gVgppMTqdWq0SMujnr16tGid2+qtmqFk8mUZXHs9nO0ITb5GX73naXjnMkEzs6WL37Ozvf+e65cliv7J54wJis2+hnaGUf4DC9cuMCaRYvYungxuw8e5GR8PG2BtWnbu7u745k/P+VKl6Zo+fIUrVaNCs2bU65GDcByK/pxFiCy9meY2aLu/Nhnur86wCmtdWRaoIVAEHAk3T5BwLi0vy8FpimllLbyPYF7ffD/vq61JjExETc3N5RSXL9+nZs3b5Kamnr7obXGN60T2eXLl4m7Y85iJycnyvj6gtnM3+fOcfPqVcuVTVISJCfjnJpKqYIFITGRv06dIvbCBW5ducKtq1e5dfUq7klJtMyfH65dY9KBAxz+5x9ibt4kNj6ev7WmMrAq7VwjleKM1jzh6UmlkiUJ8vWlRfv28NJLOHl5sdhqn+KDP0eROTb5GQ4YYHnYCZv8DO2MI3yGxYoVY+CbbzLwzTcB+CcmhiuHDsHNm6RGRBAzZQo7oqO5cOHC7bHvLwNfFytGio8PbgcO4OXiQh43N/LkykUeT09eaNqUgV26EO/pydtz5uCVL9/tR+58+ahTrx6Vq1QhJSWFSxfOY/LMhzL9t/Uzuz9Daxb14sC5dM+jgLr320drnaKUugYUAGKsmOu2muePUfHnCfyQEMctrUnQGtOnGg0klCuHm1KMiY5m6h0zbJmAlKJFwWxm1LVrzElM/M/2vMBVpUBrXgOW3HHeEvz/BzMIWH/H9kpKcaRwYcibl3UXL3I0MZECuXJRoGRJAnx8qFOjBrz2GpQuzaarVylStChubm5Z9KkIIYT9y1+wIPmbNAHAqX17No0cCUBCfDx/79nDpd27yXfpEvzzD+ZLl3g3Jobr8fFci4/n+o0bXI+OxjR3Lsydy3VgPhAHpJ9IeSJQ2c2NMy4unL9xw3Je5cTHuQswdcB0brh5ZONPbGHNon6vdow7r8Azsw9KqUFY6h+lSpV6/GRp4l3dIV8RKtxwxc3JhKtywq9oXjycnaFiRXB2pnt0NJWvXcNJKZycnCwPk8kyFaaTEwMuX6ZJXBw4OVk6Zjg54eLsbLkf6erKkKgoOsTF/X+TpbMzHh4eEBgIrq6MiYxkiNmMu7c3uQoUIFeBAuTJnx+efBK4u+DfqbSd3ycSQojs5O7hQelGjSjdqNHt19yA8el3Sk2Fa9cgNhb++YcisbFciY1F37hBwo0b3IiLIy4ujrxKgclE/mvX6L92B9eTbnE9+RbJeQqR7JR1tzofhjWLehRQMt3zEsCF++wTpZRyxnKR+8+dB9Jafwt8C5Z76lkV8LiPL8f7fP6f11Z80v4/z59Ke9xPg7TH/7V3/6F213Ucx5+vLZdjU3Juf+ZypS4AAAUZSURBVMw5dH+MSELqD/dHiIZajQg1MFCEBv3lHyMlhEmCoiUiggT9pTDBYEuEFQws8ldiEdPZWDa9KkMYW4pGITki5Oq7P853eZd33jOu9/s5fs/zAV843+89cF+fN+ee9/3+/JzK5QtkuHSBn0uSerZsGZxzzmiZI8DKblk3Z/sa4OnbHv/f+v09RDyVpXyi3H5gc5JNSVYA1/PRqeAT9gLbutfXAc8s9fl0SZKGasmaelXNAtuB3wEzwGNV9XKSu5OcmN5rJ3BuksPAj4DblirPXGtXrzit7ZqfdVw8a7h41nDxrOHiTUoNffiMJEkTbtxb2pzQRZKkgbCpS5I0EDZ1SZIGwqYuSdJA2NQlSRqIz9zV70n+Dhxp8KvX0tPjayfUtI8frMG0jx+sAViDVuO/oKrWLfSmz1xTbyXJi+PcTjBU0z5+sAbTPn6wBmANJn38Hn6XJGkgbOqSJA2ETX18D7UO0Ni0jx+swbSPH6wBWIOJHr/n1CVJGgj31CVJGgib+piS/CTJS0kOJnkiyXmtM/Utyf1JXu3q8OskX2idqW9Jvpfk5SQfJpnYK2A/bUm2JnktyeEkvcymOEmSPJzknSSHWmdpIcnGJL9PMtN9/m9unalvSc5M8kKSv3Q1uKt1pvl4+H1MSc6uqn91r38IXFRVNzWO1ask32Q05/1skvsAqmpH41i9SvIl4EPgQeDWqhr8lIFJlgOvA98AjgH7gRuq6pWmwXqU5DLgOPCLqvpy6zx9S7IeWF9VB5KcBfwZuHbKPgMBVlXV8SRnAH8Ebq6qfY2jncQ99TGdaOidVcDU/TdUVU9U1Wy3ug84v2WeFqpqpqpea52jZ1uAw1X1RlW9DzwKXNM4U6+q6jngn61ztFJVb1XVge71e8AMsKFtqn7VyPFu9Yxumbg+YFM/DUnuSXIUuBG4o3Wexn4A/LZ1CPViA3B0zvoxpuwLXR9JciHwVeD5tkn6l2R5koPAO8CTVTVxNbCpz5HkqSSH5lmuAaiq26tqI7AL2N427dJYqAbde24HZhnVYXDGqcGUyTzbJm4PRUsvyWpgD3DL/x29nApV9UFVfYXRUcotSSbuVMznWgeYJFV11Zhv3Q08Dty5hHGaWKgGSbYB3wGurIFekHEan4NpcQzYOGf9fODNRlnUSHceeQ+wq6p+1TpPS1X1bpJnga3ARF086Z76mJJsnrN6NfBqqyytJNkK7ACurqp/t86j3uwHNifZlGQFcD2wt3Em9ai7SGwnMFNVD7TO00KSdSfu+EmyEriKCewDXv0+piR7gC8yuvL5CHBTVf2tbap+JTkMfB74R7dp3xTeAfBd4OfAOuBd4GBVfattqqWX5NvAz4DlwMNVdU/jSL1K8kvg64xm6HobuLOqdjYN1aMklwJ/AP7K6DsQ4MdV9Zt2qfqV5GLgEUZ/A8uAx6rq7rapPs6mLknSQHj4XZKkgbCpS5I0EDZ1SZIGwqYuSdJA2NQlSRoIm7okSQNhU5ckaSBs6pI+UZJLkrzUzSe9qptLeuKeeS3Jh89IGkOSnwJnAiuBY1V1b+NIkuZhU5e0oO6Z7/uB/wBfq6oPGkeSNA8Pv0saxxpgNXAWoz12SRPIPXVJC0qyF3gU2ASsr6rtjSNJmofzqUv6REm+D8xW1e4ky4E/Jbmiqp5pnU3SydxTlyRpIDynLknSQNjUJUkaCJu6JEkDYVOXJGkgbOqSJA2ETV2SpIGwqUuSNBA2dUmSBuK/tRYNh9x86LcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe3169e7128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def Uniform_nodes(a, b, N):\n",
    "    '''Computes uniformly spaced nodes'''\n",
    "    h = (b-a) / N\n",
    "    return [a + i*h for i in range(N+1)]\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def init_mesh(Omega, N, n, **kwargs):\n",
    "    ''' Function initializes a 2D array which defines N+1 global \n",
    "        1D finite elements, with n+1 nodes inside each element\n",
    "        in the domain Omega = (x0, x1) '''\n",
    "    node_spacing = kwargs.get('node_spacing', 'uniform')\n",
    "    if (node_spacing == 'uniform'):\n",
    "        nodes = Uniform_nodes(Omega[0], Omega[1], n + N*n)\n",
    "    elif (node_spacing == 'chebyshev'):\n",
    "        nodes = Chebyshev_nodes(Omega[0], Omega[1], n + N*n)\n",
    "    # Initialize elements\n",
    "    elements = [[i-j%(N+1) for i in range(j*(n+1), (j+1)*(n+1))] for j in range(N+1)]\n",
    "    return (nodes, elements)\n",
    "\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def assemble(nodes, elements, X, phi, f):\n",
    "    ''' Function to assemble the global linear system of \n",
    "    equations A*x=b for a given mesh which is defined \n",
    "    through nodes and elements\n",
    "    '''\n",
    "    N_n, N_e = len(nodes), len(elements)\n",
    "    A = np.zeros((N_n, N_n))\n",
    "    b = np.zeros(N_n)\n",
    "    \n",
    "    for e in range(N_e):\n",
    "        Omega_e = (nodes[elements[e][0]], nodes[elements[e][-1]])\n",
    "        A_e = element_matrix(X, phi, Omega_e)\n",
    "        b_e = element_rhs(X, f, phi, Omega_e)\n",
    "        \n",
    "        for r in range(len(elements[e])):\n",
    "            for s in range(len(elements[e])):\n",
    "                A[elements[e][r], elements[e][s]] += A_e[r, s]\n",
    "            b[elements[e][r]] += b_e[r]\n",
    "        \n",
    "    return (A, b)\n",
    "\n",
    "def u_glob(nodes, elements, C, resolution):\n",
    "    x_patches = []\n",
    "    u_patches = []\n",
    "    for e in elements:\n",
    "        (x_L, x_R) = (nodes[e[0]], nodes[e[-1]])\n",
    "        d = len(e)-1\n",
    "        X = np.linspace(-1, 1, resolution)\n",
    "        x = 0.5*(x_L + x_R) + 0.5*(x_R - x_L)*X\n",
    "        x_patches.append(x)\n",
    "        \n",
    "        u_element = np.zeros(resolution)\n",
    "        for r in range(len(e)):\n",
    "            i = e[r] # global node indice\n",
    "            u_element += C[i] * phi_r(r, X, d)\n",
    "        u_patches.append(u_element)\n",
    "    x = np.concatenate(x_patches)\n",
    "    u = np.concatenate(u_patches)\n",
    "    return (x, u)\n",
    "    \n",
    "# Test assembly of linear equation system\n",
    "#-------------------------------------------------------\n",
    "from numpy import tanh, pi, sin, cos, exp\n",
    "d = 3\n",
    "N = 5\n",
    "\n",
    "Omega = (-pi, pi)\n",
    "f = lambda x: -tanh(x) * sin(x-pi) * exp(1. - 0.75* (x-0.25)**2)\n",
    "\n",
    "#Omega = (0, 1)\n",
    "#f = lambda x: x*(1-x)**8\n",
    "\n",
    "\n",
    "elem_resolution = 40\n",
    "x   = np.linspace(Omega[0], Omega[1], elem_resolution)\n",
    "X   = np.linspace(-1, 1, d+1)\n",
    "phi = basis(X, d)\n",
    "\n",
    "# Initialize mesh\n",
    "(nodes, elements) = init_mesh(Omega, N-1, d)\n",
    "# Assemble the global linear equation system \n",
    "(A, b) = assemble(nodes, elements, X, phi, f)\n",
    "# Solve the global linear equation system\n",
    "C = np.linalg.solve(A, b)\n",
    "# Compute the approximation u(x) to f(x)\n",
    "(x_u, u) = u_glob(nodes, elements, C, 40)\n",
    "\n",
    "# Plot results\n",
    "#-------------------------------------------------------\n",
    "element_edges = [nodes[elements[i][0]] for i in range(0,N)]\n",
    "element_edges.append(nodes[-1])\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_size_inches(8,4)\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "e_plt = ax.plot(element_edges, [0]*len(element_edges), ls='None', marker='s')\n",
    "\n",
    "f_plt = ax.plot(x, f(x), c='r', ls='-', label='$f(x)$')\n",
    "u_plt = ax.plot(x_u, u, c='k', ls='--', label='$u(x)$')\n",
    "ax.set_xlabel('x') \n",
    "ax.set_ylabel('f(x), u(x)')\n",
    "ax.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finite Element Method in 1D using Collocation Method\n",
    "\n",
    " \n",
    "For the collocating finite element method, the entries in the element matrices and element right hand sides are not evaluated through integration. \n",
    "Instead the are calculated through:\n",
    "$$\n",
    "\\tilde{A}^{(e)}_{r,s}  = \\phi_r(x_s) \\text{det}(J)\n",
    "$$\n",
    "$$\n",
    "\\tilde{b}^{(e)}_{r} =  f(x_r) \\text{det}(J)\n",
    "$$\n",
    "since $\\phi_r(x_s) = 1$ for $r = s$ and $0$ for all other nodes.\n",
    "\n",
    "Therefore, we do not need no longer the trapezoidal rule for the numerical integration of terms like $\\int_{\\Omega^{(e)}} \\phi_{q(e,r)} \\phi_{q(e,s)} \\text{d}x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAEKCAYAAAALjMzdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xd0FNX7x/H3zaYXeqgBglIDUiO9d0IJoFIVpKMgoIhSFEG+FhDBLr2oKE2kSJcSQGqQHkAxUkINnUDq7v39sYEfIpAAmUx287zO2UN2Z3buZwmcZ+fOnXuV1hohhBBCOD4XswMIIYQQIm1IURdCCCGchBR1IYQQwklIURdCCCGchBR1IYQQwklIURdCCCGchGFFXSk1Qyl1QSl18AHbOyul9ic/tiqlyhmVRQghhMgMjDxTnwU0fcj2f4A6WuuywBhgioFZhBBCCKfnatSBtdablFKBD9m+9a6n24EAo7IIIYQQmYFhRf0R9QBWPmijUqo30BvAx8enUsmSJdMrlxBCCGG63bt3X9Ra+6e0n+lFXSlVD3tRr/mgfbTWU0jung8ODtbh4eHplE4IIYQwn1LqRGr2M7WoK6XKAtOAZlrrS2ZmEUIIIRydabe0KaUKAYuAl7TWf5qVQwghhHAWhp2pK6V+AuoCuZRSUcB7gBuA1noSMBLICXyjlAJI0loHG5VHCCGEcHZGjn7vmML2nkBPo9oXQgjhXBITE4mKiiIuLs7sKIbx9PQkICAANze3x3q/6QPlhBBCiNSIiorCz8+PwMBAknt4nYrWmkuXLhEVFUWRIkUe6xgyTawQQgiHEBcXR86cOZ2yoAMopciZM+cT9URIURdCCOEwnLWg3/akn0+634UQ4hGsXr2aY3/8QdmSJakRGoqLi5wbiYxD/jUKIcRD/PHHHzRq1IizZ89CUhIHvv6a/sOHU7ttWwoHBDB16lSsVqvZMYUApKgLIcR9aa0ZNWoUzz77LAcOHOD46tVQtSqvL1vG2bp1mevlReD16/Tu3ZtmzZqhtTY7shBS1IUQ4l5JSUl06dKF0aNH82KHDhxp145qPXtCVBSW+fPJu3497RcuZNPNm8yuXJluL7/s9Nd6xb8NHjyYoKAgevXqRZ06dR7aW5OQkEDt2rVJSkoyPJcUdSGEuMe4ceP44Ycf+F/XrszaupVsX34J3bvD4cPwwgugFISEoMaOpcvOnXT85x8Ali9fTlRUlMnphdEiIyP5/fffiYiIoHz58rRt2xaLxfLA/d3d3WnQoAHz5s0zPJsMlBNCiHsM6NiRwsuX03n2bChRAsLCoHbt/+44ZAjs3w/vvMO1IkV4qX9/ChUqxLZt2/Dy8kr/4JnJoEGwd2/aHrN8efjss4fucvToURo2bEhSUhIVKlQA4JdffrmzvV69egwfPpxGjRrxzjvvcP36db744gtat27NsGHD6Ny5c9pmvocUdSGESHbwwAGKbN2K7/DhdL5xA957D4YNAw+P+79BKZg6FY4eJWufPnw/diwt+vVjwIABTJ06NX3Di3RRokQJunbtSmBgIF26dKFQoUIEBgbe2T569GhGjhzJhQsX2LNnD0uXLgWgTJky7Nq1y/B8UtSFEAK4fPYsTZ59lmfj41lcsyZMngxBQSm/0csLFi+G4GCaf/opwwYN4qPPPqNZs2a0bdvW+OCZVQpn1EY6cOAAoaGhXLx4kWzZsv1rW+3atdFaM2HCBDZu3HinW95iseDu7s6NGzfw8/MzLJtcUxdCCGDoiy9yPj6ekQMH2rvbU1PQbytQAH75BaKieH/vXiqUL0///v25du2acYGFaQ4dOkTp0qXx8vL6z+xvBw4c4OzZs3h4ePyneMfHx+Pp6WloNinqQohM7/fff2fq+vUM8vOj4oQJ8DgTylStClOm4LpxI9NKlGDUqFGGnpEJc9y4cQM3Nze8vb3Jnj07Vqv1TmE/e/YsnTt3ZsmSJfj4+LB69eo777t06RL+/v6PvVBLaklRF0JkajabjYGvvEJBYNSrrz5eQb+ta1d44w0qzptHbxcXmW3OCR08eJAyZcrced64cWO2bNnCrVu3aNu2LZ9++imlSpXi3XffZdSoUXf227BhAyEhIYbnk2vqQohM7fLly3jfuMEHgG/v3k9+wLFj4eBBePVVZvz9NxtPn+a777578uOKDKFatWosWLDgzvP+/fszYcIEGjZsyLZt2+68Xrt27X89//HHH/noo48MzydfI4UQmVqunDkJs1h4sXZteOqpJz+gqyvMnQuBgVz+8ku+//57wsLCnvy4IkOqUKEC9erVS3HymdatW1OiRAnD80hRF0JkWjt37uTMkiWov/9Gde+edgfOnh2WLqWfiwt5XF0ZM3p02h1bZDjdu3dPcfKZLl26pEsWKepCiEzJarXy0ksv8Vzv3uDrC88/n7YNlCyJ1+zZDElKYt2GDf/qihXCKFLUhRCZ0s8//8yff/7JGzduQLt24OOT9o20bk2fokXJ5erKhx9+mPbHF+IeUtSFEJnShAkTKJYnD23j4qBbN2MaUQrf3r35PCmJ/q1aGdOGEHeRoi6EyHR27NjBjh07GODnh6VYMahRw7jGunalk5sbTY4cMa4NIZJJURdCZDo7duwge9asdD12DF5+2T6Hu1Fy54bQUM7MnMnQIUO4evWqcW2JTE+KuhAi0xkwYAAn+vTBz8UF0mNUcq9enLlyhbHjx/Pjjz8a357ItKSoCyEyldjYWLDZ8Js7Fxo1goAA4xtt2JBKhQtT3teXqVOnorU2vk2RKUlRF0JkGlprgoODef2FF+DkSeMGyN3LxQXVsye9YmLYu3cvu3fvTp92RbqKjY2lTp06KU5EU7t2bZKSkgzJIEVdCJFpbNu2jYiICMqcPg3ZskFoaPo13q0bnZXCy9WVadOmpV+7It3MmDGDtm3bpjgRTYMGDZg3b54hGQwr6kqpGUqpC0qpgw/YrpRSXyiljiml9iulKhqVRQghAGbOnImPjw/t9u6FTp3A4GUw/6VAAbI2b86Lbm7oh5zJidSrW7fufx7ffPMNALdu3brv9lmzZgFw8eLF/2xLrapVq3L8+HEATp8+TXBwMABz5swh9K4vivXq1WPt2rUAvPPOOwwYMACA1q1bM2fOnCf89Pdn5IIus4CvgAetZNAMKJb8qAJ8m/ynEEKkufj4eBYuXEjbMmXw27Ej/bre79arF5N//RXVvHn6ty3ShNaakydPUrhwYQD279/PM888Q0JCApGRkQQGBt7Zd/To0YwcOZILFy6wZ88eli5dCkCZMmXYtWuXIfkMK+pa601KqcCH7BIKfKftI0a2K6WyKaXyaa3PGpVJCJF5rVixgqtXr9Lp6lUoUwYqVUr/ECEhqPz5YepUzlapQr58+dI/gxPZuHHjA7d5e3s/dHuuXLkeuv1Bjh07RpEiRVDJt0HeLuoXL14kW7Zs/9q3du3aaK2ZMGECGzduvNMtb7FYcHd358aNG/j5+T1yhocx85p6AeDUXc+jkl/7D6VUb6VUuFIqPDo6Ol3CCSGcS40aNfhqxAgaHj1qP0s38t70B3F1hW7dGL9yJYGBgVy5ciX9M4gncuDAAZ555pk7z8PDwylbtixeXl7ExcX9Z9+zZ8/i4eHxn+IdHx+PpwGXf8ws6vf7H3Xf+zy01lO01sFa62B/f3+DYwkhnFHu3LnpFx+Pq6srvPiieUF69KCO1iQkJLBkyRLzcojHcvnyZby8vAA4fPgwy5cv55lnniF79uxYrdY7hf3s2bN07tyZJUuW4OPjw+rVq+8c49KlS/j7++Pm5pbm+cws6lFAwbueBwBnTMoihHBiYWFhzJg6lcTvvoPmze2zvJmlSBGCGzaksMXCgvnzzcshHkuTJk1Yt24d7dq1Y8GCBeTMmZM8efIA0LhxY7Zs2cKtW7do27Ytn376KaVKleLdd99l1KhRd46xYcMGQkJCjAmotTbsAQQCBx+wrTmwEvsZe1VgZ2qOWalSJS2EEI+ibdu2Ol+OHNoKWi9ebHYcrefP14NBu7m66itXrpidxmFERESYHeGh/vjjD/3iiy+muF+bNm30kSNHHrj9fp8TCNepqJFG3tL2E7ANKKGUilJK9VBK9VVK9U3eZQUQCRwDpgKvGpVFCJF53bx5k5UrV9I2WzZccucGo86QHkVoKC9ky0ZiUtKdEdHC8VWoUIF69eqlOPlM69atKVGihCEZjBz93jGF7RroZ1T7QggBsHLlSmJjY3n++HEYNAgMuI75yNzdqdyjBz9MnEhIFbmT15l07979odvd3d3pYuB6AzKjnBDCqS1ZsoScPj7UstnMuTf9AVSvXnS22ci+eLHZUYQTkaIuhHBqZ8+eJcTNDUtwsP3+9IyiRAlu1ajBF+PHs23rVrPTOAzt5IvhPOnnk6IuhHBqv40bx/SrVzPUWfptrj178u7Fi0z/4AOzozgET09PLl265LSFXWvNpUuXnuj+dSOniRVCCFNprVHz5tnvB+740GE+pnBv355mvXuzbN06rFbrQxcCERAQEEBUVBTOPAmZp6cnAU+wHLAUdSGE06pbty61jh3jf5UrQ/bsZsf5Ly8vWtWty7y1a9m9fj2VGzUyO1GG5ubmRpEiRcyOkaFJ97sQwimdOXOGTZs24XvuHNSqZXacB2o0fDgKWDNxotlRhBOQoi6EcEorV64EoLnNlqGLun/dujzr48Pp7dvBSa8Vi/QjRV0I4ZSWL19OQJYslAGoXt3sOA+1Zfx4vr1yBbZtMzuKcHBS1IUQTichIYG1a9fS3M8PVa4c3LMkZkbj9uKL4OsLs2aZHUU4OCnqQginExsby6ABA+h06RLUrGl2nJT5+tIma1YG/fADPGSKUSFSIkVdCOF0smbNypjWrakdF5ehr6ffzZYvH0tjY9G//252FOHApKgLIZzOtm3buLVunf2JgxT1xh068A/wt3TBiycgRV0I4VSuXLlCjRo1+GT2bHjqKcif3+xIqdIkNBSA1UuXyih48dikqAshnEpYWBhaa+qfPesY19OTFS1alKf8/Vl96RLs3292HOGgpKgLIZzKhg0b8PL0pPK1aw7T9X7bgAEDaKgU/PKL2VGEg1KONjF+cHCwDg8PNzuGECKDKlu2LHmSklh7+DAcOQIlSpgd6dHUrg3XrsG+fWYnERmIUmq31jo4pf3kTF0I4TSio6M5cOAA9Tw8IHduKF7c7EiP7FrTphzevx8iI82OIhyQFHUhhNPIkSMHO3fu5KXb96crZXakRxa6dCkvgnTBi8ciRV0I4TQsFgvP5s9PwVOnHGqQ3N0atmjBH0D0ggVmRxEOSIq6EMJpjB49mq3Tp9ufONgguduaNGkCwNodO+D8eZPTCEcjRV0I4RTOnDnDqFGj2LpqlX0e9fLlzY70WCpWrEjObNlYBbB0qdlxhIORoi6EcAobN24EoH50NFSrBq6u5gZ6TBaLhSYhIax0ccG6aJHZcYSDkaIuhHAKmzdvJoufH+WOHXPYrvfbRowYQViXLrisWwfXr5sdRzgQKepCCKewdetWqhUrhgUcdpDcbUFBQQT16IFKTIQVK8yOIxyIFHUhhMOLj4/n0qVLVPf2Bjc3qFLF7EhPbH1sLCO9vWHxYrOjCAdiaFFXSjVVSh1VSh1TSg29z/ZCSqkNSqk9Sqn9SqkQI/MIIZyTh4cHp06dYmhSElSqBN7eZkd6Ytt37WLMrVuc/fVXiI83O45wEIYVdaWUBfgaaAYEAR2VUkH37PYOMF9rXQHoAHxjVB4hhHNTcXG4797t8NfTb2vevDkAq27ehNvLyAqRAiPP1CsDx7TWkVrrBGAuEHrPPhrIkvxzVuCMgXmEEE6qe/fujOzbFxITHf56+m1ly5alQIECLHd1ldnlRKoZWdQLAKfueh6V/NrdRgEvKqWigBXAa/c7kFKqt1IqXCkVHh0dbURWIYSDslqtLFy4kItHjthfqFHD3EBpRClFSEgIa4DExYvBajU7knAARhb1+026fO+ScB2BWVrrACAE+F4p9Z9MWuspWutgrXWwv7+/AVGFEI7q0KFD3Lhxg+qxsVC6NOTMaXakNNO8eXP8c+TgxMWLsHWr2XGEAzCyqEcBBe96HsB/u9d7APMBtNbbAE8gl4GZhBBOZmtysaseGek019Nva9myJcf+/JOi7u4yCl6kipFFfRdQTClVRCnljn0g3L1zHp4EGgAopUphL+rSvy6ESLWtW7eSN1cuity86XRF3cXFBZU1K7YGDezX1fW9nZ1C/JthRV1rnQT0B1YDh7GPcj+klHpfKdUqebfBQC+l1D7gJ+BlreVfrRAi9YoUKUKnMmXs1/ucZJDc3VatWkXuzZv5559/YP9+s+OIDM7QyZG11iuwD4C7+7WRd/0cATjHqBYhhClGjx4NL7wAhQrZH07m6aef5lJMDMuB/r/8AuXKmR1JZGAyo5wQwmHdvHkTm9UKmzc7Xdf7bcWKFaNYsWIsz55dbm0TKZKiLoRwWGPGjCF/3rxYz5932qIOEBISwoYbN7i5fz9ERpodR2RgUtSFEA5r165dFPT1tS/i4sRFvVWrVsQnJbEWZBS8eCgp6kIIh2Sz2di9ezfBHh6QIweULGl2JMPUqlWLwYMHU7R4cemCFw9l6EA5IYQwyrFjx7h27Zq9qNesCS7Oe47i5ubG+PHjwc8PRo+G8+chTx6zY4kMyHn/FwghnFp4eDgAz1644NRd77dZrVa2FC7MMa1h6b1TfghhJ0VdCOGQypQpw4i2bQmCTFHUb968Sf3evZmUJYt0wYsHkqIuhHBIZcuW5X8BAbh6eUGFCmbHMVyWLFmoX78+S1xc0L/9Btevmx1JZEBS1IUQDsdqtbJ582ZubtwIVauCu7vZkdJFaGgox65e5XBiIqxaZXYckQFJURdCOJzDhw9Tu3ZtFh04kCm63m9r1co+w/YSHx/pghf3JUVdCOFw7gyS0zpTFfUCBQoQHBzMb1mzwvLlEB9vdiSRwUhRF0I4nPDwcHzd3Snu4mLvfs9EFi5cyMqvv4YbN2DjRrPjiAxGiroQwuHs2rWLSl5euFSsCL6+ZsdJV4ULF8a9aVOQLnhxH1LUhRAOJSEhgX379hEcE+OUS62mxmeTJtE/b15YsgRsNrPjiAxEiroQwqFYLBbWT5pEL6sVqlQxO44pTpw4wbSTJ4k5dw527jQ7jshAUizqSilPpdTzSqnPlVILlFLfKaXeUkqVTo+AQghxN4vFQvWEBEoAPPus2XFMERoaSnxiIqtdXGSBF/EvDy3qSqlRwO9ANWAHMBmYDyQBHyul1iqlyhodUgghblu6dCnLFi60L+Ly1FNmxzFFzZo1yZkzJ7/kzm2/rq612ZFEBpHSgi67tNajHrBtglIqN1AobSMJIcSDjR07Fsu+fbSsWROUMjuOKVxdXWndujXz58wh/tw5PI4cgVKlzI4lMoCHnqlrrZeDvQv+3m1KqVxa6wta63CjwgkhxN1sNhv79++n/M2bULmy2XFM1bFjR5o1asRlkC54cUdqB8rtUkrduRlUKfUcsNWYSEIIcX+RkZHExMRQDjLt9fTbGjRowLylS8lXubLc2ibuSO166p2AGUqpjUB+ICdQ36hQQghxP3v37gWgPGT6on7bsVq1KPTpp7hHRUFAgNlxhMlSdaautT4AfAD0BeoB/bXWUUYGE0KIex06dAiLUpQOCIC8ec2OY7r169dT7NNPWQ+yxroAUlnUlVLTgUFAWaAbsEwp1c/IYEIIca+RI0dyvGBBPDPp/en3qlGjBn5+fizMmlWuqwsg9dfUDwL1tNb/aK1XA1WBisbFEkKI/1KXLxNw8qR0vSfz8PCgVatW/BIfT+L69XDlitmRhMlS2/0+Uev/vxFSa31Na93DuFhCCPFvly9fplvHjuwGKep3ee6557gcF0eY1QorVpgdR5gspclnlimlWiql3O6z7Sml1PtKqe7GxRNCCLu9e/cya+1a+y1clSqZHSfDaNq0KT4+Piz09pYueJHimXovoBZwRCm1Sym1Qim1XikViX12ud1a6xkPerNSqqlS6qhS6phSaugD9mmnlIpQSh1SSv342J9ECOHU9u3bB0C5okUha1aT02QcXl5e/Prrr3zUvj2sXAmxsWZHEiZ66C1tWutzwFvAW0qpQCAfEAv8qbW+9bD3KqUswNdAIyAK+73uS7XWEXftUwwYBtTQWl9JnqFOCCH+Y+/eveRzcSF39epmR8lw6tatC/HxMHMmrFsHLVqYHUmYJNWrtGmtj2utt2mt96ZU0JNVBo5prSO11gnAXCD0nn16AV9rra8kt3EhtXmEEJnLvvBwyttscj39Aab8/TdfeHhIF3wml9pb2m4opa4nP+KUUlal1PUU3lYAOHXX86jk1+5WHCiulPpdKbVdKdX0Ae33VkqFK6XCo6OjUxNZCOFEbDYbtpgY+y03mXx62AdZvW4dHymFdckSsFrNjiNMktrR735a6yzJD0/gOeCrFN52v5UW7l1KyBUoBtQFOgLTlFLZ7tP+FK11sNY62N/fPzWRhRBOxMXFhf0dOzLG1RXKlTM7TobUvn17zsXFseniRdgqs3hnVqnufr+b1noxKU8TGwUUvOt5AHDmPvss0Vonaq3/AY5iL/JCCPFvO3eiypcHDw+zk2RILVq0wNfXlx9ljfVMLbXd723vejyvlPqY/55132sXUEwpVUQp5Q50AO6dx3Ax9mlnUUrlwt4dH/lIn0AI4fQ++vBDQjdvRgcHmx0lw/L29qZ169b8bLEQL2usZ1qpXdCl5V0/JwHH+e+gt3/RWicppfoDqwELMENrfUgp9T4QrrVemrytsVIqArACQ7TWlx7xMwhhmkuXLhFz7Bhex4/j36IFysfH7EhOadOqVZxLSkLJ9LAP1blzZ07s2MG5v/6i8IEDULas2ZFEOlPawb7NBQcH6/BwWcJdmCM6Opr58+fz6quvopSiV5MmTFuzBgAfoFru3IQ0bky/iRNxz5XL3LBOpGDOnNS9fJnvDx6E0qXNjpOxnTsH+fPDqFEwcqTZaUQaUUrt1lqn2FX1WNfUhchszp8/T79+/QgICKB///78deQIvPcefdesYXrBgnzevTvdypTh3OXLTP7hB9wKFIBWrUiaMUPm435CV69eJeryZUq7u0PJkmbHyfjy5uVScDCxixaZnUSYILXd70JkSlprZs+ezcCBA7l16xY9evRgQJcuFH/jDVi1ikovv0ylb74BLy/7G6xWrq1di1q9musLFhC8bBkje/Wic4MGqOefh9atIbfMsfQoIiLs81WVKV4cLBaT02R8Bw8epMLu3cy22eh04gQULmx2JJGO5ExdiIe4fPkyw4YNo3z58hw8eJBJvXsT1LkzrF8PkyfDjBn/X9ABLBayNm0KEydyY/t2/MuW5SWbjW7btxPXpw8UKADLl5v3gRxRYiL1laJsjRpmJ3EIQUFB5M2dmx9BRsFnQo9V1JVSoUopGbEinNb169fRWpMzZ0527NjBhg0bKLF1K1Svbp/YY/Nm6N0b1P2mY7ArEBDApj/+4L333mP2jRvUL1eOK8WKQd++cONGOn4ax1bdx4d1WlOoQQOzozgEFxcXOr70EquBSwsWmB1HpLPHPVOvAryjlFqZlmGEyAjOnz9PjRo1ePvttwEolDs3Ln37QvfuULMm7N6d6lnNLBYLo0aNYsGCBew+fJi3ixeH06fh3XeN/AhOxbZjh/0HmR421Tp27EgSsHDrVrgkNxRlJjL6XYi7XL16lTp16nDs2DGWLVtG/aefhuefh/BwGDYMxox57Ou6mzdvpmzZsmQdPhwmTYIdO0Duu05RQR8fugL/i4l5aM+I+H9aa4Keeorcx48TNn26/QupcGhpOvpdKZVbKdVGKdVPKdVdKVVZKSXX44VTiY+PJzQ0lMOHD7N48WLq22z2dbv//NN+bfLDD59ooFatWrXImjUrN0eMYKS3N4m9ekFSUhp+AucTHR1N1K1b5CxcWAr6I1BKMWX2bGYULAg//WR2HJGOHlqYlVL1lFKrgeVAM+xLrwYB7wAHlFKjlVJZjI8phPEGDBjApk2bmD17No3y5oXmzSFfPvtZeuhD51p6JOvCwxkTE8OgvXvhiy/S7LjO6NCuXQCUrlTJ5CSOp1bt2jzdrZt9KdYz987QLZxVSmfbIUAvrfWzWuveWut3tNZvaq1bAeWAPdjXSxfC4bVs2ZL//e9/dHzuOejaFbJlgw0boFjaLkfQqlUrhrz5Jt8AM4cNgxMn0vT4zuTg6tUAlGnc2OQkjmlT8eL01xr9449mRxHp5KFFXWs9RGt98gHbkrTWi7XWPxsTTYj0kZiYCNgXxBgxYgR88AHs2QNTp4JBs8J99PHH1K9enf4JCRzu0kXm6X6AQzt3kg3I16SJ2VEc0pGbN/kaCJ861ewoIp2k9pr690qprHc9D1RKrTMulhDpIy4ujsqVKzNp0iT7C7t324t6ly7QqpVh7VosFn5YuBAfX1/6btoEP8t34/uprTVvZMuGkgl7Hku7du3wdHNj9p9/wqFDZscR6SC1g922ADuUUiFKqV7AGuAz42IJkT5Gjx7N3r17CQwMhLg4ezHPmxc+/9zwtvPly8fi5cuZU7o0DBgA164Z3qaj6Xj+PO9K1/tjy5YtG21atOBHIH72bLPjiHSQqqKutZ4M9ASWAO8DtbXWy4wMJoTRdu3axbhx4+jRowdNmza1L34REQHTp9uvp6eD6rVrEzBrFvrcOa688Ua6tOkoYk+e5PTx47Lc6hN6uW9frgDLZs4Em83sOMJgqe1+fwmYAXQBZgErlFLlDMwlhKFsNhv9+vUjT548fPrpp7B1K4wfD336QHpfvw0Opn3RorScMQO9dWv6tp2Bbf7uOwKATR4eZkdxaA0aNKB2yZJYL16E3383O44wWGq7358Damqtf9JaDwP6Yi/uQjikHTt2sHv3bsaNG0dWV1f7aPfAQPjkE1PyNBk4kN+BOe3bQ/LAvcxub1gYAGXS8HbCzMhisRAWHk57Hx/44Qez4wiDpbb7vbXW+sJdz3dinypWCIdUrVo1IiIi6Ny5MwwdCseOwcyZ4OdnSp5ur7zCs8WK8VZUFDc++siUDBnNH4fFHR89AAAgAElEQVQOUdjNzT7xjHgyPj7Et2zJsZ9+gvh4s9MIA6U0+cw7Sqkc99umtU5QStVXSrUwJpoQxjiRfF94iRIlUBs2wFdfwaBBUKeOaZlcXFz48vvvOQuMGTMGIiNNy5IhaM3uCxeolD+/2UmcRuhff9Hmxg30ihVmRxEGSulM/QCwTCm1Tin1iVLqLaXUyORb3A4ALYEdxscUIm0cPnyYokWLMn36dLh+Hbp1g+LF7VPAmqxKlSp0b9+eBVYrsX36ZOp7168dOMAxq5WKZcuaHcVptOnenYPAH199ZXYUYaCUJp9ZorWugf0a+iHAAlwHfgAqa61f11pHGx9TiLTx7rvv4uXlRatWrWDwYIiKgtmz/70muokmTJ5MxCef4PXbbzB3rtlxTOO6bx9zgLZdupgdxWm079QJD4uFWRs3wtWrZscRBkmp+/375B9DtNaztNYfaa0/01qv1lrHpkM+IdLMvn37+Pnnnxk0aBD+u3bBtGnw1ltQtarZ0e7ImjUrXoMGkVCuHBffeSfTnq377N9PJ3d3Shk4AVBmky1bNto0aMAcm41YWeTFaaXU/V5JKVUY6K6Uyq6UynH3Iz0CCpFWRo0aRdasWXn95ZehZ08oUwZGjTI71n/YlCL40iX6R0bCli1mxzHF+jVr2F+iBLi7mx3FqfR++22uAMulC95ppVTUJwGrgJLA7nsesqi5cBjR0dGsX7+e119/newjR0J0NHz3HWTAe6BdXFxo2aED84AD48ebHSf9Wa30O3CAd2NizE7idOrWq8fWnj15LiICTp0yO44wQErX1L/QWpcCZmitn9JaF7nr8VQ6ZRTiifn7+xMZGckbJUrAnDnw7rtQoYLZsR5o8LBhZHFzY9Svv2a66WNj/viDo1pTKQP/fhyVUopqQ4eiQNZZd1KpvU/9FaODCGGUmJgYtNbk9PHBb/hwKF0ahg0zO9ZD5ciRg9e7dmWRzca+THa2vnfRIjRQsWFDs6M4p6efZnC+fIzIZP+uMovUzignhMPq2bMn9erVQ3/6Kfzzj32xFjc3s2OlaODYsfi4uPDjtGlmR0lXf2zeDEAlGSRnmHMFCvB1dDS3dsgdyc7G0KKulGqqlDqqlDqmlBr6kP2eV0pppZSs3CDS1IkTJ1i4cCHPliyJ+ugjaNMGGjQwO1aqZM+Rg/C33+bjc+dg/36z46Sb3YcPk9fdnXwFCpgdxWn1efddrgHzMuBAUfFkDCvqSikL8DXQDAgCOiqlgu6znx8wAJnERhjgyy+/BOC18+chKQk+/dTkRI+m5ODBKHd3rFOnmh0lfcTH89n166zq3NnsJE6tVsuWlPL1ZdK6dbJym5Mx8ky9MnBMax2ptU4A5gL3W5lhDDAOiDMwi8iErl+/ztSpU3mhXj0KLV4Mb74JRYqYHevR5MzJnAoVKP7NN9y8fNnsNMbbt4/sSUmUayGzTxtJKUWftm3ZmZjI3unTzY4j0pCRRb0AcPc9E1HJr92hlKoAFNRa//qwAymleiulwpVS4dHRMoGdSJ3vv/+e69ev83pUFOTPb1+4xQEV6dCBSJuNWW++aXYUwx1etIjRwDlH+/LlgLp8+CF93dzwW7PG7CgiDRlZ1NV9XrszPZZSygWYCAxO6UBa6yla62CtdbC/v38aRhTOrGfPnix97TUqHzkC48aBr6/ZkR5Ltddeo4q7OxPnzsVqtZodx1C/rVjBKCBJ/p8bLnuBAnzbqRNPr1kDcdJR6iyMLOpRQMG7ngcAZ+567geUATYqpY4DVYGlMlhOpBWPuDhazpsH1atDp05mx3lsymLhjVat+Ds2lpWzZpkdxzg2G9uPHqWAtzcBAQFmp8kcXnyRHdev81sGWNBIpA0ji/ouoJhSqohSyh3oACy9vVFrfU1rnUtrHai1DgS2A6201jJTnXhir7zyCtPatbPPHPf556Du13HkONp8/DF5gKnjxpkdxTiHDrEtIYGqQf8ZTysMouvWpZ+bG69NnIhNBsw5BcOKutY6CegPrAYOA/O11oeUUu8rpeQGVGGY48ePM3nyZE799pt9adVgx+/8cXv6ab4pX54RV6+Ck3bBn1+2jH+Aao0bmx0l01Curgxq0IAjMTGsXbTI7DgiDRh6n7rWeoXWurjW+mmt9QfJr43UWi+9z7515SxdpIXp06eD1vTw8soQ66SnlbbDh1P5wgVYu9bsKIb4a+1avJWiakiI2VEylXajR5MX+EzuWXcKMqOccCpJSUnM+PZbmgGFRo2CPHnMjpR2WrUiImtWhgwc6HwD5rSmZkQE1zp1omoGWgo3M3B/9ln6+fuz6tAhDh8+bHYc8YSkqAunsmLJEs5cukSv/PlhwACz46QtDw8O1ajB+D//ZM2CBWanSVtHj8KFC7jWrYvFYjE7TeaiFH0GDSIfcGTJErPTiCckRV04lVwbN9IRaP7tt065Fnfo//6HPzDlgw/MjpKmrBs2UAtYcOuW2VEyJf8BAziZLRttdu0yO4p4QlLUhfM4f57q333Hj82a4eaki4G4V6hAt3z5WHbwIOfPnTM7Tpo5tHQpW4CEHDnMjpI5+fri+uqr2BYt4viGDWanEU9AirpwGj+0a8e5mzdh4kSzoxjq+W7dsAJhzrJ6m9Zs37YNgKrVqpkcJhN77TV6urhQq2VLEhMTzU4jHpMUdeEUDs6dS9dNm/ikfHkoUcLsOIaqMHgwxZXiprOMgo+MZPu1a+Ty9eWpp54yO03mlTcvbRo0IOrmTX5x5kmOnJwUdeHwtNXKkFdewU8phjvbALL7cM2Rg6Mvv0y3P/6AmBiz4zy5sDC2A1UrVUI5+CRBjq75Z59RFJgwejRa6xT3FxmPFHXh8JYPGsSqq1cZ2b49OTPLQiA9ekBMDHrePLOTPDHrhg2U8fCg2QsvmB0l03MJCuKNZ55hx+nTbJKFXhyScrRvY8HBwTo8XOaoEXbxp05ROjAQd09P9l25gpsTjni/n4MHDtC8YkWmFCtGk4gIs+M8mcBA+6x/CxeanUQAsevWEdiwIc2qVGHW9u1mxxHJlFK7tdYpTo8pZ+rCod0aMoSqWvPZ559nmoIOUDgwkNM2G1sOHwZHnjDkxAlunTgBdeqYnUQk86pfnw1lyzItOtpppyR2ZlLUhePauJHs8+bxw9ChNO7Z0+w06crPz4/yZcqwRSmYPt3sOI9v0yYaA20XLzY7ibhNKYLeew/XyEis0nvicKSoC8eUkMDYdu3Ynz8/vPOO2WlMUat+fbYrRcLs2ZCQYHacx3Jr3Tp2AsUqVTI7irhbaChhBQpQ+KWX+PvYMbPTiEcgRV04pC39+zM0Opr5tWqBt7fZcUxRs2ZN4mw2/rh4EVatMjvOY9m+di2JQJ26dc2OIu5msVC8f38uJibyyeDBZqcRj0CKunA41r/+4rVp0wjw8mKYI3c9P6GaNWvySp8+ZMmeHX74wew4j+7MGcLOnMFFKWrWrGl2GnGPfAMH0s3Tk5m//sqpU6fMjiNSSYq6cCxa823r1uzVmvHjx+Pj42N2ItPkyZOHbyZNIujFF2HpUrh61exIjyYsjDCgQsmSZMmSxew04l5eXgzr2xdts/HxW2+ZnUakkhR14VCipk5lWEQETUqWpN0rr5gdx3RWq5V9lSuj4+Ph55/NjvNoNm3iVU9P3h450uwk4gEKjRhBd4uFafPnc86J1hpwZlLUheO4cYOco0YxKE8evl26VGYfA77//nvKv/QSEYUKOV4XfFgY7erV44UOHcxOIh4kVy5GdOrEChcX8thsZqcRqSBFXTiOkSPxOneOMUuWUKRYMbPTZAj16tUD4LfSpWHjRjh50txAqXXhAuGHD3OoZEmzk4gUFHzvPRrYbKivvjI7ikgFKerCIVwJC6P6Z5+xpVUrqFLF7DgZRuHChSlWrBi/xcfbX/jxR3MDpdamTQwFOv/6q9lJREqefhratmXohAm8OWCA2WlECqSoi4zPauXtF15gJ+A9aJDZaTKchg0bsnHnThKrVYPvvwcHmPo5Yf16tgJ1mjQxO4pIjSFDuBYfzxfffENkZKTZacRDSFEXGd6mt95ianQ0r4eEUFHuZ/6Phg0bEhMTw47q1SEiAvbuNTtSisLXrCEWqJN8+UBkcJUr826VKrjabIzMpJM9OQop6iJDuxkZSc/PPiPQ05NRTrAimREaNmzIihUrqDhwILi5ZfwBc5cvs+bvv3FRirryJc1h5B8xgoFa8+Pcuezbt8/sOOIBpKiLDG1G+/b8ZbMxY9IkfHx9zY6TIWXJkoVmzZrhXbAgNG9uv66ekRfi2LyZNUDloCBy5MhhdhqRWs2b81axYmR1cWHEiBFmpxEPIEVdZFzLltEvPJywrl2p17Wr2WkytH/++YfRo0dzvW1bOHcO1q0zO9KDhYWxwsODmY4yqE/YubiQfehQZlmtjG/c2Ow04gGkqIsM6dqBA5x+6SVcKlSg9qRJZsfJ8I4fP86oUaMI8/aGrFkzdhd8WBjZqlWjZNmyZicRj6pLF0IrVKDkxx873gyGmYQUdZHxJCQwoH59yl+/zrUZM8DT0+xEGV61atXw8vLit02boF07WLQIbt40O9Z/XbvGJ3v28I2Xl9lJxONwdYWpU7l27hzPVajAnDlzzE4k7mFoUVdKNVVKHVVKHVNKDb3P9jeUUhFKqf1KqXVKqcJG5hGOYU5oKN9dvMirzz1H1vLlzY7jEDw9PalVqxa//fYbvPiivaBnwDXKbZs386nWbI6LMzuKeFyVKuH3xhscP36ct19/nZsZ8ctjJmZYUVdKWYCvgWZAENBRKRV0z257gGCtdVlgITDOqDzCMfw9aRKvrFpFjXz5ePenn8yO41AaNWpEREQEUYGBULhwhuyC379wIeeBpjI1rENzef99PsuXj9PR0Yz/+GOz44i7GHmmXhk4prWO1FonAHOB0Lt30Fpv0FrfSn66HQgwMI/I4BKOHqVj//5YLBbmbNyIq6ur2ZEcSvPmzfH29mb/wYPQuTOsWQPnz5sd619WrV8PQOOWLU1OIp6Itze1Zs/mBWDsxx8TFRVldiKRzMiiXgC4exHeqOTXHqQHsPJ+G5RSvZVS4Uqp8Ojo6DSMKDKMhAQSO3WiuIsLUz//nMLFi5udyOGULFmSixcvEhISYu+Ct9kgI/V2xMSw6tQpyufJQ758+cxOI55Uo0aMe+45bElJvC8zPWYYRhb1+y2hdd/5K5VSLwLBwCf32661nqK1DtZaB/v7+6dhRJFR6CFD8PnjD36YO5fn+/UzO45DUkrhlTwATZcsCZUqZagueNvvv+MJtGrY0OwoIo0ETp7MgixZ+CgyMmPPjZCJGFnUo4CCdz0PAM7cu5NSqiEwAmiltY43MI/IoP766itqfvEFf3XpAm3bmh3HoUVFRVGhQgUWLFhgP1vfvRsOHzY7FgAumzezymJh9Lffmh1FpJWcOWk5aRI59+zB+sUXJCQkmJ0o0zOyqO8Ciimliiil3IEOwNK7d1BKVQAmYy/oFwzMIjKoW4cO8dzAgRyxWHCXOaWfWL58+Th9+jSLFy+Gjh3BYskwZ+ux69fbew/8/MyOItJShw7ENGpElTffZOywYWanyfQMK+pa6ySgP7AaOAzM11ofUkq9r5RqlbzbJ4AvsEAptVcptfQBhxNOSMfF0at2bQ7abPw4fTqFZY30J2axWGjRogUrVqwgMUcOaNQI5syxX183UeL16xTato2PPDxMzSEMoBS+U6dSVCk++Owzjv31l9mJMjVD71PXWq/QWhfXWj+ttf4g+bWRWuulyT831Frn0VqXT360evgRhTP5sF49frx8mf917kwTmQY2zYSGhnLt2jXCwsLsXfAnTsCWLaZmCpsyhYtAUO3apuYQBilcmAkjR+Jus9GnTRtsJn+JzMxkRjlhioSffmLJ9u10LlmSYd9/b3Ycp9KoUSO8vLxYsmQJtG4NPj6md8Evnj8fL6DRq6+amkMYJ/+IEXwSGMj6Q4eYMmGC2XEyLaX1fQekZ1jBwcE6PDzc7BjiSfz9N1SsyK1ixXBZvx7PLFnMTuR0PvzwQ4oXL87zzz8PXbrA0qX2hV5MmHLXarUS4OVFNW9vFsl84U5N799P4/LluZEzJ1vPn8fFRc4b04pSarfWOjil/eRvXKSr47//Ttfy5YlxccH755+loBtk+PDh9oIO9i74a9dg+XJTsmyePJlziYl0aNPGlPZF+lFly/LjwIGEXbyIS0ZeKdCJSVEX6ebi7t00rVuXpTExRE2aZJ/KVBjmzJkz9uvqDRpA3rymdcGXWruWTz08aP7RR6a0L9KX/0cf4VG8ODd69WLDihVmx8l0pKiLdHFz/35aVK/O8aQklk6dSsn27c2O5PT69etHx44dsQJ06mQ/U798OX1DnDpFnmXLeKN/f3zy5k3ftoU5PD1h6lQGnThBSKtWHDxwwOxEmYoUdWG4uP37ea5yZXYlJDB3wgRq9expdqRMoUOHDpw9e5ZNmzbBSy9BYiLMnZuuGXYMH873NhsJffqka7vCZLVr8+HAgWS1WulQty6xt26l/B6RJqSoC2NFRHCmfn32JyYyZfRoWr/+utmJMo2WLVvi6+vLDz/8AOXKQeXKMHYspNeypzExTJw3j9fd3ORSSyaUZ+JEvmvblkOXL/NG9ergYIOyHZUUdWEY65496Dp1eMrdnaM7d9Jj5EizI2Uq3t7etGvXjvnz5xNz8yZ89BGcPAnffJMu7V/59lsWJybSqVUr3N3d06VNkYEoReOFC3mrYkUm7dvHouefl8KeDqSoC0PEb99OuypVeCsuDsLC8KtUyexImVK3bt24efMmW7Zsgfr1oXFj+OAD+2h4I9lszPvkE+KBrm+/bWxbIuNSijFbt9I3KIjKixbB+++bncjpSVEXae7mpk20rFWLRYmJFHz9dZDpX01To0YNTp06RdOmTe0vfPyxfbDcuHHGNrxiBdOjoykdEEBF+UKXqbl7ePDtgQMEdOtGzKhRnHzrLbMjOTUp6iJNXVm9mkb167MuKYlZ48czQL6Zm0opRYECBQDsU3dWqGBf6GXiRDh71rB2r4wbR7ybG33ffBOl7rcKs8hUXFxg6lQGFy1K+U8+YcnLL5udyGlJURdpxrpxIw1CQthts7FgyhS6Dh5sdiQBaK1p06YNAwcOtL8wZox9JPzo0cY0uG8f2TdvZt+YMfTt18+YNoTjsVgYsmwZRbJnp/Xs2bxet64s1WoAKeriySUlwdixWJo2ZWjevKyeP5+2vXqZnUokU0qRJUsWZs2axfXr1+Hpp6FvX5g2Df78M83bi/nkE657eaF698bV1TXNjy8cV9GSJdl68iSvPf00n4WFUaNYMU6cOGF2LKciRV08mf37mVW0KHOGDoXmzWm3dy91b09PKjKM1157jZiYGGbOnGl/4Z137JOEjBiRtg2dP88Xc+cSkJTEufj4tD22cAoevr58ERHBouBgLp48iWXePLMjORUp6uLxJCQQO3w4PcqXp9uJE8ypUAG9cCH4+5udTNxHcHAw1atXZ8KECSQmJkKePPDmm7BwIezalWbtxH35JV9YrVSvUoW8MoOceBB3d9ps2cKfTZoQMHQoeuZMRowYIWftaUFr7VCPSpUqaWGyHTv04aJFdVnQgH5n8GCdlJRkdiqRgl9//VUDevbs2fYXrl/X2t9f63r1tLbZnryB2Fg9xddXA3rdunVPfjzh/GJjtW7USEeA9nV11b4+Pnry5MnaarWanSzDAcJ1KmqknKmL1Lt1C4YM4XTVqlT8+29OZ8nCihUrGDN+PBaLxex0IgUhISF89dVXtGjRwv6Cn5+9G37DBliz5omPn/Ddd3wYE0NwiRLUq1fviY8nMgFPT1i6lFIjR3JAKSrHxdGnTx9qVK/Onj17zE7nkGQ9dZE6mzYR/fLL+P/zD/TuzZSgIFq1by9drI4uPh5KlYKsWWH3bvutR49Da1YUKULzEydYsXw5zUJC0jancH5Hj6JffZXv169niKsrvnny8OeJE3LCkEzWUxdp4+JFrvfqxZt16lDw+HH2TJ4MkyfTe+BAKegOKiwsjJCQEOLi4sDDw36L2969T7bYy7p1hJw4wYH336dps2ZpF1ZkHiVKoH77jS5z5nA0e3YWnjmD5fXXuXX2LNOnT7ePBREpkqIu/ktrCAsjrn17vs6bl+LTpjEBeKlrVwq0bm12OvGEkpKSWLlyJV9++aX9hY4d7Qu+vPMOPOZ9wxfHjYPcuSkzZIhMNiMen1LQqRPZ/vyTCq++Cl99xcJSpejZsyelSpXip59+sk+iJB4sNRfeM9JDBsoZ6OJFrSdM0LpkSZ0EuqSLiwZ0zYoV9c6dO81OJ9JQy5Yttbe3tz5x4oT9hZUrtQatv/jikY91bM0a7Q16duvWaZxSZHo7d2pbxYp6OeiyyYMwy5Urp5ctW6ZtaTG404GQyoFyck09s9MatmyByZM5On8+CxITGVG1KqpvX6bfukVg8eLUr19fzr6czIkTJwgKCqJBgwYsWbIEBdCgARw8CH//bR9Elwo2m42mhQuzLSqKo/v2kb9sWUNzi0zIaoVvv8U2fDhzY2N519OTgMKFCdu9Gzw8sNlsuDzuWBAHktpr6qafeT/qQ87U00h0tLZNnKiPPPWU/hB0cPJZuZurqz548KDZ6UQ6GD9+vAb0smXL7C/s2GE/W3/vvVQf47MPP9SA/rZaNWNCCnHbmTNa9+mjE7Jl01Ggddas+ky7drpI3rx67Acf6IsXL5qd0FDImbq4Q2tshw9zesUK/lq3jqf+/JPAyEg2ArdvPKocHMzz7drRpUsX8uTJY2JYkV5sNhvz58+nXbt2/3+m88ILsGoVHDtmn6DmIXbu3EntGjVonJTEkn37UHKWLtJDQgKsWwfz5nHk5595LSaG3wB3Fxda16xJj6FDadC4sdONmk/tmboUdSditVo5c+YMxMZS8OxZbmzYQLfJk/kzOppjViuxyfu9W6wY73frRkzduny3Zw8tW7akYMGCpmYX5vrrr7/w8PCgUFwcBAVB69bQtKl9Xv/ERPuftx/Jz3/Yu5cxK1awpXp1/DdtMvsjiMwoPh7WrGH/5MlMX7OGHxITuQKc7NSJgHbtuFywINnLlUM5QYHPEEVdKdUU+BywANO01h/fs90D+A6oBFwC2mutjz/smGlR1IP/t5aLMf8d5ZvL153wdxo90bGNdPPmTc6dO4fNaqWYvz+cOEGfoUP5KzKS4xcucOr6dZK0pqdSTNUaGxBkccPHOysqZwAxBUpxJaA0eQsXZd/Y9mZ/nDTjqL/PR2XU50xKSiIoKAilFKtWraLI+PHwzTcP3D8W8AJwcyPe0xOPxYuhfv3Hbv9e8vuUz/lY4uKIX7aMbZMmUXfHDrh5k6bAIaUIzZWL0LJlqVm3Ll7ly0Pp0lC48OPPy3AfRn/O1BZ1w5ZQUkpZgK+BRkAUsEsptVRrHXHXbj2AK1rrokqpDsBYwPBqc7+/+Ie9bqQL585x4Z9/uHL6NFfOnuXS2bN4JCXRqVIluHaNzpMmsevECc7euEFM8n2aTSwWVlmtABxMPk41Fxc6ZstG4Tx5qFitGjz3HBXW3CDOy4+4u9rzAq6l6yc0Xkb6fRrJqM/p6urKzJkzadmyJVWqVOG7mTNp+tZbYLGAq6v94eYGrq4sXrGCfoMGMWP6dJo0bYrHE7V8f/L7lM/5WDw98XjhBeq+8ALExsLevXSZMYMFGzYw459/+HrdOtzXreMV4DMAHx9uFC+OX+nS9ktNOXPaHzly/P/Pt597eaXYfEb5fRq5LmJl4JjWOhJAKTUXCAXuLuqhwKjknxcCXymllE6nawJPXYqi3IbpHL5xEVdccHVRjFk/GncXF15/5hncLRZ2nD9PxOXLdwYh2Gw2sNnoVbQoWK2sO3uWQ1eukGSzEZuURGxSEhatGR0YCAkJjD9+nE1XrxJrtRJjtXIlKYkcSrHVzw/i4+lw8yYb7skVBHRK/jkHUBHI5+lJXn9/8ubNS7Gnn4aqVaFwYX4vVAgKFYLcue33eN7l2ublhv8dCudQo0YNtm/fTmhoKM1atCAkJITFixfj5ubG33//TXh4ONOmTeO3336jfPnyFHnqKbMjC/FgXl5QrRqdqlWjE3Dr1i3Wr1/PpjVrKOnrC0WKcGPPHrJPmkSe/fspozWlbTaCgNpAccCGfXELy+3j5cwJvr72CZs8PMDdHe3uTrybG2e1pveRc5y1JlIqax68PX0ZW+dlYt090/2jG1nUCwCn7noeBVR50D5a6ySl1DUgJ3Dx7p2UUr2B3gCFChVKs4C5bl3F+8Q+fk9KIBGIB5acs28bEBkJLi78eP06X8TG/ut9LkCvCxfs26OjmXHt2r+2+bu6Mjp3bnB355xSRAFenp74ublR0MODQlmyQJ064O7O0PPn6eviQvZcucieOzc58+UjT+HC9iKdNStfZskCPj5p2k0kxP0UL16cvXv3MnHiRHbu3ImbmxsAbdq04cCBAwQEBDBu3DgGDRp0Z5sQjsDb25sWLVr8/7oHQOLly3wYGMihQ4c4dOgQkyIiiI2N5cuhQyneuDERe/ZQfsgQ/NzdcVcK98uX8bhyhU+DggjNnp2N587RcN8+rPe09d3Fk7RUik9qd0nfD5nMyKJ+vxub7z0DT80+aK2nAFPAfk39yaPZ7SxYhp2DF5Hz/9vh7w+akpCQgKeXFyjFqCtXeP3aNZRSuLi44OLiYr9nO39+AD6P+b/27ihErquO4/j3tzubZJPZ2F0SQ+yGbB4WUUQsoX1oRYLWGoq0FixEBAs+tVCsD5Ks1rYYs4hUROKTlgZUWosQhYWmEEsUFYmmlliraUuQll1brBhKDFKWdH8+zCTdjZvsBDL3zN75fWBgzuUu+zuHmfufe++5957j0fl5BgYGWL9+PUNDQ0uu6f7OChluu1adibgG1q5dy9TU1JJlBw8eZOpMz54AAAVcSURBVHh4mJ07d9JodHOTEVGdsbEx9u7de7G9sLDA7OwsIyMjMDbGdZOTTJ05w9mzZ5mfn7/42nTffXDLLUy8+ir7HnuMZrPJli1b+OrRN2hs3MzXN76XhxrlfvR28xs6ByyeUj0OvH6ZdeYkNYD3AGe6mOmKJNFoNJZsuEZHRxkdHb3s3zSbzSqiRRSza9eu0hEium5gYIDt27dfbI+Pj3PgwIHLrj8xMcH09PTF9v5XeuN0ZzeP6Z4AJiXtkLQG2APMXLLODHBP+/1ngWNVnE/f1FxzVctXq/Qz/VyN0s/0czXqlX52+5K222lNNBwEDtmelrSf1p1xZiStA34C3EBrD33PhYl1l5Pr1CMiot8Uv6QNwPYR4Mglyx5e9P5t4O5uZoiIiOgXmVIdERFREynqERERNZGiHhERURMp6hERETWRoh4REVETq+7Rq5L+BbxW4F9v4pLb1/aZfu8/ZAz6vf+QMYCMQan+b7e9eaWVVl1RL0XSc51cI1hX/d5/yBj0e/8hYwAZg17vfw6/R0RE1ESKekRERE2kqHfuh6UDFNbv/YeMQb/3HzIGkDHo6f7nnHpERERNZE89IiKiJlLUIyIiaiJFvUOSvinpBUknJR2V9L7Smaom6VFJL7XH4ReSriudqWqS7pb0V0kLknr2spZrTdJuSS9LOi1pqnSeqkk6JOlNSS+WzlKCpG2SfiXpVPvz/0DpTFWTtE7SHyX9uT0G3yidaTk5p94hSRttn22//xLwQdv3Fo5VKUm3Acdsn5f0bQDb+wrHqpSkDwALwA+Ar9h+rnCkrpM0CLwCfBKYA04An7P9t6LBKiTpY8A54Me2P1Q6T9UkbQW22n5e0gjwJ+AzffYZELDB9jlJQ8DvgAdsHy8cbYnsqXfoQkFv2wD03a8h20dtn283jwPjJfOUYPuU7ZdL56jYTcBp23+3PQ88BdxZOFOlbP8GOFM6Rym237D9fPv9f4BTwPVlU1XLLefazaH2q+fqQIr6VZA0LWkW+DzwcOk8hX0ReKZ0iKjE9cDsovYcfbZBj3dJmgBuAP5QNkn1JA1KOgm8CfzSds+NQYr6IpKelfTiMq87AWw/aHsb8ARwf9m03bHSGLTXeRA4T2scaqeTMegzWmZZz+2hRPdJagKHgS9fcvSyL9h+x/ZHaB2lvElSz52KaZQO0Ets39rhqk8CTwOPdDFOESuNgaR7gE8Dn3BNJ2RcxeegX8wB2xa1x4HXC2WJQtrnkQ8DT9j+eek8Jdl+S9Kvgd1AT02ezJ56hyRNLmreAbxUKkspknYD+4A7bP+3dJ6ozAlgUtIOSWuAPcBM4UxRofYksceBU7a/WzpPCZI2X7jiR9IwcCs9WAcy+71Dkg4D76c18/k14F7b/yibqlqSTgNrgX+3Fx3vwysA7gK+D2wG3gJO2v5U2VTdJ+l24HvAIHDI9nThSJWS9FNgF63Hbv4TeMT240VDVUjSR4HfAn+htQ0E+JrtI+VSVUvSh4Ef0foODAA/s72/bKr/l6IeERFREzn8HhERURMp6hERETWRoh4REVETKeoRERE1kaIeERFREynqERERNZGiHhERURMp6hFxRZJulPRC+3nSG9rPku65e15HRG4+ExEdkHQAWAcMA3O2v1U4UkQsI0U9IlbUvuf7CeBt4Gbb7xSOFBHLyOH3iOjEGNAERmjtsUdED8qeekSsSNIM8BSwA9hq+/7CkSJiGXmeekRckaQvAOdtPylpEPi9pI/bPlY6W0QslT31iIiImsg59YiIiJpIUY+IiKiJFPWIiIiaSFGPiIioiRT1iIiImkhRj4iIqIkU9YiIiJr4Hwu5QnLsP6EYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f124ccdef98>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# https://github.com/hplgit/INF5620/blob/master/src/fem/fe_approx1D.py\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "#=======================================================\n",
    "# Implementation of the Finite Element Method in 1D\n",
    "# for the approximateion of a function f(x)\n",
    "#=======================================================\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def Uniform_nodes(a, b, N):\n",
    "    '''Computes uniformly spaced nodes'''\n",
    "    h = (b-a) / N\n",
    "    return [a + i*h for i in range(N+1)]\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def init_mesh(Omega, N, n, **kwargs):\n",
    "    ''' Function initializes a 2D array which defines N+1 global \n",
    "        1D finite elements, with n+1 nodes inside each element\n",
    "        in the domain Omega = (x0, x1) '''\n",
    "    node_spacing = kwargs.get('node_spacing', 'uniform')\n",
    "    if (node_spacing == 'uniform'):\n",
    "        nodes = Uniform_nodes(Omega[0], Omega[1], n + N*n)\n",
    "    elif (node_spacing == 'chebyshev'):\n",
    "        nodes = Chebyshev_nodes(Omega[0], Omega[1], n + N*n)\n",
    "    # Initialize elements\n",
    "    elements = [[i-j%(N+1) for i in range(j*(n+1), (j+1)*(n+1))] for j in range(N+1)]\n",
    "    return (nodes, elements)\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def Lagrange_polynomial(x, i, points):\n",
    "    ''' Compute lagrangian interpolation polynomial \n",
    "    over domain x at specific points '''\n",
    "    p = 1\n",
    "    for k in range(len(points)):\n",
    "        if k != i:\n",
    "            p *= (x - points[k]) / (points[i] - points[k])\n",
    "    return p\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def phi_r(r, X, d):\n",
    "    ''' Defines a Lagrangian basis function for the  \n",
    "    local node r of an reference element of \n",
    "    degree d '''\n",
    "    nodes = np.linspace(-1., 1., d+1)\n",
    "    return Lagrange_polynomial(X, r, nodes)\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def basis(d=1):\n",
    "    ''' Create all basis functions for a reference \n",
    "    element of degree d\n",
    "    '''\n",
    "    X = np.linspace(-1, 1, d+1)\n",
    "    phi = [phi_r(r, X, d) for r in range(d+1)]\n",
    "    return phi\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def element_matrix(phi, Omega_e):\n",
    "    ''' Creates the element matrix for an element over\n",
    "    its local domain Omega_e=[x0_e, x1_e] and its \n",
    "    basis functions stored in phi \n",
    "    '''\n",
    "    n = len(phi)\n",
    "    A_e = np.zeros((n,n))\n",
    "    h = Omega_e[1] - Omega_e[0]\n",
    "    detJ = 0.5*h\n",
    "    for r in range(n):\n",
    "        for s in range(n):\n",
    "            A_e[r, s] = phi[r][s]*detJ\n",
    "            A_e[s, r] = A_e[r, s]\n",
    "    return A_e\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def element_rhs(f, phi, Omega_e):\n",
    "    n = len(phi)\n",
    "    b_e = np.zeros(n)\n",
    "    h = Omega_e[1] - Omega_e[0]\n",
    "    detJ = 0.5*h\n",
    "    # Compute mapping from f(x) to f(x(X))\n",
    "    x = np.linspace(Omega_e[0], Omega_e[1], n)\n",
    "    f_map = f(x)\n",
    "    for r in range(n):\n",
    "        b_e[r] = f_map[r]*detJ\n",
    "    return b_e\n",
    "\n",
    "\n",
    "#-------------------------------------------------------\n",
    "def assemble(nodes, elements, phi, f):\n",
    "    ''' Function to assemble the global linear system of \n",
    "    equations A*x=b for a given mesh which is defined \n",
    "    through nodes and elements\n",
    "    '''\n",
    "    N_n, N_e = len(nodes), len(elements)\n",
    "    A = np.zeros((N_n, N_n))\n",
    "    b = np.zeros(N_n)\n",
    "    \n",
    "    for e in range(N_e):\n",
    "        Omega_e = (nodes[elements[e][0]], nodes[elements[e][-1]])\n",
    "        A_e = element_matrix(phi, Omega_e)\n",
    "        b_e = element_rhs(f, phi, Omega_e)\n",
    "        \n",
    "        for r in range(len(elements[e])):\n",
    "            for s in range(len(elements[e])):\n",
    "                A[elements[e][r], elements[e][s]] += A_e[r, s]\n",
    "            b[elements[e][r]] += b_e[r]\n",
    "        \n",
    "    return (A, b)\n",
    "\n",
    "\n",
    "\n",
    "def u_glob(nodes, elements, C, resolution):\n",
    "    x_patches = []\n",
    "    u_patches = []\n",
    "    for e in elements:\n",
    "        (x_L, x_R) = (nodes[e[0]], nodes[e[-1]])\n",
    "        d = len(e)-1\n",
    "        X = np.linspace(-1, 1, resolution)\n",
    "        x = 0.5*(x_L + x_R) + 0.5*(x_R - x_L)*X\n",
    "        x_patches.append(x)\n",
    "        \n",
    "        u_element = np.zeros(resolution)\n",
    "        for r in range(len(e)):\n",
    "            i = e[r] # global node indice\n",
    "            u_element += C[i] * phi_r(r, X, d)\n",
    "        u_patches.append(u_element)\n",
    "    x = np.concatenate(x_patches)\n",
    "    u = np.concatenate(u_patches)\n",
    "    return (x, u)\n",
    "    \n",
    "    \n",
    "# Test assembly of linear equation system\n",
    "#-------------------------------------------------------\n",
    "from numpy import tanh, pi, sin, cos, exp\n",
    "d = 2\n",
    "N = 7\n",
    "\n",
    "Omega = (-pi, pi)\n",
    "f = lambda x: -tanh(x) * sin(x-pi) * exp(1. - 0.75* (x-0.25)**2)\n",
    "\n",
    "#Omega = (0, 1)\n",
    "#f = lambda x: x*(1-x)**8\n",
    "\n",
    "\n",
    "elem_resolution = 40\n",
    "x   = np.linspace(Omega[0], Omega[1], elem_resolution)\n",
    "phi = basis(d)\n",
    "\n",
    "# Initialize mesh\n",
    "(nodes, elements) = init_mesh(Omega, N-1, d)\n",
    "# Assemble the global linear equation system \n",
    "(A, b) = assemble(nodes, elements, phi, f)\n",
    "# Solve the global linear equation system\n",
    "C = np.linalg.solve(A, b)\n",
    "# Compute the approximation u(x) to f(x)\n",
    "(x_u, u) = u_glob(nodes, elements, C, 40)\n",
    "\n",
    "# Plot results\n",
    "#-------------------------------------------------------\n",
    "element_edges = [nodes[elements[i][0]] for i in range(0,N)]\n",
    "element_edges.append(nodes[-1])\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_size_inches(8,4)\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "e_plt = ax.plot(element_edges, [0]*len(element_edges), ls='None', marker='s')\n",
    "\n",
    "f_plt = ax.plot(x, f(x), c='r', ls='-', label='$f(x)$')\n",
    "u_plt = ax.plot(x_u, u, c='k', ls='--', label='$u(x)$')\n",
    "ax.set_xlabel('x') \n",
    "ax.set_ylabel('f(x), u(x)')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}