oned_fem_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Functions Required to Solve Simple Elliptic Equation with Dirichlet Conditions Using Finite Element Method\n",
    "In this module, we will solve the equation\n",
    "\n",
    "\\begin{align}\n",
    "-\\frac{d}{dx} & \\left(q(x) \\frac{d u(x)}{dx} \\right)  = f(x) \\qquad\n",
    "\\Longleftrightarrow \\qquad         &&\n",
    "-\\frac{d q(x)}{dx} \\frac{d u(x)}{dx} - q(x)\\frac{d^2 u(x)}{d^2x} = f(x)\\\\\n",
    "& u(a) = u0,  && \\qquad \\qquad u(a) = u0, \\\\\n",
    "& u(b) = u1  && \\qquad \\qquad u(b) = u1\n",
    "\\end{align}\n",
    "\n",
    "Since this article will try to explain the functions and method of finite element solution for this problem, we will occassionally throw in diagrams where necessary so we'll embed images using markdown embedding syntax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step in solving ODE/PDE using finite element method is to formulate the _weak form_ of the problem. For the elliptic equation above, we do this as follows:  \n",
    "For a _test_ function $v \\in H^1_0(\\Omega) = \\{ v \\in L^2(\\Omega) : v'\\in L(\\Omega), v|_{\\partial \\Omega} = 0 \\}$ where $\\Omega = [a,b]$,\n",
    "\\begin{align}\n",
    "\\int_{\\Omega}-\\frac{d}{dx} \\left(q(x) \\frac{d u(x)}{dx} \\right)v(x) dx  &= \\int_{\\Omega} f(x)v(x) dx\\\\\n",
    "-q(x)u'(x)v(x) |_{\\partial \\Omega} + \\int_{\\Omega} q(x)u'(x)v'(x)dx &=  \\int_{\\Omega} f(x)v(x) dx\\\\\n",
    " \\int_{\\Omega} q(x)u'(x)v'(x)dx &=  \\int_{\\Omega} f(x)v(x) dx\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$u$ which solves the given differential equation is infinite dimensional. However, we approximate $u$ with a  $u^h$ where \n",
    "$$\n",
    "u^h(x) = \\sum_{i=1}^n c_i \\phi_i(x)\n",
    "$$\n",
    "and \n",
    "$\n",
    "\\phi_1, \\phi_2, \\cdots, \\phi_n\n",
    "$\n",
    "are finite number of basis functons which could be _linear, quadratic, cubic, etc_ over the compact domain $\\Omega$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plugging $u^h$ into the weak form and also approximating $u$ by $u^h = \\phi_i$ for $i = 1,2,\\cdots,n$, we get"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    " \\int_{\\Omega} q(x)u'(x)v'(x)dx &=  \\int_{\\Omega} f(x)v(x) dx\\\\\n",
    " \\int_{\\Omega} q(x)\\left( \\sum_{j=1}^n c_j \\phi_j'(x) \\right) \\phi_i'(x)dx &=  \\int_{\\Omega} f(x)\\phi_i(x) dx, \\qquad i = 1,2,3,\\cdots,n\\\\\n",
    " \\sum_{j=1}^n \\left(\\int_{\\Omega} q(x)\\phi_j'(x) \\phi_i'(x)dx \\right)c_j  &=  \\int_{\\Omega} f(x)\\phi_i(x) dx, \\qquad i = 1,2,3,\\cdots,n\\\\\n",
    "\\end{align}\n",
    "\n",
    "The integral\n",
    "$\n",
    "A_{ij} = \\int_{\\Omega} q(x)\\phi_j'(x) \\phi_i'(x)dx\n",
    "$\n",
    "is called the bilinear form."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we write out the last equation explicitly for each $i$, we would get an $n\\times n$ system of equations in $c_j$'s. Below, we have $6$ python functions that computes the bilinear form and the RHS over each element:  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "oned_mesh = generates computational nodes on $\\Omega$ and their corresponding indices  \n",
    "oned_gauss = stores Gaussian quadrature nodes (on $[-1,1]$, the reference element) and weights for approximating the bilinear form  \n",
    "oned_shape = evaluates $q(x)$ (called kernel in the code), $\\phi_j$ (and it's 1st and 2nd derivatives p_x, p_xx) on the Gaussian nodes and on each of the the physical elements   \n",
    "oned_linear = computes the right hand side integral for each $i = 1, 2, \\cdots, n$.  \n",
    "oned_bilinear = computes the bilinear form $\\int_{\\Omega} q(x)\\phi_j'(x) \\phi_i'(x)dx$ for each $i,j$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### oned_mesh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function oned_mesh generates finite elements nodes and the connectivity matrix e_conn which is an $n$ by $m$ indices matrix where $n$ is the number of elements and $m$ is the number of nodes per elements. For linear bases functions $m = 2$, for quadratic $m = 3$, for cubic $m = 4$ and it continues in that manner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "\"\"\"\n",
    "Module for solving one dimensional finite element problems\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "def oned_mesh(x_min, x_max, n_elements, element_type):\n",
    "    \"\"\"\n",
    "    Description: Construct a one dimensional finite element mesh. \n",
    "\n",
    "    Usage: x, e_conn = oned_fem_mesh(x_min, x_max, n_elements, element_type)\n",
    "\n",
    "    Inputs:\n",
    "\n",
    "        x_min: double, left interval endpoint\n",
    "\n",
    "        x_max: double, right interval endpoint\n",
    "\n",
    "        n_elements: int, number of elements\n",
    "\n",
    "        element_type: str, polynomial order 'linear', 'quadratic', or 'cubic'\n",
    "\n",
    "\n",
    "    Outputs: \n",
    "\n",
    "        x: double, vector of finite element nodes\n",
    "\n",
    "        e_conn: int, (n_elements, n_dof) connectivity matrix whose ith row\n",
    "           contains the indices of the nodes x(i) contained in the ith element\n",
    "\n",
    "    Last modified: \n",
    "\n",
    "        02/05/2021, Hans-Werner van Wyk\n",
    "\n",
    "    \"\"\"\n",
    "    if element_type == 'linear':\n",
    "        # Define mesh nodes\n",
    "        x = np.linspace(x_min,x_max,n_elements+1);\n",
    "        \n",
    "        # Define connectivity matrix\n",
    "        e_conn = (np.arange(n_elements),np.arange(1,n_elements+1))\n",
    "        e_conn = np.stack(e_conn,axis=1)\n",
    "        \n",
    "    elif element_type == 'quadratic':\n",
    "        # Define the mesh nodes\n",
    "        x = np.linspace(x_min, x_max, 2*n_elements + 1)\n",
    "        \n",
    "        # Define the connectivity matrix\n",
    "        e_conn = (np.arange(0,2*n_elements-1,2),\n",
    "                  np.arange(1,2*n_elements,2),\n",
    "                  np.arange(2,2*n_elements+1,2))\n",
    "        e_conn = np.stack(e_conn,axis=1)\n",
    "        \n",
    "    elif element_type == 'cubic':\n",
    "        # Define the mesh nodes\n",
    "        x = np.linspace(x_min,x_max,3*n_elements+1)\n",
    "        \n",
    "        # Define the connectivity matrix\n",
    "        e_conn = (np.arange(0,3*n_elements-2,3),\n",
    "                  np.arange(1,3*n_elements-1,3),\n",
    "                  np.arange(2,3*n_elements,3),\n",
    "                  np.arange(3,3*n_elements+1,3))\n",
    "        e_conn = np.stack(e_conn,axis=1)\n",
    "     \n",
    "    else: \n",
    "        # Error\n",
    "        raise Exception('Use only \"linear\", \"quadratic\", '+\\\n",
    "                        'or \"cubic\" for element_type')\n",
    "    \n",
    "    return x,e_conn\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose $a = -1$ and $b = 1$ and number of elements is $10$ (i.e., there are $11$ nodes). Then"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_min = -1\n",
    "x_max = 1\n",
    "n_elements = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# With linear bases elements\n",
    "#\n",
    "element_type = 'linear'\n",
    "x, e_conn = oned_mesh(x_min, x_max, n_elements, element_type)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1. , -0.6, -0.2,  0.2,  0.6,  1. ])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1],\n",
       "       [1, 2],\n",
       "       [2, 3],\n",
       "       [3, 4],\n",
       "       [4, 5]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e_conn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These outputs can be visualized as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Linear Nodes and Elements](./nodes_and_elements1.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# With quadratic bases elements\n",
    "#\n",
    "element_type = 'quadratic'\n",
    "x, e_conn = oned_mesh(x_min, x_max, n_elements, element_type)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1. , -0.8, -0.6, -0.4, -0.2,  0. ,  0.2,  0.4,  0.6,  0.8,  1. ])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2],\n",
       "       [ 2,  3,  4],\n",
       "       [ 4,  5,  6],\n",
       "       [ 6,  7,  8],\n",
       "       [ 8,  9, 10]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e_conn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that we still have $10$ elements but the number of nodes have changed to accommodate for the number of points needed to compute each basis function over an element. Two points are needed to compute a linear function, three points are need to compute quadratic function, etc.\n",
    "Again, we visualize as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Quadratic Nodes and Elements](./nodes_and_elements2.png)\n",
    "\n",
    "Cubic bases functions follow similarly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# With cubic bases elements\n",
    "#\n",
    "element_type = 'cubic'\n",
    "x, e_conn = oned_mesh(x_min, x_max, n_elements, element_type)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.        , -0.86666667, -0.73333333, -0.6       , -0.46666667,\n",
       "       -0.33333333, -0.2       , -0.06666667,  0.06666667,  0.2       ,\n",
       "        0.33333333,  0.46666667,  0.6       ,  0.73333333,  0.86666667,\n",
       "        1.        ])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2,  3],\n",
       "       [ 3,  4,  5,  6],\n",
       "       [ 6,  7,  8,  9],\n",
       "       [ 9, 10, 11, 12],\n",
       "       [12, 13, 14, 15]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e_conn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that in each case ('linear', 'quadratic', 'cubic'),  \n",
    "total # of nodes = (n_dof * n_elements) + 1  \n",
    "where n_dof is the degree of freedom (i.e., # of points required to compute a basis function over each element: n_dof = 2 for linear, 3 for quadratic and 4 for cubic bases) and n_elements is the number of elements.  \n",
    "We could draw a numbered line representation for cubic elements just like we did for linear and quadratic, but hopefully, we understand the concept clearly from those of linear and quadratic as shown previously."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### oned_gauss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "oned_gauss function takes 'rule' as input and returns a vector of Gaussian points (r) in $[-1,1]$ and a vector of Gaussian weights (w) in $[-1,1]$. So\n",
    "$$\n",
    "\\int_{-1}^1 f(x)dx \\approx w_0 f(r_0) + w_1 f(r_1) + w_2 f(r_2) + \\cdots + w_{n-1} f(r_{n-1})\n",
    "$$\n",
    "assuming we have have chosen rule $=n$. In this case, the integral would be exact (accurate) up to degree $2n-1$ polynomial approximation of $f$.  \n",
    "Gaussian points and weights are pre-computed and oned_gauss simply records them for call and use as can be clearly seen from the function definition.\n",
    "\n",
    "A complete and illustrative example on mathematics of the Gaussian Quadrature Rule and how the Gaussian points and weights are derived can be found here\n",
    "[Gaussian Quadrature Rule](http://webhome.auburn.edu/~cae0027/resources/notes/gauss-quadrature-rule.pdf)\n",
    ".\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reset element type to quadratic for illustration purposes\n",
    "#\n",
    "element_type = 'quadratic'\n",
    "x, e_conn = oned_mesh(x_min, x_max, n_elements, element_type)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def oned_gauss(rule):\n",
    "    \"\"\"\n",
    "    Calculate Gauss integration points on (-1,1)\n",
    "\n",
    "    \n",
    "    Inputs: \n",
    "    \n",
    "        rule: int, number of Gauss points (between 1 and 11). The precision of\n",
    "            the Gauss rule is 2*rule - 1\n",
    "            \n",
    "            Int[-1,1] f(x) dx ~= sum_{i=1}^rule wi*f(ri)\n",
    "        \n",
    "    Outputs:\n",
    "    \n",
    "        r: double, (rule,) vector of Gauss points in (-1,1]\n",
    "        \n",
    "        w: double, (rule,) vector of Gauss weights\n",
    "    \"\"\"\n",
    "\n",
    "    r = np.zeros(rule)\n",
    "    w = np.zeros(rule)\n",
    "\n",
    "    if rule == 1:      # up to order 1 polynomials exact\n",
    "        r[0] = 0\n",
    "        w[0] = 2\n",
    "        \n",
    "    elif rule == 2:   # up to order 3 polynomials exact\n",
    "        r[0] =-1.0 / np.sqrt(3.0);\n",
    "        r[1] =-r[0];\n",
    "        w[0] = 1.0;\n",
    "        w[1] = 1.0;\n",
    "        \n",
    "    elif rule == 3:  # up to order 5 polynomials exact\n",
    "        r[0] =-np.sqrt(3.0/5.0);\n",
    "        r[1] = 0.0;\n",
    "        r[2] =-r[0];\n",
    "        w[0] = 5.0 / 9.0;\n",
    "        w[1] = 8.0 / 9.0;\n",
    "        w[2] = w[0];\n",
    "        \n",
    "    elif rule == 4:  # up to order 7 polynomials exact\n",
    "        r[0] =-np.sqrt((3.0+2.0*np.sqrt(6.0/5.0))/7.0);\n",
    "        r[1] =-np.sqrt((3.0-2.0*np.sqrt(6.0/5.0))/7.0);\n",
    "        r[2] =-r[1];\n",
    "        r[3] =-r[0];\n",
    "        w[0] = 0.5 - 1.0 / ( 6.0 * np.sqrt(6.0/5.0) );\n",
    "        w[1] = 0.5 + 1.0 / ( 6.0 * np.sqrt(6.0/5.0) );\n",
    "        w[2] = w[1];\n",
    "        w[3] = w[0];\n",
    "        \n",
    "    elif rule == 5:  # up to order 9 polynomials exact\n",
    "        r[0] =-np.sqrt(5.0+4.0*np.sqrt(5.0/14.0)) / 3.0;\n",
    "        r[1] =-np.sqrt(5.0-4.0*np.sqrt(5.0/14.0)) / 3.0;\n",
    "        r[2] = 0.0;\n",
    "        r[3] =-r[1];\n",
    "        r[4] =-r[0];\n",
    "        w[0] = 161.0/450.0-13.0/(180.*np.sqrt(5.0/14.0));\n",
    "        w[1] = 161.0/450.0+13.0/(180.*np.sqrt(5.0/14.0));\n",
    "        w[2] = 128.0/225.0;\n",
    "        w[3] = w[1];\n",
    "        w[4] = w[0];\n",
    "        \n",
    "    elif rule == 6:\n",
    "        r[0] = -0.2386191860831969;\n",
    "        r[1] = -0.6612093864662645;\n",
    "        r[2] = -0.9324695142031521;\n",
    "        r[3] = - r[0];\n",
    "        r[4] = - r[1];\n",
    "        r[5] = - r[2];\n",
    "        w[0] = 0.4679139345726910;\n",
    "        w[1] = 0.3607615730481386;\n",
    "        w[2] = 0.1713244923791704;\n",
    "        w[3] = w[0];\n",
    "        w[4] = w[1];\n",
    "        w[5] = w[2];\n",
    "        \n",
    "    elif rule == 7:\n",
    "        r[0] = -0.9491079123427585;\n",
    "        r[1] = -0.7415311855993945;\n",
    "        r[2] = -0.4058451513773972;\n",
    "        r[3] =  0.0000000000000000;\n",
    "        r[4] = - r[2];\n",
    "        r[5] = - r[1];\n",
    "        r[6] = - r[0];\n",
    "        w[0] = 0.1294849661688697;\n",
    "        w[1] = 0.2797053914892766;\n",
    "        w[2] = 0.3818300505051189;\n",
    "        w[3] = 0.4179591836734694;\n",
    "        w[4] = w[2];\n",
    "        w[5] = w[1];\n",
    "        w[6] = w[0];\n",
    "        \n",
    "    elif rule == 8:\n",
    "        r[0] = -0.9602898564975363;\n",
    "        r[1] = -0.7966664774136267;\n",
    "        r[2] = -0.5255324099163290;\n",
    "        r[3] = -0.1834346424956498;\n",
    "        r[4] = - r[3];\n",
    "        r[5] = - r[2];\n",
    "        r[6] = - r[1];\n",
    "        r[7] = - r[0];\n",
    "        w[0] = 0.1012285362903763;\n",
    "        w[1] = 0.2223810344533745;\n",
    "        w[2] = 0.3137066458778873;\n",
    "        w[3] = 0.3626837833783620;\n",
    "        w[4] = w[3];\n",
    "        w[5] = w[2];\n",
    "        w[6] = w[1];\n",
    "        w[7] = w[0];\n",
    "\n",
    "    elif rule == 9:\n",
    "        r[0] = -0.9681602395076261;\n",
    "        r[1] = -0.8360311073266358;\n",
    "        r[2] = -0.6133714327005904;\n",
    "        r[3] = -0.3242534234038089;\n",
    "        r[4] =  0.0000000000000000;\n",
    "        r[5] = - r[3];\n",
    "        r[6] = - r[2];\n",
    "        r[7] = - r[1];\n",
    "        r[8] = - r[0];\n",
    "        w[0] = 0.0812743883615744;\n",
    "        w[1] = 0.1806481606948574;\n",
    "        w[2] = 0.2606106964029354;\n",
    "        w[3] = 0.3123470770400029;\n",
    "        w[4] = 0.3302393550012598;\n",
    "        w[5] = w[3];\n",
    "        w[6] = w[2];\n",
    "        w[7] = w[1];\n",
    "        w[8] = w[0];\n",
    "      \n",
    "    elif rule == 10:\n",
    "        r[0] = -0.9739065285171717;\n",
    "        r[1] = -0.8650633666889845;\n",
    "        r[2] = -0.6794095682990244;\n",
    "        r[3] = -0.4333953941292472;\n",
    "        r[4] = -0.1488743389816312;\n",
    "        r[5] = - r[4];\n",
    "        r[6] = - r[3];\n",
    "        r[7] = - r[2];\n",
    "        r[8] = - r[1];\n",
    "        r[9] = - r[0];\n",
    "        w[0] = 0.0666713443086881;\n",
    "        w[1] = 0.1494513491505806;\n",
    "        w[2] = 0.2190863625159820;\n",
    "        w[3] = 0.2692667193099963;\n",
    "        w[4] = 0.2955242247147529;\n",
    "        w[5] = w[4];\n",
    "        w[6] = w[3];\n",
    "        w[7] = w[2];\n",
    "        w[8] = w[1];\n",
    "        w[9] = w[0];\n",
    "\n",
    "    elif rule == 11:\n",
    "        r[0] = -0.9782286581460570;\n",
    "        r[1] = -0.8870625997680953;\n",
    "        r[2] = -0.7301520055740494;\n",
    "        r[3] = -0.5190961292068118;\n",
    "        r[4] = -0.2695431559523450;\n",
    "        r[5] =  0.0000000000000000;\n",
    "        r[6] = - r[4];\n",
    "        r[7] = - r[3];\n",
    "        r[8] = - r[2];\n",
    "        r[9] = - r[1];\n",
    "        r[10] = - r[0];\n",
    "        w[0] = 0.0556685671161737;\n",
    "        w[1] = 0.1255803694649046;\n",
    "        w[2] = 0.1862902109277343;\n",
    "        w[3] = 0.2331937645919905;\n",
    "        w[4] = 0.2628045445102467;\n",
    "        w[5] = 0.2729250867779006;\n",
    "        w[6] = w[4];\n",
    "        w[7] = w[3];\n",
    "        w[8] = w[2];\n",
    "        w[9] = w[1];\n",
    "        w[10] = w[0];\n",
    "\n",
    "    else:\n",
    "        raise Exception('Quadrature rule not supported')\n",
    "        #Return computed quantities\n",
    "        \n",
    "    return r,w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "r,w = oned_gauss(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.90617985, -0.53846931,  0.        ,  0.53846931,  0.90617985])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$r$ is computed on $[-1,1]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.23692689, 0.47862867, 0.56888889, 0.47862867, 0.23692689])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$w$ is computed on $[-1,1]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  oned_shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To apply the Gaussian quadrature rule to the integral \n",
    "$\n",
    "\\int_a^b f(x)dx\n",
    "$,\n",
    "we make the following transformation\n",
    "$$ \n",
    "x = \\frac{b-a}{2}r + \\frac{b+a}{2} \n",
    "$$\n",
    "which is gotten by finding the parametric line between $(-1,a)$ and $(1,b)$. Thus,\n",
    "\n",
    "\\begin{align}\n",
    "\\int_{a}^b f(x)dx &\\approx \\frac{b-a}{2} \\int_{-1}^1 f \\left(\\frac{b-a}{2}r + \\frac{b+a}{2} \\right)dr\\\\\n",
    "                  &= \\frac{b-a}{2} \\left(w_0 f(x_0) + w_1 f(x_1) + w_2 f(x_2) + \\cdots + w_{n-1} f(x_{n-1}) \\right)\n",
    "\\end{align}\n",
    "where \n",
    "$\n",
    "x_i = \\frac{b-a}{2}r_i + \\frac{b+a}{2} \n",
    "$\n",
    "for $i = 0,1,2,\\cdots,n-1$. Again, we have assumed rule $= n$ and for this, the integral is exact up to a polynomial of degree $2n-1$ approximation of $f$. i.e.,\n",
    "$$\n",
    "f(x) \\approx a_0 + a_1x + a_2x^2 + \\cdots + a_{2n-1}x^{2n-1}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\frac{b-a}{2} =$ Jacobian of the transformation from $x$ to $r$  \n",
    "$w = (w_0, w_1,\\cdots, w_{2n-1}) =$ Gauss weights on the reference element $[-1,1]$   \n",
    "$w = (\\frac{b-a}{2}w_0, \\frac{b-a}{2}w_1,\\cdots, \\frac{b-a}{2}w_{2n-1}) =$ Gauss weights on the physical element. In oned_shape, $w\\_g := w$.     \n",
    "$r = (r_0, r_1,\\cdots, r_{2n-1}) =$ Gauss points/nodes on the reference element $[-1,1]$  \n",
    "$x = (x_0, x_1,\\cdots, x_{2n-1}) =$ Gauss points/nodes on the physical element. In oned_shape, $x\\_g := x$\n",
    "\n",
    "\n",
    "In the function oned_shape, $x$ is a vector of nodal values indexed by each row of"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2],\n",
       "       [ 2,  3,  4],\n",
       "       [ 4,  5,  6],\n",
       "       [ 6,  7,  8],\n",
       "       [ 8,  9, 10]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e_conn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "assuming element_type is 'cubic'. So the call of oned_shape must loop over the total number of elements to get each $x$ which now represents .\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def oned_shape(x,r,w):\n",
    "    \"\"\"\n",
    "    computes test functions and derivatives for a Lagrange C0 element given\n",
    "    element coordinates and Gauss points. (assumes all nodes are uniformly \n",
    "    distributed in the element).\n",
    "\n",
    "\n",
    "    Usage:    x_g,w_g,phi,p_x,p_xx = oned_shape(x,r,w)\n",
    "\n",
    "    Inputs:     \n",
    "    \n",
    "        x: double, coordinates of the element nodes\n",
    "                 \n",
    "        r: double, coordinates of Gauss points in (-1,1)\n",
    "        \n",
    "        w: double, Gauss weights associated with r\n",
    "        \n",
    "        \n",
    "    Outputs:\n",
    "\n",
    "        x_g: double, coordinates of Gauss points in the element\n",
    "                 \n",
    "        w_g: double, Gauss weights scaled by the element Jacobian\n",
    "                 \n",
    "        phi: double, value of element shape functions at x_g\n",
    "        \n",
    "        p_x: double, first spatial derivatives of phi\n",
    "        \n",
    "        p_xx: double, second spatial derivatives of phi\n",
    "    \n",
    "    \n",
    "    Modified: \n",
    "    \n",
    "        02/05/2021, Hans-Werner van Wyk \n",
    "    \"\"\"\n",
    "    n_dof = len(x)  # number of dofs\n",
    "    n_g = len(r)  # number of Gaussian quadrature points\n",
    "    \n",
    "    if n_dof==2:\n",
    "        # Transform coordinates for linear elements\n",
    "        c0 = ( x[-1]-x[0] )/2\n",
    "        c1 = ( x[-1]+x[0] )/2\n",
    "    \n",
    "        # Gaussian nodes\n",
    "        x_g = c0*r + c1\n",
    "\n",
    "        # Evaluate basis function at Gauss points\n",
    "        phi = np.empty((n_g,n_dof))\n",
    "        phi[:,1] = ( 1+r )/2\n",
    "        phi[:,0] = ( 1-r )/2\n",
    "        \n",
    "        # Evaluate first derivative at Gauss points\n",
    "        p_x = np.empty((n_g,n_dof))\n",
    "        p_x[:,1] = 0.5*np.ones(n_g)/c0\n",
    "        p_x[:,0] =-p_x[:,1]\n",
    "\n",
    "        # Jacobian\n",
    "        djac = c0\n",
    "    \n",
    "        # Compute physical weights\n",
    "        w_g = djac*w\n",
    "\n",
    "        # Second derivative \n",
    "        p_xx = np.zeros((n_g,n_dof))\n",
    "    \n",
    "    elif n_dof==3:\n",
    "        \n",
    "        # Transform coordinates for quadratic elements\n",
    "        c0 = ( x[-1]-x[0] )/2\n",
    "        c1 = ( x[-1]+x[0] )/2\n",
    "        \n",
    "        # Physical Gauss nodes\n",
    "        x_g = c0*r + c1\n",
    "\n",
    "        # defined backwards to help Matlab create the proper sized array\n",
    "        phi = np.empty((n_g,n_dof))\n",
    "        phi[:,2] = .5*r*( r+1 )\n",
    "        phi[:,1] =-( r+1 )*( r-1 )\n",
    "        phi[:,0] = .5*r*( r-1 )\n",
    "\n",
    "        # First derivative \n",
    "        p_x = np.empty((n_g,n_dof))\n",
    "        p_x[:,2] = ( r+.5 )/c0\n",
    "        p_x[:,1] =-2*r/c0\n",
    "        p_x[:,0] = ( r-.5 )/c0\n",
    "\n",
    "        # Jacobian \n",
    "        djac = c0\n",
    "\n",
    "        # Gaussian weight on physical element\n",
    "        w_g = djac*w\n",
    "        \n",
    "        # Second derivative\n",
    "        p_xx = np.empty((n_g,n_dof))\n",
    "        p_xx[:,2] = np.ones(n_g)/c0**2\n",
    "        p_xx[:,1] =-2*p_xx[:,2]\n",
    "        p_xx[:,0] = p_xx[:,2]\n",
    "\n",
    "    elif n_dof==4:\n",
    "        # Transform coordinates for (nonconforming) cubic elements\n",
    "        c0 = ( x[-1]-x[0] )/2\n",
    "        c1 = ( x[-1]+x[0] )/2\n",
    "\n",
    "        x_g = c0*r + c1\n",
    "\n",
    "        r2  = r*r\n",
    "        r3  = r*r2\n",
    "\n",
    "        # defined backwards to help Matlab create the proper sized array\n",
    "        phi = np.zeros((n_g,n_dof))\n",
    "        phi[:,3] =  9*( r3+r2-r/9-1/9 )/16\n",
    "        phi[:,2] =-27*( r3+r2/3-r-1/3 )/16\n",
    "        phi[:,1] = 27*( r3-r2/3-r+1/3 )/16\n",
    "        phi[:,0] =- 9*( r3-r2-r/9+1/9 )/16\n",
    "\n",
    "        # Derivative on reference element\n",
    "        p_r = np.zeros((n_g,n_dof))\n",
    "        p_r[:,3] =  9*( 3*r2+2*r-1/9 )/16\n",
    "        p_r[:,2] =-27*( 3*r2+2*r/3-1 )/16\n",
    "        p_r[:,1] = 27*( 3*r2-2*r/3-1 )/16\n",
    "        p_r[:,0] =- 9*( 3*r2-2*r-1/9 )/16\n",
    "        \n",
    "        # Second derivative on reference element\n",
    "        p_rr = np.zeros((n_g,n_dof))\n",
    "        p_rr[:,3] =  9*( 6*r+2   )/16\n",
    "        p_rr[:,2] =-27*( 6*r+2/3 )/16\n",
    "        p_rr[:,1] = 27*( 6*r-2/3 )/16\n",
    "        p_rr[:,0] =- 9*( 6*r-2   )/16\n",
    "\n",
    "        # Jacobian \n",
    "        dxdr = p_r.dot(x)\n",
    "        drdx = 1/dxdr\n",
    "        \n",
    "        # Quadrature weights on physical element\n",
    "        w_g = dxdr*w\n",
    "        \n",
    "        # Derivative on physical element\n",
    "        p_x = np.zeros((n_g,n_dof))\n",
    "        p_x[:,3] = p_r[:,3]*drdx\n",
    "        p_x[:,2] = p_r[:,2]*drdx\n",
    "        p_x[:,1] = p_r[:,1]*drdx\n",
    "        p_x[:,0] = p_r[:,0]*drdx\n",
    "        \n",
    "        # Second derivative on physical element\n",
    "        p_xx = np.zeros((n_g,n_dof))     \n",
    "        p_xx[:,3] = p_rr[:,3]*drdx**2\n",
    "        p_xx[:,2] = p_rr[:,2]*drdx**2\n",
    "        p_xx[:,1] = p_rr[:,1]*drdx**2\n",
    "        p_xx[:,0] = p_rr[:,0]*drdx**2\n",
    "    else:\n",
    "        raise Exception('Elements higher than cubic not currently supported')\n",
    "    \n",
    "    # Return computed quantities\n",
    "    return x_g, w_g, phi, p_x, p_xx    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since oned_shape is computed over each element, we print the function returns over the last element for the sake of illutsration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# Recompute mesh for this illustration\n",
    "#\n",
    "x, e_conn = oned_mesh(x_min, x_max, n_elements, element_type)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# Ilustration of how the function is used\n",
    "for i in range(n_elements):\n",
    "    y = x[e_conn[i]]\n",
    "    x_g, w_g, phi, p_x, p_xx  = oned_shape(y,r,w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.61876403, 0.69230614, 0.8       , 0.90769386, 0.98123597])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.04738538, 0.09572573, 0.11377778, 0.09572573, 0.04738538])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w_g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.86367088,  0.17883809, -0.04250897],\n",
       "       [ 0.41420925,  0.7100508 , -0.12426006],\n",
       "       [-0.        ,  1.        ,  0.        ],\n",
       "       [-0.12426006,  0.7100508 ,  0.41420925],\n",
       "       [-0.04250897,  0.17883809,  0.86367088]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first column of phi represents the the blue curve (below) evaluated at the Gaussian nodes; the second column, the orange curve; and the third column, the green curve."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Quadratic Nodes and Elements](./quadratic_elements1.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-7.03089923,  9.06179846, -2.03089923],\n",
       "       [-5.19234655,  5.3846931 , -0.19234655],\n",
       "       [-2.5       , -0.        ,  2.5       ],\n",
       "       [ 0.19234655, -5.3846931 ,  5.19234655],\n",
       "       [ 2.03089923, -9.06179846,  7.03089923]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first column of p_x is the derivative of the blue quadratic curve above evaluated at the Gaussian nodes; the second column, the orange curve; the third column, the green curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 25., -50.,  25.],\n",
       "       [ 25., -50.,  25.],\n",
       "       [ 25., -50.,  25.],\n",
       "       [ 25., -50.,  25.],\n",
       "       [ 25., -50.,  25.]])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_xx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first column of p_xx is the derivative of the blue quadratic curve above evaluated at the Gaussian nodes; the second column, the orange curve; the third column, the green curve.  \n",
    "Clearly, we can see p_xx is constant for each of the curves since they are quadratics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### oned_linear"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "oned_linear computes approximation of the integral \n",
    "$\n",
    "\\int_{\\Omega} f(x)\\phi_i(x) dx\n",
    "$\n",
    "for each $i$, using Gaussian quadrature rule.  \n",
    "kernel = f(x_g) where x_g is a vector of Gaussian nodes on the physical element  \n",
    "w_g = Gaussian weights on the physical element  \n",
    "test = matrix of size # of gaussian nodes by n_dof representing basis function an element"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "    \n",
    "def oned_linear(kernel, test, w_g):\n",
    "    \"\"\"\n",
    "    Compute the one-dimensional local linear form (f,phi)\n",
    "    \n",
    "    Inputs:\n",
    "    \n",
    "        kernel: double (n_gauss,) kernel evaluated at the Gauss nodes.\n",
    "        \n",
    "        test: double, (n_gauss, n_dofs) array of element test functions\n",
    "            evaluated at the Gauss nodes.\n",
    "            \n",
    "        w_g: double, (n_gauss,) Gauss quadrature weights on element.     \n",
    "        \n",
    "        \n",
    "    Output:\n",
    "    \n",
    "        L: double, (n_dofs,) vector with entries \n",
    "        \n",
    "            L_i = int_Ke f(x) phi_i(x) dx ~= sum_j w_g(j) f(x_j) phi_i(x_j) \n",
    "    \"\"\"\n",
    "    b = np.dot(test.T, kernel*w_g)\n",
    "    return b\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since oned_linear loops over the total number of elements, we call the function over the last element for illustration. We will take $f(x) = x^2$ for simplicity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = lambda x: x**2\n",
    "f_g = f(x_g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.02293333, 0.1728    , 0.0656    ])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "oned_linear(f_g,phi,w_g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These entries correspond to the integral \n",
    "$\n",
    "\\int_{\\Omega} f(x)\\phi_i(x) dx\n",
    "$\n",
    "when $\\phi_i$ is blue, orange, and green respectively as shown."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Quadratic Nodes and Elements](./quadratic_elements1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### oned_bilinear"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "oned_bilinear computes approximation of the integral \n",
    "$\n",
    "\\int_{\\Omega} q(x)\\phi_j'(x) \\phi_i'(x) dx\n",
    "$\n",
    "for each $i,j$, using Gaussian quadrature rule. Just like the previous functions, oned_bilinear computation is done over an element and then a loop is performed over all possible elements. In the code:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "kernel = a vector of $q(x)$ evaluated at the Gauss nodes    \n",
    "phi = matrix whose columns represent $\\phi_j'(x)$ evaluated at the Gauss nodes. For linear shape/hat/basis functions, the matrix has 2 columns; for quadratic, 3 columns; etc.  \n",
    "test = matrix whose columns represent $\\phi_i'(x)$ evaluated at the Gauss nodes. For linear shape/hat/basis functions, the matrix has 2 columns; for quadratic, 3 columns; etc.  \n",
    "w_g = Gauss weights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the sake of illustration, we take $q(x) = x$ and retain w_g, phi above. Furthermore, we take test = phi (Galerkin method)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = lambda x: x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.61876403, 0.69230614, 0.8       , 0.90769386, 0.98123597])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q(x_g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def oned_bilinear(kernel,phi,test,w_g):\n",
    "    \"\"\"\n",
    "    Compute the local bilinear form (kernel, phi, test)     \n",
    "    \n",
    "    \n",
    "    Inputs:\n",
    "    \n",
    "        kernel: double (n_gauss,), kernel function in the integral evaluated at the\n",
    "            Gauss points.\n",
    "    \n",
    "    \n",
    "        phi: double, (n_gauss,n_dofs) matrix of trial functions evaluated at \n",
    "            Gauss points. \n",
    "            \n",
    "            \n",
    "        test: double, (n_gauss,n_dofs) matrix of test functions evaluated at \n",
    "            Gauss points. \n",
    "            \n",
    "        w_g: double, (n_gauss,) Gauss quadrature weights.\n",
    "        \n",
    "        \n",
    "    Output:\n",
    "    \n",
    "        B: double, (n_dofs, n_dofs) local blinear form over element Ke with entries\n",
    "        \n",
    "            B_ij = int_Ke f(x) phi_j(x) test_i(x) dx \n",
    "                 \n",
    "                 ~= sum_{k} w_g(k) f(x_k) phi_j(x_k) test_i(x_k) \n",
    "        \n",
    "    \"\"\"\n",
    "    M = np.dot(np.dot(test.T,np.diag(kernel*w_g)),phi)\n",
    "    return M\n",
    "    \n",
    "    \n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_loc = oned_bilinear(q(x_g),p_x,p_x,w_g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 4.        , -4.66666667,  0.66666667],\n",
       "       [-4.66666667, 10.66666667, -6.        ],\n",
       "       [ 0.66666667, -6.        ,  5.33333333]])"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_loc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output of oned_bilinear, A_loc is the local assembly of the stiffness matrix, A over the specific element."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Quadratic Nodes and Elements](./quadratic_elements1.png)\n",
    "\n",
    "Since we ran the code on the last element which is indexed by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 8,  9, 10])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e_conn[-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "it follows that index $8$ represents the blue parabola; index $9$, the orange parabola; and index $10$, the green parabola. Thus \n",
    "\n",
    "A_loc$_{8,8} = \\int_{\\Omega} q(x)\\phi_8'(x) \\phi_8'(x) dx = 4$  \n",
    "A_loc$_{8,9} = \\int_{\\Omega} q(x)\\phi_8'(x) \\phi_9'(x) dx = -4.66666667$  \n",
    "A_loc$_{8,10} = \\int_{\\Omega} q(x)\\phi_8'(x) \\phi_{10}'(x) dx = 0.66666667$  \n",
    "A_loc$_{9,8} = \\int_{\\Omega} q(x)\\phi_9'(x) \\phi_8'(x) dx = -4.66666667$  \n",
    "A_loc$_{9,9} = \\int_{\\Omega} q(x)\\phi_9'(x) \\phi_9'(x) dx = 10.66666667$  \n",
    "A_loc$_{9,10} = \\int_{\\Omega} q(x)\\phi_9'(x) \\phi_{10}'(x) dx = -6$  \n",
    "A_loc$_{10,8} = \\int_{\\Omega} q(x)\\phi_{10}'(x) \\phi_8'(x) dx = 0.66666667$  \n",
    "A_loc$_{10,9} = \\int_{\\Omega} q(x)\\phi_{10}'(x) \\phi_9'(x) dx = -6$  \n",
    "A_loc$_{10,10} = \\int_{\\Omega} q(x)\\phi_{10}'(x) \\phi_{10}'(x) dx = 5.33333333$  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coupling the functions to solve specific problems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exact solution\n",
    "ue = lambda x: x**2  + x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Diffusion coefficient\n",
    "q = lambda x: 1 + 0.5*x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Forcing term\n",
    "f = lambda x: -(2.5 + 2*x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Interval end points\n",
    "a = -1\n",
    "b = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy\n",
    "from scipy.sparse import linalg\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Generate the computational mesh\n",
    "n_elements = 10      # specify number of elements\n",
    "\n",
    "# Compute nodes and connectivity matrix\n",
    "x, e_conn = oned_mesh(a,b,n_elements,'linear')\n",
    "\n",
    "n_nodes = len(x)            # number of nodes\n",
    "n_dofs = len(e_conn[1,:])   # degrees of freedom per element\n",
    "\n",
    "# Index to keep track of equation numbers\n",
    "ide = np.zeros(n_nodes, dtype=int)  # initialize\n",
    "\n",
    "# Mark Dirichlet nodes by -1\n",
    "i_dir = [0,n_nodes-1]  # indices of dirichlet nodes\n",
    "ide[i_dir] = -1        # dirichlet nodes are marked by -1\n",
    "\n",
    "# Number remaining nodes from 0 to n_equations-1\n",
    "count = 0\n",
    "for i in range(n_nodes):\n",
    "    if ide[i] == 0:\n",
    "        ide[i] = count\n",
    "        count = count + 1\n",
    "\n",
    "n_equations = count   # total number of equations\n",
    "\n",
    "\n",
    "# Initialize sparse stiffness matrix\n",
    "nnz = n_elements*n_dofs**2        # estimate the number of non-zero elements\n",
    "\n",
    "rows = np.zeros(nnz, dtype=int)   # row index\n",
    "cols = np.zeros(nnz, dtype=int)   # column index\n",
    "vals = np.zeros(nnz)              # matrix entries\n",
    "\n",
    "# Initialize the RHS\n",
    "c = np.zeros(n_equations)\n",
    "\n",
    "# Assembly\n",
    "r,w = oned_gauss(11) # Gauss rule accurate to degree (2n-1), n is the input of oned_gauss\n",
    "count = 0\n",
    "for i in range(n_elements):\n",
    "    # local information\n",
    "    i_loc = e_conn[i,:] # local node indices\n",
    "    x_loc = x[i_loc]    # local nodes\n",
    "    \n",
    "    #compute shape function on element\n",
    "    x_g,w_g,phi,phi_x,phi_xx = oned_shape(x_loc,r,w)\n",
    "    \n",
    "    # compute local stiffnes matrix\n",
    "    q_g = q(x_g)\n",
    "    A_loc = oned_bilinear(q_g,phi_x,phi_x,w_g)\n",
    "    \n",
    "    # local RHS\n",
    "    f_g = f(x_g)\n",
    "    c_loc = oned_linear(f_g,phi,w_g)\n",
    "    \n",
    "    # global\n",
    "    for j in range(n_dofs):\n",
    "        # for each row\n",
    "        j_test = i_loc[j]    # global node number\n",
    "        j_eqn = ide[j_test]  # equation number\n",
    "\n",
    "        \n",
    "        if j_eqn >= 0 :\n",
    "            #update RHS\n",
    "            c[j_eqn] = c[j_eqn] + c_loc[j]\n",
    "            \n",
    "            for m in range(n_dofs):\n",
    "                # for each column\n",
    "                i_trial = i_loc[m]    #global node number\n",
    "                i_col = ide[i_trial]  # equation number\n",
    "                \n",
    "                if i_col >= 0:\n",
    "                    # interio node: fill column\n",
    "                    rows[count] = j_eqn\n",
    "                    cols[count] = i_col\n",
    "                    vals[count] = A_loc[j,m]\n",
    "                    count = count + 1\n",
    "                else:\n",
    "                    # Dirichlet node: apply dirichlet condition\n",
    "                    u_dir = ue(x_loc[m])\n",
    "                    c[j_eqn] = c[j_eqn] - A_loc[j,m]*u_dir\n",
    "                    \n",
    "# Delete entries that weren't filled\n",
    "noz = len(rows) - (count + 1)    # read # of rows not affected by count\n",
    "rows = rows[:-noz]               # delete the last 'noz' entries\n",
    "cols = cols[:-noz]               # delete the last 'noz' entries\n",
    "vals = vals[:-noz]               # delete the last 'noz' entries\n",
    "\n",
    "# Assemble sparse stiffness matrix\n",
    "A = scipy.sparse.coo_matrix((vals,(rows,cols)), shape=(n_equations,n_equations)).tocsr()\n",
    "\n",
    "# Compute finite element solution\n",
    "ua = np.zeros(n_nodes)\n",
    "ua[i_dir] = ue(x[i_dir])           # apply prescribed values at Dirichlet nodes\n",
    "ua[ide >=0] = linalg.spsolve(A,c)  # solve the system at interior nodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4MAAAFNCAYAAAC66roXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABdjklEQVR4nO3dd3hUVf7H8fc3IZTQq0KABBASOkpAsIG9oYBtVRSxIba1LLvrytrlp2uXVUFsgGJdFbFi7ygEpPdepRN6SXJ+f9ygAQIkkJkz5fN6nnlm5t47934yDLnzzTn3HHPOISIiIiIiIvElwXcAERERERERCT8VgyIiIiIiInFIxaCIiIiIiEgcUjEoIiIiIiISh1QMioiIiIiIxCEVgyIiIiIiInFIxaCIiIiIiEgcUjEoIiIiIt6Z2QIz22pmmwrcnvGdSySWlfIdQEREREQk3znOuS/3t4GZlXLO5eyxLNE5l1vUgxR3e5FYpZZBEREREYlYZtbLzH4ysyfNbC1wr5kNMbOBZvaJmW0GTjSzpmb2rZmtN7OpZnZugX3stb23H0gkgqgYFBEREZFIdzQwD6gF9M9fdmn+44rAr8CHwOf529wMDDez9AL7KLj9j+GJLRLZVAyKiIjEEDN72cxWmtmUEtrfZ/ktLR/tsdzMrL+ZzTKz6Wb215I4nsS9Efmft123a/OXL3PO/dc5l+Oc25q/7APn3E/OuTygDVABeNg5t8M59zXwEXBJgX3/sb1zblvYfiKRCKZiUEREJLYMAc4owf09ClxeyPJeQD0gwznXFHizBI8p8aubc65KgdsL+csXF7JtwWV1gMX5heEuC4GUfWwvIqgYFBERiSnOue+BtQWXmVmj/Ba+cWb2g5llFGN/XwEbC1l1PXD/ri/fzrmVh5Jb5ADcAZYtA+qZWcHvtvWBpQfYh0hcUzEoIiIS+wYDNzvn2gJ9gedKYJ+NgL+YWZaZfWpmjUtgnyIH61dgM/APM0sys87AOajFWmS/NLWEiIhIDDOzCsAxwDtmtmtxmfx15wH3F/Kypc650w+w6zLANudcZv5+XgaOL5nUEsc+NLOCUz58AXxwoBc553bkjx76HPAvghbBns65GaGJKRIbzDm1mIuIiMQSM0sDPnLOtTCzSsBM51ztQ9hfZ6Cvc65LgWUzgDOccwssqDLXO+cqH1pyEREJJ3UTFRERiWHOuQ3AfDO7EP4YBbR1Cex6BHBS/uNOwKwS2KeIiISRWgZFRERiiJm9AXQGagArgHuAr4GBQG0gCXjTOVdY99DC9vcDkEEwbP8a4Grn3CgzqwIMJxikYxPQxzk3sUR/GBERCSkVgyIiIiIiInFI3URFRERERETikIpBERERERGROBR1U0vUqFHDpaWl+Y4hIiJhMG7cuNXOuZq+c0QLnSNFROJDSZ0fo64YTEtLIysry3cMEREJAzNb6DtDNNE5UkQkPpTU+VHdREVEREREROKQikEREREREZE4pGJQREREREQkDkXdNYOF2blzJ0uWLGHbtm2+o4gHZcuWpW7duiQlJfmOIiIiIiISNWKiGFyyZAkVK1YkLS0NM/MdR8LIOceaNWtYsmQJDRo08B1HRERERCRqxEQ30W3btlG9enUVgnHIzKhevbpahUVEREREiikmikFAhWAc07+9iIiIiEjxhawYNLN6ZvaNmU03s6lmdksh25iZDTCzOWY2ycyOClWeaDZixAimTZt2yPv59ttv6dKly363Wb9+Pc8999wfz5ctW8YFF1xwyMcWESmW4cMhLY220NZ3FBERkVgVypbBHOBvzrmmQAfgRjNrtsc2ZwKN82+9gYEhzBO1SqoYLIo9i8E6derwv//9LyzHFhEBgkKwd29YqPnmRUREQilkA8g455YDy/MfbzSz6UAKULCq6QoMc8454Bczq2JmtfNfG1Vee+01BgwYwI4dOzj66KN57rnnGD9+PFdffTVjxowhNzeX9u3b89Zbb5GWlkbXrl1Zt24dO3fu5MEHH6Rr164ADBs2jMceewwzo1WrVlx//fWMHDmS7777jgcffJB3332XRo0a/XHcd955h/vuu4/ExEQqV67M999/z7Zt27j++uvJysqiVKlSPPHEE5x44om75b333nupUKECffv2BaBFixZ89NFH3HHHHcydO5c2bdpw6qmncuONN9KlSxemTJmyz/0OGTKEkSNHsmXLFubOnUv37t155JFHwvfmi0hs6dcPtmzxnUJEJHa99hps3br7skaN4KSTgsdDhsDOnbuvz8iA448PHr/wwt77bNkSOnSAHTtg6NC91x91FLRtC5s3w+uv773+6KOhVStYvx7eeWfv9ccdB02bwqpVMGLE3utPPBGOOAKWLYOPP957/WmnQWpq8IfGzz/fe/3ZZ0OdOjBnDnzzzZ/LmzeHY47Ze/sYEZbRRM0sDTgS+HWPVSnA4gLPl+Qv260YNLPeBC2H1K9fP2Q5D9b06dN56623+Omnn0hKSuKGG25g+PDh9OzZk3PPPZd///vfbN26lcsuu4wWLVqQk5PD+++/T6VKlVi9ejUdOnTg3HPPZdq0afTv35+ffvqJGjVqsHbtWqpVq8a5555Lly5dCu2uef/99zNq1ChSUlJYv349AM8++ywAkydPZsaMGZx22mnMmjWrSD/Lww8/zJQpU5gwYQIACxYs+GPd/vY7YcIEfvvtN8qUKUN6ejo333wz9erVO8h3VETi2qJFvhOIiMS2v/8dfv9992WXXPJnMXjzzbBp0+7re/f+sxjs3Xvvfd5+e1AMbt9e+Pp77gmKwezswtc/9lhQDK5YUfj6558PisFFiwpf//rrQTE4a1bh60eODIrBSZMKX//NN0ExOGbM7uurVYPVqyFGx6gIeTFoZhWAd4FbnXMb9lxdyEvcXgucGwwMBsjMzNxrfUH3fTiVacv2PMyhaVanEvec03yf67/66ivGjRtHu3btANi6dSu1atUC4O6776Zdu3aULVuWAQMGAMF0CHfeeSfff/89CQkJLF26lBUrVvD1119zwQUXUKNGDQCqVat2wGzHHnssvXr14qKLLuK8884D4Mcff+Tmm28GICMjg9TU1CIXg/uzv/2efPLJVK5cGYBmzZqxcOFCFYMiclDy6tUjQQWhiEjJ+vhjWLwYrrkGxo+HvLzd15cr9+fjmTPB7fGVOzn5z8dLluy9/woVgvvy5QtfX7FicH/YYYWvr1QpuG/YsPD1VaoE9y1bFr6+atXgvmPHwtfv+l59yimFr8///k23bn+uf+UVuOsuWLkyyB2DQloMmlkSQSE43Dn3XiGbLAEKVgx1gWWhzBQKzjmuuOIKHnroob3WrV27lk2bNrFz5062bdtG+fLlGT58OKtWrWLcuHEkJSWRlpbGtm3bcM4Ve2TMQYMG8euvv/Lxxx/Tpk0bJkyYgNvzP28hSpUqRV6BXwJFmZphf/stU6bMH48TExPJyck54P5ERArzwYU3cvrTd5Ocs913FBGR2JCXB//8Z1Dg9e4NtWvvf/s6dfa/PiVl3+sSEva/PjFx/+uTkva/vnTp/a8vU2b/68uV2//65OQ/C98LL4R69XYvlGNMyIpBC6qal4Dpzrkn9rHZSOAmM3sTOBrIPtTrBffXghcqJ598Ml27duW2226jVq1arF27lo0bN5Kamkrv3r154IEHmD9/Pv/85z955plnyM7OplatWiQlJfHNN9+wMH+QhJNPPpnu3btz2223Ub169T+6iVasWJGNGzcWeuy5c+dy9NFHc/TRR/Phhx+yePFiTjjhBIYPH85JJ53ErFmzWLRoEenp6YwePfqP16WlpfHRRx8BMH78eObPnw+w32Pta7/jx48vybdTROLYtGUb6Fu6BdtvvJeLRwzSIDIiIiXh/fdh6tSgK2VCzMwsF3rp6cEthoWyZfBY4HJgsplNyF92J1AfwDk3CPgEOAuYA2wBrgxhnpBp1qwZDz74IKeddhp5eXkkJSXx7LPP8t1331GqVCkuvfRScnNzOeaYY/j666/p0aMH55xzDpmZmbRp04aMjAwAmjdvTr9+/ejUqROJiYkceeSRDBkyhIsvvphrr72WAQMG8L///W+3AWT+/ve/M3v2bJxznHzyybRu3ZqMjAz69OlDy5YtKVWqFEOGDNmt5Q7g/PPPZ9iwYbRp04Z27drRpEkTAKpXr86xxx5LixYtOPPMM7nxxhv/eM0NN9xwwP2KiBysvDzHv0dMpkq5JM749+3w1B2MMxvnO5eISFTLy4P774cmTeCii3yniT6//Ra0qB4VmzPgWVG6FEaSzMxMl5WVtduy6dOn07RpU0+JJBLoMyAS/d4Ys4h/vTeZJy5qzXlH1QXAzMY55zI9R4sahZ0jRSTOffBBcB3csGFw+eW+00Sf1q2DrqL5PeoiRUmdH9VOLCIi3q3etJ2HP51Bh4bV6H7kfq7lEBGR4qlSBS64IBgtVIovIwNmzPCdImRUDIqIiHf/98l0tuzI4cFuLYo9kJaIiOxHp07BvH2lwjKjXOxJT4f584MpM2KQikEREfFq9Nw1vDd+Kded0IgjalX0HUdEJDY4B88+C2vW+E4S3dLTg+su58zxnSQkVAyKiIg3O3Ly+PeIydSrVo6bTjrCdxwRkdgxahTcdBO8V9jsblJk+QM9MnOm3xwhovZiERHx5oUf5jF31WZe6dWOskmJvuOIiMQG5+C++4KBT664wnea6NasGXz1VcyOJqpiUEREvFi0ZgsDvprNWS0P58SMWr7jiIjEjq++gl9+geeeCyZpl4NXrhycdJLvFCGjbqIxYuTIkTz88MMH9dq0tDRWr159SMfv1asX//vf//a7zbfffsvPP//8x/NBgwYxbNiwQzquiEQn5xx3j5xCqQTj7i7NfccREYkdu1oFU1Lgqqt8p4kNP/0EQ4f6ThESahmMATk5OZx77rmce+65vqPs17fffkuFChU45phjAOjTp4/nRCLiy2dTfufbmau4q0szDq9c1nccEZHYsWULVKsGF18MZcr4ThMbXn8dhg+Hnj0hxka8js+WweHDIS0NEhKC++HDD2l3CxYsoGnTplx77bU0b96c0047ja1btwLQuXNndk0AvHr1atLS0gAYMmQI3bp145xzzqFBgwY888wzPPHEExx55JF06NCBtWvXAjB37lzOOOMM2rZty/HHH8+M/HlOevXqxe23386JJ57IP//5T4YMGcJNN90EwIoVK+jevTutW7emdevWf7TGdevWjbZt29K8eXMGDx68358pNzeXXr160aJFC1q2bMmTTz4JwIQJE+jQoQOtWrWie/furFu3bq/XFmxpzMrKonPnzixYsIBBgwbx5JNP0qZNG3744QfuvfdeHnvssf3ut3Pnzvzzn/+kffv2NGnShB9++OHg/pFEJGJs2p7DfR9Oo1ntSlzRMdV3HBGR2FK+fDDR/A03+E4SO9LTITsbVq70naTExV8xOHw49O4NCxcGzegLFwbPD7EgnD17NjfeeCNTp06lSpUqvPvuuwd8zZQpU3j99dcZM2YM/fr1Izk5md9++42OHTv+0X2yd+/e/Pe//2XcuHE89thj3FDgP/asWbP48ssvefzxx3fb71//+lc6derExIkTGT9+PM2bB12wXn75ZcaNG0dWVhYDBgxgzX6GGp4wYQJLly5lypQpTJ48mSuvvBKAnj178p///IdJkybRsmVL7rvvviK9P2lpafTp04fbbruNCRMmcPzxx++2fn/7zcnJYcyYMTz11FNFPp6IRK4nv5jFio3b6N+9BaUS4+80JCISMtOmwezZweMYa8HyKoZHFI2/s3C/fkHzeUFbtgTLD0GDBg1o06YNAG3btmXBggUHfM2JJ55IxYoVqVmzJpUrV+acc84BoGXLlixYsIBNmzbx888/c+GFF9KmTRuuu+46li9f/sfrL7zwQhIT9x597+uvv+b6668HIDExkcqVKwMwYMAAWrduTYcOHVi8eDGzd/2yKETDhg2ZN28eN998M5999hmVKlUiOzub9evX06lTJwCuuOIKvv/++yK9P/tzoP2ed955QNHfVxGJXFOWZvPKT/O5tH19jqxf1XccEZHYcsst0Lkz5OT4ThJb0tOD+/weerEk/q4ZXLSoeMuLqEyBPtmJiYl/dBMtVaoUeXl5AGzbtm2fr0lISPjjeUJCAjk5OeTl5VGlShUmTJhQ6DHLly9f5HzffvstX375JaNHjyY5OZnOnTvvlaegqlWrMnHiREaNGsWzzz7L22+//UdX0QPZ3898MHa9L4mJieTol5tI1MrLc/x7xBSqlS/NP07P8B1HRCS2/PwzfPklPPYYlIq/r/ghVa9eMKqoWgZjQP36xVt+iNLS0hg3bhzAAUfb3FOlSpVo0KAB77zzDhCMvjdx4sQDvu7kk09m4MCBQHDt34YNG8jOzqZq1aokJyczY8YMfvnll/3uY/Xq1eTl5XH++efzwAMPMH78eCpXrkzVqlX/uG7v1Vdf/aM1b18/c8HushUrVmTjxo17bV/U/YpIdHtj7CImLF5Pv7ObUjk5yXccEZHY8sADUKMGaIC+kpeQANOnw0GO3B/J4q8Y7N8fkpN3X5acHCwPgb59+zJw4ECOOeaYg5q+Yfjw4bz00ku0bt2a5s2b88EHHxzwNU8//TTffPMNLVu2pG3btkydOpUzzjiDnJwcWrVqxV133UWHDh32u4+lS5fSuXNn2rRpQ69evXjooYcAGDp0KH//+99p1aoVEyZM4O67797rtffccw+33HILxx9//G7dWM855xzef//9PwaQKago+xWR6LVq43b+8+kMOjasTrc2Kb7jiIjEljFj4LPPoG/fYAAZKXmpqZAUe3/INOec7wzFkpmZ6XaNzrnL9OnTadq0adF3Mnx4cI3gokVBi2D//tCjRwknlXAq9mdARMLqtrcm8PGk5Xx66/E0qlmhyK8zs3HOucwQRosphZ0jRSQODB4czC04YwZUrOg7TWz69Vd49VV4/PGImLKjpM6P8dcyCEHht2AB5OUF9yoERURC5ue5q3n/t6Vc16lhsQpBCZhZopn9ZmYf+c4iIhGqd2+YO1eFYCjNnQvPPgtz5vhOUqLisxgUEZGw2J6Ty79HTKF+tWRuPPEI33Gi1S3AdN8hRCRCTZsW3Jct6zdHrIvREUVVDIqISMgM/m4e81Zt5v6uzSmbtPdUOLJ/ZlYXOBt40XcWEYlAEyZA8+bw8su+k8S+XcVgjI0oGjPFYLRd+yglR//2IpFp4ZrNPPPNHM5uWZvO6bV8x4lWTwH/API85xCRSPTgg1CpEnTv7jtJ7KtQAVJSVAxGorJly7JmzRoVBXHIOceaNWsoq64RIhHFOcfdH0wlKTGBu7o08x0nKplZF2Clc27cAbbrbWZZZpa1atWqMKUTEe+mTIF334W//hWqVvWdJj40bQpr1/pOUaJiYkbKunXrsmTJEnQSjE9ly5albt26vmOISAGfTP6d72at4p5zmnF4Zf2x5iAdC5xrZmcBZYFKZvaac+6yghs55wYDgyEYTTT8MUXEiwcfDFqrbrvNd5L48emnUComyqc/xMRPk5SURIMGDXzHEBERYOO2ndz/0VSa16nE5R1SfceJWs65fwH/AjCzzkDfPQtBEYlTGzbA11/DzTdDtWq+08SPGCsEIUa6iYqISOR44otZrNy4nf7dW1IqUacZEZESV6kSzJsHd9zhO0l8mT0bunWDsWN9JykxOkuLiEiJmbI0m6E/L+Cyo1NpU6+K7zgxwzn3rXOui+8cIhIBsrMhNzfoIlqpku808aV0afjgAxg/3neSEqNiUERESkRunqPf+5OpVr4MfU9P9x1HRCQ23XILtG8PeRpkOOzq1YNy5WJqRFEVgyIiUiJeH7OIiUuyuatLUyqXS/IdR0Qk9sydC6+9Bp07Q4K+xoddQgI0aRJTE8/rUyQiIods5cZtPPLZDI49ojrntq7jO46ISGz6v/+DpCTo29d3kviVnq6WQRERkYL+7+PpbN+ZxwNdW2BmvuOIiMSe+fNh2DDo3Rtq1/adJn61awd16sRMN10VgyIickh+mrOaEROW0adzIxrWrOA7johIbHrxxaCb4j/+4TtJfOvbF374IWa66cbGTyEiIl5sz8nlrhFTSK2ezA2dG/mOIyISu+6/H378EVJSfCeRGKJiUEREDtrz381j3urNPNC1BWWTEn3HERGJTTk5kJgYdFEUv3buDP4dnnzSd5ISoWJQREQOyoLVm3nmmzl0aVWbE5rU9B1HRCQ2LVkCaWnw2We+kwgEA/j8/jv89pvvJCVCxaCIiBSbc467PphCmcQE7urSzHccEZHY9cgjsGIFZGT4TiK7xNCIoioGRUSk2D6evJwfZq/mb6c14bBKZX3HERGJTcuXw+DB0LNn0DookSEjI5hr0DnfSQ6ZikERESmWDdt2cv+H02iZUpnLO6b5jiMiErsefTS4XvDOO30nkYLS02HDhqDFNsqV8h1ARESiyxOfz2LVpu28eEUmiQmaU1BEJCTWroVBg6BHD2ik0ZojSmYmXHghbNvmO8khUzEoIiJFNnlJNsNGL+DyDqm0qlvFdxwRkdhVrRp8+inUq+c7ieypY8fgFgNUDIqISJHk5jn6jZhM9Qpl6Ht6uu84IiKxr1Mn3wlkf3bsgNKlfac4JLpmUEREimT4rwuZtCSbu7o0o1LZJN9xRERi1/33wy23xMQAJTHrlFOgWzffKQ6ZikERETmglRu28ehnMznuiBqc06q27zgiIrFr7dpg4JgVK8B0XXbEqlEjJqaXUDEoIiIH9ODH09mem8cD3Vpg+nIiIhI6Tz4JmzbBXXf5TiL7k5EBCxZE/SAyKgZFRGS/fpi9ipETl3FD50Y0qFHedxwRkdi1bh0MGAAXXADNm/tOI/uTng55eTBnju8kh0TFoIiI7NO2nbncNWIKadWT6dNJQ5uLiITUgAHB/HX//rfvJHIgGRnBfZR3FdVooiIisk+DvpvLgjVbePXq9pRNSvQdR0Qktl1yCVStCq1b+04iB9KkCdx2GzRo4DvJIVExKCIihZq/ejPPfTOXc1vX4fjGNX3HERGJfU2aBDeJfOXLwxNP+E5xyNRNVERE9uKc464RUyhTKoF/d2nqO46ISGzbuBEuvxymT/edRIpj2zZdM7gvZvayma00syn7WN/ZzLLNbEL+7e5QZRERkeL5cNJyfpyzmr+fkU6timV9xxERiW3PPguvvRYUhRI9+vaFtm2jej7IULYMDgHOOMA2Pzjn2uTf7g9hFhERKaLsrTt54KNptKpbmR5Hp/qOIyIS2zZtgscfhzPOgPbtfaeR4khPDwb8+f1330kOWsiKQefc98DaUO1fRERC4/HPZ7Jm03b6d2tJYoLmFBQRCalBg2D1arhbneSiTnp6cB/FI4r6vmawo5lNNLNPzUyTqYiIeDZpyXpe/WUhPTum0bJuZd9xRERi25Yt8OijcMop0LGj7zRSXDEwvYTP0UTHA6nOuU1mdhYwAmhc2IZm1hvoDVC/fv2wBRQRiSe5eY47359MzQpluP00jWYnIhJyublwzTVw1lm+k8jBqFsXypWDGTN8Jzlo3opB59yGAo8/MbPnzKyGc251IdsOBgYDZGZmRu8VmiIiEezV0QuYsnQD/73kSCqVTfIdR0Qk9lWsCP37+04hByshIejm2zR6R932Vgya2eHACuecM7P2BF1W1/jKIyISz1Zs2MZjn8/i+MY16NKqtu84IiKx7913oXRp6NIFTNdnR62ePX0nOCShnFriDWA0kG5mS8zsajPrY2Z98je5AJhiZhOBAcDFzkXxuKwiItFo+HBIS6NW5WRGPdWTJ3ZOxfSlREQktLZvh1tuCa4XlOi2ejV8+mnwbxqFQtYy6Jy75ADrnwGeCdXxRUTkAIYPh969YcsWDEjZsBL+djNULAM9evhOJyISu155BZYuhaFD1SoY7b74Ai69FCZNgpYtfacpNt+jiYqIiC/9+gUj2RW0ZUuwXEREQmPHDnjoITjmGDjpJN9p5FBF+YiiPkcTFRERnxYtKt5yERE5dEOHBr9nBw9Wq2AsaJI/+naUjiiqlkERkTi19fCUwldoCh8RkdCpWBHOPx9OO813EikJ5csHU0xEacugikERkTj0e/Y27u/Yg21JZXZfkZysYc5FRELp4ovhf/9Tq2AsychQMSgiItEhL8/xt3cmMKJZZ7IHPAepqcGXktTUoNuSBo8RESl5OTkwZEjUjjop+/HII8GgQFFI1wyKiMSZF36Yx09z1vDweS05rP0Z0Ocq35FERGLf66/DlVdC9epwzjm+00hJOvJI3wkOmloGRUTiyKQl63l01EzObHE4f2lXz3ccEZH4kJsLDz4IrVsHk8xLbFm3Dl54AWbP9p2k2FQMiojEic3bc7jlzQnUrFiGh85rqcnlRUTC5a23gkLhrrt0rWAs2rQpmLf3q698Jyk2dRMVEYkT9304lQVrNvPGtR2oklzadxwRkfiwq1WweXPo3t13GgmFlJRgALYoHERGxaCISBz4aNIy3s5awk0nHkGHhtV9xxERiR8rVwbTD/TtCwnqlBeTEhIgPV3FoIiIRJ4l67bwr/cm06ZeFW45pbHvOCIi8aV2bRgzBpzznURCKT0dfv3Vd4pi058nRERiWG6e47a3JuAcDLj4SJIS9WtfRCRspkyBNWuC6wTVKhjb0tNhwQLYts13kmLRp1JEJIY9+80cxi5YxwPdmlO/erLvOCIi8SMvDy67DE4/3XcSCYebb4bVq6FsWd9JikXdREVEYtS4hWt5+qvZdGtTh+5H1vUdR0Qkvnz4IUycCMOG+U4i4VA9Oq/HV8ugiEgM2rBtJ7e8OYE6Vcpyf7cWvuOIiMQX5+D++6FRI7jkEt9pJBycg3vvhfff952kWFQMiojEGOcc/35/Csuzt/HUX46kUtkk35FEROLLxx/D+PHQrx+UUke8uGAGL78M773nO0mxqBgUEYkx7/+2lJETl3HryY1pm1rVdxw5SGZW1szGmNlEM5tqZvf5ziQiRTR7NqSlBdcMSvyIwuklVAyKiMSQhWs2c9eIKbRPq8YNJx7hO44cmu3ASc651kAb4Awz6+A3kojs18qVwf2tt8LXX0OSembElYwMmDEjqqYRUTEoIhIjdubm8dc3J5CYYDx5cRsSE8x3JDkELrAp/2lS/i16vmGIxJPcXPj734NiYOHCoMtggwa+U0m4pafDxo3w++++kxSZikERkRjx1JezmLh4PQ+f34qUKuV8x5ESYGaJZjYBWAl84ZyLvhmNRWLdunVw1lnw2GPQowfUqeM7kfiSkQHJybB4se8kRaYrWkVEYsDouWt47tu5/CWzHme1rO07jpQQ51wu0MbMqgDvm1kL59yUgtuYWW+gN0D9+vXDH1Iknk2bBl27Bq2BL7wA11zjO5H4dOKJQctgQvS0t0VPUhERKdT6LTu47a0JNKhenrvPaeY7joSAc2498C1wRiHrBjvnMp1zmTVr1gx3NJH49sQTwZf/b75RISiQmBhVhSCoGBQRiWrOOe54dzJrNm/n6YuPpHwZdfiIFWZWM79FEDMrB5wCzPAaSkSCwUFWrw4eDxgAWVlw7LF+M0nkeOAB6NvXd4oiUzEoIhLF3hy7mM+m/s7fT0+nZd3KvuNIyaoNfGNmk4CxBNcMfuQ5k0h827QJLroIOnWCLVuC68Pq1vWdSiLJtGnw7ru+UxSZ/oQsIhKl5qzcxH0fTuW4I2pwzXENfceREuacmwQc6TuHiOSbPx+6dYMpU+CRR6CcBuqSQqSnw1tvwdatUfEZUcugiEgU2p6Ty1/f+I3k0qV44qLWJGgaCRGR0Pn6a2jXDhYtgk8+gb/9LZg+QmRPGRlBV+LZs30nKRIVgyIiUejRz2YybfkGHjm/FbUqlfUdR0QkdjkHd94JtWrBmDFw+um+E0kkS08P7mfO9JujiNRNVEQkynw/axUv/jifyzukckqzw3zHERGJTdu3w86dUKECvPdecF+pku9UEumaNAlueXm+kxSJikERkSiyetN2bn97Io1rVaDf2U19xxERiU3Ll8P550PNmjBihCaSl6IrXz5qWgVBxaCISNRwzvGP/01iw7advHp1e8omJfqOJCISe8aMge7dYf16GDpU1wZKTNM1gyIiUWLY6IV8PWMld56ZQdPa6qokIlLihg2DE06A0qVh9Gi44ALfiSQaPfsstGgRXG8a4VQMiohEgRm/b6D/J9M5Mb0mVxyT5juOiEjsyc6Gf/wDjjkGxo6FVq18J5JoNnUq/P677xQHpG6iIiIRbtvOYBqJSmWTePTC1pi6LImIlJz166FiRahcGX74AdLSICnJdyqJZrtGFJ0xA2rX9pvlANQyKCIS4f7vk+nMWrGJxy9qTY0KZXzHERGJHVOmQNu2cO+9wfPGjVUIyqHLyAjuZ8zwm6MIVAyKiESwL6etYNjohVxzXAM6NanpO46ISOx4/33o0AG2bIGzzvKdRmJJSkrUjCqqbqIiIhFqxYZt/P1/E2lWuxJ/PyPddxw5SGbWs4ibTnDOTQppGBEJ5n978EG45x5o1y4oClNSfKeSWGIGl10WtDRHOBWDIiIRKC/P8be3J7J1Zy4DLjmSMqU0jUQUa1DE7RaEMoSI5JsxA/r3h5494fnnoWxZ34kkFg0a5DtBkagYFBGJQC/+OI8f56zmofNackStCr7jyCFwzt3nO4OIEAwUU6UKNGsG48cH9xqQS0IpJwcSEoJbhIrcZCIicWrykmweHTWTM5ofzsXt6vmOIyXEzGqZWXczu9HMrjKz9mam87BIOHz1FTRqBG+9FTxv3lyFoITWJ59AcnIwSFEE00lIRCSCbN6ew1/f/I3q5cvw8PktNY1EDDCzE81sFPAxcCZQG2gG/BuYbGb3mVklnxlFYpZz8PTTcPrpwRD/mZm+E0m8qFMHdu6M+EFk1E1URCSC3P/hNBas2czr13SgSnJp33GkZJwFXOucW7TnCjMrBXQBTgXeDXcwkZi2bRtcfz0MGQLdusGwYcF8giLh0KRJcK9iUEREiuLjSct5K2sxN57YiI6NqvuOIyXEOff3/azLAUaEL41IHPnss6AQvPdeuOuuiL5uS2JQcjLUrx/xcw3qf4WISARYun4r/3pvEq3rVeHWU5r4jiMhYGavmlnlAs/TzOwrn5lEYlJ2dnDfrRtMnBhMIaFCUHxIT4/4lkH9zxAR8Sw3z3HbmxPIzXMMuLgNSYn61RyjfgR+NbOzzOxa4HPgKb+RRGLMK69AWhr89lvwvFUrr3EkzvXsCT16+E6xX+omKiLi2XPfzGHMgrU8cVFrUquX9x1HQsQ597yZTQW+AVYDRzrnfvccSyQ25ORA377BYDGnnAKpqb4TiQQTz0c4/flZRMSjcQvX8dRXs+napg7dj0zxHUdCyMwuB14GegJDgE/MrLXXUCKxYM2aYLTQp5+GW2+FTz+FatV8pxIJRrNdtgzWrfOdZJ9C1jJoZi8TjJC20jnXopD1BjxNMMraFqCXc258qPKIiESaDdt2csubv1G7clke6NZC00jEvvOB45xzK4E3zOx9gqLwSK+pRIpi3jyYPn33ZdWrQ4cOweNvvoEtW3Zff9hhf07l8OWXsH377uvr1IEj8z/+n30Gubm7r69fH1q2DL5Qf/LJ3pkaNoSmTeG55+DHH4Muor16HdSPJxISS5dCvXrBZ/T6632nKVQou4kOAZ4Bhu1j/ZlA4/zb0cDA/HsRkbhw94gpLM/extvXdaRS2STfcSTEnHPd9ng+xsx03pPIl5MDp50Gc+fuvvzEE+Hrr4PHvXvDnDm7rz/3XPjgg+DxZZfBihW7r7/0Uhg+PHh8/vl7F5PXXQeDBgXFYJcue+fq2xcefRT+9a9gsJiWLQ/qxxMJmZQUKF8+okcUDVkx6Jz73szS9rNJV2CYc84Bv5hZFTOr7ZxbHqpMIiKR4v3fljBiwjJuP7UJbVOr+o4jIWRm/waec86t3XOdc26HmZ0EJDvnPgp/OpH92NVSV6oUvPde0NUtOfnP9QXn7Hv33b1b/qpU+fPxp58GRWVBBbtyfv895OXtvr5mzeDeDMaM2Tvf4Yf/mU+FoEQis2C+wQgeUdTnADIpwOICz5fkL1MxKCIxbeGazdw1Yirt06px44lH+I4joTcZ+NDMtgHjgVVAWYKeMW2AL4H/85ZOpDDr1sEllwRF1qOPHnhUzgOtP/IAvaHbtt33OjNo127/rxeJVBkZMHq07xT75HMAmcIujnGFbmjW28yyzCxr1apVIY4lIhI6O3PzuOXNCZjBkxe3ITFB1wnGOufcB865Y4E+wFQgEdgAvAa0d87d5pzTyU0ix7Rp0L590AW0ieY9FTkk6emwcCFs3eo7SaF8tgwuAeoVeF4XWFbYhs65wcBggMzMzEILRhGRiDZ8OPTrR6lFi3imYg1W3nE3KVVO951Kwsg5NxuYbWblnXObfecRKdTIkcH1fcnJwaAwxx7rO5FIdOvePRjsKEL5bBkcCfS0QAcgW9cLikhMGj48GFxh4ULMOepuWMVRD/7zz4ETJC6YWUczmwZMz3/e2sye8xxL5E8rV8LFFwctGVlZKgRFSkKrVnD55VCunO8khQpZMWhmbwCjgXQzW2JmV5tZHzPrk7/JJ8A8YA7wAnBDqLKIiHjVr9/eo+Rt2RIsl3jyFHA6sAbAOTcROMFnIBEAduwI7mvVglGjgsFc6tb1m0kklmRlwaRJvlMUKpSjiV5ygPUOuDFUxxcRiRRu0aJCL5Jm0aJwRxHPnHOL95hPMndf24qExfz50LVrME1Dz55w/PG+E4nEnvPPh+OOi8geQcVqGTSzqmbW3MwampnPLqYiIlEhL8+xtvphha+sXz+8YcS3xWZ2DODMrLSZ9SW/y6iIF19/HYzSuXhxMEG8iIRGRkbETi9xwILOzCqb2Z1mNhn4BXgeeBtYaGbvmNmJoQ4pIhKtHvt8Jvd16MHOMmV3X5GcDP37+wklvvQh6BGTQjCIWhvUQ0Z8cA7++99gIvlatYI5/E7XgFYiIZOeHhSDLvLGwSxK697/COYDPN45l+6cO845l+mcqwc8DHQ1s6tDmlJEJAq9OWYRz307l/JX9qTUSy9CamowX1ZqKgweDD16+I4oYeScW+2c6+GcO8w5V8s5d5lzbo3vXBKHfv0V/vpXOPts+OUXaNzYdyKR2JaeDps2wbJCJ07w6oDXDDrnTt3PunHAuBJNJCISA36YvYp+I6ZwfOMa3N+1BZbYSsVfnDOzVyhkPl3n3FUe4kg82rkTkpKgQ4dgoJhTToEEXfUjEnIZGcH9zJmQkuI3yx6K/Btgz9Y/M0s0s3tKPpKISHSb+ftGbnhtPEfUrMCzPY4iKVFftgSAj4CP829fAZWATV4TSfwYMyb4Qvrrr8Hz005TISgSLu3aBdfoZmb6TrKX4vwWONnMPjGz2mbWguD6wYohyiUiEpVWbtjGVUPGUrZ0Ii9f2Y5KZZN8R5II4Zx7t8BtOHAR0MJ3LokDw4bBCSdAXl7EznUmEtMqVYITTwzuI0yRp5Zwzl1qZn8BJgNbgEuccz+FLJmISJTZsiOHq4dmsXbzDt6+riMpVfSlS/arMaAhZSV0cnLgH/+AJ58Mvoi+/TbUqOE7lUh8+uYbWLECLr7Yd5LdFKebaGPgFuBdYAFwuZklhyiXiEhUyc1z3PLmBKYsy+a/lxxJy7qVfUeSCGNmG81sw64b8CHwT9+54tLvv8OWLcHj3NyIHOGvRAwZEhSCf/1rcI2gCkERf158Ef71L98p9lKcbqIfAnc5564DOgGzgbEhSSUiEmX6fzydL6at4O4uzTilmebrkr055yo65yoVuDVxzr3rO1fcGTs2uG7nppuC58OGQcuW8MwzkJ3tN1tJyckJ7q+8Ej75BJ5+Ohg4RkT8SU+HhQth61bfSXZTnGKwvXPuKwAXeBzoFpJUIiJRZOjPC3j5p/n0OiaNK49t4DuORAEzu9d3hrj06qtw/PFBYXTLLcGymjWDeT9vvhnq1IGrrw4KxmhtLRwxApo1g6VLITERzjzTdyIRgWAAJ+dg9mzfSXZTlEnnjwNwzm3Yc51zbraZVcofUEZEJO58NX0F9304lVOaHsZdXZr5jiPR41zfAeJKTg707Qs9e0LHjkGx17p1sK5Ll2CkzaysYPqXN9+EPn12f200yMuD+++H7t2hShXfaURkT+npwf2MGX5z7KEoLYPnm9nPZna3mZ1tZu3N7AQzu8rMXiUYKlujJIhI3JmyNJubXv+NZnUqMeCSNiQmmO9IEj30YQmn5cvh5ZeDrqGff174tXNt28LgwcGk0K++Cmawfj3Uqwc33giTJoU9dpFt3AgXXAD33BMUvN9/H3FzmYnEvcaNg/uZM/3m2ENRJp2/zcyqAhcAFwKHA1uB6cAgjSgqIvFo2fqtXDVkLFWTk3j5inYkly7y4MwiAG19B4gLCxdC/fpBQTd1KtSufeDXVK4c3AA2bw7m43vpJXjuuaBVsU8fuPDCyJqi4a674IMPgsFibrklKGRFJLIkJ8PcucHvpAhSpGsGnXPrgDrAfGA0MAHYBpwcsmQiIhFq47adXDVkLFt25PLyle2oVams70gS4cysrJldYGZPm9k7wBAz+4eZNfedLWZ98AG0aAFPPRU8L0ohuKeUFBg6NLj+7oknYM0auOIKmD8/WL9zZ4nFPSi5ucH9fffBV1/BrbeqEBSJZA0bQqnI+uNxcQaQ2VTglgOcCaSFIJOISMTKyc3jxtd/Y/bKTTzX4ygyDo+8CWQlsuQPFvMT0BH4FXgeeJvgXPqwmX1hZq38JYwxu66d69YtGLDhwgsPfZ/Vq8NttwXX+owdGwzQAnD55XDSScH8fTt2HPpxisq5YITQ444LRiasXBk6dw7f8UXk4Hz3XTDVSwQNUFWcSecfL/jczB4DRpZ4IhGRCOWc4+6RU/l+1ioeOq8lJzSp6TuSRIexzrl797HuCTOrhSafLxmbNgUtd++9B5ddFlwDWJLdOc2CaSl2ycwMpqT4y1+gVq1gJNLevSEtreSOuadt24KuqkOHBgXvrtZBEYl8U6fCf/8L//xnxFzXW5yWwT0lAw1LKoiISKQb/P08Xv91EX06NeKS9vruLkXjnPsYgq6ie64zsxrOuZXOuaxC1tUzs2/MbLqZTTWzW8KRN6pNmAAffQSPPx7MHxjq6/r69g2uAfrkE+jQAf7zn6AAhaCFsqRHIl22DDp1CgrBe++Fd9+FChVK9hgiEjq7RhSNoEFkilwMmtlkM5uUf5sKzASeDl00EZHI8cnk5Tz06QzOblWbf5ye7juORKexZtZh1xMzOx/4eT/b5wB/c841BToAN5qZ5i8pzOLFwf1xx8G8eXD77eG7dm7XXH4ffAALFgTX7QF8+ik0aBB0WV26tGSO1bNn0LLw3nvByKEJh/I3fREJu4yM4D6CppcozhWMXQo8zgFWOOeiZPIdEZGDN37ROm57awJH1a/C4xe2JkFTSMjBuRR42cy+JRiUrTpw0r42ds4tB5bnP95oZtOBFGBa6KNGCeeCLld9+wYtgqed5rfrVb16fz6uWhWaNw+Ktvvvh3POCbp3nnpq8Yu43Nyg6Bw4MOgm2rJlyeYWkfCoUydozY+glsHiXDO4MJRBREQi0aI1W7h2aBaHVSrLCz0zKZuU6DuSRCnn3GQz6w+8CmwETnDOLSnKa80sDTiSYAAaAdi+Ha6/Hl55Bbp2DaZ9iCTHHAOffRZ0I33hhWCew/Hjg5ZLCPKXKbP/feTkBIXuqlXw2mt/zlMmItHJLGgdXLPGd5I/qH+BiMg+ZG/ZyZVDxpCT53jlynZUr3CAL24i+2FmLwG3Aq2AK4EPzezGIryuAvAucKtzbkMh63ubWZaZZa1ataqEU0eoXdfOvfIK3H130G2yYkXfqQrXqBE8/HDQlfXTT4MWvh074Igj4JJLgtEFCxtZcM0aOP30YNTQww4LrkEUkeg3enTwx50IoWJQRKQQO3LyuO61LBat3cLzl7elUU0N0iCHbApwonNuvnNuFMF1gEft7wVmlkRQCA53zr1X2DbOucHOuUznXGbNmnEywu3HH8OUKcEAKvfdFx3XzpUp8+eUFFu2wPnnBy2HnTsHy59+GtatC9ZPngzt2sFPP8GQIcEch4nqlSASE6J4nkERkbjgnOOO9ybxy7y1PHJBKzo0rO47ksQA59yTzv3ZBOScy3bOXb2v7c3MgJeA6c65J8KRMeItye9Ve801MH06nHee3zwHq0oVeOqpYGCZV14J5gm89VaYNCloNTz77ODawO++C6bKEJHYMWkSnHtu8DssAqgYFBHZw4Cv5vDe+KXcdkoTuh9Z13cciXJm9qGZnZPfyrfnuoZmdr+ZXVXIS48FLgdOMrMJ+bezQh44EuXkBMVS06YwZ05w3U3BwVqiVXIy9OoFv/wSfEE84QQoXRreeAOysuDoo30nFJGS5hx8+GHwfz4CRFY7pYiIZ+//toQnv5zFeUel8NeTj/AdR2LDtcDtwFNmthZYBZQF0oC5wDPOuQ/2fJFz7kdAQ9euWRNM6v7VV3DLLaGd0N2ngiOEHnusvxwiElq7BoKKkBFFVQyKiOT7Zd4a/vG/SXRoWI2Hz2uFhWueMolpzrnfgX8A/8gfFbQ2sBWY5Zzb4jNbxJs8ORgpdOnSYDTOK6/0nUhE5NAkJ0NqqopBEZFIMnfVJq57dRz1qyXz/GWZlC6lXvRS8pxzC4AFnmNEj8GD/7x2rkMH32lEREpGenrETDyvbzsiEvfWbNrOla+MpVSC8Uqv9lRO3uvSLpFDZmYdzWysmW0ysx1mlmtme00VEffy8mD58uDx44/DuHEqBEUktrRvDzVq+E4BqGVQROLctp25XDssixUbtvFG7w7Ur57sO5LErv8CFwPvAJlAT0AXpha0cSNcfnnQPXTChGDuwNq1facSESlZDzzgO8EfVAyKSNzKy3P87e2JjF+0nud6HMVR9av6jiQxzjk3x8wSnXO5wCtm9rPvTBFj7tzg+sAZM4IWwQqa21NEJNTUTVRE4tajn8/k48nL+deZGZzVUq0PEnJbzKw0MMHMHjGz24DyvkNFhC++CCZZX74cRo0KRg3VAE4iEquys+Goo4KBsTxTMSgicemNMYsY+O1cLj26Pr1PaOg7jsSHy4FE4CZgM1APON9rokjgHPTvDykpMHYsnHyy70QiIqFVqRLMng0TJ/pOom6iIhJ/vp+1in+PmMIJTWpy/7nNNYWEhIVzbmH+w63AfT6zRIRt24JblSrwzjtQrpy6hopIfDCLmBFFVQyKSFyZ+ftGbhg+nsa1KvDspUdSKlEdJCQ8zGw+4PZc7pyLv6bppUuhe/egEBw1CmrW9J1IRCS80tPhxx99p1AxKCLxY+WGbVz5yhjKl0nk5V7tqFhWU0hIWGUWeFwWuBCo5imLP6NHw3nnwaZN8NprujZQROJTRga8/jps2RJMRO+J/iQuInFhy44crh6axfqtO3npinbUqVLOdySJM865NQVuS51zTwEn+c4VVi+9BJ07Q/ny8MsvweihIiLxqEMH+Mtfgil1PFLLoIjEvNw8x1/f+I2py7J5oWcmLVIq+44kccjMjirwNIGgpbCipzjht2kT3HcfdOoEb74J1eKvUVRE5A+nnhrcPFMxKCIx78GPp/Hl9JXcd25zTm56mO84Er8eL/A4B5gPXOQpS/isXRtMHl+hAnz/PdStC6X09UNEBICdOyHJ32Ur+m0sIjFtyE/zeeWnBVx1bAOuOCbNdxyJQ2ZWP/9hrwKLXSHr1zvnNoQrV1jk5MCxx8JZZwUTyael+U4kIhI5MjODgWSGD/cWQcWgiMSsL6et4P6PpnFqs8Pod3ZT33Ekfg2lkFFE8+0aPcUBQ4Bh4QgUNh9+GAydfp9m0hAR2UuNGt6nl1AxKCKxZfhw6NcPt2gRTSvV5Iau13HDfXeTmKARC8UP59yJvjN4M3Bg0C30vPN8JxERiTwZGcH0Es55G1lZo4mKSOwYPhx694aFCzHnSMleyd/eeZzkd97ynUwk/syeDV98Efyf1DWCIiJ7S0+HzZuDuVc9UTEoIrGjX79gvp4CbOuWYLmIhNcLLwRF4DXX+E4iIhKZ0tOD+5kzvUXQn+pEJGa4RYsotJPFokXhjiIid90FJ50EtWv7TiIiEplatoS//c3r70kVgyISExat2ULpyjU5fP3KvVfWr7/3MhEJrYoV4YwzfKcQEYlcNWvCY495jaBuoiIS9Wb8voELBv3MUyddSW65cruvTE6G/v39BBOJV716wTvv+E4hIhL5tm6F+fO9HT6kxaCZnWFmM81sjpndUcj6zmaWbWYT8m93hzKPiMSecQvXcdGg0ZjBVc/1I/GFFyA1NRiVKzUVBg+GHj18xxSJH7/9BkOHwrJlvpOIiES+a6+Fzp29HT5k3UTNLBF4FjgVWAKMNbORzrlpe2z6g3OuS6hyiEjs+n7WKq57dRy1KpXhtauPpl615KDwU/En4s/AgVCuHFxxhe8kIiKRb9ek85s3Q/nyYT98KFsG2wNznHPznHM7gDeBriE8nojEkY8nLefqoWNJrZ7MO306BoWgiPiVnR18qbnkEqhSxXcaEZHIl5ER3M+e7eXwoSwGU4DFBZ4vyV+2p45mNtHMPjWz5iHMIyIx4o0xi7jpjfG0rluFt67rSK2KZX1HEhGAV18Npne5/nrfSUREooPn6SVCOZpoYSO8uz2ejwdSnXObzOwsYATQeK8dmfUGegPU16iAInFt4Ldz+c9nM+icXpOBPdpSrnSi70giskuDBsH1L5mZvpOIiESHxo2DcQ5mzPBy+FC2DC4B6hV4XhfY7Wpy59wG59ym/MefAElmVmPPHTnnBjvnMp1zmTVr1gxhZBGJVM45Hvp0Ov/5bAbntK7D4MszVQiKRJqzzw4GbRIRkaIpVy74vdnVz9V0oSwGxwKNzayBmZUGLgZGFtzAzA43M8t/3D4/z5oQZhKRKJSb5/jXe5N5/rt5XNahPk/9pQ2lS2lmHJGIMnIkrF7tO4WISPS55hpo08bLoUP2bco5lwPcBIwCpgNvO+emmlkfM+uTv9kFwBQzmwgMAC52zu3ZlVRE4tj2nFxufmM8b45dzM0nHcEDXVuQmFBYL3QR8WbFCrjgAvi///OdREQk+qxcCZ98Anl5YT90KK8Z3NX185M9lg0q8PgZ4JlQZhCR6LV5ew59XhvHD7NX8++zm3LN8Q19RxKRwrz0EuzcCddd5zuJiEj0efdduOEGWLQI6tU78PYlSP2sRCQird+ygx4v/spPc1bzyAWtVAiKRKrcXHj+eTjppD9HxRMRkaLbNb2EhxFFVQyKSMRZsWEbFz0/mmnLNvBcj7ZclBnev5KJSDF8+mnw12xNJyEicnA8Ti8R0m6iIiLFtWD1Zi576VfWbd7BK1e249gj9hpgWEQiyZgxUKeOt5HwRESiXu3aUKGCl+kl1DIoIhFj+vINXDBoNJu35/D6tR1UCIpEg/vvD77AJCX5TiIiEp3Mgq6iahkUkXiVtWAtVw4ZS/nSpXjj2o40Pqyi70giciBbtkByMlTU/1cRkUPy7LNQqVLYD6uWQRHx7tuZK7nspV+pUaEM/7tehaBIVNi+HY44Av7zH99JRESiX/v2fw4kE0YqBkXEq5ETl3HN0Cwa1qjA29d1pG7VZN+RRKQo3n0Xli/3NlGyiEhMWbkyGJl58eKwHlbFoIh489ovC7nlzd84qn5V3ryuAzUrlvEdSUSKauBAaNQITj3VdxIRkei3YgX06QM//xzWw+qaQREJO+ccz307l0dHzeSkjFo8e+lRlCud6DuWiBTV5Mnw44/wyCOQoL8ri4gcsiOOCAaSCfOIoioGRSSsnHP83yfTeeGH+XRtU4fHLmxNUqK+TIpElUGDoEwZuPJK30lERGJDuXKQmhr2EUVVDIpI2OTk5nHn+5N5O2sJV3RM5Z5zmpOQYL5jiUhx3XknnHwy1ND0LyIiJSYjQy2DIhKbtu3M5ZY3f2PU1BX89eTG3HZKY8xUCIpEpZQUOO883ylERGJLejr88AM4F3QZDQP1zRKRkNu0PYerhoxl1NQV3N2lGbef2kSFoEg0cg6uuw6++sp3EhGR2HPnncFoomH8jqRiUERCat3mHfR44Rd+nb+Wxy9szVXHNfAdSUQO1ujRMHgwzJ3rO4mISOypVQuqVg3rIVUMikjILM/eyoXPj2b67xsZdFlbzm9b13ckETkUAwdCxYpw6aW+k4iIxJ7t2+Guu2DUqLAdUtcMikhIzF+9mcte/JXsrTsZemV7Ojaq7juSiByK1avh7bfhmmugQgXfaUREYk/p0vD007B+PZx+elgOqWJQRErc1GXZXPHyGPIcvHFtB1rWrew7kogcqldegR074PrrfScREYlNZsEgMmGcXkLdREWkRI2Zv5aLn/+FpMQE3r6uowpBkYNkZi+b2Uozm+I7CwCHHQY9e0KLFr6TiIjErjBPL6FiUERKzNczVnD5S79Ss2IZ/nf9MRxRS13JRA7BEOAM3yH+0LMnDB3qO4WISGxLTw9GFN28OSyHUzEoIiXigwlL6T1sHI0Pq8DbfTqSUqWc70giUc059z2w1ncOIJhKYutW3ylERGJfRgaULw9LloTlcCoGRaT4hg+HtDRISIC0NH649ylufWsCbVOr8sa1HahRoYzvhCJSUhYvhtNOg4ce8p1ERCT2de8OGzcGLYRhoGJQRIpn+HDo3RsWLgwmoF64kLb97+DOteMZelV7KpZN8p1QJK6YWW8zyzKzrFWrVpX8AV54Ifi/ftVVJb9vERHZXWKiJp0XkQjWrx9s2bLbouSc7Vzz2YuUTUr0FEokfjnnBjvnMp1zmTVr1izZne/cGRSDZ54Z9AYQEZHQu+MOuOeesBxKxaCIFM+iRYUutsWLwxxEREJuxAj4/XdNJyEiEk4TJ8LIkWE5lIpBESky5xxbDq9T+Mr69cMbRiTGmdkbwGgg3cyWmNnVYQ/x1VeQmhq0DIqISHikp8OsWZCXF/JDqRgUkSJZv2UHN7/xG3e0vZhtSXsMEJOcDP37+wkmEqOcc5c452o755Kcc3Wdcy+FPcTAgfDrr8E1LCIiEh4ZGcElOUuXhvxQKgZF5IB+mL2K05/6ns+m/E767X1IeunFoLXALLgfPBh69PAdU0RKUk5O8H/8sMN8JxERiS+7RhINw+TzpUJ+BBGJWlt35PKfz2Yw5OcFHFGrAi9d0Y4WKZWBI+Dyy3zHE5FQ2bIl+DJy771wdfh7p4qIxLWMjOB38I4dIT+UikERKdTkJdnc+tZvzF21mV7HpHHHmRkaLVQkXrz5ZjDhcePGvpOIiMSf2rXD0ioIKgZFZA85uXkM/HYuT381m+oVSvPq1e05vnEJD1cvIpFt4EBo3hyOP953EhERCSEVgyLyh4VrNnPbWxMYv2g9XVrV5sFuLaiSXNp3LBEJp7FjISsL/vvfsE58LCIiBTz8MHzwAYweHdLDqBgUEZxzvDl2MQ98NI3EBOPpi9vQtU2K71gi4sPAgcEIwZdf7juJiEj8ys2FX36BzZuhfPmQHUbFoEicW7VxO3e8O4mvZqzkmEbVeezC1tSpUs53LBHx5fbb4ZRToHJl30lEROLXrhFFZ82CI48M2WFUDIrEsc+n/s6/3pvMxu053NWlGVcek0ZCgrqFicS1Fi2Cm4iI+LOrGJw5U8WgiJSsTdtzeODDabyVtZhmtSvxxsVtaHJYRd+xRMQn5+Bvf4PLLoOjjvKdRkQkvjVuHFy3HeJRRVUMisSZrAVrue3tCSxdt5UbOjfi1lOaULpUgu9YIuLbN9/Ak09CmzYqBkVEfCtbFnr0gNTUkB5GxaBInNiRk8dTX85i0HdzSalajreu60i7tGq+Y4lIpBg4EKpVg4su8p1EREQAXn015IdQMSgSB2av2Mitb01g6rIN/CWzHned04wKZfTfX0TyLV8OI0bALbcEf40WEZHIsHMnlCoVsql+1DdMJIbl5Tle/nE+Z//3R37P3sbgy9vynwtaqRAUkd29+CLk5MB11/lOIiIiu7zxRjDVz5IlITuEvhGKxKjl2Vvp+85EfpqzhpMzavHw+a2oWbGM71giEonKlYNLLgkGLBARkchQu3bwh7oZM6BevZAcQsWgSAz6YMJS7hoxhZw8x8PnteQv7ephIepeICIxoG9f3wlERGRPGRnB/cyZcOqpITmEikGRGJK9ZSf//mAKH05cxlH1q/DERW1Iq1HedywRiWRjxkBmJiToyhERkYhy2GFQqVJIp5fQb36RGPHj7NWc/tT3fDp5OX1Pa8Lb13VUISgi+zdnDhx9dDClhIiIRBazYPL5mTNDdgi1DIpEuW07c/nPZzN45acFNKpZnhd6HkvLupV9xxKRaPD885CYGFwvKCIikefaayEvL2S7VzEoEsWmLM3m1rcmMGflJnodk8Y/z8igXOlE37FEJBps2wavvALdukGdOr7TiIhIYa69NqS7VzEoEoVy8xyDvpvLk1/MonqF0gy7qj0nNKnpO5aIRJN33oE1a+D6630nERGRfXEOli0Lrh2sWLHEdx/SawbN7Awzm2lmc8zsjkLWm5kNyF8/ycyOCmUekViwaM0WLnp+NI+OmsnpLQ5n1K0nqBAUkeJ7/31o0gROOsl3EhER2ZepU6FuXfjoo5DsPmTFoJklAs8CZwLNgEvMrNkem50JNM6/9QYGHnDH48ZBWhoMH16ieUUi0vDhwec9IQGXmsroBwZw5tPfM2vFRp6+uA3PXHIkVZJL+04pItHo7bfhs8+CAQpERCQyHXFE8Hs6RIPIhLJlsD0wxzk3zzm3A3gT6LrHNl2BYS7wC1DFzGofcM8LF0Lv3ioIJbYNHx58zhcuBOewRYtoff8/uH7Zr4y69QS6tknR3IEicnCcg1KloEED30lERGR/ypYNfleHaHqJUBaDKcDiAs+X5C8r7jaF27IF+vU7lHwika1fv+BzXkByznZu/OJl6lQp5ymUiES9DRuCocrff993EhERKYoQTi8RymKwsCYLdxDbYGa9zSzLzLJ223DRInLz9tpcJOqt3bwDt2hRoets8eJCl4uIFMlrr8Hs2ZBStL+9ioiIZxkZQTEYgikmQlkMLgHqFXheF1h2ENvgnBvsnMt0zmUWXL60Yg1OeeI7Xv91Edt25pZQbBF/Fq/dwj0fTOGYh79iacUahW9Uv354Q4lI7HAOnnsOjjoK2rXznUZERIrikktg8GDILfl6J5TF4FigsZk1MLPSwMXAyD22GQn0zB9VtAOQ7ZxbXpSdu+RkVt1xDxXKlOLO9ydz3H++4dlv5pC9dWfJ/hQiYTBt2QZuefM3Oj/2La+PWcQ5reqQ+NBDkJy8+4bJydC/v5+QIhL9fvwxGJnu+us1cIyISLRo1w4uuwySkkp81yGbZ9A5l2NmNwGjgETgZefcVDPrk79+EPAJcBYwB9gCXFmknaemYv37c2SPHox0jtFz1zDwu7k8Omomz30zh0uPrs/VxzXk8MplQ/KziZQE5xyj561h0Hfz+H7WKsqXTuSqY9O46rgG1K5cDmgNlcsG1w4uWhS0CPbvDz16+I4uItFq4ECoXDn4K7OIiEQH5+DXX4Pf302bluiuzbnouuYuMzPTZWVlFbpuytJsBn8/j48mLSMxwejWJoXrOjXkiFolP0GjyMHKzXN8PvV3Bn03l4lLsqlRoQxXHpvGZUenUjm55P/iIxLNzGzcnpcIyL7t7xwJwLffwoIF0KtXmBKJiMghcw6qVAlaB599Fii582PIWgZ9aJFSmQGXHMnfT0/nhR/m8XbWYt4Zt4RTmh7G9Z0b0ja1mu+IEse27czlvfFLeeGHecxfvZm06sn0796C84+qS9mkRN/xRCQedO7sO4GIiBSXWchGFI2pYnCXetWSub9rC245uTFDRy9k2OgFnD9wBe3SqtKnUyNOTK9FQoKulZDwyN66k9d+WcgrPy1g9abttEypzLOXHsUZLQ4nUZ9DEQmH3Fy47z644gpo1Mh3GhERKa6MDPjmmxLfbUwWg7tUr1CG209twnUnNOStsYt58Yd5XD00iyaHVaD3CY04t3UdSpcK5Rg6Es9+z97Gyz/N5/VfF7Fpew7HN67B9Z3a0LFRdU0WLyLhNWoUPPAAtGypYlBEJBqlp8Orr8KmTVChQontNqaLwV3KlynFVcc14PKOqXw4cRnPfzePvu9M5PHPZ3L1cQ24pH19ypeJi7dCwmDOyo08/908RkxYSm6eo0urOvQ+oSEtUir7jiYi8WrgQDj8cOjWzXcSERE5GBkZwf2sWcH0QCUkriqgpMQEzjuqLt2PTOHbmasY+N1cHvx4Ov/9eg6Xd0il17Fp1KhQxndMiVLjFq5j0Hdz+WLaCsomJXBJ+/pce3xD6lVLPvCLRURCZcEC+PjjYGTiEAxLLiIiYdCpUzAIWHp6ie42rorBXcyMEzNqcWJGLcYvWsfz383l2W/n8MIP87gwsy7XHt+Q1OrlfceUKJCX5/hm5koGfTeXsQvWUSU5ib+e3JgrOqZSXX9YEJFIMHhwMPhA796+k4iIyMGqUSMoCEtYXBaDBR1VvyrPX57J3FWbGPzdPN4eu4TXf13EmS1rc32nRuraJ4XakZPHyInLGPz9XGat2ERKlXLc3aUZf2lXT12ORSSy7NwJF14I9er5TiIiIofi889h82bo3r3EdhlT8wyWhBUb8gf9+GURG7fncNwRNbiuU0OOO6KGBv0QNm/P4Y0xi3jpx/ksz95G+mEV6dO5IV1a1SEpUYMRiZQ0zTNYPPs8RzoXtA6KiEj0OvtsWLYMfvtN8wyGymGVyvKvM5ty44lHMPyXRbz803wuf2kMLVIqcd0JjTizxeGU0pf+uLN603aG/LSAV39ZSPbWnbRvUI3/696Szuk19UcCEYlcM2YEgw7o95SISPRLTw+ml8jLK7Fdqhjch0plk7i+cyOuOi6N98cvZfD387j5jd+oXy2Za49vwIWZ9TRReBxYtGYLg3+YyztZS9iRm8dpzQ7juk6NOKp+Vd/RRET2b+pUaNECXn4ZrrzSdxoRETlU6emwdSssWVJiu1QT1wGUKZXIxe3r88XtnRh0WVuqlS/NXR9M5diHv+a/X81m8ytDIS0NEhKC++HDfUeWgzF8+G7/joufeZGbXh9P58e+4e2xS+jWJoUvb+/E85dnqhAUkegwaBCULg1duvhOIiIiJWHX9BIzZ5bYLtUyWESJCcYZLQ7n9OaH8ev8tQz6bi6zn3qehFHPwM7twUYLF/45WluPHv7CSvEMHx78u23ZEjxfuJDqt91E+XNu4dpre3HVsQ04rFJZvxlFRIpj0yYYNgwuughq1vSdRkRESsKuaSVmzCixXWoAmUOwo259Si9dvNfyzYensHTCdI6oWYGEBF2nEam27MhhwuL1tDimNZVWLttrfV69+iQsWughmYjsogFkiuePc+TgwXDddfDTT3DMMb5jiYhISXAO5s2DtDSsVCkNIONb6WWF99ct9/syTnvyeyqXS6JtalUy06qSmVqNVnUr6zpDj1Zu3Ma4BesYu2Ad4xauZeqyDeTkOeatXF7o9glL9i70RUSiwvDh0KoVdOzoO4mIiJQUM2jUqER3qWLwUNSvH3QN3UNu3bo8ckErxi1YR9bCtXw9YyUASYlGy5TKZKZVIzO1Km1Tq2pi8hDJy3PMWbWJrPx/g6wF61i0NugGWqZUAq3rVaH3CQ1pl1YNXq8HixftvZP69cOcWkSkhHzyCSxapFFERURizahR8PXXJbY7FYOHon//3a81A0hOJunhh7gosx4XZQYT/K7dvINxC9eRtWAtWQvXMeSnBQz+fh4ADWuU/6PlsG1aVRrWKK+pCg7Ctp25TFqS/UfhN27hOrK37gSgevnSZKZV5fIOqbRNq0qLOpUpXarA2EkP/V+h/4707x/mn0JEpISULw9Nm/pOISIiJW3cOHjkkRLbnYrBQ7FrkJh+/YK/wNavHxQQewweU618aU5tdhinNjsMCAqXyUuz84uWtXw+bQVvZwVdTquXL81RqVVpl1aVtqnVaJFSiTKl1LV0T2s2bQ8K7Pwie/LSbHbmBte/NqpZnjOaH07btKq0S6tGWvXk/RfYRfx3FBGJeDk50LYtPP44dO7sO42IiJS0XYPIlBANIBMB8vIc81ZvYuyCdX8UiAvWBK1UpUsl0KZulfzCpipt61ejcnKS58Th5Zxj/urNu3X5nLd6MwClExNoWbfyn62rqVWpVr6058QiUlI0gEzxZNar57KWLIFJk6BlS99xRESkpOXPIWtQIudHFYMRauXGbYxfGAx2krVwHVOXZpOTF/xbNa5V4Y/rDtulVaNetXJ7t3wNHx6ZLV1FyLU9J5cpSzcwbuFaxi5Yx/iF61izeQcAVZKT8q+3rEZmWlVapmhQHpFYFu/FoJmdATwNJAIvOuce3t/2mWXLuqzMTPjxx7DkExGRMNu+HZKTsbw8jSYay2pVLMsZLWpzRovaAGzdkcuExesZtzC47vCjSct4Y0ww6EnNimXITK36R4HY/NuPKNXnut3mzYuI+Q8Lmc+P3r3ZsiOXXzqeHozyuWAdE5asZ0dOHgBp1ZPpnF6LdmnBqKwNa2i6DhGJD2aWCDwLnAosAcaa2Ujn3LR9vmj7drj++jAlFBGRsCtTJhhRdPbsEtmdWgajVF6eY9bKjUHXyfyBaZas2wrAT4OuIiV75V6v2Vq7Lj9+6e+9O+6UTMot33s6jiWVanLc9a9QKsFonlKZdvnTcbRNrUbNihptVSSexXPLoJl1BO51zp2e//xfAM65h/b1msykJJe1aVPwZUFERGLTzp1Y6dJqGYxnCQlGxuGVyDi8Epd1SAXg9+xtZC1cS51HVhX6mjLLl3LtMH/F4LzlSwtdnrJxNW/27kDrulUoV1pdPkVE8qUABSc8XQIcvedGZtYb6A3QpFo1FYIiIrEuqeTGD1ExGEMOr1yWLq3q7HP+w5yUunx083EekuUf/7W6lF6690TuVr8+HRpW95BIRCSiFdYnfq/uPM65wcBgCHrPhDqUiIjEDhWDsWgf8x+W/s9DtEip7C/Xfx7SfH4iIkW3BKhX4HldYJmnLCIiEoMSDryJRJ0ePWDwYEhNBbPgfvBg/6OJRmouEZHINBZobGYNzKw0cDEw0nMmERGJIWoZjFU9ekRmkRWpuUREIoxzLsfMbgJGEUwt8bJzbqrnWCIiEkNUDIqIiEQo59wnwCe+c4iISGxSN1EREREREZE4pGJQREREREQkDqkYFBERERERiUMqBkVEREREROKQikEREREREZE4pGJQREREREQkDqkYFBERERERiUPmnPOdoVjMbCMw03eOQ1ADWO07xEGK5uwQ3fmjOTsov0/RnB0g3TlX0XeIaBHl58ho/6xGc/5ozg7K71M0Z4fozl8i58donHR+pnMu03eIg2VmWdGaP5qzQ3Tnj+bsoPw+RXN2CPL7zhBlovYcGQuf1WjNH83ZQfl9iubsEN35S+r8qG6iIiIiIiIicUjFoIiIiIiISByKxmJwsO8Ahyia80dzdoju/NGcHZTfp2jODtGfP9yi+f2K5uwQ3fmjOTsov0/RnB2iO3+JZI+6AWRERERERETk0EVjy6CIiIiIiIgcoogsBs3sQjObamZ5ZrbPEX7M7Awzm2lmc8zsjgLLq5nZF2Y2O/++aniSF+3YZpZuZhMK3DaY2a356+41s6UF1p0VruxFzZ+/3QIzm5yfMau4rw+VIr7/9czsGzObnv85u6XAurC///v6HBdYb2Y2IH/9JDM7qqivDbUiZO+Rn3mSmf1sZq0LrCv0MxRORcjf2cyyC3we7i7qa8OhCPn/XiD7FDPLNbNq+eu8vv9m9rKZrTSzKftYH7Gfe58sis+PRT2+6RzpLbvp/FiiipBf58gQKUJ2nR93cc5F3A1oCqQD3wKZ+9gmEZgLNARKAxOBZvnrHgHuyH98B/CfMGYv1rHzf47fgdT85/cCfT2+90XKDywAahzqz+8jP1AbOCr/cUVgVoHPTljf//19jgtscxbwKWBAB+DXor42ArIfA1TNf3zmruz7+wxFWP7OwEcH89pIyL/H9ucAX0fQ+38CcBQwZR/rI/Jz7/tGFJ8fD+b46BwZ1uzo/Bju/DpHesq+x/ZxfX6MyJZB59x059yBJs1tD8xxzs1zzu0A3gS65q/rCgzNfzwU6BaSoIUr7rFPBuY65xaGMlQxHOp75/O9L9LxnXPLnXPj8x9vBKYDKeEKuIf9fY536QoMc4FfgCpmVruIr/Wa3Tn3s3NuXf7TX4C6Ycx3IIfy/vl+7w8mwyXAG2FJVgTOue+BtfvZJFI/915F+fnxYI6vc2TJ0fkxvHSOjODvJ3uI6/NjRBaDRZQCLC7wfAl//sI6zDm3HIJfbECtMOYq7rEvZu8P4E35zb4vh7sLD0XP74DPzWycmfU+iNeHSrGOb2ZpwJHArwUWh/P939/n+EDbFOW1oVTc419N8JesXfb1GQqXoubvaGYTzexTM2tezNeGUpEzmFkycAbwboHFvt//A4nUz300iNTz48EcX+fIkqPzY3jpHHng14aKzo+7L9+vUiUarRjM7Evg8EJW9XPOfVCUXRSyLCxDo+4vezH3Uxo4F/hXgcUDgQcIfpYHgMeBqw4u6T6PWxL5j3XOLTOzWsAXZjYj/y8ZIVeC738Fgv/8tzrnNuQvDvn7v2eMQpbt+Tne1zbe/g/kK/LxzexEghPdcQUWe/sM7YpVyLI9848n6J62Kf/6mBFA4yK+NtSKk+Ec4CfnXMG/NPp+/w8kUj/3IRfN50fQOTKfl/9fOj/+sU0k/J7QOXL/rw0lnR93X75f3opB59wph7iLJUC9As/rAsvyH68ws9rOueX5zaYrD/FYu9lfdjMrzrHPBMY751YU2Pcfj83sBeCjkshcUEnkd84ty79faWbvEzRNf0+I3/uSym9mSQQnuuHOufcK7Dvk7/8e9vc5PtA2pYvw2lAqSnbMrBXwInCmc27NruX7+QyFywHzF/gShHPuEzN7zsxqFOW1YVCcDHu1rkTA+38gkfq5D7loPj+CzpH5+/ByjtT5MaJ+T+gcuZ/XhpjOj7sv369o7iY6FmhsZg3y/3p4MTAyf91I4Ir8x1cARflLakkpzrH36qOc/wt6l+5AoSMJhdAB85tZeTOruOsxcBp/5vT53hfp+GZmwEvAdOfcE3usC/f7v7/P8S4jgZ4W6ABk53fxKcprvWY3s/rAe8DlzrlZBZbv7zMULkXJf3j+5wUza0/wO3NNUV4bBkXKYGaVgU4U+L8QIe//gUTq5z4aROr5sbjH1zmyZOn8GF46R0bw95P8zDo/QsSOJtqdoOrdDqwARuUvrwN8UmC7swhGuppL0H1m1/LqwFfA7Pz7amHMXuixC8meTPAfpvIer38VmAxMyv8HrB3m9/6A+QlGKZqYf5saKe99MfIfR9BsPgmYkH87y9f7X9jnGOgD9Ml/bMCz+esnU2AEwX39Hwjj+32g7C8C6wq8z1kH+gxFWP6b8vNNJLi4/5hIee+Lkj//eS/gzT1e5/39J/iSvxzYSfD7/upo+dz7vBHF58f9Hb+Q/DpH+smu82N48+sc6Sl7/vNe6PyI5b9QRERERERE4kg0dxMVERERERGRg6RiUEREREREJA6pGBQREREREYlDKgZFRERERETikIpBERERERGROKRiUEREREREJA6pGBQREREREYlDKgZFIpiZtTOzSWZW1szKm9lUM2vhO5eIiIhPOj+KlAxNOi8S4czsQaAsUA5Y4px7yHMkERER73R+FDl0KgZFIpyZlQbGAtuAY5xzuZ4jiYiIeKfzo8ihUzdRkchXDagAVCT4C6iIiIjo/ChyyNQyKBLhzGwk8CbQAKjtnLvJcyQRERHvdH4UOXSlfAcQkX0zs55AjnPudTNLBH42s5Occ1/7ziYiIuKLzo8iJUMtgyIiIiIiInFI1wyKiIiIiIjEIRWDIiIiIiIicUjFoIiIiIiISBxSMSgiIiIiIhKHVAyKiIiIiIjEIRWDIiIiIiIicUjFoIiIiIiISBxSMSgiIiIiIhKH/h9JkHt8OtKFPQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(15,5))\n",
    "ax1 = fig.add_subplot(121)\n",
    "ax2 = fig.add_subplot(122)\n",
    "\n",
    "ax1.plot(x,ue(x),label='exact solution')\n",
    "ax1.plot(x,ua,'r o',label='numerical solution')\n",
    "ax1.legend(loc='upper left')\n",
    "ax1.set_xlim([-1,1])\n",
    "ax1.set_ylabel('u(x)')\n",
    "ax1.set_xlabel('x')\n",
    "\n",
    "ax2.plot(x,np.abs(ua-ue(x)),'r--')\n",
    "ax2.set_xlim([-1,1])\n",
    "ax2.set_ylabel('|ua(x)-ue(x)|')\n",
    "ax2.set_xlabel('x')\n",
    "ax2.set_title('Error')\n",
    "plt.show()"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}