{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Importing Snapshots from OpenFOAM\n",
    "\n",
    "**Aim of this notebook:** learn how to import parametric snapshots from OpenFOAM (version 6 from org-version) using *fluidfoam* package.\n",
    "\n",
    "The snapshots are related to the Buoyant Cavity problem in fluid dynamics, governed by the Navier-Stokes equations, including energy, under the Boussinesq approximation. In particular, the snapshots have been generated using the case reported in [ROSE-ROM4FOAM tutorials](https://ermete-lab.github.io/ROSE-ROM4FOAM/Tutorials/BuoyantCavity/problem.html).\n",
    "\n",
    "*Disclaimer*: the OpenFOAM snapshots are not included in the Zenodo repository for storage issues, but they can be generated using the case reported in the link above. Therefore, to generate the OpenFOAM case and later execute this notebook, you need to have OpenFOAM-v6 installed in your machine.\n",
    "For interested readers, please contact stefano.riva@polimi.it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from mpi4py import MPI\n",
    "from dolfinx.io import gmshio\n",
    "import gmsh\n",
    "from dolfinx.fem import Function, FunctionSpace\n",
    "import ufl\n",
    "\n",
    "from pyforce.tools.functions_list import FunctionsList\n",
    "from pyforce.tools.backends import LoopProgress\n",
    "\n",
    "path = './Snapshots/OpenFOAM/'\n",
    "\n",
    "importing = ['TrainSet', 'TestSet']\n",
    "\n",
    "var_names = ['p', 'T', 'U']\n",
    "is_vector = [False, False, True]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The snapshots are dependent on two different parameters: the Reynolds and the Richardson number, split into train and test set\n",
    "\\begin{equation*}\n",
    "\\begin{array}{cc}\n",
    "Re_{train} = [15:5:150] & Ri_{train} = [0.2:0.4:5] \\\\\n",
    "Re_{test} = [17.5:10:147.5] & Ri_{test} = [0.4:0.8:44] \n",
    "\\end{array}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fc379e6d7e0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhkAAAIfCAYAAADdQ86WAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABw7klEQVR4nO3deXwTdf4/8Nekd0spR4UWWs7lvg9RFCiosIIH3ooXIB6wgGBZRNejsN9V1FXX3XWLgoqAIrgU/IkcUrGUVlYohXKJHAKCUMCqpIXSNmnfvz+6zfZKmmQmmSTzej4eeQDJfOb9+sykzZtkZqKIiICIiIhIYya9AxAREVFgYpNBREREHsEmg4iIiDyCTQYRERF5BJsMIiIi8gg2GUREROQRbDKIiIjII9hkEBERkUcE6x1ADxUVFThz5gyio6OhKIrecYiIiPyGiKCoqAitWrWCyeT4vQpDNhlnzpxBYmKi3jGIiIj81qlTp5CQkOBwGUM2GdHR0QAqN1Djxo11TlPJYrFg06ZNGDVqFEJCQvSO41VGnbtR5w0Yd+5GnTfAuQfS3AsLC5GYmGh7LXXEkE1G1UckjRs39qkmIzIyEo0bNw6IJ6ErjDp3o84bMO7cjTpvgHMPxLk7c7gBD/wkIiIij2CTQURERB7BJoOIiIg8gk0GEREReQSbDCIiIvIINhlERETkEYY8hZWIiJwnIigtLUVFRYXb67BYLAgJCUFxcXFAncbpDF+fe3BwMEJDQz2zbo+slYiI/F5paSl++uknFBUVoby8XPX6WrZsiaNHj2qQzP/4+twjIiIQFxeHZs2aabpeNhlERFTHxYsXcfToUQQFBaFFixZo1KgRgoKC+H1PAUZEUFZWhoKCAhw/fhwANG002GQQEVEdZ86cQWhoKDp37ozgYL5UBLKoqCg0adIER48exdmzZzVtMnjgJxER1WCxWFBUVISWLVuywTAIRVEQGxuLy5cvo6ysTLP1sskgIqIaLBYLACA8PFznJORNVQd/Wq1WzdbJJoOIiOrF4y+MxRP72y/fB5s7dy7mzZtX476WLVvi7Nmz+gQqLweysoD8fCA+Hhg6FAgKcmm8kpmJ1lu3QomKAkaMcHm82/U1yK7b3H0gu5rtrts+1yA797nB9jmRu8QPpaSkSI8ePSQ/P992O3/+vNPjzWazABCz2aw+TFqaSEKCCPC/W0JC5f2+Pp7Zmd2fsnO7eW38pUuXZOfOnXLp0iXn1k8Bwdn97sprqN82GX369HF7vGZNRlqaiKLU/MEFKu9TlIZ/Aeg5ntmZ3Z+yc7t5dTybDGNik/FfKSkpEhkZKfHx8dKuXTu599575YcffnB6vCZNhtVa938GtX+AExMrl/O18czO7P6UndvN6+PZZNQEwKVb27ZtNc+QlJQkAOT48eOar7uKJ5oMvzwm46qrrsLSpUvRuXNnnDt3Dn/5y19wzTXX4MCBA2jevHmd5UtLS1FaWmr7d2FhIYDKI6irjqJ2lZKZieCffrK/gAhw6hSsGRmQpCSfGs/szO5P2bndvD/e3d+LgWr8+PF17svOzsYPP/yAPn36oG/fvjUei42N9VIyz2jotdGV54dfNhmjR4+2/b1Xr14YPHgwOnbsiCVLliA5ObnO8vPnz69zoCgAbNq0CZGRkW5laL11KwY6sVzehg04femST41ndn3GM7v/1VY73l+zh4SEoGXLlk6MNIYPP/ywzn0TJkzADz/8gNtuuw1z5871eIalS5eiuLgYrVu39nitb775xmEjUVxc7PzKtHqbRW833HCDTJ48ud7HSkpKxGw2226nTp0SAFJQUCBlZWVu3Szp6fbfgqx2s6Sn+9x4Zmd2f8rO7eb98RcuXPD+xyVWq0hGhsjy5ZV/2vsIyEeMHz9eAEhKSoreUTRT9XHJhQsXHL7+FRQUCBDAx2TUVlJSIq1bt5Z58+Y5tbymx2TUd0AV4PxnpXqMZ3Zm96fs3G5eH+/1YzLUnj2jA3tNRkZGhgCQ8ePHS35+vkyaNElat24tQUFB8re//U1ERM6cOSOvvvqqDBs2TFq1aiUhISHSsmVLuf3222XHjh311rN3TAZQeQyI1WqVV199VTp16iShoaGSkJAgTz/9tJSUlDg9Jx74+V+zZs2SLVu2yLFjx+Tbb7+Vm2++WaKjo+XEiRNOjdf87JLaP8CuHvWtx3hmZ3Z/ys7t5tXxXm0y1J49o5OGmowxY8ZIQkKCxMXFyV133SU333yzvPvuuyIismDBAgEgv/vd7+T3v/+93H333dKvXz8BICEhIfLll1/WqddQk3HvvfdKVFSUjBgxQm6++WaJiYkRAPLAAw84PSc2Gf917733Snx8vISEhEirVq3kjjvukAMHDjg93uPXyUhMVHf+urfGMzuz+1N2bjevjfdak6H27BkdNdRkAJDbb79dLl++XGfs3r17Zc+ePXXu37hxo4SGhkrHjh2loqKixmOOmgwA0q1btxqPHTt2TJo2bSoA5OjRo07NiU2GRjRtMkTUf5ZotYolPV1ykpMrP0N1Y7zb9TXIrtvcfSC7mu2u2z5XO577POD3udeajIwM+w1G9VtGhmdzuKGhJiMsLEx++uknl9f7wAMPCADZu3dvjfsbajK++uqrOuuaPn26AJDFixc7VZunsPqqoCBg+HBV4yUpCacvXUKfpCTXL/Wrpr4G2XWbuw9kV7Pdddvnasdzn7s91m/3uafk52u7nA/p37+/wzNBSktLsXHjRuzYsQM///yz7ZtP9+3bBwA4cuQIevXq5VStkJAQDK9n33bu3BkAkK/j9mOTQURE+oiP13Y5H9KmTRu7j+3btw+33norTpw4YXeZoqIip2vFx8cjqJ6mtVGjRgBQ4zpR3sZvYSUiIn0MHQokJAD2vv1TUYDExMrl/Ex4eHi994sI7rnnHpw4cQKTJ09GXl4eCgsLUVFRARHBs88+a1vOWb78bbl8J4OIiPQRFAT8/e/AXXdVNhTVX1irXjjfeiugvi32+++/x/fff4+BAwdiwYIFdR4/duyYDqk8h+9kEBGRfu64A1i1Cqh9/EJCQuX9d9yhTy4P+e233wAACQkJ9T6Wnp7u7UgexXcyiIhIX3fcAYwdC2RlVR7kGR9f+RFJAL2DUeV3v/sdTCYTvv76axw5cgSdOnUCAJSUlGDy5Mn49ddfdU6oLTYZRESkP188+8UDWrRogUmTJmHRokXo06cPrrvuOkRERCArKwvl5eWYMGFCvd+V4q/4cQkREZEXLViwAG+88Qbat2+PzZs3IysrCzfccAN27tyJtm3b6h1PU3wng4iIyEUffvhhve84DB8+vMEzQ4KCgpCcnFzvt4bPnTu33m913bJlS73rclRrwoQJmDBhgsMsnsZ3MoiIiMgj2GQQERGRR7DJICIiIo9gk0FEREQewSaDiIiIPIJNBhEREXkEmwwiIiLyCDYZRERE5BFsMoiIiMgj2GQQERGRR7DJICIiIo9gk0FEREQewSaDiIiIPIJNBhERkQOKorh0a9eund6RfQa/6p2IiMiB8ePH17kvOzsbP/zwA/r06YO+ffvWeCw2NtYjOU6cOIH27dsjKSnJ7le/+xo2GVooLweysoD8fCA+Hhg6FAgKcmm8kpmJ1lu3QomKAkaMcHm82/U1yK7b3H0gu5rtrts+1yA797nB9rnBffjhh3XumzBhAn744QfcdtttmDt3rtcz+Q0xILPZLADEbDarX1lamkhCggjwv1tCQuX9vj6e2Zndn7Jzu3lt/KVLl2Tnzp1y6dIl59ZvQOPHjxcAkpKS4rWax48fFwCSlJTkkfU7u99deQ1lk6FGWpqIotT8wQUq71OUhn8B6Dme2Zndn7Jzu3l1vDebjB8v/Ci5Z3Lt3n688KPHM7jDUZNRWloqb731lgwcOFAaNWokkZGRcuWVV8p7770nFRUVdZY/efKk/OEPf5DOnTtLRESENG3aVLp37y6PP/64fP/99yIikpKSIgDqvY0fP16TObHJ0IgmTYbVWvd/BrV/gBMTK5fztfHMzuz+lJ3bzevjvdVk/HjhRwn/S7hgLuzewv8S7pONhr0m4+LFizJ06FABILGxsXLjjTfKmDFjpGnTpgJAnnjiiRrLnzp1SmJjYwWA9O7dW+655x659dZbpU+fPqIoiixevFhERNasWSN33nmnAJCWLVvK+PHjbbdFixZpMic2GRrRpMnIyLD/g1v9lpHhe+OZndn9KTu3m9fHe6vJyD2T67DBqLrlnsn1aA532GsypkyZIgDkoYcekqKiItv958+fl6uuukoAyBdffGG7v+odijfeeKNOjRMnTsjRo0dt//bHj0t4Cqu78vPVLafneGbXZzyz+19tteP9OTu57Pz583jvvffQvn17LFq0CI0aNbI9dsUVV+Ddd98FANufVWMA4LrrrquzvrZt26Jjx44eTu1ZbDLcFR+vbjk9xzO7PuOZ3f9qqx3vz9nJZZmZmbBYLLjxxhsRFhZW5/E+ffogOjoaOTk5tvsGDBgAAJg6dSoyMjJgtVq9ltcrtHqbxZ9oekxGfQdUAc5/VqrHeGZndn/Kzu3m9fH8uKRh9X1c8tprrwlQ/8GZ1W/BwcG2MVarVe655x7bY5GRkTJs2DCZP3++nDt3rkZNflxiJEFBwN//Xvl3Ran5WNW/33rL/nnoeo5ndmb3p+zcbvqNJ5eUl5cDAPr164fx48fbvT3wwAO2MUFBQVi5ciV27dqFlJQUDBw4EN9++y2effZZdOrUCd9++61e09GGVh2QP/H4dTISE9Wdv+6t8czO7P6UndvNa+P5TkbD6nsnY9myZQJAnnrqKVXrNpvNkpycLABk0KBBtvv98Z0MRUREzyZHD4WFhYiJiYHZbEbjxo3Vr1CDK/FZMzKQt2ED+o4ejWCDXfHT7bn7QHY12123fa5Bdu7zwN7nxcXFOHjwILp164bIyEjna7hoV/4uDFg4oMHlch/PRf/4/h7L4Y4JEyZgyZIlSElJsV3x8/Tp02jbti06dOiAgwcPIkjFO0SlpaWIiIhAeHg4iouLAQBnzpxB69atce211yI7O1uLadTg7H535TWUlxXXQlAQMHy4qvGSlITTly6hT1KS629dqqmvQXbd5u4D2dVsd932udrx3Oduj/Xbfe4hsZGxCA8OR4m1xO4y4cHhiI30zHeBaK1169aYMGEC3n//fTz00EP4xz/+Ued7TLZt24YLFy5gzJgxAIBly5ahX79+6NmzZ43lNm7cCBFBmzZtbPfFxsYiJCQEP/zwA8rLy1U1Md7CJoOIiHTRJqYNDk07hILiArvLxEbGok1MG7uP+5p//OMfOHbsGD755BN88cUX6Nu3L1q1aoWzZ8/i6NGjOH36NGbMmGFrMtLS0vDwww+jY8eO6NWrFyIiInDixAl8++23CAoKwssvv2xbd2hoKG688UasXbsWffr0Qf/+/REaGoprr70WEydO1GvKDrHJICIi3bSJaeNXTURDIiMjsWnTJixZsgTLli3D3r17sX37drRo0QIdO3bEjBkzMG7cONvyycnJSEhIwDfffIOsrCxcunQJrVu3xrhx4/DHP/4R/fr1q7H+9957D3/84x+Rnp6O5cuXo7y8HFarlU0GERFRoPjwww/r/XZWAAgODsakSZMwadKkBtczbNgwDBs2zOm6LVq0wNKlS51eXm88hZWIiIg8gk0GEREReQQ/LvGgk+aTuh3QpGdtvesbtbbe9Y1aW+/6es+dyBE2GR5y0nwSXd7u0uCpWYemHdL8F4CetfWub9Taetc3am1X68dHavsdIXrPnagh/LjEQwqKCxz+4ANAibXE4f9A/LG23vWNWlvv+katrXd9vedO1BA2GUREROQRbDKIiIjII9hkEBFRvQz41VaG5on9zSaDiIhqCAkJAQCUlDg+3oMCS1lZGYDKi4lphU0GERHVEBISgujoaJw7dw5Wq1XvOOQFIoKCggJEREQgNDRUs/XyFFYiIqqjVatWOHr0KA4ePIjmzZujUaNGCAoKgqIoekcjDYkIysrKUFBQgMLCQrRv317T9TvVZLz66qsYOHAgBgwYgCZNmmgaIFDp+RXGen99slHnzu1uvNp61/dk7UaNGqFbt244ffo0zp8/j/z8fDVRycdFRESgffv2aNasmabrdarJePbZZ6EoChRFqfHW2b///W/06dMHnTt31jRUINDzK4z1/vpko86d2914tV2tb7FYdKvtjrCwMHTo0AEigtLSUlRUVNR4PL8oHxdKLtgd3yS8CeKjKy9AZrFY8M033+Daa6+1HfOhliv1tabn3LWed3BwsKYfkdRYt7MB6vtc7t5774WiKIiIiECvXr3Qt29f9OnTB3379kWvXr0QFRWleWB/oudXGOv99clGnTu3u/Fq613fG7UVRUF4eHid+ztGdnR6HRaLBRaLBZGRkZo1Ga7U15qec9dz3q5yqsm4ePEi8vLykJubW+cxEUFxcTG2b9+OHTt22O5XFAUdO3a0NR1Vf7Zu3Vq79EREROSznGoyQkNDMWjQIAwaNKjG/UeOHEFeXh727Nlj+/PUqVMAKpuPI0eO4OjRo0hLS7ONadasma3heP311zWcChEREfkSVWeXdOzYER07dsSdd95pu++3336r0XTk5eXh4MGDtvNvf/nlF3z99dfIyMhgk0FERBTAND+FtWnTphg+fDiGDx9uu89qteK7776zNR15eXnYu3ev1qWJiIjIh3jlOhnBwcHo3bs3evfujYceesgbJYmIiEhnvOInEREReQSv+KmF8nIgKwvIzwfi44GhQ4GgIJfGK5mZaL11K5SoKGDECJfHu11fg+y6zd0HsqvZ7rrtcw2yc59zn/N3nPPj/fb5rgVx0ZkzZ+SNN96QRx55RB5++GFJSUmRXbt2ubSOn3/+WT788ENXS2vGbDYLADGbzepXlpYmkpAgAvzvlpBQeb+vj2d2Zven7NxuzM7s3svugCuvoS41GStXrpSoqCgxmUx1bvfee69cvHjR7th9+/bJyy+/LNdcc40EBwdLUFCQK6U1pVmTkZYmoig1dyJQeZ+iNLwz9RzP7MzuT9m53Zid2b2XvQEeaTIOHjwo4eHhoihKvTeTySRJSUlSXl5uG7Nr1y5JTk6Wdu3a1WhIqpbXiyZNhtVat0usvTMTEyuX87XxzM7s/pSd243Zmd172Z3gkSZj8uTJtuZg8ODBsnTpUsnMzJRPP/1U7r//fgkKChKTySTvvPOOnDhxQpKSkuo0FlXjBwwYIC+++KLbE1RLkyYjI8P+Tqx+y8jwvfHMzuz+lJ3bjdmZ3XvZneDKa6jTB35u2bIFiqLgmmuuwdatW2t83e/dd99tu73//vv4+9//ju+//972eKNGjXDDDTfgpptuwpgxYxAf75kvrPEqZ7+R0N5yeo5ndn3GM7v/1VY7ntn1Gc/s7o/XmNNNRtXlwmfOnFmjwahy2223Yfr06Xjrrbdsj/fo0QPPP/88xo4dW++X6/g1Zxsle8vpOZ7Z9RnP7P5XW+14ZtdnPLO7P15rzr49UvVRx+7du+0us2/fvhofiZSUlDi7eq/S9JiM+g6uAZz/3EyP8czO7P6UnduN2Znde9md4JFjMqqahwMHDthdprS01Lbcxx9/7OyqvU7zs0tq70xXjwDWYzyzM7s/Zed2Y3Zm9172BujWZFRfLicnx9lVe53Hr5ORmKjuXGZvjWd2Zven7NxuzM7s3svugCuvoYqIiDMfq5hMJiiKgv3796Nbt24NLrdv3z50795dgw90tFdYWIiYmBiYzWY0btxY/Qo1uCqbNSMDeRs2oO/o0Qg22NXw3J67D2RXs9112+caZOc+5z7n7zjnx/vt890OV15DXW4yYmJi0KtXL/Tu3dv2Z8+ePdGoUaMayxmqydCAxWLB+vXrMWbMGISEhOgdx6uMOnejzhsw7tyNOm+Acw+kubvyGurSd5eICC5cuIDs7GxkZ2fXeKxdu3bo1auX7d+//vqrK6smIiKiAON0k5Gamoq8vDzk5eVh//79KC4urvH48ePHceLECdvpq0lJSYiNjUW/fv3Qt29f25+dO3eu9xRYIiIiCixONxmTJ0+2/V1EcOjQIVvTkZeXhz179uDcuXM1xvz8889IT09Henq67b7IyEj06tUL/fr1w7/+9S8NpkBERES+yK2velcUBV27dkXXrl1x33332e4/d+4cdu/eXaP5OHr0KCoqKmzLXLp0Cd9++y22b9/OJoOIiCiAudVk2NOyZUvceOONuPHGG233Xb58GXv27KnReOzfvx+XL1/WsjQRERH5GE2bjPpERETg6quvxtVXX227T0Rw+PBhT5cmIiIiHZn0KKooCrp06aLJuubPnw9FUTBz5kxN1kdERETa0KXJ0EpOTg4WLlyI3r176x2FiIiIavHbJuPixYt44IEHsGjRIjRt2lTvOERERFSLx4/J8JSpU6fipptuwg033IC//OUvDpctLS1FaWmp7d+FhYUAKq/CZrFYPJrTWVU5fCWPNxl17kadN2DcuRt13gDnXv1Pf+fKPJy+rLgvWbFiBV566SXk5OQgPDwcw4cPR9++ffHWW2/Vu/zcuXMxb968OvcvX74ckZGRHk5LREQUOIqLi3H//fdr+90lvuLUqVMYOHAgNm3ahD59+gBAg01Gfe9kJCYmoqCgwKe+uyQ9PR0jR44MiGvbu8KoczfqvAHjzt2o8wY490Cae2FhIWJjY7X/7hJfkJubi/Pnz2PAgAG2+8rLy7F161a8/fbbKC0tRVCtb5kLCwtDWFhYnXWFhIT43A73xUzeYtS5G3XegHHnbtR5A5x7IMzdlTn4XZNx/fXXY9++fTXumzhxIrp27Yo5c+bUaTCIiIhIH6qajKVLlwIAunTpgquuukqTQA2Jjo5Gz549a9wXFRWF5s2b17mfiIiI9KPqFNYJEyZg4sSJ+PHHH7XKQ0RERAFC1TsZMTExKCwsRKdOnbTK45YtW7boWp+IiIjqUvVORvv27QEAv/32myZhiIiIKHCoajJuv/12iAjWrl2rVR4iIiIKEKqajBkzZqBt27ZYsGABvv76a60yERERUQBQ1WQ0btwY6enp6Nq1K37/+9/j8ccfx5YtW/Drr7/Cz67xRURERBpTdeBn9WtSiAjef/99vP/++06NVRQFVqtVTXkiIiLyYaqajNrvVvDdCyIiIqqiqslISUnRKod/Ky8HsrKA/HwgPh4YOhRw5cqj5eVQMjPReutWKFFRwIgRLo93u74G2XWbuw9kV7PdddvnGmTnPuc+5+8458f77fNdC2JAZrNZAIjZbFa/srQ0kYQEEeB/t4SEyvt9fTyzM7s/Zed2Y3Zm9152B1x5DWWToUZamoii1NyJQOV9itLwztRzPLMzuz9l53Zjdmb3XvYGsMlogCZNhtVat0usvTMTEyuX87XxzM7s/pSd243Zmd172Z3AJqMBmjQZGRn2d2L1W0aG741ndmb3p+zcbszO7N7L7gRXXkM1+6r3iooKbNmyBf/5z39w9uxZFBcX4y9/+Qvi4+Nty5SVlcFqtSIoKAhhYWFaldZHfr665fQcz+z6jGd2/6utdjyz6zOe2d0frzFNmox169bhySefxIkTJ2rcP2vWrBpNxvvvv49p06ahUaNGOHPmDKKiorQor49q83JrOT3HM7s+45nd/2qrHc/s+oxndvfHa83t90v+a9GiRWIymURRFFEURa644gpRFEVMJpMcOHCgxrKlpaXSrFkzMZlMsmzZMrWl3abpMRn1HVwDOP+5mR7jmZ3Z/Sk7txuzM7v3sjvBa8dkHDlyREJDQ8VkMsn1118vBw8eFBGx22SIiDz22GOiKIo89NBDakqrovnZJbV3pqtHAOsxntmZ3Z+yc7sxO7N7L3sDvNZkTJ06VRRFkV69eklpaantfkdNxtKlS21j9OLx62QkJqo7l9lb45md2f0pO7cbszO797I74MprqCIi4u5HLd26dcPhw4exaNEiPPLII7b7TSYTFEXBvn370L179xpjtm3bhiFDhqBx48a4cOGCu6VVKSwsRExMDMxmMxo3bqx+hRpclc2akYG8DRvQd/RoBBvsanhuz90HsqvZ7rrtcw2yc59zn/N3nPPj/fb5bocrr6GqmoxGjRrh8uXLyMnJQf/+/W33O2oy9uzZg379+iE4OBhlZWXullZF8yZDAxaLBevXr8eYMWMQEhKidxyvMurcjTpvwLhzN+q8Ac49kObuymuoqq96VxQFAOBKn/Lzzz8DgM+8uBMREZFnqGoyWrVqBQA4fPiw02MyMzMBAO3atVNTmoiIiHycqiZj2LBhEBEsX77cqeULCgrw7rvvQlEUXHfddWpKExERkY9T1WQ8/vjjAID169dj8eLFDpf96aefMGbMGBQUFCAoKMg2loiIiAKTqibjyiuvxOTJkyEiePTRR3H33Xfj008/tT2+d+9erFy5EpMmTUKXLl2Qm5sLRVEwa9Ys/O53v1MdnoiIiHyX6suK//Of/8SlS5ewbNkyrF69GqtXr7YdEPrAAw/Ylqs6OHTChAl4+eWX1Zb1CyfNJ1FQXGD38djIWLSJaRNwtfWub9Taetc3am296xu1tt71jVrbVaqbjKCgICxZsgS33nor5s+fj127dtW7XPfu3fH888/jvvvuU1vSL5w0n0SXt7ugxFpid5nw4HAcmnZI8yeDnrX1rm/U2nrXN2ptV+vHR2r7fRHc7sabu97b3VWqPi6p7s4778TOnTvx008/4bPPPsPChQuxYMECfPrppzhy5Aj2799vmAYDAAqKCxw+CQCgxFrisBv1x9p61zdqbb3rG7W23vWNWlvv+kat7Q7Nvuq9SqtWrXDrrbdqvVoiIiLyM5o3GQBgtVrx22+/AQCaNm2K4GCPlCEiIiIfptnHJQcOHMD06dPRrVs3hIeHIy4uDnFxcQgPD0e3bt0wffp07N+/X6tyRERE5ONUNxkVFRWYOXMm+vbti9TUVBw6dAgVFRWQym94RUVFBQ4dOoTU1FT069cPTz31FCoqKrTITkRERD5M9ecY9913H9LS0mynqPbo0QODBg1Cy5YtISI4f/48cnJysH//fpSXl+Mf//gHzpw5g5UrV6oOT0RERL5LVZOxfPlyrFq1CoqioE+fPli4cCGuvPLKepfduXMnnnjiCezevRurVq3CihUrDHW2CRERkdGo+rhk0aJFAIDOnTsjOzvbboMBAAMHDsTWrVvRpUsXiAjeffddNaV9XmxkLMKDwx0uEx4cjtjI2ICqrXd9o9bWu75Ra+td36i19a5v1NruUPVOxt69e6EoCubMmYOoqKgGl4+KisKcOXPwyCOPYM+ePWpK+7w2MW1waNohXa7Kpmdtvesbtbbe9Y1a29X6FotFt9pa86ftztr6UdVklJWVAQB69+7t9JiqZbX+YfNFbWLa6Laj9aytd32j1ta7vlFr613fqLX1rm/U2q5S9XFJ27ZtAQBms9npMYWFhTXGEhERUWBS1WTceeedEBGkpaU5PabqQNHbb79dTWkiIiLycaqajOTkZHTo0AHvvvtuja94t2fVqlV499130b59e/zxj39UU5qIiIh8nKomIyYmBl999RX69++PcePG4bbbbsNnn32G06dPw2KxwGq14vTp0/jss89w++23495770X//v2xefNmxMTEaDUHIiIi8kFOHfgZFBTU4DIigrVr12Lt2rUOl9m5cyc6dOgARVFgtVqdT0pERER+xakmo+pqnlos5+y6iIiIyL851WSkpKR4Ood/Ky8HsrKA/HwgPh4YOhRw4t2f6uOVzEy03roVSlQUMGKEy+Pdrq9Bdt3m7gPZ1Wx33fa5Btm5z7nP+TvO+fF++3zXghiQ2WwWAGI2m9WvLC1NJCFBBPjfLSGh8n5fH8/szO5P2bndmJ3ZvZfdAVdeQ9lkqJGWJqIoNXciUHmfojS8M/Ucz+zM7k/Zud2Yndm9l70BbDIaoEmTYbXW7RJr78zExMrlfG08szO7P2XndmN2ZvdediewyWiAJk1GRob9nVj9lpHhe+OZndn9KTu3G7Mzu/eyO8GV11BV311SxWq1Yt26dcjKysKxY8dQVFSE8vJyh2MURcHmzZu1KK+P/Hx1y+k5ntn1Gc/s/ldb7Xhm12c8s7s/XmOqm4wtW7Zg4sSJOHnypO0+EbG7vKIoEBEoiqK2tL7i49Utp+d4ZtdnPLP7X22145ldn/HM7v54rbn9fomI7N69W8LDw8VkMomiKBIRESG9e/eWYcOGyfDhwxu86UXTYzLqO7gGcP5zMz3GMzuz+1N2bjdmZ3bvZXeC147JGDt2rCiKIuHh4ZKamiqXL19Wszqv0fzskto709UjgPUYz+zM7k/Zud2Yndm9l70BXmsymjdvLiaTSf7v//5PzWq8zuPXyUhMVHcus7fGMzuz+1N2bjdmZ3bvZXfAlddQRUTE3Y9aGjVqhMuXL2P79u0YOHCgNp/feEFhYSFiYmJgNpvRuHFj9SvU4Kps1owM5G3YgL6jRyPYYFfDc3vuPpBdzXbXbZ9rkJ37nPucv+OcH++3z3c7XHkNVdVk9OzZEwcPHkR2djYGDx7s7mq8TvMmQwMWiwXr16/HmDFjEBISonccrzLq3I06b8C4czfqvAHOPZDm7sprqKqver/tttsAAFu3blWzGiIiIgpAqpqMGTNmID4+Hq+//jpOnDihUSQiIiIKBKqajCuuuALr169HREQErrrqKrz33nswm81aZSMiIiI/pvpiXL1798bWrVtx1VVX4YknnsDkyZMRGxuLyMhIh+MURcEPP/ygtjwRERH5KNVNRlpaGiZNmoSioiJI5SmxOH/+fIPj/P6Kn0REROSQqibjP//5D+677z7b95S0bdsWvXv3RpMmTWAyqfokhoiIiPycqibjL3/5C8rLyxETE4Ply5dj9OjRWuUiIiIiP6fq7Ybc3FwoioJ58+axwSAiIqIaVDUZly5dAgAMGTJEkzBEREQUOFQ1Ge3btwcAFBcXaxKGiIiIAoeqJuOOO+6AiODLL7/UKg8REREFCFVNxqxZs9CpUye89dZb2Llzp1aZiIiIKACoajKio6OxefNm9OzZE8OGDcNzzz2HvXv3oqSkRKt8RERE5KdUncIaVO0rY0UEr7zyCl555RWnxiqKAqvVqqY8ERER+TBVTUbtb4lX8a3xREREFGBUNRkpKSla5SAiIqIAwyaDiIiIPIJfMEJEREQewSaDiIiIPIJNBhEREXmEqmMyTp48qap4mzZtVI0nIiIi36Wqyaj67hJ38DoZREREgU3T62QQERERVVHVZCxevLjBZS5duoRDhw4hLS0NZ86cwTXXXIPHHntMTVnfU14OZGUB+flAfDwwdChQ7WqozoxXMjPReutWKFFRwIgRLo93u74G2XWbuw9kV7PdddvnGmTnPuc+5+8458f77fNdC+IlZWVlMnnyZDGZTDJ79mxV60pNTZVevXpJdHS0REdHy9VXXy3r1693erzZbBYAYjabVeUQEZG0NJGEBBHgf7eEhMr7fX08szO7P2XndmN2ZvdedgdceQ31WpNRZfjw4WIymWTjxo1ur+Pzzz+XdevWyaFDh+TQoUPypz/9SUJCQmT//v1OjdesyUhLE1GUmjsRqLxPURremXqOZ3Zm96fs3G7Mzuzey94An24yVq5cKYqiyE033aTpeps2bSrvvfeeU8tq0mRYrXW7xNo7MzGxcjlfG8/szO5P2bndmJ3ZvZfdCT7dZOzatUsURZGWLVtqsj6r1SqffPKJhIaGyoEDB+pdpqSkRMxms+126tQpASAFBQVSVlbm1s2Snm5/J1a7WdLTfW48szO7P2XndmN2ZvdedmduBQUF4myToerAT3eYzeYaf7pr3759GDx4MEpKStCoUSOsWbMG3bt3r3fZ+fPnY968eXXu37RpEyIjI92q33rrVgx0Yrm8DRtw+tIlnxrP7PqMZ3b/q612PLPrM57Z3R/vjOLiYucX1uTtBBdMmDBBFEWR9u3bq1pPaWmpHDlyRHJycuSZZ56R2NhYvpPhJ50yszO7v9RmdmY3WnZnbq68kwFVr/QuOHz4sDzxxBOiKIqYTCaZOnWqpuu//vrr5fHHH3dqWU2Pyajv4BrA+c/N9BjP7MzuT9m53Zid2b2X3QleOyajffv2Dd7atm0rMTExYjKZxGQyiaIoEhcXJ2fOnFFTuo7rrrtOxo8f79Symp9dUntnunoEsB7jmZ3Z/Sk7txuzM7v3sjfAa02Goigu3wYPHizff/+9mrLy7LPPytatW+X48eOyd+9e+dOf/iQmk0k2bdrk1HiPXycjMVHduczeGs/szO5P2bndmJ3ZvZfdAVdeQxUREbeO/AAwceLEBpcxmUyIjo5G+/btkZSUhL59+7pbzmbSpEnYvHkz8vPzERMTg969e2POnDkYOXKkU+MLCwsRExMDs9mMxo0bq86jxVXZrBkZyNuwAX1Hj0awwa6G5/bcfSC7mu2u2z7XIDv3Ofc5f8c5P95vn+92uPIaqqrJ8FeaNxkasFgsWL9+PcaMGYOQkBC943iVUedu1HkDxp27UecNcO6BNHdXXkNNXspEREREBsMmg4iIiDyCTQYRERF5BJsMIiIi8ginLysepPF30CuKAqvVquk6iYiIyHc43WQY8CQU1U6aT6KguMDu47GRsWgT0ybgautd36i19a5v1Np61zdqbb3rG7W2q5xuMsaPH6+qkIhg/fr1+OWXXwzRsJw0n0SXt7ugxFpid5nw4HAcmnZI8yeDnrX1rm/U2nrXN2ptV+vHR8brVtvI2z2Q5q73dneV003G4sWL3S7y2WefISUlBb/88ovtvsTERLfX5w8KigscPgkAoMRagoLiAs2fCHrW1ru+UWvrXd+otV2tr3WTwe1uvLnrvd1d5dEDP9evX4+BAwfizjvvxP79+yEiiIuLwz//+U8cPnzYk6WJiIhIZ06/k+GK9PR0vPjii9ixYweAyo9KWrRogTlz5mDKlCkIDw/3RFkiIiLyIZo2GZmZmXjhhRfwzTffAKhsLpo3b47Zs2dj2rRpiIyM1LIcERER+TBNmoxt27bhxRdfREZGBoDK5qJJkyZITk7GzJkz0ahRIy3KEBERkR9R1WTk5OTgxRdfxKZNmwBUNhfR0dGYOXMmkpOTERMTo0lIIiIi8j9uNRl5eXlISUnBF198AaCyuYiKisK0adMwe/ZsNGvWTNOQRERE5H9cOrvkwIEDuOuuuzBgwAB88cUXEBGEh4cjOTkZx48fx/z589lg/FdsZCzCgx0f4BoeHI7YyNiAqq13faPW1ru+UWvrXd+otfWub9Ta7nD6nYz7778fn376KUQEIoKwsDA88cQTeOaZZxAXF+fJjH6pTUwbHJp2SJersulZW+/6Rq2td32j1na1vsVi0a221vxpu7O2fpxuMlasWGH7+xVXXIEZM2YgISHBdjyGOx5++GG3x/qDNjFtdNvRetbWu75Ra+td36i19a5v1Np61zdqbVe5dEyGoigAgIKCArzwwguqCiuKEvBNBhERkZG51GQY4TtHiIiISBtONxlV18AgIiIicobTTUZSUpIncxAREVGA8egXpBEREZFxsckgIiIij/DIt7AaTnk5kJUF5OcD8fHA0KFAUJBL45XMTLTeuhVKVBQwYoTL492ur0F23ebuA9nVbHfd9rkG2bnPuc/5O8758X77fNeCGJDZbBYAYjab1a8sLU0kIUEE+N8tIaHyfl8fz+zM7k/Zud2Yndm9l90BV15D2WSokZYmoig1dyJQeZ+iNLwz9RzP7MzuT9m53Zid2b2XvQFsMhqgSZNhtdbtEmvvzMTEyuV8bTyzM7s/Zed2Y3Zm9152J7DJaIAmTUZGhv2dWP2WkeF745md2f0pO7cbszO797I7wZXXUJ5d4q78fHXL6Tme2fUZz+z+V1vteGbXZzyzuz9eY2wy3BUfr245Pcczuz7jmd3/aqsdz+z6jGd298drze33S/yYpsdk1HdwDeD852Z6jGd2Zven7NxuzM7s3svuBB6T0QDNzy6pvTNdPQJYj/HMzuz+lJ3bjdmZ3XvZG6Bbk1FYWCh79uyR7OxsyczMbPCmF49fJyMxUd25zN4az+zM7k/Zud2Yndm9l90BV15DFRERtR+5LFq0CKmpqdi7d6/TYxRFgdVqVVvaLYWFhYiJiYHZbEbjxo3Vr1CDq7JZMzKQt2ED+o4ejWCDXQ3P7bn7QHY12123fa5Bdu5z7nP+jnN+vN8+3+1w5TVUVZNRXl6OO++8E2vXrgUAuLIqRVFQXl7ubmlVNG8yNGCxWLB+/XqMGTMGISEhesfxKqPO3ajzBow7d6POG+DcA2nurryGqvruknfeeQeff/45AKBly5aYOHEiBgwYgGbNmsFk4okrRERERqaqyVi6dCkAoHv37sjKykLTpk01CUVERET+T9XbDQcPHoSiKHjhhRfYYBAREVENmnym0aVLFy1WQ0RERAFEVZPRqVMnAMCvv/6qSRgiIiIKHKqajPvuuw8igi+++EKrPERERBQgVDUZTz75JHr37o0FCxYgKytLq0xEREQUAFQ1GWFhYdi0aRMGDBiAkSNH4umnn0ZeXh5KSkq0ykdERER+StUprEHVrhwmInjjjTfwxhtvODVWzyt+EhERkeepajJqX+FTgyuUExERUYBQ1WSkpKRolYOIiIgCDJsMIiIi8gh+wQgRERF5BJsMIiIi8ghVH5fU59y5c9i/f7/tKqDNmjVDz5490bJlS61LERERkQ/TpMkQESxcuBBvv/02vvvuu3qX6d69O6ZPn47HHnsMiqJoUZaIiIh8mOqPS3777TcMHToUf/jDH/Ddd99BROq9fffdd5gyZQqGDRuGCxcuaBCdiIiIfJnq62SMHTsW27ZtAwA0b94c99xzD6666irExcVBRHDu3Dns2LEDn376KQoKCrBt2zaMHTsWmZmZmkyAiIiIfJOqJmP58uXIzs6Goii4//77kZqaiujo6DrLPfzww3jllVcwdepULFu2DNnZ2fjkk08wbtw4NeWJiIjIh6n6uGT58uUAgKSkJCxbtqzeBqNKo0aNsGTJEiQlJUFE8NFHH6kpTURERD5OVZOxa9cuKIqCadOmOT1m+vTpAIDdu3erKU1EREQ+TlWTUXWaavv27Z0eU7Vs1VgiIiIKTKqajJiYGADAmTNnnB5TtWzjxo3VlCYiIiIfp+rAz549eyIzMxOLFy/GTTfd5NSYDz74wDY2YJSXA1lZQH4+EB8PDB0KBAW5NF7JzETrrVuhREUBI0a4PN7t+hpk123uPpBdzXbXbZ9rkJ37nPucv+OcH++3z3ctiAr/+te/RFEUMZlMkpKSIhUVFXaXraiokJSUFNvyqampakqrYjabBYCYzWb1K0tLE0lIEAH+d0tIqLzf18czO7P7U3ZuN2Zndu9ld8CV11BVTUZZWZl069bN1jj07NlTXn/9dcnKypLDhw/LkSNHJCsrS15//XXp1auXmEwmURRFunfvLhaLRU1pVTRrMtLSRBSl5k4EKu9TlIZ3pp7jmZ3Z/Sk7txuzM7v3sjfAa02GiMiJEyekQ4cOtkbD0U1RFOnYsaP8+OOPasuqokmTYbXW7RJr78zExMrlfG08szO7P2XndmN2Zvdedid4tckQEbl48aL88Y9/lKZNm4qiKPXemjZtKrNnz5aioiItSqqiSZORkWF/J1a/ZWT43nhmZ3Z/ys7txuzM7r3sTnDlNVSTL0iLiorCX//6V7z00kvIzc2t91tYBwwYgNDQUC3K+Yb8fHXL6Tme2fUZz+z+V1vteGbXZzyzuz9eY5p+1XtoaCgGDx6MwYMHa7la3xQfr245Pcczuz7jmd3/aqsdz+z6jGd298drze33S/yYpsdk1HdwDeD852Z6jGd2Zven7NxuzM7s3svuBK8fk9GQkpISOXv2rJSXl3ujXIM0P7uk9s509QhgPcYzO7P7U3ZuN2Zndu9lb4DXmoyioiJZt26drFu3rt4DOn/++We54447JDQ0VEwmkzRu3FhmzZolpaWlasqq5vHrZCQmqjuX2VvjmZ3Z/Sk7txuzM7v3sjvgymuoIiLi7kctS5YswcSJE9GmTRscO3YMJtP/rlJeUVGBq666Crt27UL1Eoqi4M4778Snn37qblnVCgsLERMTA7PZrM3lzTW4Kps1IwN5Gzag7+jRCDbY1fDcnrsPZFez3XXb5xpk5z7nPufvOOfH++3z3Q6XXkPVdDPjxo0TRVEkOTm5zmPLly+XqmtnDBgwQJKTk2XAgAG2+zZs2KCmtCqavpOhkbKyMvnss8+krKxM7yheZ9S5G3XeIsadu1HnLcK5B9LcvXYK6/79+6EoSr1nkyxbtgwAMGDAAGzbtg3BwcGwWCwYOnQocnJysHTpUtx4441qyhMREZEPU/UtrD///DMAoG3btjXut1gsyMzMhKIo+MMf/oDg4MpeJiQkBJMnT4aIYPv27WpKExERkY9T1WRUXXArJCSkxv07d+7E5cuXAQCjR4+u8Vjnzp0BAGfPnlVTmoiIiHycqiYjIiICAHD+/Pka92dmZgIAOnbsiJYtW9Y7hoiIiAKbqiajY8eOAIAtW7bUuH/NmjVQFAVJSUl1xlR9xNKiRQs1pYmIiMjHqTrwc+TIkdi9ezdSU1MxdOhQDB06FIsXL0ZOTg4URcEtt9xSZ8zevXsBAK1atVJT2i+cNJ9EQXGB3cdjI2PRJqZNwNXWu75Ra+td36i19a5v1Np61zdqbVepajJmzJiBd955B0VFRbj55ptrPNatW7d6m4x169bZPSPFWfPnz8fq1avx/fffIyIiAtdccw1effVVdOnSxe11au2k+SS6vN0FJdYSu8uEB4fj0LRDmj8Z9Kytd32j1ta7vlFru1o/PlLb74vgdjfe3PXe7q5S9XFJfHw81q5di7i4OEjl1UMhIujQoQNWrVoFRVFqLP/DDz8gKysLQOW7IO7KzMzE1KlT8e233yI9PR1WqxWjRo3CpUuX1ExHUwXFBQ6fBABQYi1x2I36Y2296xu1tt71jVpb7/pGra13faPWdofqb2EdOnQojh8/jm+++QZnz55FfHw8hgwZYjtttbr8/Hy88MILAFDv8RrO2rhxY41/L168GC1atEBubi6GDRvm9nqJiIhIO5p81XtoaChGjBjR4HJDhgzBkCFDtChZg9lsBgA0a9as3sdLS0tRWlpq+3dhYSGAyut5WCwWzfMAgNVqdXq56jm0yONqba0Zde56ztud+loy6j53tz73uffra8nIz3fAtXlo0mToSUSQnJyMIUOGoGfPnvUuM3/+fMybN6/O/Zs2bUJkZKRHcv1Q/INTy2VnZyM/Mt/27/T0dN1qa8Woc9dz3mrq61nb3/e5u/W5z/Wrr2ftQHi+A0BxcbHTyzrVZJw8edLtMI60aaP+oJRp06Zh7969yM7OtrvMs88+i+TkZNu/CwsLkZiYiFGjRmnzBWn12H12N3C44eWGDBmCfnH9YLFYkJ6ejpEjR9a5uJmna2vNqHPXc97u1NeSUfe5q/V7Nu/Jfc65e7W2J1R9GuAMp5qM9u3bux3GHkVRnH7bx57p06fj888/x9atW5GQkGB3ubCwMISFhdW5PyQkRJMne33qOybF3nLVM2iRyd3aWjHq3PWct5r6WjDqPne3Pvc55+7t2lpyZb1OpRX3vw3eI0QE06dPx5o1a7BlyxaPNEFERESkjlNNxuLFix0+npqaipycHISEhGDUqFEYNGgQWrZsCRHB+fPnkZOTg02bNsFiseDKK6/ElClTVIWeOnUqli9fjv/3//4foqOjbd+DEhMT4zOXLY+NjEV4cHiD5zLHRsYGVG296xu1tt71jVpb7/pGra13faPWdodTTcb48ePtPvboo49i586dGDVqFN5//320bt263uVOnz6Nxx57DF9++SV69eqFRYsWuZcYwIIFCwAAw4cPr3H/4sWLMWHCBLfXq6U2MW1waNohXa7Kpmdtvesbtbbe9Y1a29X6Wh/tz+1uvLnrvd1dpersklWrVuGDDz7AlVdeiXXr1iEoKMjusq1bt8batWsxePBgfPDBBxg5ciTuuecet+r62sc39rSJaaPbjtaztt71jVpb7/pGra13faPW1ru+UWu7StUVP999910oioLk5GSHDUaVoKAgzJo1CyKChQsXqilNREREPk5Vk1H1ZWedO3d2ekzVsvv27VNTmoiIiHycqiajqKgIAHD+/Hmnx1QtWzWWiIiIApOqJqNt27YAgKVLlzo9pmpZLS7ERURERL5LVZMxduxYiAhWrFiB1157rcHlX3/9dXzyySdQFAW33367mtJERETk41SdXfLMM89g6dKlOHfuHJ599ll88sknGD9+PK688kq0aNECiqLg3LlzyMnJwbJly5CXlwcAiIuLw5w5c7TIT0RERD5KVZPRpEkTfPXVV/j973+P06dPY+/evZg1a5bd5UUECQkJ2LhxI5o0aaKmtG8pLweysoD8fCA+Hhg6FHDibJvq45XMTLTeuhVKVBQwYoTL492ur0F23ebuA9nVbHfd9rkG2bnPuc/5O8758X77fNeCaMBsNktycrI0a9ZMFEWp99asWTOZNWuWmM1mLUqqzgtAmyxpaSIJCSLA/24JCZX3+/p4Zmd2f8rO7cbszO697A648hqqSZNRpbS0VLZt2ybvvvuuvPLKKzJ//nx59913Zdu2bVJaWqplKVU0azLS0kQUpeZOBCrvU5SGd6ae45md2f0pO7cbszO797I3QLcmw19o0mRYrXW7xNo7MzGxcjlfG8/szO5P2bndmJ3ZvZfdCV5rMpYsWSJLliyRb7/9Vs1qvE6TJiMjw/5OrH7LyPC98czO7P6UnduN2Znde9md4MprqKpTWCdMmICJEyfixx9/VHtoiP/Jz1e3nJ7jmV2f8czuf7XVjmd2fcYzu/vjNaaqyYiJiQEAdOrUSZMwfiU+Xt1yeo5ndn3GM7v/1VY7ntn1Gc/s7o/Xmtvvl4hIv379xGQyyebNm9Wsxus0PSajvoNrAOc/N9NjPLMzuz9l53Zjdmb3XnYneO2YjD//+c+iKIrMnDlTzWq8TvOzS2rvTFePANZjPLMzuz9l53Zjdmb3XvYGeK3JMJvN0q5dOwkLC/OrdzM8fp2MxER15zJ7azyzM7s/Zed2Y3Zm9152B1x5DVVERNR83HL06FHcddddOHDgACZOnIj7778fvXv3RtOmTaEoivrPczygsLAQMTExMJvNaNy4sfoVanBVNmtGBvI2bEDf0aMRbLCr4bk9dx/Irma767bPNcjOfc59zt9xzo/32+e7Ha68hqpqMoKqhRURl5oKRVFgtVrdLa2K5k2GBiwWC9avX48xY8YgJCRE7zheZdS5G3XegHHnbtR5A5x7IM3dlddQVd9dUrs/UfmmCBEREQUQVU1GSkqKVjmIiIgowLDJICIiIo9QdTEuIiIiInvYZBAREZFHsMkgIiIij1B1TEZ1IoK8vDzs2bMHBQUFuHz5coNnm7z44otalSciIiIfo0mTsWTJEsybN8/lb2Nlk0FERBS4VDcZzz33HF555RWnrpGhKAqvpUFERGQQqo7J2L59O+bPnw8AGDlyJPLy8rBr1y4AlQ1FeXk5CgoKsHHjRowdOxYigiFDhiA/Px8VFRXq0xMREZHPUtVkLFiwAADQtm1brFu3Dr17965xyVRFUdCsWTOMGjUKa9aswb/+9S9kZ2fjxhtvRFlZmbrkRERE5NNUNRnbtm2Doih48sknERzc8CcvU6ZMwZ133om9e/ciNTVVTWkiIiLycaqajPz8fABAjx49/rdC0/9WabFY6ox56KGHICJYuXKlmtJERETk41Q1GVVNRIsWLWz3NWrUyPb3n3/+uc6YxMREAJVfEU9ERESBS1WTccUVVwCo/NrXKi1btrR9BfzBgwfrjKl696OoqEhNaSIiIvJxqpqMqo9Jvv/+e9t9oaGhtvvr+0jk448/BgC0atVKTWkiIiLycaqajKFDh0JEkJGRUeP+e++9FyKCDz74AC+++CIOHDiAnJwcTJs2DZ988gkURcHo0aNVBSciIiLfpqrJuO222wAAX3zxRY2PTGbMmIF27dqhoqICL730Enr37o2rr77adspr06ZN8eyzz6opTURERD5O9cclGRkZWLNmDaxWq+3+yMhIZGRk4Nprr4WI1Lj17NkTmzdvRkJCgurwRERE5LtUX1Y8KSmp3vvbtm2LrKwsHDp0CAcOHIDVakWnTp3Qr18/tSWJiIjID2j2Laz2dOnSBV26dPF0GSIiIvIxqj4uISIiIrJH03cyioqKcPz4cRQVFaG8vLzB5YcNG6ZleSIiIvIhmjQZixYtQmpqKvbu3ev0GEVRahws6tfKy4GsLCA/H4iPB4YOBf57QTJnxyuZmWi9dSuUqChgxAiXx7tdX4Psus3dB7Kr2e667XMNsnOfc5/zd5zz4/32+a4FUcFqtcrYsWPFZDKJyWQSRVGcvplMJjWlVTGbzQJAzGaz+pWlpYkkJIgA/7slJFTe7+vjmZ3Z/Sk7txuzM7v3sjvgymuoqibj7bfftjUNcXFx8uyzz8qqVavk66+/li1btjR404tmTUZamoii1NyJQOV9itLwztRzPLMzuz9l53Zjdmb3XvYGeK3JGDRokCiKIj169JBff/1Vzaq8SpMmw2qt2yXW3pmJiZXL+dp4Zmd2f8rO7cbszO697E7wWpMRHR0tJpNJVqxYoWY1XqdJk5GRYX8nVr9lZPjeeGZndn/Kzu3G7MzuvexOcOU1VJNTWA15HYz/fpus28vpOZ7Z9RnP7P5XW+14ZtdnPLO7P15jqpqMTp06AQB+/fVXTcL4lfh4dcvpOZ7Z9RnP7P5XW+14ZtdnPLO7P15rbr9fIiKvvfaaKIoiTz31lJrVeJ2mx2TUd3AN4PznZnqMZ3Zm96fs3G7Mzuzey+4Erx2TUVJSIn369JHw8HDZunWrmlV5leZnl9Tema4eAazHeGZndn/Kzu3G7MzuvewN8FqTISJy7tw5ufbaayUsLExmz54tu3fvlsuXL6tdrUd5/DoZiYnqzmX21nhmZ3Z/ys7txuzM7r3sDrjyGqqIiDT0kUqQE1cIExEoiuL0xzR6XvGzsLAQMTExMJvNaNy4sfoVanBVNmtGBvI2bEDf0aMRbLCr4bk9dx/Irma767bPNcjOfc59zt9xzo/32+e7Ha68hjrVZJhM2n+PmqIoTn2/iSdo3mRowGKxYP369RgzZgxCQkL0juNVRp27UecNGHfuRp03wLkH0txdeQ116rtLUlJSNAlGRERExsEmg4iIiDxC+89BiIiIiMAmg4iIiDzEqY9L7CkqKsLf/vY3AMDjjz+OuLg4h8vn5+dj0aJFAIDZs2cjIiJCTXkiIiLyYaqajM8++wxz585Fp06d8OKLLza4fFxcHD7++GMcPXoUXbt2xT333KOmvM87aT6JguICu4/HRsaiTUybgKutd32j1ta7vlFr613fqLX1rm/U2q5S1WSsXr0aiqI43SwoioL77rsP//d//4d///vfAd1knDSfRJe3u6DEWmJ3mfDgcByadkjzJ4OetfWub9Taetc3am1X68dHavt9Edzuxpu73tvdVaqOyfj+++8BANdcc43TYwYPHgwA+O6779SU9nkFxQUOnwQAUGItcdiN+mNtvesbtbbe9Y1aW+/6Rq2td32j1naHqibjp59+AgDEu/BtblXHbZw+fVpNaSIiIvJxqpqMqiuBFhcXOz2malm9LilORERE3qGqyah6B2Pnzp1Oj6latqEzUYiIiMi/qWoyhg4dChFBamoqLBZLg8tbLBakpqZCURQMGTJETWkiIiLycaqajIkTJwIAjhw5gvvvv9/hxybFxcUYN24cDh8+XGMsERERBSZVp7Bec801uO+++7BixQqsXr0a27dvx2OPPYZhw4YhPj4eiqLgzJkz2Lp1K9577z389NNPUBQFd911F5KSkrSaAxEREfkgVU0GAHzwwQcoKCjAV199hdOnT2Pu3Ln1Llf1jfIjR47EkiVL1Jb1ebGRsQgPDm/wXObYyNiAqq13faPW1ru+UWvrXd+otfWub9Ta7lDdZISHh+PLL7/EP/7xD7z++ut2T01NTEzE7NmzMXXqVCiKorasz2sT0waHph3S5apsetbWu75Ra+td36i1Xa3vzLFrnqqtNX/a7qytH9VNBlB5Jc8ZM2bgySefRF5eHnbv3o2CgsoNEBsbi/79+6NPnz6GaC6qaxPTRrcdrWdtvesbtbbe9Y1aW+/6Rq2td32j1naVJk1GFUVR0K9fP/Tr10/L1RIREZEf4le9ExERkUewySAiIiKPcOrjkpMnT9r+3qZNm3rvd0f1dREREVFgcarJaN++PYDKYy6qf+dI1f3uqL0uIiIiCixONRlV17hw9n4iIiIip5qMxYsXu3S/4ZSXA1lZQH4+EB8PDB0KBAW5NF7JzETrrVuhREUBI0a4PN7t+hpk123uPpBdzXbXbZ9rkJ37nPucv+OcH++3z3ctiB/KzMyUm2++WeLj4wWArFmzxqXxZrNZAIjZbFYfJi1NJCFBBPjfLSGh8n5fH8/szO5P2bndmJ3ZvZfdAVdeQ/2yyVi/fr0899xzkpaWpm+TkZYmoig1dyJQeZ+iNLwz9RzP7MzuT9m53Zid2b2XvQF+0WSsWrVKk/Xo1mRYrXW7xNo7MzGxcjlfG8/szO5P2bndmJ3ZvZfdCT7bZFRUVMhHH30k3bt3l6CgIE3W6UyTUVJSImaz2XY7deqUAJCCggIpKytz62ZJT7e/E6vdLOnpPjee2Zndn7JzuzE7s3svuzO3goICcbbJ0PSy4vZYrVYsWbIEr7zyCo4dOwYR8er3mMyfPx/z5s2rc/+mTZsQGRnp1jpbb92KgU4sl7dhA05fuuRT45ldn/HM7n+11Y5ndn3GM7v7451RXFzs/MKuvnNw8eJFefPNN2X06NHSs2dP6du3r9x2223y8ccfi7XW2y8VFRXywQcfSNu2bcVkMonJZBJFUSQsLEyeeOIJV0vXC+A7Gf7WKTM7s/tLbWZndqNld+bmyjsZcOUFfc+ePdK6dWtbw1D7NmTIECkqKhIRkcOHD8vAgQNrNBcREREybdo0OXXqlCtlHU8AOh+TUd/BNYDzn5vpMZ7Zmd2fsnO7MTuzey+7EzxyTMbFixelXbt2oiiK3ZvJZJLx48fL/v37pXnz5rbmolGjRjJr1iw5e/as25OyOwH4wNkltXemq0cA6zGe2Zndn7JzuzE7s3svewM80mSkpqbaGonrr79esrKy5NKlS1JSUiK7du2ScePGiaIoEhoaKj179hRFUSQoKEimTp0q58+fVzWh2oqKimT37t2ye/duASBvvvmm7N69W3788Uenxnv8OhmJierOZfbWeGZndn/Kzu3G7MzuvewOuPIaqoiIOHPsxpgxY7Bx40Z06dIF+/btQ3Bw3WNGq5YBgNDQUKxevRpjxoxx/gARJ23ZsgUjRoyoc//48ePx4YcfNji+sLAQMTExMJvNaNy4sfpAGlyVzZqRgbwNG9B39GgEG+xqeG7P3Qeyq9nuuu1zDbJzn3Of83ec8+P99vluh0uvoc52LgkJCWIymWThwoV2l8nOzra92zFlyhRnV+11mr6ToZGysjL57LPPpKysTO8oXmfUuRt13iLGnbtR5y3CuQfS3F15DTU527n88ssvAIAuXbrYXaZbt262v996663OrpqIiIgCkNNNRklJCQAgNjbW7jLNmjWz/b1Vq1YqYhEREZG/c7rJcFV9x2wQERGRcXisySAiIiJjc/nthtTUVLRo0UKT5V588UVXyxMREZGfcLnJWLBggcPHq76TpKHlADYZREREgcylJkOcu6SGU7z5BWlERETkfU43GRkZGZ7MQURERAHG6SYjKSnJkzmIiIgowPDsEiIiIvIINhlERETkEWwyiIiIyCPYZBAREZFHsMkgIiIij2CTQURERB7BJoOIiIg8gk0GEREReQSbDCIiIvIINhlERETkEWwyiIiIyCPYZBAREZFHsMkgIiIij2CTQURERB7h9Fe9kwPl5UBWFpCfD8THA0OHAkFBLo1XMjPReutWKFFRwIgRLo93u74G2XWbuw9kV7PdddvnGmTnPuc+5+8458f77fNdC2JAZrNZAIjZbFa/srQ0kYQEEeB/t4SEyvt9fTyzM7s/Zed2Y3Zm9152B1x5DWWToUZamoii1NyJQOV9itLwztRzPLMzuz9l53Zjdmb3XvYGsMlogCZNhtVat0usvTMTEyuX87XxzM7s/pSd243Zmd172Z3AJqMBmjQZGRn2d2L1W0aG741ndmb3p+zcbszO7N7L7gRXXkN5dom78vPVLafneGbXZzyz+19tteOZXZ/xzO7+eI2xyXBXfLy65fQcz+z6jGd2/6utdjyz6zOe2d0frzW33y/xY5oek1HfwTWA85+b6TGe2Zndn7JzuzE7s3svuxN4TEYDND+7pPbOdPUIYD3GMzuz+1N2bjdmZ3bvZW8Am4wGePw6GYmJ6s5l9tZ4Zmd2f8rO7cbszO697A648hqqiIh454MZ31FYWIiYmBiYzWY0btxY/Qo1uCqbNSMDeRs2oO/o0Qg22NXw3J67D2RXs9112+caZOc+5z7n7zjnx/vt890OV15D2WRo0WRowGKxYP369RgzZgxCQkL0juNVRp27UecNGHfuRp03wLkH0txdeQ3l2SVERETkEWwyiIiIyCPYZBAREZFHsMkgIiIij2CTQURERB4RrHeAQHbSfBIFxQV2H4+NjEWbmDYBV1vv+katrXd9o9bWu75Ra+td36i1XcUmw0NOmk+iy9tdUGItsbtMeHA4Dk07pPmTQc/aetc3am296xu1tqv14yO1/b4IbnfjzV3v7e4qflziIQXFBQ6fBABQYi1x2I36Y2296xu1tt71jVpb7/pGra13faPWdgebDCIiIvIINhlERETkEWwyiIiIyCPYZBAREZFHsMkgIiIij2CTQURERB7BJsNDYiNjER4c7nCZ8OBwxEbGBlRtvesbtbbe9Y1aW+/6Rq2td32j1nYHL8blIW1i2uDQtEO6XJVNz9p61zdqbb3rG7W2q/UtFotutbXmT9udtfXDJsOD2sS00W1H61lb7/pGra13faPW1ru+UWvrXd+otV3Fj0uIiIjII9hkEBERkUewySAiIiKPYJNBREREHsEmg4iIiDyCTQYRERF5BE9h1UJ5OZCVBeTnA/HxwNChQFCQS+OVzEy03roVSlQUMGKEy+Pdrq9Bdt3m7gPZ1Wx33fa5Btm5z7nP+TvO+fF++3zXghiQ2WwWAGI2m9WvLC1NJCFBBPjfLSGh8n5fH8/szO5P2bndmJ3ZvZfdAVdeQ9lkqJGWJqIoNXciUHmfojS8M/Ucz+zM7k/Zud2Yndm9l70BbDIaoEmTYbXW7RJr78zExMrlfG08szO7P2XndmN2ZvdediewyWiAJk1GRob9nVj9lpHhe+OZndn9KTu3G7Mzu/eyO8GV11CeXeKu/Hx1y+k5ntn1Gc/s/ldb7Xhm12c8s7s/XmNsMtwVH69uOT3HM7s+45nd/2qrHc/s+oxndvfHa83t90v8mKbHZNR3cA3g/Odmeoxndmb3p+zcbszO7N7L7gQek9EAzc8uqb0zXT0CWI/xzM7s/pSd243Zmd172RvAJqMBHr9ORmKiunOZvTWe2Zndn7JzuzE7s3svuwOuvIYqIiLe+WDGdxQWFiImJgZmsxmNGzdWv0INrspmzchA3oYN6Dt6NIINdjU8t+fuA9nVbHfd9rkG2bnPuc/5O8758X77fLfDlddQNhlaNBkasFgsWL9+PcaMGYOQkBC943iVUedu1HkDxp27UecNcO6BNHdXXkN5dgkRERF5BJsMIiIi8gg2GUREROQRbDKIiIjII9hkEBERkUewySAiIiKP8NsmIzU1Fe3bt0d4eDgGDBiArKwsvSMRERFRNX7ZZKxcuRIzZ87Ec889h927d2Po0KEYPXo0Tp48qXc0IiIi+i+/bDLefPNNTJo0CY8++ii6deuGt956C4mJiViwYIHe0YiIiOi/gvUO4KqysjLk5ubimWeeqXH/qFGjsG3btnrHlJaWorS01PbvwsJCAJVXYbNYLJ4L64KqHL6Sx5uMOnejzhsw7tyNOm+Ac6/+p79zZR5+d1nxM2fOoHXr1vjmm29wzTXX2O5/+eWXsWTJEhw6dKjOmLlz52LevHl17l++fDkiIyM9mpeIiCiQFBcX4/7773fqsuJ+905GFUVRavxbROrcV+XZZ59FcnKy7d+FhYVITEzEqFGjfOq7S9LT0zFy5MiAuLa9K4w6d6POGzDu3I06b4BzD6S5V30a4Ay/azJiY2MRFBSEs2fP1rj//PnzaNmyZb1jwsLCEBYWVuf+kJAQn9vhvpjJW4w6d6POGzDu3I06b4BzD4S5uzIHv2syQkNDMWDAAKSnp+P222+33Z+eno6xY8c6tY6qT4hc6cY8zWKxoLi4GIWFhQHxJHSFUedu1HkDxp27UecNcO6BNPeq105njrbwuyYDAJKTk/HQQw9h4MCBGDx4MBYuXIiTJ09i8uTJTo0vKioCACQmJnoyJhERUcAqKipCTEyMw2X8ssm499578csvv+DPf/4z8vPz0bNnT6xfvx5t27Z1anyrVq1w6tQpREdH2z2Ow9uqjhM5deqUzxwn4i1GnbtR5w0Yd+5GnTfAuQfS3EUERUVFaNWqVYPL+t3ZJYGqsLAQMTExTh2tG2iMOnejzhsw7tyNOm+Aczfq3P3yYlxERETk+9hkEBERkUewyfARYWFhSElJqfdU20Bn1Lkbdd6Acedu1HkDnLtR585jMoiIiMgj+E4GEREReQSbDCIiIvIINhlERETkEWwyiIiIyCPYZHjR/PnzceWVVyI6OhotWrTAbbfdVuer6UUEc+fORatWrRAREYHhw4fjwIEDOiX2jPnz50NRFMycOdN2XyDP+/Tp03jwwQfRvHlzREZGom/fvsjNzbU9Hqhzt1qteP7559G+fXtERESgQ4cO+POf/4yKigrbMoEy961bt+KWW25Bq1atoCgKPvvssxqPOzPP0tJSTJ8+HbGxsYiKisKtt96Kn376yYuzcJ2jeVssFsyZMwe9evVCVFQUWrVqhYcffhhnzpypsQ5/nDfQ8D6v7oknnoCiKHjrrbdq3O+vc3cFmwwvyszMxNSpU/Htt98iPT0dVqsVo0aNwqVLl2zLvPbaa3jzzTfx9ttvIycnB3FxcRg5cqTt+1b8XU5ODhYuXIjevXvXuD9Q5/3bb7/h2muvRUhICDZs2IDvvvsOb7zxBpo0aWJbJlDn/uqrr+Kdd97B22+/jYMHD+K1117DX//6V/zzn/+0LRMoc7906RL69OmDt99+u97HnZnnzJkzsWbNGqxYsQLZ2dm4ePEibr75ZpSXl3trGi5zNO/i4mLs2rULL7zwAnbt2oXVq1fj8OHDuPXWW2ss54/zBhre51U+++wzbN++vd5LcPvr3F0ipJvz588LAMnMzBQRkYqKComLi5NXXnnFtkxJSYnExMTIO++8o1dMzRQVFUmnTp0kPT1dkpKSZMaMGSIS2POeM2eODBkyxO7jgTz3m266SR555JEa991xxx3y4IMPikjgzh2ArFmzxvZvZ+Z54cIFCQkJkRUrVtiWOX36tJhMJtm4caPXsqtRe9712bFjhwCQH3/8UUQCY94i9uf+008/SevWrWX//v3Stm1b+dvf/mZ7LFDm3hC+k6Ejs9kMAGjWrBkA4Pjx4zh79ixGjRplWyYsLAxJSUnYtm2bLhm1NHXqVNx000244YYbatwfyPP+/PPPMXDgQNx9991o0aIF+vXrh0WLFtkeD+S5DxkyBJs3b8bhw4cBAHv27EF2djbGjBkDILDnXp0z88zNzYXFYqmxTKtWrdCzZ8+A2hZmsxmKotjeyQvkeVdUVOChhx7C7Nmz0aNHjzqPB/Lcq/PLb2ENBCKC5ORkDBkyBD179gQAnD17FgDQsmXLGsu2bNkSP/74o9czamnFihXYtWsXcnJy6jwWyPM+duwYFixYgOTkZPzpT3/Cjh078OSTTyIsLAwPP/xwQM99zpw5MJvN6Nq1K4KCglBeXo6XXnoJ48aNAxDY+706Z+Z59uxZhIaGomnTpnWWqRrv70pKSvDMM8/g/vvvt31JWCDP+9VXX0VwcDCefPLJeh8P5LlXxyZDJ9OmTcPevXuRnZ1d57HaXz8vIj7zlfTuOHXqFGbMmIFNmzYhPDzc7nKBNm+g8n8zAwcOxMsvvwwA6NevHw4cOIAFCxbg4Ycfti0XiHNfuXIlPvroIyxfvhw9evRAXl4eZs6ciVatWmH8+PG25QJx7vVxZ56Bsi0sFgvuu+8+VFRUIDU1tcHl/X3eubm5+Pvf/45du3a5PA9/n3tt/LhEB9OnT8fnn3+OjIwMJCQk2O6Pi4sDgDpd7Pnz5+v8L8if5Obm4vz58xgwYACCg4MRHByMzMxM/OMf/0BwcLBtboE2bwCIj49H9+7da9zXrVs3nDx5EkDg7nMAmD17Np555hncd9996NWrFx566CE89dRTmD9/PoDAnnt1zswzLi4OZWVl+O233+wu468sFgvuueceHD9+HOnp6TW+6jxQ552VlYXz58+jTZs2tt95P/74I2bNmoV27doBCNy518Ymw4tEBNOmTcPq1avx9ddfo3379jUeb9++PeLi4pCenm67r6ysDJmZmbjmmmu8HVcz119/Pfbt24e8vDzbbeDAgXjggQeQl5eHDh06BOS8AeDaa6+tc5ry4cOH0bZtWwCBu8+ByrMLTKaav2KCgoJsp7AG8tyrc2aeAwYMQEhISI1l8vPzsX//fr/eFlUNxpEjR/DVV1+hefPmNR4P1Hk/9NBD2Lt3b43fea1atcLs2bPx5ZdfAgjcudeh1xGnRjRlyhSJiYmRLVu2SH5+vu1WXFxsW+aVV16RmJgYWb16tezbt0/GjRsn8fHxUlhYqGNy7VU/u0QkcOe9Y8cOCQ4OlpdeekmOHDkiH3/8sURGRspHH31kWyZQ5z5+/Hhp3bq1fPHFF3L8+HFZvXq1xMbGytNPP21bJlDmXlRUJLt375bdu3cLAHnzzTdl9+7dtrMonJnn5MmTJSEhQb766ivZtWuXXHfdddKnTx+xWq16TatBjuZtsVjk1ltvlYSEBMnLy6vxO6+0tNS2Dn+ct0jD+7y22meXiPjv3F3BJsOLANR7W7x4sW2ZiooKSUlJkbi4OAkLC5Nhw4bJvn379AvtIbWbjECe99q1a6Vnz54SFhYmXbt2lYULF9Z4PFDnXlhYKDNmzJA2bdpIeHi4dOjQQZ577rkaLzCBMveMjIx6f7bHjx8vIs7N8/LlyzJt2jRp1qyZREREyM033ywnT57UYTbOczTv48eP2/2dl5GRYVuHP85bpOF9Xlt9TYa/zt0V/Kp3IiIi8ggek0FEREQewSaDiIiIPIJNBhEREXkEmwwiIiLyCDYZRERE5BFsMoiIiMgj2GQQERGRR7DJICIiIo9gk0FEREQewSaDyA99+OGHUBQFiqLgxIkTuq+H3DN37lzb9icKRGwyyGO2bNli+wVa+xYREYGEhASMHj0aqampuHjxot5xyctqPz/uu+++BsdMmDCBL8p+ZsGCBXZ/D0RFRaFr166YPHkyDh48qHdU8gA2GaSLkpISnD59Ghs3bsTUqVPRs2dP7NmzR+9YpKNPP/0U+/bt0zsGaSwvL8/uY8XFxTh06BDeffdd9OvXDytXrvReMPIKNhnkFVOmTMG+fftst4yMDLzzzjvo3LkzAODHH3/E6NGjUVRUpHNS0ouIICUlRe8YpLGqJiMmJqbG74Dc3Fx8+umnuOaaawAApaWlmDhxIk6fPq1jWtIamwzyihYtWqBnz5622/Dhw/HEE09g3759uO666wAA+fn5WLhwoc5JSQ+xsbEAgDVr1mDXrl06pyGtVFRUYP/+/QCAXr161fgd0L9/f9x9993IysrC1VdfDQC4fPkyPvnkEz0jk8bYZJCuQkNDMXfuXNu/09PT9QtDupkxYwbCwsIAgO9mBJBDhw6huLgYANC7d+96lzGZTJgyZYrt3wcOHPBKNvIONhmku/79+9v+furUqQaX37FjBx577DF07twZjRo1sh08NnXqVBw5cqTO8haLBXFxcVAUBaNHj25w/fv377cdmPbyyy/b7q99JkBJSQn++te/on///oiOjkZ0dDQGDRqEt99+G1artcE6ZWVlSE1NxYgRI3DFFVcgNDQUcXFxGDNmDD766CNUVFQ0uI6G/Pbbb3jmmWfQtWtXREREoEWLFrjhhhvw73//26nxZ86cwTPPPIP+/fsjJibGlrFXr14YN24cPvzwQxQWFqrOmZCQgMcffxwA8MUXX2D79u0ur6PqoNB27do5XK6hM2pq7+fCwkLMnTsXvXr1QqNGjdCyZUuMGTMG27ZtqzHu/PnzeP7559GjRw9ERUWhefPmGDt2LHbv3u30HC5cuICUlBT06NEDjRo1QrNmzTB8+HB8/PHHTq/D1Z8Pe/M2m834v//7P/Tr1w9NmjSBoij48MMPnc4B1Dweo1evXnaXa9u2re3vzvzsfP/995g1axb69euHpk2bIiwsDO3bt8ekSZPw/fffu5SRPEyIPCQjI0MACABJSUmxu9zly5dty/Xu3dvuchaLRaZMmWJbtr5bSEiILFy4sM7Y2bNnCwAxmUzy008/Ocz91FNPCQAJCgqqsWxKSoqtztmzZ6VPnz52c9xyyy1SXl5ut8aJEyekW7duDucyZMgQ+eWXX+odv3jxYttyx48fr3eZAwcOSHx8vN31P/LIIw7Xs3XrVmncuLHDjABk7dq1DrenPdWfH4sXL5YzZ85IRESEAJBRo0bVO2b8+PG2MfYea9u2rcO6DW276vv55MmT0rlz53rnHRQUJJ9++qmIiOzZs0dat25d73JhYWGyefPmerNUr3Xs2DHp2LGj3e181113icVisTsvNT8ftbMcPnxY2rVrV2f84sWLHW7b2p5++mnb2G+++cbucitXrrQt98ILL9hdrrS0VKZPny4mk8nuHMPDwyUtLc2lnOQ5fCeDdPfdd9/Z/u7of6GTJk3CggULAACjR4/GRx99hB07diAnJweLFi1Cjx49YLFY8Pjjj2Pt2rU1xj766KMAKj8jXrp0qd0aFosFH330EQBg1KhRaN26db3L3XHHHTh48CCefPJJpKenIzc3F8uXL0e3bt0AAGvXrsWiRYvqHXvx4kVcd911tlP2brvtNnz++efYuXMn/v3vfyMpKQkAkJ2djZtvvhnl5eV289pjNpvx+9//Hvn5+QCAe++9F+vXr8fOnTuxfPlyDBw4EB988AFSU1PrHV9aWor77rsPhYWFiI6OxtNPP40NGzYgNzcX3377LVauXImZM2ciMTHR5Wz2xMfH294237RpE7KzszVbt7vuvvtu/PTTT3j22WeRmZmJnJwc/O1vf0Pjxo1RXl6OSZMm4fjx47j55ptx+fJlvPTSS8jOzsb27dsxb948hIaG2g5oLCsrc1jr3nvvxfHjxzF58mR89dVXyMnJwfvvv287OHrVqlVITk62O17Nz0dtd911F06fPo3p06cjPT0dO3fuxCeffIIuXbq4tP2qv5PRs2dPu8t99tlntr+PHTu23mXKy8txxx134J///CcqKiowZswYLF26FNu2bUN6ejrmzZuHmJgYlJSU4MEHH+R1X3yF3l0OBS5n38kYN26cbbmlS5fWu8yqVatsyyxatKjeZS5fvizXXXedAJB27drV+V/fsGHDBIB06tTJbpbVq1fb6qxatarGY9X/pxcSEiIZGRl1xv/yyy/SsmVLh+/K/PGPf7St5/nnn6/zeEVFhTzwwAO2ZVJTU+ss09D/xpOTk22Pv/zyy3UeLysrk1GjRtX4H2D19WzevNmpdyosFouYzWa7jztS+50MEZHz589LVFSUAJARI0bUGePtdzLCwsLk22+/rbPMunXrbMtcccUVEhsbK0ePHq2z3L/+9S/bcqtXr3ZYC4AsX768zjKFhYW2d81MJpPs3bu3zjJa/HxUz2IymWTTpk31rscVVT8L7dq1s7vMmjVrbO9M3HXXXXaXmzlzpgCQqKgoSU9Pr3eZ3bt3S3h4uACQWbNmqc5P6rHJII9x1GQUFBRIVlaWjB492rbM4MGDpaysrN51DRgwQADI7bff7rDmd999Z1tf7V9ES5cutT2WnZ1d7/hbbrlFAEhsbKyUlpbWeKz6L+Hk5GS7GZ555hnbchcuXKjxWElJiTRp0kQASPfu3cVqtda7DrPZLM2bN7ctV5ujF8qSkhJp2rSprdGx97HNqVOnJCQkpN71fPzxx7b73W0iGlJfkyEiMmfOHNv9X3/9dY0x3m4y5syZY3c9bdu2tS33zjvv1LtMcXGx7UXvqaeecljr5ptvtltr+/bttuX+8Ic/1Hlci5+P6lkeeeQRh+txRn5+vm19t9xyS43HSkpKZO/evfLUU09JUFCQAJUfDxYWFta7rpycHFEURQDIV1995bDuTTfdJADkmmuuUT0HUo8fl5BXzJs3r8aV/mJjYzF06FBs2LABwcHBePDBB7Fx40aEhITUGXv69Gnk5uYCAO655x6Hdbp162Y7HfI///lPjcfuuusuNGnSBACwePHiOmPPnTuHDRs2AAAefPBBhIaG2q3zwAMP2H1swIABtr8fP368xmO5ubm4cOECgMoDFYOCgupdR+PGjW1z/e6772wfezgjNzcXv/32GwBg/PjxMJnq/zFPSEjAqFGj6n0sPj7e9vf6tpUnzZ49G9HR0QCAF154wau1a3N0FdKqsyUURbH7vIyIiECnTp0AAMeOHXNYa+LEiXYfGzRoEHr06AEA+Oqrr2o8ptXPR3WOnt/Oqn7A69q1a2v8/IeHh6N3797429/+hr59++Jf//oXMjIybPu9tmeffRYigrvuugvXX3+9w7pt2rQBANvPAOmLTQbprnPnznj66afRuHHjeh/fuXOn7e/jxo2ze4niqltBQQEA4OzZszXWExERgfvvvx9A5dUlL126VOPxZcuW2Y5sf+SRRxxm7tq1q93HmjVrZvt77YuLVV0zAACuuuoqhzWqP159XEOqXzXzyiuvdLjsoEGD6r1/yJAh6NChAwBg5syZGDRoEObPn49t27Y1eGyBWs2bN8fMmTMBAN988w2+/PJLj9ZzpOp4iPpUNayxsbFo2rRpg8s1dKE5Z/fVkSNHauwDrX4+qrN3uqkrHF3ps7qLFy9i9OjRCA4OrvfxM2fOYPPmzQAqj0tpaH5Vx6VUbXfSF5sM8orqV/zcvXs31q1bhyeeeAIhISH47rvvMHz4cBw6dKjesefPn3erZtX5+dU99thjACp/4aelpdV4rOp/7FdeeaXD0+0AIDIy0u5j1d85qH3Q5q+//mr7e8uWLR3WiIuLq3dcQ6r/D65FixYOl7WXISQkBGvXrrUdyJqTk4M//elPuPbaa9GkSROMHj0ay5cvd+ugVGckJyfbXiRefPFFj9RwhjP72dEy1ZdraFs5u69EpMY+1vLno4qjpslZ1ZuMLVu22H7+t2/fjmXLlqFfv34AKq+lMX78eLvrWbduHUTE5fodO3Z0eQxpr/7WkUhjVVf8rNK3b1+MGTMGt9xyC2699Vb8+uuvuP/++7Fjx446HyFU/+X88ccfO/2/rPp+Ufbt2xcDBgxAbm4uFi9ejIcffhgAsH37dttZLg29i6GVhr7ky51frLXHqanRvXt37Nu3D2vXrsXatWuRmZmJH374AZcvX8bGjRuxceNGvPnmm1i/fn2DL5CuatKkCZKTk/Hiiy9ix44d+OKLL3DzzTdrWsPXuLuvtPz5qGLvYzxXVDUZsbGxtjOmqgwaNAh33XUXBg4ciAMHDiArKwu5ubk1PmqsUvUuXuvWrbFx40an619xxRXuhyfNsMkgXd10002YPHkyUlNTsWvXLnz44YeYNGlSjWWaN29u+7uiKA5PhXPGo48+itzcXGRmZuLYsWPo0KGD7V2MiIgIjBs3TtX6Han+UcrZs2cdvh1/7ty5ese5UuPcuXMOazT0v+CgoCDcdtttuO222wBUXvp9w4YNSE1NRW5uLnJzc/HEE09gzZo1Tudz1syZM/H3v/8dv/zyC1588UXcdNNNDpevesegoYuY1f6YzFecO3fO4SnBVftKUZQaDYLWPx9aKC4uxtGjRwHA9o5FbeHh4Xj++edtP28fffRRvU1G1c9BXFycT8yNXMOPS0h3KSkpiIqKAlB5gGjtz/yr/5LatGmT6nr3338/IiMjISJYsmQJLl++jBUrVgCovP5FTEyM6hr2VP8l2dBVLXfs2FHvuIZU/6gnJyfH4bINPV5bfHw8HnnkEfznP/+xXan1iy++wOXLl11ajzOio6Mxe/ZsAJUHEVa/loK95QHYDqy1x97Hcnpzdl916tSpxkHJWv98aGHPnj22Zq9v3752lxs7diwaNWoEAFi9enW9y1Q9t3ggp39ik0G6a9GiBZ544gkAlZcVX7JkSY3Hf/e736F79+4AgBUrVuDkyZOq6lU/c2PJkiVYtWoVzGYzANR5F0VrAwYMsB1rsGTJEruf0xcVFeHTTz8FUPmxRfWzPZypUfU/3WXLltl9m/306dNuvyiFhITY3gK3Wq0NvrC7a9q0abaPYlJSUhx+vNO+fXsAldvOXiNRVlZW51gcX1H7eV/dzp07bR8b3HDDDTUe0/rnQwvVj8ew904GUPnOYdV8Tp48ib1799ZZpuqCeMeOHbO9O0L+g00G+YTZs2cjPDwcAPDKK6/UefF9/vnnAVR+X8gdd9yBn3/+2e66SktLkZqaipKSErvLVF0B9Mcff8TTTz8NoPJFavjw4Wqm0aCwsDBb7QMHDmDevHl1lhERTJs2zXYWwLRp01yuUXU6ZF5eHv7617/WWcZqteKxxx6ze6ZIVlaWw1/oZWVlyMzMBAA0atTIY59/R0VFYc6cOQAqz5pZv3693WWrf+7/xhtv1HlcRDBjxgycOXNG+6Aa+Pzzz22NZXUXL160fa+LyWSyNeTVaf3zoVb1JsPROxkAanwMVt+VSEeOHGn7+5QpUxzmvnTpku3sEvIR+lyeg4zA2St+Vpk6darDK39WvxBTbGysPPfcc7Jp0ybZvXu3ZGdny5IlS+TRRx+VZs2aCQApKipyWK/2d4f8+c9/drh89YsVOTvv+q4KWlhYKB06dLAtc/vtt8vatWslNzdXVq1aJcOHD69xgbL6LtjV0AWlLly4IAkJCbZlxo0bJxs2bJDc3Fz55JNP5MorrxQAtj9rryclJUVMJpMkJSXJa6+9Jhs3bpTc3FzJzs6WDz74QAYNGmQbN3PmTIfbw5nt5Og7MYqLi+v9Dpb6XH311bbHx48fL19//bXk5ubKihUrbNt18ODBTl+MyxFnL/6VlJQkACQpKclhrYEDB0pQUJD84Q9/kK+//lp27twpH3zwgXTp0sW2zPTp0xvM4+7Ph7PzdkbV8yMyMtLuBeeqnD592nahrauvvrrO4+Xl5dKvXz9btq5du8rbb78t27Ztk127dklGRoa888478uCDD0p0dLQMHjxYdX7SDpsM8hhXm4yTJ09KaGio7RdJ7StVWq1Wefrpp21XCHR0i4qKkuLiYof1Xn/9ddvyJpNJTp486XB5rZoMEZHjx49L165dHc7h2muvVfUFafv375e4uDi76584caLd9dS+3LW92x133CGXL192uD2c2U4NffHWP//5T6eajIMHD0qLFi3s5k1OTnbpip+OaN1kHDt2TNq3b283+5133unwC9LU/nxo1WSUl5dLZGSkAJCrrrrKqTH9+/e3/RyeO3euzuOnTp2S3r17O/WcfPLJJ1XlJ23x4xLyGYmJibbz5b///vs6n50HBQXh1VdfxXfffVfja56DgoIQHR2NHj164IEHHsCSJUuQn5+PiIgIh/Ueeugh299Hjhyp6Zd9NaRdu3bYs2cP3n77bSQlJaF58+YICQlBy5YtceONN2LZsmXYunWrS2eV1NajRw8cOHAATz/9NDp16oSwsDDExsZixIgRWL58OT744AO7Y59++mmsX78eTz31FK6++mq0adMG4eHhCA8PR7t27XDvvfdi3bp1SEtLs33M5UmPPfaYU/una9eu2LVrF6ZMmYK2bdsiNDQUV1xxBW688UasW7eu3o9RfEX79u2Rm5uLP/3pT+jWrRsiIyMRExODYcOG4aOPPsKqVavsXrAK0P7nw12HDx+2XYOjoY9KqlR9ZFJRUYF169bVeTwhIQE7duzAokWL8Pvf/x4tW7ZESEgIwsLC0KpVKyQlJeGZZ55BdnY2/v73v2s2F1JPEXHzZHwiP7d582bbQWcrV65s8JLMRETkGr6TQYZV9T/55s2b2/16aSIich+bDDKkEydO4N///jeAyi+mCgsL0zkREVHg4cclZBinT59GcXExjh8/jmeeeQa7d+9GeHg4jh49ajsXn4iItMPLipNhPPDAA7ZrO1T585//zAaDiMhD2GSQ4URGRqJz586YOXOmw29/JCIidfhxCREREXkED/wkIiIij2CTQURERB7BJoOIiIg8gk0GEREReQSbDCIiIvIINhlERETkEWwyiIiIyCPYZBAREZFHsMkgIiIij2CTQURERB7BJoOIiIg84v8DI1sqizPj3X8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 600x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dRe = 5.\n",
    "dRi = 0.4\n",
    "\n",
    "# Train/Test Parameters\n",
    "Re_train_test = [np.arange(15,  150+dRe/2, dRe), np.arange(15+dRe/2,  150+dRe/2, dRe*2)]\n",
    "Ri_train_test = [np.arange(0.2,   5+dRi/2, dRi), np.arange(0.2+dRi/2,   5+dRi/2, dRi*2)  ]\n",
    "\n",
    "mu_train_test = [np.meshgrid(Re_train_test[kk], Ri_train_test[kk]) for kk in range(len(importing))]\n",
    "\n",
    "fig = plt.figure(figsize=(6,6))\n",
    "\n",
    "plt.plot(mu_train_test[0][0].flatten(), mu_train_test[0][1].flatten(), 'ro', label='Train')\n",
    "plt.plot(mu_train_test[1][0].flatten(), mu_train_test[1][1].flatten(), 'gs', label='Test')\n",
    "plt.xlabel('Reynolds Number $Re$', fontsize=20)\n",
    "plt.ylabel('Richardson Number $Ri$', fontsize=20)\n",
    "plt.grid()\n",
    "plt.legend(framealpha=1, fontsize=15, loc='upper right')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us generate the mesh for importing OpenFOAM dataset into *dolfinx*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Info    : Reading 'cavity.geo'...\n",
      "Info    : Done reading 'cavity.geo'\n",
      "Info    : Meshing 1D...\n",
      "Info    : [  0%] Meshing curve 1 (Line)\n",
      "Info    : [ 30%] Meshing curve 2 (Line)\n",
      "Info    : [ 50%] Meshing curve 3 (Line)\n",
      "Info    : [ 80%] Meshing curve 4 (Line)\n",
      "Info    : Done meshing 1D (Wall 0.000607034s, CPU 0.000692s)\n",
      "Info    : Meshing 2D...\n",
      "Info    : Meshing surface 1 (Transfinite)\n",
      "Info    : Done meshing 2D (Wall 0.00732976s, CPU 0.007374s)\n",
      "Info    : 16384 nodes 32770 elements\n",
      "Info    : Optimizing mesh (Netgen)...\n",
      "Info    : Done optimizing mesh (Wall 3.26009e-06s, CPU 1e-05s)\n"
     ]
    }
   ],
   "source": [
    "mesh_comm = MPI.COMM_WORLD\n",
    "model_rank = 0\n",
    "\n",
    "# Initialize the gmsh module\n",
    "gmsh.initialize()\n",
    "\n",
    "# Load the .geo file\n",
    "gmsh.merge('cavity.geo')\n",
    "gmsh.model.geo.synchronize()\n",
    "\n",
    "# Set algorithm (adaptive = 1, Frontal-Delaunay = 6)\n",
    "gmsh.option.setNumber(\"Mesh.Algorithm\", 6)\n",
    "gdim = 2\n",
    "\n",
    "# Linear Finite Element\n",
    "gmsh.model.mesh.generate(gdim)\n",
    "gmsh.model.mesh.optimize(\"Netgen\")\n",
    "\n",
    "# Import into dolfinx\n",
    "model_rank = 0\n",
    "domain, ct, ft = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, model_rank, gdim = gdim )\n",
    "gmsh.finalize()\n",
    "\n",
    "########################################################################################################\n",
    "\n",
    "tdim = domain.topology.dim\n",
    "fdim = tdim - 1\n",
    "domain.topology.create_connectivity(fdim, tdim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us define the functional space onto which the OpenFOAM data are projected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "fun_spaces = [FunctionSpace(domain, ('Lagrange', 1)), FunctionSpace(domain, ('Lagrange', 1)), FunctionSpace(domain, ufl.VectorElement(\"CG\", domain.ufl_cell(), 1))]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us import the pressure $p$, temperature $T$ and velocity $\\mathbf{u}$ fields. The mapping between OpenFOAM and *dolfinx* is performed using $N$-dimensional interpolation implemented in *scipy*.\n",
    "\n",
    "`pyforce` comes with a class called `ImportOF`: this class requires the path to the OpenFOAM case, there is also an option to specify if the centroids of the cells must be imported. The class uses two important methods:\n",
    "\n",
    "- `import_field`: imports the field from OpenFOAM into numpy arrays\n",
    "- `foam_to_dolfinx`: interpolates the OpenFOAM field into the dolfinx mesh using *scipy*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Import U - TrainSet: 364.000 / 364.00 - 1.324 s/it\n"
     ]
    }
   ],
   "source": [
    "from pyforce.tools.write_read import ImportOF, StoreFunctionsList as store\n",
    "\n",
    "impo_i = 0 # for TrainSet\n",
    "impo_i = 1 # for TestSet\n",
    "\n",
    "dolfinx_path = './Snapshots/'+importing[impo_i]+'_'\n",
    "snaps = {var_names[field_i]: FunctionsList(fun_spaces[field_i]) for field_i in range(len(var_names))}\n",
    "\n",
    "for field_i in range(len(var_names)):\n",
    "    bar = LoopProgress('Import '+var_names[field_i]+' - '+importing[impo_i], final = len(Re_train_test[impo_i]) * len(Ri_train_test[impo_i]))\n",
    "\n",
    "    caseI = 0\n",
    "    for Re_i in range(len(Re_train_test[impo_i])):\n",
    "        Re = Re_train_test[impo_i][Re_i]\n",
    "        for Ri_i in range(len(Ri_train_test[impo_i])):\n",
    "            Ri = Ri_train_test[impo_i][Ri_i]\n",
    "\n",
    "            path_ = path+importing[impo_i]+'/'+f'Case_{caseI+0:03}_Re{Re:.2f}_Ri{Ri:.2f}'\n",
    "\n",
    "            oF = ImportOF(path=path_, extract_dofs=True)\n",
    "            of_snaps = oF.import_field(var_names[field_i], vector=is_vector[field_i], verbose=False)[0]\n",
    "\n",
    "            # Projection in dolfinx\n",
    "            dolfinx_snap = oF.foam_to_dolfinx(fun_spaces[field_i], of_snaps, variables=['x', 'y'],  cut_value = oF.of_dofs[2,0])\n",
    "\n",
    "            for uu in dolfinx_snap._list:\n",
    "                snaps[var_names[field_i]].append(uu)\n",
    "            \n",
    "            bar.update(1)\n",
    "            caseI += 1\n",
    "    del bar\n",
    "    \n",
    "    store(domain, snaps[var_names[field_i]], var_names[field_i], dolfinx_path+var_names[field_i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us normalise the temperature field, using `min-max` to have it scaled between $(0,1)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyforce.tools.write_read import ImportH5\n",
    "\n",
    "field_i = 1\n",
    "field = var_names[field_i]\n",
    "\n",
    "for imp in importing:\n",
    "    dolfinx_path = './Snapshots/'+imp+'_'\n",
    "\n",
    "    T_snaps = ImportH5(fun_spaces[field_i], dolfinx_path+field, field)[0]\n",
    "\n",
    "    if imp == 'TrainSet':\n",
    "        _min = min([np.min(snap) for snap in T_snaps._list])\n",
    "        _max = max([np.max(snap) for snap in T_snaps._list])\n",
    "\n",
    "    T_norm_snaps = FunctionsList(T_snaps.fun_space)\n",
    "    for snap in T_snaps._list:\n",
    "        T_norm_snaps.append((snap - _min) / (_max - _min))\n",
    "        \n",
    "    store(domain, T_norm_snaps, 'norm_'+field, dolfinx_path+'norm_'+field)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mp",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}