FEM-exercise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "## Introduction to Finite Element Methods - Exercise Class\n",
    "\n",
    "### DST Module 4 - November 2018\n",
    "\n",
    "### James Kermode\n",
    "\n",
    "<b>j.r.kermode@warwick.ac.uk</b>\n",
    "\n",
    "Warwick Centre for Predictive Modelling / School of Engineering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup code\n",
    "\n",
    "Startup code for examples and visualisations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 345,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "using QuadGK\n",
    "using Plots\n",
    "using LaTeXStrings\n",
    "pyplot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 356,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 356,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L = 10.0\n",
    "E = 100.0\n",
    "A = 1.0\n",
    "n_cells = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Main solver routine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 357,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "solver (generic function with 1 method)"
      ]
     },
     "execution_count": 357,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function solver(L, f, n_cells; quad_degree=3)\n",
    "    n_nodes = n_cells + 1\n",
    "    mesh = linspace(0.0, L, n_nodes)\n",
    "    h = L/n_cells\n",
    "    k_e = (E*A/h)*[ 1 -1; -1 1]\n",
    "    K = zeros(n_nodes, n_nodes)\n",
    "    for element in 1:n_cells\n",
    "        K[element:element + 1, element:element + 1] += k_e\n",
    "    end\n",
    "    quad_points, quad_weights = gauss(Float64, quad_degree)\n",
    "    b = zeros(n_nodes)\n",
    "    for element in 1:n_cells\n",
    "        x_mid = (mesh[element + 1] + mesh[element])/2.0    \n",
    "        for (zeta, weight) in zip(quad_points, quad_weights)    \n",
    "            x = x_mid + zeta*h/2.0\n",
    "            w = weight*(h/2.0)\n",
    "            b[element] += w*(0.5 - zeta/2.0)*f(x) \n",
    "            b[element + 1] += w*(0.5 + zeta/2.0)*f(x)\n",
    "        end\n",
    "    end\n",
    "    K[1, :], K[:, 1], K[1, 1] = 0.0, 0.0, 1.0\n",
    "    b[1] = 0.0\n",
    "    return K \\ b, mesh\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualisation of solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 358,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAYAAAByNR6YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd4VFX+BvB3eq/phST0JiUgSJNmQV0URRRWXRu2tbvrD9e1BMW1LcvaVte67mJdBRu6dkVUmnTpHZKQhCTTMr2c3x8xIyFlJiSTTJL38zw8mjvn3vlm6ptzzz1HIoQQICIiIqI2I+3oAoiIiIi6Ggasbio3Nxd9+vRp1TH27NkDiUSCa665po2qSqwJEyZALpd3dBndVnu9XpLxef70008xduxYmM1mSCQSzJo1q6NLqqctPg/aw5dffgmJRIKHHnqoo0tJ6Os5GV/D1HIMWF3EgQMHIJFImv1nt9vbpRZ+OCSHZPoyakv33nsvJBIJvv/++44uJS579+7F+eefj4MHD2Lu3LkoKirCxRdf3K41dLbHLFl0luBJyYnfgl1M7969cdlllzV6m1qtjv7/8uXLIZFIWnVf+fn52L59O8xmc6uOQ9SW3njjDXi93o4uI+qLL76A3+/HE0880e7BihKHn38UCwNWF9OnTx/Mnz8/ZrvevXu3+r4UCgUGDBjQ6uMQtaW8vLyOLqGe0tJSAEB2dnYHV0JtiZ9/FAtPEXZTjXV9H3sa4bXXXsOwYcOgVquRlZWFO+64Az6fr17748cghEIhSCQS/PDDDwiHw/VOTx4/TmHjxo2YPXs2MjMzoVQqkZ+fj9tuuw3V1dUt+j3Ky8tx++23o3fv3lCpVEhLS8NFF12Ebdu2teg4S5cuxdSpU2E2m6FWqzFkyBAsWrQI4XC4XruXXnoJEokEr732Gj744AOMHj0aWq0WOTk5KCoqQiQSAQC88sorGDp0KDQaDfLz8/G3v/2t0fsVQuCll17C2LFjYTAYoNVqMWrUKPz73/9u0LYlz8+9996LM844AwBw33331XsuiouLm30sIpEIXnjhBYwaNQpWqxUajQY9evTAeeedh++++65B+5dffhmjR4+GTqeDXq/HmDFj8J///Kf5B/wYzZ2GOf5084QJE/CXv/wFAHDqqadGf6dj92/qFHUwGMTChQujz4vJZMLUqVPx8ccfN2h77PP8v//9D2PHjoVWq0Vqaiquuuoq2Gy2mL9X3ftjwYIFDeo99lTd5s2bcdFFFyE9PR0qlQq9evXCH/7wh0bfC3WPlc1mw4033ojc3FzIZDK89tprTdYRz2NWx+Vy4ZZbbkFWVhbUajWGDRuGpUuXNnpcv9+PhQsXorCwEDqdDgaDARMnTmz08Ywl3vdfU7766itcddVV6N+/P/R6PfR6PUaPHo2XXnqp0fY//fQTLrzwQvTo0SP6uTF69Gg88sgjAH597kpKSrB3795675+6U+7NjcFyOp0oKirCkCFDoNVqYTKZUFhYiPvvvx+hUCjux8Xr9WLevHnIy8uDWq3GoEGD8OyzzzZoZ7PZ8Oijj2LixInIysqCUqlETk4OrrjiCuzfv79B+2M/S15++WUUFhZCq9Xi9NNPj7s2io09WNTAE088gc8++wwzZszAaaedhk8++QRPPPEEqqurG/3iryOVSlFUVIRXXnkFxcXFuP/++6O3jRgxIvr/7733HubMmQO5XI4ZM2YgNzcXW7duxVNPPYXPPvsMq1evhslkilnn7t27MXnyZBw5cgTTpk3DzJkzUVZWhiVLluDTTz/FN998g5NPPjnmcf7v//4PCxcuRG5uLmbNmgWj0Yjly5fjj3/8I9auXYs333yzwT7vvPMOPv/8c1xwwQUYP348li1bhgcffBASiQQqlQqPPfYYZsyYgcmTJ2PJkiW48847kZWVhUsuuSR6DCEE5syZg//+97/o168fLrvsMigUCnz++ee48sorsX37djz66KMN7jue52fq1Kk4dOgQFi9ejClTpmDixInR/Y1GY8zHY9GiRejTpw8uvfRS6HQ6lJaWYsWKFfj666/rHevGG2/Ec889hx49euDaa69FJBLBu+++iyuuuAKbNm1qMlieqKuvvhpSqRQrVqzAVVddFe2tslqtze4XiUQwc+ZMLFu2DP3798fNN98Ml8uFt99+G9OnT8dTTz2FW265pcF+S5cuxccff4zzzjsP48ePx7fffotXX30V+/btw/Lly5u9T6vViqKiInz99dcN6q377/Lly3H22WcjGAzioosuQn5+Pn744Qf8/e9/x7Jly7Bq1aoGv5vP58PkyZPh8/kwY8YMyGQypKent/ox8/v9OP300+FyuTBr1iy43W689dZbmDVrFr744gucdtpp9Wo488wzsWLFChQWFmLu3LkIBAL4+OOPMX36dDz33HO44YYbmn186pzI++94jzzyCA4ePIgxY8YgJycHdrsd//vf/3Dttddi9+7deOyxx6Jt161bh/Hjx0OhUOC8885DQUEBbDYbtm3bhpdeegl333139LlbtGgR5HI5br311uj+x77+G1NeXo6JEydi165dKCwsxO9//3tEIhFs374djzzyCObNmwe9Xh/XYzNr1ixs2bIFF154ISKRCJYsWYKbbroJBw8erPc7bd26FUVFRZg6dSpmzpwJnU6Hbdu24bXXXsMnn3yC9evXo0ePHo0+bt9++y3OO+88TJs2DUqlMq66KE6CuoT9+/cLAKJ3796iqKiowb+VK1fWa5+TkyN69+5db9s999wjAAiLxSJ27doV3e52u0Xv3r2FVCoVZWVl0e27d+8WAMTcuXPrHWf8+PFCJpM1Wmd5ebnQ6/WiR48e4tChQ/VuW7x4sQAgbr/99rh+59GjRwu5XC6++OKLetu3b98udDqdKCwsjFnXJ598IgCIc845R7jd7uj2cDgsrr32WgFAvP/++9HtL774ogAgVCqVWLduXXS7w+EQqampQqfTiaysLLF///7obfv37xcKhaJBPc8++6wAIK655hoRDAaj2/1+vzj77LOFRCIRGzZsiG5v6fPzxRdfCABiwYIFzT6OxzMajaJHjx7C4/HU2x6JRERVVVX056+++koAEIMHDxYOhyO6vaqqSvTu3VsAED/88EN0e1Ovl8Zei3Uae87qHocVK1bEvc/LL78sAIjTTjtNBAKB6Pb9+/cLq9UqFAqFOHDgQHR73fOsUCjqvXdCoZCYMGGCACDWrl3b6P0fr6l6Q6GQKCgoEBKJRHz55Zf1brvjjjsEAHHdddfV256TkxN9vXq93rjuv7kajj/uzJkzhd/vj27/9NNPBQDxm9/8pl77efPmCQDigQceEJFIJLrd4XCIwsJCoVar670Wm9LS919Tr+l9+/Y1OHYgEBBTp04VcrlcFBcXR7ffeuutAoBYtmxZg30qKyvr/dzca7Op1/OMGTMEAHH//fc32Ke0tFSEQqFGj3es8ePHCwBi0KBB9d5bNptN9O3bV0gkErF+/fp626urqxsc5/PPPxdSqVTccMMN9bbXvR4MBoP4+eefY9ZDJ4anCLuYvXv34oEHHmjwb9WqVXEf44477kDfvn2jP2u1WsyZMweRSATr169vVX3//ve/UVNTg8cee6zBX1SXXXYZhg4dirfeeivmcdauXYs1a9bg6quvbtCtPWDAAMydOxcbNmzAjh07mj3OM888AwB48cUXodVqo9ulUmm096ixv6Avv/zyer1yRqMR55xzDtxuN26++WYUFBREbysoKMDYsWPx888/R08h1t23wWDAM888U++UllKpxF/+8hcIIRp9LBL5/Bxbw/Gn2SQSSb1ej7resgcffLBer5jVao32XjbX49me6up4/PHHoVAootsLCgpw2223IRgM4o033miw3+9+9zuMGTMm+rNMJsMVV1wBoPY12BrfffcdDhw4gOnTp9frHQKA+fPnw2w24/XXX2/0lNJf//rXehettJUnnniiXi/GtGnTkJOTU+93DYVC+Oc//4n+/ftHTz/XMRqNuO++++Dz+fDee+/FvL8Tff8dr2fPng22KRQKXH/99QiFQvj2228b3K7RaBpsS0lJiXlfzSkpKcGHH36Ifv364b777mtwe1ZWFmQyWdzHu//+++u9t8xmM/785z9DCFHvNLzZbIbFYmmw/xlnnIH+/fvjyy+/bPT4N9xwAwYPHhx3PdQyPEXYxUybNg2ffvppq44xcuTIBttyc3MBoNVTPdQFvZUrV2Lnzp0Nbg8EAigrK4Pdbm/26py64xw5cqTRQf27du0CAOzYsaPZgairVq2CwWDACy+80OjtarW60ZA2fPjwBtuysrKavS0YDOLo0aPIyMiA0+nEtm3b0KNHj+i4j2P5/f5o/cdL5PMDAHPmzMELL7yAk046CbNnz8aUKVMwZsyYBl9IGzZsAABMnjy5wTGmTJkCoHasXTLYsGEDDAZDvVBcZ8qUKSgqKmq01kQ+1s09fkajESNGjMDXX3+N3bt3Y+DAgdHbdDodBg0a1Kr7bkxaWlqjp5Fyc3OjtQLA9u3b4XQ6kZ+fjwceeKBB+/LycgCNv3aPd6Lvv+M5nU4sXLgQ77//Pvbt2we3213v9roLDQDg4osvxtNPP41zzz0Xc+bMwRlnnIGJEye2yUUIa9euhRACU6dObZOpak499dQmtx37nADA119/jSeffBJr1qxBZWVlvWB+bHg91ujRo1tdIzWNAYsaaGyMTt2HRbyDTptSN3D36aefbrad2+1uNmDVHeejjz7CRx991OxxmmOz2SCEaPSLorljNPcYNXdbMBiM3i8AHD58uE3vu7XPD1Dbq9C7d2+8+uqrWLBgARYsWACNRoPZs2dj4cKF0b/ynU4n5HJ5o+Of6sKm0+lsdT2tJYSAy+VqciB9c7Um8rGuu7+MjIwW1dVU+9ZqamyeXC6v97vWvfe2bNmCLVu2NHm8WO894MTff8fy+/2YOHEiNm3ahBEjRuDyyy+H1WqFXC7Hvn37sHjx4ugfLAAwfvx4fPPNN3j44Yfx2muv4ZVXXgFQGzYef/xxTJo0KWbdTXE4HACAnJycEz7GsRobW1f3/NfdFwC89dZbuOSSS2AwGDBt2jQUFBRAo9FAIpHglVdeqRcwGzsWJQYDFrWrug/x7du3t+oS57rjtGQwbWMMBgM0Gg3KyspO+Bgnoq7+U045pUWnb9uDQqHAvHnzMG/ePJSWluLbb7/FK6+8gldffRUVFRXRq8SMRiMOHjyI6urqBiGr7vGMNaAeqD0d1NSVVcd+iZwoiUQCg8EQ7Vk5XktqbUt199fSulo7f11r1dUze/bsuE7nN6ct3n9Lly7Fpk2bcP311+Of//xnvdtee+01LF68uME+kyZNwqRJk+DxeLB69Wp89NFHePbZZ3HOOedg69at9U7xt0TdH4UlJSUntP/xKioqGvSs1b1ejr0QqKioCFqtFuvXr28wBU9zV5h29Gupq+MYLGpzMpkMQgiIRtYRP+WUUwDUniJsjbY8Tnl5eaOXMieSxWJB3759sXXr1oT18tSN9WhNT0t2djYuueQSfPbZZ+jZsyc+++yzaG9AYWEhADQ6vqVuW2OnS49nsVhQVlbWoE6Xy4U9e/Y0aH8iv1dhYSFcLlejY9RaUmtbau7xc7lcWLduHXQ6Xb3xdieqLV4LdQYPHgydToe1a9e2aMqBxrTF+2/v3r0AgBkzZjS4bcWKFc3uq9VqMWXKFCxatAh33XUXPB4Pvvrqq+jtMpmsRY/ZqFGjIJFI8PXXX7f6sQEar79uW93rB6h9DAYPHtwgXBUXF+PAgQOtroNODAMWtTmr1YpIJNLoX3FXX301dDod7r77bmzfvr3B7XV/UcYybtw4nHzyyVi8eDHefffdBrdHIpGYl9EDiF5+fdVVVzU679CRI0carbMt3HrrraipqcF1110Hj8fT4PZ9+/bh4MGDJ3z8ul6lw4cPx72P1+tttEfN7XbD7XZDoVBEv6zrBnvPnz8fNTU10bY2mw0PPvhgvTbNGTVqFPx+f73ekEgkgrvuuqvB3Gsn+nvV1fGnP/2p3hffoUOH8MQTT0ChUNSbQqM9TJw4EQUFBVi2bFmDkPXggw/Cbrfj0ksvbZOxPCfymDVFqVTihhtuwL59+3DXXXc1GiS2bNmCysrKmMdqi/dffn4+ADRYBujrr7+Onv471o8//ljvlGGdup6hYy8esFqtqKioaLR9Y7Kzs3H++edj165d0fnPjtXYHxLNefDBB+v9AWa32/Hwww9DIpHg8ssvj27Pz8/Hrl27cPTo0eg2r9eLG264oU1CNZ0YniKkNjd16lS8//77uOCCCzBt2jSo1WoUFhbiN7/5DTIzM/HGG29g9uzZGDp0KM466yz0798fPp8PBw8exLfffotJkyZh2bJlMe/nrbfewpQpU3DRRRdh7NixGDFiBNRqNQ4ePIiVK1fCbrfX++JvzPTp03H33XfjkUceQZ8+fXDWWWchLy8PVVVV2L17N77//ns8+uij9QYZt5WbbroJq1atwuuvv44VK1bgtNNOQ3Z2NsrLy7F9+3asXr0a77zzTvQLpKUGDRqEjIwMvP7665DL5cjJyYFEIsFtt90Gg8HQ6D5utxtjx45F//79MWLECOTl5cHlcmHZsmWoqKjAn/70p+gX/tSpU/H73/8ezz33HAYPHoyZM2dG58EqLS3FH/7wB4wbNy5mnbfccgv+85//4KqrrsKnn36K1NRUfPfdd6ipqcFJJ53U4At2ypQpkEgkuOuuu7Bp0yaYTCZYrVbceOONTd7HlVdeGZ3TaujQoZg+fXp0HiybzYYnn3zyhB/nEyWTyfDqq6/i7LPPxplnnomLL74YeXl5+OGHH/Ddd9+hb9++ePjhh9vkvk7kMWvOQw89hA0bNmDRokX48MMPMXHiRKSlpaGkpASbN2/G5s2bsXbtWqSmpjZ7nLZ4/82YMQN5eXl45JFHsHnzZgwcOBA7d+7EsmXLcMEFF2DJkiX12j/88MP4/vvvceqpp6Jnz55Qq9X46aef8M0336BPnz71esKmTp2KjRs34pxzzsGECROgUCgwefJkTJgwocl6/vnPf2Lbtm148MEHsWzZMkydOhXhcBg7d+7EF198gerq6rjmwZLJZMjPz8dJJ52ECy+8EEIIvPvuuygpKcG8efPq9WDdcsstuOOOOzB8+HBceOGFCAaD+PzzzyGXyzFkyJAWT7xMbaTDJoigNlU3D9a0adPiat/cPFiNzZVTNy/Q4sWLo9uamgcmEAiIO++8U+Tl5Qm5XN5om23btomrr75a5OXlCaVSKaxWqxg6dKi47bbbxE8//RTvry2qqqrEPffcIwYPHiw0Go3Q6/Wib9++4tJLL603f44Qzc/P9dlnn4np06eLtLQ0oVAoRGZmphg3bpx46KGHxOHDh5t9HOo09/hdeumlAkC9Y9V58803xWmnnSYsFotQKBQiJydHTJkyRSxatKjevDwtfX6EEOLHH38Up556qtDr9QJAkzXU8fv94tFHHxVnnHGGyMnJEUqlUmRkZIhJkyaJt956q0H7SCQiXnzxRXHyyScLjUYjtFqtGD16tHj11VcbtG3q9SKEEF9++aUYPXq0UKlUIiUlRVxxxRWioqKiyefs5ZdfFieddJJQqVTR+d/qNLVPIBAQjz/+eHQ/g8EgJk+e3Oh8SM09zy2dXyzWHFQbN24UF154oUhJSREKhUIUFBSI22+/vcGcTEI0Py9TLM09Zi2di0wIIYLBoHjuuefEuHHjhNFoFCqVSuTn54uzzz5bPP/88/XmtYol3vdfU4/9nj17xMyZM0VaWlr0NfjOO+802v6TTz4Rv/vd70S/fv2EXq8Xer1eDBo0SNx7770NHnOn0ynmzp0rMjMzhVQqrXes5l7Pdrtd3HPPPWLAgAFCpVIJs9ksCgsLRVFRUdzzYMlkMuHxeMSdd94pcnNzhVKpFAMGDBDPPPNMg/aRSET84x//EIMGDRJqtVpkZWWJa6+9Vhw9evSE5pKjtiERopGBMkRERER0wjgGi4iIiKiNMWARERERtTEGLCIiIqI2xoBFRERE1MYYsIiIiIjaGAMWERERURtLioDl8Xiwfv36RmezJiIiIupskiJg7dixAyNHjsSOHTs6uhQiIiKiVkuKgEVERETUlTBgEREREbUxBiwiIiKiNsaARURERNTGGLCIiIiI2pi8owsgIiIiSjQhBPzhSO2/0K//3b5zB775/DPIpRIM6t8PV1w0ExKJpNX3x4BFREREnVIo8mtY8tULTuGGYSocabD/vq2b8cDVc+D3eiGTy+HzuLF+cxGeXFDU6pDFgEVERERJQQiBwPFh6ZjA5Ku3LYKwECd8X36fFw/O/S2y8nvivpffhFZvwNtPL8TTf3kApxQOx6UXnt+q34UBi4iIiBImFBHwh8P1TsvV9jg13NZYL1OieGtq4LLbcNPDf4fOYAQAzL7lTnzw8nPYvGMnLm3l8RmwiIiIKG61vUy/hKZfwtGvPU7hBmOcQq3oZWoPTltV9P9rHHZEImEopK2/BpABi4iIiJoUEQI2XxBHPQE8/fe/Ydl/34QAkJGbh5sfeQIGs6WjSzwhppRUDJ8wGS8/dC+EEDCYLXj32Seg0elx4YzprT6+RIiOj5br16/HyJEjsW7dOowYMaKjyyEiIuq2IkLA5g3iqDeASo8fVd4gwkLg3eeewJtPPo4JvzkfRmsKVny0FGnZuSj619vQm8wdXfYJ8fu8ePzmudj4/bcAAHNKKt5d9glOGzOq1cdmDxYREVE39mug8qPSE0CVN4DwcV0v3poavPXUX3H+NTfid3feCwA4bdZvMe/Cs/DdR0txzmVXd0DlLSOVACqZFCqZDCq5FCqZFGqrDm+//yE8Dhs0cjny061Qq9Vtcn8MWERERN1IOFJ7yq/S68dRTwDVjQSq44VCQQgh0G/YyOi2gv6DoNJoEQoEElxx05RSyS9hSVYbnn4JTnXb1Mdsk0slTU+9kGJo89oYsIiIiLqwukB11ONHpTe+QHU8hVIFmVyONV99hlGnTYNUKsWmH5bD566BWqdvs1ob62VSyaRQR4NT/SAlbYMJQROFAYuIiKgLqQ1UARz1BFDprT3lF2nlaGu1VosbH/obnrn7dpTu3wuDxYItK7/H8FOnYPL5s5rdV3FML5O6QS/Tr2FKHauXqZNhwCIiIurEEhGoGjP5/IugUmvw7QfvQCqRYPoll+OmPxfBqNPWC0zqul6mX7Ylcy9TIjFgERERdSLhiEC1L4BKT22oqvYlJlAdSymVIFWrxPWXX4J7r78SRpW8y/Q0JQoDFhERURLrkEAlkyBVo0KaVolUrRJGJQNVSzFgERERJZG6QHXUUxuqGKg6JwYsIiKiDtQxgUqKVI2SgSqBGLCIiIjaUTgiUOWtHZB+1BOArZ0CVZqmNkylaZUwMFAlHAMWERFRAoUiAtXtHKhUMilStcpoLxUDVftjwCIiImpjgXAE//34U7z5+hvwhUJIy+mBmdfdArlCkZD7Y6BKPgxYREREbUCI2lN/BxxefLBsGR69aS7Sc3vAYDJj+YdLcGj3Ttyx8B+QyVv/1VsXqOpO+zFQJR8GLCIiolbwh8I46PTigMODmkAYwYAff73tOgyfMAl3PvkiFEol1nz1Kf56yzX4asmpOHP2ZS2+D5VMWjsgnYGq02DAIiIiaiEhBCo8ARywe1Ba48OxQ6pCgSD8Xi8mnjsTCqUSADD6tLOgN5lR47DFdfxooNIqkaZRQa+UMVB1MgxYREREcfIEwzjo8OCgwwtPKNxoG6lMCgA4tHsnxv+yraL4MHweD2Syxr926wUqrQp6BQNVZ8eARURE1IyIEChz+3HA7kGZ2x+zvUqjxfnX3Ih3n3sCNQ4bzKnp+OrdN2HNzMLEc2fWtmGg6vIYsIiIiBrhDoRw4JfeKl840qJ9L/vjPZArlPjq3TchIJCR0wN/ffFVnNSnF1K1SgaqbkAihEjwbByxrV+/HiNHjsS6deswYsSIji6HiIi6qXBE4EiND/sdHhz1BFp1LKkEyNarUWDSIk2rZKDqZtiDRURE3Z7TH8QBhxeHnF4EWthbdTyDUo4CkwZ5Ri1UcmkbVUidDQMWERF1S6GIQImrdnqFKm+wVceSSYAcgwY9zVpY1Qr2VhEDFhERdS92XxAHHB4cdnoRbOWaNWaVHAVmLXoYNFDI2FtFv2LAIiKiLi8YjuCwy4sDdi/s/tb1VsmlEvQwalBg0sKiTszSN9T5MWAREVGXJIRAtS+IA3YPil0+hFt5TZdVrUBPsxY5BjXkUvZWUfMYsIiIqEvxhyM47KgdW+UMhFp1LKVMgjyjFgUmDYwq9lZR/BIesM4880yUlZVBKpXCYDDgqaeeQmFhYaLvloiIuhEhBI56AjjgqF26ppVDq5CmVaLApEW2Xg2ZlAPWqeUSHrD++9//wmw2AwDee+89XHnlldi0aVOi75aIiLoBXyiMg7/0VrmDjS9dEy+VTIp8U+3YKr2SJ3iodRL+CqoLVwDgcDh46SoREbWKEALlbj/2Ozwoq/GjtbNlZ+pUKDBpkalXQcrvKGoj7RLRL7/8cnzzzTcAgE8++aTJdjU1NXA6ndGfVSoVVCpVwusjIqLk5wmGf1m6xgNvqHWTgWrkMhSYNMg3aaFVyNqoQqJftetSOf/+97/x9ttvNwhZdUvlHK+oqAjz589vp+qIiCjZRETt0jUHHF6Ux7HQcnMkALL0avQ0a5HOpWsowdp9LUKNRoPi4mKkpKREt9UFrOXLl2P48OHR7ezBIiLqnmqOWWjZ38qla/QKGQrMWuQZNVDL2VtF7SOhpwjtdjs8Hg+ys7MBAO+//z5SUlJgtVobba/X62E0GhNZEhERJalwRKC0xof9dg8qva1faDnHoEGBSYNUDXurqP0lNGA5HA5cdNFF8Hq9kEqlSEtLw7Jly/hCJyKiqEA4goXPvYjX/vUyIkIgNSsb189/DDqjqcXHMirl6GnWoodRAyWXrqEOlNCAlZ+fjzVr1iTyLoiIqJMKRiLYU+3GM8/9E8/ePw/DJ0yGNT0Da776HA9cPQdFr7wVV8iSSSToYVRHl67hH/GUDDiPgnRCAAAgAElEQVTRBxERtatwRGCf3Y2d1W54fH68sOAeTL1wDm586G+QSCQ4e9sW/HnOefjinddx/twbmzyORa1AgUmLXKMaCi5dQ0mGAYuIiNpFRAgcsHuwo7oGvl+mWYhEwggFgzjplPHRnqdeg4ZAbzIj4PM1OIbimIWWzVxomZIYAxYRESWUEAKHnF5sr6qB57jZ1qVSGeQKJTb/uAKnTr8AUqkUe7ZsQo3DDpVaE22XqlGiwKRBjkHDpWuoU2DAIiKihBBCoKTGh+2VLrgCjS9jI1cocP38R/HsvX9EVVkpLOkZWPftl+g58CRM/+3l6GvRocCshYFL11Anw1csERG1KSEEytx+bKt0weEPxWw/9cI5UCiV+OKdN1B1pBTjTjsTTz71FPrnZHDpGuq0GLCIiKjNHPX4sfWoC9W+YIv2O/XcmTjvotkYnGpAGmdZpy6AAYuIiFqt2hvA1koXjnpaPkGoSSXHoFQDMnUqBivqMhiwiIjohDl8QWytdKHsBNYJ1CtkGJhqQK5BzWBFXQ4DFhERtZgrEML2SheKXQ2nUohFI5dhYKoeeUYNx1hRl8WARUREcfMEQ9heWYODTm+L91XJpBiQokeBScupFqjLY8AiIqKYvKEwdlbVYL/dA9HCfZVSCfpZ9ehl0UHOYEXdBAMWERE1yR+OYFdVDfbZ3Qi3MFnJpRL0sejQ16KDggsvUzfDgEVERA0EwxHstrmxx+ZGKNKyZCWTAL3MOvRL0UPFYEXdFAMWERFFhX5ZiHlXVQ0CLQxWEgA9zVr0T9FDI5clpkCiToIBi4iIEI4IHHB4sKOqBv5wpMX75xs1GJCqh07BrxUigAGLiKhbi/yyEPOOyhp4Qo2vF9icXIMaA1MNXCuQ6Dh8RxARdUNCCJS4fNhW6UJNsOXBKlOnwqBUA8xqRQKqI+r8GLCIiLqRli7EfLw0rRKDUw2wapQJqI6o62DAIiLqJircfmytdMHWwoWYAcCqVmBwmgFpWlUCKiPqehiwiIi6uCpvANu4EDNRu2LAIiLqouy+ILad4ELMBmXtQsw5ei7ETHQiGLCIiLoYlz+EbVUulJzAQsxahQwDU/TowYWYiVqFAYuIqItwB0LYXlWDQyewELNaLsUAqx4FZi2DFVEbYMAiIurkvKEwdlTV4MCJLMQsk6C/VY9eZh1kXIiZqM0wYBERdVL+UAS7qmuw1+5GC1e1gVwqQV+LDn2sOiikXC+QqK0xYBERdTIutwfvf/M9ip0eGFLSkF3QK+59ZRIJelu06GfVQ8mFmIkShgGLiKiTEELg50NHMOu86di1eQMAQCqV4pbHnsLEc2c2u69UAvQ01S7ErOZCzEQJx4BFRNQJ1ARC2FBmx/UXXoDSg/tx/ytvIS0rB++9+AyevutWGC1WDJ8wucF+EgD5Jg0GpBigVTBYEbUX9g8TESWxiBDYWVWDLw8cRVmND9vXrcFvb78Lw8ZNRHbP3vj9Q3+DMSUV29etabBvD4MaZ/RMw4hMM8MVUTtjDxYRUZKq8gawocwBZ6D+moFSya9/G0skEkhQ/+q/LH3tQswmFRdiJuooDFhEREkmEI5g61EX9js89bZLpVL0HVqIN596HPn9ByI1OwdLn38atqPl6DdsBNK1SgziQsxESYEBi4goSQghUOLyYVOFE/5wpMHtEokE9zy/GPOvuhh//u150e13ProId/zuYqRyIWaipMGARUSUBNzBEDaWO1EeY91Ag8WKh17/ADvWr4VCCowd0Bunjx3F9QKJkgwDFhFRB4oIgT02N7ZX1iAs4pstVKPTYea552BwqoFzWRElKQYsIqIOUu0NYEO5Aw5/KHbjXxiVchRmmpDCcVZESY0Bi4ionQXDEWyrdGGv3RO78S+kEmBgigF9rTouxkzUCTBgERG1o1KXDxsrHPCFGg5ib0q6VonhGSbolfzIJuos+G4lImoHnmAYmyocOFLT/CD2Y6lkUgxJN6KHQc1B7ESdDAMWEVECCSGw1+7BtqMuhOIcxA7ULm8zJM3IQexEnRQDFhFRgth9Qawvd8DuC8a9j14pw4gME+e0IurkGLCIiNpYKBLB9soa7LG5EW+flVQC9Lfq0c+qh0zK04FEnR0DFhFRGyqr8WFjuROeUDjufVI1ShRmmmDgIHaiLoPvZiKiNuANhbG5wokSly/ufZRSCYakG5Fn1HAQO1EXw4BFRNQKQgjsd3iw9agLwUj8g9jzjBoMSTNAJZclsDoi6igMWEREJ8jhD2JDmQPVLRjErlPIUJhhQrqOg9iJujIGLCKiFgpHBLZXubC7Ov5B7BIA/ax6DEjhIHai7oABi4ioBcrdfmwsd8AdjH8Qu1WtwIhME4wqRQIrI6JkwoBFRBQHXyiMLRVOHG7BIHaFVILBaQb0NGk5iJ2om2HAIiJqhhACB51ebKlwtmgQe65BjaHpRqg5iJ2oW2LAIiJqgssfwoZyByq9gbj30cplGJ5hRKZencDKiCjZMWARER0nHBHYWV2DXdU1iLfTSgKgj0WHgal6yKVcP5Cou2PAIiI6xlGPHxvKHKhpwSB2i1qBwgwTzGoOYieiWgxYREQA/OEIfq5w4qDTG/c+cqkEg1MN6GXmIHYiqo8Bi4i6NSEEDju92HzUhUA4Evd+2XoVhqaboFVwEDsRNcSARUTdVk0ghI3lDlR44h/ErpFLMSzdhGwDB7ETUdMYsIio24kIgV3VNdhRFf8gdgDobdFiUKoBCg5iJ6IYGLCIqFup8gSwodwBZyAU9z4mlRwjMk2wqJUJrIyIuhIGLCLqFgLhCH4+6sIBhyfufWQSCQal6tHbooOUg9iJqAUYsIioSxNCoMTlw6YKJ/wtGMSeqVNheIYRWgU/Jomo5TiQgIi6LG8whEt/fyt6ppoxY2Aurp9yMg7s2NrsPmqZFKdkmzE2x8JwRUQnjAGLiLqkIzU+XHbjbXjz+Wcw/YprcPU9C2C0pGD+lRejeO/uRvfpZdbijJ5pyDFoOK8VEbUKAxYRdSnhiMDGcgdWltjwwasv4IJrb8bv7rwXZ/32ChT9621IJBIs/+DdevsYlXJMykvB8AwTFDJ+LBJR67H/m4i6DIc/iLWl9ugVgiISQXpuXvR2vckMrdGEiKgdiyWVAANTDOhr5SB2Impb/FONiDo9IQT22tz45mBlvekXzKnp+OrdN+B2OgAAKz/7GBXFh2BJTUe6VokzCtLQP0XPcEVEbY49WETUqflDYawrc6DM7W9w293PvYoHrp6Dm8+aAL3RhLLDBzHhrHNx3513oMCi5zgrIkoYBiwi6rTK3X78dMTe5PQLvQYPxYLFS7H8w3cRiURQkJuDR/48D1oVJwwlosRiwCKiTiccEdhW6cJumztm27x+A3Dl/92HYelG5Jt4dSARtQ8GLCLqVFz+ENYcscHhj2+pG7NKjlHZFhiU/LgjovbDTxwi6hSEEDjg8GJzhRNhEd8KzX0tOgxOM3AQOxG1OwYsIkp6gXAE68scKK3xxdVeLZNiZJYZGTpVgisjImocAxYRJbWjHj/WHrHDF4pvHcFMnQojM01QyWUJroyIqGkMWESUlCJCYHulCzurYw9kB2onDR2SZkQvs5YD2YmowzFgEVHSqQmEsPaIHTZfMK72RqUco7PNMKoUCa6MiCg+DFhElDSEEDjk9GJTuROhOAey9zZrcVKaETIpe62IKHkwYBFRUgiGI9hQ7kCxK76B7CqZFCMzTcjUqxNcGRFRyzFgEVGHq/IGsPaIHZ5gOK726VoVTs4yQc2B7ESUpBiwiKjDRITAzqoabK+qiau9VAKclGpEbwsHshNRcmPAIqIO4QnWDmSv8sY3kN2glGFUlgVmNQeyE1HyY8AionZ32OnFxnIHgpH4BrL3NGkxJN0IOQeyE1EnwYBFRO0mGIlgU7kTh5zeuNorpRKMyDQj28CB7ETUuTBgEVG7qP5lILs7zoHsaVolTs40Q6PgQHYi6nwYsIgooYQQ2FXtxrZKF+I5ISgBMCjVgH5WHQeyE1GnxYBFRAnjCYaxrsyOo55AXO11ChlGZZlh1SgTXBkRUWIxYBFRQpS4vFhfFv9A9nyjBsMyjJBLpQmujIgo8RiwiKhNhSIRbK5w4YDDE1d7hVSCwgwTco2aBFdGRNR+GLCIqM3YfUGsPWKDKxDfQPYUjQKjsszQKvhRRERdS0I/1Xw+H+bMmYNt27ZBo9EgPT0dzz33HPr06ZPIuyWidiaEwB6bG1srXYjnjKAEwIAUPfqn6CHlQHYi6oISPtjhuuuuw86dO7Fp0ybMmDED11xzTaLvkojakS8Uxg/F1dhyNL5wpVXIMDEvBQNTDQxXRNRlJTRgqdVqnHPOOdFLrceMGYMDBw4k8i6JqB0dqfHhqwOVqIjzKsEeBjVOy09FCq8SJKIurl0HPjz55JOYMWNGk7fX1NTA6XRGf1apVFCpVO1RGhG1QDgi8PNRJ/ba4xvILpdIMDzDiDyTNsGVERElh3YLWA8//DD27NmDr776qsk2kyZNqvdzUVER5s+fn+DKiKglHP4g1pba4QyE4mpvUdcOZNcrOZCdiLqPdvnEW7hwIZYuXYovv/wSWm3Tf8EuX74cw4cPj/7M3iui5CGEwD67B1uOOuMaawUA/a16DEzlQHYi6n4SHrAWLVqEN998E19++SXMZnOzbfV6PYxGY6JLIqIW8ofCWFfmQJnbH1d7jVyKk7PMSNPyjyQi6p4SGrCKi4vxxz/+Eb169cKUKVMA1PZKrV69OpF3S0RtqNztx09H7PCHI3G1z9arMSLTBKWMM7ITUfeV0ICVm5sLIeI8l0BESSUcEdhW6cJumzuu9jKJBMPSjcg3abhIMxF1exx1SkQNuAIhrCm1weGPbyC7WSXHqGwLDBzITkQEgAGLiI4hhMABhxebK5wIx9n73Neiw+A0ThpKRHQsDpIgIgBAIBTGDXfdi5PysvDbEX1x07TxOLRrR5Pt1TIpxudaMSTdyHBFRHQcBiwigtMfxLX/dw9e+OvDmHjuTFx88x+h1mhQdMUsFO/b3aB9pk6F0wpSkaHjVYJERI1hwCLq5opdXnxzsApvP/8UzvndXFx7/8M476rrMf9f/wUkEny95K1oW6kEGJZuxNgcC1RyWQdWTUSU3DgGi6ibiojaqwR3VddeJRgJR5Bd0Ct6u8FihdFiReSX6RmMSjlGZ5thVCk6pF4ios6EPVhE3ZA/FMEPxdXRcAUARqsVyz9cAk+NCwCwYcU3KDt0EEarFb3NWkzJT2W4IiKKE3uwiLoZmy+IVSU2eEPhetvnPf0KHrh6Nm47ZxKMViuK9+7BiFMn4y/3/An5Vq6wQETUEgxYRN3IQYcHG8odja4l2GfIMCxYvARfLXkLIhLB1N+ch0UP3g+TVtP+hRIRdXIMWETdQEQIbKpwYr/d02y7ggGDMfeeBehl1mIop18gIjphDFhEXZw3FMbqEhuqfcGYbaUSoDDDhHyTth0qIyLquhiwiLqwSk8Aq0ttcS3UrFXIMCbbArOaA9mJiFqLAYuoCxJCYJ/dg80VTsSz4E26VolR2RaoZLywmIioLTBgEXUxoYjAxnIHDjm9cbXvZ9VhcKoBEo63IiJqMwxYRF2IOxDCqlIbHP5QzLZyiQQjs0zIMfAqQSKitsaARdRFlLl9WFtqR7CxORiOo1fWjrfixKFERInBgEXUyQkhsLPajW2VrrjaZ+lVODnTDAXHWxERJQwDFlEnFgxH8FOZHUdq/HG1H5xqQD+rjuOtiIgSjAGLqJNy+oNYVWpDTSAcs61CKsHobAsydKp2qIyIiBiwiDqhEpcX6444EBKxx1uZVHKMybZAp+TbnYiovfATl6gTEUJga6ULu6rdcbXvYdSgMMMEuZSnBImI2hMDFlEn4Q9HsLbUhgpPIGZbCYCh6Ub0Mms53oqIqAMwYBF1AjZfEKtLbPCEYo+3UsmkOCXbglStsh0qIyKixjBgESW5gw4PNpQ7EMf0VrCqFTglxwKNXJb4woiIqEkMWERJKiIENlc4sc/uiat9T7MWw9KNkPKUIBFRh2s0YH333Xcxd5w4cWKbF0NEtbyhMNaU2lDlDcZsK5UAhRkm5Ju07VAZERHFo9GANXnyZEgkEohjLgE/fqBsOBx7LAgRtVyVJ4BVpTb4w5GYbTVyGcbkWGBRc8kbIqJk0mjA2rBhQ4NtNpsNn332GZYsWYLnn38+4YURdTdCCOyze7C5wok4hlshTavE6CwLVHIueUNElGwaDVjDhg1rtPHkyZOh1Wrx/PPPY8qUKQktjKg7CUcENpQ7cMjpjat9P6sOg1INHG9FRJSkWvyn77hx4/DJJ58kohaibskdDOHbQ5VxhSuZRIJTss04KY2D2YmIklmLryJ8//33YbVaE1ELUbdT7vZjbakNgTjmYNArasdbGVUcb0VElOwaDVjnnXdeg22BQAA7d+7EoUOH8Pjjjye8MKKuTAiBXdVubK10xdU+S6/CyZlmKGQcb0VE1Bk0GrCcTmeDqwbVajVOP/10zJo1C9OmTWuX4oi6omAkgnVHHCit8cXVflCqHv2tei55Q0TUiTQasL799tt2LoOoe3D5Q1hVWg1XIPY0JwqpBKOyzMjUq9uhMiIiakucyZ2onZS6fPipzI5QHOOtTCo5Tsm2QK/kW5SIqDPipzdRggkhsK3ShZ3V7rja9zCoUZhpglzK8VZERJ0VAxZRAgXCEawptaPC44/ZVgJgSLoRvc1ajrciIurkGLCIEsTuC2JVqQ2eYOzxViqZFKOzzUjTqtqhMiIiSjQGLKIEOOTwYEO5A+E41ryxqBU4JdsCrUKW+MKIiKhdMGARtaGIENhS4cReuyeu9j1NWgxNN0Im5SlBIqKuhAGLqI34QmGsLrWjyhuI2VYqAYanm1Bg1rZDZURE1N4YsIjaQJU3gNUlNvjCkZhtNXIpTsm2wKpRtkNlRETUERiwiFpBCIH9dg82VTgRx3ArpGmVGJ1lhkrO8VZERF0ZAxbRCQpHBL7Ztg+HHG6oNFpodLpm2/e16DA4zQApp2AgIuryOJMh0Qlw+QL4zexLccaQPpg7YRiuHj8E67/7utG2MokEo7PMGJJuZLgiIuomGLCIWsjmC2D21dfii6Vv48o/zcfdz/0bQ8eeisduuhpb166q11ankGFyfgpyjZoOqpaIiDoCTxEStUCJy4ufjtjxzQdLcP61N+PcK68DAAwbPwm3nDUBKz/9CINHjQEAZOpUODnLDKWMf8cQEXU3DFhEcRBCYFe1G1srXdFtBrMl+v8KpRIavSH688AUPQak6LnkDRFRN8U/rYliCEcE1pU56oWrrIKe+GTxy6goPgwA+GrJmzi8ewdye/bC2BwLBqYaGK6IiLox9mARNcMfimBVaTWqvMF62//0j3+h6IqLcNO0cZArFAj4fJh+6ZVYeM88GNWc34qIqLtjwCJqgtMfxMoSG9yNLNacmpWDh15/H6s+/xiRcBj52Vm46/oroZLzLUVERAxYRI0qd/uxptSGYKTp6UMtaek4+9Kr0Meiw5A0nhIkIqJfMWARHWefzR3XzOwSAMMzTOjJ9QSJiOg4DFhEv4gIgS0VTuy1e2K2VUglOCXbgnSdqh0qIyKizoYBiwhAMBzBmiN2lLv9MdvqFDKMy7XCoOTbh4iIGsdvCOr23MEQVhbb4AyEYrZN1ShxSo4FKk4eSkREzWDAom6tyhvAqhIb/OFIzLb5Rg0KM01cT5CIiGJiwKJu67DTi3VldjRzoWDUSakG9LXqeKUgERHFhQGLuh0hBLZX1WBHVU3MtjKJBKOyzMg2qNuhMiIi6ioYsKhbqV32xo5ily9mW7VcirE5VljUinaojIiIuhIGLOo2fKEwVpbYYPMFY7Y1qxQYm2OBRiFrh8qIiKirYcCibsHhC+LHEhu8oYbL3hwvW6/GyVlmyKUcb0VERCeGAYu6vCM1PqwttSMkYo9m72fVYXAql70hIqLWYcCiLksIgb12DzZXOGO2lQAYkWlCvonL3hARUesxYFGXFBECm8qd2O+IveyNUibBmGwrUrXKdqiMiIi6AwYs6nIC4QjWlNpQ4QnEbGtQyjA2xwo9l70hIqI2xG8V6lJqAiGsLKmGKxB7MHu6VonR2RYouewNERG1MQYs6jIqPX6sKrUhEI49mL2nSYthGUYue0NERAnBgEVdwkGHB+vLHIhj1RsMTTeit1nLKwWJiChhGLCoUxNCYGulC7uq3THbyiUSjMo2I0vPZW+IiCixGLCo0wpFIvjpiAOlNbGXvdHKZRiba4FJxWVviIgo8RiwqFPyBsNYWVINuz8Us61VrcCYHAvUci57Q0RE7YMBizodmy+IlSXV8IUiMdvmGtQYmWmGjMveEBFRO2LAok6l1OXD2iN2hONY9mZgih4DUvQczE5ERO2OAYs6BSEEdle78XOlK2ZbqQQYmWlGD6OmHSojIiJqiAGLkl5ECGwoc+Cg0xuzrUomxdgcC6waLntDREQdhwGLkpo/HMHqEhsqvbGXvTEq5RiXa4FWwZc1ERF1LH4TUdJyBUL4sbga7mDsZW8ydCqMzjJDwWVviIgoCTBgUVKqcPuxutSGYCT2YPY+Fh2GpBk4mJ2IiJIGAxYlnf12DzaWx172RgJgWIYRvcy69iiLiIgobgxYlDSEENhy1IU9ttjL3iikEpySbUG6TtUOlREREbUMAxYlhWAkgrWldpS5/THb6hQyjMuxwqDiy5eIiJITv6Gow3l+WfbGEceyNykaJcZkW6CSczA7ERElLwYs6lDV3gBWltjgD8de9ibPqEFhhonL3hARUdJjwKIOU+z04qcyO+K4UBCDUw3oZ9XxSkEiIuoUGLCo3QkhsKOqBturamK2lUmAk7PMyDFw2RsiIuo8GLCoXYUjAuvL7Djs8sVsq5ZLMTbHCota0Q6VERERtZ2EjxS+9dZbUVBQAIlEgo0bNyb67iiJ+UJhrDhcFVe4MqnkmJKXynBFRESdUsJ7sGbNmoV58+ZhwoQJib4rSmIHyyrwzOvvwBsIwmixYuTk05scT5WlV2FUlhlyKa8UJCKizinhAWvixImJvgtKctsOluCcM8/AwV3bo9umX3EtrvzT/AYhq59Vh8GpXPaGiIg6t6Qag1VTUwOn0xn9WaVSQaXiTN2d2UGHB7NmXQh71VH8/cOvkdOrDz5/ezFeWnAPMnrk45zLrgZQu+xNYaYJBSZtxxZMRETUBpLqHMykSZNgMpmi/x555JGOLolOkBAC2ytdWFfmwOHdOzH9imuR128AZHI5zr70KvQcdBIO794JAFBKJZjQw8pwRUREXUZS9WAtX74cw4cPj/7M3qvOKSIENpY7cMDhjW7zun9dXzAcDsPv8QAA9MraZW/0yqR6KRIREbVKUn2r6fV6GI3Gji6DWiEUiWB1qR3lx6wpOPr0s/DBy88iLTsXOb364Mt3XkfZoQO4bf5fMDkvFUpZUnWkEhERtVrCA9b111+Pjz/+GGVlZZg2bRoMBgP27NmT6LulDuALhfFjcTXsx60peMMDj8Pv9eL5onkAAKlMhgXPPI87L58NKQezExFRFyQRQsSxUElirV+/HiNHjsS6deswYsSIji6HToDTH8SPxTZ4QuFGbxdCoPzwQUTCYRQWZGNc/568UpCIiLqspDpFSJ1TpSeAlSXVCDazqKBEIkFWXgFGZpqQx8HsRETUxTFgUasUu7z46UjsBZvlUgnGZFuQruOFC0RE1PUxYNEJEUJgj82NLUddMduq5VKMz7HCxGVviIiom2DAohYTQmBzhRN77Z6YbY1KOcblWqFVyNqhMiIiouTAgEUtEo4IrD1iQ2mNP2bbNK0SY7ItUHAaBiIi6mYYsChu/lAEK0uqUe0Lxmzbw6DGyCwzp2EgIqJuiQGL4lITCOHH4mrUBBufhuFYXLCZiIi6OwYsiqnaG8DKEhv84UjMtsMzjOhl1rVDVURERMmLAYuadaTGhzWlNoRjTMMgk0gwOtuMLL26fQojIiJKYgxY1KR9djc2ljtjtlPJpBibY4FVo2yHqoiIiJIfAxY1IITA1koXdlW7Y7bVK2QYn2uFTsmXEhERUR1+K1I9ESGw7ogdh12+mG2tagXG5lqh4jQMRERE9TBgUVQgHMHqUhuOegIx22br1RiVZYZMyisFiYiIjseARQAATzCMH4ur4QyEYrbtbdFiaJqR0zAQERE1gQGLYPcF8WNJNXyh2NMwDEkzoK9V3w5VERERdV4MWN1chduPVaU2hCLNz8MglQAnZ5mRa9C0U2VERESdFwNWN3bQ4cH6MgdiTHEFhVSCsTlWpGo5DQMREVE8GLC6ISEEdlbXYFtlTcy2WoUM43OsMKj4UiEiIooXvzW7mYgQ2FjuwAGHN2Zbs0qOcblWqOWydqiMiIio62DA6kZCkQhWl9pR7vbHbJuhU+GUbDPkUs5xRURE1FIMWN2EL1Q7DYPdH3sahgKTBsMzTJByGgYiIqITwoDVDbj8IfxQUg1PMByz7aBUPfpb9ZzjioiIqBUYsLq4Sk8AK0uqEYwxDYMEwIhME/JN2vYpjIiIqAtjwOrCil1e/HTEjhjZCnKpBKdkW5ChU7VPYURERF0cA1YXtbu6BluOumK2U8ulGJdjhVmtaIeqiIiIugcGrC5GCIHNR53Ya/PEbGtU1k7DoFVwGgYiIqK2xIDVhYQjAmuP2FFa44vZNk2rxCnZFihlnIaBiIiorTFgdRH+UAQrS6pR7QvGbNvDoMaITDNkUl4pSERElAgMWF2AOxDCD8XVqIljGoZ+Vh0Gpxo4DQMREVECMWB1ctXeAFaW2OAPR2K2HZ5hRC+zrh2qIiIi6t4YsDqxIzU+rCm1Iyyan4dBJpFgdLYZWXp1O1VGRETUvTFgdVL77G5sLHfGbKeSSTE2xwKrRtkOVRERERHAgNXpCCGwtdKFXdXumG31ChnG5VqhV/JpJluJsq4AAA9+SURBVCIiak/85u1EIkJg3RE7DrtiT8NgVSswNscKlZzTMBAREbU3BqxOIhiOYFWpDUc9gZhts/UqjMqycBoGIiKiDsKA1Ql4gmH8WFwNZyAUs21vsxZD042choGIiKgDMWAlOYcviB9KquELxZ6GYUiaAX0sOoYrIiKiDsaAlcQq3H6sKrUhFGl+GgapBDg504xco6adKiMiIqLmMGAlqUMOD9aVOdB8tAIUUgnG5liRquU0DERERMmCASvJCCHwytIP8dSivyEcCkNvMuOaex+CNSOzQVutXIZxuRYYVYoOqJSIiIiawmv4k4gQAv94cyl+f8lseFwuWNMzsHvTetx/+SxUl5fVa2tWyTEpP4XhioiIKAkxYCWJcERgVakNC+bdgYEnj8ZDr7+H2xf+AwteWwq3046lLzwdbZuhU2FiXgo0clkHVkxERERN4SnCJBAIR7CypBpV3iD8Pi8Gjx4HhVIFAMjMK0BmXk/4fV4AQL5Jg8IME6S8UpCIiChpsQerg3mCYSw/VIUqbxAAoFSpsf2n1QgFa3+uKClG2eEDUKrUGJiixwiGKyIioqTHHqwO5PQH8X1x/Tmu5t77EBbedi2KrpiF7J69sXnlCqi1Otx3958wMNXQgdUSERFRvNiD1UEqPQEsP1TVYALRUVPPxLynX4ZEIkHxnt3I7zsAn3zxJcad1L+DKiUiIqKWYg9WByhx+bD2iA1NzR86cvLpGDn5dKhkUozL/f/27jU26nrP4/hn6LTTljK9Uhla2m4L5Sa0FMqWshF8gGfDIpDF3bBGkYUV2cSoTzbug9UYo7AJYuJqVEywXtAe0DaLJx7WoOeI7Cmr1MJWWVigtnQoxQK90dJOL/PbB2C32MIMdC6dmfcr4UFnvg7f2D8/PvnP8GmykmPpuAIAIJQQsALsp/ZuHf+50+PcxOgoLc1MUUIM3yIAAEINf3sHiDFGJ6906dSVLo+zSbZolWYmK5YaBgAAQhIBKwDcxujYzx0619HjcTY9PkYlGcmyTuDjcQAAhCoClp8NuI2+u9Cmi90uj7PT7HFaOIUaBgAAQh0By49cg24dOd+q1t5+j7P5KRM1N22SLIQrAABCHgHLT7r7B/Sn863q6hv0ODt/sl3TUyYGYCsAABAIBCw/aO/tV9X5VvUOum87Z5G0yJGkafa4wCwGAAACgoDlY5euuXSkqU0Dtyq5usE6waKSqclKn2gL0GYAACBQCFg+dL6zR9UX229ZIPoLW9QELc1MUVJsdGAWAwAAAUXA8pG6tm79d4vnAtGEGwWiEykQBQAgbPG3/BgZY3Ti8lWdbu32OJscG63SjGTZKBAFACCsEbDGwG2Mai52qLHTc4HoPRNt+vOpSRSIAgAQAQhYd2nA7da3F9r1sxcFoln2OBVRIAoAQMQgYN0F18Cgqpra1OZFgejMlATNSUugQBQAgAhCwLpD3X0D+s/zreru91wgWpBuV14yBaIAAEQaAtYdaLtRIOryUCA6wSIVO5KUMYkCUQAAIhEBy0st3S79V1ObBsztS66iJ1hUkpGsyfEUiAIAEKkIWF5wdvaourldHvpDFWu9XiCaaKNAFACASEbA8uBMa5d+uHTV49ykmCgtzUxVfDQdVwAARDoC1i0YY/Tjpas60+a5QDQlNlpLMlNki6LjCgAAELBG5TZG3ze3y3m11+OsI8GmYkeyrBOoYQAAANcRsH6l3+3Wt01tarnW53E2JzFehffYKRAFAAA3IWAN0zswqKrzrWp3DXicnZ2aoFmpFIgCAICRCFg3dN0oEL3mRYHognsS9WdJ8QHYCgAAhCIClqTWnj5VNbWpz4sC0cWOZE2dFBugzQAAQCiK+IB1sbtX3za1a9CLAtHSjBSlxscEaDMAABCqIjpgneu4ppqLHR4LRONuFIjaKRAFAABeiMiAZYzR6dZunbjsuUDUHmNVaWYKBaIAAMBrERewjDGqbelUXfs1j7OpcTFakpGsGApEAQDAHYiogDXoNqq+2K4mLwpEpybEqtiRpCgKRAEAwB2KmIDVP+jWkaY2Xe7xXCCamxSvgnQ7HVcAAOCuRETA6rlRINrhRYHonLQEzUyhQBQAANy9sA9YV/sG9Cdnq64N3L5A1CJpwZRE5SRSIAoAAMYmrAPW9QLRVvUN3r6IIcoiLZ6aLEcCBaIAAGDswjZgNXf16rsLbfKQrRQTdb1ANCWOAlEAAOAbYRmwGjqu6ZgXBaLx1igtnZaiSTFh+b8BAAAESVglC2OM/re1S/9zucvjbKLteoFonJUCUQAA4FthE7CMMTre0ql6LwpEJ8fHqGRqsqIpEAUAAH4QFgFr0G10tLldF7o8F4hmTIrVoikUiAIAAP/x+y2cM2fOqLS0VPn5+SouLtaJEyd8+vqugUGVf/Odvj1eqys/N992Ni8pXotpZwcAAH7m94D1xBNPaMuWLTp9+rSeffZZbdy40Wev3dvXr1V/83d69P4SPbNqubYsW6iD+z4adfbetEmaTzs7AAAIAL8GrJaWFlVXV+uRRx6RJK1bt05Op1Nnz571yetv2PKP+sNnFdrywr/q5Y/364H1G/T28/+kqv/43dCMRdKiKYnKT6WdHQAABIZfP4PldDrlcDhktV7/bSwWi7KystTY2Kjp06ePmO/q6lJnZ+fQ1zabTTab7Zav/4cvfq+Vj/6DfrN+gyRp5oJFOlNbo+o/HlTpXz6oKItFJRnJumfirV8DAADA18bVP6NbtmyZEhMTh35t3779lrPGXG+5ip34/z/axmKxyBYbJ0myRU3QfVkphCsAABBwfr2DNW3aNDU3N2tgYEBWq1XGGDU2NiorK2vU+UOHDqmwsHDo69vdvbJYLMqfM1cH9pRp0fIVyp45W4f2V+hUzVH9xQMrtSwrVQkUiAIAgCDwawJJT09XUVGR9uzZo40bN6qiokKZmZmjvj0oSQkJCbLb7V6//r7f/lYrVjygf/7bvxp6bM2jf69/e/FfFBdNuAIAAMHh9xSya9cubdy4Udu2bZPdbldZWZnPXjtzcqqOHPqj9lT+u3p6++SYnKqH/3qNJkwYV+98AgCACGMxv3yYKYhqamq0cOFCff/99yoqKgr2OgAAAGMSFrd6XC6XXnjhBblcrmCvgnGCawLDcT3g17gmMJw/roewuIPV2dmpxMREdXR03NFnuBC+uCYwHNcDfo1rAsP543oIiztYAAAA4wkBCwAAwMfGRZdBT0+PJOnkyZN39d93dXVJko4fP66EhASf7YXQxTWB4bge8GtcExjOV9fDrFmzFB9/vQB9XHwG66OPPhr6eYUAAAChaPhnycdFwLp8+bK++OIL5eTkKC4uLtjrAAAA3LFxdwcLAAAgnPAhdwAAAB8jYAEAAPhYSAWsM2fOqLS0VPn5+SouLtaJEydGndu9e7dmzJihvLw8Pf744+rv7w/wpgiE3t5erV27Vvn5+SooKNCKFSt09uzZEXMNDQ2KiopSYWHh0K+6urogbAx/y8nJ0cyZM4e+z3v37h11jjMiMly5cuWmP/f5+fmyWq1qbW29aY4zInw99dRTysnJkcVi0fHjx4ce9zZPSGM4L0wIuf/++01ZWZkxxphPPvnELFq0aMTMTz/9ZBwOh2lubjZut9s8+OCD5o033gjwpgiEnp4e8/nnnxu3222MMeb11183y5YtGzFXX19vEhMTA7wdgiE7O9scO3bstjOcEZFrx44dZtWqVSMe54wIX4cOHTJOp3PE2eBNnjBmbOdFyNzBamlpUXV19VCdw7p16+R0Okfcsfj000+1evVqTZkyRRaLRVu3blV5eXkwVoafxcbGauXKlbJYLJKkkpISNTQ0BHcpjHucEZFr9+7d2rx5c7DXQADdd999yszMvOkxb/OENLbzImQCltPplMPhkNV6vRvVYrEoKytLjY2NN801NjYqOzt76OucnJwRMwhPr732mtasWTPqc93d3SouLlZRUZFefPFFDQ4OBng7BMqGDRs0b948bd68WZcuXRrxPGdEZKqqqlJbW5tWrVo16vOcEZHD2zwhje28CJmABdzOtm3bdPbsWW3fvn3Ecw6HQ01NTTp69Ki+/PJLHT58WDt37gzClvC3b775RrW1taqpqVFaWpoee+yxYK+EcWL37t3asGHD0F+qw3FGwB9CJmBNmzZNzc3NGhgYkCQZY9TY2KisrKyb5rKysnTu3LmhrxsaGkbMILy88sorqqys1IEDB4YK3oaz2WxKT0+XJKWkpGjTpk06fPhwoNdEAPzyZz06OlrPPPPMqN9nzojI09XVpX379mnTpk2jPs8ZEVm8zRPS2M6LkAlY6enpKioq0p49eyRJFRUVyszM1PTp02+aW7dunT777DNdvHhRxhi9/fbbWr9+fTBWRgC8+uqrKi8v18GDB5WUlDTqTEtLy9C/+nC5XKqsrNSCBQsCuSYCoLu7W+3t7UNfl5eXj/p95oyIPHv37lVBQYFmzZo16vOcEZHF2zwhjfG8GMun8wPt1KlTpqSkxMyYMcMsXLjQ1NbWGmOM2bx5s9m/f//Q3DvvvGNyc3NNbm6u2bRpk+nr6wvWyvAjp9NpJJnc3FxTUFBgCgoKzOLFi40xxjz33HPmrbfeMsYYU1FRYebOnWvmz59v5syZY5588knT29sbzNXhB3V1daawsNDMmzfP3HvvvWb16tWmvr7eGMMZEemWLFli3n333Zse44yIDFu2bDEZGRkmKirKpKenm7y8PGPMrfOEMb47L/hROQAAAD4WMm8RAgAAhAoCFgAAgI8RsAAAAHyMgAUAAOBjBCwAAAAfI2ABAAD4GAELAADAxwhYAAAAPkbAAgAA8DECFoCw1NnZqezsbD300EM3Pb5161alpqbqwoULQdoMQCQgYAEIS3a7XWVlZaqsrNSHH34oSTpw4IB27dqlN998U1OnTg3yhgDCGT+LEEBYe/rpp/X+++/r66+/1sqVK7V8+XJ9/PHHwV4LQJgjYAEIaz09PSoqKlJ9fb3S0tL0ww8/KDk5OdhrAQhzvEUIIKzFxcVp7dq1crlcevjhhwlXAAKCO1gAwlptba2Ki4s1d+5cnTx5UjU1NZo9e3aw1wIQ5ghYAMJWX1+fiouLZbfb9dVXX2np0qWSpCNHjshqtQZ5OwDhjLcIAYSt559/XnV1dXrvvfcUExOjDz74QD/++KNeeumlYK8GIMwRsACEpaqqKu3YsUM7d+5UXl6eJGn27Nnavn27Xn75ZVVXVwd5QwDhjLcIAQAAfIw7WAAAAD5GwAIAAPAxAhYAAICPEbAAAAB87P8A5YmooHdfnzwAAAAASUVORK5CYII="
     },
     "execution_count": 358,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u, mesh = solver(L, identity, n_cells)\n",
    "plot(mesh, u, marker=:o, xlabel=:x, lw=5, ylabel=:u, leg=false, \n",
    "     grid=false, title=\"Finite element solution for the elastic bar\", c=[:lightblue])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Mesh convergence - $h$-refinement\n",
    "\n",
    "Use the solver with the loading function $f(x) = \\sin(x)$ to demonstrate that the error in the displacement at the right hand end of the rod scales as $h^{2p}$ when the quadrature is performed to order $p$.\n",
    "\n",
    "Adapt the code given below, replacing the `...` sections.\n",
    "\n",
    "Plot a graph of the displacement error vs. $h$ to illustrate your results for $p=1,2,3$.\n",
    "\n",
    "Now try changing $f(x)$ to something which can be integrated exactly with piecewise linear or quadratic elements. What do you notice?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "n_cells = ... # range of sizes to use for convergence test\n",
    "n_cell_ref = ... # size to use for reference results\n",
    "\n",
    "plot(title=\"Mesh convergence\")\n",
    "for (p, color) in zip([1, 2, 3], \n",
    "                       [:red, :green, :blue])\n",
    "    \n",
    "    # Compute solutions at a range of sizes\n",
    "    solutions = [solver(L, sin, n, quad_degree=p) for n in n_cells]    \n",
    "    \n",
    "    # Compute a reference solution\n",
    "    ref_u, ref_mesh = solver(L, sin, n_cell_ref, quad_degree=p)\n",
    "        \n",
    "    hs = [... for (u, mesh) in solutions ] # find the element size h for each solution\n",
    "    err = [... for (u, mesh) in solutions] # find the error in displacement at end for each\n",
    "\n",
    "    # Plot results\n",
    "    plot!(hs, err, lw=2, legend=true,\n",
    "          marker=:o, ms=5, yscale=:log10, xscale=:log10,\n",
    "          xlabel=\"log h\", ylabel=\"log error\", c=[color], label=\"p=$p\")\n",
    "\n",
    "    # Plot dashed lines to show expected rates\n",
    "    C = err[end] / (hs[end]^(2*p))\n",
    "    plot!(hs, C*hs.^(2.0*p), label=\"h^$(2*p)\", \n",
    "          line=(1, :dash, color), c=[color])\n",
    "end\n",
    "current()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.4",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}