permafrost.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "<div style=\"display: flex; align-items: center;\">\n",
    "    <img src=\"https://github.com/nagelt/Teaching_Scripts/raw/9d9e29ecca4b04eaf7397938eacbf116d37ddc93/Images/TUBAF_Logo_blau.png\" width=\"500\" height=\"auto\" height=\"auto\" style=\"margin-right: 100px;\" />\n",
    "    <div>\n",
    "        <p><strong>Prof. Dr. Thomas Nagel</strong></p>\n",
    "        <p>Chair of Soil Mechanics and Foundation Engineering<br>Geotechnical Institute<br>Technische Universität Bergakademie Freiberg.</p>\n",
    "        <p><a href=\"https://tu-freiberg.de/en/soilmechanics\">https://tu-freiberg.de/en/soilmechanics</a></p>\n",
    "    </div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#Some plot settings\n",
    "import plot_functions.plot_settings\n",
    "\n",
    "#FEM routines\n",
    "%run plot_functions/fem_routines.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Liquid-solid phase transition in soils: permafrost melting\n",
    "\n",
    "## Governing differential equation\n",
    "\n",
    "We consider the problem as one dimensional with piecewise constant properties:\n",
    "\n",
    "$$\n",
    "\t\\varrho \\dot{h} = \\lambda T_{,zz}\n",
    "$$\n",
    "\n",
    "$[\\lambda] = \\text{W m}^{-1} \\text{K}^{-1}$: heat conductivity\n",
    "\n",
    "$[\\varrho] = \\text{kg m}^{-3}$: density\n",
    "\n",
    "$[h] = \\text{J kg}^{-1}$: specific enthalpy\n",
    "\n",
    "Let the ice volume fraction be given as an equilibrium function\n",
    "\n",
    "$$\n",
    "    \\phi_\\text{I} = \\phi \\left[1 + e^{k(T - T_\\text{m})} \\right]^{-1}\n",
    "$$\n",
    "\n",
    "so that the apparent ice density is\n",
    "\n",
    "$$\n",
    "    \\varrho_\\text{I} = \\phi \\varrho_\\text{IR} \\left[1 + e^{k(T - T_\\text{m})} \\right]^{-1}\n",
    "$$\n",
    "\n",
    "With $h(T,\\varrho_\\text{I})$ we then find the energy balance\n",
    "\n",
    "$$\n",
    "    \\left[ C - L \\frac{\\partial \\varrho_\\text{I}}{\\partial T} \\right] \\dot{T} = \\lambda T_{,zz}\n",
    "$$\n",
    "\n",
    "with \n",
    "\n",
    "$[L] = \\text{J kg}^{-1}$: specific latent heat of fusion\n",
    "\n",
    "$[C] = \\text{J m}^{-3} \\text{K}^{-1}$: volumetric heat capacity.\n",
    "\n",
    "We further use\n",
    "\n",
    "\\begin{align}\n",
    "    C &= (1 - \\phi) c_{p\\text{S}} \\varrho_\\text{SR} + \\phi_\\text{I} c_{p\\text{I}} \\varrho_\\text{IR} + (\\phi - \\phi_\\text{I}) c_{p\\text{W}} \\varrho_\\text{WR} \n",
    "    \\\\\n",
    "    \\lambda &= (1 - \\phi) \\lambda_\\text{SR} + \\phi_\\text{I} \\lambda_\\text{IR} + (\\phi - \\phi_\\text{I}) \\lambda_\\text{WR} \n",
    "\\end{align}\n",
    "\n",
    "\n",
    "## Weak form\n",
    "\n",
    "The pore pressure can have (the essential/Dirichlet) boundary conditions in the form:\n",
    "\n",
    "$$\n",
    "    p = \\bar{p}\\ \\forall z \\in \\partial \\Omega_\\mathrm{D}\n",
    "$$\n",
    "\n",
    "We now introduce a test function $\\eta$ which vanishes where the pore pressure is given\n",
    "\n",
    "$$\n",
    "    \\eta = 0\\ \\forall z \\in \\partial \\Omega_\\mathrm{D}\n",
    "$$\n",
    "\n",
    "and construct the weak form (using integration by parts):\n",
    "\n",
    "\\begin{align}\n",
    "   0 &= \\int \\limits_0^H \\eta \\left[(C - L \\partial_T \\varrho_\\text{I}) \\dot{T} - \\lambda T_{,zz} \\right] \\text{d}z\n",
    "   \\\\\n",
    "   &= \\int \\limits_0^H \\left[\\eta (C - L \\partial_T \\varrho_\\text{I}) \\dot{T} - \\left( \\eta \\lambda T_{,z} \\right)_{,z} + \\eta_{,z} \\lambda T_{,z} \\right] \\, \\text{d}z\n",
    "   \\\\\n",
    "   &= \\int \\limits_0^H \\left[\\eta (C - L \\partial_T \\varrho_\\text{I}) \\dot{T} + \\eta_{,z} \\lambda T_{,z} \\right] \\, \\text{d}z + \\left[ \\eta q_z \\right]^H_0\n",
    "\\end{align}\n",
    "\n",
    "where the natural/Neumann boundary conditions have appeared."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{k e^{k \\left(T - T_{m}\\right)}}{\\left(e^{k \\left(T - T_{m}\\right)} + 1\\right)^{2}}$"
      ],
      "text/plain": [
       "-k*exp(k*(T - T_m))/(exp(k*(T - T_m)) + 1)**2"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import sympy as sp\n",
    "k,T,Tm = sp.symbols('k T T_m')\n",
    "f = 1/(1+sp.exp(k*(T-Tm)))\n",
    "sp.diff(f,T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def phi_I(T,phi=0.35,T_m=273.15):\n",
    "    return phi / (1 + np.exp(k*(T-(T_m-4/k))))\n",
    "\n",
    "def dphi_I_dT(T,phi=0.35,T_m=273.15):\n",
    "    temp = np.exp(k*(T-(T_m-4/k)))\n",
    "    return -phi*k * temp / (1 + temp)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0sAAAIrCAYAAAA+xwKYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACLqUlEQVR4nOzdeXxU13n/8c+ZRfu+ICSxG4MBG7PI+747MWA3jussTtwlQJM0vyRtmr0F0jSN2zRO3WYB0ixO0zTx0hhInNiJbbxjS2I1NjY7SAIkFu3bzJzfHzMjBBagZWbundH3/XopQnfuPecZIEaPznOeY6y1iIiIiIiIyKk8TgcgIiIiIiLiRkqWREREREREBqBkSUREREREZABKlkRERERERAagZElERERERGQAPqcDSAa33367/d3vfud0GCIiIiIiEh9moItaWRqEpqYmp0MQEZEEqaiooKKiwukwRETEBbSyJCIi0k9DQ4PTIYiIiEsoWRIREemnrq7O6RBERMQllCyJiIj0oxI8ERGJ0p4lERERERGRAShZEhER6WfJkiUsWbLE6TBERMQFjLXW6Rhcr6qqylZXVzsdhoiIJIAx4e6x+vdRRGRUGbB1uPYsiYiI9LNy5UqnQxCRJNHc3ExTUxM9PT1OhyID8Hq95ObmUlRURHp6+rDG0MrSIGhlSURERET66+rqYv/+/YwbN47MzMy+VWlxB2stvb29tLS0cPz4cSZMmHCuhEmH0oqIiIiIxEJjYyOlpaVkZWUpUXIhYwxpaWmUlJRQWFjIsWPHhjWOkiUREZF+1q5dy9q1a50OQ0Rcrquri5ycHKfDkEHIy8ujtbV1WM9qz5KIiEg/ixYtAtTgQUTOLhAI4PPpW+lk4Pf7CQaDw3pWf8IiIiL9LFiwwOkQRCRJqPwuOYzkz0nJkoiISD8qwRMRkSjtWRIRERERERmAkiUREREREZEBKFkSERHpxxijfQgiMmpt374dYwxPP/103Od69NFHufvuu5k4cSKZmZlMnz6dL33pS8PuXBcPju9ZMsaMAR4EqiKXtgKfsdYePMdzOcBfAwuAzMhHCPgR8KDt18bIGPMcMAY4/Xjlb1trH47B2xARERERSXq1tbUAzJ8/P+5zfetb32LChAl84xvfYNy4cWzcuJHly5fz7LPP8vLLL+PxOL+u42iyZIxJA54G3gZmAZZwsvOsMWautbbtLI9PAr4O3GWtXRcZbxHwayAr8lp/77XW7o1l/CIiknrUMlxERrOamhomT55MUVFR3Odau3YtpaWlfV9fd911FBUVcf/99/Pcc89x4403xj2Gc3E6XbsfmA18wVobsNYGgS8AU4CPn+PZduAH0UQJwFq7BtgCvC9O8YqIiIiIpKyamhqqqqr6vm5paeH9738/Y8eOZf369TGdq3+iFHXJJZcAUFdXF9O5hsvpMry7gf3W2t3RC9baQ8aY7ZHX/vVMD1pr9xAuwztdHvBOrAN12omjh+lobXY6jLgzJHifwFmmO1ss73rFvPuLUy+deq1vbGNOvWZO/vrUr8P/68H0PWOMJ/LZhJ/QFgtHDW2PS797Tf8/ZHNyrPAfbuTXHown/Gcf/bXxeDBeL8brjdE7EBGRkVr4t084HQIAa//tzmE9Z61l06ZNfPWrXwVg69at3H333RQXF1NTU0NlZeW77h/MYa/GGLyD/PcqmpDNmDFjiNHHh9PJ0mzCJXin2wPcNJSBjDHZwJeA6OfT/Y0x5nLCydQR4MfW2h8PLVznvPD9b1P0+k6nwxARtzEmnDT5fBivF4/Piyc9HW9mJt6sLHxZmXgzs/Bmhb9OKywks6KczMpK0svG4NHp8++ycOFCQOcticjo8/bbb9Pa2kpVVRUPP/wwH//4x7n//vv5zne+Q1pa2rvuX79+PTfccMM5x73uuut47rnnznlfXV0d//AP/8DNN998yuqWk5z+V7IEqBngeguQZYzJtNZ2nmsQY8yLwKXALuB91tra0245AewknER1AX8C/LcxZpa19nNnGHMJsARgwoQJg3s3cdTYBb6sM2fk5ixfnbo6MfDDbliUSPguAce3JdgBfjWEJ5NkX4UxkXUyY/AYD57I576v8eDxePBGrrnh7+JQDemP4pSbbd+nk/tkbPiytYCNfAphQxasxYZCEAphQyFsMBi+FghgAwEAzv3zvZOM10t6WRmZleHkKXfaNIourcLj9w9hlNSzbt26c98kIjKA4a7ouEVNTfjb8pUrV7J27Vp+8IMf8NGPfvSM98+fP5/XX3/9nOPm5uae8562tjbuvPNOfD4fP/6xe9YznE6WzmRI3y9Za682xviBDwJ/NMZ83lr7UL/X7zrtkUeNMTcAnzXGPGSt3T/AmKuAVQBVVVWOf1faO+/DfK9zT9zGT0/zkpPpJzvTT06mn5zMNHKy/ORmpZGfk0ZedvQjve/XuVlpeDzJ+K1t6rHWYrHhz9YSwhKyofCvbajfR/jrYOTrYChIMBQkZEME+n0OhAKRjyC9oV56g+Gve4MBekO9dAd66A720hPooTsY+Qh00xXopqOnk47ekx9BGxr0+8j0ZzAmu4Qx2cWUZZcwJqeEcXljmVQ4npy07Dj+DiYvGwxig0FCgSA2GE6aQj09BDs6CXR0EOzsJNjeEf51Rwc9R4/SWVdPZ3093Y1NdNXX01Vfz/HXw/9A+vPzKbvlJspuvYWMsjEOvztnrFmzxukQREQcUVtbS0FBAY899hif+tSnzpooAeTk5DBnzpxzjnuuUvWuri4WLVrE7t27Wb9+PePGjRtK2HFlnOz6Y4ypB9621l5/2vU1wE3W2iF/d2SMWQn8OTDBWnvoLPd9FPgp8H5r7WNnG7OqqspWV1cPNZSYs9bSGwjR0xukO/LR0xv5uif8dVdPgK7uyOeeU7/u6ArQ2R2go6s38jn80d7VSyg09L8HXo+hIDedwrwMCnPTKcrLoDA3g6K8dEoKMiktzKK0IJPszNH9U+rRzFpLT7CXjt5OWrvbaO5upbmrheau1sivw183dRznSHsTXYHuM45Vml3M5MLxTC4Yz+TCCUwpHE9BZn4C303qCXZ309VwiM76ejoP1tH04kt07Iv87MgYCufPY+x7bqNw7hztjRIROc2bb77pmn01sXLDDTdQVlbGtddey6c+9SkeeeQR3ve+M/dNe+6550Zchtfb28tdd93F+vXr+cMf/sDll18+3PDPahB/XgNmdE6vLG0BLhjg+mTC5y2dUaTteMhaGzjtpc2AH7gIOBS5L9Nae3p3hGi1StJ8B2CMIc3vJc3vJSeG41pr6eoJ0tbRS1tnD+2dvbR19tLW0UtLew8t7d2Rzz19X59oC993tLmLo81dZx0/K8NHaTR5KsykvDibipJsykuyGVucTZo/af4IZIiMMaT70kj3pVF4jsTGWktrdxtH2o9ypL2JI+1HOdzWxP4TB9nXXEdj+1Ea24/y2sFNfc+UZZcwt/xC5lbMYlbpNNJ8766nljPzpqeTPWki2ZMmAjDunrtpfWsHDb/9HUdffoXj1TUcr64hfcwYKt93J2Nvv02HtYqIpChrLRs3buSLX/win/jEJ9ixYwf33Xcf69ev7+tQd7qRluGFQiE+/OEP88c//pHf/OY3cUuURsLpZOlxYKUxZlL0DCRjTBkwg9OaNESuN1rbV9PzZaAb+OfTxpwU+Xw08vnKyFi3nXZf9KStjSN7C8nPGENmuo/MdB+lhZmDfq6nN8iJ1m6OtXZxvKWb45HPx1q6aDrRyZHjHTSe6KSjK8C+Q63sO/Tu05iNgdKCTMpLsqkoyWHC2FwmlucxcWweedn6xnc0McaQl5FLXkYuU4snnfJaMBSkvvUwe44fiHzsZ8+JAxxub+J3O5/jdzufw+/1c+GYaeHkqXwWZTnvbkcqZ2eMIW/GBeTNuIDej/05h//wDId+/xTdh4+w+werCXZ2Me59dzkdZtytWrUKgCVLljgciYhI4uzatYvm5ua+sroHH3yQ3bt3s2jRIjZs2DDgHv7c3NwRNWL45Cc/ySOPPMJXvvIVsrOzefXVV/teGzdunCvK8Zwuw0sDqoE3gQ8DIeC/gKuBvkNpjTFXAc8Dq6y1H49cWw58DLjVWrs9cu1q4EmgFrjeWmuNMdcDfwDutNb+JnLf9cBvgMestWcvxsQ9ZXjJyFpLa0cvjZHE6cixDhqa2qlvaqehqZ3DxzvOWAJYlJfOxLF5fcnT+eMLGFeWi1f7pITwT6N2HtvLxoY32Niwjd3HT916OLlwPLdPvZ6rJl5CmleloMNlQyEO/+EZdn33+wCc/+m/ZsyN5y65SGbR1TMdTisiZ5NqZXi//OUv+cAHPkBDQwNjx44Fwk0Xrr76aoLBIC+99BJ5eXkxnXPSpEns27dvwNeWLVvG8uXLYzbXcMvwHE2WoG/F6EGginAPqG3AZ6y1B/rdczHwLPAda+3XItcmA39FeMXIEF4lCwL/CzzUL9HKI5xU3Q3kEG4t3gP8DPjXAcr43kXJUvz0BkIcOR5OoA4eaWP/oRb2HWph36FWunve3dcrM93L1HGFTJtQwPkTCpk+oZDi/AyVBgknOpvZdGg7GxveYPOh7XT0hhtp5qZlc9N5V3PreddSkh3/08hTVf3a37Dnhz8Cj4cZX/4CRZe4o6VrPERXlKIrTCIiA0m1ZCnVJW2ylAyULCVeKGQ5fKwjnDg1tLCnvoV3DhznyPF3d5Ivyktn1pQSLj6/hNlTSxlbnKXkaZTrCfby8v5qnnznWfYcD//cxRjDJZUX857zb2Bm6fn6OzIM+372cw4++jietDRmfW0ZeTMG2nIqIjI6KFlKLkqW4kjJknscb+3inf0neHv/8fDHgRO0d/aeck9pYSazp4YTp9lTSygpGPw+LEkt1lrePrqbJ995jg0HavvamJ9fPJklVR9iYoHztdDJxFrLzv/8Pkf+8Ed8OTlc9M//SJYLzqETEXGCkqXkomQpjpQsuVcoZKlrbGPLzia27Gxk684mWjtOTZ4mV+Rx+YXlXH5hOZMr8rSiMEod6zzB0ztf4Oldz9PS3YbXeFgw/WbeP+sO0tVFb9BsMMhbD/wrxza8TlpxEbMf+AbppanVTKO+vh6AiooKhyMRETdTspRclCzFkZKl5BEKWfY2tLBlZyOb32nijd1NdHaf3Ps0pjCzL3GaObkIr9fjYLTihI6eTn6x9Qme2vk8FsuY7GI+Nv9DzCmf6XRoSSPY3c325f9Iy/Y3yRxXyUX//HX8Md706yQ1eBCRwVCylFyULMWRkqXk1RsIsmVnE69uO8SGbQ0cbz156Glulp+rL67k5ksncP74Aq04jTJvN+1mVfX/sL+5DoCrJlRx/9x7KMhInW/64ynQ1sbWL/89Hfv2k3vBdC76569jPKnxw4foilJ0hUlEZCBKlpKLkqU4UrKUGkIhy9sHjvPq1gZe3dZAXWN732vjy3K5+ZIJ3DB/HIV5GQ5GKYkUCAX5zY4/8sgb6+gJ9pKdlsVfzL2XayZd6nRoSaH76DE2/83f0XviBBd+42vkz5rldEgiIgmjZCm5KFmKIyVLqWlPfTN/fP0Az9UeoLmtBwCPxzD/gjHcfMkELp01Fp/K9EaFw22N/LDmf9l8aDsAH7zoTu6acZtWGwdh709/Rt3jv6bstluZ+omlTocjIpIwSpaSi5KlOFKylNoCwRDVbx7mD6/tp/rNwwQjh+SWFGRy57XncetlE8jK0KGmqc5ay+93rufHtb/CYrlj2k18ZM778BglzGfTvncfmz79N/hyc7jkxz/E49f/V0RkdFCylFyGmyz54hOOSPLweT19TR9OtHbzXO0BfvfKPuoa2/ivNdv436fe4j1XTmbhNVMoUoleyjLGcPv515OXnst/bPgxv3n7j7R2t/FXl34En8frdHiulT1pIlkTJ9Cxbz8nNm6i6NJLnA5pxObPnw9ATU2Nw5GIiIjT9CNTkX4KctO567qpfO/zN/L3f3EZs6YU094V4NFn3uEvv/40D/1yIwcOtzodpsTRlRPm86VrPkm6L53n923gWy/+gO5Aj9NhuVrptdcA0Pj8Cw5HEhu1tbXU1tY6HYaIiLiAkiWRAXg8hktnjeWbn7yab/2/a7hydjnBUIinX9vPJ/7lGR54+HUamtrPPZAkpdljZ7Ds+s+Qm5ZNbcM2vr7+Idp69Od9JiXXXA3AsQ2vE+zsdDiakauurkal1yIiAkqWRM5p+sQivnT/pfzgizfxnisnkebz8OLmej7xL39k9a+30tKuVYdUNLV4Eitu+luKMwvZ0bSL5c88yLHOE06H5UoZZWPInXEBoZ4ejm54zelwRmz+/Pl9pXgiIjK6KVkSGaSKkhw+cffFrPzSzdx0yXiCIcuaF3az5BtP89gz79DTGzz3IJJUxuWV8483f47K3LHsb67ja89+h87eLqfDcqVoKV5TipTiiYiMVtu3b8cYw9NPP53wuW+//XaMMXz1q19N+NxnomRJZIhKCjL5zAfm8e9/cz1zp5XS3hXgJ7/Zzl898EeeqT5AKKQOk6mkJKuIFTf9LePzK6hvPcyq6p+jLqLvVnL1lRivl+MbN9Pb3Ox0OCOyfPlyli9f7nQYIiKOiO7ZTPQK+y9+8Qs2b96c0DkHQ8mSyDBNrsjna0uvZMWSK5hUnkfj8U4e/EUtX/zuixw8oiYQqSQvPYe/vXIxGb50XtpfzR92veh0SK7jz8ujYM7FEArR9OLLToczIitWrGDFihVOhyEi4oiamhomT55MUVFRwuY8ceIEn/3sZ/n2t7+dsDkHS8mSyAjNmz6G7/zN9Xz63rkU5qbz5t5jfPrfnuPxZ3f2ndkkya8ibyxLqj4MwE82/oo9xw84HJH7lKRIV7xly5axbNkyp8MQEXFETU0NVVVVfV+3tLTw/ve/n7Fjx7J+/fq4zPn5z3+eWbNm8cEPfjAu44+EzlkSiQGvx3DzpRO4/MKxrH5iG89UH+DH697g5S31/L975zBhbJ7TIUoMXD3xErY3vsMfdr3At19ezQO3fImstEynw3KN4ssuYVd6Oq1v7aDr8GEyysqcDmlYVIInIsP1p7/8uNMhAPCre78/rOestWzatKlvz9DWrVu5++67KS4upqamhsrKynfdHwyee8+2MQavd+AzC1988UUefvhhV5bggVaWRGIqJyuNz35wHss+djnF+Rns2H+cT397PY/88W2CwZDT4UkM/Nnce5hUMI7DbY384PX/1v6lfryZmRRdFj6UtumFlxyORkREhurtt9+mtbWVqqoqHn74YS6//HJuvvlm1q9f/65ECWD9+vX4/f5zftx0000Dztfb28vSpUv53Oc+x/Tp0+P99oZFK0sicVA1o4zv/t2N/NeabTz92n4e/u2bvLylns9+cJ5WmZJcmtfP31y5mC889c+8erCW373zHO+ZdoPTYblG6bXX0PT8izSuf57Ku/8EY4zTIQ1ZTU0NkPjNzSKS/Ia7ouMW0f/+rVy5krVr1/KDH/yAj370o2e8f/78+bz++uvnHDc3N3fA6w888ACdnZ185StfGV7ACaBkSSROsjP9/L9753L1nEr+85FN7DzYzN/++/N8+gNzufrid/90RpLH2NwxfPzSj/Dtl1fz8ObHOL94MlOLJzkdlisUzJ2DLzeXjv0H6Ni3j+xJk5wOaciitfpaNRSR0aa2tpaCggIee+wxPvWpT501UQLIyclhzpw55xx3oB+c7d+/n3/6p3/ihz/8Id3d3XR3d/e91t3dzYkTJ8jNzT1j+V6iqAxPJM7mTR/Df37uBq6bO46uniAPPFzNj9e+obK8JHf5+HncPvV6gqEgD77yQ9p62p0OyRU8Ph8lV10BQOP65Gz0MG/ePObNm+d0GCIiCVdTU8Ntt93GQw89xEMPPcTjjz9+1vtHUoa3e/duurq6uO+++ygsLOz7APjWt75FYWEhW7dujcv7HAqtLIkkQFaGn7/98DymTSzgR2ve4PHndrLz4Ak+/5Eq8nPSnQ5Phukjc97HO0f3sOv4Pla+/nP+9qolTofkCiXXXsOh3z1F4/MvMvEjH8Z4kuvnctEyFBGR0cRay8aNG/niF7/IJz7xCXbs2MF9993H+vXrueSSSwZ8ZiRleHPmzOHZZ5991/UbbriB++67j7/8y79k6tSpQ38jMaZkSSRBjDEsuuY8plTk88DPqtmys4nPPLieL//ZJZw/vtDp8GQY/F4/n73yY3zu919nw8GNbD/yNjPHTHM6LMflzbiA9NISuhubaHnzLfJnzXQ6JBEROYddu3bR3NzcV1b34IMPsnv3bhYtWsSGDRuYMGHCu57Jzc09pc34UBQUFHD99dcP+NrEiRPP+FqiJdeP+0RSwIXnlfCdz17H9ImFNJ3o5Av/+SJPb9jndFgyTGNySlh0wS0A/M+WJ7TPBTAeDyXXXA1AU5KfuSQiMlpEV9WjyZLH4+EXv/gFZWVl3HHHHbS0tDgYnXOULIk4oDg/k3/+xNW858pJ9AZCPPSrTaz+9VZCOsQ2Kd0x7Sby0nN4++huauq3OB2OK5ReFz6gtumllwkFAg5HMzQVFRVUVFQ4HYaISELde++9WGsZO3Zs37WcnBw2bdrE1q1byctLTDdfay1f//rXEzLXYChZEnGI3+fhE3dfzKfvnYvP62HNC7v5919uVOOHJJTpz+B9M98DwC+2PEEopD/D7EmTSC8tIdDaRtehQ06HMyQNDQ00NDQ4HYaIiLiAkiURh9186QSWfewyMtK8PFN9gG8+/Dq9gXOfhi3ucst511CaVcSBlgae37fB6XBcIXPcOAA66+odjmRo6urqqKurczoMERFxASVLIi4wZ9oY/nHplWRn+nl12yG+9sMNdHYnV+nSaOf3+vnTCxcC8Ktt6+gN9jockfMyI6VsyZYsqQxPRESilCyJuMQFk4r4509cRUFuOpveaeTvV75Ma0eP02HJEFwz8VLG51fQ1HGMp3epsUFmZXImSyIiIlFKlkRcZHJFPg988mpKCzPZse84X/7eSxxv6XI6LBkkj8fDBy+6E4DHtj9JZ+/o/rOLJktd9cmVLC1ZsoQlS3RmloiIKFkScZ2K0hwe+OQ1VJbmsLehhS9890UOH+twOiwZpPkVFzG9eAqt3W2s2/EHp8NxVGZlJQCdSbb/Z/Xq1axevdrpMERExAWULIm4UGlhJg/89dVMqcynoamdv//ByxxvHd2rFMnCGMOHLr4LgLU7/kBLV6uzATkorbgIT1oavc0tBNranA5n0FauXMnKlSudDkNERFxAyZKIS+XnpPONj1/FeePyaTjazoofvkpHl5oGJIMZpeczt/xCugLdPP7m75wOxzHG40nKfUsqwxMRkSglSyIulp3pZ9nHLqe8OJtdB5v555+8Tm9AZ/gkg+jepad2Pk9j+1GHo3FORpJ2xBMREQElSyKuV5ibwYolV1CQE+6S953/rSUUsk6HJecwqXAcV0+4hEAowK+2rXM6HMecXFlKnn1La9euZe3atU6HISIiLqBkSSQJlJdks2zx5WSme3l+Yx0/WvsG1iphcrs/vWghHuPhxX2vcaKrxelwHNGXLNU3OBzJ4C1atIhFixY5HYaIiLiAkiWRJDF1XAFf/rNL8XkNTzy/i/97bqfTIck5jM0pZV75hQRtiOf3bnA6HEckY0e8BQsWsGDBAqfDEBERF1CyJJJE5kwbw2c/OA+AH6/bzjPVBxyOSM7lhilXAvDsnpdH5Wpg31lLDYewoeTYb6cyPBEZzbZv344xhqeffjrucx08eJBPfepTXHHFFWRlZWGMYe/evWe8/8CBA7z//e8nPz+fvLw83ve+97F///64xqhkSSTJXDt3HB+780IAHvrlRmrfOuJwRHI2c8svJD89l7qWQ7xzdI/T4SScLysLf2EBoZ4euhubnA5HRETOoba2FoD58+fHfa6dO3fyq1/9isLCQq655pqz3tvR0cGNN97IW2+9xU9/+lN+9rOf8c4773DDDTfQ3t4etxiVLIkkoTuvPY+7b5hKMGT5l/+u5tDR+P1HQkbG5/Fy7aTLAHh2zysOR+OMzGhHvHp1xBMRcbuamhomT55MUVFR3Oe69tprOXz4ML/97W+55557znrv6tWr2b17N7/+9a+56667uPPOO1mzZg379u2L69l4SpZEktT9d8zkslljae/s5ZsPv05Pb9DpkOQMoqV4L++vpivQ7XA0iZdsZy0ZYzDGOB2GiIgjampqqKqq6vu6paWF97///YwdO5b169fHdC6PZ/CpyJo1a7j88suZOnVq37XJkydz1VVX8cQTT8Q0rv58cRtZROLKGMNnPjiPz3z7OXYdbGbVr7fy1/fMcTosGcC4vHKmFU/h7aO7efVALddPvsLpkBIqGZs8iIgMx0t33u10CABc9cRjw3rOWsumTZv46le/CsDWrVu5++67KS4upqamhsrIf8/73x8MnvuHtcYYvF7vsGKKeuONN7jzzjvfdX3WrFk88sgjIxr7bLSyJJLEcjL9fOn+S/D7PPz+1X08Ux3fTY4yfDdEEqTRWIqXbCtL1tpR2YxDROTtt9+mtbWVqqoqHn74YS6//HJuvvlm1q9f/65ECWD9+vX4/f5zftx0000jju3YsWMUFha+63pRURHHjx8f8fhnopUlkSR33rgClv7JbP7zkU1899EtTKksYFJ5ntNhyWmumDCfn2x8hDcb36Gh9QjluWOcDilhonuWurRnSURS3HBXdNyipqYGgJUrV7J27Vp+8IMf8NGPfvSM98+fP5/XX3/9nOPm5ubGJL6BSqTj/cMtJUsiKeDWyybw5t6j/PH1A/zzT17jwc9eR1aG3+mwpJ8sfyZXjJ/Pc3tf4dk9L/Oh2Xc5HVLCpJeNwXi9dDc2Eezuxpue7nRIIiIygNraWgoKCnjsscf41Kc+ddZECSAnJ4c5c+acc9xY7AMtLCzk2LFj77p+/PjxAVecYkVleCIpwBjDX71vNpPK86hvaufff7lRZUQudMOUcCne+r2vEgyNnoYcHp+PjLFlAHTVNzgczbktXLiQhQsXOh2GiEjC1dTUcNttt/HQQw/x0EMP8fjjj5/1/kSW4c2aNYs33njjXde3b9/OzJkzRzz+mWhlSSRFZKT5+NL9l/DZ76zn5S0NPPH8bu667jynw5J+LiiZyticUg61NbL50JvMq7jQ6ZASJrOyks66ejrr6siePMnpcM5q3bp1TocgIpJw1lo2btzIF7/4RT7xiU+wY8cO7rvvPtavX88ll1wy4DOJLMNbtGgRn/vc59i9ezdTpkwBYO/evbz00kt885vfHPH4Z6JkSSSFVJTm8Ol75/LPP32dn6x7g2kTCpg5udjpsCTCGMMNk6/kF1uf4Nk9L4+qZCmjohyAziRYWVqzZo3TIYiIJNyuXbtobm7uK6t78MEH2b17N4sWLWLDhg1MmDDhXc/k5uae0mZ8OB599FHg5H6pJ598ktLSUkpLS7nuuuv67lu8eDH/+Z//yZ133snXv/51jDH8/d//PePHj2fp0qUjiuFsVIYnkmKunF3BXdedRzBk+bf/qaWzO+B0SNLPdZMvxxhDdf0WWrpanQ4nYZKpfbjK8ERkNIomK9FkyePx8Itf/IKysjLuuOMOWlpa4jLvPffcwz333MMPfvADAD7xiU9wzz33sGzZslPuy87O5plnnmHatGl85CMf4cMf/jCTJ0/mmWeeIScnJy6xgVaWRFLS/XfMZMvOJnbXNfPT32znr9432+mQJKIos4C5Y2dR27CN5/e9xoLpI6/jTgbJ1j5cRGS0uffee7n33ntPuZaTk8OmTZviOu9Q9lhPmDCBxx5LbMdBrSyJpCCf18NnPjAXr8fwm5f2sHVXk9MhST83TLkSgGd3vzRqGnGcXFmqd/17XrVqFatWrXI6DBERcQElSyIpanJFPn968zQA/uOXm+jqUTmeW8wvv4i89BwOtDSw69g+p8NJCH9+Ht7sLIIdHfQ2NzsdzlktXbo0rvXvIiKSPJQsiaSwe26axqTyPBqOtvOzJ990OhyJ8Hl9XDvxMgCe3fOyw9EkhjGGzIrk2Le0ePFiFi9e7HQYIiLiAo4nS8aYMcaYnxtjdkQ+HjXGjBvEcznGmC8aY140xtQYY7YbY7YZY/7GDHDylTHmM5F7thhjao0xd8XlDYm4iN/n4dMfmIvHY1j7wm627znqdEgSES3Fe2l/NYFRcuZSsuxbUhmeiIhEOZosGWPSgKeBNGAWMBNoB541xpyrrcUk4OvAN6218621M4EvA98CvnLaPF8EvgostNbOBr4APGKMeU8M346IK00dV8DdN0zFWnjolxvp7h0d35i73fj8Cipzx9LR28nbTbudDichkiVZEhERiXJ6Zel+YDbwBWttwFobJJzITAE+fo5n24EfWGv7Tg+01q4BtgDvi14zxhQAfw98z1q7K3Lf08BThBMrkZT3wVunM74sl7rGdn7+u7ecDkci5pTPAmDToXefSJ6KkiVZqq+vp77e3TGKiDu4vWGNhI3kz8npZOluYL+1tu/HqtbaQ8D2yGtnZK3dY6396wFeygMa+319O5AFPHvafc8AM40xFwwncJFk4vd5+cwH5uIx8MT6nezYd8zpkASYG0mWNtZvcziSxEiWZKmyspLKSPc+EZEz8fv9dHZ2Oh2GDEJnZyfp6enDetbpZGk2sGeA63uAi4YykDEm2xjzdSAb+NJpc0THPH2O/q+fPt4SY0y1Maa6sbFxoFtEksq0CYXcdd1UQhb+/Zcb6VE5nuNmlE4l3ZvGvuY6jnWccDqcuMsoLwdj6D58mFDAvd0Zy8vLKS8vdzoMEXG5MWPGUFdXR0dHh1aYXMhaS29vL8eOHePgwYMUFxcPaxynD6UtAWoGuN4CZBljMq2150zZjTEvApcCu4D3WWtrT5sDoHWAOQAG/J2z1q4CVgFUVVXp/wGSEj50+wVseOMQBw638b9P7+Cj753pdEijmt/r58Ky6dTUb2XToTe4ccpVTocUV970dNJLS+g+0kj34SN9K01uoxI8ERmMvLw8IPzfjN7eXoejkYH4fD4yMjKYMGECGRkZwxsjxjHFyru62Z2NtfZqY4wf+CDwR2PM5621D8VyDpFUkO738ul75/L5/3yB/3tuF7dcOpHykmynwxrV5pbPoqZ+KxsbUj9ZAsisqKD7SCOddXWuTZZERAYrLy+vL2mS1OR0GV4TkDvA9VygYzCrSlHW2l5r7cPAT4FvGWPG9psjOubpcwCol7KMKjMmF3Fj1XgCwRA/Xjc6Ggu42ZzyCwHYcvjNUdFCPFn2LYmIiIDzydIWwi3ATzcZ2Hq2B40xacaYgVbGNgN+Tu552hL5fPo8k097XWTU+Oh7Z5Ce5uWVrQ1s3dl07gckbsZkF1OZO5bO3i7ebtrldDhxlwzJ0vz585k/f77TYYiIiAs4nSw9Dkw0xkyKXjDGlAEzgMf632iMKTPG9I/3y8DfDTBmdKzoitHvgA7g+tPuuwHYbq1VH2UZdYrzM7nnxvMBWP3EVoIhbctzUl9XvIbUX+nLjHSZ63TxvqDa2lpqa2vPfaOIiKQ8p5OlnxBeQXrAGOOLJEPfJNyp7vvRm4wxVwH1wHdPe/6TxpiZ/e67mvD5TM8DGwGstSeAf4zcOyVy383AbcDn4vKuRJLAXddPpbQwkz31LfzhtX1OhzOq9Z23NAqSpYyKcJc5N68sVVdXU11d7XQYIiLiAo42eLDW9hhjbgEeJHy2kgW2ATdaa9v63doGNAMN/a79FMgE/scYYwi/lyDwz8BDtl8PR2vtN40xXcA6Y0wgct891ton4/fuRNwt3e/lz++Yxb/8dzU/e/JNrr64kuxMv9NhjUozSqeS7ktnX3MdRzuOU5xV6HRIcZNeUoInLY3eEycItLfjy3ZfgxGV4ImISJTTK0tYaw9baz9krZ1mrZ1urb3bWnvgtHs2W2uLrLVf63dtj7X2C9baOdbai621s6y1s6213zgt0Yre/x1r7czIPXOttb9OwNsTcbWr51Qwc3IRzW09/OoPbzsdzqjl9/q5cMw0IPVXl4zHkxSrSyIiIuCCZElEnGOMYfGdF2EMrHlhF/VN7/o5gyTI3EhXvI2HUjtZgn5NHly6b2n58uUsX77c6TBERMQFlCyJjHJTxxdEWolbfrQm9b9Rd6tok4eth95K+RbimRXu7oi3YsUKVqxY4XQYIiLiAm49lFZEEuij753JS5vr2fDGITa/3cjF00qdDmnUKc0upjJvLHUth9jRtItZkbK8VHSyfXidw5EMbNmyZU6HICIiLqGVJRGhKC+De24Kf3P+wzXbCAZDDkc0Os0dOzpaiPe1D3fpypLK8EREJErJkogAcNd15zGmKIu9DS08tUGtxJ0wWlqIR8vwuuobsCEl5iIi4l5KlkQEgDS/l79YEP5m/RdP7aC7N7X3zbhRtIX4/kgL8VTly8nGl5dHqKeHnuMnnA7nXWpqaqipqXE6DBERcQElSyLS58rZ5Zw3Lp/jrd387pW9Tocz6vi9fi4aMx1I/VK8tIJ8AAKtLQ5H8m5VVVVUVVU5HYaIiLiAkiUR6WOM4UO3XgDAY8+8o9UlB4yWUjxfXh4Avc3uS5bmzZvHvHnznA5DRERcQMmSiJzikpllTI2sLj358l6nwxl1+lqIH36LQDDgcDTx43dxsqQyPBERiVKyJCKnMMbwwdsiq0vPvkNXT+p+w+5G0RbinYEudhzd7XQ4cePPjyRLLe5LlkRERKKULInIu1wyo4yp4ws4ob1LjphbfiEAGxu2ORxJ/ERXlgJKlkRExMWULInIu4T3LoUbDTz2zE66urW6lEjRUrxUbvLg5pWliooKKiLtzUVEZHRTsiQiA6qaUca0CQWcaOvmt9q7lFAXlJyH3+vnQHM9Ld1tTocTF75c9+5ZamhooKGhwekwRETEBZQsiciAjDF8MNIZ7/Hn3tHqUgL5vX7OK5wAwDtH9zgcTXy4eWWprq6Ouro6p8MQEREXULIkImc0/4IxTJ9QSHNbD799OTW/aXeraSXnAfB2U2o2eYgmS27cs6QyPBERiVKyJCJnFO6MF9m79OxOOrW6lDDTiicD8HaKdsRzcxmeiIhIlJIlETmredPDq0st7T389iWtLiXKtJIpAOw8to9gKPUOB/bn5QLQ29qKDYUcjuZUS5YsYcmSJU6HISIiLqBkSUTOyhjDh26L7l3S6lKiFGTkMSa7mO5ANwea650OJ+Y8fj/e7CwIhQi0tzsdzilWr17N6tWrnQ5DRERcQMmSiJzT3OmlTJ8YXl1a92JqloW50bTi8OpSqpbi+V1airdy5UpWrlzpdBgiIuICSpZE5Jz6ry6teX43Pb2pVxbmRtFSvLebUrP80a1NHlSGJyIiUUqWRGRQ5k4rZXJFHifaullfe9DpcEaFlG/ykOfe9uEiIiKgZElEBskYw13XTQXg18/vwlrrcESpb0LBONK8fg61NdLS1ep0ODHnd2mytHbtWtauXet0GCIi4gJKlkRk0K6ZU0lRXjr7D7Wy8e1Gp8NJeT6Pl/OKJgHwdgoeTtt3MK3L9iwtWrSIRYsWOR2GiIi4gJIlERk0v8/DgqvD+2ieWL/L4WhGh1QuxYuuLLltz9KCBQtYsGCB02GIiIgLKFkSkSG5/YpJpKd5qd1xhH0N7vomNxVN72vykILJUr7K8ERExN2ULInIkORmpXFT1XgAnnheq0vxdn5kZWlXCh5O29fgwWVleCIiIlFKlkRkyO689jyMgWdrDnK8tcvpcFJafkYeZTmldAd72HeizulwYupkg4fUa14hIiKpQcmSiAxZRWkOl84cSyAY4rcv7XU6nJSXqvuWTp6z1OxwJKcyxmCMcToMERFxASVLIjIsd113HgC/fXkP3TqkNq6mFUf2LaVYRzxfrsrwRETE3ZQsiciwzJpSzNRx+bS09/BczQGnw0lp0yJNHt5JsSYP3swMjN9PqKeHYJd7yjmttTpHTEREACVLIjJMxhjujBxS+8TzuwiF9M1lvEzIryDdl87h9iaau1JnFcYY49qDaUVEREDJkoiMwNUXV1Ccn8GBw23U7jjidDgpy+vxMrVoIpB6pXh+dcQTEREXU7IkIsPm83pYGDmk9tfrdzocTWqLthBPtfOWTjZ5cE+ytHDhQhYuXOh0GCIi4gJKlkRkRG67YhIZaV42v9PEnnp3dTVLJSnb5MGFZXjr1q1j3bp1TochIiIuoGRJREYkJ9PPLZeFS8R+vV6H1MZLtMnDrmN7CaTQ4bTRlSU3JUtr1qxhzZo1TochIiIuoGRJREYsWor3wqY6Wjt6HI4mNeWl51CeM4aeYC/7Thx0OpyYceOeJZXhiYhIlJIlERmx8pJs5k4rpTcQ4plqtRGPl/NLUm/fkhuTJRERkSglSyISE7dfMQmA372yV2fUxMnJfUsplCxFGzy0uidZWrVqFatWrXI6DBERcQElSyISE5fOGkthbjoHj7Txxu6jToeTklKxyYMvLxdw18rS0qVLWbp0qdNhiIiICyhZEpGY8Hk93HzpBAB+98o+h6NJTRPyK8jwpdPYfpQTnanRedCflw+4q8HD4sWLWbx4sdNhiIiICyhZEpGYue3ySRgDL22pp6VdjR5izePxMLVoEpA6q0snz1lqdTiSk1SGJyIiUUqWRCRmyoqymDt9DIFgiGeq9zsdTkqaFm3ykCL7lnzZ2WAMgbY2QoGA0+GIiIicQsmSiMTU7ZdPAsKleGr0EHt9+5ZSpCOe8Xrx5Yb3LQVa3bG6VF9fT319vdNhiIiICyhZEpGYunRmGUV5GdQ1trFtlxo9xNr5xeGVpV3H9xNMkcNp+9qHu6QUr7KyksrKSqfDEBERF1CyJCIx5fV6uOWyaKOHvc4Gk4Jy03MozS6mN9hLQ+sRp8OJiei+pd5mdzStKC8vp7y83OkwRETEBZQsiUjM3XrZRIyBl7c20NzW7XQ4KWdiwTgA9p446HAkseGPtA8PuKQjnsrwREQkSsmSiMTcmMIs5l9QRiAY4o+vH3A6nJQzqSBcIrYvRZIlX7R9uIvOWhIREQElSyISJ7dfPhGA37+6V40eYiy6spQqyVJfGZ5LGjyIiIhEKVkSkbiomlFGcX4G9U3tbNnZ5HQ4KSVVy/Dcsmdp/vz5zJ8/3+kwRETEBZQsiUhceL0ebr0svLqkRg+xNSa7mAxfOie6WmjuSv7SNb/LyvBqa2upra11OgwREXEBJUsiEje3XjYRj4FXtzVwolWNHmLFYzxMzI/uW6pzOJqRi5bhueWcperqaqqrq50OQ0REXMDxZMkYM8YY83NjzI7Ix6PGmHGDeK7cGLPCGPOaMWajMeYtY8zjxpiLBrj3OWPMdmPMptM+PhqfdyUiACUFmVTNGEsgaPnj6/udDielpFIpnk9leCIi4lKOJkvGmDTgaSANmAXMBNqBZ40xOed4fBnwQeBua+1cYA4QBDYMlDAB77XWzjnt4+FYvRcRGdhtV4RL8Z7asE+NHmIolZo89JXhuaR1uIiISJTTK0v3A7OBL1hrA9baIPAFYArw8UE8/y/W2gMA1tou4ItAJrAkTvGKyBDNnz6Gwtx06pva2bH/uNPhpIxJhdFkKYXK8FpaXZFQL1++nOXLlzsdhoiIuIDTydLdwH5r7e7oBWvtIWB75LWz+WvgR6ddi54iWBizCEVkRLxeD9fPHw/AM9U6cylWxudXYDDUtTTQG+x1OpwR8fj9eDMzscEgwfYOp8NhxYoVrFixwukwRETEBZxOlmYDewa4vgcYqJSuT2QlKnTa5WmRz88N8MjfRPY3vWWMed4Y8+dDjlZEhuXGqnCy9MLGOnoDQYejSQ0ZvnTG5pYStCHqWg45Hc6I9e1banF+39KyZctYtmyZ02GIiIgLOJ0slQADtT9qAbKMMZlDHG8J8Abws9OunwB2AjcQ3hv1EPB9Y8y3zjSQMWaJMabaGFPd2Ng4xDBEpL9J5XlMqcinrbOX17cfdjqclJFKTR7c1D5cZXgiIhLldLJ0JmbIDxhzI3Av8KfW2lN6FFtr77LWPmStbbfWBq21jwL/BXzWGDNhoPGstaustVXW2qrS0tLhvAcR6eeGKpXixdqkgtTbt9Tb4o724SIiIuB8stQE5A5wPRfosNZ2DmYQY8zFwMPAImvt9kHOvYHw+79kkPeLyAhcN7cSj8dQ/eZhmtt05lIspFZHvPA/BQEXlOHV1NRQU1PjdBgiIuICTidLW4BJA1yfDGwdzADGmNnAr4EPWGtfHuD1NGNM/gCPRjdOeAcVqYiMSGFeBvOmjyEYsrywKflXQtxgYkH4YNq9Jw66oovcSPjyIitLLijDq6qqoqqqyukwRETEBZxOlh4HJhpjJkUvGGPKgBnAY/1vNMaUGWM8p12bDTwBfMRa+2LkWrkxZmW/264EfjXA3NETBzeO9E2IyODcqK54MVWcWUh2WhZtPe0c6zzhdDgj4s+P7Flqdb4Mb968ecybN8/pMERExAWcTpZ+QngF6QFjjC+SDH2TcDe870dvMsZcRbgt+Hf7XbsI+CPwO2CSMeY+Y8x9hPctTT9tnpuMMXf0e/Z6YCnwM2vtO7F/WyIykEsvHEtWho93DpzgwGHnvylOdsaYfvuWkrsUL1qG54aVJZXhiYhIlKPJkrW2B7iFcEncduBNIA+40Vrb1u/WNqAZaOh3bQXhbnp/Rbj7XfTjwdOmqQU+D3zZGLPZGLMT+B7wdeAvYv2eROTM0v1err44XDqm1aXYSJWOeP686MG0zu9ZEhERifI5HYC19jDwoXPcsxkoOu3a+wY5fgvw7ciHiDjsxqrxPLVhH8/VHOC+98zA6xly80vpJ1U64vWV4akbnoiIuIjTZXgiMsrMmFREWVEWTc1dbNvZ5HQ4SS9VOuK5qcFDRUUFFRUVTochIiIuoGRJRBLK4zHcGD1zqUaleCM1Lm8sXuOhoe0I3YEep8MZtr49Sy3OJ0sNDQ00NDSc+0YREUl5SpZEJOFuiHTFe3lLPZ3dAYejSW5+r5+KvLFYaznQXO90OMPmzcrC+HyEuroIdjt7DlddXR11dcld1igiIrGhZElEEq68JJsZk4ro6gnyylb9BH+kUqHJgzGmX5MHZ/ctqQxPRESilCyJiCOipXjPqiveiE2KHE6b/PuW3FOKJyIiAkqWRMQhV19cgd/nYfPORppOdDodTlJLlSYP/r4mD862D1+yZAlLlixxNAYREXEHJUsi4oicrDQunTUWa+G52uT+Jt9pE/u1Dw/ZkMPRDJ8/P5IsOVyGt3r1alavXu1oDCIi4g5KlkTEMX2leOqKNyIFGXkUZOTRGeiisf2o0+EM28k9S86W4a1cuZKVK1c6GoOIiLiD44fSisjoNW/6GHKz/Ow/1Mq+hhYmluc5HVLSmlgwjhOHtrPvRB1lOaVOhzMsPheV4YmIiIBWlkTEQT6vhytnh7uOvbBJrZpHIhU64vWV4bU6W4YnIiISpWRJRBx1zZxwJ7fnN9VhrXU4muSVCh3x/Hn5APQ2O1uGt3btWtauXetoDCIi4g4qwxMRR114XgkFuek0NLWzq66ZqeMKnA4pKaVCRzx/pHW403uWFi1aBKDkXUREtLIkIs7yegxXR0vxNqoUb7gqcsvwe3wcaT9KR29ytmLvK8NzeM/SggULWLBggaMxiIiIOyhZEhHHXR0pxXths0rxhsvr8TIuvxyA/SeSM+n0RcvwHG4drjI8ERGJUrIkIo6bMamIkvwMGo93smPfcafDSVrJ3uTBn5sDQKCtDRsMOhyNiIiIkiURcQGPx5xcXVJXvGGb1O9w2mRkvF58uTlgLb2tbU6HIyIiomRJRNwh2hXvxc11BEMqxRuOSSnR5MH5g2mNMRhjHJtfRETcQ8mSiLjC+eMLGFucxbGWbrbvPup0OElpQqR9+P7mOkKhkMPRDE/fwbQOd8QTEREBJUsi4hLGmL7VJZXiDU9OWjZFmQX0BHtp7EjOhDO6suTkWUvWWjUaERERQMmSiLhINFl6aUs9wWByrow4rTJvLAAHWw45HMnw9LUP18qSiIi4gJIlEXGNSeV5jBuTQ0t7D5t3NjkdTlIalxduH17X0uBwJMPjhj1LIiIiUUqWRMQ1TinF0wG1w9K3stScnCtLPheU4S1cuJCFCxc6Nr+IiLiHkiURcZVosvTKtgZ6AzprZ6iiK0sHk3Vlqa8Mr9mxGNatW8e6descm19ERNzD53QAIiL9jS/LZVJ5HnsbWti4o5FLZ411OqSkMi6yslTXcghrbdK1wO5r8NDS6lgMa9ascWxuERFxF60siYjrXDs3vLr0vErxhiwvI5fc9Bw6A10c6zzhdDhD5oY9SyrDExGRKCVLIuI60VK817Y30NUTcDia5HOyyUPy7VvqK8NzcM+SiIhIlJIlEXGdscXZnD++gM7uIDVvHnE6nKRzsn148u1b6n8orVNnHa1atYpVq1Y5MreIiLiLkiURcaXo6tLzmw46HEnyGZfEZy1509PxpKdjAwGCnZ2OxLB06VKWLl3qyNwiIuIuavAgIq509cWV/GjtG1S/eYSu7gAZ6frP1WAl/1lLuXQ3dhNobcWXlZXw+RcvXpzwOUVExJ303YeIuFJpYSbTJxSyY/9xanYc4arZFU6HlDT62oc3J2ey5M3OhsYmAu3tjsyvEjwREYlSGZ6IuNaVs8Pf9L+8pd7hSJJLYWY+mb4MWnvaaelyrgX3cPmyswEItDmTLImIiEQpWRIR17riovBq0uvbD+uA2iEwxvTbt5R8q0vRZCnY3uHI/PX19dTXK0EXERElSyLiYuUl2UyuyKOzO8DGtxudDiepVOZHSvGSsMmDLyeysuRQGV5lZSWVlZWOzC0iIu6iZElEXO3KyF4lleINTTKvLHmznU2WysvLKS8vd2RuERFxFyVLIuJqV14U/qZ1w7ZDBIIhh6NJHsncEe9kGZ4zyZLK8EREJErJkoi42oSxeYwbk0NbZy9bdzY5HU7SqEzis5Z8Dq8siYiIRClZEhHX6yvF25p8qyROKc0qJs3r53hnMx09zhzuOlze7PDZSgGHGjyIiIhEKVkSEdeLluK9urWBYMg6HE1y8Hg8VOSWAcm3b+lk6/A2R+afP38+8+fPd2RuERFxFyVLIuJ6UyrzKSvK4kRbN2/uOep0OEmj73DaJCvFc3rPUm1tLbW1tY7MLSIi7qJkSURczxijUrxhGNfXPjy5fs98OTmAc3uWqqurqa6udmRuERFxFyVLIpIUrpwd/sb/lS31hFSKNyjRJg91Sbay5PSeJZXhiYhIlJIlEUkK08YXUpyfQVNzF+8cOO50OEnhZBlekq0sOVyGJyIiEqVkSUSSgsdjuCLS6OHlLcn1zb9TynJK8RoPTe3H6Ap0Ox3OoHkzM8EYgp2d2GAw4fMvX76c5cuXJ3xeERFxHyVLIpI0Tu5bqsdaleKdi8/jpTy3DIulvuWw0+EMmvF48GZFSvE6El+Kt2LFClasWJHweUVExH18TgcgIjJYMycXU5CTzqGjHeypb2FKZb7TIbleZd5YDrY0cLClgSlFE5wOZ9B82dkE29sJtrfjz81N6NzLli1L6HwiIuJeSpZEJGl4PYbLLhzL71/dx8tb6pUsDcK4vHI2sDHpmjz4srPpBgJtid+3pBI8ERGJUhmeiCSV/qV4cm7RjnjJ1uThZEc8NXkQERHnKFkSkaQye2oJOZl+Dhxu48DhVqfDcb1oR7ykW1ly8KylmpoaampqEj6viIi4j5IlEUkqPq+HS2eFV0te3qLVpXOpyB2DMYZDbY0EggGnwxk0J9uHV1VVUVVVlfB5RUTEfZQsiUjSuSpaiqcW4ueU5ktjTHYJIRuioe2I0+EMms/Bg2nnzZvHvHnzEj6viIi4j5IlEUk6c6aVkpHmZXd9M0eOJf6b6WQzLgn3LXkjK0sqwxMRESc5niwZY8YYY35ujNkR+XjUGDNuEM+VG2NWGGNeM8ZsNMa8ZYx53Bhz0Rnu/4wxZrsxZosxptYYc1fM34yIJESa38v8C8oAePWN5EkAnFIZ2bd0sDl5fq+cLMMTERGJcjRZMsakAU8DacAsYCbQDjxrjMk5x+PLgA8Cd1tr5wJzgCCw4fSEyRjzReCrwEJr7WzgC8Ajxpj3xPDtiEgCXX5ROAF4dWtyNS5wQnRlKZmaPPgcXFkSERGJcnpl6X5gNvAFa23AWhsknMhMAT4+iOf/xVp7AMBa2wV8EcgElkRvMMYUAH8PfM9auyty79PAU8C3YvdWRCSRqmaU4fUY3tjdREt7j9PhuFq0I97BJEqW+srwHDhnqaKigoqKioTPKyIi7uN0snQ3sN9auzt6wVp7CNgeee1s/hr40WnXoq2xCvtdux3IAp497d5ngJnGmAuGGrSIOC8n089FU0sIWXh9e/IkAU6InrXU0HqYYCjocDSD48txbmWpoaGBhobkKVkUEZH4cTpZmg3sGeD6HmDAvUdRkZWo0GmXp0U+P3faHNExT5+j/+sikmQuvzBSirdN39ieTaY/g+KsQnpDAY60H3U6nEFxcs9SXV0ddXV1CZ9XRETcx+lkqQQY6FTJFiDLGJM5xPGWAG8APzttDgaYpyXyuXiggYwxS4wx1caY6sbGxiGGISKJcFnkvKXaHY109STPGUJOSLaOeE7uWVIZnoiIRDmdLJ2JGfIDxtwI3Av8qbW2e6RzWGtXWWurrLVVpaWlQw1HRBKgpCCT88cX0NMbZNPb+qHG2UQ74iVLk4eTrcPVGl5ERJzjdLLUBOQOcD0X6LDWdg5mEGPMxcDDwCJr7fYB5oiOefocAMlRkyIiA7oi0hXvla3JsWLilL6VpSRpH+7NzACPh1BXF6FAYlcNlyxZwpIlS859o4iIpDynk6UtwKQBrk8Gtg5mAGPMbODXwAestS+fYQ4GmGfyaa+LSBKK7lt6ffshgsHTtzFK1LgkW1kyxuDLzgIg2JHY1aXVq1ezevXqhM4pIiLu5HSy9Dgw0RgzKXrBGFMGzAAe63+jMabMGOM57dps4AngI9baFyPXyo0xK/vd9jugA7j+tLlvALZba9+KzVsRESeMG5NDZWk2rR29bN9zzOlwXCvaEe9g6yFC7+qN405O7VtauXIlK1euPPeNIiKS8pxOln5CeAXpAWOML5IMfZNwp7rvR28yxlxFuC34d/tduwj4I+FkaJIx5j5jzH2E9y1Nj95nrT0B/CPwSWPMlMizNwO3AZ+L55sTkfgzxqgr3iDkpueQl55Dd6CbY50nnA5nUJw6a0lleCIiEuVosmSt7QFuAYKEz1Z6E8gDbrTWtvW7tQ1oBvp/J7SCcKe7vyLc/S768eAA83wT+CdgnTFmC/CvwD3W2idj/Z5EJPH6J0vWWoejca/o6lJ9y2GHIxkcJ9uHi4iIAPicDsBaexj40Dnu2QwUnXbtfUOc5zvAd4YYnogkgWkTCinMTefI8U721LcwpTLf6ZBcqTy3jDcbd1LfepjZY2c4Hc45OVWGt3btWgAWLlyY0HlFRMR9nC7DExEZMY/HcJlK8c6pMje8spQsTR68DiVLixYtYtGiRQmdU0RE3EnJkoikhMsvDCcCaiF+ZpV5ZQDUtyZHsuTLcWbP0oIFC1iwYEFC5xQREXdyvAxPRCQWZk8tITPdx96GFg4dbWdscbbTIblORW4kWWo54nAkg+PUnqVoGZ6IiIhWlkQkJfh9XqpmhJOBV7clx8pJopVmF+Pz+DjaeZyu3i6nwzmn6DlLgfbEnrMkIiISpWRJRFJGtBRP+5YG5vV4Kc8pBaC+1f2rS07tWRIREYlSsiQiKWP+BWX4vIY39xylua3b6XBcqTyJ9i2dLMNrO8edsWWMwRiT0DlFRMSdlCyJSMrIzvQz+/xSQhZe3+7+ZMAJJzviuf+sJV9ODgCBNpXhiYiIM5QsiUhKiR5Q+8pWJUsD6Wvy0JoEyVLfnqXEluFZa3W4sYiIAEqWRCTFXDYrvHKy6e0jdHUHHI7GfSrzwr8/9Ulw1pL2LImIiNOULIlISinKy2D6hEJ6AiE2vu3+JgaJ1rey1HaEkA05HM3ZOdU6XEREJErJkoiknMv6uuK5f/Uk0bLSMinIyKM32EtTx3GnwzkrT3o6xusl1NNDqLc3YfMuXLiQhQsXJmw+ERFxLyVLIpJyoqV4r28/TDDo7tUTJyRLKZ4xpl8pXuKaPKxbt45169YlbD4REXEvJUsiknLGl+VSXpJNa0cPb+495nQ4rpOMTR4SWYq3Zs0a1qxZk7D5RETEvZQsiUjKMcb0rS5teMPdqydOiCZLdS5fWYJ+7cMTmCypDE9ERKKULIlISoq2EN+w7ZDaQJ+mrwwvKVaWImV4bYk9mFZERASULIlIirpgUhF52Wk0HG1n/+FWp8NxlYq+PUvuT5a8fWctJW7P0qpVq1i1alXC5hMREfdSsiQiKcnrMVwyM1xutkFd8U5RklWI3+vneFczHT2dTodzVk60D1+6dClLly5N2HwiIuJeSpZEJGVdNitSivdGg8ORuIvHeKjIGQO4vxTP58DBtIsXL2bx4sUJm09ERNzL53QAIiLxMndaKWk+D2/vP8HR5k6K8zOdDsk1yvPK2NdcR33rYaYWT3I6nDPyOpAsqQRPRESihp0sGWNmn+11a+2W4Y4tIhILGek+5kwbw2vbD/Ha9sO854pJTofkGpW54X1Lbu+I50QZnoiISNRIVpaeOMtrFpgygrFFRGLisgvH8tr2Q2zY1qBkqZ/KvOQ4a8mJMrz6+noAKioqEjaniIi407CTJWvt5FgGIiISD5fMLMMY2PxOEx1dvWRl+J0OyRX6DqZ1+8pSTjRZSlw3vMrKSgC1nBcRkaE3eDDGfNQY83/GmKeMMd8zxtxijPHGIzgRkZEqzM3ggolFBIIhNu5odDoc14gmSw1tjYRCIYejOTOvA+cslZeXU15enrD5RETEvYaULBlj/gH4CXAtUArcDfwe2GqMmRnz6EREYuCyWeH9Oa+qK16fDH8GxZmFBEIBjnQcdTqcM3Jiz1J9fX1fKZ6IiIxu50yWIitJ0yJffgL4IVBqrZ1rrS0DLgX2Aq8aYy6IW6QiIsN02YXhZKl6+2ECQfeuoiRaRV6kfbiLS/FO7llKXBmeiIhI1GBWln4MvGmMOU54NSkTuNsYMxXAWlttrX0v8DvggbhFKiIyTOPG5FJZmkNbZy/b97h3FSXRKvo64rm3yYM3OwtIbIMHERGRqMEkS0XArcC/EO5y917gl8AOY0yzMeZ5Y8x3gCPADfEKVERkJC6PrC5t2ObeVZRE62vy4OKOeJ60NIzPh+3tJdTTk5A558+fz/z58xMyl4iIuNs5kyVrbbO19o/W2n8mXG73RcIrTO8FvgE0AHcAHweyjTGtxpgXjDHfjl/YIiJDc/mF4Q37r75xSF3OIirzwglkfat7E0hjTMLbh9fW1lJbW5uQuURExN2G2g1vNfBPwHnW2t9bax+w1t5rrT0f+AgQAP4ROAQsim2oIiLDN21CIQW56Rw51sHehhanw3GFiuhZSy4uw4P+7cMTkyxVV1dTXV2dkLlERMTdhnrO0r8Cs4GXjDG/IbxPqQGYDHwe2GCt/ZfYhigiMnIej+HSmWN5asM+Xt12iMkV+U6H5LiizALSvWk0d7fS1t1OTnq20yENyNvXES8xTR5UgiciIlFDWlmy1oastR8GlgAVwPeA/wO+DbQBi2MeoYhIjES74m1QC3EAPMaTFPuWfA6ctSQiIgLDOJQWwFr7Y2vtpUAlcA0wD7jAWrsjlsGJiMTSxeeXkp7mZdfBZhqPdzodjiv0leIlQ7KUoDK85cuXs3z58oTMJSIi7jasZCnKWttgrX3JWrvJWqvDS0TE1dL9XuZND58t9JpWl4CTHfHqXHzWkjfBydKKFStYsWJFQuYSERF3G+qeJRGRpHb5hWN5ZWsDr247xB1XT3E6HMdFO+LVuXplKXzWUqL2LC1btiwh84iIiPspWRKRUaVqxlg8HsPWXU20dfaSk+l3OiRHRVeWGlzcEc+JMjwREREYYRmeiEiyyctOY9bkYoIhS/Wb7k0QEqU8kiwdajtCIBR0OJqBJboMT0REJErJkoiMOpdHu+Jt076ldF8apVlFBG2II+1NToczIF9ODgCBtsQkSzU1NdTU1CRkLhERcTclSyIy6lx2YTkANW8dpjfgztWURDp5OK07mzyc3LOUmGSpqqqKqqqqhMwlIiLupmRJREadsqIsJlfk0dkdZMtOd66mJFJFbqTJg0v3LSV6z9K8efOYN29eQuYSERF3U7IkIqPSZbPCq0uvbnPnakoiVbr8rKVE71lSGZ6IiEQpWRKRUSm6b+m1NxoIhazD0Tgr2hHPtWV4OeFkKVFleCIiIlFKlkRkVJpSmU9pYSbHWrp558Bxp8NxVIXLz1ryZYX3LAXaO7B2dCe2IiKSWEqWRGRUMsZw2axIV7w33LmikiiFGflk+jJo62mnpavV6XDexZOWhictDRsIEOrpift8FRUVVFRUxH0eERFxPyVLIjJqXX5hdN/S6G4hbozp64h30KWleN5IR7xEtA9vaGigoWF0/50QEZEwJUsiMmrNmlJMdqafA4fbqGtsczocR1VGS/Fcmiz5ssNnLSVi31JdXR11dXVxn0dERNxPyZKIjFo+r4dLZoZXVEb7AbXj8sKrbHUt7vx9SGT7cJXhiYhIlJIlERnVLlcLcaDfylKrO38ffDnRJg/qiCciIomjZElERrW500vx+zy8te8Yx1u7nA7HMdFkyb17liIrSwnYs7RkyRKWLFkS93lERMT9lCyJyKiWleHn4vNLsRZe3+7O1tmJUJZdgs/j42jHcbp63Zc0RsvwErFnafXq1axevTru84iIiPspWRKRUS96QO1o7orn9XgpzykFoN6F5y0lcs/SypUrWblyZdznERER9/M5HYCIiNMunTkWYzaz6e1GOrsDZKaPzv80VuaVc6ClgYMth5hSNNHpcE7hTWCypBI8ERGJSvqVJWNMujHmW8aYkDHmeqfjEZHkU5iXwfQJhfQGQmzcccTpcBzj5vbhvpxoGV6Hw5GIiMho4niyZIwZY4z5uTFmR+TjUWPMuEE+eyGwAbgZMGe57zljzHZjzKbTPj4ao7chIknussgBtRvecF+ikCiVkYNpXZksJXBlae3ataxduzbu84iIiPs5WmtijEkDngbeBmYBFvgR8KwxZq619lynRH4Z+AJQDvz4HPe+11q7d2QRi0iquvzCsfz0N9t57Y1DBIIhfF7Hf5aUcJV9Zy2N7mRp0aJFAFhr4z6XiIi4m9PfDdwPzAa+YK0NWGuDhJOfKcDHB/H8R621v49ngCIyOowbk8v4shzaOnt5Y9dRp8NxREVuGQbDobYjBEJBp8M5RSJbhy9YsIAFCxbEfR4REXE/p5Olu4H91trd0QvW2kPA9shrZ2WtDcQxNhEZZa64qAKAl7fWOxyJM9J9aZRkFxG0IQ61uWvvli87fChtIlqHqwxPRESinE6WZgN7Bri+B7goxnP9jTHmNWPMW8aY540xf362m40xS4wx1caY6sbGxhiHIiJudEVk39Kr2w4RCo3OEqxxLm3ykMgyPBERkSink6USoHWA6y1AljEmM0bznAB2AjcQ3hv1EPB9Y8y3zvSAtXaVtbbKWltVWloaozBExM3OG5dPaWEmx1q6ePvAcafDcURlrjuTJW9WeGUp0N6uvUQiIpIwTidLZ3LGznbDYa29y1r7kLW23VobtNY+CvwX8FljzIRYziUiycsYwxUXhVeXXtkyOg+ojbYPP+iyZMnj9+NJT4dQiFBXV1znMsZgTEz/GRIRkSTldLLUBOQOcD0X6LDWdsZx7g2E3/8lcZxDRJJMtBTvlW0No3IFI9oRr95lyRKcPGspoLOWREQkQZxOlrYAkwa4PhnYGosJjDFpxpj8AV6KtnryxmIeEUkNMyYXk5+TRkNTO/sODVQlnNr69iy1HiZkQw5Hc6pE7Vuy1o7KRFlERN7N6WTpcWCiMWZS9IIxpgyYATzW/0ZjTJkxZjjxXgn8aoDr8yOfNw5jTBFJUV6P4fLo6tKW0dcVLyc9m/z0XLoD3RzrOOF0OKeItg9PREc8ERERcD5Z+gnhFaQHjDG+SDL0TcLd8L4fvckYcxVQD3x3mPPcZIy5o9941wNLgZ9Za98Z5pgikqIu71eKNxq5dd+SOuKJiEiiOZosWWt7gFsIl8RtB94E8oAbrbVt/W5tA5qBU75zMcbca4zZBHwtcumHxphNxpgr+91WC3we+LIxZrMxZifwPeDrwF/E/l2JSLK7+PwSsjJ87KlvoaFp9H1jXtnXPtxdyWJfstTWdo47R2bhwoUsXLgwrnOIiEhy8DkdgLX2MPChc9yzGSga4PovgV+e49kW4NuRDxGRc/L7vFwyYyzrNx7kla0NvO+GqU6HlFCVLj1ryZsdbR8e3wYP69ati+v4IiKSPJwuwxMRcaVoC/FXR2EpXl+y1OquZMmXoD1La9asYc2aNXGdQ0REkoPjK0siIm4074IxpPk8vLn3GMdauijKy3A6pIQZ7XuWVIInIiJRWlkSERlAZrqPudPHAKNvdak4s5AMXzqt3W20dMd3f9BQnDxnafTtIxMREWcoWRIROYNoKd4rW0dXsmSMoTLXfU0eEtU6fNWqVaxatSquc4iISHJQsiQicgaXzhqLx2PYurOJ1o4ep8NJKDc2eThZhhffBg9Lly5l6dKlcZ1DRESSg/YsiYicQW5WGhedV8zmd5p4ffshbqya4HRICXMyWTrscCQnJWrP0uLFi+M6voiIJA8lSyIiZ3HFRRVsfqeJV7Y2jNJkyX1lePE+Z0kleCIiEqUyPBGRs7j8wnDSUPvWEbq6Aw5HkzjjXFyGF4xzGZ6IiEiUkiURkbMozs9k+sRCegIhancccTqchCnLKcXr8dLYcYyuQLfT4QDgix5K29GBDYXiNk99fT319fVxG19ERJKHkiURkXO4chR2xfN6vJTnhFun17tk35LxevFmZkIoRLCrK27zVFZWUllZGbfxRUQkeShZEhE5h8sjydJr2w/RGwg6HE3iVOSVAe4qxUvEvqXy8nLKy8vjNr6IiCQPJUsiIudQUZLDlIp8OroCbNzR6HQ4CdO3b6nVPStq/rxcAAKt8UuWVIYnIiJRSpZERAbhqosrAHhxc53DkSROZW54deWgi1aWfLnhZKm3pcXhSEREZDRQsiQiMghXR5KlDW+MnlI8Nx5M27ey1NLqcCQiIjIaKFkSERmEitLRV4oX3bN0qPUIgZA7EkR/Xh4Q35Wl+fPnM3/+/LiNLyIiyUPJkojIIF09J7y69MIoKcXL8KVTmlVE0IY43OaOBNGXgGSptraW2trauI0vIiLJQ8mSiMggRfctvfbGIXp63bHSEm9uK8Xz50YbPMSvDK+6uprq6uq4jS8iIslDyZKIyCBVlOQwpTJaijc6DqitzAs3eXBLspSIlSWV4YmISJSSJRGRIYg2enhxy+hoLe26lSU1eBARkQRSsiQiMgTRUrwN20ZHKd441yVL8V9ZWr58OcuXL4/b+CIikjyULImIDEG0FK+ze3SU4vWtLLUewlrrcDTgy4uesxS/laUVK1awYsWKuI0vIiLJw+d0ACIiyebqiyvYXdfMi5vruezCcqfDiavc9Bzy03Np7m6lqeMYpdnFjsbTv8GDtRZjTMznWLZsWczHFBGR5KSVJRGRIbr64kogfEDtaCjFm1AQLj3cd8L5lumetDQ8GRnYQIBgZ2dc5lAZnoiIRClZEhEZovKSbM4bFy7Fqx0FpXgT8scBsL/Z+WQJTu5bCsRx35KIiAgoWRIRGZarZodXW17anPpd8SYWhFfS9rtgZQlOdsSL176lmpoaampq4jK2iIgkFyVLIiLDMJpK8SbkR8rwXLKyFO+zlqqqqqiqqorL2CIiklyULImIDMNoKsUbl1eOMYaG1iP0BHudDifuZy3NmzePefPmxWVsERFJLkqWRESGKbq69OKm1C7FS/OlUZFTRsiGXHHeki83srLUGp+VJZXhiYhIlJIlEZFhujpyQO1r2xvoTvFSvPF9HfEOOhxJ/FeWREREopQsiYgM09jibKaOy6ezO0jtW6ldijcx3z1NHvxx3rMkIiISpWRJRGQEroqU4qV6V7y+jnjNzr9PX5y74VVUVFBRURGXsUVEJLkoWRIRGYFoKd6GNxro6gk4HE38TIisLLmhI168z1lqaGigoaEhLmOLiEhyUbIkIjICY4uzmT6hkK6eIK+94Xzzg3gpyS4i05dBc1cLzV3Olr/F+5yluro66uqcTwpFRMR5SpZEREbounnjAHiu1vnmB/HiMZ6+85acLsWLnrMUiFM3PJXhiYhIlJIlEZERumZOJR6PofatIzS3dTsdTtyMj+xb2udwkwdfTg4Ava1t2FDI0VhERCS1KVkSERmhgtx05kwrJRiyvJjCjR7c0hHP4/Phzc6GUIhAe3vMx1+yZAlLliyJ+bgiIpJ8lCyJiMTADZFSvPUpXIoX7Yi3r9n59xjPs5ZWr17N6tWrYz6uiIgkH5/TAYiIpILLLiwnPc3Lm3uPcehoO2OLs50OKeaiHfEOthwiGAri9Xgdi8Wfl0dXwyF6W1rIrIzt/qKVK1fGdDwREUleWlkSEYmBzHQfl88qB1J3dSkrLZOSrCJ6g70camt0NBZfbvw64qkMT0REopQsiYjEyPXzT3bFs9Y6HE18THBJk4e+MrzW+LQPFxERASVLIiIxM2daKfk5aRw80sauumanw4mLviYPDh9OG20f3huHg2nXrl3L2rVrYz6uiIgkHyVLIiIx4vN6uObicDKRqqV4EwoiZy05vbKUG7+VpUWLFrFo0aKYjysiIslHyZKISAxdFynFe37jQYKh1CvFm5gffn/7UnhlacGCBSxYsCDm44qISPJRNzwRkRiaPqGQ8uJsGo62s3VnI3OmjXE6pJgqzx2Dz+Ojsf0oHb2dZPkzHYkjnq3DVYInIiJRWlkSEYkhYwzXzguX4j2XgqV4Xo+XcXljATjQ7NwBvP44riyJiIhEKVkSEYmx6yMH1L68pYHu3qDD0cSeGzri+fLi1zpcREQkSsmSiEiMjRuTy9TxBXR2B3h9+yGnw4m56L4lJzviRVeWAq2xX1kyxmCMifm4IiKSfJQsiYjEQXR16bma1CvFc0NHPF92NhhDoK0dG0y91TsREXEHJUsiInFw7ZxKPAZq3jpMa0eP0+HE1MmzluodO3zXeL34cnLAWgJtbTEd21qbsocKi4jI0ChZEhGJg8K8DC4+v5RA0PLiZucaIcRDfkYeeek5dPR2crTjuGNx+LVvSURE4izpkyVjTLox5lvGmJAx5nqn4xERibp+frQU74DDkcSWMYYJkdUlJ89biudZSyIiIuCCZMkYM8YY83NjzI7Ix6PGmHGDfPZCYANwM3DW3bjGmM8YY7YbY7YYY2qNMXeNPHoRkTO7/MJy0tO8bN9zjPqm2JaKOe1kRzzn9mTF66ylhQsXsnDhwpiOKSIiycnRZMkYkwY8DaQBs4CZQDvwrDEmZxBDfBn4AvCdc8zzReCrwEJr7ezIM48YY94z/OhFRM4uK8PPVbPDzRD+8Np+h6OJrf77lpziy42sLMW4I966detYt25dTMcUEZHk5PTK0v3AbOAL1tqAtTZIOJGZAnx8EM9/1Fr7+7PdYIwpAP4e+J61dheAtfZp4CngWyOIXUTknG69bCIAf3x9P8FgyOFoYie6suRkR7x4rSytWbOGNWvWxHRMERFJTk4nS3cD+621u6MXrLWHgO2R187KWhsYxBy3A1nAs6ddfwaYaYy5YPDhiogMzczJRVSWZnOspZuaHUecDidmxuWVY4yhvvUwvcFeR2Lwx2nPksrwREQkyulkaTawZ4Dre4CLYjhHdMzT5+j/uohIzBljuPnS8OrS0xv2ORxN7KT70ijPGUPIhqhrcebgXZ+64YmISJw5nSyVAAP9K9cCZBljMmM0BwPME/1RZPFADxljlhhjqo0x1Y2NjTEIQ0RGq5uqxuPxGF7ffpjjrV1OhxMzfR3xHCrFi64sBWK8srRq1SpWrVoV0zFFRCQ5OZ0snclZO9slYg5r7SprbZW1tqq0tDQB4YhIqirMy+CSGWUEQ5Znq1OnjXhfRzyH2oefLMOL7crS0qVLWbp0aUzHFBGR5OR0stQE5A5wPRfosNZ2xmiO6JinzwFwNAZziIic1S2XTgDgqQ37sdY6HE1sTHS4yUO0DC8Q4254ixcvZvHixTEdU0REkpPP4fm3AAM1WJgMbI3hHACTgL2nzdH/dRGRuKmaUUZhbjp1jW28ufcYMycPWAGcVCbkh9ui73dqZSk3PitLKsETEZEop1eWHgcmGmMmRS8YY8qAGcBj/W80xpQZY4YT7++ADuD6067fAGy31r41jDFFRIbE6/VwY9V4AJ7ekBpnLpVmF5Ppy+BEVwsnOpsTPr83Ows8HoIdHYR6nenIJyIiqc3pZOknhFeQHjDG+CLJ0DcJd6r7fvQmY8xVQD3w3aFOYK09Afwj8EljzJTIeDcDtwGfG2H8IiKDdkvkzKUXN9fR0ZX839x7jIcpReHywp3H9iZ8fmPMySYPrW0xG7e+vp76eucO2xUREfdwNFmy1vYAtwBBwmcrvQnkATdaa/v/y9cGNAMN/Z83xtxrjNkEfC1y6YfGmE3GmCtPm+ebwD8B64wxW4B/Be6x1j4Z+3clIjKwytIcZk0ppqsnyAubUuOb8alFkwB45+heR+b397UPj92+pcrKSiorK2M2noiIJC+n9yxhrT0MfOgc92wGiga4/kvgl4Oc5zvAd4YeoYhI7Nxy6QTe2H2Up1/bx22XT3Q6nBE7vzi8/dOJlSUAXxwOpi0vL4/ZWCIiktycLsMTERlVrppdQWa6jx37jrP/UGy7uDkhurK089heQjaU8Pn9udGOeLFr8qAyPBERiVKyJCKSQBnpPq6dGy7xevq15G/0UJRVQFFmAZ29XdS3Hk74/L44lOGJiIhEKVkSEUmwWyONHp6pPkBvIPGrMbE2tXgSADsd2LfUt7IU4/bhIiIioGRJRCThzh9fwMSxubS09/Da9kNOhzNifaV4DiRL8dizNH/+fObPnx+z8UREJHkpWRIRSTBjTF8b8ac27HM4mpGLNnl459iehM99shte7FaWamtrqa2tjdl4IiKSvJQsiYg44Pp54/B5PWzccYRDR9udDmdEziucgDGG/Sfq6An0JHTuvnOWYriyVF1dTXV1dczGExGR5KVkSUTEAfk56VwzpwJr4TcvJX5FJpYy/BmMyysnaEPsOXEgoXP3leHFsBueyvBERCRKyZKIiEMWXD0FCHfF6+oOOBzNyJzv0OG00TK8WLYOFxERiVKyJCLikGkTCpk+sZD2zl6erT3odDgjMtWhw2l9udEGD7FLlpYvX87y5ctjNp6IiCQvJUsiIg6Kri6te3E31lqHoxm+8/vahye2pNCbmYHx+Qh1dRHs7o7JmCtWrGDFihUxGUtERJKbkiUREQddNbuCwtx09h9qZcvOJqfDGbZxeeWke9M40n6Ulq7ElcQZY042eWhti8mYy5YtY9myZTEZS0REkpuSJRERB/l9Ht5zxSQgvLqUrLweL1OKwu3QE16K19c+PDYd8VSGJyIiUUqWREQcdvsVk/B5Da+9cYjDxzqcDmfYoqV4iW/yEPv24SIiIqBkSUTEcYV5GVx9cSUhC79N4jbiUyMd8XYm+HBaX25sD6atqamhpqYmJmOJiEhyU7IkIuICC64Od5N7asM+unqSs4341L4mD3sT2qzi5J6l2KwsVVVVUVVVFZOxREQkuSlZEhFxgekTizh/fAFtnb2sr61zOpxhKc4spDAjn/beThrajiRsXn9ebFeW5s2bx7x582IyloiIJDclSyIiLrHwmuRuI26MOWV1KVF8edGzlmKzsqQyPBERiVKyJCLiEldfXEFBTjp7G1rYtvuo0+EMS9++pQQmS9GVpUAMD6YVEREBJUsiIq7h93m57Ypw++1kbSPe1xEvgU0e/DFeWRIREYlSsiQi4iLvuWISXo/h1W2HOHI8+dqITymaiMGw98RBeoK9CZkzes5SoDU2K0sVFRVUVFTEZCwREUluSpZERFykOD+Tq2ZXEApZnnx5r9PhDFmWP5PKvLEEQ0H2nTiYkDn9ubFdWWpoaKChoSEmY4mISHJTsiQi4jILrg43evj9q3vp6k6+NuJT+w6nTUwpnq9fN7xYNMaoq6ujri45OxKKiEhsKVkSEXGZCyYVMn1iIa0dvfzu1X1OhzNk5xeFz4xKVJMHb3o6nvR0bG8voa6uEY+nMjwREYlSsiQi4jLGGP70pmkA/N9zO+kNBB2OaGj62ocf25uwOWN91pKIiAgoWRIRcaWqGWVMKs/jWEsXz1QfcDqcIZmQX0Ga18+htkZau9sSMmcsz1pasmQJS5YsGfE4IiKS/JQsiYi4kMdjuOem8wF49Jl3CAZDDkc0eF6PlymFE4DErS75c2PXEW/16tWsXr16xOOIiEjyU7IkIuJSV11cSXlJNoeOdvDC5nqnwxmSqcWJ3bcUy5WllStXsnLlyhGPIyIiyU/JkoiIS3k9hrtvCK8uPfLHtwmFRt7pLVGmFk0CEriyFD1rKQZ7llSGJyIiUUqWRERc7Maq8ZTkZ7D/UCuvbT/kdDiDdn60ycPRvTFp530uvtxog4fYnLUkIiICSpZERFzN7/PwJ9dPBcKrS4lIPGKhJKuIwox8WnvaqWuJf5Ln7yvDG/nK0tq1a1m7du2IxxERkeSnZElExOVuvXwiedlpvL3/BJvfaXQ6nEExxnBh2XQAthx+M+7znSzDG/nK0qJFi1i0aNGIxxERkeSnZElExOUy0nzcee15ADzyx3ccjmbwZpfNAGDzofgnS30NHmLQDW/BggUsWLBgxOOIiEjyU7IkIpIE7rhqMlkZPrbsbOKtvcecDmdQLhp7AQDbG98hEAzEda5YriypDE9ERKKULImIJIHsTD93XBVux/2rP77tcDSDU5RZwPj8CroD3bx9dHdc5/Llxm7PkoiISJSSJRGRJHHnteeR5vfy+vbD7KlvdjqcQUlUKV50Zam3pSVpmmCIiIj7KVkSEUkS+Tnp3H75RCB59i5dPDacLMW7yYPH78ebmQmhEMH2jhGNZYzBGBOjyEREJJkpWRIRSSJ/cv1UfF7Di5vrkmJ1aUbp+fg8PnYf209rd1tc5/Ln5wPQc+J4XOcREZHRQ8mSiEgSKSnI5D1XTsZa+Mlvtjsdzjml+9KYXjIFi2XbkR1xnStjbBkAXQ0jO9fJWqtSPhERAZQsiYgknXtvnkZmuo/at44kxblLF4+dCcR/31JGeTkw8mRJREQkSsmSiEiSyc9J5+4bpwLwk3VvEAq5exVkdlm4hfjWQ2/GdcUmo3wsAJ31DXGbQ0RERhclSyIiSejOa86jKC+dnQebeWlzvdPhnNWkwvHkpmXT2HGMhrYjcZsnsyK6sjSyZGnhwoUsXLgwFiGJiEiSU7IkIpKEMtJ9fOi28IrNw09upzcQcjiiM/MYDxdFVpe2xLEUr29laYTJ0rp161i3bl0sQhIRkSSnZElEJEndfMkExo3J4dDRDn73yl6nwzmr2X0txN+K2xwZZWXg8dDd2ESot3fY46xZs4Y1a9bEMDIREUlWSpZERJKU1+vh/jvCzRP+9+kddHQNP0GIt+jhtG8c3kEgFIzLHB6/n/SSEgiF6Do8/HI/leGJiEiUkiURkSR22ayxzJhUREt7D48/u9PpcM6oJLuIitwyOgNd7Dy6N27z9O1bOqSOeCIiMnJKlkREkpgxhj9fMAuAXz+/i2MtXQ5HdGYnS/Hidz5UdN9S1wg64q1atYpVq1bFKiQREUliSpZERJLcjMlFXHFROd09Qf7n9/HbEzRS0VK8LYfiuG8pctbSSJo8LF26lKVLl8YqJBERSWJKlkREUsBH3jMDj8fw9Gv7OXC41elwBjRrzDS8xsM7x/bQ3tMRlzkyY7CytHjxYhYvXhyrkEREJIkpWRIRSQHjy3K55dIJhEKWn/4mfmVuI5Hpz2BayRSstbxx5O24zBFdWRrJniWV4YmISJSSJRGRFPGh2y4gPc3LhjcOUfPWYafDGdBFkVK8zYfik9BljC0DY+g60jii9uEiIiKgZElEJGUU5WXwoVunA/D9x7bQ1RNwOKJ3uzjO5y15/H7SS8Ptw7uPNA5rjPr6eurr62McmYiIJCPHkyVjzBhjzM+NMTsiH48aY8YN8lm/MeYfjTFvGWO2GWNeNsZcPcB9zxljthtjNp328dHYvyMREecsuvY8JpXncfhYB798Oj6lbiNxXuFEsv2ZHG5r5HDb8JKZcxlpk4fKykoqKytjGZKIiCQpR5MlY0wa8DSQBswCZgLtwLPGmJxBDPEfwL3ANdbaC4EfAU8bY+YMcO97rbVzTvt4OCZvRETEJXxeD5+852KMgf97bif7GlqcDukUHo+HWWXh1a94dcXra/IwzGSpvLyc8kjCJSIio5vTK0v3A7OBL1hrA9baIPAFYArw8bM9aIyZDiwBvmmtbQSw1v4Q2A38U1yjFhFxsQsmFnH7FZMIhizffXQzoZB1OqRTXFw2E4Ath9+My/gZkYNpO4fZEU9leCIiEuV0snQ3sN9auzt6wVp7CNgeee1s/gQwwLOnXX8GuHWQK1MiIinpo++dSWFuOm/uPcbTr+1zOpxTzB57AQDbDr9FMBSM+fgZYyMd8RqG3xFPREQEnE+WZgN7Bri+B7hoEM+GgP0DPOsjXNLX398YY16L7G963hjz58MJWEQkGeRk+vnYnRcC8ON12zne2uVwRCeV5ZRSnjuG9t5Oth3ZEfPxMyuiydLwz1oSEREB55OlEmCg0xNbgCxjTOY5nu2IlO6d/ixAcb9rJ4CdwA2E90Y9BHzfGPOtMw1ujFlijKk2xlQ3NsZnE7KISDxdM6eSedPH0N7Zy4/WvOF0OKe4esIlALyw97WYj51RNmZE7cPnz5/P/PnzYx6XiIgkH6eTpTMxsXzWWnuXtfYha227tTZorX0U+C/gs8aYCQMNYq1dZa2tstZWlZaWjiAcERFnGGP4+N2zSfN5eK72IBt3HHE6pD7XTLwUgA11m+gKdMd0bE9a2ojah9fW1lJbWxvTmEREJDk5nSw1AbkDXM8lvGrUeY5ns4wx3gGeBTh6jrk3EH7/lwwmUBGRZDS2OJsPRM9eenwL3b2x3yM0HGNzxzCteArdgW5eP7g55uNnjA13xBtO+/Dq6mqqq6tjHZKIiCQhp5OlLcCkAa5PBrYO4lkPMH6AZwPAmxBuT26MyR/g+eh3DKcnWyIiKeWu66YyviyXhqZ2HvmDe85eiq4uvbBvQ8zHPrlvaehNHlSGJyIiUU4nS48DE40xk6IXjDFlwAzgsf43GmPKjDH94/0/wALXnzbmDcBT1troXqgrgV8NMHf0X8KNww1eRCQZ+H0ePvn+iwF49Jl3eHv/cYcjCrtiwny8xsPmw29yoiu250FFD6ZVkwcRERkJp5OlnxBeQXrAGOOLJEPfJNzR7vvRm4wxVwH1wHej16y1O4BVwJeMMSWR+/4COA/4ymnz3GSMuaPfeNcDS4GfWWvfifm7EhFxmVlTill4zRSCIcu3/ruGjq6hNz6Itbz0HOaWX4i1lpf2vR7TsaPJUucwVpaWL1/O8uXLYxqPiIgkJ0eTJWttD3AL4ZK47YRL5/KAG621bf1ubQOagdN/RPgp4BHgJWPMNmAxcKu1dlO/e2qBzwNfNsZsNsbsBL4HfB34i5i/KRERl/qzO2YyqTyPhqPtrPy/c1U6J8Y1k6KleLHtipdZHt6z1DWMg2lXrFjBihUrYhqPiIgkJ5/TAVhrDwMfOsc9m4GiAa73Al+NfJzp2Rbg25EPEZFRK83v5e/um89nv/M8z1QfYO70MVw/b5yjMc0vv4hMfwa7j+/nYEsD4/LKYzJuxtiySPvwI4QCATy+wf9zt2zZspjEICIiyc/pMjwREUmgCWPz+g6r/d6jmzl0tN3ReNJ8aVwxbh4Q2zOXPGlppJcUR9qHD61lusrwREQkSsmSiMgoc/vlE7nionI6uwN8679rCARDjsZzzaTLAHhx32uEbOxiOdnkYej7lkREREDJkojIqGOM4VN/OoeS/Ax27D/O//z+LUfjmVE6leKsQho7jrGjaVfMxu1r8jDEfUs1NTXU1NTELA4REUleSpZEREah3Kw0/ubD8zEm3E58y85Gx2LxGM/JM5diWIo33LOWqqqqqKqqilkcIiKSvJQsiYiMUhedV8Kf3jQNa+Hb/1NLS3uPY7FEk6VXDtTQE4xNW/OMSEe8ziGetTRv3jzmzZsXkxhERCS5KVkSERnFPnjrdC6YWMjR5i4e+uVGQiHrSBzj8yuYVDCO9t5ONjZsi8mYmcM8mFZleCIiEqVkSURkFPN6PXzuviqyM3xseOMQ//OUc/uXro00eohVKV5f+/DD4fbhIiIiQ6VkSURklCsryuLvPlKFx8Avn36b52oOOBLHlROqMMZQ07CVtu6RtzQ/pX14o3N7skREJHkpWRIREeZfUMbiuy4C4N9/uYnte44mPIaizAIuGnMBwVCQVw7UxmTMvvbhQ+iIV1FRQUVFRUzmFxGR5KZkSUREAFhw9RTuuGoygWCIf/rxa44cWNtXirdvQ0zG62sfPoSOeA0NDTQMcZ+TiIikJiVLIiLSZ/GdFzJv+hha2nv42n9toL0zNp3pBuvSyotJ96bxVtMu9p+oG/F4mZGOeENZWaqrq6OubuRzi4hI8lOyJCIifbxeD5//SBXjy3I5cLiVBx5+nWAwlLD5M/wZXD/5CgAe2/7kyMeLJkuHVIYnIiJDp2RJREROkZ3p5x/+8jLyc9LY+HYjq5+ITSvvwbprxm34PD5ePVDLweaRlcNFD6btHMLKkoiISJSSJREReZexxdl85c8uw+/z8JuX9rD2hd0Jm7s4q5Abp1yJxfLo9t+OaKz0sjIAuo80Drp9+JIlS1iyZMmI5hURkdSgZElERAY0Y3IR/+/euQCsfmIrz1TvT9jcd824Da/Hyyv7a0a0uuRNTyetpAQbDA66ffjq1atZvXr1sOcUEZHUoWRJRETO6Pp54/izO2ZiLXznfzcmLGEqySrixsnh1aXHRri61NfkYZAd8VauXMnKlStHNKeIiKQGJUsiInJWd994Ph9974x+CVNiDq39kxm34/V4eXl/DXUtg2/9fbqMIe5bUhmeiIhEKVkSEZFzuuemaXzkPdGEqZZna+KfMJVkF3FD3+rS8DvjZYyNriypyYOIiAyNkiURERmUP715Gve95wKshQd/kZiE6U8ie5de2v869cNcXYp2xBtssrR27VrWrl07rLlERCS1KFkSEZFBu/fm6dx3ezhh+s4vankuzglTaXYx10+6AmuHv7qUUR4pwxvknqVFixaxaNGiYc0lIiKpRcmSiIgMyb23hBOmUGSFKd4J05/MvB2v8fDi/tepbz085Oczxkbahx8+gg0Gz3n/ggULWLBgwZDnERGR1KNkSUREhuzeW6bz4UjC9O1f1PL4s+9grY3LXGOyi7lucnh16fFhrC71bx/ecbDunPerDE9ERKKULImIyLB84JbpfW3Ff7xuO//xq030BkJxmet9M8KrSy/se42G1iNDfr5g9oUAHNvwWqxDExGRFKZkSUREhu3uG8/ni/dfQprfy9Ov7WfZqldo7eiJ+Txjckq4dtLlw15dKr7yCgCOvvxKrEMTEZEUpmRJRERG5KrZFXzzk1dRmJvO1l1NfO7fn6eusS3m87xv5u14IqtLu47tG9KzBXMuxpuZSfuevXTW15/1XmMMxpiRhCoiIilCyZKIiIzY+eML+bdPX8eUinzqm9r53L8/z9adTTGdoyynlPecfwMhG+I/Xv0x3YHBr2B5/H6KLr0EgKMvvxrTuEREJHUpWRIRkZgoLczkm399NZfOHEtbZy9/v/Jlfv/q0FaAzuWDFy1ifF459a2H+dnmx4b0bPGVlwPQdI5kyVobt2YVIiKSXJQsiYhIzGSm+/jyn1/KXdedRzBk+c9HNvGv/11NW2dvTMZP86Xxqcv/Ap/Hx1M7n6e2ftugny2YOwdPRgbtu3bRdXjoLchFRGT0UbIkIiIx5fUY/nLRhXz63rlkpHl5fmMd/+/fnmXbrtiU5U0qHMcHLloIwPdfe5jmrpbBxZWeTlHVfECleCIiMjhKlkREJC5uvnQC//4313P++AIaj3fy5e+/xMO/3U4gOPL24gum3cysMdNo7m5l5es/H3TZXLQrXtNZuuItXLiQhQsXjjhGERFJfkqWREQkbipKc/iXT13Dn948DQM88sd3+Lv/eGHE3fI8Hg+fvPR+svyZVNdv4ZndLw3qucL5c/GkpdH29jt0Nw680rVu3TrWrVs3ovhERCQ1KFkSEZG48nk9fOQ9M/jGJ65mTGEmOw+c4NPffo4nX9lLKDT8Rgol2UV8bP4HAPjJxkcGdVitNyODwvnzADj66sCleGvWrGHNmjXDjktERFKHkiUREUmIWVOKeehvb+C6uePo7gnyvUc38/n/eIF3Dhwf9phXT7yUqyZU0R3s4T9e/TGBUPCcz5w8oHbgZElleCIiEqVkSUREEiY708/n7pvP3903n8LcdHbsP87f/vvz/Ocjm2hu6x7WmH85/wMUZxWy89heHt/+5DnvL6yaj/H7aXnzLbqPHhvWnCIiMjooWRIRkYS7du44fvDFm/iT66fiMYbfv7qPv/rmH/nNi7sJDrEBRE5aNp+89H4Mhse2/5YX97121vt9WZkUzpsD1nLs1Q3ven3VqlWsWrVqSDGIiEhqMjp479yqqqpsdXW102GIiKSkA4dbWfXrrWx6uxGAyRV5LLnrIi48r2RI4zy+/Un+d+sajDF88tL7uXbSZWe898hz63nnwYfIu3AWF/3T1055zRgDoINpRURGFzPQRV+ioxAREelvfFkuX1tyBa9ua+CHT2xjT30LX/reS8w5v5QP3DqdWVOKBzXO+2a+h5C1/GrbWr674aeEbIjrJ18x4L1Fl1RhfD5atr9Jz4kTpBUU9L22ePHiWLwtERFJAVpZGgStLImIJEZXT4D/e3Ynv35+Fx1dAQAuOq+ED9w6jYvOK+lb9TmbvhUmDEsvuY8bp1w54H3b//EbHK+u4byPL2Xs7bfG9H2IiEjSGfAfGO1ZEhER18hI8/HB2y7gv75yCx+6dTrZmX627mriK99/mS9+90U27jhyzvK49818Dx+afRcWyw9e/xl/2PXigPcVX3k5cPYDakVEZHTTytIgaGVJRMQZ7Z29rHtpN0+s30VrRy8A0yYUcMdVk7nq4krS/d4zPrvmraf5782PA/Cx+R/k1qnXnvJ6b2srr9//l1hrufSn/4U/Lw+A+vp6ACoqKuLxlkRExJ0GXFlSsjQISpZERJzV0dXLb17aw6/X76KlvQeA3Cw/N10ygduvmERlac6Az63b8Uce3vQoAH8x715uP//6U15/Y/k/cmLjJqb+9ccpu+VmQA0eRERGKSVLw6VkSUTEHbq6Azy/qY4nX9nLzgMn+q5ffH4J77liMpddOBaf99QK89++/Qw/2fgIAFdPuIS/mHcvOenZABx66g/s+u73KZg3l1nLvgqcXFGKrjCJiMiooGRpuJQsiYi4zzsHjvPky3tZv7GOnt4gAAU56Vwxu5xr5lQyc3IxXk/4375ndr/Mj2t/SXewh8KMfJZe8mHmVVxEb0sLr93/lxhjmPf9/yCjrMzJtyQiIs5RsjRcSpZERNyrrbOXZ6sP8OQrezlwuLXvelFeOlfOruCaOZVcMLGII+2NfPe1h9nRtAuAGyZfyf1z38/+B7/L0ZdeIWvCeC564Bv4srKceisiIuIcJUvDpWRJRMT9rLXsrmvmxc31vLCpjsPHOvpeK87P4MrZFcydXsKB0BYe276O3lCAkqwi/urC9xP41k/oPHiQwvnzmPGVL2K8Z24cISIiKUnJ0nApWRIRSS7WWnYePMELm+p5cXMdjcc7+17z+zxMneqlpeg1jgUOA/Dewrlc8NOXCLa1sXRTDeklxdTU1DgVvoiIJJ6SpeFSsiQikrystezYf5zXtx+mdseRk40hTAhf+W78lbvAWMYdCfAnzxzn2l8/1veciIiMGkqWhkvJkohI6mhu62bj241s3HGE2h1HaA404avYhbfoEDP3dDLhd/uxBtLuu5+517+XaROLznqek4iIpAQlS8OlZElEJDWFQpa9DS1s3dXEpr17eaujmsv2vsElb7bT5Tf87zUTaGqfTnnaZM6rKGZKRT5TKvKZXJlPXnaa0+GLiEjsKFkaLiVLIiKjg7WWnQ2HefNfvk7xngaO53r55a2FdPl8hFpKCB4fQ/DEGAikUVKQyaTyPMaX5VJZmsO4MeGP/Jx0p9+GiIgM3YDJki/RUYiIiLiVMYafr/oBoex87pqcBnv2cc/LXfxmXhrHC4/gLTwC1mDbCzl+dAxHd5dS/WYW/f+Nzc1KY9yYHCpLcygrzmJMYRZlReGPorwMPJ4B/z0WEREX0srSIGhlSURk9DAmnMx0NTax+XNfoPf4cQACMyaxbVY+L6YfJmhDffenmQwyQyWE2vJpbcyk80QuBAZeXfJ5PZQWZjKmMJPi/EyK8zMozsugKD+j7+uCnHS8Xk/836iIiPTnzpUlY8wY4EGgKnJpK/AZa+3BQTzrB/4BuAcIAC3A5621Lw5w72eAJZH7AsDXrLW/jsFbEBGRFLJs2TIA0kuKmf3ANzj42P/R+Oxz+N7cy5w34fKJE+i6+kKqxwbYemwnLd1t9JiDkHsQciETyEvLp8A7Bl8gl2BnJp0t6Zxo8tF8wtLQ1E5DU/sZ5zcmvDqVn5NOQU46+TnhX+dHfp2blUZulp+crJO/zkz39SV5IiISO46uLBlj0oDXgbeBDwIW+BFwJTDXWtt2jud/ANwIXGWtbTTGfAz4D+AKa+2mfvd9EfgccJm1dpcx5hbgt8Aia+2T54pTK0siIqNbb0sLh373FA2/fZLe4ycA8BcUMOamGwiNG8OhnBC7va3sbD7A7uP76Qp0DziO3+OjKKOIbG8+fpsFgQyC3X66O3x0tHlpPWFobgZCXs7wQ84BeTyGnEw/WRk+sjL8ZGdEf+0jO8NPZoaPzHQfGWk+MtO9ZPT92kdGmpf0NC/pfh/paV7S/B7SfF6VC4rIaOO+Bg/GmMXAKuA8a+3uyLWxQB3wRWvtv57l2enAm8DHrLU/6nf9DWCvtfaOyNcFkfH+zVr7D/3u+w0wyVo761xxKlkSERGAUG8vTS+8SN0Ta+nYu++U14zXS0Z5OZnjKgmMKaA518cJ081R28HhYCsNgWaabAc9foM9RyLiNV4yfBmkezLwkYbXpkHIjw34CAY8BHo99PYYenqgpwt6egw26AXrgZAHrAcb+dz3tTXhr62JfER+fYakLM3vJd3vIc3vJc3nxe/3kObz4PeFEyq/z4vf58Hn9ZzyOfprr9fg93rwesNf+7wGr9eD32vweCJfe8L3eT0nf+3xhL+OfvZ6PHg8Bo8JJ4XhX/e7L/Jrj8dgTOQ+YzCefr82BmPQ6puInI0ry/DuBvZHEyUAa+0hY8z2yGtnTJaAPyH8pp497fozwF8ZY3IiK1O3A1lnuO9bxpgLrLVvjfB9iIhIiqipqQFg/vz573rN4/cz5sYbKL3hepq3buNE7UY6Dh6k88BBug4fofPgQToPhqvIDVAY+Zh6+kB+H9brIeT1EPJA0AMBAwETohdLCAuRnCZkAGOwRHIbwmUYfZMM8PUp1xjgxbOx/e8zp45zyhdnGy/8WgAIWBh4nW0oEpvkmIHmHGYI8Yp88HmfEkRxr2Canw//52qnwzgrp5Ol2YRL8E63B7hpEM+GgP0DPOsDZgKvRe6LXj/9vug470qWjDFLCO9xYsKECecIRUREUkVVVXgL7dkqL4wxFMy+iILZF/VdC3Z301lXR8f+g3QeOEB3YxOBjnaC7R0EOjoItreHP3d0Qm8A0wtewh/+OL8nERE36kh3fzMbp5OlEqBmgOstQJYxJtNa23mWZzustcEBngUo7ncfQOs57juFtXYV4RJBqqqq1DJQRGSUmDdv3rCe86ankzNlCjlTppz1PhsKEeruxgaDhAIBbCCIDfQSCgSxgQA2GARrw8laKIQNhSJfhyBkT03irA1/nD7H6dfOkPhZIGRDWGvDn4l8tmBtiFAoFFlMCq922UhckbUvsOHnT04RfcWeMq3lZJwnX42+YvtdsP3+9/TQz/xP8WC2FPTNa23k/UXm5tTfxuhI1kbisrbv2eh7tvbUOPu+siefPdOKnH33pVPGfHfUEDrb2zvLH7U92+/ZWYYc7ovD3dlxtjhdJUnCjId4vXVfmvsP93Y6WTqTkawZD/ZZrUuLiMi7RMvw4sV4PHgzM+M6h4iIxIbTa19NQO4A13MJrxqdaVUp+myWMcY7wLMAR/vd1//6me4TERERERHp43SytAWYNMD1yYTPWzrXsx5g/ADPBgh3yovexwDzTD7tdRERERERkT5OJ0uPAxONMZOiF4wxZcAM4LH+Nxpjyowx/eP9P8IllNefNuYNwFPW2ugepd8BHWe4b7s64YmISH8VFRVUVFQ4HYaIiLiA08nSTwivID1gjPFFkqFvEu5U9/3oTcaYq4B64LvRa9baHYQbMHzJGFMSue8vgPOAr/S77wTwj8AnjTFTIvfdDNxG+KBaERGRPg0NDTQ0NDgdhoiIuICjDR6stT3GmFuAB4HthFeKtgE3Rs5IimoDmoHT//X6FLAMeMkY00u4492t1tpNp83zTWNMF7DOGBMAgsA91ton4/C2REQkidXV1TkdgoiIuIQZTLvN0a6qqspWV1c7HYaIiIiIiMTHgJ2ynS7DExERERERcSUlSyIiIv0sWbKEJUuWOB2GiIi4gMrwBkFleCIio4cx4UoM/fsoIjKqDFiG52iDBxEREbdZuXKl0yGIiIhLKFkSERHpRyV4IiISpT1LIiIiIiIiA1CyJCIi0s/atWtZu3at02GIiIgLqAxPRESkn0WLFgFq8CAiIkqWRERETrFgwQKnQxAREZdQsiQiItKPSvBERCRKe5ZEREREREQGoGRJRERERERkAEqWRERE+jHGYMyAB7mLiMgoo2RJRERERERkAEatUc/NGNMI7HM6DjmjEqDJ6SAkKejvigyW/q7IYOnvigyG/p64X5O19vbTLypZkqRnjKm21lY5HYe4n/6uyGDp74oMlv6uyGDo70nyUhmeiIiIiIjIAJQsiYiIiIiIDEDJkqSCVU4HIElDf1dksPR3RQZLf1dkMPT3JElpz5KIiIiIiMgAtLIkIiIiIiIyACVLIiIiIiIiA1CyJCnLGPMxY4w1xix3OhYREREZ3YwxL0S+L5nkdCwyeEqWJCUZY7KBrzkdh7iTMabQGPM5Y8zLxphaY8ybxpinjDFXOx2bOMcYM8YY83NjzI7Ix6PGmHFOxyXuYoyZY4xZ/f/bu9dQy8o6juPfH46ioqbmjaIo6TJljaMppKIEqYQNIyUkoWmFVyqoFzGal5xelCVpWZQmkYVCyRRhFwyjptDyMsZoOmEvGkussdJMm1HMmX8v9ho5jut0zplzPM86e74f2Oyzn/2svX4HFnvv/36e9awk9yS5N8m6JFcn2b91Ng1XklMAP2MWIIsljasVwO9ah9BgnQtcAJxTVYcDbwHWAb9KcmLTZGoiyS7ArcAuwCHAm4GNwC+T7NEymwbnu8C+wHFVdShwAnAicHuS3Zom0yB17y+fA37aOotmzmJJYyfJK4GzgJWts2jQrqmq+wGqajOj4mkz8LGmqdTKmcASYEVVPdcdEyuAg4HzmybTEK2oqo0AVfUIcAXweuCkpqk0VB8B1gB3tw6imVvUOoD0EvgscDXwj9ZBNFhXbNtQVc8k+RewT4M8au8U4C9V9aetDVW1Icm67rkXHTPaYS2pqme3aftrd+/7h14gyb7AJ4GjgQ+2TaPt4ciSxkqSw4DjgC81jqIBq6rN3cjB85LsA+wHrG4SSq0tAdb3tK8H3jrPWTRgPYUSwBuAAn49z3E0fJcCN1TVQ62DaPs4sqRx80Xgkm6UoHUWLSxnAX8HrmodRE3sB9zT0/4ksHuS3arq6XnOpAUgyU7Ah4FvVtUfW+fRcCR5HfA+4E2ts2j7ObKkwUpyfLfE5lS31V3/5cBewI1Ng2vezfRY6dn+EEbnLL2/qh6b1/AaOn910VQuAZ4DPtE6iAbnC8DlVfXv1kG0/RxZ0pD9hun9GrMpySJGb0rnVVW9tLE0QNM+VrZtSPJq4EfA2VW1eo5zaeH4J7BnT/uewCZHldQnyYcYjRy8o6r+0zqPhiPJsYxWWj21dRbNTvxeqXHQjQzcymga1Va7MPoC/SiwAbijqs5rEE8DleRVjI6bi6tqVes8aifJLcDiqnrNNu2/BzZW1dubBNNgJfkAoxUTj6+qDa3zaFiSfAY4DXhqQvNBwIHAH4BngU9VlcuJD5zFksZWd4Xs9cDKqrqsbRoNTXex0Z+zTaGU5OaqWt4umVpIcg5wLfDarSdiJzkQeAS4sKpcDU/PS3I6o6m7zxdKSZYBr6iqbzQNp8FKchnwaSa8z2j4nIYnaYfTFUqrgbXArt0Xn62WtMik5q4HPgp8PslpwBbgckY/uHy9YS4NTHd8XMfoXKXjJywmdCzwt1a5JL00HFnSWEqyBtidF07Du7Sqbm4aTIOQ5EomPxn7z9tOxdKOoRtJugo4gtEy0PcDH6+qh5sG06AkeZzJr6fkTAa9SJKTGF0D8gXT8Kpqactcmh6LJUmSJEnq4dLhkiRJktTDYkmSJEmSelgsSZIkSVIPiyVJkiRJ6mGxJEmSJEk9LJYkSZIkqYfFkiRJkiT1sFiSJEmSpB4WS5IkSZLUw2JJkrQgJVmcpKZx+8n/eY3dkzyd5Ngp9nVGkruT3JXktEn6LElyQ5L7kqxN8kCS1UkuSPLG2f6/kqT5l6pqnUGSpBlLsjeweELTEcBXgIuAX0xo31BVD03yGsuB64H9q2rzJH2OAlYBb2P0I+Ma4D1VdeeEPmcCXwbOAVZV1ZYkAd4LfAvYUlV7z/iflCQ1tah1AEmStkdVPQHcsfVxkiO7P39YVeum+TLLgFsmK5Q6RwO3VdWGbj+3AccAd3aPDweuA86vqpsm5Cvg+0leBlw5zTySpAFxGp4kaVwsBZ4BHpzBNicBP56izxrgmCQHJjmICYVS52JgI/CdSbb/HqORKUnSAuM0PEnSWEiyhtGAzpFTdh71Pwy4Gzigqh6fou+5jKbYbQG+WlXf7tp3Ap4A7qqqd84iviRpgJyGJ0la8Lqi5RDgxhlstgz47VSFEkBVXQtc2/PUy4E9gEdnsF9J0gLhNDxJ0jhYDOwKrJ3BNsuYegredDlNQ5LGkMWSJGkcLO3u751O5yQHMFo9b7bF0mOMzlc6aJavI0kaIIslSdI4OLS7n1axxGhhh4er6oHZ7LRbRe9nwOFJdu7rk2TfJO9Kstds9iVJmn8WS5KkcbAUWF9VT06z/1xOwVsJ7AacPsnzFwFfAzbN0f4kSfPEBR4kSePgUOD26XTsRoBOAE6dix1X1X1JzgCuSfIU8IPuorQ7A+cCZwPLq+q5udifJGn+WCxJkha07tpHBzD9KXjHMfr8Wz1XGarqpiQPAhcCK5P8l9HsjbXAUbOd7idJasPrLEmSdihJrgIOrqqTW2eRJA2b5yxJknY072buzleSJI0xR5YkSZIkqYcjS5IkSZLUw2JJkiRJknpYLEmSJElSD4slSZIkSephsSRJkiRJPSyWJEmSJKmHxZIkSZIk9fgf7lrvUw8yTCUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "T = np.linspace(-5,5,100)\n",
    "fig, ax = plt.subplots()\n",
    "k=2\n",
    "ax.plot(T,phi_I(T+273.15),label='$k=2$')\n",
    "k=4\n",
    "ax.plot(T,phi_I(T+273.15),label='$k=4$')\n",
    "k=10\n",
    "ax.plot(T,phi_I(T+273.15),label='$k=10$')\n",
    "#ax.plot(T,dphi_I_dT(T+273.15),label='$k=10$')\n",
    "ax.axvline(0,ls=':')\n",
    "ax.set_xlabel('$T$ / °C')\n",
    "ax.set_ylabel('$\\\\phi_\\\\mathrm{I}$')\n",
    "ax.legend()\n",
    "fig.tight_layout();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def heat_capacity(T,phi=0.35,T_m=273.15):\n",
    "    c_pW = 4186. #J/kgK\n",
    "    c_pI = 2052. #J/kgK\n",
    "    c_pS = 500. #J/kgK\n",
    "    rho_W = 1000. #kg/m³\n",
    "    rho_I = 900. #kg/m³\n",
    "    rho_S = 2600. #kg/m³\n",
    "    L = 334.e3 * 1000/900 #J/kg renormalized to ice\n",
    "    phi_ice = phi_I(T,phi,T_m)\n",
    "    rho_cp = rho_S * (1-phi) * c_pS + rho_W * (phi - phi_ice) * c_pW + rho_I * phi_ice * c_pI\n",
    "    #rho_cp = rho_S * (1-phi) * c_pS + rho_W * phi * c_pW\n",
    "    return rho_cp - dphi_I_dT(T,phi,T_m)*L*rho_I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enhalpy of fusion is 338.127 MJ/m³\n"
     ]
    }
   ],
   "source": [
    "T = np.linspace(-1.,.5,1000)\n",
    "H = 0.\n",
    "k=50\n",
    "for i in range(len(T)-1):\n",
    "    H += (T[i+1]-T[i])*heat_capacity((T[i+1]+T[i])/2.+273.15,phi=1)\n",
    "print(\"Enhalpy of fusion is %.3f MJ/m³\" %(H/1e6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def heat_conductivity(T,phi=0.35,T_m=273.15):\n",
    "    lam_W = 0.3 #W/mK\n",
    "    lam_I = 2.5 #W/mK\n",
    "    lam_S = 2.2 #W/mK\n",
    "    phi_ice = phi_I(T,phi,T_m)\n",
    "    #return lam_S * (1-phi) + lam_W * phi\n",
    "    return lam_S * (1-phi) + lam_W * (phi - phi_ice) + lam_I * phi_ice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time and space discretization, Picard iterations\n",
    "\n",
    "Using a backward Euler approach for simplicity we find the time-discrete weak form as ($p_{n+1} \\equiv p$)\n",
    "\n",
    "$$\n",
    "    0 = \\int \\limits_0^H \\left[\\eta (C - L \\partial_T \\varrho_\\text{I}) \\frac{T - T_n}{\\Delta t} + \\eta_{,z} \\lambda T_{,z} \\right] \\, \\text{d}z + \\left[ \\eta q_z \\right]^H_0\n",
    "$$\n",
    "\n",
    "Introducing standard FE approximations\n",
    "\n",
    "$$\n",
    "    T \\approx N_i \\hat{T}_i, \\quad \\eta \\approx N_i \\hat{\\eta}_i, \\quad \\frac{\\partial T}{\\partial z} \\approx \\nabla  N_i \\hat{T}_i, \\quad \\frac{\\partial \\eta}{\\partial z} \\approx \\nabla N_i \\hat{\\eta}_i\n",
    "$$\n",
    "\n",
    "This yields\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    " 0 &= \\int \\limits_0^H \\left[ N_i \\hat{\\eta}_i (C - L \\partial_T \\varrho_\\text{I}) \\frac{N_k \\hat{T}_k - N_k \\hat{T}_{n,k}}{\\Delta t} + \\nabla N_i \\hat{\\eta}_i \\lambda \\nabla N_k \\hat{T}_k \\right] \\, \\text{d}z + \\left[ N_i \\hat{\\eta}_i q_z \\right]^H_0\n",
    " \\end{align}\n",
    "$$\n",
    "\n",
    "Now we bring all quantities associated with the unknown pressure to the left-hand side (LHS) and all known quantities to the RHS:\n",
    "\n",
    "$$\n",
    "    \\hat{\\eta}_i \\int \\limits_0^H \\left[ N_i \\frac{(C - L \\partial_T \\varrho_\\text{I})}{\\Delta t} N_k + \\nabla N_i \\lambda \\nabla  N_k \\right] \\, \\text{d}z\\ \\hat{T}_k =  \\hat{\\eta}_i \\left[  \\bar{q}_z|_{z=H} \\delta_{n_\\text{n}i} -  \\bar{q}_z|_{z=0} \\delta_{i0} \\right] +  \\hat{\\eta}_i \\int \\limits_0^H N_i \\frac{(C - L \\partial_T \\varrho_\\text{I})}{\\Delta t} N_k \\, \\text{d}z\\ \\hat{T}_{n,k}\n",
    "$$\n",
    "\n",
    "can be simplified by realizing that the nodal test function values are arbitrary and thus\n",
    "\n",
    "$$\n",
    "    \\int \\limits_0^H \\left[ N_i \\frac{(C - L \\partial_T \\varrho_\\text{I})}{\\Delta t} N_k + \\nabla N_i \\lambda \\nabla N_k \\right] \\, \\text{d}z\\ \\hat{T}_k =  \\bar{q}_z|_{z=H} \\delta_{n_\\text{n}i} -  \\bar{q}_z|_{z=0} \\delta_{i0} +  \\int \\limits_0^H N_i \\frac{(C - L \\partial_T \\varrho_\\text{I})}{\\Delta t} N_k \\, \\text{d}z\\ \\hat{T}_{n,k}\n",
    "$$\n",
    "\n",
    "which leaves us with $n_\\text{n}$ equations for the $n_\\text{n}$ unknown nodal pressures $\\hat{p}_k$. If any coefficients in the above are taken as pressure-dependent, the system could be solved repeatedly usind Picard iterations. For strong non-linearities, a Newton linearization would typically be used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we require now is the local assembler to calculate the stiffness matrix and the local right-hand side. Local integration is performed by Gauss quadrature:\n",
    "\n",
    "$$\n",
    "    \\int \\limits_{-1}^1 f(\\xi)\\,\\text{d}\\xi \\approx \\sum \\limits_{i=1}^{n_\\text{gp}} f(\\xi_i) w_i \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Local assember"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def local_assembler(elem,dt,prev_sol,sol,mass_lumping=False):\n",
    "    element_order = elem._line_element__nnodes\n",
    "    K_loc = np.zeros((element_order,element_order))\n",
    "    M_loc = np.zeros((element_order,element_order))\n",
    "    b_loc = np.zeros(element_order)\n",
    "    z_nodes = elem._line_element__coords\n",
    "    for i in range(elem._line_element__quad_degree):\n",
    "        #local integration point coordinate\n",
    "        xi = elem._line_element__quad_points[i]\n",
    "        #shape function\n",
    "        N = shape_function(element_order,xi)\n",
    "        #gradient of shape function\n",
    "        dN_dX = grad_shape_function(elem,xi)\n",
    "        #determinant of Jacobian\n",
    "        detJ = np.abs(element_jacobian(elem,xi))\n",
    "        #integration weight\n",
    "        w = elem._line_element__quad_weights[i]\n",
    "        \n",
    "        #global integration point coordinate (for spatially varying properties)\n",
    "        #z_glob = np.dot(N,z_nodes)\n",
    "        #evaluation of local material/structural properties\n",
    "        #E = Stiffness(z_glob)\n",
    "        #evaluation of local body force\n",
    "        #CV = ConsolidationCoeff(z_glob)\n",
    "        \n",
    "        T_prev = np.dot(N,prev_sol)#T_prev in integration point\n",
    "        T = np.dot(N,sol)#T in integration point\n",
    "        \n",
    "        C = heat_capacity(T)\n",
    "        lam = heat_conductivity(T)\n",
    "        #assembly of local stiffness matrix\n",
    "        M_loc = np.outer(N,N) * C / dt\n",
    "        if (mass_lumping):\n",
    "            M_loc = np.diag(M_loc.sum(0)) #diagonal of column sum\n",
    "           \n",
    "        K_loc += (np.outer(dN_dX,dN_dX) * lam + M_loc)* w * detJ\n",
    "        \n",
    "        #assembly of local RHS\n",
    "        b_loc += N * C * T_prev/dt * w * detJ\n",
    "    return K_loc,b_loc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time loop and problem solution\n",
    "\n",
    "We now establish the time loop and in each time step perform the global assembly, apply a vanishing traction on the top and constrain the displacement at the bottom to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def time_loop(dt,nodes,elements,solution,mass_lumping=False):\n",
    "    #Startwerte\n",
    "    t_end = 366*24*60*60 * 500#s\n",
    "    absolute_tolerance = 1.e-6\n",
    "    max_iter = 100\n",
    "    iteration_counter = np.array([0])\n",
    "    apply_initial_conditions(solution,268.15)\n",
    "    y = [solution] #create a list that will hold the solution vectors at all time points\n",
    "    times = np.array([0.])\n",
    "    #\n",
    "    while times[-1]+dt < t_end: #repeat the loop as long as the final time step is below the end point\n",
    "        times = np.append(times,times[-1]+dt) #here define the next time point as the previous time point plus the time increment dt\n",
    "        y_old = y[-1] #Starting value for recursive update\n",
    "        i = 0\n",
    "        #\n",
    "        while True:\n",
    "            K, f = global_assembler(nodes,elements,y[-1],y_old,dt,mass_lumping)\n",
    "            #heat flux bc\n",
    "            f = apply_Neumann_bc(f,0,0.2)\n",
    "            #alternative: variable surface temperature\n",
    "            #value = 278.\n",
    "            #if (times[-1] > 366*24*60*60*300):\n",
    "            #    value = 268.\n",
    "            #K, f = apply_Dirichlet_bc(K, f, 0, value)#fixed temperature top\n",
    "            solution = np.linalg.solve(K,f)\n",
    "            i += 1\n",
    "            if (np.abs(np.linalg.norm(solution) - np.linalg.norm(y_old)) < absolute_tolerance or i > max_iter): #if change is below tolerance, stop iterations\n",
    "                break\n",
    "            y_old = solution #preparation of next recursion\n",
    "\n",
    "        y.append(solution) #append the new found solution to the solution vector\n",
    "        iteration_counter = np.append(iteration_counter,i) #store how much iterations this time step took to converge\n",
    "    return times, y,iteration_counter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "#spatial discretization\n",
    "H = 150.\n",
    "nel = 20\n",
    "n_per_el = 3\n",
    "nodes,elements,solution=generate_mesh(H,nel,n_per_el)\n",
    "k = 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "times, sols, iters = time_loop(24*60*60*365,nodes,elements,solution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(ncols=2,figsize=(10,10))\n",
    "ax[0].set_ylabel('$z$ / m')\n",
    "ax[0].set_xlabel('$T$ / °C')\n",
    "ax[0].plot(sols[100]-273.15, -nodes, marker='o', label='$t = %i$ a' %(times[100]/60/60/24/365))\n",
    "ax[0].plot(sols[300]-273.15, -nodes, marker='o', label='$t = %i$ a' %(times[300]/60/60/24/365))\n",
    "ax[0].plot(sols[500]-273.15, -nodes, marker='o', label='$t = %i$ a' %(times[500]/60/60/24/365))\n",
    "ax[0].axvline(0,ls=':')\n",
    "\n",
    "ax[1].set_ylabel('$z$ / m')\n",
    "ax[1].set_xlabel('$\\\\phi_\\\\mathrm{W} / \\\\phi$')\n",
    "ax[1].plot((0.4-phi_I(sols[100]))/0.4, -nodes, marker='o')\n",
    "ax[1].plot((0.4-phi_I(sols[300]))/0.4, -nodes, marker='o')\n",
    "ax[1].plot((0.4-phi_I(sols[500]))/0.4, -nodes, marker='o')\n",
    "fig.legend(loc='center right',fontsize=12);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}