JuliaFEM.jl/notebooks/2015-06-14-data-structures.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data structures\n",
    "\n",
    "Author(s): Jukka Aho (<jukka.aho@kapsi.fi>)\n",
    "\n",
    "####Abstract\n",
    "A proposal how to work with element, nodal etc. data. Main point is that we have one simple principle: *everything is field*. That is, variables are recorded to nodes, element nodes, element gauss points dynamically and there's nothing \"fixed\". This allows us to write very elegant \"microsolvers\" and very easy model structure. Concept will be very easy to unit test. Moreover, when we keep structures as simple as possible, the concept is easy to understand."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we talk about this \"model\". It could be something very easy, like struct of dicts. We don't care *what* we are recording to field or *how many* fields exists for single node/element/model, but we are caring about *where* (in topological sense) we need to save values. For example:\n",
    "\n",
    "### Model fields\n",
    "Global variables like \n",
    "- mesh file name\n",
    "- results file name\n",
    "- time\n",
    "\n",
    "### Nodal fields\n",
    "- primal variable after solution\n",
    "\n",
    "### Element fields\n",
    "- body id\n",
    "- connected nodes\n",
    "- shape function family\n",
    "- dimension\n",
    "- everything else we don't have to interpolate\n",
    "\n",
    "### Element node fields\n",
    "Defined in element local nodes\n",
    "- density\n",
    "- young's modulus\n",
    "- thickness\n",
    "- boundary load\n",
    "- if shape functions are defined, we can interpolate these inside element\n",
    "\n",
    "### Element gauss fields\n",
    "Defined in element gauss points\n",
    "- stress tensor\n",
    "- pretty much everything we need during the solution\n",
    "- these can be extrapolated to element local nodes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, what this means in practice?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "type Model\n",
    "    mo  # For global variables\n",
    "    no  # For nodes\n",
    "    el  # For elements\n",
    "    elno  # For element nodes\n",
    "    elga  # For element gauss points\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "new_model (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function new_model()\n",
    "    \"\"\"Initialize empty model.\n",
    "    \"\"\"\n",
    "    return Model(Dict(), Dict(), Dict(), Dict(), Dict())\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "new_field (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function new_field(field_type, field_name)\n",
    "    d = Dict()\n",
    "    setindex!(field_type, d, field_name)\n",
    "    return d\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And that's it. We can do a lot of things using that simple construction when we are dealing with dynamic languages. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Structure is following:\n",
    "\n",
    "    model.mo[parameter] = value\n",
    "    model.no[field][node_id] = value\n",
    "    model.el[field][element_id] = value\n",
    "    model.elno[field][element_id][local_node_id] = value\n",
    "    model.elga[field][element_id][local_gauss_id] = value\n",
    "\n",
    "`Value` can be scalar, vector, tensor, composite type or whatever is needed. \n",
    "\n",
    "`Field` can also be whatever needed: temperature, displacement, velocity, pressure, density, etc.\n",
    "\n",
    "Examples:\n",
    "\n",
    "    model.mo[\"simulation_time\"] = datetime.datetime.now()\n",
    "    model.no[\"coords\"][node_id] = [1, 2, 3]\n",
    "    model.el[\"color\"][element_id] = \"red\"\n",
    "    model.elno[\"thickness\"][element_id][1] = 5  #  1 here stands for \"first local node of element\"\n",
    "    model.elno[\"thickness\"][element_id][2] = 10\n",
    "    model.elga[\"stress\"][element_id] = [1 2; 3 4]\n",
    "\n",
    "We apply *solvers* to our model. Solver is generic name for a function that uses model data and updates it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1, normal directions in global nodes\n",
    "\n",
    "Next we can start writing our microsolvers. Our goal is to calculate normal direction in nodal points, averaged from element nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: Loading help data...\n"
     ]
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAqMAAAIUCAYAAAAuUCWlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3X10VOW9//1PIEUEhGUxxUpBKUrkyYOiEWiqNZFA58AgKEWEUhLSWgwYIwb4IUIQEYkoIinFVaJFgYA8isYDmATFcCwBouekNgGkVAwKEmixMSBNmPuPHr0JYZQN2bmumf1+rZU/zuRhf2e9z7Zf5mlHBAKBgAAAAAADGpkeAAAAAN7FMgoAAABjWEYBAABgDMsoAAAAjGEZBQAAgDEsowAAADCGZRQAAADGsIwCAADAGJZRAAAAGMMyCgAAAGMuaBn98ssvNX36dPXv31/f//731ahRIy1ZsuScP1taWqr+/fvrsssuU+vWrTVq1ChVVFRc1NAAAAAIDxe0jB45ckQzZ87U7t271aNHD0lSREREnZ8rLy/Xbbfdpr/+9a+aPXu2HnnkEeXm5qpv377617/+dXGTAwAAIORFXsgvXXXVVTp06JB+8IMfaNeuXbrlllvO+XNPPvmkTpw4offff18/+tGPJEkxMTHq27ev/vjHP+rXv/71hU8OAACAkHdBj4w2adJEP/jBDyRJgUAg6M+tWbNGAwYM+GYRlaT4+Hh16tRJr7766oUcGgAAAGHEtTcwHTx4UEeOHNHNN99c53u33HKL3n//fbcODQAAgBDh2jL62WefSZJ++MMf1vneD3/4Qx07dozXjQIAAHjcBb1m9HycOHFCknTJJZfU+V7Tpk2/+Znvfe97tb5XUVGhTZs26ZprrtGll17q1ngAAAC4QCdOnNDf/vY39evXT1dcccVF/S3XltGvF8mvvvqqzvdOnjxZ62fOtGnTJo0cOdKtsQAAAFBPli5dqhEjRlzU33BtGf366fmvn64/02effabWrVvXeVRUkq655hpJ/75znTt3dms8XKB77rlHq1evNj0GzoE2dqOPvWhjL9rYq7S0VCNHjvxmb7sYri2jbdu2VVRUlHbs2FHne0VFRd98PunZvn60tHPnzrrpppvcGg8X6Hvf+x5dLEUbu9HHXrSxF23sVx8vqXT1cqB333233njjDZWXl39zW35+vvbu3auhQ4e6eWi4JDo62vQICII2dqOPvWhjL9p4wwU/MpqVlaV//OMf+vTTTyVJGzZs0IEDByRJDz74oFq2bKkpU6Zo1apVuuOOO5Samqp//vOfevrpp3XDDTcoMTGxfu4BAAAAQtYFL6PPPPOMPv74Y0n/vhTounXrtHbtWkVERGjUqFFq2bKlfvSjH+mdd97Rww8/rMmTJ+uSSy7RgAED9Mwzz5zz9aIAAADwlgteRvfv339eP9elSxdt3LjxQg8DywwYMMD0CAiCNnajj71oYy/aeIOrrxlF+HnjjTdMj4AgaGM3+tiLNvaijTewjMKRjIwM0yMgCNrYjT72oo29aOMNLKNwhI/YsBdt7EYfe9HGXrTxBpZRAAAAGMMyCgAAAGNYRuFIdna26REQBG3sRh970cZetPEGllE4UlxcbHoEBEEbu9HHXrSxF228ISIQCARMD3Gm4uJi9ezZU7t27eKFywAAABaqz32NR0YBAABgDMsoAAAAjGEZBQAAgDEso3DE7/ebHgFB0MZu9LEXbexFG29gGYUj48aNMz0CgqCN3ehjL9rYizbewDIKRxISEkyPgCBoYzf62Is29qKNN7CMAgAAwBiWUQAAABjDMgpH1q9fb3oEBEEbu9HHXrSxF228gWUUjuTk5JgeAUHQxm70sRdt7EUbb+ByoAAAAHCEy4ECAAAgLLCMAgAAwBiWUQAAABjDMgpHEhMTTY+AIGhjN/rYizb2oo03sIzCEa6GYS/a2I0+9qKNvWjjDbybHgAAAI7wbnoAAACEBZZRAAAAGMMyCkcKCwtNj4AgaGM3+tiLNvaijTewjMKRzMxM0yMgCNrYjT72oo29aOMNLKNwZMWKFaZHQBC0sRt97EUbe9HGG1hG4UizZs1Mj4AgaGM3+tiLNvaijTewjAIAAMAYllEAAAAYwzIKR9LT002PgCBoYzf62Is29qKNN7CMwpH27dubHgFB0MZu9LEXbexFG2/gcqAAAABwhMuBAgAAICywjAIAAMAYllE4UlZWZnoEBEEbu9HHXrSxF228gWUUjkycONH0CAiCNnajj71oYy/aeAPLKBzJysoyPQKCoI3d6GMv2tiLNt7AMgpH+JgNe9HGbvSxF23sRRtvYBkFAACAMSyjAAAAMIZlFI7MmTPH9AgIgjZ2o4+9aGMv2ngDyygcqaqqMj0CgqCN3ehjL9rYizbewOVAAQAA4AiXAwUAAEBYYBkFAACAMSyjcKSiosL0CAiCNnajj71oYy/aeAPLKBxJSkoyPQKCoI3d6GMv2tiLNt7AMgpHMjIyTI+AIGhjN/rYizb2oo03sIzCET7hwF60sRt97EUbe9HGG1hGAQAAYAzLKAAAAIxhGYUj2dnZpkdAELSxG33sRRt70cYbWEbhSHFxsekREARt7EYfe9HGXrTxBi4HCgAAAEe4HCgAAADCAssoAAAAjGEZBQAAgDEso3DE7/ebHgFB0MZu9LEXbexFG29gGYUj48aNMz0CgqCN3ehjL9rYizbewDIKRxISEkyPgCBoYzf62Is29qKNN7CMAgAAwBiWUQAAABjj+jK6c+dODRo0SFdddZWaN2+uzp07a+bMmTpx4oTbh4YL1q9fb3oEBEEbu9HHXrSxF228wdVltKSkRLGxsSopKdFDDz2k+fPnq3fv3po+fbqGDx/u5qHhkpycHNMjIAja2I0+9qKNvWjjDZFu/vEVK1bo1KlTys3NVefOnSVJycnJOn36tF5++WUdP35crVq1cnME1LOVK1eaHgFB0MZu9LEXbexFG29w9ZHRSy+9VJL0gx/8oNbtV155pRo3bqwmTZq4eXgAAABYztVlNCkpSW3atNGYMWP0P//zP/rkk0+0cuVKLVq0SA8++OA3yyoAAAC8ydWn6a+66ipt27ZNPp9PN9544ze3T506VY8//ribhwYAAEAIcPWR0cOHD+vnP/+5AoGA/vCHP2jt2rVKSkrSrFmz9Lvf/c7NQ8MliYmJpkdAELSxG33sRRt70cYbXF1GZ86cqYMHD2rLli0aM2aM7rrrLi1evFi/+tWvNGnSJB07dizo7/p8Pvn9/lpfvXv3rvMxD5s3bz7ntWtTUlKUnZ1d67bi4mL5/X5VVFTUun369OmaM2dOrdsOHDggv9+vsrKyWrcvWLBA6enptW6rqqqS3+9XYWFhrdtzcnLOeSINGzYsZO/HmVfDCOX7caZwuR9lZWVhcT/CpcfZ9yMhISEs7ocUHj3OvB9nX+UnVO/H2cLhfiQkJITF/ZBCu0dOTs43u1iHDh3Uo0cPpaWl1fk7FyoiEAgE6u2vnaVHjx667LLL9O6779a6fd26dbr77ruVl5enuLi4Wt8rLi5Wz549tWvXLt10001ujQYAAIALVJ/7mquPjP7rX/9STU3NOW+XpOrqajcPDwAAAMu5uozedNNNKi4u1t69e2vdnpOTo8aNG+uGG25w8/AAAACwnKvLaHp6uho1aqSf/vSneuKJJ7Rw4UL5fD699tprSkpK0pVXXunm4eGCs19/AnvQxm70sRdt7EUbb3B1Gb3hhhv09ttv68Ybb9TTTz+ttLQ07d+/X08++aR+//vfu3louCQzM9P0CAiCNnajj71oYy/aeIOrb2C6ELyByW5VVVVq1qyZ6TFwDrSxG33sRRt70cZeIfMGJoQf/qNgL9rYjT72oo29aOMNLKMAAAAwhmUUAAAAxrCMwpGzr/IAe9DGbvSxF23sRRtvYBmFI+3btzc9AoKgjd3oYy/a2Is23sC76QEAAOAI76YHAABAWGAZBQAAgDEso3CkrKzM9AgIgjZ2o4+9aGMv2ngDyygcmThxoukREARt7EYfe9HGXrTxBpZROJKVlWV6BARBG7vRx160sRdtvIFlFI7wMRv2oo3d6GMv2tiLNt7AMgoAAABjWEYBAABgDMsoHJkzZ47pERAEbexGH3vRxl608QaWUThSVVVlegQEQRu70cdetLEXbbyBy4ECAADAES4HCgAAgLDAMgoAAABjWEbhSEVFhekREARt7EYfe9HGXrTxBpZROJKUlGR6BARBG7vRx160sRdtvIFlFI5kZGSYHgFB0MZu9LEXbexFG29gGYUjfMKBvWhjN/rYizb2oo03sIwCAADAGJZRAAAAGMMyCkeys7NNj4AgaGM3+tiLNvaijTewjMKR4uJi0yMgCNrYjT72oo29aOMNXA4UAAAAjnA5UAAAAIQFllEAAELckSNHlJiYqK5duyo6Olpdu3ZVYmKijhw5Yno04DtFmh4AAABcuM8//1x9+vTRvn37at3+l7/8Re+++67ee+89RUVFGZoO+G48MgpH/H6/6REQBG3sRh97hXqbSZMm1VlEv7Zv3z5NnDixgSeqP6HeBueHZRSOjBs3zvQICII2dqOPvUK9TVFR0UV932ah3gbnh2UUjiQkJJgeAUHQxm70sVeot6murr6o79ss1Nvg/LCMAgAQwiIjv/3tH9/1fcA0llEAAELYlVde+a3fj4mJaaBJgAvDMgpH1q9fb3oEBEEbu9HHXqHcZuPGjXr33XfVvHnzc36/Y8eOyszMbOCp6k8ot8H5YxmFIzk5OaZHQBC0sRt97BWqbQoKCjR48GD169dPu3fv1ujRo9WlSxd16tRJXbp00ejRo0P+Y51CtQ2c4XKgAACEmMLCQvXr10+xsbF67bXX1LRpU9MjwWO4HCgAAB61fft2+Xw+xcTEaN26dSyiCHksowAAhIji4mL169dP3bt31+uvv65mzZqZHgm4aCyjAACEgJKSEiUkJKhTp05688031aJFC9MjAfWCZRSOJCYmmh4BQdDGbvSxVyi0KSsr05133ql27dpp06ZNatWqlemRGkQotMHFYxmFI1wNw160sRt97GV7m48++khxcXGKiorS5s2bdfnll5seqcHY3gb1g3fTAwBgqY8//li33XabmjZtqnfeeec7P+AeaCi8mx4AgDBXXl6uuLg4RUZGqqCggEUUYYsL1gIAYJlDhw4pPj5e1dXV2rp1q9q2bWt6JMA1PDIKRwoLC02PgCBoYzf62Mu2NkeOHFF8fLwqKyuVn5+vq6++2vRIxtjWBu5gGYUjoXyN43BHG7vRx142tTl27JgSEhJUUVGh/Px8XXvttaZHMsqmNnAPT9PDkRUrVpgeAUHQxm70sZctbY4fP67+/fvrk08+0dtvv63rr7/e9EjG2dIG7mIZhSNc7cNetLEbfexlQ5vKykr5fD7t3btX+fn56tatm+mRrGBDG7iPZRQAAIOqqqo0cOBAlZSU6K233uJjDeE5LKMAABhy8uRJDR48WEVFRdq0aZNuvfVW0yMBDY43MMGR9PR00yMgCNrYjT72MtXm1KlTGjp0qLZu3arXX39dsbGxRuawGeeNN/DIKBxp37696REQBG3sRh97mWhTXV2t4cOHa/PmzXrttdcUFxfX4DOEAs4bb+ByoAAANKCamhr98pe/1KpVq7RmzRr5/X7TIwGO1ee+xiOjAAA0kNOnTys5OVkrV67UypUrWUQBsYwCANAgAoGAHnjgAS1ZskSvvPKK7rnnHtMjAVbgDUxwpKyszPQICII2dqOPvRqiTSAQUFpaml544QUtXrxYI0aMcP2Y4YDzxhtYRuHIxIkTTY+AIGhjN/rYy+02gUBAkydP1vz587Vw4UIlJSW5erxwwnnjDSyjcCQrK8v0CAiCNnajj73cbjNjxgxlZmbq2Wef1dixY109VrjhvPEGllE4wsds2Is2dqOPvdxsM3v2bM2YMUOzZ89WWlqaa8cJV5w33sAyCgCAC+bNm6cpU6Zo+vTpmjx5sulxAGuxjAIAUM8WLlyohx9+WJMmTdL06dNNjwNYjWUUjsyZM8f0CAiCNnajj73qu82LL76olJQUpaamavbs2YqIiKjXv+8lnDfewDIKR6qqqkyPgCBoYzf62Ks+2yxbtkzJycn67W9/q3nz5rGIXiTOG29okMuBFhcXKyMjQ9u2bdPJkyf14x//WL/5zW80fvz4c/4slwMFAISaVatW6d5779WoUaOUnZ2tRo14vAfhK6QuB7p582YNHDhQPXv21LRp09SiRQt99NFHOnjwoNuHBgCgQWzYsEH33Xefhg0bpsWLF7OIAg64uox+8cUXGjVqlAYOHKjVq1e7eSgAAIzYuHGjhg4dqkGDBunll19W48aNTY8EhBRX/+m2fPlyff7555o1a5Yk6csvv9Tp06fdPCRcVlFRYXoEBEEbu9HHXhfTpqCgQIMHD1ZCQoKWL1+uyEjXn3D0FM4bb3B1Gc3Ly1PLli31ySefKDo6WpdddplatWqlBx54QF999ZWbh4ZLuIydvWhjN/rY60LbFBYWauDAgbrtttu0atUqNWnSpJ4nA+eNN7i6jO7du1fV1dW666679POf/1xr165VUlKSFi1apMTERDcPDZdkZGSYHgFB0MZu9LHXhbTZvn27fD6fYmJitG7dOjVt2rT+BwPnjUe4+nxCZWWlqqqqNHbsWD333HOSpLvuukunTp3SCy+8oMcff1zXXnutmyOgnvEJB/aijd3oYy+nbYqLi9WvXz91795dr7/+upo1a+bSZOC88QZXHxm99NJLJUnDhw+vdfvX//ef/vSnoL/r8/nk9/trffXu3Vvr16+v9XObN2+W3++v8/spKSnKzs6udVtxcbH8fn+d16BMnz69zgfrHjhwQH6/X2VlZbVuX7BggdLT02vdVlVVJb/fr8LCwlq35+TknPMR4GHDhnE/uB/cD+4H9yME78fWrVuVkJCgTp066c0339TTTz8dkvcjXHpwPxrmfuTk5Hyzi3Xo0EE9evRQWlpanb9zoVz9nNGEhATl5eVp9+7duu666765vaysTF26dNH8+fPrfNYonzMKALBRWVmZbr/9dl111VUqKCjQ5ZdfbnokwJj63NdcfWT05ptvliSVl5fXuv3TTz+VJEVFRbl5eLjg7H/VwR60sRt97HU+bT766CPFxcUpKipKmzdvZhFtIJw33uDqMvqLX/xCUt3/Z1q8eLG+973v6Wc/+5mbh4cLiouLTY+AIGhjN/rY67vafPzxx4qPj9dll12mvLw8HkhpQJw33uD65UCTk5P14osv6he/+IVuu+02vf3221q9erWmTJmiJ554os7P8zQ9AMAW5eXluv322yVJW7duVdu2bQ1PBNghpC4HumjRIrVv314vvfSS1q1bp2uuuUbPPfecHnzwQbcPDQDABTt06JDi4+NVXV3NIgq4yPVlNDIyUtOmTdO0adPcPhQAAPXiyJEjio+PV2Vlpd555x1dffXVpkcCwhbXLQMA4AzHjh1TQkKCKioq9M477/B52IDLXH0DE8LPuT47DXagjd3oY68z2xw/flz9+/fXJ598ovz8fF1//fUGJwPnjTfwyCgcGTdunOkREARt7EYfe33dprKyUj6fT3v37lV+fr66detmeDJw3ngDyygcSUhIMD0CgqCN3ehjr4SEBFVVVWngwIEqKSnRW2+9xae5WILzxhtYRgEAnnby5EkNHjxYRUVF2rRpk2699VbTIwGewjIKAPCsU6dOaejQodq6datyc3MVGxtreiTAc3gDExxZv3696REQBG3sRh/7VFdXa/jw4dq4caPWrVunuLg40yPhLJw33sAyCkdycnJMj4AgaGM3+tilpqZGo0aN0oYNG3Trrbeqf//+pkfCOXDeeAPLKBxZuXKl6REQBG3sRh97nD59WsnJyVq5cqVycnJUWFhoeiQEwXnjDbxmFADgGYFAQA888ICWLFmiV155Rffcc4/pkQDPYxkFAHhCIBBQWlqaXnjhBWVnZ2vEiBGmRwIgnqYHAHhAIBDQ5MmTNX/+fC1cuFBJSUmmRwLwf1hG4UhiYqLpERAEbexGH7NmzJihzMxMPfvssxo7dmyt79HGXrTxBpZROMLVMOxFG7vRx5zZs2drxowZmj17ttLS0up8nzb2oo03RAQCgYDpIc5UXFysnj17ateuXVyODQBwUebNm6eHH35Y06dPV0ZGhulxgLBRn/saj4wCAMLSwoUL9fDDD2vSpEmaPn266XEABMEyCgAIOy+++KJSUlKUmpqq2bNnKyIiwvRIAIJgGYUjfDi0vWhjN/o0nGXLlik5OVm//e1vNW/evO9cRGljL9p4A8soHMnMzDQ9AoKgjd3o0zBWrVqlUaNG6Ve/+pV+97vfndcjorSxF228gWUUjqxYscL0CAiCNnajj/s2bNig++67T8OGDdPixYvVqNH5/U8cbexFG29gGYUjzZo1Mz0CgqCN3ejjro0bN2ro0KEaNGiQXn75ZTVu3Pi8f5c29qKNN7CMAgBCWkFBgQYPHqyEhAQtX75ckZFc6RoIJSyjAICQVVhYqIEDB+q2227TqlWr1KRJE9MjAXCIZRSOpKenmx4BQdDGbvSpf9u3b5fP51NMTIzWrVunpk2bXtDfoY29aOMNLKNwpH379qZHQBC0sRt96ldxcbH69eun7t276/XXX7+o1xbSxl608QYuBwoACCklJSW644479OMf/1hvvfWWWrVqZXokwHO4HCgAwJPKysp05513ql27dtq0aROLKBAGWEYBACHho48+UlxcnKKiorR582ZdfvnlpkcCUA9YRuFIWVmZ6REQBG3sRp+L8/HHHys+Pl6XXXaZ8vLyFBUVVW9/mzb2oo03sIzCkYkTJ5oeAUHQxm70uXDl5eWKi4tTZGSkCgoKdOWVV9br36eNvWjjDXwyMBzJysoyPQKCoI3d6HNhDh06pPj4eFVXV2vr1q1q27ZtvR+DNvaijTewjMIRPmbDXrSxG32cO3LkiOLj41VZWal33nlHV199tSvHoY29aOMNLKMAAOscO3ZMCQkJqqio0DvvvKNrr73W9EgAXMIyCgCwyvHjx9W/f3998sknevvtt3X99debHgmAi3gDExyZM2eO6REQBG3sRp/zU1lZKZ/Pp71792rz5s3q1q2b68ekjb1o4w08MgpHqqqqTI+AIGhjN/p8t6qqKg0cOFAlJSV66623GuwqfLSxF228gcuBAgCMO3nypAYNGqTCwkJt2rRJsbGxpkcC8C3qc1/jkVEAgFGnTp3S0KFDtXXrVuXm5rKIAh7DMgoAMKa6ulrDhw/X5s2b9dprrykuLs70SAAaGG9ggiMVFRWmR0AQtLEbfeqqqanRqFGjtGHDBq1atUr9+/c3Mgdt7EUbb2AZhSNJSUmmR0AQtLEbfWo7ffq0kpOTtXLlSuXk5Mjv9xubhTb2oo038DQ9HMnIyDA9AoKgjd3o8/8LBAJ64IEHtGTJEr3yyiu65557jM5DG3vRxhtYRuEIn3BgL9rYjT7/FggElJaWphdeeEHZ2dkaMWKE6ZFoYzHaeANP0wMAGkQgENDkyZM1f/58LVy4kKdgAUhiGQUANJAZM2YoMzNTzz77rMaOHWt6HACWYBmFI9nZ2aZHQBC0sZvX+8yePVszZszQ7NmzlZaWZnqcWrzexma08QaWUThSXFxsegQEQRu7ebnPvHnzNGXKFE2fPl2TJ082PU4dXm5jO9p4A5cDBQC4ZuHChUpJSdHkyZP15JNPKiIiwvRIAOpBfe5rPDIKAHDFiy++qJSUFD300EMsogCCYhkFANS7ZcuWKTk5WWPHjtWzzz7LIgogKJZRAEC9WrVqlUaNGqXRo0crKyuLRRTAt2IZhSMmL9mHb0cbu3mlz4YNG3Tffffp3nvv1R/+8Ac1amT//8x4pU0ooo032P9fCVhl3LhxpkdAELSxmxf6bNy4UUOHDtWgQYO0ZMkSNW7c2PRI58ULbUIVbbyBZRSOJCQkmB4BQdDGbuHep6CgQIMHD1a/fv20fPlyRUaGztWmw71NKKONN7CMAgAuSmFhoQYOHKjbb79dq1atUpMmTUyPBCCEsIwCAC7Y9u3b5fP5dOutt2rt2rW65JJLTI8EIMSwjMKR9evXmx4BQdDGbuHYp7i4WP369dMNN9ygDRs2qFmzZqZHuiDh2CZc0MYbWEbhSE5OjukREARt7BZufUpKSpSQkKDo6Gi9+eabatGihemRLli4tQkntPEGLgcKAHCkrKxMt99+u6666ioVFBTo8ssvNz0SgAbG5UABAEZ89NFHiouLU1RUlN566y0WUQAXjWUUAHBePv74Y8XHx6tly5bKz8/XFVdcYXokAGGAZRQA8J3Ky8sVFxenyMhI5efnq02bNqZHAhAmWEbhSGJioukREARt7BbKfQ4dOqT4+HhVV1eroKBAbdu2NT1SvQrlNuGONt4QOpfIgBW4Goa9aGO3UO1z5MgRxcfHq7KyUlu3btXVV19teqR6F6ptvIA23sC76QEA53Ts2DHFx8frs88+0zvvvKPo6GjTIwGwRH3uazwyCgCo4/jx4+rfv78++eQTvf322yyiAFzToK8ZnTVrlho1aqTu3bs35GEBAA5UVlbK5/Np7969euutt9StWzfTIwEIYw22jJaXl+vJJ59U8+bNFRER0VCHRT0rLCw0PQKCoI3dQqVPVVWVBg4cqJKSEm3atEk33nij6ZFcFyptvIg23tBgy+gjjzyiPn366Oabb5ZlL1OFA5mZmaZHQBC0sVso9Dl58qQGDx6sHTt26L/+678UExNjeqQGEQptvIo23tAgy+jWrVu1Zs0aPffccwoEAjwyGsJWrFhhegQEQRu72d7n1KlTGjp0qLZu3arXX39dP/nJT0yP1GBsb+NltPEG19/AVFNTo/Hjx+vXv/61unbt6vbh4LJmzZqZHgFB0MZuNveprq7W8OHDtXnzZm3YsEF33HGH6ZEalM1tvI423uD6Mrpo0SIdOHBABQUFbh8KAOBQTU2NRo0apQ0bNmjt2rXq16+f6ZEAeIyry+jRo0c1bdo0TZs2Ta1bt3bzUAAAh06fPq3k5GS9+uqrWrkzT5XNAAAgAElEQVRypQYOHGh6JAAe5OprRqdOnaorrrhC48ePd/y7Pp9Pfr+/1lfv3r21fv36Wj+3efNm+f3+Or+fkpKi7OzsWrcVFxfL7/eroqKi1u3Tp0/XnDlzat124MAB+f1+lZWV1bp9wYIFSk9Pr3VbVVWV/H5/nXf95eTknPNSZsOGDQvZ+3Hm74Ty/ThTuNyP66+/PizuR7j0OPt+pKenW3U/3n33XT3wwANasmSJXn75ZZ06dcpTPc68H2f/7VC9H2cLh/uRnp4eFvdDCu0eOTk53+xiHTp0UI8ePZSWllbn71ywgEv27NkTaNy4cWDBggWB/fv3f/N16623BqKjowN/+9vfAseOHavze7t27QpICuzatcut0XARnn/+edMjIAja2M2mPqdPnw6kpqYGJAVefPFF0+MYZ1Mb1EYbe9Xnvuba5UDffvttxcXFfevPPPTQQ3r22Wdr3cblQAHAPYFAQJMnT1ZmZqZ+//vf67e//a3pkQCEoJC4HGj37t21bt26Wh/jFAgENHXqVFVWVmr+/Pnq2LGjW4cHAJzDjBkzlJmZqXnz5rGIArCCa8to69atNWjQoDq3z5s3T5LO+XoKAIB7Zs+erRkzZuipp57SQw89ZHocAJDUwNeml6SIiAg+9D6Enf3CadiDNnYz3WfevHmaMmWKMjIyNGnSJKOz2MZ0GwRHG29o8GV0y5Yt+t///d+GPizqycSJE02PgCBoYzeTfRYuXKiHH35YkydP1rRp04zNYSvOHXvRxhsafBlFaMvKyjI9AoKgjd1M9XnxxReVkpKihx56SE8++STPTJ0D5469aOMNLKNwpH379qZHQBC0sZuJPsuWLVNycrLGjh2rZ599lkU0CM4de9HGG1hGASAMrVq1SqNGjdLo0aOVlZXFIgrAWiyjABBmNmzYoPvuu0/33nuv/vCHP6hRI/5TD8Be/BcKjpx9eTLYgzZ2a6g+Gzdu1NChQzVo0CAtWbJEjRs3bpDjhjLOHXvRxhtYRuFIVVWV6REQBG3s1hB9CgoKNHjwYPXr10/Lly9XZKRrHyUdVjh37EUbb3DtcqAXisuBAoBzhYWF6tevn37605/qtdde0yWXXGJ6JABhrD73NR4ZBYAQt337dvl8Pt16661au3YtiyiAkMIyCgAhrLi4WP369dMNN9ygDRs2qFmzZqZHAgBHWEbhSEVFhekREARt7OZGn5KSEiUkJCg6OlpvvvmmWrRoUe/H8ALOHXvRxhtYRuFIUlKS6REQBG3sVt99ysrKdOedd6pdu3bauHGjWrZsWa9/30s4d+xFG29gGYUjGRkZpkdAELSxW332+eijjxQXF6eoqCi99dZbuvzyy+vtb3sR5469aOMNLKNwhE84sBdt7FZffT7++GPFx8erZcuWys/P1xVXXFEvf9fLOHfsRRtvYBkFgBBRXl6uuLg4RUZGKj8/X23atDE9EgBcND4RGQBCwKFDhxQfH6/q6mpt3bpVbdu2NT0SANQLHhmFI9nZ2aZHQBC0sdvF9Dly5Iji4+NVWVmpgoICXX311fU4GTh37EUbb2AZhSPFxcWmR0AQtLHbhfY5duyYEhISdPToURUUFKhjx471PBk4d+xFG2/gcqAAYKnjx4+rb9+++utf/6q3335b3bp1Mz0SAEiq332N14wCgIUqKyvl8/m0d+9eFRQUsIgCCFssowBgmaqqKg0cOFAlJSXKy8vTjTfeaHokAHANyygAWOTkyZMaPHiwduzYoU2bNikmJsb0SADgKt7ABEf8fr/pERAEbex2Pn1OnTqloUOHauvWrXr99df1k5/8pAEmA+eOvWjjDTwyCkfGjRtnegQEQRu7fVef6upqDR8+XJs3b9aGDRt0xx13NNBk4NyxF228gWUUjiQkJJgeAUHQxm7f1qempkajRo3Shg0btHbtWvXr168BJwPnjr1o4w0sowBg0OnTp5WcnKxXX31VK1eu1MCBA02PBAANimUUAAwJBAJ64IEHtGTJEi1dulR333236ZEAoMHxBiY4sn79etMjIAja2O3sPoFAQGlpaXrhhReUnZ2t++67z9Bk4NyxF228gWUUjuTk5JgeAUHQxm5n9gkEApo8ebLmz5+v3//+90pMTDQ4GTh37EUbb+ByoADQwDIyMjRjxgzNmzdPDz30kOlxAMAxLgcKACHgyJEjmjhxooqKilRdXa3IyEg1bdpUxcXFeuqpp1hEAUAsowDgis8//1x9+vTRvn376nzv8ssvV1JSkoGpAMA+vGYUAFwwadKkcy6ikvT3v/9dEydObOCJAMBOLKNwhDda2Is2dikqKrqo76PhcO7YizbewDIKR7gahr1oY5fq6uqL+j4aDueOvWjjDSyjcGT48OGmR0AQtLFHZWWljh079q0/ExnJS/ZtwbljL9p4A8soANSTQCCgZcuW6frrr//OZTQmJqaBpgIAu7GMAkA92Llzp2JjYzVy5Ej16tVL27dvV8eOHc/5sx07dlRmZmYDTwgAdmIZhSOFhYWmR0AQtDHj8OHDGjNmjGJiYvTFF18oPz9fq1ev1s0336z33ntPo0ePVpcuXdSuXTt16dJFo0eP1nvvvaeoqCjTo+P/cO7YizbewDIKR3g0x160aVinTp3S3Llzdd1112ndunVasGCB3n//fcXFxX3zM1FRUXrppZf04YcfqkePHvrwww/10ksvsYhahnPHXrTxBl5BD0dWrFhhegQEQZuGk5ubq7S0NO3bt09jx47VjBkz1Lp162/9HfrYizb2oo038MgoHGnWrJnpERAEbdy3e/du+Xw+DRgwQO3atdMHH3ygrKys71xEJfrYjDb2oo03sIwCwHc4fvy4JkyYoG7duqm0tFRr1qxRXl6eunfvbno0AAh5PE0PAEHU1NToj3/8o6ZMmaLKykplZGRowoQJatq0qenRACBs8MgoHElPTzc9AoKgTf3atm2bYmJilJycrL59+2rPnj169NFHL3gRpY+9aGMv2ngDyygcad++vekREARt6kd5eblGjBih2NhYRUREaNu2bVq6dKnatm17UX+XPvaijb1o4w0RgUAgYHqIMxUXF6tnz57atWuXbrrpJtPjAPCIEydO6JlnntHs2bPVokULzZ49W6NHj1ajRvybHQDOVp/7Gq8ZBeBpgUBAa9eu1SOPPKLy8nKlpqbqscceU6tWrUyPBgCewDIKwLNKSkqUmpqqLVu2yOfzaePGjYqOjjY9FgB4Cs8/wZGysjLTIyAI2py/o0ePKiUlRT169NDBgweVm5ur3NxcVxdR+tiLNvaijTewjMKRiRMnmh4BQdDmu1VXVysrK0vXXXedli5dqszMTJWUlMjn87l+bPrYizb2oo038DQ9HMnKyjI9AoKgzbcrKChQamqqPvzwQyUlJWnWrFlq06ZNgx2fPvaijb1o4w08MgpH+JgNe9Hm3Pbv368hQ4YoPj5eLVu21I4dO7R48eIGXUQl+tiMNvaijTewjAIIS5WVlZo6dao6d+6soqIiLVu2TIWFherZs6fp0QAAZ+BpegBhJRAIaPny5Zo0aZIqKiqUnp6uyZMnq3nz5qZHAwCcA4+MwpE5c+aYHgFB0EbauXOnYmNjNXLkSPXq1UulpaWaOXOmFYsofexFG3vRxhtYRuFIVVWV6REQhJfbHD58WGPGjFFMTIy++OIL5efna/Xq1erQoYPp0b7h5T62o429aOMNXA4UQMg6deqUnn/+eT3++OOKjIzUzJkzdf/99ysyklcgAYCbuBwoAM/Lzc1VWlqa9u3bp7Fjx2rGjBlq3bq16bEAAA7xND2AkLJ79275fD4NGDBA7dq10wcffKCsrCwWUQAIUSyjcKSiosL0CAgi3NscP35cEyZMULdu3VRaWqo1a9YoLy9P3bt3Nz3aeQn3PqGMNvaijTewjMKRpKQk0yMgiHBtU1NTo+zsbHXq1EmLFi1SRkaGSktLNWTIEEVERJge77yFa59wQBt70cYbWEbhSEZGhukREEQ4ttm2bZtiYmKUnJysvn37as+ePXr00UfVtGlT06M5Fo59wgVt7EUbb2AZhSN8woG9wqlNeXm5RowYodjYWEVERGjbtm1aunSp2rZta3q0CxZOfcINbexFG29gGQVgjRMnTuiJJ55QdHS08vLylJ2draKiIvXp08f0aAAAl/DRTgCMCwQCWrt2rR555BGVl5crNTVVjz32mFq1amV6NACAy3hkFI5kZ2ebHgFBhGqbkpISxcfH65577lGXLl305z//WXPnzg27RTRU+3gBbexFG29wdRndsWOHxo0bp65du6pFixa6+uqrNWzYMO3du9fNw8JFxcXFpkdAEKHW5ujRo0pJSVGPHj108OBB5ebmKjc3V9HR0aZHc0Wo9fES2tiLNt7g6uVA77nnHr333nsaOnSobrjhBn322WfKyspSZWWl/vSnP6lr1651fofLgQLhrbq6WosWLdK0adNUU1OjadOmafz48WrSpInp0QAA5ylkLgc6YcIE3XLLLbWuEz1s2DB1795dTz31lF555RU3Dw/AMgUFBUpNTdWHH36opKQkzZo1S23atDE9FgDAIFefpu/du3etRVSSrr32WnXp0kVlZWVuHhqARfbv368hQ4YoPj5eLVu21I4dO7R48WIWUQBAw7+BKRAI6PDhw7riiisa+tAAGlhlZaWmTp2qzp07q6ioSMuWLVNhYaF69uxpejQAgCUafBldtmyZPv30Uw0bNqyhD4164Pf7TY+AIGxqEwgEtGzZMl1//fWaO3eu0tPTtXv3bt13330hdQnP+mRTH9RGG3vRxhsadBktKytTSkqK+vTpo1/96lcNeWjUk3HjxpkeAUHY0mbnzp2KjY3VyJEj1atXL5WWlmrmzJlq3ry56dGMsqUP6qKNvWjjDQ22jB46dEj/+Z//qcsvv1yrV6/+zkdHfD6f/H5/ra/evXtr/fr1tX5u8+bN5/yXU0pKSp3PJysuLpbf71dFRUWt26dPn645c+bUuu3AgQPy+/11Xtu6YMECpaen17qtqqpKfr9fhYWFtW7PyclRYmJindmGDRsWsvcjISEhLO7HmcLlfmRnZxu9H4cPH1bv3r11yy236IsvvlB+fr5Wr16tNm3aeLLH2fcjISEhLO6HFB49zrwfZ/53LZTvx9nC4X4kJCSExf2QQrtHTk7ON7tYhw4d1KNHD6WlpdX5OxfK1Y92+trx48f1s5/9TOXl5Xr33Xd1/fXXB/1ZPtoJCC2nTp3S888/r8cff1yRkZGaOXOm7r///jpvXgQAhI+Q+WgnSTp58qQGDhyojz76SHl5ed+6iAIILbm5uUpLS9O+ffs0duxYzZgxQ61btzY9FgAghLj6NH1NTY2GDRum7du3a9WqVbr11lvdPBwawNkP48MeDdlm9+7d8vl8GjBggNq1a6cPPvhAWVlZLKLfgnPHXrSxF228wdVldMKECXr99dfVv39/VVRUaOnSpbW+EHpycnJMj4AgGqLN8ePHNWHCBHXr1k2lpaVas2aN8vLy1L17d9ePHeo4d+xFG3vRxhtcfc3oHXfcoa1bt+pch4iIiFBNTU2d23nNKGCfmpoa/fGPf9SUKVNUWVmpKVOmaMKECWratKnp0QAABoTMa0a3bNni5p8H0AC2bdumBx98UMXFxRoxYoTmzJmjtm3bmh4LABAmGvxD7wGEhvLyco0YMUKxsbGKiIjQtm3btHTpUhZRAEC9YhkFUMuJEyf0xBNPKDo6Wnl5ecrOzlZRUZH69OljejQAQBhiGYUj5/qwXNjhYtsEAgGtWbNGXbp00YwZMzR27Fjt2bNHSUlJatSI/1RcLM4de9HGXrTxBj6VGo6cfaUS2ONi2pSUlCg1NVVbtmyRz+fTxo0bFR0dXY/TgXPHXrSxF228gYc74Mjw4cNNj4AgLqTN0aNHlZKSoh49eujgwYPKzc1Vbm4ui6gLOHfsRRt70cYbeGQU8KDq6motWrRI06ZNU01NjTIzMzV+/Hg1adLE9GgAAI9hGQU8pqCgQKmpqfrwww+VlJSkWbNmqU2bNqbHAgB4FE/Tw5HCwkLTIyCI72qzf/9+DRkyRPHx8WrZsqV27NihxYsXs4g2EM4de9HGXrTxBpZROJKZmWl6BAQRrM2XX36pqVOnqnPnzioqKtKyZctUWFionj17NvCE3sa5Yy/a2Is23sDT9HBkxYoVpkdAEGe3CQQCWr58uSZNmqSKigqlp6dr8uTJat68uaEJvY1zx160sRdtvIFHRuFIs2bNTI+AIM5ss3PnTsXGxmrkyJHq1auXSktLNXPmTBZRgzh37EUbe9HGG1hGgTBy+PBhjRkzRjExMfriiy+Un5+v1atXq0OHDqZHAwDgnHiaHggDp06d0vPPP6/HH39ckZGRWrBgge6//35FRnKKAwDsxiOjcCQ9Pd30CDhLbm6uunXrpvT0dI0aNUp79+5VSkoKi6hlOHfsRRt70cYbWEbhSPv27U2PgP+ze/du+Xw+DRgwQO3atdPkyZOVlZWl1q1bmx4N58C5Yy/a2Is23sBDJ3Bk/PjxpkfwvOPHj+vxxx/X888/rx/96Edas2aNBg8erIiICNOj4Vtw7tiLNvaijTewjAIh4vTp03rppZc0ZcoUVVZWKiMjQxMmTFDTpk1NjwYAwAXjaXogBGzbtk0xMTFKTk5W3759tWfPHj366KMsogCAkMcyCkfKyspMj+Ap5eXlGjFihGJjYyX9eyldunSp2rZtW+dnaWM3+tiLNvaijTewjMKRiRMnmh7BE06cOKEnnnhC0dHRysvLU3Z2toqKitSnT5+gv0Mbu9HHXrSxF228gdeMwpGsrCzTI4S1QCCgtWvX6pFHHlF5eblSU1P12GOPqVWrVt/5u7SxG33sRRt70cYbWEbhCB+z4Z6SkhKlpqZqy5Yt8vl82rhxo6Kjo8/792ljN/rYizb2oo038DQ9YNjRo0eVkpKiHj166ODBg8rNzVVubq6jRRQAgFDFI6OAIdXV1Vq0aJGmTZummpoaZWZmavz48WrSpInp0QAAaDA8MgpH5syZY3qEsFBQUKAbb7xRDz74oIYMGaI9e/ZowoQJF7WI0sZu9LEXbexFG29gGYUjVVVVpkcIafv379eQIUMUHx+vli1baseOHVq8eLHatGlz0X+bNnajj71oYy/aeENEIBAImB7iTMXFxerZs6d27dqlm266yfQ4QL348ssvNXv2bM2dO1dXXHGFMjMzNXz4cC7hCQAISfW5r/GaUcBFgUBAy5cv16RJk1RRUaH09HRNnjxZzZs3Nz0aAABW4Gl6wCU7d+5UbGysRo4cqV69eqm0tFQzZ85kEQUA4Awso3CkoqLC9AjWO3z4sMaMGaOYmBh98cUXys/P1+rVq9WhQwdXj0sbu9HHXrSxF228gWUUjiQlJZkewVqnTp3S3Llzdd1112ndunVasGCB3n//fcXFxTXI8WljN/rYizb2oo038JpROJKRkWF6BCvl5uYqLS1N+/bt09ixYzVjxgy1bt26QWegjd3oYy/a2Is23sAjo3CETziobffu3fL5fBowYIDatWunDz74QFlZWQ2+iEq0sR197EUbe9HGG1hGgQtw/PhxTZgwQd26dVNpaanWrFmjvLw8de/e3fRoAACEFJ6mBxw4ffq0XnrpJU2ZMkWVlZXKyMjQhAkT1LRpU9OjAQAQknhkFI5kZ2ebHsGYbdu2KSYmRsnJyerbt6/27NmjRx991JpF1MttQgF97EUbe9HGG1hG4UhxcbHpERpceXm5RowYodjYWEn/XkqXLl2qtm3bGp6sNi+2CSX0sRdt7EUbb+ByoEAQJ06c0DPPPKPZs2erRYsWmj17tkaPHq1Gjfg3HADA27gcKOCiQCCgtWvX6pFHHlF5eblSU1P12GOPqVWrVqZHAwAg7LCMAmcoKSlRamqqtmzZIp/Pp40bNyo6Otr0WAAAhC2ebwQkHT16VCkpKerRo4cOHjyo3Nxc5ebmsogCAOAyllE44vf7TY9Qr6qrq5WVlaXrrrtOS5cuVWZmpkpKSuTz+UyP5li4tQk39LEXbexFG2/gaXo4Mm7cONMj1JuCggKlpqbqww8/VFJSkmbNmqU2bdqYHuuChVObcEQfe9HGXrTxBh4ZhSMJCQmmR7ho+/fv15AhQxQfH6+WLVtqx44dWrx4cUgvolJ4tAln9LEXbexFG29gGYVnfPnll5o6dao6d+6soqIiLVu2TIWFherZs6fp0QAA8CyepkfYCwQCWr58uSZNmqSKigqlp6dr8uTJat68uenRAADwPB4ZhSPr1683PYIjO3fuVGxsrEaOHKlevXqptLRUM2fODMtFNNTaeA197EUbe9HGG1hG4UhOTo7pEc7L4cOHNWbMGMXExOiLL75Qfn6+Vq9erQ4dOpgezTWh0sar6GMv2tiLNt7A5UARVk6dOqXnn39ejz/+uCIjIzVz5kzdf//9iozkFSkAANQXLgcKnENubq7S0tK0b98+jR07VjNmzFDr1q1NjwUAAL4FT9Mj5O3evVs+n08DBgxQu3bt9MEHHygrK4tFFACAEMAyipB1/PhxTZgwQd26dVNpaanWrFmjvLw8de/e3fRoAADgPLGMwpHExETTI+j06dPKzs5Wp06dtGjRImVkZKi0tFRDhgxRRESE6fGMsaENgqOPvWhjL9p4A8soHDF9NYxt27YpJiZGycnJ6tu3r/bs2aNHH31UTZs2NTqXDUy3wbejj71oYy/aeAPLKBwZPny4keOWl5drxIgRio2NlfTvpXTp0qVq27atkXlsZKoNzg997EUbe9HGG1hGYbUTJ07oiSeeUHR0tPLy8pSdna2ioiL16dPH9GgAAKAe8NFOsFIgENDatWv1yCOPqLy8XKmpqXrsscfUqlUr06MBAIB6xCOjcKSwsND1Y5SUlCg+Pl733HOPunTpoj//+c+aO3cui+h3aIg2uHD0sRdt7EUbb2AZhSOZmZmu/e2jR48qJSVFPXr00MGDB5Wbm6vc3FxFR0e7dsxw4mYbXDz62Is29qKNN/A0PRxZsWJFvf/N6upqLVq0SNOmTVNNTY0yMzM1fvx4NWnSpN6PFc7caIP6Qx970cZetPEGllE40qxZs3r9ewUFBUpNTdWHH36opKQkzZo1S23atKnXY3hFfbdB/aKPvWhjL9p4A0/Tw4j9+/fr7rvvVnx8vFq2bKkdO3Zo8eLFLKIAAHgMyyga1JdffqmpU6eqc+fO2r59u5YtW6bCwkL17NnT9GgAAMAAllE4kp6efkG/FwgEtGzZMkVHR2vu3LlKT0/X7t27dd9993n6Ep716ULboGHQx160sRdtvIFlFI60b9/e8e/s3LlTsbGxGjlypHr16qXS0lLNnDlTzZs3d2FC77qQNmg49LEXbexFG2+ICAQCAdNDnKm4uFg9e/bUrl27dNNNN5keBxfh8OHDmjJlil566SV17dpV8+fPV1xcnOmxAADARarPfc31R0a/+uorTZo0SVdddZWaNWumXr16KS8vz+3DwqBTp05p7ty5uu6667Ru3TotWLBA77//PosoAACow/VldPTo0Zo3b55++ctf6vnnn1fjxo3l8/m0bds2tw8NA3Jzc9WtWzdNmjRJo0aN0t69e5WSkqLISD5FDAAA1OXqMlpUVKSVK1fqqaee0pw5c5ScnKyCggJdffXVmjhxopuHhkvKysrOefvu3bvl8/k0YMAAtWvXTh988IGysrLUunXrBp7Qu4K1gR3oYy/a2Is23uDqMrp69WpFRkbqN7/5zTe3XXLJJRozZozee+89HTx40M3Do54cOXJEiYmJ6tq1q26++WZ17dpViYmJOnLkiI4fP64JEyaoW7duKi0t1Zo1a5SXl6fu3bubHttz+Aee3ehjL9rYizbe4Opzp++//746deqkFi1a1Lr9lltukSR98MEHatu2rZsj4CJ9/vnn6tOnj/bt2/fNbX/5y1/0l7/8Rbm5uQoEAqqqqlJGRoYmTJigpk2bGpzW27KyskyPgG9BH3vRxl608QZXl9HPPvtMP/zhD+vc/vVtn376qZuHRz2YNGlSrUX0TEeOHNGPf/xj/lFhCT4CxW70sRdt7EUbb3D1afoTJ07okksuqXP714+enThxws3Dox4UFRV96/ebNm3KIgoAAC6Yq8vopZdeqq+++qrO7SdPnvzm+8H4fD75/f5aX71799b69etr/dzmzZvl9/vr/H5KSoqys7Nr3VZcXCy/36+Kiopat0+fPl1z5sypdduBAwfk9/vrvHh6wYIFda4IUVVVJb/fr8LCwlq35+TkKDExsc5sw4YNC5n7UV1dXWemM/3jH/8IifshhUcP7gf3g/vB/eB+cD8a+n7k5OR8s4t16NBBPXr0UFpaWp2/c8ECLrrzzjsDXbp0qXN7Xl5eICIiIvDGG2/U+d6uXbsCkgK7du1yczScpy5dugQkBf06V1+Y8dRTT5keAd+CPvaijb1oY6/63NdcfWT0xhtv1J49e/TPf/6z1u3bt2+XJPXo0cPNw6MexMTEXNT30XCqqqpMj4BvQR970cZetPEGVy8HWlRUpF69eunpp5/WhAkTJP37ikzdunVTVFSU/vu//7vO73A5ULscOXJEvXv3PuebmDp27Kj33ntPUVFRBiYDAACm1Oe+5uq76WNiYjR06FD9v//3//T555+rY8eOWrJkiQ4cOKCXXnrJzUOjnkRFRem9997TxIkTVVRUpOrqakVGRiomJkaZmZksogAA4KK4fo3Gl19+WY899pheeeUV/f3vf9d//Md/6I033lBsbKzbh0Y9iYqK4h8PAADAFa5fm/6SSy5RZmamPv30U504cUJ/+tOf1LdvX7cPC5ec/Y4/2IM2dqOPvWhjL9p4g+vLKMJLUlKS6REQBG3sRh970cZetPEGllE4kpGRYXoEBEEbu9HHXrSxF228gWUUjvAJB+OPxFQAAA+JSURBVPaijd3oYy/a2Is23sAyCgAAAGNYRgEAAGAMyygcOfu6urAHbexGH3vRxl608QaWUThSXFxsegQEQRu70cdetLEXbbzB1cuBXgguBwoAAGC3+tzXeGQUAAAAxrCMAgAAwBiWUQAAABjDMgpH/H6/6REQBG3sRh970cZetPEGllE4Mm7cONMjIAja2I0+9qKNvWjjDSyjcCQhIcH0CAiCNnajj71oYy/aeAPLKAAAAIxhGQUAAIAxLKNwZP369aZHQBC0sRt97EUbe9HGG1hG4UhOTo7pERAEbexGH3vRxl608QYuBwoAAABHuBwoAAAAwgLLKAAAAIxhGQUAAIAxLKNwJDEx0fQICII2dqOPvWhjL9p4A8soHOFqGPaijd3oYy/a2Is23sC76QEAAOAI76YHAABAWGAZBQAAgDEso3CksLDQ9AgIgjZ2o4+9aGMv2ngDyygcyczMND0CgqCN3ehjL9rYizbewDIKR1asWGF6BARBG7vRx160sRdtvIFlFI40a9bM9AgIgjZ2o4+9aGMv2ngDyygAAACMYRkFAACAMSyjcCQ9Pd30CAiCNnajj71oYy/aeAPLKBxp37696REQBG3sRh970cZetPEGLgcKAAAAR7gcKAAAAMICyygAAACMYRmFI2VlZaZHQBC0sRt97EUbe9HGG1hG4cjEiRNNj4AgaGM3+tiLNvaijTewjMKRrKws0yMgCNrYjT72oo29aOMNLKNwhI/ZsBdt7EYfe9HGXrTxBpZRAAAAGMMyCgAAAGNYRuHInDlzTI+AIGhjN/rYizb2oo03sIzCkaqqKtMjIAja2I0+9qKNvWjjDVwOFAAAAI5wOVAAAACEBZZRAAAAGMMyCkcqKipMj4AgaGM3+tiLNvaijTewjMKRpKQk0yMgCNrYjT72oo29aOMNLKNwJCMjw/QICII2dqOPvWhjL9p4A8soHOETDuxFG7vRx160sRdtvIFlFAAAAMawjAIAAMAYllE4kp2dbXoEBEEbu9HHXrSxF228gWUUjhQXF5seAUHQxm70sRdt7EUbb+ByoAAAAHCEy4ECAAAgLLCMAgAAwBiWUQAAABjDMgpH/H6/6REQBG3sRh970cZetPEGllE4Mm7cONMjIAja2I0+9qKNvWjjDSyjcCQhIcH0CAiCNnajj71oYy/aeAPLKAAAAIxhGQUAAIAxri6j+fn5SkpKUqdOndS8eXN17NhRv/71r3Xo0CE3DwsXrV+/3vQICII2dqOPvWhjL9p4g6vL6KRJk7R161bdfffdWrBgge699169+uqruvHGG3X48GE3Dw2XzJkzx/QICII2dqOPvWhjL9p4Q6Sbf/y5555TbGxsrdv69++v22+/XVlZWZo5c6abh4cLoqKiTI+AIGhjN/rYizb2oo03uPrI6NmLqCT99Kc/1fe//32VlZW5eWgAAACEgAZ/A1NlZaX++c9/6oorrmjoQwMAAMAyDb6MPvfcc/rXv/6lYcOGNfShAQAAYJnzfs1oIBDQV199dV4/27Rp03PevnXrVs2YMUPDhg3Tz372s3P+zIkTJyRJpaWl5zsaGlBRUZGKi4tNj4FzoI3d6GMv2tiLNvb6ek/7em+7KIHztGXLlkBERMR5fe3evbvO75eWlga+//3vB2666aZAZWVl0OMsXbo0IIkvvvjiiy+++OKLL8u/li5der6rZFARgUAgoPNw+PBhbdq06Xx+VHfddZdatmz5zf/9ySef6Cc/+YmaNGmibdu2qU2bNkF/t6KiQps2bdI111yjSy+99LyOBwAAgIZz4sQJ/e1vf1O/fv0u+n1A572MXqijR48qNjZW//jHP1RYWKiOHTu6eTgAAACEEFc/Z/TLL7+Uz+fTZ599pi1btrCIAgAAoBZXHxm96667tGHDBiUlJdV5w9Jll12mQYMGuXVoAAAAhABXl9EOHTrowIEDOtchrrnmGv31r39169AAAAAIAa6/ZhQAAAAIpsE/9B4AAAD4mtXLaH5+vpKSktSpUyc1b95c/1979xrS5PvGAfy7lYeartwyc0VqHjqYkqVoR6TC1LIssyF0AM3ouOpFBCGGmVR0wMoXagf1RUalKCERstK0UlPKikxaQeZhWVmWVtqa9/+FbeRvqwb99N7/91wf2Itdz8OeL4xnu57TfXt6eiIxMRGvX7/mHU1Qent7sXfvXigUCowcORIhISFQq9W8YwlebW0ttm/fDl9fXzg4OMDNzQ1KpRIajYZ3NGJGWloaxGIx/Pz8eEchP9y/fx/Lly+HXC6HRCKBn58fTp8+zTuW4NXV1WHFihVQKBSQSCSYOnUqUlNT/53B1YnFPn/+jP379yM8PBwymQxisRh5eXlm13369CnCw8Ph6OgIuVyO9evX4927dxZvy6ov0wcGBqKzsxOxsbHw9vbGixcvkJGRgZEjR6K+vv6345WSf09cXBwKCwuxe/dueHt7IycnB7W1tSgrK8PcuXN5xxOs1atXo6qqCrGxsfD394dWq0VGRga6u7tRXV0NX19f3hHJDy0tLZg8eTLEYjE8PDzw6NEj3pEEr7S0FFFRUZg1axaUSiUcHBzw/PlzMMZw+PBh3vEE6/HjxwgKCoJCocDmzZshk8lw9+5d5ObmYvny5SguLuYdUTBevnyJSZMmwc3NDR4eHigvL0dubi7Wr18/YL2WlhYEBATAyckJKpUKXV1dOHbsGCZOnIh79+7Bxsbmzxv762HzB1FlZaVJraKigolEIpaUlMQhkfDU1NQwkUjEjh8/bqz19PQwLy8vNmfOHI7JyN27d5lOpxtQ02g0zN7enq1du5ZTKmKOUqlkixcvZqGhoWz69Om84wjex48fmYuLC4uJieEdhfzDvn37mEgkYg0NDQPqGzZsYCKRiHV2dnJKJjy9vb2svb2dMcZYXV0dE4lELC8vz2S9LVu2MIlEwpqbm401tVrNRCIRy87OtmhbVn2Zft68eSa1+fPnQyaTobGxkUMi4SkoKMDw4cOxadMmY83Ozg4JCQmoqqpCa2srx3TCNnv2bAwfPnCoYC8vL0ybNo32DytSUVGBwsJCpKengzEGkUjEO5Lg5efn482bN0hLSwPQfzmyr6+PcyoCwDjz4tixYwfUx40bh2HDhsHW1pZHLEGytbU1fg/sNxfRCwsLsWzZMkyYMMFYW7RoEXx8fHD58mWLtmXVzag53d3d6Orq+uupp4hlHjx4AB8fHzg4OAyoBwUFAQDq6+t5xCK/wBhDe3s77R9WQq/XY8eOHUhMTKTbJqyIWq2GVCpFc3MzJk+eDEdHR4waNQpbt25Fb28v73iCFh8fDxcXFyQkJODhw4dobm7GpUuXkJmZCZVKRdOEW5nW1la8ffsWgYGBJsuCgoLw4MEDiz5nUGdgGgzp6enQ6XRQKpW8owiCVquFq6urSd1Qa2trG+pI5DcuXLiAtrY2HDx4kHcUAiAzMxOvXr3CzZs3eUchP9FoNPj+/Tuio6OxceNGHDlyBGVlZTh9+jQ6OzuRn5/PO6JgKRQK3LlzB5GRkQgICDDWk5KScODAAY7JiDlarRYAftknvH//Hjqd7o/3jQ5ZM8oYs/iI097e3my9oqICKSkpUCqVJjM6kcHx9etX2NnZmdQN3xE93Wg9GhsbsW3bNsyZMwcbNmzgHUfwOjo6kJycjOTkZMjlct5xyE+6u7vx5csXbNmyBenp6QD6Zwz89u0bsrKycODAAXh5eXFOKUzt7e2IiIgAAJw5cwZyuRwlJSVIS0uDi4sLtm3bxjkh+ZmhB/hTn2A1zeitW7ewcOFCi9ZtbGyEj4+PSW3lypXw9/fH2bNnByMiMWPEiBFmDyJ6enqMywl/r1+/xtKlS+Hk5ISCggK6L9EKJCUlYcyYMdixYwfvKOQfDL9bcXFxA+pxcXHIyspCdXU1NaOcpKamorW1Fc+ePYNCoQDQf6DQ19eHvXv3Ii4uDjKZjHNKYmDYl/62TxiyZnTq1KnIzc21aN1x48YNeN/c3IywsDA4OTnh2rVrkEgkg5CQmOPq6mr2Urzh1Lzhx4Lw8/HjR0RERODTp0+orKw02X/I0NNoNDhz5gzS09PR0tJirPf09ODbt29oamqCVCqFk5MTx5TCpVAo0NDQYDI8oOFhjQ8fPvCIRQDcvn0bAQEBJv8tUVFRyM3NRX19vcUntsjgM1yeN/QEP9NqtZDL5RYN7TRkzaiLi4vJ2FSW6OjoQFhYGHQ6HcrKymhs0SEWEBCA8vJydHV1wdHR0VivqakBAMyYMYNXNIL+5iYqKgrPnz+HWq3GlClTeEci6L+pv6+vDyqVCiqVymS5h4cHdu3ahRMnTnBIRwIDA6FWq9HS0gJvb29j3XDg7ezszCua4Ol0Ouj1erN1APj+/ftQRyK/MX78eDg7O6O2ttZk2b179yzuEaz6afrPnz8jMjISWq0W165dg6enJ+9IgrN69Wro9XpkZ2cba729vcjJyUFISAjGjx/PMZ2w6fV6KJVK1NTU4MqVKwgODuYdifzg5+eHoqIiFBcXG19FRUXw9fWFm5sbiouLkZCQwDumYK1ZswYAcO7cuQH1s2fPwsbGhp5J4GjmzJm4f/++yUxyFy9exLBhw+Dv788pGfmVmJgYlJSUDLgKdOPGDWg0GsTGxlr0GVY9A1N0dDSuXr2K+Ph4kx8HR0dHrFixgk8wgVEqlSgqKsLu3bvh6emJvLw81NXV4caNG2bHgiVDY9euXTh16hSioqLM7vBr167lkIr8TmhoKDo6OvD48WPeUQRv48aNOH/+PNasWYMFCxagvLwcBQUF2LdvH41GwdGjR48QEhICqVSK7du3QyaToaSkBNevX0diYiKysrJ4RxSUjIwMdHZ2oq2tDZmZmVi1apXxbKdKpYJUKjXOwDR69Gjs3LkTXV1dOHr0KCZOnIja2tr//xmY3N3dmVgsZiKRyOTl4eHBO55g9PT0sD179jBXV1dmb2/PgoODWWlpKe9YghcaGvrL/UMsFvOOR8wIDQ1lfn5+vGMQxphOp2MpKSnM3d2d2draMh8fH3by5EnesQjrn/kvPDycSaVSZmtry6ZMmcIOHTrE9Ho972iC4+7uPuB/xfCfIxaLWVNTk3G9J0+esCVLljCJRMJkMhlbt24de/PmjcXbseozo4QQQggh5L/Nqu8ZJYQQQggh/23UjBJCCCGEEG6oGSWEEEIIIdxQM0oIIYQQQrihZpQQQgghhHBDzSghhBBCCOGGmlFCCCGEEMINNaOEEEIIIYQbakYJIYQQQgg31IwSQgghhBBuqBklhBBCCCHcUDNKCCGEEEK4+R/1Av42w/Zt6gAAAABJRU5ErkJggg=="
      ],
      "text/plain": [
       "Figure(PyObject <matplotlib.figure.Figure object at 0x10dd04390>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "(-2,10)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using PyPlot\n",
    "plot([0, 4, 7], [0, 3, 7], \"-ko\")\n",
    "grid()\n",
    "xlim(-2, 10)\n",
    "ylim(-2, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First obvious thing is to generate or read geometry to model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model(Dict{Any,Any}(),{\"connects_to\"=>{2=>[1,2],3=>[2],1=>[1]},\"coords\"=>{2=>[4.0,3.0],3=>[7.0,7.0],1=>[0.0,0.0]}},{\"connectivity\"=>{2=>[2,3],1=>[1,2]},\"dimension\"=>{2=>1,1=>1}},Dict{Any,Any}(),Dict{Any,Any}())\n"
     ]
    }
   ],
   "source": [
    "function dummymeshreader(model)\n",
    "    \"\"\"Read mesh to model. Rewrite older mesh. This microsolver is for demonstration purposes.\"\"\"\n",
    "\n",
    "    # Create new field for node coordinates. Coordinates is nodal field.\n",
    "    coords = new_field(model.no, \"coords\")\n",
    "    reverse_connectivity = new_field(model.no, \"connects_to\")\n",
    "\n",
    "    # Create new field for connectivity and dimension. This is element property.\n",
    "    connectivity = new_field(model.el, \"connectivity\")\n",
    "    dimension = new_field(model.el, \"dimension\")\n",
    "\n",
    "    # Do actual reading, call another microsolver which reads geometry e.g. ABAQUS\n",
    "\n",
    "    # Here we just create three nodes and two linear segments to demonstrate the concept\n",
    "    coords[1] = [0.0, 0.0]\n",
    "    coords[2] = [4.0, 3.0]\n",
    "    coords[3] = [7.0, 7.0]\n",
    "    connectivity[1] = [1, 2]\n",
    "    connectivity[2] = [2, 3]\n",
    "    dimension[1] = 1\n",
    "    dimension[2] = 1\n",
    "    reverse_connectivity[1] = [1]\n",
    "    reverse_connectivity[2] = [1, 2]\n",
    "    reverse_connectivity[3] = [2]\n",
    "    return model\n",
    "end\n",
    "\n",
    "function test_dummymeshreader()\n",
    "    m = new_model() # Create new model\n",
    "    dummymeshreader(m) # Apply microsolver \"dummymeshreader\" to model m\n",
    "    println(m)\n",
    "    @assert m.no[\"coords\"][1] == [0.0, 0.0]\n",
    "    @assert m.el[\"connectivity\"][1] == [1, 2]\n",
    "end\n",
    "\n",
    "test_dummymeshreader()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our next microsolver assigns some function space to elements. We use linear Lagrange polynomials. We test routine with some simple interpolation. We want to find midpoint of first element [0.0, 0.0] -> [4.0, 3.0] but because we can interpolate only from element nodes, we need some utility functions to move field variables from nodes to element nodes. So we need to write another microsolver which makes such transformation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{2=>[4.0,3.0],1=>[0.0,0.0]}\n"
     ]
    }
   ],
   "source": [
    "function no2elno(model, field)\n",
    "    \"\"\"Move field from nodes to element nodes.\"\"\"\n",
    "    f = new_field(model.elno, field)\n",
    "    for (elid, node_ids) in model.el[\"connectivity\"]  # loop through element connectivity list\n",
    "        f[elid] = Dict()\n",
    "        for (i, c) in enumerate(node_ids)\n",
    "            f[elid][i] = model.no[field][c]\n",
    "        end\n",
    "    end\n",
    "end\n",
    "\n",
    "function test_no2elno()\n",
    "    m = new_model()\n",
    "    dummymeshreader(m)\n",
    "    no2elno(m, \"coords\")\n",
    "    @assert \"coords\" in keys(m.elno)  # new element node field\n",
    "    el = m.elno[\"coords\"][1]  # first element ..\n",
    "    println(el)\n",
    "    @assert el[1] == [0.0, 0.0]\n",
    "    @assert el[2] == [4.0, 3.0]\n",
    "end\n",
    "\n",
    "test_no2elno()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.0 0.0]\n",
      "[0.0 1.0]\n",
      "{2=>[4.0,3.0],1=>[0.0,0.0]}\n",
      "[2.0,1.5]\n"
     ]
    }
   ],
   "source": [
    "function get_shape_functions(family, dim, order)\n",
    "    if family == \"Lagrange\"\n",
    "        if dim == 1\n",
    "            if order == 1\n",
    "                return xi -> [0.5*(1-xi) 0.5*(1+xi)]\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    throw(\"Unknown shape function requested.\")\n",
    "end\n",
    "\n",
    "function interpolate(shape_functions, elnofield, xi)\n",
    "    result = zeros(size(elnofield[1]))\n",
    "    sf = shape_functions(xi)\n",
    "    for i=1:length(sf)\n",
    "        result += sf[i]*elnofield[i]\n",
    "    end\n",
    "    return result\n",
    "end\n",
    "\n",
    "function test_shape_functions()\n",
    "    sf = get_shape_functions(\"Lagrange\", 1, 1)\n",
    "    println(sf(-1))\n",
    "    println(sf(1))\n",
    "    @assert sf(-1) == [1 0]\n",
    "    @assert sf(1) == [0 1]\n",
    "end\n",
    "\n",
    "function test_define_fuction_space()\n",
    "    m = new_model()\n",
    "    dummymeshreader(m) # Read mesh\n",
    "    no2elno(m, \"coords\") # Move nodal field (coordinates) from nodes to element nodes\n",
    "    sf = get_shape_functions(\"Lagrange\", 1, 1)\n",
    "    midpoint = interpolate(sf, m.elno[\"coords\"][1], 0.0)\n",
    "    println(m.elno[\"coords\"][1])\n",
    "    println(midpoint)\n",
    "    @assert midpoint == [2.0, 1.5]\n",
    "end\n",
    "\n",
    "test_shape_functions()\n",
    "test_define_fuction_space()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our next microsolver introduces new field, namely normal direction in *element nodes*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.6,-0.8]\n",
      "{2=>[0.6,-0.8],1=>[0.6,-0.8]}\n"
     ]
    }
   ],
   "source": [
    "function get_shape_function_derivatives(family, dim, order)\n",
    "    if family == \"Lagrange\"\n",
    "        if dim == 1\n",
    "            if order == 1\n",
    "                return xi -> [-0.5 0.5]\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    throw(\"Unknown shape function requested.\")\n",
    "end\n",
    "\n",
    "function calculate_element_normals(model)\n",
    "    \"\"\"A solver to calculate boundary element normals.\"\"\"\n",
    "    normals = new_field(model.elno, \"normals\")\n",
    "    coords = model.elno[\"coords\"]\n",
    "    dsf = get_shape_function_derivatives(\"Lagrange\", 1, 1)\n",
    "    for (elid, node_ids) in model.el[\"connectivity\"]\n",
    "        normals[elid] = Dict()\n",
    "\n",
    "        t = interpolate(dsf, coords[elid], -1)\n",
    "        n = [t[2], -t[1]]\n",
    "        n /= norm(n)\n",
    "        normals[elid][1] = n/norm(n)\n",
    "\n",
    "        t = interpolate(dsf, coords[elid], 1)\n",
    "        n = [t[2], -t[1]]\n",
    "        n /= norm(n)\n",
    "        normals[elid][2] = n/norm(n)\n",
    "\n",
    "    end\n",
    "end\n",
    "\n",
    "function test_element_normal_field_calculation()\n",
    "    m = new_model()\n",
    "    dummymeshreader(m) # Read mesh\n",
    "    no2elno(m, \"coords\") # Move nodal field (coordinates) from nodes to element nodes\n",
    "    calculate_element_normals(m)\n",
    "    n = [3, -4]\n",
    "    n /= norm(n)\n",
    "    println(n)\n",
    "    element_normals = m.elno[\"normals\"][1] # Normals of first element\n",
    "    println(element_normals)\n",
    "    @assert element_normals[1] == n\n",
    "    @assert element_normals[2] == n\n",
    "end\n",
    "\n",
    "test_element_normal_field_calculation()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our last solver reads element node values and transforms them back to nodal field values. We now have different values for each element node so we have to somehow combine them. Let's take simple average for now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.7,-0.7]\n",
      "{\"connects_to\",\"normals\",\"coords\"}\n",
      "{2=>[0.7,-0.7],3=>[0.8,-0.6],1=>[0.6,-0.8]}\n"
     ]
    }
   ],
   "source": [
    "function elno2no(model, field)\n",
    "    \"\"\"Move field from element nodes to global nodes.\"\"\"\n",
    "    f = new_field(model.no, field)\n",
    "    for (nid, elids) in model.no[\"connects_to\"]  # loop through nodes \"reverse\" connectivity list\n",
    "        # find variable dimension from first element first local variable\n",
    "        vardim = size(model.elno[field][elids[1]][1])\n",
    "        f[nid] = zeros(vardim)\n",
    "        for elid in elids\n",
    "            local_node_id = findfirst(model.el[\"connectivity\"][elid], nid)\n",
    "            f[nid] += model.elno[field][elid][local_node_id]\n",
    "        end\n",
    "        f[nid] /= length(elids)\n",
    "    end\n",
    "end\n",
    "\n",
    "function test_elno2no()\n",
    "    m = new_model()\n",
    "    dummymeshreader(m)\n",
    "    no2elno(m, \"coords\") # Move nodal field (coordinates) from nodes to element nodes\n",
    "    calculate_element_normals(m) # Calculate normals in elements\n",
    "    #println(m.elno[\"normals\"])\n",
    "    elno2no(m, \"normals\") # Move element node normals back to global nodes\n",
    "    n1 = [3, -4]\n",
    "    n1 /= norm(n1)\n",
    "    n2 = [4, -3]\n",
    "    n2 /= norm(n2)\n",
    "    navg = 0.5*(n1+n2)\n",
    "    println(navg)\n",
    "    println(keys(m.no))\n",
    "    println(m.no[\"normals\"])\n",
    "    @assert m.no[\"normals\"][1] == n1  # Normal of first node\n",
    "    @assert m.no[\"normals\"][2] == navg  # Normal of middle node\n",
    "    @assert m.no[\"normals\"][3] == n2  # Normal of third node\n",
    "end\n",
    "\n",
    "test_elno2no()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So our final solver combines all above to find correct normal directions in boundary nodes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{2=>[0.7071067811865475,-0.7071067811865475],3=>[0.8,-0.6],1=>[0.6,-0.8]}\n"
     ]
    }
   ],
   "source": [
    "function calculate_node_normals(model)\n",
    "    no2elno(model, \"coords\") # Move nodal field (coordinates) from nodes to element nodes\n",
    "    calculate_element_normals(model) # Calculate normals in elements\n",
    "    elno2no(model, \"normals\") # Move element node normals back to global nodes\n",
    "    for (nid, ndir) in model.no[\"normals\"]\n",
    "        model.no[\"normals\"][nid] /= norm(model.no[\"normals\"][nid])\n",
    "    end\n",
    "end\n",
    "\n",
    "function test_calculate_node_normals()\n",
    "    m = new_model()\n",
    "    dummymeshreader(m)\n",
    "    calculate_node_normals(m)\n",
    "    n1 = [3, -4]\n",
    "    n1 /= norm(n1)\n",
    "    n2 = [4, -3]\n",
    "    n2 /= norm(n2)\n",
    "    navg = 0.5*(n1+n2)\n",
    "    navg /= norm(navg)\n",
    "    println(m.no[\"normals\"])\n",
    "    @assert m.no[\"normals\"][1] == n1  # Normal of first node\n",
    "    @assert m.no[\"normals\"][2] == navg  # Normal of middle node\n",
    "    @assert m.no[\"normals\"][3] == n2  # Normal of third node\n",
    "end\n",
    "    \n",
    "test_calculate_node_normals()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary for example 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2 {2=>[7.0,7.0],1=>[4.0,3.0]}\n",
      "1 {2=>[4.0,3.0],1=>[0.0,0.0]}\n",
      "[4.0,3.0]\n",
      "[0.7071067811865475,-0.7071067811865475]\n"
     ]
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAqMAAAIUCAYAAAAuUCWlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XtcVAX+//E3CuIlWs3wq7ZZpqJCGUqaZIGJ4gUb2SyhvAGmUqLgBTBEIW8ImghaCyYmWg6aEiil4g3R9S6o+8AIL5tUmkmrZaLIZX5/+MuN0BI8M5+Zc97Px8M/dkDOZ16c2k8DZ46VwWAwgIiIiIhIQD3pAYiIiIhIu7iMEhEREZEYLqNEREREJIbLKBERERGJ4TJKRERERGK4jBIRERGRGC6jRERERCSGyygRERERieEySkRERERiuIwSERERkZg6LaPXr19HVFQUBgwYgEceeQT16tVDamrqXT/3q6++woABA2BnZ4fmzZtj1KhRKCkpeaChiYiIiEgd6rSMXr58GXPmzMHXX38NZ2dnAICVlVWNz/vuu+/g5uaGc+fOISYmBtOmTcMXX3yBfv36oby8/MEmJyIiIiKLZ12Xv9S6dWv88MMPaNGiBY4dO4bu3bvf9fPmz5+PGzduID8/H3//+98BAD169EC/fv2watUqjB07tu6TExEREZHFq9Mrow0aNECLFi0AAAaD4Z6ft3HjRgwePPjOIgoAHh4ecHBwwPr16+tyaCIiIiJSEaNdwPT999/j8uXLeO6552p8rHv37sjPzzfWoYmIiIjIQhhtGb148SIAoFWrVjU+1qpVK/z3v//l740SERERaVydfmf0fty4cQMAYGtrW+NjDRs2vPM5NjY21T5WUlKCbdu24cknn0SjRo2MNR4RERER1dGNGzfwzTffoH///nj00Ucf6GsZbRn9bZEsKyur8bGbN29W+5zf27ZtG0aMGGGssYiIiIhIIZ988gmGDx/+QF/DaMvobz+e/+3H9b938eJFNG/evMarogDw5JNPArj95Dp37mys8Yzmtddew4YNG6TH0Cz2l8P2sthfDtvLYXs5X331FUaMGHFnb3sQRltGH3vsMdjb2+PIkSM1Pnb48OE770/6R7+9Wtq5c2d069bNWOMZjY2NjUXOrRbsL4ftZbG/HLaXw/bylPiVSqPeDnTo0KHIysrCd999d+exnTt34vTp03j99deNeWgxHTt2lB5B09hfDtvLYn85bC+H7dWhzq+MLlu2DFevXsWFCxcAAJs2bUJxcTEAYNKkSXj44YcRERGBzz77DC+//DKCg4Nx7do1LFy4EF26dIG/v78yz4CIiIiILFadl9H3338f58+fB3D7VqCff/450tPTYWVlhVGjRuHhhx/G3//+d+zZswdTpkzB9OnTYWtri8GDB+P999+/6++LEhEREZG21HkZ/c9//nNfn+fo6IitW7fW9TAWZ/DgwdIjaBr7y2F7Wewvh+3lsL06GPV3RrUoKytLegRNY385bC+L/eWwvRy2VwcuowqLjo6WHkHT2F8O28tifzlsL4ft1YHLqML4FhOy2F8O28tifzlsL4ft1YHLKBERERGJ4TJKRERERGK4jCosJSVFegRNY385bC+L/eWwvRy2VwcuowrLy8uTHkHT2F8O28tifzlsL4ft1cHKYDAYpIf4vby8PLi4uODYsWP8xWQiIiIiM6TkvsZXRomIiIhIDJdRIiIiIhLDZZSIiIiIxHAZVZhOp5MeQdPYXw7by2J/OWwvh+3VgcuowoKCgqRH0DT2l8P2sthfDtvLYXt14NX0RERERFQrvJqeiIiIiFSByygRERERieEyqrCMjAzpETSN/eWwvSz2l8P2ctheHbiMKkyv10uPoGnsL4ftZbG/HLaXw/bqwAuYiIiIiKhWeAETEREREakCl1EiIiIiEsNllIiIiIjEcBlVmL+/v/QImsb+ctheFvvLYXs5bK8OXEYV5unpKT2CprG/HLaXxf5y2F4O26sDr6YnIiIiolrh1fREREREpApcRomIiIhIDJdRhe3bt096BE1jfzlsL4v95bC9HLZXBy6jCouLi5MeQdPYXw7by2J/OWwvh+3VgRcwKay0tBSNGzeWHkOz2F8O28tifzlsL4ft5fACJjPGfyhksb8ctpfF/nLYXg7bqwOXUSIiIiISw2WUiIiIiMRwGVVYaGio9Aiaxv5y2F4W+8thezlsrw5cRhXWpk0b6RE0jf3lsL0s9pfD9nLYXh14NT0RERER1QqvpiciIiIiVeAySkRERERiuIwqrLCwUHoETWN/OWwvi/3lsL0ctlcHLqMKCwsLkx5B09hfDtvLYn85bC+H7dWBy6jCli1bJj2CprG/HLaXxf5y2F4O26sDl1GF8W0mZLG/HLaXxf5y2F4O26sDl1EiIiIiEsNllIiIiIjEcBlVWGxsrPQImsb+ctheFvvLYXs5bK8OXEYVVlpaKj2CprG/HLaXxf5y2F4O26sDbwdKRERERLXC24ESERERkSpwGSUiIiIiMVxGFVZSUiI9gqaxvxy2l8X+ctheDturA5dRhQUEBEiPoGnsL4ftZbG/HLaXw/bqwGVUYdHR0dIjaBr7y2F7Wewvh+3lsL06cBlVGN8BQBb7y2F7Wewvh+3lsL06cBklIiIiIjFcRomIiIhIDJdRhaWkpEiPoGnsL4ftZbG/HLaXw/bqwGVUYXl5edIjaBr7y2F7Wewvh+3lsL068HagRERERFQrvB0oEREREakCl1EiIiIiEmP0ZfTo0aMYMmQIWrdujSZNmqBz586YM2cObty4YexDExEREZGZM+oy+u9//xsvvvgi/v3vfyMkJAQJCQlwdXVFVFQU3njjDWMeWoxOp5MeQdPYXw7by2J/OWwvh+3VwdqYXzwtLQ23bt3CF198gc6dOwMA3nrrLVRVVWH16tX4+eef8be//c2YI5hcUFCQ9Aiaxv5y2F4W+8thezlsrw5GfWW0UaNGAIAWLVpUe7xly5aoX78+GjRoYMzDi/D09JQeQdPYXw7by2J/OWwvh+3VwajLaEBAAP7v//4PY8aMwYkTJ/Dtt99i3bp1SEpKwqRJk+4sq0RERESkTUb9MX3r1q3xr3/9C4MGDULXrl3vPB4ZGYnZs2cb89BEREREZAGM+sropUuXMHDgQBgMBnz00UdIT09HQEAA5s2bhw8++MCYhxaTkZEhPYKmsb8ctpfF/nLYXg7bq4NRl9E5c+bg+++/x+7duzFmzBh4e3tjxYoVGD16NMLDw/Hf//73nn930KBB0Ol01f64urrWOPGys7PvejXdhAkTatyzNi8vDzqdDiUlJdUej4qKQmxsbLXHiouLodPpUFhYWO3xpUuXIjQ0tNpjpaWl0Ol02LdvH/R6/Z3H9Xo9/P39a8zm4+Nj9s/j9yzpeej1elU8D8Dyvh96vV4VzwOwzO/HhAkTVPE8LPH7MWvWLFU8D0v8fqxatUoVz8Pcvx96vf7OLta2bVs4Oztj8uTJNb5OXRn1dqDOzs6ws7PD3r17qz3++eefY+jQodixYwf69OlT7WO8HSgRERGRebOY24GWl5ejsrLyro8DQEVFhTEPT0RERERmzqjLaLdu3ZCXl4fTp09Xe1yv16N+/fro0qWLMQ9PRERERGbOqFfTh4aGYuPGjXjppZcQFBSERx55BFlZWdi6dSvGjh2Lli1bGvPwRERERGTmjPrKaJcuXZCTk4OuXbti4cKFmDx5Mv7zn/9g/vz5+Oc//2nMQ4u52y8Gk+mwvxy2l8X+ctheDturg1FfGQWAHj16YMuWLcY+jNng3SBksb8ctpfF/nLYXg7bq4NRr6avC15NT0RERGTeLOZqeiIiIiKiP8NllIiIiIjEcBlV2B/vdkCmxf5y2F4W+8thezlsrw5cRhUWFxcnPYKmsb8ctpfF/nLYXg7bqwMvYFJYaWkpGjduLD2GZrG/HLaXxf5y2F4O28vhBUxmjP9QyGJ/OWwvi/3lsL0ctlcHLqNEREREJIbLKBERERGJ4TKqsNDQUOkRNI395bC9LPaXw/Zy2F4duIwqrE2bNtIjaBr7y2F7Wewvh+3lsL068Gp6IiIiIqoVXk1PRERERKrAZZSIiIiIxHAZVVhhYaH0CJrG/nLYXhb7y2F7OWyvDlxGFRYWFiY9gqaxvxy2l8X+ctheDturA5dRhS1btkx6BE1jfzlsL4v95bC9HLZXBy6jCuPbTMhifzlsL4v95bC9HLZXBy6jRERERCSGyygRERERieEyqrDY2FjpETSN/eWwvSz2l8P2ctheHbiMKqy0tFR6BE1jfzlsL4v95bC9HLZXB94OlIiIiIhqhbcDJSIiIiJV4DJKRERERGK4jCqspKREegRNY385bC+L/eWwvRy2VwcuowoLCAiQHkHT2F8O28tifzmW2P7y5cvw9/eHk5MTOnbsCCcnJ/j7++Py5cvSo9WKJbanmqylB1Cb6Oho6RE0jf3lsL0s9pdjae1//PFHvPDCCzh79my1x0+dOoW9e/fiwIEDsLe3F5qudiytPd0dXxlVGN8BQBb7y2F7Wewvx9Lah4eH11hEf3P27FmEhYWZeKK6s7T2dHdcRomIiDTk8OHDD/RxIqVxGSUiItKQioqKB/o4kdK4jCosJSVFegRNY385bC+L/eVYWntr6z+/XKR+/fommuTBWVp7ujsuowrLy8uTHkHT2F8O28tifzmW1N5gMKCqqupPP+fKlSs4efKkiSZ6MJbUnu6NtwMlIiLSAIPBgGnTpmHx4sWwt7e/69s4tWrVCg899BDOnj2L8ePHY86cOWjevLnAtGTueDtQIiIiqpWZM2di8eLFWLp0KQoKCuDn5wdHR0c4ODjA0dERfn5+OHHiBAoKCvD+++9j7dq16NChA5YuXcrfIyWj4vuMEhERqdzcuXMxb948LFq0CEFBQQCAjz/++J6fHxISgjfffBORkZEIDg5GcnIylixZgr59+5pqZNIQvjJKRESkYgsXLsTMmTMxd+5cTJ069b7/XosWLbB8+XIcPXoUzZo1Q79+/fCPf/wD586dM+K0pEVcRhWm0+mkR9A09pfD9rLYX445t1+6dCnCwsIQGRmJGTNm1OlrdOvWDbm5udDr9Th69Cg6d+6MiIgI/PrrrwpPW3vm3J7uH5dRhf324w+Swf5y2F4W+8sx1/bLly/HpEmTMG3aNMyePfuBvpaVlRV8fX1RWFiI6dOnIz4+Hh07dsQnn3wCyeugzbU91Q6vpiciIlKZ1NRU+Pv7Y8KECUhMTISVlZWiX/+bb75BWFgYPvvsM7i6uiIhIQHdu3dX9Bhk3ng1PREREd1VWloaAgIC8NZbbyEhIUHxRRQAnnzySaxfvx67d+/G9evX0aNHDwQEBOCHH35Q/FikflxGiYiIVOLzzz/HiBEjMHz4cCQlJaFePeP+33zv3r1x7Ngx/POf/8SmTZvg4OCAhQsX4tatW0Y9LqkLl1GFZWRkSI+gaewvh+1lsb8cc2n/xRdfwMfHB6+99hpWrlxp9EX0N9bW1ggMDERRURH8/Pzw7rvv4umnn0ZWVpbRf5/UXNrTg+EyqjC9Xi89gqaxvxy2l8X+csyh/fbt2zF06FB4eXlhzZo1f3n/eWN45JFHkJiYiBMnTuCJJ57AK6+8gkGDBqGwsNBoxzSH9vTgeAETERGRBduzZw8GDhyIl19+Genp6bC1tZUeCQaDAZmZmZgyZQq+/fZbTJw4EbNmzULTpk2lRyOF8AImIiIiwv79++Hl5YVevXph48aNZrGIArffCsrb2xunTp3C7NmzsXz5cjg4OOCjjz5CZWWl9HhkZriMEhERWaCjR49i4MCBcHFxQUZGBho2bCg9Ug0NGzbEu+++i6+//hoDBgzAuHHj0L17d+zbt096NDIjXEaJiIgszIkTJ+Dp6QlHR0dkZWWhSZMm0iP9qcceewyrV6/G/v37YW1tjZdeeglvvvkmvv32W+nRyAxwGVWYv7+/9Aiaxv5y2F4W+8sxdftTp06hb9++eOqpp7BlyxbY2dmZ9PgPwtXVFQcPHsTHH3+MXbt2oVOnTpgzZw5u3LhRp6/H814duIwqzNPTU3oETWN/OWwvi/3lmLJ9UVERPDw80Lp1a2zbts0iLwiqV68e/Pz8UFRUhAkTJmDOnDno3LkzNmzYUOu3guJ5rw68mp6IiMgCnDt3Dm5ubnj44YeRk5ODFi1aSI+kiKKiIkydOhVZWVno3bs3EhIS0KVLF+mx6C/wanoiIiINKS4uhoeHBxo1aoSdO3eqZhEFAAcHB2zevBlbtmzBxYsX0bVrV7zzzjv46aefpEcjE+EySkREZMYuXLgADw8PAMCuXbvQqlUr4YmMY8CAAfj3v/+NRYsW4dNPP8VHH30kPRKZCJdRhfHtKmSxvxy2l8X+cozZ/scff4SHhwdu3ryJXbt24fHHHzfascyBjY0NJk+ejNOnTyM4OPgvP5/nvTpwGVVYXFyc9Aiaxv5y2F4W+8sxVvuffvoJffv2xdWrV7Fr1y60bdvWKMcxRy1atECjRo3+8vN43quD6W9eq3JpaWnSI2ga+8the1nsL8cY7a9evQpPT0/88MMPyMnJQYcOHRQ/hhrwvFcHLqMKa9y4sfQImsb+ctheFvvLUbr9tWvXMGDAAHzzzTfYvXs3HB0dFf36asLzXh24jBIREZmJ69evw8vLC1999RV27drFtzgiTeAySkREZAZu3LgBnU6H/Px8ZGdnw8XFRXokIpPgBUwKCw0NlR5B09hfDtvLYn85SrQvKyvD0KFDceDAAXzxxRdwdXVVYDL143mvDiZZRvPy8qDT6dC8eXM0adIEzzzzDJYuXWqKQ5tcmzZtpEfQNPaXw/ay2F/Og7YvLy+Hj48Pdu3ahU2bNsHNzU2hydSP5706GP12oNnZ2XjllVfg4uICHx8fPPTQQzhz5gwMBgMWLFhQ4/N5O1AiItKKiooKvPnmm8jIyEBGRgYGDRokPRLRfVFyXzPq74z+8ssvGDVqFF555RVs2LDBmIciIiKyKJWVlfD390d6ejo2bNjARZQ0y6g/pl+7di1+/PFHzJs3D8DtqwSrqqqMeUgiIiKzV1VVhcDAQKxduxaffvopvL29pUciEmPUZXTHjh14+OGH8e2336Jjx46ws7PD3/72N7zzzjsoKysz5qHFFBYWSo+gaewvh+1lsb+c2rY3GAyYNGkSUlJSsGrVKvj4+BhpMvXjea8ORl1GT58+jYqKCnh7e2PgwIFIT09HQEAAkpKS4O/vb8xDiwkLC5MeQdPYXw7by2J/ObVpbzAYMG3aNHzwwQdITk7GyJEjjTiZ+vG8Vwej/s7or7/+itLSUrz99ttYsmQJAMDb2xu3bt1CcnIyZs+ejfbt2xtzBJNbtmyZ9Aiaxv5y2F4W+8upTfuZM2di8eLFWLp0KcaOHWvEqbSB5706GPWV0UaNGgEA3njjjWqP//a/Dx48eM+/O2jQIOh0ump/XF1dkZGRUe3zsrOzodPpavz9CRMmICUlpdpjv73FVElJSbXHo6KiEBsbW+2x4uJi6HS6Gj8CWLp0aY33NSstLYVOp8O+ffuqvc2EXq+/6yvAPj4+Zv88fs+SnkebNm1U8TwAy/t+tGnTRhXPA7DM70dQUJAqnoclfj/0ev19PY/Bgwdj3rx5WLRoEYKCgszueVji9+OPtwO11Odh7t8PvV5/Zxdr27YtnJ2dMXny5Bpfp66M+tZOnp6e2LFjB77++mt06NDhzuOFhYVwdHREQkICJk6cWO3v8K2diIhIbRYuXIiwsDDMnTsXM2bMkB6H6IEpua8Z9ZXR5557DgDw3XffVXv8woULAAB7e3tjHp6IiEjc0qVLERYWhsjISC6iRHdh1GV02LBhAFDjZekVK1bAxsYGvXv3NubhRfzxZXUyLfaXw/ay2F/On7Vfvnw5Jk2ahGnTpmH27NkmnEobeN6rg1EvYHJ2dkZAQABWrlyJiooKuLm5IScnBxs2bEBERARatmxpzMOLKC0tlR5B09hfDtvLYn8592qfmpqKwMBABAUFIS4uDlZWViaeTP143quD0W8HWlFRgfnz5+Pjjz/GhQsX8OSTT2LChAmYNGnSXT+fvzNKRESWLi0tDcOHD8eYMWOQlJSEevWM+oNIIpOzmNuBAoC1tTVmzZqFWbNmGftQRERE4tLT0zFixAgMHz6ciyjRfeA/IURERArJysqCr68vXnvtNaxcuZKLKNF94D8lCvvje4WRabG/HLaXxf5yfmufnZ2NoUOHwsvLC2vWrIG1tdF/+Kh5PO/VgcuowgICAqRH0DT2l8P2sthfTkBAAHJycuDt7Y2+ffsiLS0NNjY20mNpAs97deAyqrDo6GjpETSN/eWwvSz2l/OPf/wDgwcPRq9evbBx40bY2tpKj6QZPO/VgcuowvgOALLYXw7by2J/GUeOHEFISAhcXFyQkZGBhg0bSo+kKTzv1YHLKBERUR0cP34c/fv3h6OjI7KystCkSRPpkYgsEpdRIiKiWiooKEC/fv3w1FNPYcuWLbCzs5MeichicRlV2B9vfUqmxf5y2F4W+5tOUVERPDw80Lp1a2zbtg0bN26UHkmzeN6rA5dRheXl5UmPoGnsL4ftZbG/aZw7dw59+vTBI488gu3bt6N58+ZsL4jt1cHotwOtLd4OlIiIzFFxcTHc3NxgY2OD3NxctGrVSnokIjFK7mt8ZZSIiOgvXLhwAX369IGVlRV27drFRZRIQbw9BBER0Z+4dOkSPDw8UFZWhtzcXDz++OPSIxGpCpdRIiKieygpKUHfvn1x9epV5Obmom3bttIjEakOf0yvMJ1OJz2CprG/HLaXxf7Ku3r1Kjw9PXHp0iXs3LkTHTp0uOvnsb0ctlcHvjKqsKCgIOkRNI395bC9LPZX1i+//IIBAwbg/Pnz2L17NxwdHe/5uWwvh+3VgVfTExER/c7169cxYMAAnDx5Ert27YKLi4v0SERmR8l9ja+MEhER/X83btyATqfD8ePHkZ2dzUWUyAS4jBIREQEoKyvDq6++igMHDmDr1q1wdXWVHolIE3gBk8IyMjKkR9A09pfD9rLY/8GUl5dj2LBh2L17NzZt2gQ3N7f7/rtsL4ft1YHLqML0er30CJrG/nLYXhb7111FRQWGDx+OLVu2ID09HX379q3V32d7OWyvDryAiYiINKuyshKjR49GWloaNmzYAG9vb+mRiCwCL2AiIiJ6QFVVVRg/fjz0ej3Wrl3LRZRICJdRIiLSHIPBgIkTJ2LlypVITU2Fj4+P9EhEmsVllIiINMVgMGDq1Kn48MMPsXz5cowcOVJ6JCJN4wVMCvP395ceQdPYXw7by2L/+2MwGDBjxgzEx8dj6dKlGDt27AN/TbaXw/bqwGVUYZ6entIjaBr7y2F7Wex/f+bOnYuYmBgsWrRIsVtJsr0ctlcHXk1PRESaEBcXh/DwcMydOxczZsyQHofIoim5r/GVUSIiUr3ExESEh4cjMjKSiyiRmeEySkREqpacnIzg4GBMmzYNs2fPlh6HiP6Ay6jC9u3bJz2CprG/HLaXxf53l5qaisDAQAQFBSEuLg5WVlaKH4Pt5bC9OnAZVVhcXJz0CJrG/nLYXhb716TX6xEQEICxY8ciISHBKIsowPaS2F4deAGTwkpLS9G4cWPpMTSL/eWwvSz2ry49PR3Dhg3Dm2++iVWrVqFePeO99sL2ctheDi9gMmP8h0IW+8the1ns/z9ZWVnw9fXFa6+9hpUrVxp1EQXYXhLbqwOXUSIiUo3s7GwMHToUXl5eWLNmDayteaNBInPHZZSIiFQhJycH3t7e6Nu3L9LS0mBjYyM9EhHdBy6jCgsNDZUeQdPYXw7by9J6//3792Pw4MHo1asXNm7cCFtbW5MdW+vtJbG9OnAZVVibNm2kR9A09pfD9rK03P/IkSMYOHAgXFxckJGRgYYNG5r0+FpuL43t1YFX0xMRkcU6fvw4+vTpg44dOyI7Oxt2dnbSIxFpAq+mJyIizSsoKEC/fv3w1FNPYcuWLVxEiSwUl1EiIrI4RUVF8PDwQOvWrbFt2zY0bdpUeiQiqiMuoworLCyUHkHT2F8O28vSUv9z586hT58+eOSRR7B9+3Y0b95cdB4ttTc3bK8OXEYVFhYWJj2CprG/HLaXpZX+xcXF6NOnDxo1aoSdO3eiRYsW0iNppr05Ynt14DKqsGXLlkmPoGnsL4ftZWmh/4ULF9CnTx9YWVlh165daNWqlfRIALTR3lyxvTrw1hQK49tMyGJ/OWwvS+39L126BA8PD5SVlSE3NxePP/649Eh3qL29OWN7deAySkREZq2kpAR9+/bF1atXkZubi7Zt20qPREQK4jJKRERm6+rVq/D09MSlS5eQk5ODDh06SI9ERArj74wqLDY2VnoETWN/OWwvS439f/nlFwwYMADnz5/Hjh074OjoKD3SXamxvaVge3XgK6MKKy0tlR5B09hfDtvLUlv/69evw8vLC1999RV27dqFLl26SI90T2prb0nYXh14O1AiIjIrN27cwODBg3H48GFkZ2fD1dVVeiQi+gMl9zW+MkpERGajrKwMr776Kg4cOICtW7dyESXSAC6jRERkFsrLyzFs2DDs3r0bWVlZcHNzkx6JiEyAFzAprKSkRHoETWN/OWwvy9L7V1RUYPjw4diyZQvS09PRt29f6ZHum6W3t2Rsrw5cRhUWEBAgPYKmsb8ctpdlyf0rKyvh5+eH9PR0rF+/HoMGDZIeqVYsub2lY3t14I/pFRYdHS09gqaxvxy2l2Wp/auqqjB+/Hjo9XqsXbsW3t7e0iPVmqW2VwO2VwcuowrjOwDIYn85bC/LEvsbDAZMnDgRK1euRGpqKnx8fKRHqhNLbK8WbK8OXEaJiMjkDAYDpk6dig8//BDLly/HyJEjpUciIiH8nVEiIjIpg8GAGTNmID4+HkuXLsXYsWOlRyIiQVxGFZaSkiI9gqaxvxy2l2VJ/efOnYuYmBgsWrQIQUFB0uM8MEtqrzZsrw5cRhWWl5cnPYKmsb8ctpdlKf3j4uIwa9YszJ07F1OnTpUeRxGW0l6N2F4deDtQIiIyicTERAQHByPzjQCWAAAgAElEQVQyMhJz5syRHoeIHoCS+5pJXxmdN28e6tWrh2eeecaUhyUiImHJyckIDg7GtGnTMHv2bOlxiMiMmGwZ/e677zB//nw0adIEVlZWpjosEREJS01NRWBgIIKCghAXF8f/DyCiakz21k7Tpk3DCy+8gIqKCt6+i4hII/R6PQICAjB27FgkJCRwESWiGkzyymhubi42btyIJUuWwGAwqPpfRjqdTnoETWN/OWwvyxz7p6enY+TIkRg+fDiSkpJQr546r5k1x/ZawfbqYPR/M1RWVmLixIkYO3YsnJycjH04cWp4mxJLxv5y2F6WufXPysqCr68vXnvtNaxcuVK1iyhgfu21hO3Vweg/pk9KSkJxcTF27dpl7EOZBU9PT+kRNI395bC9LHPqn52djaFDh8LLywtr1qyBtbW6b/ZnTu21hu3Vwaj/qfrTTz9h1qxZmDVrFpo3b27MQxERkRnIycmBt7c3+vbti7S0NNjY2EiPRERmzqjLaGRkJB599FFMnDix1n930KBB0Ol01f64uroiIyOj2udlZ2ff9XdGJkyYUOPODHl5edDpdDUuoIqKikJsbGy1x4qLi6HT6VBYWFjt8aVLlyI0NLTaY6WlpdDpdNi3b1+1x/V6Pfz9/WvM5uPjw+fB58HnweehuucxdOhQ9O/fH7169cLGjRtha2trkc9DLd8PPg8+D6Weh16vv7OLtW3bFs7Ozpg8eXKNr1NXRnvT+9OnT6Nz585YsmQJBg8efOdxX19fXL16Fdu2bcPDDz+MZs2aVft7lv6m9xkZGfD29pYeQ7PYXw7by5Luf+TIEfTt2xfOzs748ssv0aRJE7FZTE26vZaxvRyLeNP777//HlVVVZg0aRKeeuqpO38OHz6MoqIitG3bVpV34NDr9dIjaBr7y2F7WZL9jx8/jv79+8PR0RFZWVmaWkQBnvuS2F4djPbK6E8//YR9+/ZVexsng8GAyMhI/Prrr0hISEC7du1qXGFv6a+MEhFpSUFBAXr37o0nnngCO3bsQNOmTaVHIiITUHJfM9oljs2bN8eQIUNqPB4fHw+A7w1GRGTpioqK4OHhgdatW2Pbtm1cRImoTkz+xm9WVlaqftN7IiItOHfuHPr06YNHHnkE27dv5zumEFGdmfzN33bv3m3qQxIRkYKKi4vRp08fNGrUCDt37kSLFi2kRyIiC6beW2IIudtbJpDpsL8ctpdlqv4XLlxAnz59YGVlhV27dqFVq1YmOa4547kvh+3VQd23xRDAu0HIYn85bC/LFP0vXboEDw8PlJWVITc3F48//rjRj2kJeO7LYXt1MNrV9HXFq+mJiMxPSUkJXn75ZZSUlCA3NxcdOnSQHomIBFnE1fRERKQOV69ehaenJy5duoScnBwuokSkKC6jRER0T7/88gsGDBiA8+fPY/fu3XB0dJQeiYhUhhcwKeyP94El02J/OWwvyxj9r1+/Di8vLxQWFiI7OxtdunRR/BhqwHNfDturA5dRhcXFxUmPoGnsL4ftZSnd/8aNG9DpdDh+/Di2bt0KFxcXRb++mvDcl8P26sAf0yssLS1NegRNY385bC9Lyf5lZWV49dVXcfDgQWzduhU9e/ZU7GurEc99OWyvDlxGFda4cWPpETSN/eWwvSyl+peXl2PYsGHYvXs3vvjiC7z00kuKfF0147kvh+3VgcsoEREBACoqKjB8+HBs2bIFmZmZ8PDwkB6JiDSAyygREaGyshJ+fn74/PPPsWHDBgwcOFB6JCLSCF7ApLDQ0FDpETSN/eWwvawH6V9VVYXx48dDr9fj008/xZAhQxScTP147sthe3XgK6MKa9OmjfQImsb+ctheVl37GwwGTJw4EStXrsTq1asxbNgwhSdTP577ctheHXg7UCIijTIYDJg6dSri4+Px0Ucf4a233pIeiYgshJL7Gn9MT0SkQQaDATNmzEB8fDyWLVvGRZSIxHAZJSLSoLlz5yImJgbvv/8+JkyYID0OEWkYl1GFFRYWSo+gaewvh+1l1aZ/XFwcZs2ahXnz5mHKlClGnEobeO7LYXt14DKqsLCwMOkRNI395bC9rPvtn5iYiPDwcMycORMRERFGnkobeO7LYXt14DKqsGXLlkmPoGnsL4ftZd1P/+TkZAQHByM0NBTvvfeeCabSBp77ctheHbiMKoxvMyGL/eWwvay/6p+amorAwEBMnDgRsbGxsLKyMtFk6sdzXw7bqwOXUSIildPr9QgICMC4ceOQkJDARZSIzAqXUSIiFUtPT8fIkSMxcuRI/POf/+QiSkRmh8uowmJjY6VH0DT2l8P2su7WPysrC76+vnj99deRkpKCevX4r3xj4Lkvh+3Vgf9mUlhpaan0CJrG/nLYXtYf+2dnZ2Po0KEYPHgwVq9ejfr16wtNpn489+WwvTrwdqBERCqTk5ODQYMGoU+fPkhPT0eDBg2kRyIileHtQImI6K7279+PwYMHo1evXtiwYQMXUSIye1xGiYgs1FdffYUOHTrA1tYWDRo0QIMGDfDiiy+iU6dOyMzMRMOGDaVHJCL6S1xGFVZSUiI9gqaxvxy2N62CggJ06dIFZ86cwa1bt1BeXo7y8nIYDAacOHEC58+flx5RM3juy2F7deAyqrCAgADpETSN/eWwvWl5e3ujoqICAPDHN2uqqKiATqcz/VAaxXNfDturA5dRhUVHR0uPoGnsL4ftTau4uBgA0AbAcQDP3uPjZHw89+WwvTpwGVUY3wFAFvvLYXvTqqqqAgCUAbgC4CiAkN993MzeKEXVeO7LYXt14DJKRGRBrly5gpCQkDs/or8EwAPAYgDxANYDeAjgnZaIyGJwGSUisgCVlZVITk6Gg4MDUlJS0Lx58/99DEA4gFcB9AdwBECHW7eQmJiI8vJymYGJiO4Tl1GFpaSkSI+gaewvh+2NJzc3Fy4uLggMDISXlxeKioqwd+9eWFtbV/u8zwE8B6AcwNH69XEwOBjOzs7Yvn27xNiawXNfDturA5dRheXl5UmPoGnsL4ftlXf+/Hn4+PjA3d0dtra2OHjwIFatWoVWrVqhc+fOOHnyJNq3b48GDRqgXr16aNCgAW498QT+0bIltjVujLUAZl+5Ai9PT3h7e+Ps2bPST0mVeO7LYXt14O1AiYjMTGlpKeLi4hAbG4umTZsiNjYWI0aMQL169/f6wdmzZ+Hu5oZxVVWY+dNP+OnJJzHo119x4qefMGXKFERERMDOzs7Iz4KI1Iy3AyUiUiGDwYB169ahU6dOiImJQUhICIqKijBq1Kj7XkQBoF27dti5axc+NBgw8okn0Ky0FIfKy7HcxwdLlixBx44dsXr16jtX5BMRSeIySkRkBo4fP47evXvD19cX3bp1Q0FBAWJiYur8CmbHjh2xc+dObLt6Ff0ffRQVzzyD0Z9+iosTJ8LtxRcxevRovPDCCzh8+LDCz4SIqHa4jBIRCbp8+TICAwPh4uKCy5cvIzs7GxkZGWjfvv0Df20nJyds374dx86fx8s3bqBs2jQ0XbgQaTdv4l9ZWbh58yaef/55+Pv74+LFiwo8GyKi2uMyqjDegk8W+8th+9opLy9HQkICHBwckJaWhsWLF+PEiRPo169fnb7evfo7OzsjOzsbJwsK0P/QIdzcuBHYtw8vTJqEYykpSEpKwubNm+Hg4IC4uDiUlZU9yNPSJJ77ctheHbiMKiwoKEh6BE1jfzlsf/+ys7Px7LPPYvLkyfDx8cHp06cRHBwMGxubOn/NP+vfvXt3bNmyBUePHoUuKQll+/cDTZui/osvYnyDBjh9+jTGjBmDiIgIPP3009i8eTPv4FQLPPflsL06cBlVmKenp/QImsb+ctj+r505cwZDhgxB//79YW9vj7y8PCQlJcHe3v6Bv/Zf9e/VqxeysrKwd+9eDJ02Dbd27wZGjAACAtAsLAxLFizAyZMn0bZtW+h0OgwcOBBfffXVA8+lBTz35bC9OnAZJSIysmvXrmH69OlwcnJCfn4+1q1bh5ycHDg7O5t0jt69eyMzMxPbt2+Hr58fyj/8EEhJAT75BOjVC46NGmHbtm3IzMzEmTNn0KVLF0yePBlXr1416ZxEpC1cRomIjKSqqgqpqalwcHBAYmIiIiIiUFhYiGHDhondO97T0xMbN27E5s2bMWrUKFSOHg3s3w9cvQq4uMBqyxbodDoUFBRg7ty5WLFiBTp06IDly5ejsrJSZGYiUjcuowrLyMiQHkHT2F8O21d36NAhuLq6ws/PD+7u7igsLERUVBQaN25slOPVpv/gwYORlpaGzz77DGPGjEHVs88CR48CvXoBXl7ArFmwtbZGeHg4vv76awwaNAjjx4/Hc889h7179xplfkvGc18O26sDl1GF6fV66RE0jf3lsP1tFy9exOjRo9GzZ0/cunULe/bsQVpaGtq0aWPU49a2/9ChQ7FmzRqsXr0ab7/9NgxNmwKZmcC8ebf/DBwIlJSgdevWSE1NxcGDB9GgQQO4ubnB19cXxcXFRnomlofnvhy2VwmDmTl27JgBgOHYsWPSoxAR3bebN28aYmJiDA899JDh0UcfNSQnJxsqKiqkx/pLH3/8sQGAYdKkSYaqqqrbD27fbjA8+qjB8PjjBsPBg3c+t7Ky0pCammpo2bKloVGjRobo6GhDaWmp0OREJEnJfY2vjBIRPQCDwYBNmzbByckJM2fOxFtvvYWioiKMGzcO9evXlx7vL/n5+SEpKQmJiYkIDw+//ZZOffsC+fnAY48BL70EfPghYDCgXr16GDVqFIqKijBp0iTMnz8fnTp1wmeffca3giKiOuMySkRUR6dOncKAAQMwZMgQtGvXDidPnkR8fDyaNWsmPVqtjB8/HkuWLMHChQsRHR19+8G//x3YswcIDAQmTABGjQKuXwcA2NnZYcGCBSgoKICzszOGDRuGl19+GSdOnJB7EkRksbiMEhHV0pUrVxASEoIuXbrg7NmzyMzMxNatW9G5c2fp0eosODgYsbGxmD17NubPn3/7wQYNgMREYO1aID0d6NkTKCq683fat2+PzMxMbNu2DT/++CO6deuGt99+GyUlJULPgogsEZdRhfn7+0uPoGnsL0cL7SsrK5GcnAwHBwekpKRg3rx5KCgogE6nE3urpt8o0T8sLAzvvfceZsyYgcWLF//vA2+8ARw+DNy6BTz33O3F9Hc8PT1x4sQJLF68GHq9Hh06dEBiYiLKy8sfeCZLoIVz31yxvTpwGVUY7wYhi/3lqL19bm4uXFxcEBgYCC8vLxQVFSE8PBy2trbSowFQrv/MmTPx7rvvYurUqfjggw/+9wEnJ+DIEaB/f2DoUCAsDKiouPNhGxsbBAcH4/Tp0/Dx8UFISAicnZ2xf/9+ReYyZ2o/980Z26sDl1GFvfHGG9IjaBr7y1Fr+/Pnz8PHxwfu7u6wtbXFwYMHsWrVKrRq1Up6tGqU6m9lZYV58+Zh8uTJCAoKQkpKyv8++PDDwPr1wOLFt/94eAA//FDt79vb2yMpKQl5eXl49NFHNXFhk1rPfUvA9upgLT0AEZE5Ki0tRVxcHGJjY9G0aVOkpqZixIgRqFdP/f8Nb2Vlhffffx83b97E2LFjYWtrixEjRvz2QWDy5Ns/rh82DOja9faC+tJL1b6Gs7MzcnJyxH99gYjMn/r/rUpEVAsGgwHr1q1Dp06dEBMTg5CQEBQVFWHUqFGaWER/Y2VlhWXLlsHf3x+jR4/GZ599Vv0TXnrp9ts/OTgAL78MxMcDf3gVlIsoEd0P7fyb1UT27dsnPYKmsb8cNbTPz8+Hu7s7fH190a1bNxQUFCAmJgZ2dnbSo/0lY/SvV68eli9fjjfeeANvvvkmMjMzq39Cy5bAzp3AlCm3/wwbBvzyi+JzmDs1nPuWiu3VgcuowuLi4qRH0DT2l2PJ7S9fvozx48fDxcUFJSUlyM7ORkZGBtq3by892n0zVv/69etj1apV8Pb2xuuvv44tW7ZU/wRrayAuDti4Edi2DejRAygoMMos5sqSz31Lx/bqwGVUYWlpadIjaBr7y7HE9uXl5UhISICDgwPWrVuH+Ph4nDhxAv369ZMerdaM2d/a2hpr167FwIED8eqrr2Lnzp01P+nVV4GjRwEbm9sLqYbuGW6J575asL06cBlVWOPGjaVH0DT2l2Np7bOzs/Hss89i8uTJ8PHxwenTpxEcHAwbGxvp0erE2P1tbGywfv169O7dGzqdDnv37q35SQ4OwMGDtxfTN98EJk68/d6kKmdp576asL06GHUZPXLkCIKCguDk5ISHHnoITzzxxJ1/6RMRSThz5gyGDBmC/v37w97eHnl5eUhKSoK9vb30aGbP1tYW6enpeP755xEQEICK373P6B1NmgCrV9++n31yMuDuDnz7remHJSKLYdRlNDY2Fp9//jn69euHxMREjBs3Drm5uXcuDCAiMpVr165h+vTpcHJyQn5+PtatW4ecnBw4OztLj2ZRGjVqhM2bN2PLli2wtr7HuwNaWQFvvw3s2wd8/z3QrdvtC52IiO7CqMvo1KlTcf78eSxZsgQBAQGYMWMG9u7di4qKCixYsMCYhxYTGhoqPYKmsb8cc21fVVWF1NRUODg4IDExERERESgsLMSwYcNU9dZDpuzfpEmT+7u4q0cPIC/v9nuRenoC8+cDVVXGH9DEzPXc1wK2VwejLqOurq41/su5ffv2cHR0RGFhoTEPLaZNmzbSI2ga+8sxx/aHDh2Cq6sr/Pz84O7ujsLCQkRFRany98zMsT8A4NFHgS1bgMhIYMYMYMgQ4MoV6akUZbbtNYDt1cHKYOJ7tRkMBjz++ON45plnar5FCIC8vDy4uLjg2LFj6NatmylHIyKVuHjxIqZPn47Vq1fD2dkZCQkJcHNzkx6LvvwSGDECaNr09ltBde0qPRER1ZGS+5rJr6b/9NNPceHCBfj4+Jj60ESkcmVlZViwYAEcHBzw5ZdfIjk5GUePHuUiai4GDQKOHQOaNQNeeAH4+GPpiYjIDJh0GS0sLMSECRPwwgsvYPTo0aY8NBGpmMFgwKZNm+Dk5ISZM2firbfeQlFREcaNG4f69etLj0e/17Yt8K9/ASNHAgEBwNixwM2b0lMRkSCTLaM//PADvLy80KxZM2zYsEFVFw78nlp/F9ZSsL8cqfanTp3CgAEDMGTIELRr1w4nT55EfHw8mjVrJjKPFIs69xs2BJYvB1auBD75BOjVC/jPf6SnqjOLaq8ybK8OJllGf/75ZwwcOBC//PILtm7dipYtW/7l3xk0aBB0Ol21P66ursjIyKj2ednZ2dDpdDX+/oQJE5CSklLtsby8POh0OpSUlFR7PCoqCrGxsdUeKy4uhk6nq3GiL126tMbVe6WlpdDpdNi3bx/CwsLuPK7X6+Hv719jNh8fH7N/Hr9nSc8jLCxMFc8DsLzvR1hYmEmfx5UrV/Dss8/i6aefxtmzZ5GZmYmtW7fixo0bmvx+9O7d2/Kex48/Avv3A1evAi4uKF61yiK/H6+++mq1x9R0Xpn78wgODlbF8zD374der7+zi7Vt2xbOzs6YPHlyja9TV0a/gOnmzZvw9PREfn4+duzYgeeff/5PP9/SL2AqLi7m1X2C2F+OqdpXVlZixYoViIyMxM2bNxEZGYmQkBDY2toa/djmzKLP/StXgNGjgc2bgZkzgagowIJ+vcKi21s4tpdjMRcwVVZWwsfHB4cOHcJnn332l4uoGvAfClnsL8cU7XNzc+Hi4oLAwEB4eXmhqKgI4eHhml9EAQs/95s1AzIygHnzbv8ZOBD4wytF5syi21s4tleHe9w+QxlTp07F5s2b8corr6CkpASffPJJtY+PGDHCmIcnIpU4f/48wsLCsH79evTo0QMHDx7UxH/cakq9ekBExO03yn/jjdt3bfrsM4DfZyLVM+oyeuLECVhZWWHz5s3YvHlztY9ZWVlxGSWiP1VaWoq4uDjExsaiadOmSE1NxYgRI1CvnsnflY5MpW9fID8feP114KWXgCVLbt9aVKUXvRKRkX9Mv3v3blRWVqKqqqrGn8rKSmMeWswff+GYTIv95SjZ3mAwYN26dejUqRNiYmIQEhKCoqIijBo1iovoPajq3P/734E9e4DAQGDCBGDUKOD6demp7klV7S0M26sD/62usNLSUukRNI395SjVPj8/H+7u7vD19UW3bt1QUFCAmJgY2NnZKfL11Up1536DBkBiIrB2LZCeDvTsCRQVSU91V6prb0HYXh1MfjvQv2LpV9MTUd1cvnwZkZGR+Oijj9CpUyckJCSgX79+0mOROSgoAIYOBS5cAFatAv7wVkpEZHoWczU9EdFfKS8vR0JCAhwcHLBu3TrEx8fjxIkTXETpf5ycgMOHgf79by+loaFARYX0VESkEC6jRCQmOzsbzz77LCZPngwfHx+cPn0awcHBsLGxkR6NzM3DDwPr1wOLFwPx8YCHB/DDD9JTEZECuIwq7I93USDTYn85tWl/5swZDBkyBP3794e9vT3y8vKQlJQEe3t7I06obpo4962sgMmTgd27gdOnga5dgb17pafSRnszxfbqwGVUYQEBAdIjaBr7y7mf9teuXcP06dPh5OSE/Px8rFu3Djk5OXB2djbBhOqmqXP/pZeAvDygY0fg5Zdvv1oqePmDptqbGbZXBy6jCouOjpYeQdPYX86fta+qqkJqaiocHByQmJiIiIgIFBYWYtiwYbDi+0cqQnPnfsuWwI4dwJQpwNSpwLBhwC+/iIyiufZmhO3VgcuowvgOALLYX8692h86dAiurq7w8/ODu7s7CgsLERUVhcaNG5t4QnXT5LlvbQ3ExQEbNwLbtt2+e1NBgcnH0GR7M8H26sBllIiM4uLFixg9ejR69uyJW7duYc+ePUhLS+O9pEl5r74KHD0K2NjcXkjXrpWeiIhqgcsoESmqrKwMCxYsgIODA7788kskJyfj6NGjcHNzkx6N1MzBATh48PZieviw9DREVAtcRhWWkpIiPYKmsb+cFStWIDMzE05OToiMjMSYMWNQVFSEcePGoX79+tLjqR7PfQBNmgCrVwPvv2/Sw7K9HLZXBy6jCsvLy5MeQdPYX8apU6fw3nvvwdvbG+3atcPJkyexZMkSNGvWTHo0zeC5//9ZWQEm/o8ftpfD9urA24ESUZ1duXIF7733HpYtW4Ynn3wSixcvxiuvvMIr5ImIVE7Jfc1aoZmISEMqKyuxYsUKREZG4ubNm5g3bx5CQkJga2srPRoREVkY/pieiGolNzcXLi4uCAwMhJeXF4qKihAeHs5FlIiI6oTLKBHdl/Pnz8PHxwfu7u6wtbXFwYMHsWrVKrRq1Up6NCIismBcRhWm0+mkR9A09ldeaWkpoqOj0alTJ+Tm5iI1NRUHDhzA888/X+3z2F4W+8thezlsrw78nVGFBQUFSY+gaeyvHIPBgPXr1yM0NBSXLl3ClClTEBERATs7u7t+PtvLYn85bC+H7dWBr4wqzNPTU3oETWN/ZeTn58Pd3R2+vr7o1q0bCgoKEBMTc89FFGB7aewvh+3lsL06cBklojsuX76M8ePHw8XFBSUlJcjOzkZGRgbat28vPRoREakUf0xPRCgvL8cHH3yA6OhoAEB8fDzeeecd2NjYyA5GRESqx1dGFZaRkSE9gqaxf+1lZ2fj2WefxZQpU+Dr64vTp08jODi41oso28tifzlsL4ft1YHLqML0er30CJrG/vfvzJkzGDJkCPr37w97e3vk5eUhKSkJ9vb2dfp6bC+L/eWwvRy2VwfeDpRIY65du4Z58+YhPj4e//d//4dFixbh9ddf5y08iYjovvF2oERUa1VVVVizZg2mT5+On3/+GREREQgNDUXjxo2lRyMiIg3jj+mJNODQoUNwdXWFn58f3N3dUVhYiKioKC6iREQkjssokYpdvHgRo0ePRs+ePXHr1i3s2bMHaWlpaNOmjfRoREREALiMKs7f3196BE1j/9vKysqwYMECODg44Msvv0RycjKOHj0KNzc3ox2T7WWxvxy2l8P26sDfGVUY7wYhS+v9DQYDNm3ahKlTp+Kbb75BUFAQoqKi0KxZM6MfW+vtpbG/HLaXw/bqwKvpiVTi1KlTCAkJwfbt2+Hp6Yn4+Hg4OjpKj0VERCqk5L7GH9MTWbgrV64gJCQEXbp0wblz55CZmYmtW7dyESUiIovAH9MTWajKykqsWLECkZGRuHnzJubNm4eQkBDY2tpKj0ZERHTf+Mqowvbt2yc9gqZppX9ubi5cXFwQGBgILy8vFBUVITw8XHQR1Up7c8X+ctheDturA5dRhcXFxUmPoGlq73/+/Hn4+PjA3d0dtra2OHjwIFatWoVWrVpJj6b69uaO/eWwvRy2VwcuowpLS0uTHkHT1Nq/tLQU0dHR6NSpE3Jzc5GamooDBw7g+eeflx7tDrW2txTsL4ft5bC9OvB3RhXGO9rIUlt/g8GA9evXIzQ0FJcuXcKUKVMQEREBOzs76dFqUFt7S8P+ctheDturA18ZJTJT+fn5cHd3h6+vL7p27YqCggLExMSY5SJKRERUV1xGiczM5cuXMX78eLi4uKCkpATbtm1DZmYm2rdvLz0aERGR4riMKiw0NFR6BE2z5P7l5eVYsmQJOnTogHXr1iE+Ph4nTpywmDuMWHJ7NWB/OWwvh+3Vgb8zqrA2bdpIj6Bplto/OzsbISEhKCwsxLhx4zBnzhzY29tLj1UrltpeLdhfDtvLYXt14O1AiQSdOXMGU6dOxaZNm+Dm5oaEhAQ4OztLj0VERPSneDtQIgt37do1TJ8+HU5OTsjPz8e6deuQk5PDRZSIiDSHP6YnMqGqqiqsWbMG06dPx88//4yIiAiEhoby7UmIiEiz+MqowgoLC6VH0DRz7n/o0CG4urrCz88P7u7uKCwsRFRUlGoWUXNurwXsL4ft5bC9OnAZVVhYWJj0CJpmjisdVpoAABukSURBVP0vXryI0aNHo2fPnrh16xb27NmDtLQ01f3ivTm21xL2l8P2ctheHbiMKmzZsmXSI2iaOfUvKyvDggUL4ODggC+//BLJyck4evQo3NzcpEczCnNqr0XsL4ft5bC9OvB3RhWmtle7LI059DcYDNi0aROmTp2Kb775BkFBQYiKikKzZs2kRzMqc2ivZewvh+3lsL068JVRIgWdOnUK/fv3h7e3N9q1a4eTJ09iyZIlql9EiYiI6orLKJECrly5guDgYHTp0gXnzp1DZmYmtm7dCkdHR+nRiIiIzBqXUYXFxsZKj6Bppu5fWVmJ5ORkODg4YOXKlZg3bx4KCgqg0+lgZWVl0lmk8dyXxf5y2F4O26sDl1GFlZaWSo+gaabsn5ubCxcXFwQGBsLLywtFRUUIDw+Hra2tyWYwJzz3ZbG/HLaXw/bqwNuBEtXS+fPnERYWhvXr16NHjx5ITEzE888/Lz0WERGRyfB2oEQCSktLER0djU6dOiE3Nxepqak4cOAAF1EiIqIHwLd2IvoLBoMB69evR2hoKC5duoQpU6YgIiICdnZ20qMRERFZPL4yqrCSkhLpETRN6f75+flwd3eHr68vunbtioKCAsTExHARvQue+7LYXw7by2F7deAyqrCAgADpETRNqf6XL1/G+PHj4eLigpKSEmzbtg2ZmZlo3769Il9fjXjuy2J/OWwvh+3VgT+mV1h0dLT0CJr2oP3Ly8vxwQcf3Pk68fHxeOedd2BjY/Pgw6kcz31Z7C+H7eWwvTpwGVUY3wFA1oP0z87ORkhICAoLCzFu3DjMmTMH9vb2Ck6nbjz3ZbG/HLaXw/bqwB/Tk+adOXMGQ4YMQf/+/WFvb4+8vDwkJSVxESUiIjIBLqOkWdeuXcP06dPh5OSE/Px8rFu3Djk5OXB2dpYejYiISDOMvoyWlZUhPDwcrVu3RuPGjdGzZ0/s2LHD2IcVk5KSIj2Cpt1P/6qqKqSmpsLBwQGJiYmIiIhAYWEhhg0bprlbeCqJ574s9pfD9nLYXh2Mvoz6+fkhPj4eI0eORGJiIurXr49BgwbhX//6l7EPLSIvL096BE37q/6HDh2Cq6sr/Pz84O7ujsLCQkRFRaFx48YmmlC9eO7LYn85bC+H7dXBqLcDPXz4MHr27IlFixZhypQpAG6/Uvr000+jRYsWd11IeTtQMoaLFy9i+vTpWL16NZydnZGQkAA3NzfpsYiIiCySxdwOdMOGDbC2tsa4cePuPGZra4sxY8bgwIED+P777415eCKUlZVhwYIFcHBwwJdffonk5GQcPXqUiygREZGZMOpbO+Xn58PBwQEPPfRQtce7d+8OADh+/Dgee+wxY45AGnD58mWEhYXh8OHDqKiogLW1Nbp3744+ffpg9uzZ+OabbxAUFISoqCg0a9ZMelwiIiL6HaMuoxcvXkSrVq1qPP7bYxcuXDDm4UkDfvzxR7zwwgs4e/ZstcdPnTqF1NRU9O7dG5s2bYKjo6PQhERERPRnjPpj+hs3bsDW1rbG4w0bNrzzcbXR6XTSI2hKeHh4jUX095544gkuoibCc18W+8thezlsrw5GXUYbNWqEsrKyGo/fvHnzzsfVJigoSHoETTl8+PCffvzIkSMmmoR47stifzlsL4ft1cGoy2irVq3u+qP4ixcvAgBat259z787aNAg6HS6an9cXV2RkZFR7fOys7Pv+l9GEyZMqPH+Y3l5edDpdCgpKan2eFRUFGJjY6s9VlxcDJ1Oh8LCwmqPL126FKGhodUeKy0thU6nw759++Dp6Xnncb1eD39//xqz+fj4mP3z+D1zfh6XL1+ucazfq6iosIjnoYbvh6enpyqeB2CZ349ly5ap4nlY4vcjPz9fFc/DEr8ff7yK21Kfh7l/P/R6/Z1drG3btnB2dsbkyZNrfJ26MupbO4WFhSE+Ph7//e9/YWdnd+fx+fPnIzIyEt9++22NC5j41k5UG05OTjh16tQ9P+7o6IiCggITTkRERKR+FvPWTq+99hoqKyuxfPnyO4+VlZXh448/Rs+ePXklPT2wHj16PNDHiYiISJZRl9EePXrg9ddfx7vvvovw8HAsX74cffr0QXFxMeLi4ox5aDF/fJmbjCsuLg7t2rW758fbt29vwmm0jee+LPaXw/Zy2F4djH470NWrVyMkJARr1qxBcHAwKisrkZWVhRdffNHYhxah1+ulR9AUe3t7HDhwAH5+fnB0dMRDDz0ER0dHjB49GkFBQYiMjERSUpL0mJrAc18W+8thezlsrw5G/Z3R/9fevQdFdd9vHH8WbwSEFNCimFEoSrRGJladmERbGqy31MVEkVqtaTB2mmqo/mGcyTi5Wac1vQyNTEcTNdqpoW10YFKGRApKiDXehhBTEyLJNAYRSbAhQRFEPL8/WviFYCPqLh895/2a4Y98WfY8++S4+4Hdc87V4DOjCBTHcbR8+XI9++yz2rJlyyU/tA0AAK5cIOe1oJ70HrDk8/mUnZ2tlpYWLV68WP369dMPf/hD61gAAOALGEbhaj6fT3/4wx/U0tKiRYsWqV+/fpozZ451LAAA8F9B/8woYC0kJESbNm1Senq6fvCDH+hvf/ubdSQAAPBfDKMBxucSbf2v/nv16qU//vGP8vv9mjt3roqKino4mfux79uifzt0b4fu3YFhNMC+eAUm9Lyv6r9Pnz7Kzc3V9773PaWlpam0tLTngnkA+74t+rdD93bo3h04mh6e09zcLL/fr3379mnXrl26++67rSMBAHBDuWGuwARcj0JDQ5Wfn69x48ZpxowZOnjwoHUkAAA8i2EUnhQWFqaCggLddtttmjZtmioqKqwjAQDgSQyjAbZ3717rCJ52Jf1HRETolVde0fDhwzVlyhT985//DGIy92Pft0X/dujeDt27A8NogD3zzDPWETztSvu/+eabtWvXLg0ZMkRTpkzRe++9F6Rk7se+b4v+7dC9Hbp3Bw5gCrCmpiaFhYVZx/Csq+3/k08+UUpKihoaGlRWVqbExMQgpHM39n1b9G+H7u3QvR0OYLqO8Y/C1tX2P3DgQBUXFys8PFz33HOPjh8/HuBk7se+b4v+7dC9Hbp3B4ZR4L8GDx6s3bt3q1evXkpNTVVNTY11JAAAXI9hFPiCW265Rbt379b58+eVmpqquro660gAALgaw2iArVy50jqCpwWi//j4eJWUlOjzzz/XlClTVF9fH4Bk7se+b4v+7dC9Hbp3B4bRABs6dKh1BE8LVP8jRoxQSUmJ6urqNHXqVH366acBuV83Y9+3Rf926N4O3bsDR9MDX+HIkSP67ne/q+HDh+vvf/+7IiMjrSMBAGCOo+mBHpKcnKyioiK99957uvfee3X27FnrSAAAuArDKHAZ48aN06uvvqqKigr5/X6dO3fOOhIAAK7BMBpglZWV1hE8LVj9T5w4UYWFhdq/f7/uv/9+tbS0BGU7NzL2fVv0b4fu7dC9OzCMBtijjz5qHcHTgtn/5MmT9fLLL2vPnj2aN2+eWltbg7atGxH7vi36t0P3dujeHRhGAywnJ8c6gqcFu//U1FTl5eXplVde0YIFC3ThwoWgbu9Gwr5vi/7t0L0duncHhtEA4zQTtnqi/xkzZuill15SXl6efvzjH6utrS3o27wRsO/bon87dG+H7t2BYRS4Cmlpadq+fbtKS0u5bCgAANeAYRS4SvPmzdOxY8f4zRwAgGvAMBpg69ats47gaT3df1hYWI9u73rGvm+L/u3QvR26dweG0QBramqyjuBp9G+H7m3Rvx26t0P37sDlQAEAAHBFuBwoAAAAXIFhFAAAAGYYRgOsvr7eOoKn0b8durdF/3bo3g7duwPDaIBlZmZaR/A0+rdD97bo3w7d26F7d2AYDbAnn3zSOoKn0b8durdF/3bo3g7duwPDaIBxBgBb9G+H7m3Rvx26t0P37sAwCgAAADMMowAAADDDMBpgmzdvto7gafRvh+5t0b8durdD9+7AMBpg5eXl1hE8jf7t0L0t+rdD93bo3h24HCgAAACuCJcDBQAAgCswjAIAAMAMwygAAADMMIwGmN/vt47gafRvh+5t0b8durdD9+7AMBpgy5Yts47gafRvh+5t0b8durdD9+7A0fQAAAC4IhxNDwAAAFdgGAUAAIAZhtEAy8/Pt47gafRvh+5t0b8durdD9+7AMBpgubm51hE8jf7t0L0t+rdD93bo3h04gAkAAABXhAOYAAAA4AoMowAAADDDMAoAAAAzDKMB9uCDD1pH8DT6t0P3tujfDt3boXt3YBgNsKlTp1pH8DT6t0P3tujfDt3boXt34Gh6AAAAXBGOpgcAAIArMIwCAADADMNogO3du9c6gqfRvx26t0X/dujeDt27A8NogD3zzDPWETyN/u3QvS36t0P3dujeHYI6jJaUlCgzM1NJSUkKDw9XYmKilixZolOnTgVzs6b+/Oc/W0fwNPq3Q/e26N8O3duhe3foHcw7X7VqlRoaGpSenq4RI0bogw8+UE5OjgoKClRRUaHY2Nhgbt5EWFiYdQRPo387dG+L/u3QvR26d4egDqPZ2dmaNGlSp7Xp06frO9/5jnJycrRmzZpgbh4AAADXuaC+Tf/lQVSSJk+erOjoaFVWVgZz0wAAALgB9PgBTGfOnFFjY6MGDBjQ05vuEStXrrSO4Gn0b4fubdG/Hbq3Q/fu0OPDaHZ2tlpbW5WRkdHTm+4RQ4cOtY7gafRvh+5t0b8durdD9+7Q7cuBOo6jlpaWbt1paGjoJdfLysqUmpqquXPnKjc395K34XKgAAAA1zeTy4G+9tprCgsL69bXsWPHuvx8ZWWl7rvvPiUnJ2vTpk3XFBoAAADu0O1hdNSoUdq6dWu3vgYNGtTpZ6urqzV16lRFRUWpsLBQ4eHhl93ezJkz5ff7O33deeedys/P73S7oqIi+f3+Lj+/dOlSbd68udNaeXm5/H6/6uvrO60/8cQTWrduXae1jz76SH6/v8uBVuvXr+/yGZWmpib5/f4uV4LIzc3Vgw8+2CVbRkYGj4PHwePgcfA4eBw8Dh7HDfE4cnNzO2axhIQE3X777VqxYkWX+7la3X6b/mqdPn1akyZNUkNDg/bu3avExMSvvP2N/jZ9ZWWlRo4caR3Ds+jfDt3bon87dG+H7u2YvE1/Nc6ePauZM2eqtrZWhYWFlx1E3eDRRx+1juBp9G+H7m3Rvx26t0P37hDUk94vWLBAhw4dUmZmpo4ePaqjR492fC8iIkJpaWnB3LyJnJwc6wieRv926N4W/duhezt07w5BHUbfeust+Xw+bdmyRVu2bOn0vfj4eFcOo5xmwhb926F7W/Rvh+7t0L07BHUY/de//hXMuwcAAMANrsdPeg8AAAC0YxgNsC+figE9i/7t0L0t+rdD93bo3h0YRgOsqanJOoKn0b8durdF/3bo3g7du0PQzzN6pW7084wCAAC43Q1znlEAAADgqzCMAgAAwAzDaIB9+fqy6Fn0b4fubdG/Hbq3Q/fuwDAaYJmZmdYRPI3+7dC9Lfq3Q/d26N4dGEYD7Mknn7SO4Gn0b4fubdG/Hbq3Q/fuwDAaYJwBwBb926F7W/Rvh+7t0L07MIwCAADADMMoAAAAzDCMBtjmzZutI3ga/duhe1v0b4fu7dC9OzCMBlh5ebl1BE+jfzt0b4v+7dC9Hbp3By4HCgAAgCvC5UABAADgCgyjAAAAMMMwCgAAADMMowHm9/utI3ga/duhe1v0b4fu7dC9OzCMBtiyZcusI3ga/duhe1v0b4fu7dC9O3A0PQAAAK4IR9MDAADAFRhGAQAAYIZhNMDy8/OtI3ga/duhe1v0b4fu7dC9OzCMBti6deusI3ga/duhe1v0b4fu7dC9OzCMBtjAgQOtI3ga/duhe1v0b4fu7dC9OzCMAgAAwAzDKAAAAMwwjAIAAMBMb+sAX3bu3DlJ0rvvvmuc5OocPHhQ5eXl1jE8i/7t0L0t+rdD93bo3k77nNY+t12L6+4KTNu3b9fChQutYwAAAOAy/vSnP2nBggXXdB/X3TBaX1+vXbt2KT4+XjfddJN1HAAAAHzJuXPn9OGHH2ratGkaMGDANd3XdTeMAgAAwDs4gAkAAABmGEYBAABghmEUAAAAZhhGAQAAYIZhFAAAAGYYRoOopKREmZmZSkpKUnh4uBITE7VkyRKdOnXKOpqrtLS0aNWqVYqLi1NYWJgmTpyo4uJi61iud+jQIS1btkyjR49W//79NWzYMGVkZKiqqso6mietXbtWISEhGjNmjHUUzygvL5ff71dMTIzCw8M1ZswYrV+/3jqW6x0+fFhpaWmKi4tTeHi4Ro0apTVr1gTk5Ov4f2fPntUTTzyh6dOnKzo6WiEhIdq2bdslb/vuu+9q+vTpioiIUExMjBYtWqT6+vpub4tTOwXR+PHj1dDQoPT0dI0YMUIffPCBcnJyFBYWpoqKCsXGxlpHdIX58+dr586dWrFihUaMGKEXXnhBhw4d0p49e3T33Xdbx3OtuXPn6o033lB6erqSk5NVW1urnJwcnTlzRvv379fo0aOtI3rGiRMndOuttyokJEQJCQk6cuSIdSTXKyoq0qxZszRu3DhlZGSof//+ev/99+U4jn71q19Zx3Ott99+WxMmTFBcXJx++tOfKjo6Wvv27dPWrVvl9/uVn59vHdE1PvzwQ33jG9/QsGHDlJCQoNLSUm3dulWLFi3qdLsTJ05o7NixioqKUlZWlhobG/Wb3/xGQ4cO1cGDB9WnT5/Lb8xB0Lz++utd1srKyhyfz+esXr3aIJH7HDhwwPH5fM5vf/vbjrXm5mZn+PDhzl133WWYzP327dvntLa2dlqrqqpyQkNDnYULFxql8qaMjAxnypQpTkpKinPbbbdZx3G9zz77zImNjXXmzJljHcVzHnvsMcfn8znvvPNOp/UHHnjA8fl8TkNDg1Ey92lpaXHq6uocx3Gcw4cPOz6fz9m2bVuX2z388MNOeHi4U11d3bFWXFzs+Hw+57nnnuvWtnibPogmTZrUZW3y5MmKjo5WZWWlQSL32bFjh3r37q2f/OQnHWv9+vXT4sWL9cYbb6impsYwnbvdeeed6t27d6e14cOH65vf/Cb7dw8qKyvTzp07lZ2dLcdx5PP5rCO53osvvqiPP/5Ya9eulfSftzMvXrxonMob2q/M+PWvf73T+qBBg9SrVy/17dvXIpYr9e3bt6Nn5yveRN+5c6e+//3v65ZbbulYS01NVVJSkv761792a1sMoz3szJkzamxsvOZLZ+E/3nzzTSUlJal///6d1idMmCBJqqiosIjlWY7jqK6ujv27h7S1temRRx7RkiVL+FhEDyouLlZkZKSqq6t16623KiIiQjfffLN+9rOfqaWlxTqeq2VmZio2NlaLFy/WW2+9perqav3lL3/Rhg0blJWVxWXEe1hNTY0++eQTjR8/vsv3JkyYoDfffLNb99P78jdBIGVnZ6u1tVUZGRnWUVyhtrZWgwcP7rLevnby5MmejuRp27dv18mTJ/WLX/zCOoonbNiwQR999JF2795tHcVTqqqqdOHCBc2ePVsPPfSQ1q1bpz179mj9+vVqaGjQiy++aB3RteLi4vSPf/xDM2fO1NixYzvWV69eraefftowmTfV1tZK0v98Hf73v/+t1tbWy35ulGG0mxzH6fZvvKGhoZdcLysr01NPPaWMjAylpKQEMJ13nTt3Tv369euy3v7/gKMre05lZaWWLl2qu+66Sw888IB1HNc7ffq0Hn/8cT3++OOKiYmxjuMpZ86cUVNTkx5++GFlZ2dLkmbPnq3z589r48aNevrppzV8+HDjlO5UV1enGTNmSJKef/55xcTEqKCgQGvXrlVsbKyWLl1qnNBb2l9jL/c6zDAaIK+99pruueeebt22srJSSUlJXdbuu+8+JScna9OmTcGI6Ek33XTTJX9JaG5u7vg+gu/UqVO69957FRUVpR07dvC5xR6wevVqDRgwQI888oh1FM9pf16ZP39+p/X58+dr48aN2r9/P8NokKxZs0Y1NTU6duyY4uLiJP3nF4GLFy9q1apVmj9/vqKjo41Tekf7v4VrfR1mGO2mUaNGaevWrd267aBBgzr9d3V1taZOnaqoqCgVFhYqPDw8CAm9afDgwZd8K779rYP2JysEz2effaYZM2bo888/1+uvv95l/0fgVVVV6fnnn1d2drZOnDjRsd7c3Kzz58/r+PHjioyMVFRUlGFK94qLi9M777zT5fR87Qd7fPrppxaxPGHv3r0aO3Zsl+f2WbNmaevWraqoqOj2H45w7drfnm9/zf2i2tpaxcTEdOvUTgyj3RQbG9vl3Frdcfr0aU2dOlWtra3as2cP5xYNsLFjx6q0tFSNjY2KiIjoWD9w4IAk6fbbb7eK5gnNzc2aNWuW3n//fRUXF2vkyJHWkTyhpqZGFy9eVFZWlrKysrp8PyEhQcuXL9fvfvc7g3TuN378eBUXF+vEiRMaMWJEx3r7L8YDBw60iuZ6ra2tamtru+S6JF24cKGnI3nakCFDNHDgQB06dKjL9w4ePNjt12COpg+is2fPaubMmaqtrVVhYaESExOtI7nO3Llz1dbWpueee65jraWlRS+88IImTpyoIUOGGKZzt7a2NmVkZOjAgQN66aWXdMcdd1hH8owxY8YoLy9P+fn5HV95eXkaPXq0hg0bpvz8fC1evNg6pmvNmzdPkrR58+ZO65s2bVKfPn04JiCIvvWtb6m8vLzLld5yc3PVq1cvJScnGyXzrjlz5qigoKDTuzQlJSWqqqpSenp6t+6DKzAF0ezZs/Xyyy8rMzOzy5NTRESE0tLSbIK5TEZGhvLy8rRixQolJiZq27ZtOnz4sEpKSi55rlcExvLly/Xss89q1qxZl3zCWbhwoUEqb0tJSdHp06f19ttvW0dxvYceekhbtmzRvHnz9O1vf1ulpaXasWOHHnvsMc4mEURHjhzRxIkTFRkZqWXLlik6OloFBQV69dVXtWTJEm3cuNE6oqvk5OSooaFBJ0+e1IYNG3T//fd3/LUzKytLkZGRHVdg+trXvqaf//znamxs1K9//WsNHTpUhw4d4gpM1uLj452QkBDH5/N1+UpISLCO5xrNzc3OypUrncGDBzuhoaHOHXfc4RQVFVnHcr2UlJT/uX+HhIRYx/OklJQUZ8yYMdYxPKG1tdV56qmnnPj4eKdv375OUlKS8/vf/946liccOHDAmT59uhMZGen07dvXGTlypPPLX/7SaWtrs47mOvHx8Z2e19uf80NCQpzjx4933O7o0aPOtGnTnPDwcCc6Otr50Y9+5Hz88cfd3g5/GQUAAIAZPjMKAAAAMwyjAAAAMMMwCgAAADMMowAAADDDMAoAAAAzDKMAAAAwwzAKAAAAMwyjAAAAMMMwCgAAADMMowAAADDDMAoAAAAzDKMAAAAw83+1AAVrL/tGMAAAAABJRU5ErkJggg=="
      ],
      "text/plain": [
       "Figure(PyObject <matplotlib.figure.Figure object at 0x10df93f50>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "(-3,10)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = new_model()\n",
    "dummymeshreader(m)\n",
    "calculate_node_normals(m)\n",
    "for i in keys(m.el[\"connectivity\"])\n",
    "    c = m.elno[\"coords\"][i]\n",
    "    println(i, \" \", c)\n",
    "    nor1 = m.elno[\"normals\"][i][1]\n",
    "    nor2 = m.elno[\"normals\"][i][2]\n",
    "    plot([c[1][1], c[2][1]], [c[1][2], c[2][2]], \"-ko\")\n",
    "    plot([c[1][1], c[1][1]+nor1[1]], [c[1][2], c[1][2]+nor1[2]], \"-k\")\n",
    "    plot([c[2][1], c[2][1]+nor2[1]], [c[2][2], c[2][2]+nor2[2]], \"-k\")\n",
    "end\n",
    "\n",
    "node2_normal = m.no[\"normals\"][2]\n",
    "node2_coords = m.no[\"coords\"][2]\n",
    "println(node2_coords)\n",
    "println(node2_normal)\n",
    "plot([node2_coords[1], node2_coords[1]+node2_normal[1]*2],\n",
    "     [node2_coords[2], node2_coords[2]+node2_normal[2]*2], \"-r\")\n",
    "\n",
    "grid()\n",
    "xlim(-3, 10)\n",
    "ylim(-3, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here small black lines are element *local* normal directions. Red line is *global* normal direction, defined in node point. Element local directions are in field `model.elno[\"normals\"][elid][id]` while global direction defined in nodes are in `model.no[\"normals\"][node_id]`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2, poisson equation in 1d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we solve PDE $ku'' = 1$ in domain $I=[0, 1]$ with boundary conditions $u(0)=0$, $u'(1)=1$. For this we need to define field variable for $k$, let's call it to stiffness and it must be defined in element nodes because we have to integrate over it. $k$ not necessarily have to be constant that way. After standard FE assembly and solution discrete $u$ values are in nodal field, let's call it to displacement. \n",
    "\n",
    "Accurate solution is $u=\\frac{1}{2}x^2$ for $k=1$.\n",
    "\n",
    "First we must again write some solver to generate geometry. In 1d this is easy enough (even for me ;))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "generate_1d_mesh (generic function with 1 method)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function generate_1d_mesh(model, N, I)\n",
    "    \"\"\"Generate 1D mesh with N nodes for interval I.\"\"\"\n",
    "    coords = new_field(model.no, \"coords\")\n",
    "    connectivity = new_field(model.el, \"connectivity\")\n",
    "    X = linspace(I[1], I[2], N)\n",
    "    for (i, x) in enumerate(X)\n",
    "        model.no[\"coords\"][i] = [x]\n",
    "    end\n",
    "    for i=1:N-1\n",
    "        model.el[\"connectivity\"][i] = [i, i+1]\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "node 1 x = [0.0]\n",
      "node 2 x = [0.25]\n",
      "node 3 x = [0.5]\n",
      "node 4 x = [0.75]\n",
      "node 5 x = [1.0]\n",
      "element 1 connects to nodes [1,2]\n",
      "element 2 connects to nodes [2,3]\n",
      "element 3 connects to nodes [3,4]\n",
      "element 4 connects to nodes [4,5]\n"
     ]
    }
   ],
   "source": [
    "N = 5  # Number of nodes in model\n",
    "mo = new_model()\n",
    "generate_1d_mesh(mo, N, [0, 1])\n",
    "no2elno(mo, \"coords\") # again move global coords to node local coords\n",
    "for i=1:N\n",
    "    println(\"node \", i, \" x = \", mo.no[\"coords\"][i])\n",
    "end\n",
    "for i=1:N-1\n",
    "    println(\"element \", i, \" connects to nodes \", mo.el[\"connectivity\"][i])\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our elements need stiffness and force. Let's define new field for them. Note that both are defined on *element nodes* (we must interpolate them in assembly)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dict{Any,Any} with 4 entries:\n",
       "  4 => {2=>[1.0],1=>[1.0]}\n",
       "  2 => {2=>[1.0],1=>[1.0]}\n",
       "  3 => {2=>[1.0],1=>[1.0]}\n",
       "  1 => {2=>[1.0],1=>[1.0]}"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "new_field(mo.elno, \"stiffness\")\n",
    "new_field(mo.elno, \"force\")\n",
    "for i=1:N-1\n",
    "    mo.elno[\"stiffness\"][i] = Dict()\n",
    "    mo.elno[\"stiffness\"][i][1] = [1.0]\n",
    "    mo.elno[\"stiffness\"][i][2] = [1.0]\n",
    "    mo.elno[\"force\"][i] = Dict()\n",
    "    mo.elno[\"force\"][i][1] = [1.0]\n",
    "    mo.elno[\"force\"][i][2] = [1.0]\n",
    "end\n",
    "mo.elno[\"stiffness\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next the solver. For simplicity boundary conditions are handled \"silly\". We also skip defining function spaces for now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "poissonsolver1d (generic function with 1 method)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function poissonsolver1d(model)\n",
    "    # Note that we could save all intermediate variables like jacobian, stiffness etc. to\n",
    "    # model gauss points. TODO.\n",
    "    ipoints = [-sqrt(1/3) sqrt(1/3)]\n",
    "    iweights = [1 1]\n",
    "    K = zeros(N, N)\n",
    "    F = zeros(N)\n",
    "    basis = get_shape_functions(\"Lagrange\", 1, 1)\n",
    "    #println(basis(-1), \" \", basis(1))\n",
    "    dbasis = get_shape_function_derivatives(\"Lagrange\", 1, 1)\n",
    "    for i=1:N-1  # Loop through elements\n",
    "\n",
    "        # create local stiffness matrix and force vector\n",
    "        k = zeros(2, 2)\n",
    "        f = zeros(2)\n",
    "        for j=1:length(iweights)\n",
    "            wi = iweights[j]\n",
    "            ip = ipoints[j]\n",
    "            N = basis(ip)\n",
    "            dNdxi = dbasis(ip)\n",
    "            J = norm(interpolate(dbasis, model.elno[\"coords\"][i], ip))\n",
    "            dNdx = dNdxi*J^-1\n",
    "\n",
    "            # stiffness and force in integration point, interpolated from element nodes\n",
    "            # (not have to be constant then necessary)\n",
    "            stiff = interpolate(basis, model.elno[\"stiffness\"][i], ip)[1]\n",
    "            force = interpolate(basis, model.elno[\"force\"][i], ip)[1]\n",
    "            k += wi * stiff * dNdx'*dNdx * J\n",
    "            # FIXME: There's a mistake somewhere ...\n",
    "            f -= wi * force * N' * J\n",
    "            #println(N')\n",
    "        end\n",
    "        \n",
    "        # assemble to global matrices\n",
    "        gdofs = model.el[\"connectivity\"][i]\n",
    "        K[gdofs, gdofs] += k\n",
    "        F[gdofs] += f\n",
    "    end\n",
    "    # add neumann boundary condition to RHS and solve\n",
    "    F[end] += 1.0\n",
    "    u = zeros(N)\n",
    "    free_dofs = 2:N\n",
    "    u[free_dofs] = K[free_dofs, free_dofs] \\ F[free_dofs]\n",
    "    #println(K[free_dofs, free_dofs])\n",
    "    #println(F[free_dofs])\n",
    "    #println(u)\n",
    "    # create new field and save results to nodes\n",
    "    new_field(model.no, \"displacement\")\n",
    "    for i=1:N\n",
    "        model.no[\"displacement\"][i] = u[i]\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "poissonsolver1d(mo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAqkAAAIUCAYAAAA5crVsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3X2cjXX+x/H3MBqGsJValEgrlQ0jdkR3W1TU6Y5kowylQkQZtEJ3m5kUuWl/a41xk8YUmUrKbRSJNSe1qclWNHKzmnRjGyYzzu+Pa802xs0c55zre53zfT0fj3lU15yb99n3nvq4br5XXCAQCAgAAADwkEqmAwAAAACHY0gFAACA5zCkAgAAwHMYUgEAAOA5DKkAAADwHIZUAAAAeA5DKgAAADyHIRUAAACew5AKAAAAz2FIBQAAgOcENaQWFRVp2LBhqlevnhITE5WcnKxly5Yd93kzZsxQpUqVjvize/fuEw4PAACA2BQfzIN79eql+fPna/Dgwfrd736nzMxMderUSe+8847atWt33Oc/8cQTatSoUZlttWrVCi4xAAAAYl5cIBAIVOSB69evV3JyssaNG6chQ4ZIcvasNmvWTKeffrrWrFlz1OfOmDFDvXv31oYNG5SUlBSe5AAAAIhZFT7cP2/ePMXHx6tv376l2xISEtSnTx+tXbtW27dvP+5rBAIB7d27VyUlJSeWFgAAAFao8JD64YcfqkmTJqpRo0aZ7a1bt5Ykbdy48bivceWVV6pWrVqqXr26brzxRn3xxRdBxgUAAIANKnxO6s6dO1W3bt1y2w9t27Fjx1GfW716daWkpOjKK69UzZo1tWHDBj333HO65JJL5Pf7deaZZ55AdAAAAMSqCg+p+/btU0JCQrntVatWLf390XTt2lVdu3Yt/Wefz6drrrlGl112mZ566in99a9/PeLzCgoKtHjxYjVs2FDVqlWraFQAAAC4ZN++fdq6dauuueYanXbaaWF73QoPqdWqVVNRUVG57fv37y/9fTDatWunP/zhD8dcwmrx4sXq0aNHUK8LAAAA97344ou64447wvZ6FR5S69ate8RD+jt37pQk1atXL+g3P/PMM7V58+aj/r5hw4aSnA99/vnnB/36iD5dunTRvHnzTMeAS+jbLvRtF/qObV26dNGWLVv0pKRTJd2v/81t4VLhIbVly5ZauXKl9u7dq5NPPrl0+7p16yRJLVq0CPrNv/rqK9WpU+eovz+0d/b8889n6SpLVKlSha4tQt92oW+70Hdsu/zyy3Xhli36s6RO/90W7lMzK3x1f5cuXVRSUqKpU6eWbisqKlJmZqaSk5NVv359SdKuXbuUl5en4uLi0sd9++235V5v0aJF8vv9uvbaa0PJjxhz3nnnmY4AF9G3XejbLvQd254ZOVJ/r1xZiyS9FaH3qPCe1DZt2qhr164aMWKEdu/ercaNG2vmzJnKz89XZmZm6eOGDx+uWbNmaevWrWrQoIEk6ZJLLlFSUpJatWqlWrVqye/3a/r06WrQoIEeeeSR8H8qAAAARMxpzzyjg1Wravl116lRbq62bNkS9vcI6raos2bN0qOPPqrZs2fr+++/V/PmzbVw4UK1b9++9DFxcXGKi4sr87zbb79db775ppYsWaLCwkLVq1dP9957r0aPHn3Mw/0AAADwmJUrpb/9TZWmTNGz/frpDr9frVq1CvvbBDWkJiQkKD09Xenp6Ud9TGZmZpk9q5L0xBNP6IknnjixhLDK9ddfbzoCXETfdqFvu9B3jCoslO65R7r0Uum++yL6VhU+JxVww8KFC01HgIvo2y70bRf6jlGjR0vbtknTpkmVIjtGBrUnFYi0MWPGmI4AF9G3XejbLvQdg/7xD+m556S//EVq0iTibxcTQ2phYaHy8vJMx8BxNG3aVImJicd8DMuV2IW+7ULfdqHvGPPLL1KfPlKLFtJDD7nyljExpObl5UXkhF2EV25uLv/SAgAgGo0dK332mbM3Nd6d8TEmhtRDuDOVN3322Wfc3hYAgGj1z39KTz4pDR/u7El1SUwNqdyZKvplZGSoT58+pmPAJfRtF/q2C33HiOJiqXdv6Xe/k0aOdPWtubofnuL3+01HgIvo2y70bRf6jhHPPSf5/dL06VJCgqtvzZAKT5kyZYrpCHARfduFvu1C3zHg88+lUaOkIUOkP/zB9bdnSAUAAEBZJSXOYf4GDaTHHzcSIabOSQUAAEAYTJkivf++9O67UrVqRiKwJxUAAAD/89VX0ogR0oABzu1PDWFIhaf4fD7TEeAi+rYLfduFvqNUICDdfbdUp4709NNGo3C4H54yYMAA0xHgIvq2C33bhb6j1N//Lr3zjrR0qVSjhtEo7EmFp3Ts2NF0BLiIvu1C33ah7yiUny89/LCzJ/Xqq02nYUgFAACwXiAg3XOPVKuWNG6c6TSSONwPAACA6dOlJUukt95yBlUPYE9qDJkxY4ZuvfVWnXPOOUpMTFStWrXUvn17zZkz54iP37Nnj/785z+rWbNmql69umrXrq0WLVpoxIgRKiwsPOHHhiInJydsrwXvo2+70Ldd6DuKfPONs2B/Sop07bWm05RiSI0h/fr107Zt23TFFVdo8ODBuv322/X111+rZ8+eGjVqVJnHbtmyRUlJSXr66aeVmJiofv36qU+fPjrzzDM1YcIEFRQUnNBjQ5WVlRW214L30bdd6Nsu9B0lAgGpb1/nIqnnnjOdpgy7D/cXFkp5ee6/b9OmUmJi2F9206ZNatSoUZltBw4c0HXXXaexY8fqvvvuU7169SRJd9xxh/Lz8/X0009r2LBhZZ6zZ88eVa9evfSfg3lsqLKzs8P2WvA++rYLfduFvqPErFnOIf6FC6XatU2nKcPuITUvT2rVyv33zc2VkpLC/rKHD6iSVKVKFfXr108rVqzQihUr1KNHD+Xm5uqDDz5Qy5Ytyw2dknTKKaf8KmrFHwsAAKLIjh3Sgw9Kd94pde5sOk05dg+pTZs6A6OJ942A/Px8paWlafny5dq2bZv27dtX5vfbt2+XJH3wwQeSpGuuuea4rxnMYwEAQJQIBKR775WqVpXGjzed5ojsHlITEyOyR9OEr776Sm3atNEPP/ygyy67TNdee61q1aqlypUra8uWLZo5c6aKiookST/88IMkqX79+sd93WAeCwAAosScOc4h/tdekzx6VJQLp2LEc889pz179mj69OlasWKFJkyYoMcee0yjRo0qt6By7f+ec/LNN98c93WDeWw4pKSkuPI+8Ab6tgt924W+PWznTmngQOlPf5I8fPtahtQY8cUXXyguLk633nprud+tWrWqzD+3bdtWkrR48WIFAoFjvm4wjw0H7lBiF/q2C33bhb496tBh/oQEadIk02mOiSE1RjRq1EiBQEDvvPNOme2LFy/WtGnTymxLSkrSJZdcoo0bNyotLa3ca3333XelpwYE89hw6N69e9heC95H33ahb7vQt0fNni298Yb0t7959jD/IXafkxpD+vXrp8zMTHXt2lVdunRR3bp19cknn2jx4sW67bbbyi0F8uKLL+qKK67QI488ovnz5+vyyy9XIBDQv/71Ly1dulSff/65GjRoEPRjAQCAR23fLg0aJPXs6enD/IcwpMaI3//+93rnnXc0cuRIvfnmmyouLlaLFi20YMEC1apVq9yQ2rBhQ/n9fqWnpysnJ0dTpkxRtWrV1LBhQz388MOqU6fOCT0WAAB40KFF+6tVk55/3nSaCmFIjSFt27bV8uXLj/i7gwcPltt2yimnaOzYsRo7duxxXzuYx4Zi9erVat++fUTfA95B33ahb7vQt8fMnCktWuQc6v/Nb0ynqRDOSYWnpKenm44AF9G3XejbLvTtId984xzmv+su6frrTaepMIZUeMrcuXNNR4CL6Nsu9G0X+vaIQEC65x6pRg1pwgTTaYLC4X54SmJioukIcBF924W+7ULfHpGZKb39tvTmm9J/1z6PFuxJBQAAiEX5+dLgwVKvXlKnTqbTBI0hFQAAINYEAlKfPlLNmlF3mP8QhlR4ytChQ01HgIvo2y70bRf6Nuz//k9atkzKyJBq1TKd5oQwpMJTuCmAXejbLvRtF/o26MsvpYcflu67T4ri29MypMJTHnjgAdMR4CL6tgt924W+DTl4UEpJkc44Q3rmGdNpQsLV/QAAALHi+eel996TVq50lp2KYjE1pH722WemI+AI6AUAABfk5UmPPOIs3H/55abThCymhtQePXqYjoAQ5eXlqWnTpqZjwCX0bRf6tgt9u6y42LmjVIMG0l/+YjpNWMTEkNq0aVPl5uaajoHjqMi/rFJTU/X666+7kAZeQN92oW+70LfLnnlG2rBBWrNGipEbKcTEkJqYmKikpCTTMRAGkydPNh0BLqJvu9C3XejbRR9/LI0eLQ0dKiUnm04TNlzdD09hyRK70Ldd6Nsu9O2SX36R7rxTOu886bHHTKcJq5jYkwoAAGClxx6TNm2S1q+XEhJMpwkr9qQCAABEow8+kMaOdQ71t2xpOk3YMaTCU9LS0kxHgIvo2y70bRf6jrDCQudq/osvloYPN50mIjjcD08pLCw0HQEuom+70Ldd6DvChg+X8vOl116T4mNznIsLBAIB0yGOxu/3q1WrVsrNzeXqfQAAAElavly6+mppwgRn4X7DIjWvcbgfAAAgWvz4o5SSIl15pfTAA6bTRBRDKgAAQLR48EHphx+kzEypUmyPcbH96RB1CgoKTEeAi+jbLvRtF/qOgNdek2bMcA7zn3226TQRx5AKT+ndu7fpCHARfduFvu1C32G2e7d0zz3SDTc4h/stwJAKTxkzZozpCHARfduFvu1C32EUCEh9+zp///e/S3FxZvO4JDbXLEDUYhUHu9C3XejbLvQdRjNmOIf6FyyQzjjDdBrXsCcVAADAq7ZudZaZ6tVLuukm02lcxZAKAADgRSUlzl2lfvMb52IpyzCkwlMyMjJMR4CL6Nsu9G0X+g6D8eOl996TZs6UatUyncZ1DKnwFL/fbzoCXETfdqFvu9B3iD75RPrzn6XBg6UrrjCdxghuiwoAAOAlv/witWkjFRdLGzZIVauaTnRMkZrXuLofAADAS8aMkT79VFq3zvMDaiQxpAIAAHjF6tVSWpr0xBNSy5am0xjFOakAAABe8NNPUs+eUtu20rBhptMYx5AKT/H5fKYjwEX0bRf6tgt9n4BBg6SCAmn2bKlyZdNpjONwPzxlwIABpiPARfRtF/q2C30H6dVXnTtLTZ8uNWpkOo0nsCcVntKxY0fTEeAi+rYLfduFvoOwc6fUt690883OnaUgiSEVAADAnEBA6tNHqlJFmjpVioszncgzONwPAABgyl//Kr31lrRokXTaaabTeAp7UuEpOTk5piPARfRtF/q2C31XwOefSw8/LPXrJ113nek0nsOQCk/JysoyHQEuom+70Ldd6Ps4DhyQevSQzjpLeuYZ02k8icP98JTs7GzTEeAi+rYLfduFvo9jzBhp40bp/felxETTaTyJIRUAAMBN770nPf209OSTUuvWptN4Fof7AQAA3PLjj85dpdq3565Sx8GeVAAAALf07y99/720ciV3lToO9qTCU1JSUkxHgIvo2y70bRf6PoKsLGnOHGnKFKlhQ9NpPI8hFZ7CHUrsQt92oW+70Pdhvv5auv9+6fbbpTvuMJ0mKjCkwlO6d+9uOgJcRN92oW+70PevlJRId94p1azpLN7PXaUqhHNSAQAAIumZZ5wr+t95R6pd23SaqMGeVAAAgEjJzZUefVRKTZUuv9x0mqjCkApPWb16tekIcBF924W+7ULfkn7+WfrTn6SLLpIef9x0mqjDkApPSU9PNx0BLqJvu9C3Xehb0uDB0jffSC+9JJ10kuk0UYdzUuEpc+fONR0BLqJvu9C3Xazve8EC6e9/l6ZOlc47z3SaqMSeVHhKIvcvtgp924W+7WJ139u3S3ffLd18s/NXnBCGVAAAgHA5eNBZbqpaNWdPKstNnbCghtSioiINGzZM9erVU2JiopKTk7Vs2bKg3/See+5RpUqVdMMNNwT9XAAAAM969llnqalZs6RTTzWdJqoFNaT26tVL48ePV8+ePTVx4kRVrlxZnTp10po1ayr8Ghs2bNDMmTNVtWpVxfGnCxxm6NChpiPARfRtF/q2i5V9+/3Sn/8sDR0q/fGPptNEvQoPqevXr1d2drbGjh2rtLQ03X333VqxYoXOPvtspaamVug1AoGABg4cqLvuuktnnHHGCYdG7GrQoIHpCHARfduFvu1iXd8//yx17+4sN/XEE6bTxIQKD6nz5s1TfHy8+vbtW7otISFBffr00dq1a7V9+/bjvsbs2bP16aef6sknn1QgEDixxIhpDzzwgOkIcBF924W+7WJd3w8+6Cw3NWcOy02FSYWH1A8//FBNmjRRjRo1ymxv3bq1JGnjxo3HfP7evXs1bNgwPfLII+xFBQAAsWPePGnaNGniRJabCqMKD6k7d+5U3bp1y20/tG3Hjh3HfP7jjz+u6tWra/DgwUFGBAAA8Kj8fOmee6SuXaXevU2niSkVHlL37dunhISEcturVq1a+vuj2bx5syZOnKhnnnlGVapUOYGYsEVeXp7pCHARfduFvu1iRd/FxdIdd0g1azqL9nNBeFhVeEitVq2aioqKym3fv39/6e+PZtCgQWrXrp1uvvnmE4goderUST6fr8xP27ZtlZOTU+ZxS5Yskc/nK/f8/v37KyMjo8w2v98vn8+ngoKCMttHjx6ttLS0Mtvy8/Pl8/nKfeEmTZpU7urFwsJC+Xy+cvcszsrKUkpKSrls3bp143P86nMMGjQoJj5HrPQR6c9x6KLLaP8ch/A5jv05+vbtGxOfI1b6iPTn6NOnT0x8jmP28dRT0vvvSy+9pG733hu9n+NXjtdHVlZW6SzWqFEjtWjRImJHyeMCFbyCqUOHDtqxY4c2bdpUZvvy5cvVoUMHvfHGG+rcuXO5561YsUJXX321Xn31VbVo0aJ0e/v27dW0aVNlZGTolFNO0cknn1zuuX6/X61atVJubq6SkpKC/WyIQvn5+fZdEWox+rYLfdsl5vtevVq6/HJp1Chp9GjTaYyK1LwWX9EHtmzZUitXrtTevXvLDJTr1q2TpDID6K/l5+dLkm655ZZyv9uxY4caNWqkCRMmaODAgUEFR2yK6X+hoRz6tgt92yWm+/7+e+cw/yWXOOuiIiIqPKR26dJF48aN09SpU/XQQw9Jcu5AlZmZqeTkZNWvX1+StGvXLv3www8699xzFR8fr6uuuqrcbuNAIKC+ffuqYcOG+vOf/6xmzZqF8SMBAABESCAg3Xuv9NNPznJT8RUepRCkCv8v26ZNG3Xt2lUjRozQ7t271bhxY82cOVP5+fnKzMwsfdzw4cM1a9Ysbd26VQ0aNNBZZ52ls846q9zrDRo0SGecccYRz9UAAADwpOnTpVdecX5ieW+xBwR1W9RZs2bpwQcf1OzZszVo0CCVlJRo4cKFat++felj4uLiKnS7U26JiiM5/ORwxDb6tgt92yUm+/7sM2ngQOnuu6UuXUyniXlB7aNOSEhQenq60tPTj/qYzMzMMntWj2bLli3BvDUsUVhYaDoCXETfdqFvu8Rc3/v3S7ff7uw9nTDBdBorVPjqfhO4uh8AAHjCAw9If/+7tH69dNFFptN4ivGr+wEAAKyUkyNNnuz8MKC6JqhzUgEAAKyybZtzu9ObbpL69TOdxioMqfCUw++0gdhG33ahb7vERN+HbntavbqUkcFtT13GkApP6d27t+kIcBF924W+7RITfT/5pLRmjfTSS9Ipp5hOYx2GVHjKmDFjTEeAi+jbLvRtl6jve9Uq6YknnNueXnqp6TRWYkiFp7CKg13o2y70bZeo7rugwDnM3769NHKk6TTWYkgFAAA4JBCQUlKcdVHnzJEqVzadyFosQQUAAHDIhAnSwoXSG29IZ55pOo3V2JMKT8nIyDAdAS6ib7vQt12isu8NG6Rhw6TBg6XrrzedxnoMqfAUv99vOgJcRN92oW+7RF3fP/4odesmNW8ujR1rOg3E4X54zJQpU0xHgIvo2y70bZeo6jsQkO6917lgaulS6aSTTCeCGFIBAIDtpk2TsrOdn3POMZ0G/8XhfgAAYK9//lMaOFDq21e67TbTafArDKkAAMBOP//snId67rnOVf3wFIZUeIrP5zMdAS6ib7vQt12iou8HHpC+/lp6+WWpWjXTaXAYzkmFpwwYMMB0BLiIvu1C33bxfN+zZkmZmdKMGdL555tOgyNgTyo8pWPHjqYjwEX0bRf6toun+/7sM+n++6VevaS77jKdBkfBkAoAAOxRWCh17So1bChNnmw6DY6Bw/0AAMAeAwdKX30l/eMfUvXqptPgGNiTCk/JyckxHQEuom+70LddPNn3iy9KGRnSlCnShReaToPjYEiFp2RlZZmOABfRt13o2y6e6zsvT7rvPunOO51zUeF5cYFAIGA6xNH4/X61atVKubm5SkpKMh0HAABEo337pD/8QTpwwDnMX6OG6UQxJVLzGuekAgCA2DZwoPTFF9L69QyoUYQhFQAAxK7Zs6Vp05xzUZs1M50GQeCcVAAAEJs+/dQ5D/Wuu6SUFNNpECSGVHhKCv8SsQp924W+7WK8759/lrp0kRo1cq7mj4szmwdB43A/PMXTdyhB2NG3XejbLkb7DgScO0rl57MeahRjSIWndO/e3XQEuIi+7ULfdjHa9/Tpzrmos2dL559vLgdCwuF+AAAQOz76SBowQOrbV+rRw3QahIAhFQAAxIaffpK6dpWaNpWef950GoSIIRWesnr1atMR4CL6tgt928X1vgMBZ+/prl3SK69IVau6+/4IO4ZUeEp6errpCHARfduFvu3iet9TpkjZ2c75qOee6+57IyIYUuEpc+fONR0BLqJvu9C3XVzte906acgQ6cEHnWWnEBMYUuEpiYmJpiPARfRtF/q2i2t9f/edcx7qxRdLaWnuvCdcwZAKAACi08GDzhX8+/Y5h/pPOsl0IoQR66QCAIDo9Je/SIsXS2+/LZ11luk0CDP2pMJThg4dajoCXETfdqFvu0S872XLpFGjpNGjJe5mFpMYUuEpDRo0MB0BLqJvu9C3XSLa9/bt0p/+JHXoII0cGbn3gVFxgUAgYDrE0fj9frVq1Uq5ublKSkoyHQcAAJh24IB05ZXS119Lfr9Up47pRNaL1LzGOakAACB6pKZK69dLq1YxoMY4hlQAABAdXn5ZmjBBmjRJatvWdBpEGOekwlPy8vJMR4CL6Nsu9G2XsPf92WdS795S9+5S//7hfW14EkMqPCU1NdV0BLiIvu1C33YJa99790q33CKdfbY0daoUFxe+14ZncbgfnjJ58mTTEeAi+rYLfdslbH0HAtLdd0vffCNt2CDVqBGe14XnMaTCU1iixi70bRf6tkvY+p440TkX9ZVXpPPOC89rIipwuB8AAHjTmjXSww9LQ4ZIXbqYTgOXMaQCAADv+fe/pdtuc67iHzvWdBoYwJAKT0lLSzMdAS6ib7vQt11C6vvAAWdAPXhQys6WqlQJXzBEDc5JhacUFhaajgAX0bdd6NsuIfU9fLj0/vvSO+9IdeuGLxSiCrdFBQAA3pGdLd1+u7No/6BBptOgAiI1r3G4HwAAeMOmTVKfPtKf/iQNHGg6DQxjSAUAAOb9+KN0883SOeewYD8kMaTCYwoKCkxHgIvo2y70bZeg+j54ULrrLmn3bunVV6Xq1SMXDFGDIRWe0rt3b9MR4CL6tgt92yWovseOlV57TZo9Wzr33MiFQlRhSIWnjBkzxnQEuIi+7ULfdqlw30uWSCNHOj833BDRTIguDKnwFFZxsAt924W+7VKhvrdskbp3lzp2lPhDDA7DkAoAANxXWCjdcotUu7b00ktS5cqmE8FjWMwfAAC4KxCQ+vaVPv9cWrtWOuUU04ngQexJhadkZGSYjgAX0bdd6Nsux+x74kRpzhwpI0Nq3ty9UIgqDKnwFL/fbzoCXETfdqFvuxy171WrpIcekoYMcc5HBY6C26ICAAB3bNsmtWolNWvmXNUfz1mHsYDbogIAgOi1f790661S1apSdjYDKo6L/4cAAIDICgSkAQOkjz+WVq+W6tQxnQhRgCEVAABE1v/9n3ORVGamdPHFptMgSnC4H57i8/lMR4CL6Nsu9G2X0r7fe08aOFB64AGpVy+jmRBdGFLhKQMGDDAdAS6ib7vQt10GDBjgXCjVpYvUrp307LOmIyHKcLgfntKxY0fTEeAi+rYLfdul46WXSpde6lwo9corUpUqpiMhyjCkAgCA8AoEpPvukzZtktas4UIpnBCGVAAAEF4TJ0qzZkkvviixzjlOEOekwlNycnJMR4CL6Nsu9G2JFSukhx5Sjs8n3XGH6TSIYgyp8JSsrCzTEeAi+rYLfVtgyxbpttukK69UVkKC6TSIchzuh6dkZ2ebjgAX0bdd6DvG/ec/0k03SbVqSXPnKvvUU00nQpRjSAUAAKE5eNBZA/Wrr6S1ayUGVIQBQyoAAAjNU09J8+dLCxZIzZqZToMYwTmpAADgxL32mjRqlPTYY87hfiBMGFLhKSkpKaYjwEX0bRf6jkGffCL16CHdeqs0cmSZX9E3QsWQCk/hjjR2oW+70HeM2bNHuvFG6ZxzpBkzpEplRwr6Rqg4JxWe0r17d9MR4CL6tgt9x5DiYqlbN+nHH6Vly6QaNco9hL4RKoZUAAAQnIcflt55R1q6VGrUyHQaxCiGVAAAUHHTpknPPy9NmSJdeaXpNIhhnJMKT1m9erXpCHARfduFvmPAe+9J/fpJ993n/PUY6BuhCmpILSoq0rBhw1SvXj0lJiYqOTlZy5YtO+7z3n33Xfl8PjVo0EDVqlXTGWecoauuukpvvfXWCQdHbEpPTzcdAS6ib7vQd5TbulW65RapXTtp4sTjPpy+EaqghtRevXpp/Pjx6tmzpyZOnKjKlSurU6dOWrNmzTGf969//Uvx8fG6//779cILLyg1NVV79uxR586dNWfOnJA+AGLL3LlzTUeAi+jbLvQdxfbulXw+6eSTpVdekapUOe5T6BuCQMARAAAgAElEQVShigsEAoGKPHD9+vVKTk7WuHHjNGTIEEnOntVmzZrp9NNPP+6gerh9+/bpnHPOUZMmTbRq1aojPsbv96tVq1bKzc1VUlJSUK8PAADC4OBBZw/q8uXSBx9IF15oOhE8JlLzWoX3pM6bN0/x8fHq27dv6baEhAT16dNHa9eu1fbt24N642rVqum0005TlQr8aQwAABgyapT0+utSVhYDKlxV4av7P/zwQzVp0kQ1DlsLrXXr1pKkjRs3qn79+sd8jZ9++km//PKLCgoKNGvWLG3evJlzVgAA8KqsLOmpp6SxY6XrrzedBpap8J7UnTt3qm7duuW2H9q2Y8eO477GbbfdptNPP10XXHCBJk2apJdfflnXXXddEHER64YOHWo6AlxE33ah7yizbp2UkiLdcYeUmhr00+kboarwkLpv3z4lJCSU2161atXS3x9PWlqali5dqoyMDF1wwQW6/fbbtXTp0iDiItY1aNDAdAS4iL7tQt9RZNs26aabpJYtnXVR4+KCfgn6RqgqPKRWq1ZNRUVF5bbv37+/9PfH07x5c1111VVKSUnR6tWrde6556p///7HfV6nTp3k8/nK/LRt21Y5OTllHrdkyRL5fL5yz+/fv78yMjLKbPP7/fL5fCooKCizffTo0UpLSyuzLT8/Xz6fT3l5eWW2T5o0qdyfFAsLC+Xz+cqtD5eVlaWUlJRy2bp168bn+NXnOPw2etH6OWKlj0h/jgceeCAmPschfI5jf4527drFxOeIlT6O+jlGjFBacrJzBX9OjlS16gl9jpYtW5r9HLHSh8c+R1ZWVuks1qhRI7Vo0UKDBw8u9zrhUOGr+zt06KAdO3Zo06ZNZbYvX75cHTp00BtvvKHOnTsH9eaPPPKIxo4dqz179qh27drlfs/V/QAAuOjgQalrV2nxYmnNGql5c9OJEAWMX93fsmVLbd68WXv37i2zfd26dZKkFi1aBP3mh04RqFSJG18BAGDcqFHSggXSSy8xoMK4Ck+HXbp0UUlJiaZOnVq6raioSJmZmUpOTi69sn/Xrl3Ky8tTcXFx6eN2795d7vV++OEHzZ8/X7///e9Vs2bNUD4DYsjhhy4Q2+jbLvTtcXPm/O9K/iMcog4WfSNUFV6Cqk2bNuratatGjBih3bt3q3Hjxpo5c6by8/OVmZlZ+rjhw4dr1qxZ2rp1a+lJ09ddd53OOusstWnTRqeffnrpc7799lvNmDEj7B8K0Ss1NVWvv/666RhwCX3bhb49bO1aqU8f6a67pDBdlU/fCFWFh1RJmjVrlh599FHNnj1b33//vZo3b66FCxeqffv2pY+Ji4tT3GFXAfbp00dz587VhAkT9MMPP+jUU09V+/btNWLECM41RRmTJ082HQEuom+70LdHff21cyX/xRdLf/vbCV3JfyT0jVBV+MIpE7hwCgCACPrpJ6ldO+nnn511UevUMZ0IUShS81pQe1IBAECMKC6Wbr9dys93DvczoMJjGFIBALDRkCHSkiXSW29JF1xgOg1QDms/wVMOX7AYsY2+7ULfHjJlijRpkjR5stShQ0Tegr4RKoZUeEphYaHpCHARfduFvj3i7belgQOlQYOk++6L2NvQN0LFhVMAANjik0+kSy6RLrtMeu01qXJl04kQA4zfcQoAAESxf/9buv56qVEjKSuLARWex5AKAECs27dPuvFGqahIeuMN6eSTTScCjoshFZ5SUFBgOgJcRN92oW93fPvtt0pJSdGFF16o8847T80uuED/OP98BT7+2BlQ/3s3yEijb4SKIRWe0rt3b9MR4CL6tgt9R97u3bvVtm1bzZgxQ59++qk2b96sHp99plZff637a9XSt2ef7VoW+kaoGFLhKWPGjDEdAS6ib7vQd+QNGzZMX375Zek/3y1puKSHJf1t1y6lpqa6loW+ESqGVHgKqzjYhb7tQt+Rt379+tK/v1rSXyW9IGn8EX4fafSNUDGkAgAQI4qLiyVJF0iaJ2mppIFH+D0QDRhSAQCIEfHx8Tpd0puSvpbUTVLJYb8HogVDKjwlIyPDdAS4iL7tQt+R175lSy2UdJKkzpL2Hvb7Nm3auJaFvhEqhlR4it/vNx0BLqJvu9B3hJWUaGJBgS6Ii1NnSd8c9uvGjRsrPT3dtTj0jVCx3x+eMmXKFNMR4CL6tgt9R1AgID34oBKWLtX+OXPUYskS/bJ+vYqLixUfH682bdooPT1dderUcS0SfSNUDKkAAES7CROkyZOlv/5Vtbp3V2b37qYTASHjcD8AANFs/nzpoYek1FTpvvtMpwHChiEVAIBotXat1KOH1LWr9PTTptMAYcWQCk/x+XymI8BF9G0X+g6zL7+UfD6pVStp5kypkrf+k07fCJW3/h8N6w0YMMB0BLiIvu1C32FUUCBdd530m99Ir70mVa1qOlE59I1QceEUPKVjx46mI8BF9G0X+g6TwkLphhukH35wDvefeqrpREdE3wgVQyoAANGipES64w7p44+llSulxo1NJwIihiEVAIBo8N+1UPX6684h/tatTScCIopzUuEpOTk5piPARfRtF/oO0bPPOmuhvvCCdP31ptMcF30jVAyp8JSsrCzTEeAi+rYLfYdg7lxp6FDpkUeke+81naZC6BuhYkiFp2RnZ5uOABfRt13o+wStWiXddZfUs6f05JOm01QYfSNUDKkAAHjVJ59IN94oXXqpNG2aFBdnOhHgGoZUAAC8aNs26dprpYYNpVdflU46yXQiwFUMqQAAeM333zsDany89NZbUs2aphMBrmNIhaekpKSYjgAX0bdd6LuC9u93DvHv2iUtXizVrWs60Qmhb4SKdVLhKdyhxC70bRf6roCSEqlHD2nDBmn5cum880wnOmH0jVAxpMJTunfvbjoCXETfdqHv4zi0WP+CBc5P27amE4WEvhEqhlQAALwgLc1ZrP9vf5N8PtNpAOM4JxUAANNmzpRGjJBGjZL69jWdBvAEhlR4yurVq01HgIvo2y70fRRvvin16SPdc480ZozpNGFD3wgVQyo8JT093XQEuIi+7ULfR/D++1LXrtINN0gvvBBTi/XTN0LFkApPmTt3rukIcBF924W+D/Ppp9L110sXXyy99JKzJmoMoW+EiiEVnpKYmGg6AlxE33ah71/Ztk265hrpzDOl11+XqlUznSjs6BuhYkgFAMBNe/Y4A2rlytLbb0u1a5tOBHhSbB1bAADAywoLnUP8334rrVkj1atnOhHgWexJhacMHTrUdAS4iL7tYn3fBw44F0l9/LFzRX+TJqYTRZT1fSNk7EmFpzRo0MB0BLiIvu1idd8HD0opKdLSpdLChVKbNqYTRZzVfSMs4gKBQMB0iKPx+/1q1aqVcnNzlZSUZDoOAADBO3S700mTpLlzpdtuM50ICKtIzWvsSQUAIJKeekqaONFZB5UBFagwzkkFACBS/u//pEcflR5/XLr/ftNpgKjCkApPycvLMx0BLqJvu1jX98svS/36SQ88II0caTqN66zrG2HHkApPSU1NNR0BLqJvu1jV99KlUo8eUvfu0oQJMXW704qyqm9EBEMqPGXy5MmmI8BF9G0Xa/r+4APp5pulq6+WZsyQKtn5n1pr+kbE2PnNgWexZIld6NsuVvT9z39K110ntWwpzZsnValiOpExVvSNiGJIBQAgHL74QurYUWrYUHrjDYl71wMhYUgFACBU27dLHTpINWtKixdLtWubTgREPYZUeEpaWprpCHARfdslZvv+7jtnD2pJiXPB1Omnm07kCTHbN1zDYv7wlMLCQtMR4CL6tktM9r13r3MO6rffSu+9J3EeZqmY7Buu4raoAACciP37pU6dpNxcaeVK52IpwELcFhUAAK84cEDq2tVZbmrxYgZUIAIYUgEACEZJidSzp7RkifT669Kll5pOBMQkLpyCpxQUFJiOABfRt11iou+DB6V773XWQJ07V7rmGtOJPCsm+oZRDKnwlN69e5uOABfRt12ivu9AQBoyRMrIkDIznbtK4aiivm8Yx5AKTxkzZozpCHARfdsl6vsePVp6/nnphRecw/04pqjvG8YxpMJTWMXBLvRtl6ju+5lnpCeekNLSpPvvN50mKkR13/AEhlQAAI7lhRek1FRp5EjnrwBcwZAKAMDRZGZK/ftLDz4oPf646TSAVRhS4SkZGRmmI8BF9G2XqOs7K0vq08e5mv+556S4ONOJokrU9Q3PYUiFp/j9ftMR4CL6tktU9b1ggXNxVM+ezuF+BtSgRVXf8CRuiwoAwK+9/bbk8zlLTM2ZI8Vz3xvgWCI1r7EnFQCAQ955xxlOr7tOevFFBlTAIIZUAAAk6f33pRtukC67TMrOlqpUMZ0IsBpDKgAA69c7e08vvtg5H7VqVdOJAOsxpMJTfD6f6QhwEX3bxbN9+/3SNddIv/+9tHChlJhoOlFM8GzfiBoMqfCUAQMGmI4AF9G3XTzZ98cfSx06SOedJy1aJNWoYTpRzPBk34gqDKnwlI4dO5qOABfRt1081/emTdJVV0kNGzpX9NesaTpRTPFc34g6DKkAAPt8/rkzoNarJy1ZItWubToRgMMwpAIA7PLFF9If/yideqq0bJnzVwCew5AKT8nJyTEdAS6ib7t4ou+vvnIG1JNPlpYvl+rUMZ0oZnmib0Q1hlR4SlZWlukIcBF928V431u2SFde6SwvtWKF9Nvfms0T44z3jajHkApPyc7ONh0BLqJvuxjte+tWZ0A96STnrlL16pnLYgm+3wgVQyoAILZ9/bUzoMbHOwNq/fqmEwGoAIZUAEDsys93BtRKlZwB9cwzTScCUEHxpgMAABAR27Y5A2ogIK1cKZ11lulEAILAnlR4SkpKiukIcBF928XVvr/5xhlQS0qcPagNGrj33pDE9xuhY0iFp3CHErvQt11c63vbNunyy6XiYmdAbdjQnfdFGXy/Eaqgh9SioiINGzZM9erVU2JiopKTk7Vs2bLjPm/58uXq3bu3mjRpourVq6tx48a65557tGvXrhMKjtjUvXt30xHgIvq2iyt95+dLV1whHTzoHOJv1Cjy74kj4vuNUAU9pPbq1Uvjx49Xz549NXHiRFWuXFmdOnXSmjVrjvm8YcOG6d1339Wtt96qSZMm6fbbb9fLL7+sli1b6t///vcJfwAAACQ5V/H/ekBlDyoQ1YK6cGr9+vXKzs7WuHHjNGTIEElSz5491axZM6Wmph5zUJ0wYYLat29fZtu1116ryy+/XJMnT9YTTzxxAvEBAND/1kGNi5NWreIcVCAGBLUndd68eYqPj1ffvn1LtyUkJKhPnz5au3attm/fftTnHj6gStKll16qU045RXl5ecHEQAxbvXq16QhwEX3bJWJ9b93q7EGtVIkB1UP4fiNUQQ2pH374oZo0aaIaNWqU2d66dWtJ0saNG4N68//85z/au3evTjvttKCeh9iVnp5uOgJcRN92iUjfX33lXCQVH+8MqCwz5Rl8vxGqoIbUnTt3qm7duuW2H9q2Y8eOoN58woQJOnDggLp16xbU8xC75s6dazoCXETfdgl73//6lzOgnnSScw4qC/V7Ct9vhCqoIXXfvn1KSEgot71q1aqlv6+od999V4899pi6deumK664IpgYiGGJiYmmI8BF9G2XsPadl+cMqDVqOHtQGVA9h+83QhXUkFqtWjUVFRWV275///7S31dEXl6ebr75Zl100UWaNm3acR/fqVMn+Xy+Mj9t27ZVTk5OmcctWbJEPp+v3PP79++vjIyMMtv8fr98Pp8KCgrKbB89erTS0tLKbMvPz5fP5yt37uykSZM0dOjQMtsKCwvl8/nKnYuTlZV1xIWNu3Xrxufgc/A5+Bx8jmA+R0qKM6Cecoq0cqUKa9eOzs8RK33wOaz6HFlZWaWzWKNGjdSiRQsNHjy43OuEQ1wgEAhU9MEdOnTQjh07tGnTpjLbly9frg4dOuiNN95Q586dj/ka27ZtU7t27XTSSSdpzZo1OuOMM476WL/fr1atWik3N1dJSUkVjQkAiFUffSRdfbVUv760dKlUp47pRID1IjWvBbUntWXLltq8ebP27t1bZvu6deskSS1atDjm87/77jt17NhRBw4c0OLFi485oMJOh/8JELGNvu0Sct9+v/THP0pnny2tWMGA6nF8vxGqoIbULl26qKSkRFOnTi3dVlRUpMzMTCUnJ6t+/fqSpF27dikvL0/FxcWlj/v555/VqVMn7dy5U4sWLVLjxo3D9BEQSxqwdIxV6NsuIfW9fr101VXS734nLVvmHOqHp/H9RqiCOtwvOecqLFiwQIMHD1bjxo01c+ZMbdiwQcuXLy9dC7VXr16aNWuWtm7dWvp/0ptuukmvv/66evfuXe5CqZNPPlk33nhjufficD8AQO+9J3XuLF10kbRokVSzpulEAH4lUvNaUHeckqRZs2bp0Ucf1ezZs/X999+refPmWrhwYZnF+uPi4hQXF1fmeR999JHi4uI0ffp0TZ8+vczvGjZseMQhFQBguWXLpBtvlJKTpddec67mB2CFoPekuok9qQBgsTfflG691TkPdf58qYIryABwlycunAIijVvk2oW+7RJU3/PnSzffLHXqJC1YwIAahfh+I1QMqfCU1NRU0xHgIvq2S4X7njNH6tZN6tJFys6WjnATGXgf32+EiiEVnjJ58mTTEeAi+rZLhfqeNk3q2VO6805p9mypSpXIB0NE8P1GqBhS4SksWWIX+rbLcfueMEG65x7p/vudYbVyZXeCISL4fiNUDKkAALMCAemJJ6TBg6Vhw6TJk6VK/OcJsF3QS1ABABA2gYCUmiqNGyf95S/SiBGmEwHwCP6oCk9JS0szHQEuom+7lOu7pMQ5tD9unDRxIgNqjOH7jVCxJxWeUlhYaDoCXETfdinT94EDUq9e0ty5Umam8/eIKXy/ESoW8wcAuGv/fun2251bnL70krPUFICo5ZnbogIAcML27nVuc7p2rZST4yzWDwBHwJAKAHDHd99J110nff65tGSJdOmlphMB8DAunIKnFBQUmI4AF9G3RbZvV8Ell0hbt0orVzKgWoDvN0LFkApP6d27t+kIcBF9W+Jf/5LatVPvbduk1aulli1NJ4IL+H4jVAyp8JQxY8aYjgAX0bcFPvrI2WtarZrGzJsnNWliOhFcwvcboeKcVHgKqzjYhb5j3Jo1UufO0rnnSm+9paQ6dUwngov4fiNU7EkFAITfwoXS1Vc7h/ZXrJAYUAEEiSEVABBes2ZJN93kXMn/1ltSzZqmEwGIQgyp8JSMjAzTEeAi+o5Bzz4r3XWXlJIivfKKVLVq6a/o2y70jVAxpMJT/H6/6QhwEX3HkEBAGjZMevhh6ZFHpKlTpcqVyzyEvu1C3wgVt0UFAISmuFi6915p+nRp/HjpwQdNJwLgIm6LCgDwnsJCqXt3adEi6cUXpTvuMJ0IQIxgSAUAnJg9e6QbbpA2bpRef925UAoAwoQhFQAQvPx86dprpW+/ld55R2rTxnQiADGGC6fgKT6fz3QEuIi+o9Qnn0iXXCLt2+cs2F/BAZW+7ULfCBVDKjxlwIABpiPARfQdhd59V2rf3lmc//33g7rNKX3bhb4RKoZUeErHjh1NR4CL6DvKLFggdewotWolrVol1a0b1NPp2y70jVAxpAIAjm/yZOnWW507SS1axF2kAEQcQyoA4OgOHpRSU6UHHpCGDJFeeklKSDCdCoAFGFLhKTk5OaYjwEX07XFFRc66p+PGSRMmOH+tdOL/2aBvu9A3QsWQCk/JysoyHQEuom8P++EHZ4mpBQukl1+WBg0K+SXp2y70jVCxTio8JTs723QEuIi+PWrbNmdh/h07pGXLnKv5w4C+7ULfCBVDKgDgfz76SOrUSapSxVliqmlT04kAWIrD/QAAx9tvO3tNf/tbae1aBlQARjGkAgCkqVOl66+XLr/8hNZABYBwY0iFp6SkpJiOABfRtwccPCiNGCHde6/zk5Mj1agRkbeib7vQN0LFOanwFO5QYhf6Nmz/fqlXL+fq/WeflQYPluLiIvZ29G0X+kaoGFLhKd27dzcdAS6ib4O++865e9SGDdIrrzh3k4ow+rYLfSNUDKkAYJvNm6XOnZ21UFeskNq2NZ0IAMrhnFQAsMmqVVJyshQfL61bx4AKwLMYUuEpq1evNh0BLqJvl82cKXXoICUlOUtMnXOOq29P33ahb4SKIRWekp6ebjoCXETfLjl4UBo50rlI6q67pLfekmrXdj0GfduFvhEqzkmFp8ydO9d0BLiIvl2wb58znL7yivTMM9JDD0X0Cv5joW+70DdCxZAKT0lMTDQdAS6i7wjbtcu5gv/jj6X586WbbzYah77tQt8IFUMqAMSijRsln08qKZHefVe6+GLTiQAgKJyTCgCxJidHatdOOv10af16BlQAUYkhFZ4ydOhQ0xHgIvoOs0BASkuTbrlF6tTJ2YNav77pVKXo2y70jVAxpMJTGjRoYDoCXETfYVRU5FwgNXy4cyV/drbksXMC6dsu9I1QxQUCgYDpEEfj9/vVqlUr5ebmKikpyXQcAPCmf//b2XuamytNny796U+mEwGwSKTmNS6cAoBo9uGH0o03SgcOSCtXOneTAoAYwOF+AIhWr7ziXCB1xhnShg0MqABiCkMqPCUvL890BLiIvk/QwYPSqFHSbbc566B67AKpo6Fvu9A3QsWQCk9JTU01HQEuou8T8J//SLfeKj35pPT009KcOVK1aqZTVQh924W+ESrOSYWnTJ482XQEuIi+g/TVV86e0y1bpNdek264wXSioNC3XegboWJPKjyFJUvsQt9BWLZMat1aKiyUPvgg6gZUib5tQ98IFUMqAHhZICCNHy9dc41z56h//EO68ELTqQAg4hhSAcCr9u2T7rxTGjJEevhhadEi6Te/MZ0KAFzBkApPSUtLMx0BLqLvY8jPl9q3l+bPl7KynNudVq5sOlVI6Nsu9I1QceEUPKWwsNB0BLiIvo/i3XelLl2cq/bXrJFatjSdKCzo2y70jVBxW1QA8IpAQJo4UXroIenSS6WXX5bq1DGdCgCOKVLzGof7AcALfv5Z6tFDevBB52fpUgZUAFbjcD8AmPbll9Itt0hffCHNnSt162Y6EQAYx55UeEpBQYHpCHCRTX1/++23SklJ0YUXXqjzzjtPF154oVJSUvTj3LnO0lKFhdK6dTE9oNrUN+gboWNIhaf07t3bdAS4yJa+d+/erbZt22rGjBn69NNPtXnzZn326ac6c8YMndy9u4ratHHWP23WzHTUiLKlbzjoG6HicD88ZcyYMaYjwEW29D1s2DB9+eWXpf/8G0kvSrpW0mOSttWtq+m1axtK5x5b+oaDvhEqhlR4Cqs42MWWvtevX1/690mS5ks6WVInSYslXfCPf5gJ5jJb+oaDvhEqDvcDQIQVFxdLkvpIWiPpWznD6uLDfg8A+B+GVACIsOqVKilD0jRJmZLaS8r/1e/j4zmoBQCHY0iFp2RkZJiOABdZ0fcXXyjn3/9Wd0l3Suon6ZfDHtKmTRv3cxlgRd8oRd8IFUMqPMXv95uOABfFfN/z50tJSapXu7a6nnmmZh/hIY0bN1Z6errr0UyI+b5RBn0jVNwWFQDC7ZdfpNRU6fnnpa5dpWnT9G1RkVJTU7V+/XoVFxcrPj5ebdq0UXp6uupwZykAUSxS8xonQgFAOH39tXTbbdKHH0qTJkn9+0txcaojKTMz03Q6AIgaDKkAEC5vvin17CnVrCmtXi1Zcq4pAEQC56QCQKgOHHAO719/vdSuneT3M6ACQIgYUuEpPp/PdAS4KCb6/vpr6bLLpPHjpWeekV57TTrlFNOpPCkm+kaF0TdCxeF+eMqAAQNMR4CLor7vnBwpJUWqVUt67z0pOdl0Ik+L+r4RFPpGqNiTCk/p2LGj6QhwUdT2XVQkDRok3XyzdOWVzkVSDKjHFbV944TQN0LFnlQACMYXX0jdukmffFLm6n0AQHixJxUAKmrOHKllS+mnn6S1a6UBAxhQASBCGFLhKTk5OaYjwEVR0/fevdJdd0k9ejiH+P1+iRuMBC1q+kZY0DdCxZAKT8nKyjIdAS6Kir5zc52B9NVXpVmznJ+TTzadKipFRd8IG/pGqBhS4SnZ2dmmI8BFnu774EHp2Weltm2dq/f9fmehfpwwT/eNsKNvhIohFQAOt3On1KmT9PDDzlX8778v/e53plMBgFWCGlKLioo0bNgw1atXT4mJiUpOTtayZcuO+7xdu3Zp+PDhuvLKK3XyySerUqVKWrVq1QmHBoCIef116aKLpI0bpbfechboP+kk06kAwDpBDam9evXS+PHj1bNnT02cOFGVK1dWp06dtGbNmmM+Ly8vT+np6dq5c6cuuugiSVIcV8QC8JLCQun++6Ubb3QO8f/zn9K115pOBQDWqvCQun79emVnZ2vs2LFKS0vT3XffrRUrVujss89WamrqMZ978cUXa8+ePcrLy9PgwYNDDo3YlZKSYjoCXOSZvj/8UGrVSpo5U/rrX51bm9apYzpVzPFM33AFfSNUFR5S582bp/j4ePXt27d0W0JCgvr06aO1a9dq+/btR31ujRo1VLt27dCSwgrcocQuxvs+eFAaN076wx+kqlWdK/nvu4+1TyPEeN9wFX0jVBUeUj/88EM1adJENWrUKLO9devWkqSNGzeGNxms1L17d9MR4CKjfX/9tXTVVdLQodLAgdIHH0jnn28ujwX4ftuFvhGqCt8WdefOnapbt2657Ye27dixI3ypACBSAgHnzlH9+ztLSy1fLv3xj6ZTAQAOU+E9qfv27VNCQkK57VWrVi39PQB42nffSd26Oeud+nzSxx8zoAKAR1V4SK1WrZqKiorKbd+/f3/p7yOlU6dO8vl8ZX7atm1b7pZrS5Yskc/nK/f8/v37KyMjo8w2v98vn8+nghTfJ+cAAB5jSURBVIKCMttHjx6ttLS0Mtvy8/Pl8/mUl5dXZvukSZM0dOjQMtsKCwvl8/m0evXqMtuzsrKOeBJ5t27d+By/+hwLFy6Mic8RK31E+nMc+qsrn2PxYun3v5eWLdPoW29VWrNm0q/OlaePyH+OjIyMmPgcsdJHpD/HCy+8EBOfI1b6CNfnyMrKKp3FGjVqpBYtWkTuovhABV199dWBCy64oNz2ZcuWBeLi4gILFy6s0Ou88sorgbi4uMCqVauO+9jc3NyApEBubm5FYyLK3XDDDaYjwEWu9L13byDQr18gIAUCHTsGAt98E/n3xBHx/bYLfdsjUvNahfektmzZUps3b9bevXvLbF+3bp0kqUWLFuGcnWGpuXPnmo4AF0W879WrpRYtpBkzpMmTpbfflurXj+x74qj4ftuFvhGqCg+pXbp0UUlJiaZOnVq6raioSJmZmUpOTlb9//6Lf9euXcrLy1NxcXH40yLmJSYmmo4AF0Ws7/37nVuaXnaZ9NvfSh995FwoxdJSRvH9tgt9I1QVvrq/TZs26tq1q0aMGKHdu3ercePGmjlzpvLz85WZmVn6uOHDh2vWrFnaunWrGjRoULr9ySeflCRt2rRJkjRr1iy9++67kqSRI0eG5cMAgDZskO68U/rySyk9XRo8WKpc2XQqAECQKjykSs5g+eijj2r27Nn6/vvv1bx5cy1cuFDt27cvfUxcXNwRb3k6atQoxcXFKRAIKC4uTtOnTy99PEMqgJD98ov05JPSX/4iNW8u+f3ShReaTgUAOEEVPtwvOXeYSk9P144dO7Rv3z598MEH6tChQ5nHZGZmqqSkpMxeVEk6ePCgSkpKyvz10N8Dhxx+VSJiW9j6zs2VLr5Yevpp6dFHnYX5GVA9h++3XegboQpqSAUi7fA/3CC2hdx3UZE0cqRzW9PKlZ1D/aNHS1WqhCcgworvt13oG6GKCwQCAdMhjsbv96tVq1bKzc1VUlKS6TgAvGTDBiklRfr8c2fv6fDhDKcAYECk5jX2pAKILvv3S488IiUnO0Pphg3OkMqACgAxJagLpwDAqDVrpD59pK++cg7rs/cUAGIWe1LhKYff4g2xrcJ9/+c/0sCB0qWXSr/5jfThh+w9jUJ8v+1C3wgVQyo8JTU11XQEuKhCfS9e7Fypn5EhjR/v3EWKK/ejEt9vu9A3QsWQCk+ZPHmy6Qhw0TH7/u47qVcv6dprpd/9TvrnP6VBg1iYP4rx/bYLfSNUnJMKT2HJErscse9AQHrpJenBB6XiYmnaNKl3b25pGgP4ftuFvhEq9qQC8I4tW6TrrpN69JD++Efps8+cC6UYUAHAOgypAMwrLpbGjXPONf30U+mNN6TsbOm3vzWdDABgCEMqPCUtLc10BLgoLS1N+sc/pDZtpGHDpHvvdYbU6683HQ0RwPfbLvSNUDGkwlMKCwtNR4BbfvxRhfPmObc0PXhQ+uAD5+r9GjVMJ0OE8P22C30jVAyp8JTHHnvMdAREWiAgzZ0rNW2qx/LypOeec+4a1bq16WSIML7fdqFvhIohFYB7vvhCuuYaqXt36ZJLnAujHnxQimehEQBAWQypACJv/37p8celZs2kzZudC6Pmz5fOPNN0MgCARzGkwlMKCgpMR0C4LVrkXLX/5JPS4MFlLoyib7vQt13oG6FiSIWn9O7d23QEhMvWrdJNN0mdO0vnnCN9/LH09NNSYmLpQ+jbLvRtF/pGqBhS4SljxowxHQGhKiqSnnpKuuAC54Kol1+WliyRmjYt91D6tgt924W+ESquVoCnJCUlmY6AULz5pnMh1NatzqH9UaOOuaQUfduFvu1C3wgVe1IBhG7zZuew/vXXS2efLX30kZSezpqnAIATxpAK4MTt3SsNH+5ctf/pp9Krr0pLlzqH+gEACAFDKjwlIyPDdARUxMGD0osvSuedJz3/vDRypDOk3nyzFBdX4Zehb7vQt13oG6FiSIWn+P1+0xFwPB984CzE37On1K6dlJfnnHtarVrQL0XfdqFvu9A3QsWQCk+ZMmWK6Qg4mm3bpDvukNq2da7gf+f/27v3qKqqPA7g33sHBQQhFRS0FEWpJAVFSRDTzNBqyDQXPsbKB0JNPpvA0mlqtGktTSdH01VksnxkWqhoptkg+SoGVHAcE8JIfAGGpQbyvLDnjz3cJIS5cO6958D5ftY6C9twvL+7vuH9rXPO3vsr4NNP5TOozcS89YV56wvzJqU4u5+IGnfrlpwE9fbbgJsbsH49MG0a8LvfqV0ZERG1YmxSiejOap87XbQIKCoCXnoJePVV2agSERHZGG/3E1F9KSlAUBDw3HPy9n5Wltwtig0qERHZCZtU0pQnn3xS7RL07exZudbpI48ATk7A11/L50579bLJyzFvfWHe+sK8SSk2qaQps2fPVrsEfbp6FXj+eaBfP9mofvIJ8M03cha/DTFvfWHe+sK8SSk+k0qaEh4ernYJ+lJcDKxcCaxYAbRtK7/+8Y+Ao6NdXp556wvz1hfmTUqxSSXSo8pKID4eWLIE+OUXYO5cuXNUx45qV0ZERASAt/uJ9KWmBti2Dbj/fmDePPn86blzcokpNqhERKQhbFJJU5KSktQuoXUSAjhwABg8GJg8GfD3B/79b2DDBuCee1Qri3nrC/PWF+ZNSrFJJU35+OOP1S6h9Tl2DBgxAhgzRm5devgwsGcP8MADalfGvHWGeesL8yal2KSSpmzfvl3tElqPzEzgiSeAYcOAmzeBzz8Hjh4FHnpI7crMmLe+MG99Yd6kFJtUotYmKwuYOBEYOBD4/ntg+3YgIwN4/HHAYFC7OiIiIouwSSVqLXJygKlT5fOmqanA+vXAt98CkZGAkb/qRETUsvCTi6il++EHYPp0OWP/0CFg7Vo5Y3/mTMCBq8wREVHLxCaVNGX69Olql9By5OUBs2YB994LfPEF8M478vb+Cy/YbTF+pZi3vjBvfWHepBQvs5CmcIcSC+TmAm+9BWzaBHToACxbJrc0bddO7cqajHnrC/PWF+ZNSrFJJU2ZPHmy2iVoV06ObE63bAE8PGRzGhMDuLioXVmzMW99Yd76wrxJKTapRFqXlQX87W/Axx8DXl7A3/8ub/M7O6tdGRERkc2wSSXSqpMn5ZXTXbuAbt2A1avlZCgnJ7UrIyIisjlOnCJNOXbsmNolqO/IEbk71KBBcuvS+Hg5IerFF1tdg8q89YV56wvzJqXYpJKmLF++XO0S1CEEsG+f3B1q+HAgP1/e3s/OBqKiWsxs/abSbd46xbz1hXmTUmxSSVO2bdumdgn2VVUFbN4M9O8vtzCtqgL27JFXUCdNavXrnOoub51j3vrCvEkpNqmkKe1a4DJKzVJcLNc19fUFnn0W6NEDOHxY7hQVEaGb7Ut1kzcBYN56w7xJqdZ9mYZIawoLgXfflbtClZQAU6YAL78M9OundmVERESawiaVyB5On5ZXTrduBdq2BaKjgfnzgXvuUbsyIiIiTeLtftKU2NhYtUuwnpoaYP9+4NFHgYAAIDlZrnd66RKwciUbVLSyvOn/Yt76wrxJKV5JJU3p3r272iUod+sW8NFHwKpVciH+QYPkTP2nnwbatFG7Ok1pFXmTxZi3vjBvUsoghBBqF9GQjIwMBAUF4eTJkxg4cKDa5RA1Li8PWLcOWL8euHkTGDsWeOklYOhQ3UyEIiIi/bFVv8YrqURKCCFn5a9eDezeDbi5yXVNX3wR8PFRuzoiIqIWi00qUXMUFwNbtsgrp2fOAH37yj9PnQq4uKhdHRERUYvHiVOkKdnZ2WqX0LhvvwVmzwa6dZNfe/cG/vlP2ajGxLBBbSLN501Wxbz1hXmTUmxSSVPi4uLs+npFRUWYPn06/P39ce+998Lf3x/Tp09HUVHRrz9UWQl88gkwYgTwwANAYiIwb558BnXXLmDUKD5z2kz2zpvUxbz1hXmTUrzdT5ry7rvv2u21fvzxR4SGhiI3N7fO+NmzZ3H06FGkf/QROu7cCSQkAEVFwEMPAdu2AePGybVOSTF75k3qY976wrxJKTappCn2XLJk4cKF9RrUtgCeAhCTm4uOQ4YAHToAzz0HzJolnzslq+ISNfrCvPWFeZNSbFJJt9LT081/9gcwA8AzADwBHAYQ160bln//PeDkpE6BREREOsZnUkm3nCoqEA0gDcAZAFMBbAZwP4ARAHa7uLBBJSIiUgmbVNKUZcuW2fYFqquBgweBZ57B1z/8gHUArgIYB+BuAH8CUDsf1cGBNxpszeZ5k6Ywb31h3qQUP4VJU0pLS23zF2dlAZs2ybVNL18G+vTBZwMGYF5GBgoaOCU4ONg2tZCZzfImTWLe+sK8SSlui0qt17Vrcjb+xo3AiRNyEtSkScCzzwIPPoiia9cQEhJSb/IUAPj6+iI1NRWenp4qFE5ERNRycFtUIkuUlMjtSbduBb78Uo49/rhc2/T3vwccHc0/6unpidTUVMTFxSE9PR0mkwkODg4IDg7G8uXL2aASERGpiE0qtXyVlbIh3bpVNqilpcDQocCqVUBkJNBIs+np6YmEhAQ7FktERESWYJNKmnLt2jV4eHj8/x80mYCvvgK2b5e7Pv38s9wN6rXX5C19Hx+b10rKWZw3tQrMW1+YNynF2f2kKTNmzGj4myaTnJkfEwN4eQHh4cChQ/K/T58G/vMf4JVX2KC2II3mTa0O89YX5k1K8Uoqacobb7xRd6CyEkhJAXbuBJKS5PakPXsCUVHyVv6AAYDBoEqtpFy9vKlVY976wrxJKTappCkDBw4Ebt0CDhyQjenevcDNm4CvLzBtmmxMg4LYmLYSXLVDX5i3vjBvUopNKmnDjz8Cn38O7NkjG9SyMqB/f2D+fGD8eKBfPzamREREOsImldQhBJCdLZvSPXuA1FQ5HhoKvPGGbEx791a1RCIiIlIPJ06R/ZSVAV98AcydC/TpA/TtCyxZAnTpAmzYABQW4sPp04G4ODaoOvHhhx+qXQLZEfPWF+ZNSrFJJdvKywPWrZML6XfqBDz2mFzL9NFHgc8+k7tC7dwpnzft3BkZGRlqV0x2xLz1hXnrC/MmpbgtKllXcbFcFurLL+WRkwM4OABhYXLnp8cfl1dQ+XwpERFRq8BtUUmbTCYgIwNITpZN6TffAFVVcq3S0aOBt94CRo0C3N3VrpSIiIhaEDap1DQ1NcCZM3Lt0pQU4PBh4JdfgPbtgZEj5Vak4eFyySheLSUiIqJmYpNKjauuljs5HTkij8OH5XOkjo7A0KFyktMjj8i1S9u0UbtaIiIiaiU4cYrqqqgA/vUv4O23gYgIwMND7uoUGyvXMo2JkVdQb9yQW5QuXgwMGWK1BvXJJ5+0yt9DLQPz1hfmrS/Mm5TilVS9u3pVrlH6zTfyOHFCNqouLnLN0j/9CXjoISA4GHBysnk5s2fPtvlrkHYwb31h3vrCvEkpNql6UlYGnDoFpKfLIy0NyM2V37v7bnn7PjISCAkBAgNVuX0fHh5u99ck9TBvfWHe+sK8SSk2qa1VZSXw7bdy5v2JE7IpPX1azsZ3dAQGDgSeeEJeLQ0NBe65R+2KiYiIiMzYpLYGxcVyxv2pU7IpzciQk52qqgCjEbjvPuDBB4FZs+Rt+379OMmJiIiINK1JE6cqKiqwcOFCdO3aFe3atcOQIUOQnJxs0bk3btxAdHQ0PD094erqipEjRyIzM7NZRetWdTVw7pzcoal2f3tfX8DNTV4NnTtXXjUNDATeeUc+Y/rLL/KK6oYNwPPPyyuoGm5Qk5KS1C6B7Ih56wvz1hfmTUo1acepyZMnY8eOHViwYAH69OmDhIQEHD9+HF999RWGDh3a4Hk1NTUYNmwYTp8+jbi4OHTq1Anr1q3DpUuXcPLkSfRuYJ923e44VVEhnxXNygLOnv31+O47+T0A6NwZCAgA+vf/9bj/fnkrvwULCQlBamqq2mWQnTBvfWHe+sK89UP1HafS09Oxfft2rFixAi+99BIA4JlnnsEDDzyAuLg4fP311w2em5iYiNTUVCQmJmL8+PEAgMjISPj5+eH111/HRx99pPBttEAmE3DpkmxGv/tObh9ae+TlyUXzAaBjR8DfX14pjYqSjWi/fkCXLqqWbyuenp5ql0B2xLz1hXnrC/MmpSxuUhMTE+Hg4IDo6GjzmKOjI2bOnIlFixbhypUr6NatW4Pnenl5mRtUAPDw8EBkZCS2bNmCqqoqtGnkFvSECRMwfPhwLF++vOX8T19TAxQWAhcvyuP8eeCHH349LlyQt+8Befvd1xfw85O38P38gD595B73np7cuYmIiIh0x+ImNTMzE35+fnB1da0zPnjwYADAqVOnGmxSMzMz73j5d/DgwYiPj0dOTg78/f0bfO3z58/j/PnzOHr0KFJTU9VvVMvLgYICeeTn/3pcuSKvjl68CFy+LCcu1XJzk41or17A00/Lr7WHjw/gwDlsRERERLUs7owKCgrg7e1db7x2LD8/v9FzR4wY0ei5jTWptXJzcxEXF4eEhAQLq7aAySR3T/r55/rHtWtyl6XfHjdv1v07nJyArl0Bb2+ge3e5A1P37nUPd3deESUiIiKykMVNallZGRzvMCnH6X+7EJWVlTV4bnl5ebPOrR3vC6A95FIEv+zfj4xVq+St8qoq2WTe/rWyUl7p/O1RVgaUlgK3bgElJfJraan83p20bSsby44dfz169QI6dJB/9vT89WjfvuEG1GSSt/fJIunp6cjIyFC7DLIT5q0vzFtfmLd+ZGVlAWi8F2wOi5tUZ2dnVNTOLL9N+f+aPGdnZ6ufm5eXBwA4e/vg1avYuWCBZUUrUVkJFBXJg+wqKChI7RLIjpi3vjBvfWHe+pKXl9foak9NZXGT6u3tfcdb+gUFBQCArl27Wv3c0aNHY8uWLfDx8Wm0CSYiIiIidZSVlSEvLw+jR4+26t9rcZM6YMAAHDp0CMXFxWjfvr15PC0tDQAQGBjY4LmBgYE4evQohBAw3HZbPC0tDS4uLvDz87vjeR4eHvjDH/5gaYlEREREpAJrXkGtZfGOUxMmTEB1dTXi4+PNYxUVFUhISMCQIUPMM/sLCwuRnZ0Nk8lU59yrV69i586d5rFr167h008/RURERKPLTxERERGR/jRpx6mJEydi165dWLBgAXx9fbFx40acOHECBw8eRFhYGABg2rRp2LRpE/Ly8tC9e3cAcsepsLAwnDlzBrGxseYdpy5fvozjx4+jT58+tnl3RERERNQiNWlxzk2bNuG1117D5s2bcf36dQQEBGDv3r3mBhUADAZDnVv6AGA0GrFv3z7ExsZi9erVKCsrQ3BwMDZt2sQGlYiIiIjqadKVVCIiIiIie7D4mVQiIiIiIntRpUmtqKjAwoUL0bVrV7Rr1w5DhgxBcnKyRefeuHED0dHR8PT0hKurK0aOHInMzEwbV0xKNDfvgwcPYsaMGfDz84OLiwt8fX0xa9YsFBYW2qFqai4lv9+3mzVrFoxGIyIiImxQJVmL0ryTk5MxcuRI3HXXXXBzc8OgQYPwySef2LBiUkJJ3snJyXjkkUfQuXNntG/fHgEBAVizZg1qampsXDU1x61bt/D6669jzJgx6NixI4xGIzZu3Gjx+Vbp14QKJk2aJNq0aSPi4uLEBx98IEJDQ0WbNm3EsWPHGj2vurpahIaGCldXV7FkyRKxdu1a4e/vL9zc3MS5c+fsVD01VXPzDgoKEr6+vuKVV14RH374oVi0aJFwc3MTXl5eorCw0E7VU1M1N+/bHT9+XLRp00Y4OzuLiIgIG1ZLSinJe8OGDcJoNIoxY8aIdevWiffff18sWLBArFy50g6VU3M0N+/9+/cLg8Eg+vXrJ1atWiXi4+PFU089JQwGg5g3b56dqqemOH/+vDAYDMLHx0c8/PDDwmAwiI0bN1p0rrX6Nbs3qWlpacJgMNT5R6i8vFz07t1bhIaGNnru9u3bhcFgEDt27DCPFRUViQ4dOogpU6bYrGZqPiV5Hz16tN7YkSNHhMFgEH/+85+tXisppyTvWjU1NSIkJERERUUJHx8fNqkapiTv8+fPC2dnZzF//nxbl0lWoiTvKVOmCCcnJ3H9+vU648OHDxfu7u42qZeUqaioEFevXhVCCHHixIkmNanW6tfsfrs/MTERDg4OiI6ONo85Ojpi5syZSE1NxZUrVxo918vLC+PHjzePeXh4IDIyErt370ZVVZVNa6emU5L37atG1Bo2bBg6duyI7Oxsm9RLyijJu9bmzZtx9uxZvPnmmxCc16lpSvJ+7733IITAkiVLAAAlJSXMW+OU5O3s7AxHR0e4u7vXGffy8kK7du1sVjM1X9u2bdG5c2cAaPLvprX6Nbs3qZmZmfDz84Orq2ud8cGDBwMATp061ei5AwcOrDc+ePBglJaWIicnx7rFkmJK8r6TkpISFBcXw8PDw2o1kvUozbu4uBgLFy7EokWL0KVLF5vVSdahJO/k5GTcd9992Lt3L+6++264ubnBw8MDf/nLX9isapSSvOfMmYOamhrExMQgOzsbFy5cwHvvvYddu3bh1VdftWndZH/W6tfs3qQWFBTA29u73njtWH5+vk3OJXVYO7NVq1ahqqoKEydOtEp9ZF1K816yZAlcXFywYMECm9RH1qUk73PnzuHixYuYMWMGoqKisGPHDjz22GN48803sXjxYpvVTM2nJO+AgACkpKTgs88+Q9++fdGzZ0/MmTMHa9aswZw5c2xWM6nDWp/9TVrM3xrKysrg6OhYb9zJycn8/YaUl5c3+1xSh5K8f+vIkSP461//iokTJ2LEiBHWKpGsSEneOTk5WL16NbZt28atklsIJXnX3t5ftmwZYmNjAQDjxo3Dzz//jH/84x9YtGhRvSt2pC4leWdnZ+OJJ55Ajx498Pbbb8PJyQlbt27F7Nmz0aVLF4wdO9ZmdZP9Watfs/uVVGdnZ1RUVNQbLy8vN3/fFueSOqyVWXZ2NsaNG4f+/ftj/fr1Vq2RrEdJ3vPmzcPQoUMxbtw4m9VH1qX033ODwYDJkyfXGZ80aRLKysqa/CgQ2Z6SvF9++WU4ODjg0KFDmDp1KiZMmICdO3ciLCwML774Iqqrq21WN9mftT777d6kent73/Eyb0FBAQCga9euNjmX1GGNzC5duoTw8HB06NAB+/btg4uLi9XrJOtobt4pKSk4cOAA5s6di7y8PPNhMplQWlqKCxcuoLi42Ka1U9Mp+f2u/d5vnz2unahx/fp1a5VJVqIk72PHjmHkyJH1JklFREQgPz8fFy5csG6xpCpr9Wt2b1IHDBiAnJyceh84aWlpAIDAwMAGzw0MDERGRka9h+rT0tLg4uICPz8/6xdMiijJGwB++uknhIeHo6qqCgcOHOBkGo1rbt4XL14EAIwfPx69evUyH/n5+UhJSUHPnj2RkJBg2+KpyZT8fg8aNAhCCFy+fLnOeO0Hm6enp5WrJaWU5G0yme54tbR2lrfJZLJipaQ2q/VrFi9WZSW166ytWLHCPFa7zlpISIh5rKCgQGRlZYmqqirzWO26W4mJieaxoqIicdddd4nJkyfb5w1QkyjJu6SkRAQHBwt3d3eRkZFh17qpeZqb98WLF8Xu3bvrHElJSaJz584iODhY7N69W+Tm5tr9/VDjlPx+JyUlCYPBIBYvXmweq66uFmFhYcLDw0NUVlba502QxZTkHRYWJjp16iR++ukn85jJZBJBQUHC3d1dmEwm+7wJapbjx483uE6qLfs1VXacioyMNO9Y8f7774vQ0FDRtm3bOou3P/fcc8JgMIgLFy6Yx6qrq0VISIho3759nR0M3N3dRU5OjhpvhSzQ3LzHjh0rDAaDmDlzpti8eXOdIykpSY23QhZobt530qNHDy7mr3FK8h41apQwGo0iJiZGrF27Vjz66KPCYDCIDz74wN5vgyzU3Lz3798vjEaj6N27t1i+fLlYvXq1CAkJEQaDQbz11ltqvBWywJo1a8TSpUvFCy+8IAwGg3j66afF0qVLxdKlS8XNmzeFELbt11RpUsvLy0VsbKzw9vYWTk5O4sEHHxRffvllnZ+ZNm2aMBqN9f5Ru379uoiKihIeHh7CxcVFPPzww+LkyZP2LJ+aqLl5+/j4CKPRKAwGQ72jZ8+e9n4bZCElv9+/xR2ntE9J3iUlJWL+/PnC29tbODo6ioCAALF161Z7lk9NpCTvL774QgwbNky4uLiY846Pj7dn+dREPj4+5s9do9Fo/ky+PV9b9msGIbhqMhERERFpi90nThERERER/T9sUomIiIhIc9ikEhEREZHmsEklIiIiIs1hk0pEREREmsMmlYiIiIg0h00qEREREWkOm1QiIiIi0hw2qURERESkOWxSiYiIiEhz2KQSERERkeawSSUiIiIizfkvSxpGwj1R8tAAAAAASUVORK5CYII="
      ],
      "text/plain": [
       "Figure(PyObject <matplotlib.figure.Figure object at 0x10e1eb890>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i=1:N\n",
    "    plot(mo.no[\"coords\"][i][1], mo.no[\"displacement\"][i], \"ko\")\n",
    "end\n",
    "x = linspace(0, 1)\n",
    "plot(x, 1/2*x.^2, \"-r\", label=\"acc\")\n",
    "legend(loc=\"best\")\n",
    "grid()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.3.8",
   "language": "julia",
   "name": "julia-0.3"
  },
  "language_info": {
   "name": "julia",
   "version": "0.3.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}