ogs-jupyter-lab-h2m-2d-liakopoulos.ipynb

{
 "cells": [
  {
   "cell_type": "raw",
   "id": "ac287b2f",
   "metadata": {},
   "source": [
    "title = \"TH2M:HM process: \"\n",
    "date = \"2022-08-16\"\n",
    "author = \"Norbert Grunwald, Olaf Kolditz\"\n",
    "notebook = \"TH2M/H2M/Liakopoulos/ogs-jupyter-lab-h2m-liakopoulos.ipynb\"\n",
    "[menu]\n",
    "  [menu.benchmarks]\n",
    "    parent = \"HM proceses\"\n",
    "<!--eofm-->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c3848074",
   "metadata": {},
   "outputs": [],
   "source": [
    "#1-modules\n",
    "from IPython.display import Image\n",
    "import os\n",
    "import pyvista as pv\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.special import exp1\n",
    "import vtk\n",
    "from vtk.util.numpy_support import vtk_to_numpy\n",
    "import matplotlib.tri as tri"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ffd2a84d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#2-settings (file handling, title, figures)\n",
    "data_dir = os.environ.get('OGS_DATA_DIR', '../../../Data')\n",
    "fig_dir = \"/home/ok/figures/\"\n",
    "ogs_bin_dir = \"/home/ok/ogs-devs/build/release/bin/\"\n",
    "out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', '_out')\n",
    "prj_name = \"liakopoulos_TH2M\"\n",
    "pvd_name = \"result_liakopoulos\"\n",
    "prj_file = f\"{data_dir}/C/{prj_name}.prj\"\n",
    "prj_file_test = \"liakopoulos_TH2M.prj\"\n",
    "vtu_mesh_file = \"domain.vtu\"\n",
    "title = \"H2M process: Liakopoulos benchmark\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "688250b6",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABQ0AAAHICAYAAADk7syNAAAAAXNSR0ICQMB9xQAAAAlwSFlzAAAXEgAAFxIBZ5/SUgAAABl0RVh0U29mdHdhcmUATWljcm9zb2Z0IE9mZmljZX/tNXEAAGKGSURBVHja7d0vVORIuwfgliORSCQSORKJRCKRuGnOiEXyKa4bORI5EolYEVBIJBKJRCL3knRmlrD86e76k0ry7DnPufd8Ozskb5Ki8+uqemf/+9//ZgAAAAAAvykCAAAAANChCAAAAABAhyIAAAAAAB2KAAAAAAB0KAIAAAAA0KEIAAAAAECHIgAAAAAAHYoAAAAAAHQoAgAAAADQoQgAAAAAQIciAAAAAAAdigAAAAAAdCgCAAAAANChCAAAAABAhyIAAAAAAB2KAAAAAAB0KAIAAAAA0KEIAAAAAECHIgAAAAAAHYoAAAAAAHQoAgAAAADQoQgAAAAAQIciAAAAAAAdigAAAAAAdCgCAAAAANChCAAAAABAhyIAAAAAAB2KAAAAAAB0KAIAAAAA0KEIAAAAAECHIgAAAAAAHYoAAAAAAHQoAgAAAADQoQgAAAAAQIciAAAAAAAdigAAAAAAdCgCAAAAANChCAAAAABAhyIAAAAAAB2KAAAAAAB0KAIAAAAA0KEIAAAAAECHIgAAAAAAHYoAAAAAAHQoAgAAAADQoQgAAAAAQIciAAAAAAAdigAAAAAAdCgCAAAAANChCAAAAABAhyIAAAAAAB2KAAAAAAB0KAIAAAAA0KEIAAAAAECHIgAAAAAAHYoAAAAAAHQoAgAAAADQoQgAAAAAQIciAAAAAAAdigAAAAAAdCgCAAAAANChCAAAAABAhyIAAAAAAB2KAAAAAAB0KAIAAAAA0KEIAAAAAECHIgAAAAAAHYoAAAAAAHQoAgAAAADQoQgAAAAAQIciAAAAAAAdigAAAAAAdCgCAAAAANChCAAAAABAhyIAAAAAAB2KAAAAAAB0KAIAAAAA0KEIAAAAAECHIgAAAAAAHYoAAAAAAHQoAgAAAADQoQgAAAAAQIciAAAAAAAdigAAAAAAdCgCAAAAANChCAAAAABAhyIAAAAAAB2KAAAAAAB0KAIAAAAA0KEIAAAAAECHIgAAAAAAHYoAAAAAAHQoAgAAAADQoQgAAAAAQIciAAAAAAAdigAAAAAAdCgCAAAAANChCAAAAABAhyIAAAAAAB2KAAAAAAB0KAIAAAAA0KEIAAAAAECHIgAAAAAAHYoAAAAAAHQoAgAAAADQoQgAAAAAQIciAAAAAAAdigAAAAAAdCgCAAAAANChCAAAAABAhyIAAAAAAB2KAAAAAAB0KAIAAAAA0KEIAAAAAECHIgAAAAAAHYoAAAAAAHQoAgAAAADQoQiM9+aezTaf7a5oU+0AGP3vyL/+3pwdV7v/8b3aVh9g8GPc92rrzTFuMc5tqRHAkuOpIjDYm3cR8u0/O33241n17PbZPxHct3/fRfv3n7Q/b0PtASj0JXmnfSk+nM2r09nx9c/Z/Kp6dvvsn7V9u3pq/57q+e88W/zd1b6Xb6DH8W5rMQ7VY1LQOPf4Z3ybX50//51H9ViqxgDteKsIDOJGXcwarAPCs2c3kYLBkEDxd5i4J0gEIPvvxb/+3nh+ud1rArz6ZXcR7P3Tk/aluw0TBYlA/PFutxljvl1dPntIPq79+bKkGdf26mNwLYBJjsGKQLE352Jm3882pPuncDftbMSvrh0A0X8n/t/fX5pArp4JM7+67zEgXCVIvGhn7Wy5hsDK49736msz5vX7pcjLIPFuNr/+YRsHYFJjsSJQ1A35b1D4MICg8D2Pz37VsxBdUwDW/p24mF1z8Pyi+quYl+aQl+16ufRxte/aAh+Me/V+q/NFQFf6mFad1MfrugGjHpcVgd5vwsXS49OBB4UfLWU+soQZgKV/LzbLjq8uBh0Sfvyy/dAEiN8rs/OB3+PebvMFyfDGs6dmNqR9EIGxjs+KQG8332y21TYweRphWPhafY7n9Tm79gD853fiYvnxUfGza+K7b2YVPZ+/+wAmOPbVM6oX2y6MYTy7sHQZGN04rQhkv+kWYeH5BILC9/ww8xCA5ndi0wG06f75OLGw8L+zD+uGA5oNwHTGv3qf1hxNTXKrx3RfhABjGasVgWw322y20XY/nmpY+Hrfw7pxig8UAFP8ndjsV/j8Yjn0vQpTNFA5rg7dIzDq8W+z6YI87i9CbgSHwCjGbEUgy402mx2MdM/CGHse2tMJYEq/ExfLkMc3uyb2C7dlfjDG8e9wMjOrj69/uubA4MdtRSDpDTabbT+rhIOfOnW/AIz+ZXn3+UXyVii4QoMB3ZZhPGNgvQXB1Max79WWaw8MeuxWBJLdXIvZhVNochLLTd1J2r0DMLLfh3WTk/n1D0HgurN1qgP3EQx8HJxiYLgYv45cf2DQ47ciEP2mms2+PPspBFx7r8M99xHASH4nfq+2J9gROf4+h2brwHDHwakGhpYoA2MYwxWBqDfUbLbZzpgTAK6vnp25734CGPjvxONqrtGJF2+Y9Dg45cBw4cJ9AAx6HFcEot1Ms9mOZidRHbqvAAb4+7DujDz2zqD59zd8cG/BwMbC4+rE+HV1614ABj2WKwJRbqRFYPgo6BMcAkz69+Fff29qdpLI92rHPQYDGQufn1fjVuPe/QAMejxXBIJvIoGh4BCAxf6Fzy+IXpJTLVHWEAUGMx768kRoCIxjPFcEgm4ggWGuPQ533W8ABf8+XMyqefSCnHRfwzP3GgxgPLQsWWgIjGdMVwTWvnlms22BYTZ1nbfcdwBFviDvaXiSxbn7DQofD79XW8ZDoSEwonFdEVjrxpnNNp7dC/Oyuqvr7v4DKOoFeccLstAQaMdETaCEhsC4xnVFYK0bZza7FOL14sL9B1DI78LFjJoHL8XZ+B0IJY+Jx9VBYWNGvWVE9YZboSHAkmO7IrDyTTObnQ0wbLt9Vr0y1JmSNoIH6Pt34V9/b4yg6cntf16mv13dFXy8lXsPCh0T/+/vL71+iVLP+P529Wt2XM1n36uvSx933cCqCTuvfzz/9zdCQ4BX46QisNINM5vtFxym1bMfT+tuw3XjkGdfVjivuqHLXvvf/yp8r8YHy5QBen85vhlOMPj8Mlw3JjiuduvZkSufb70Eu/5vj6uj5u9ahItPQkPgzzhxXO33FBY+NOPbX39H+WzcjO91iFgHkHGaWwkNgWGP74rA0jfLbLZZWJhWdxU+r2fepQjR2iDxpN1LsLTg8Kd7EqCn34fzq/OCQ8LH2fH1z+YFPtJL9Lt1WMzQ2a+7GmeZofjt6tL9B8WOi1UfHdXrkC/peS3C0IuQL27cH8Cgx3dFYOmbZTEDr5ROwvWMwM2M5/61oPP/7av7EiDz78Lj6rC4oPDfZXn7vdam3uOxno2YrhGCRihQ4rhYf4GQ+8uRFZYgxxvfrs9Wnn34PDa7R4BBj/GKwFI3ymLpbglB2ekqy44T1GH72U0p+zS6NwEyvzTGWa4WNUib/fX3ZnG1+neJ36XQEEY+Ntb7COb8kuS42u3tXJv9bKvTpX8XHF+fuUeAQY/xisCnN8ls9qWApiEPJc2sK6gZzL57FCDT2F/SPoZ1GPe92hlE3f6doXMfuhTRfQhFjo2X2ca+46qIhoBLh4eFHC/A2uOdIvDpTdJ/QHZRYuOPdsly3/sd3rhHATKM+YuXwzIamxxXe4Ot4/Oxr7332XF14l6Ewp7pRWOop0xflhT3ufdPePhWDer/rcCZ4AArjXOKwIc3yKL5yVOPjU6OCq9PPQvzR8/B4Y57FSDhWF93D+5/ZmHdIfRwRDX9unJzgZ73bATefZbzjIMFz66uw8Gmu/y/4eFjvcerewQY/DivCHx4g/QXiNXLobfVabmZmO5VgIRjfFjnzBiB4V29xHekgcPO0vUdaQ1g0M9w3fwo0yzrAdVk170BjGacVwTevTn6m2X4NMTZc+0y6r6Cw8EErAADeyHe63mWYVUvf5tAnXc/2TPy3v0IBT679ey6LGNhdareAD2M84rAuzdHf7Pn9gdar3qp8m1PNbM5PECSF+Kr2x5nGP6q9wubVL2Pq/03G6ZYmgxlPrO5mqAMpPETwOjGeUXgzRujv1mGg97kvK1bH52mH9y3AJHH9OPqsL8Zhtc/1L5ZtlzZFwwKflbr/VZzjIkaigD0M84rAm/eGP10TP45ktrt9BS4DrabZuLrUQe5u89OW/W9Xb3y88W/P2z//KRm90Sued1ZfO+Tmv928eLPnba1/6qORV3LwxfX5/LFtTt/8b/vt9dudzTnvugI+tBLYHhczd1/E3ze6oYSzXL4ulP3s+Prs7bT9Fsu/vy55s9Wu/V/r47Jrs1OU+MmzH5R9/nV+bvXqLl+f67PUXuNdkZXm0zjovuQiYw12+1Ys9cda5ptAF6OMef/+R3g9wCp7ktF4J0XxYfMgVc1ppDm+Z95D6HhL/duU/vttv51uPEYOoOzruuzg2cb6vtmvXfaev9KMG48tKHiyZjCqMKv59d2a4qbwGv32IaKB0Md2/ubZTjtGYYTC6HmzRL02OH04u+7eP77TzRkWPG6/PX3xr8v688v5v92wo3fDb1Z1tu+7A9431KhIaz57NRfTtbPfz0OLJb5P0YcY26aoPG4OqiDSPVGaEjsl8a9HgKvnZHV8EsPwetTrmCrDRbmbShQLdkNu2qDpaNnW5GPZ6udzZZ6aXgdRB5MfHz40gZBv3qaUXvRznwb7JcM7fPzctbe0xLP9u/Zfaexx8s2aP+ZcMx6GmL43tNehrdT28NwYi+HB21I+NTDvXXRBOHur7cD3Prluu5S3neX9Po4BjRTqAlZc9VnosFH3TW+eXaPr3++mGX2uPq9Vf93w7q/XtVhu/kiZHH+96ufe60J6ff6HAcX17OZRZ77M8ZjE0zWs56fj8HYj9CQ0BfIX7lDgJHWsY/ZhoeJz2krUrOXpxjH2s5y+9VDne/b6zuZl6926el5hNmbsYOonQHVMNbzU/sR6Zr20fX9Z+nhYbNMNH9g8GQ2wAjHzmYWSbOE9bHXQOrlfVYHlxNvKtEEEIvZPXdFXJe3Z4uelx7wNAFIvm0b9id1j9ZfNKTqTH1cDWoP+WbMiv18Zd4GpAkrczUNWvZLSl8kITRkzZfIjR5mD43yg2tPsw3PE59TFfl499Y8ju0Ex7Lu8tlR7x3S7lV3U0hQ+O7S/NizVwfy/JwGhO19Pz8Pqb/kCLpWiwYcufcxLLYerPWCuL9YHlZgIPWyO/fEZpy0XwhcFH1d3lpmWOjLfd2cxF6vqcaQZmbh5LfCaLZxSHdPnSR+Pjba2ZH3BY8xj80xCg8RGrLCy+RR5hfHUe/d1MNsw/uBBco/1jiGswJDq9MR3rt7EWfFZRtP6sY3E/pC5maNY/hZ2DWrl2gXNbsu60vwi6WjPoOMJizc62lpe1B4MPbOtO11qYZ1Xd6cHVVcJ3HjZKKwKc+S753ia5EycPt29ZRwzDnsrZna+uOLppoIDVnqpTLn7JOnUl/wI9azj9mGW4nOZavPmZFtiPVQcGB1M4b7ebboNn05sLDwdQOO09KWvyZ6fu5X+PkHBS0tf+t3wbyccKHZayjvB/UBN0HgRdhc1tKz1Web1Mt1R3YvNstnh31d3ltWuFtMjXMGIxNZVp9ti4zj65+Tr0Pk2daDH3OsekBoyCcvlbmXJk+iQ2QbYAx+X8M+Q8NCZxe+F1gNtlHKbNGluNRgaZ3lr8V8Y9pXaNh+cXE+kGt2WsS1yr2MaGL7dI3y9/xi+VkZexaabbK4Js1+cM2ehU8jCwxfOi8h5M28L+QkZhtm2yvy+Xkvug6LLuaDCQ2b3wVjGHMK+lICoSHlBQb7Y5gRV2Bdd8ewr2EfoWE76630/fQGH4i3XxiUsEdkihlsRYS4fYSG7d6FdwO7Zr0Gh4suqnlnDPn8MeDf74slhMNe8vrOkr260/OAA5evxTY4iX+t7vqefZf9GZjIEsps4VPB9Uz+HEdantx8STGuGc3nfscjNOS9F8yce13dTaiuXzLP3rpPdB5ZQ8Pnf74OfNbbfCD3Zx0s3Y8wMCzqWuQODdvlyE+u16ovvxlmNVgGNI7f7YuA+X7UgdQAG080TRPGPbvwvZC3t70Oo3e2XeZ8J7BMOWMIVWRAlGm2ZSUwvBrcDFSEhvT7UpszODibWG0vMr90bww5NBx44PFS0d+Gj6jOxY85OUPDdpn50K9XL2Fa1tlJ9TJQ3QqH+Tv9uDqYTDB1fD2Iz2vtrM+LSYWF/1Gd9vM8pO7y+6bH2feqqCZaCcaZw2whbImduRfbPhT9zIwyMPzN5xOEhmR6of3I7sTqOx96fXOFhj3UKvUeh0V+qK1DmYmEhZ37rZ75O9bQcDas/QuX8TXrNXp+ATWLi2Je5EvbO6/gF8g2MLyddmD4J+T92cMzcdLT+d6PecZh01wp3++j4rYjyLLsPSB4HnVgGFgbhIYItWIFKZP69qKHUHY+kHM4f/UzfowwqLqfldfJd4qB4W+XYwwN28BwyF2v35J1v7+sS5PrmR06JgsMh7Vc7bLIayIw7H25aQ97wRazNDt5bfOFUkU1mMkSmH67Ctoqaza//jHuLyA0Q0FoyH9faHPuZ/hrojXOufz7Z4LjTxoajmyG1GtVKUH5LH/DoxId9VD3ZKHhSAPD7NdqlnUz/+tBNUuiCQz3zWIrK5wRGJZzrZrtFvo954sxfhHTbIUwwS+y6vu35KXJvQblQkOEhvQYJOTsUHs40RrnDMVuExx/ktBwhEsqi9xTr72Gm7NhN5eJ5Sn3svFUoeHIA8PaQ33fZgoNHzO+nG367DGg39+LWS+Pkw+iFg0oiliy1i4NvBEQlnGtmiXs/Z/zw9hmHTb3ea6xp6DGXFlmWH6vtgLu9/F/WSE0RGjIGy+0OYOEYr7Jylzjo5yhyEBCj18jDzxeB1WbPd+DlcDwj5sRhIYPE3l+ks/Kyzxr4NbnjoH9/s46C7X4MOqmiGvST/ONoV2ru1x7URa2dP92TIFHxnu9KuJ86xnEqRtNBYxjTYf2acxWFhoiNKTzMrudc9bIhOu8l/lFeyvy8W8JmsoPPz64fnP1/4/TgYeGU7IznhdeS5OH9cI+kRfETMv6Il2TI9egrPGmXSpe2vlfzL5XXwc/Bj2fQ7aaFTALPsuS7OPqZM37fDN5oCk0RGhIoWHWfsYXvyI30h5hOFvbjXz8Qo+BzjZs772nvmf2tTMdX7ov4Hp8zXQNPD8FL+/POmvpuNr32WMwL+vbvb8gLpbgVq/cF7D09WtP1+TrZF7a4405e1muTbnLxauhByDNrNE898q8gHP9VerS5OfPCmc93cOPb/we+FeK+0NoiNCQVy+zJ1OYaVVIrXO+ZB8IPYp0nvme+zLLu2fpy6XnB7MlGsC0oea8p+XTd7MMTWo8P8GSzlLPuj+RrsnD+H3d15559QtzPdNmiWWlTai5mAlZ9XCc2Za+vromd4UGUyUvU37IMYNsAM2C6mXLB4Mcj46rk1xfUvQ+7qbewzFkaXLOL2x+/y5YMuBsandc7Tb3yiJ4Dauj0BChIa9eZs8zvvj1/g1Wz7W+y1jro8jHLvQYyFLLV9ftNON5Pc4We3duBBzvRv3lQubrsS80HIQkH2Dbl5RsL60+dwzk9/W8Os06i6RechsQKC+Wh17/yDyDbX/E12T5QK4746fvLsLvLtXNc40G0SDivgnbB/QFTtuMKU99emx2VM+KLXZpcr5l4uex7s2mnvVKinXGpREs7UdoSNyX2Zyze/YmXuucTQtOhR7Fush0v33NeU4xl17XwWrGGZJVhmvh+Qn3M8m1qWdrZXtZtZ/hIH5X59xDrN53LeIssKapT74ZklXGa7JTwFLxuya4XGIGzotZP4ezUhrpPNcwQ+hT+mzDV0s+r3+EdNHNOi7l6Cjc6G/P0lmOLtylLk2ux7eEX8Q0z+YqS7+tikBoyKuX2Zz7ig3iF3PCWv8c6jLYgYQeT20Ift7OsNtrnba1r//dQyHHup3hfrvMVPP9hOeQayb0TuJrMZTQ8H72370nH0tZopxiKXmWmQ09zcyi4Jfz1C+IOV6+MwVR7fn0M4OtmaFTnYbOvqqDinam5H1/IVmmpijDmG3436Wghc+sytIgpA3Hexx7U8/UrQLu6/vEnw/y7D1az1pdBKCPJd4DCA0pN8jK1SDhSa2zdrC9jHzspYYeT22w9HXJ8/jSLqG97/m4TxLfaztjmT2cKTj8kfgcSg4N79qxafuD46+XjR9mni2d5X7L2onVN/fl/56uZ7TlCZDTj51ZgsP0QVSPs9fOUzyzvTVTqIPqDPtQDmy24X/2uyt138Ms+/31uDS1mZmbftw9KvL3wvH1z+z1bra0aL7IeCzheBAaUnaItZHxZW/yeznN8naqvol87CWGHj9DlsS2IUhfMw+rxPdajn0BDzM+O6nDqqdZwoYohT4/N+uEcLNF45qbno45ehfljC/w9jMcwu/pHPsCHlf5xs7UsyYzBFHZZ67Vs50Sh7rtDOfH7KFYpu642Wa6pnNf4szwZn+6LPfJ9VkP55b+d/GaW0EkX5GQoVHRh2F0vcfn7yZTdXDuC06Ehrx6+dsZS0gykHrv5lxmOOLQ4y7WctL2vPqadbiZ6D7byDCD+Czzs7OZ4ZyOEh5/Sc/PY4xzrWfLZpyp/lv0/UBX2uNnIPu/dc6vXha56LxZDcKiY+RJM+skd3feeuZF6n3zMr+MN0vRkp9TlW7szD1rrX5xzvTC3Da3yL2UN8uXFyPqdH1bUhfZbPutfrt6yP+FTfKl+1XAOHSYMqCWSSA0pPQQay/jy94v9RYaxggNZgHdeQsLDg8T3Wepl8Hfp5yV90lINcjZ0AU9PzezuA1rDnJ/YZAgNMzVNOI8+zOzCCUeB/vCvpjxle2FvZltkXr2UuYgtD2vk6EGURmfz8X1ydwUo31G77M+V5meqabJVN/Na+KNRZe59u9c4pm4G9N90t4rOxnOZ+0vNxJ3br8o4b4CoSEfvfAdZnzZO1fvrN1sxxganiYOdXIvVb5IdC6pA9C9np6fL+0s05TntpXw/ur7+fmR6NxOMp9H1NAlw8br2fZ+e+PcLkfxwp5vSeV94vPoZ+zMMesrQdi2aB6S7T577CsUamaP5Q3XsoUUiWdo9REe/uq723KGLwGyf9GVOJQLXgKcdLn98/WUSSA0RGgoNOwrOHgY8LH3sn/eLO+ek0n20ctwDhc9P0OpZ1HujzQ0PEx8XXLucRh1g/Z8L5zVadZnpdnkfDQv6k+p93zKsAy237Ez9SzKBPu/ZQkSeg50ewnXMjxP3XPLtA9fVtc/+tqHrp2dmidIzzQzOvky/W9XQc0hk34BJzREaMgAQqycM0TO1DtvcDCS0OMp58y2dvlzzvPbi3z8l0MM1VY4v9R7G/4Yw7PfQ+C+k/F8ou6fljGUOMr6rAy5g2kPMzWTz8rsualC+r0N41+fjMt2i1geWB9HxvHoMO94NMLgsH6eMs2Czj5eZRy3sswoDrzfZ4s9dwc/oxOEhqz7onc6hqWlA6r3ptCw7JCqvUaPGc/v54Cu0UMhz1HKYLRKdMxbYw0Mewjco4YTGV/SD/IGEBm6AGd+QR/sS2sPTQV6CBqijp1NI5xcwU8pe9Xl2Nftz3h0/TP7+eWcOZp3bLrMPeuw/n2Sazl2hnOZJ3/GA5sbJQ0Nv13dyCQQGiI0FBq+rrnQsMDAo8fn4ibicadeulvEbOFZ2uYb0ZeMT+X5meVr9HQZ7Zhz7pmWuQNn0j2Y+vK92h7kS2vmjsm9BA31i3nEpYz5Aqb8e41+8tzmmm1428/5jTY4fMi5xL3ZpzRHk6vIz3X2QC7C0uT2C5dfY56JDkJDSgpHDtVcaDiEkDnzjNDHiMedemnyTiHP0PbQznMqz88sTzOhaDOasoaGiQKv3l7E+pDopTz5Ur9yZrJtD+U8s9y/mff2W/Ia5Ztt2EMn7za8PhpNV+UeQ+hsS74TLmXPsj9jhOPPUOv7vp5HEBqyzAveudBQaFhg6NH7/kKZ9zbcjHC8XxLv9fdY2HOU8lznAw8Nf/Z4XX4OanZuztAw9/K1cYaG0T9HtJ2FUwYYZY2dKc810v5uGa5JtqWXRYbYPc1+fiPAvh3prMNfOQKgput2puXXCQPkw8THHrw0efH7NMMM2fq5FxwiNERoKDQUGi7lNsXy0DXOdy/jtdobwPEWsVH8i/O9TXiu5wmON9fzU/X5/GR6bu4jvrDsZnuR/F5tCQ0LDA2Pq73Ex13W2Jk2qIkydmabbZd5n9EV7sl5pvOf93qezRLbke29mjkAev45d0OekZshIL8Y1DNZj88RQk4QGjLk0FBLeaHhp7PZ6p9T0LW6z3StjiIc64/Exzgv7Dn6NYTlr5mfn/p+3ej5umwJDa+KWKJab64uNFwmREscWvQczLxxX6Tcm6uK9FzOswQhhc7qSb6MvLCZlm1wfz/C8LBKHQA91+4k09h7FP3Yc8wojvTFQIYvl14+lw+5u5uD0JCSQsNTNc8aGj5GPu6tqd0jswF1g82wn+FuYdfmdAih1BSfn8RLx2vRuuhmDQ3zN0Ix03C5EO1yTNf98/si6TK7+0j37sVUlya/qEGOEO2+mPOtA6Q6LK4Dk3GNW7cpw+ksewIunpfoHX6Td4CO+MVAti0TXs9WzbwXMggN6ePFW2jYT3CQJPiYaGg4z3StLiMc62PiY9ya0Nj1IDQMOte7ocyizjajp4elkELDpev0OKZl6b2Ght+uooydWYKjQpcmv7gvz7M8U4UthayPp7lHx9Uo5TxpzXLtgRl5LEvekTjyFwPZ6vzGEuvSvnxCaIjQsOjZVEJDoWHm67WT6Vo9BB5njm7PRS3jSj12CQ2DzvViKPdj1kYomZcbJX8hG0FomGWWTmFLYFNv6D+YZ7LwGTzJG0T8+0ztF3n+zbM5ov0OEyzvfXGvHGQ6h2jbTC32s0z8hU3kLwYy7mv4/qzVuvO4PQ8RGjLy0PBcvWfbQsPhhIbteT9lul4bAceYuvlEUd0/M41dG56ftc81x7YXW4MLDedV1uszO77+KTT89CUw9T5V5Y2dqbuABr7QZgnLvl09lXZd3hibdob4TCWow1aW5eoZ7rm623GSGuXrNn43mLE3wZ6leT8vfHJu9Uxksw8RGiI0HG29dzPW+ybysU81NLzJdL22Ao7xJPVMyPbeLcn5EEIpoaHQ8N+X8+uzvKFhdSQ0/LRGJ4lf7h6afTNLknrZa+ASxixhd4L92aI/v80y3SzP1HwQn58XIWo18DHsPtUssWzL2SPN0E3+nCfas7TAAPu++SKosG0wEBoyvhDrMGOIVfSm0yMMDS+EHlHO+9es/NDwR8b7aiq+en6Kvh+jfUAey75W77xkCw0/vPYjWv5YisDZVFleyo+vfw7iM2Pq5Zs9zICOMK59HfbWC9dJtmrK1tQr0pdfyfctTbTsvvDfq1UzU7vQrvAIDREaLqtS7+TLSF/6GfnYpxoanmW6XrsDCDanZNfzs/a55pjBPrzQ8NvVZfZr8e3qTmj4YX3Gt+9j/9doN/B5rDIc43wgoeHtUEOsDOHhVjNbbWgNU+rj/evvzUT3y32OmW1Rgt+0x/iYtGN16cvlfy9f/l7tyDoQGhLr5W5faJi13kcZ630S+dinGhoelR5S1c+WkC+6Q8/PRELDHJ1a21kAPbxYfx1VF9LoMw1H2GF64NcoS9B9XO0N4jNjnlB70FsHLRqmNN2WH4bznCSabZh6v9JYXwwcX58N+Z4e2Cz+C+EhQkNivNzlXC77oN5Zl5EKPeKcd67ZoQcBx3gv5PP8CA3XfjG/yfThvZemGM2ytUX4UCWTq4bxQ8P7Ab38TSM0zHFNBtJAIEO40oQKo/h8XTcDWezjWv4znWi2YbY9egOX9ye/Rhm+FBjg1haVxikIDSn9RfalSe+zUO8zaHnl4ELDndJDKgGf0FBoGPQCk2+p0Ui/8c/3sho9NBTylRca5ngOtwbxXOWZOTa6VUD1fnazLEu7Q6SabZhh9nQ9q3PN5b91I5UhL03OXus04eG+DAShIau+3H3J/CK+PfF63w211hMODXMF60dCQ6Gh0LCP0DDjjIGB7KUmNKSPa9TMFhMa5g4Nb0f7mXvRHKTMYGcx23AjwTkfZnrO9wq9p7Mtt2+Wxg9qWXz3uQ9tWoXQkOkFWU8ZX8Qn/e1G5lpvRD52oWGB515fZwGf0FBoGPSSNc+5v5DQsIzQsH5hF/IVFhrmuo8ShDWJQq8cY9P9GMekV/fVTt2IqsBnJfo7URO859nHdq1wLvlWFpn3K223/xjuvsH1Fgi6LSM0ZMkXvJz7oc0nXOecS8EfBnr8QsOy76sp2fP8TCY03M84u2SUe/sOMjTMdczTCw33Aq5JlgYDg3mu8swaexzjmPROPfcK6yifZFZc/fdmmSm5YtiUfMzt6fdr2w36cbBjdv1MaJaC0JAlXvCqjC/iPyZc572Mdb4UGkY9d6Hh9FjeP5XQcPGBP98H9IEsjRQasub9vfbY2S4nFRrmDQ0HU4+IdZ0XEvLcJzq/3Uzj8cGKx3WUeNbczx5/B9Z7Nd4PeOx+HEpXeYSG9PeCdzHkMGtAdT7KWOczoaHQkCBfPD8TCQ3zLeeKHnwJDYWGxQlY6laHEEKy7KHhZGYa/uf5z9X1/eOQfSfJ+eUJsC5WPKYq8e+HXfdUgPpz0IpBMEJDphVmnWV8Eb+bcJ3PM9b5IMHxCw3TOlrz2LYFfNHde36mExq2LzO3Q1+SJjRc+Zi3B/tyV66gsdPMul7qMfo9DT/8wihnI6y3x7STNKFhhiY6KzRzSb6HbEFbf9TXdNj7HI7vi02EhsR5wTvI+DL+FHsGz4DqnHPvyB2hx6DOu3ZY+PFNydzzM7HQ8Pj659CXpAkNCz3mKQnsDt7sOSc0zBsafru6mfx7UK6Ow2/X/3LgY/JRETXucWnyB19KVQMey3dlJAgNef2Ct5P5hXxfjZP7kuAchIYFPhdCw+gePD9TDA0zvzSO7AO50JB6pk9oF85se7H99ffmIJ6rHLPF5leVd6EM++31ENpmCq2qJY/lYoq/U9sxbYjh4eMY919GaEjYC96XzC/lk2uGkukl+rfbAYdnJYaGuQLf3YKvy1Q8pVjaLzQcQGiYqWvrbM29oISGQsPCA8Moe2FlWzI+kJfhTKHh6LZLKLze2WaeZ/sy7JMQPvm+wQUtTR5ZeHgb+kUQQkPGF2rdDT3UKry+1dBD2QmHhruZrtvXgYT+Y3Uzi9wxWWg4nNBw8cKYuZtmog3whYZLHnOzn5nAL8ZMqZCOyb3cR0MJDY+vz9I/U2Ut6+w/OEw8I+6NGV3JziVXk69P9mVM3uDo+TkZzP01uPDwenITfRAa8vFL3kXmF/RBfGCLVNuNzLXdE3pEPe+D0p+JTEt2T0fqYN1Znp6fkYWG364ufRifTmjYBgTpl+w2s5dGqA4CIi8JzBbkDiSwr2cBpq9HVdznrp7Hstyzzv8Z/D30yRLr5Mfwvfo6uPusDg+zf+ZY69o+DWU7B4SG5HnJy9lBuXY4odqOotHMhEPDXEvLSw4NJ9td0fMzodAw976GI/owLjS8yr70cLRjZZ776GAgtbjIUIu5++513TN3VE448zXbPqHvzDZuvwh4NMZ+9LuzWRZ/b7YhQkOG8EK7mzk0nMweKvW5ZqxrlfA8phoa/spx7QKP8T71TEPjpOdn9KHhX39vZv8wPqBlVSMNDdO+qA1gr63i7qVmdmby++hkILW4G9ozNYp7MPfvgsQzX/MEUm/PWE0eWo7kd2hbq73nZ/6X2YYIDSn5hfZLO0stV7j1OJG6bmau60nCc5lqaHg7gNDwpuTjQ2g4lC0usu81FKHjrNAwKJS5GfLSw1GOlXlm3BT/xXW2pdrH1Z77rrf78Pc12E17LlkavNy/HYRd/7Q0eeXQeqPp5p3jSwPbGCA0ZI0XvZzNOmr7E6jpz8w13RF6RD/vx9KX/z7/cyk0FBoKDSMc+3F1kn+24fCXBw44NLwUGhZ2L2UIcj/bg62QZyrP3nojasg0tLEh20zDfA2GvmYOX0e//UNd02ZPyBwNbT42uSamCA3p90XvpVEPQHWAl7med0KPwV7Di8DjPB9qSCM0FBoWFRrmerl6vfQnUvdZoeGKx52jScBAOvUWcy/lmO07gGXjyTvOLkxi1c9692HGfQ0zjBF5ZtF3975rAy+z32Jcv3r24WLG6GNvwaHfZUJDaF/0cu9rOOrZhrP8HalPEp/PFEPDk0zX7izwOHOENLvGSc/P2EPD5vhzzHT6b4hxN+RlysMNDTMs20u89HCEoeG5F+BmWedZhjpcuOfeDW2PRhUa5mj09SqMT34PD/zLtqDwsI+Zh8/PhLFBaAh97Gs42tmGPcwyTD4TbKKhYa4l+/uBx5mjQ/fcOOn5mURomLuL8juzNISGGULDHLO5dKctM6wpvAFIni8v7FP2wX2YbauKLOdT75GZI2h68SVJ0j35nv/uSd+f36vtrEvoFzW/NDYIDeH3y96vHsKu/RHWMfcswyrDOU0qNMwcom8HHmuOkPqXMdLzM4nQcPFy9dBLcHhcDfL34WBDwxz7xn27MnaWdk0Kvy4Zn6dd99w71yDf8uTHfOeUZRbv+Z9QS+CdI9w+yvgZZfR7SCI0ZPmXvf0eQsNRzTZ8/mevhxoeZjivHKHHfSl752WavVd7GkjA6cOC0HASoeHig3iWpYFvvkAOcd+gwYaGeWbfGDtXvS55lt4Vu59fPTs1y/P0198b7rd3A7aLTGP+bcb7ajfL77B6XE09U3OCS5M/uK4HUx8zERqS/2WvDh8exxh6Zarfdg/1q8OiLxnObSvT+dzmOJ8lzvcy0/neRDre2wzH6gVDaDiN0PCvvzd722y83rB+YPsbDjU0bMOBW+FMcYFNlWkvuZ0Jn79uqB8H13dZ7sHMM15nabsZ/5kxn3R5/cSXJr8bHI5kKT1CQ4bzcnveQ2hYB19fB163jXamXO7a/ch0flsZz+my52uZc0/Kn5GOOcfWAnNjpNBwCqFh+3J10WtwOKCgadCh4fNLu30NSwvNci0NLW8f0XxfWAx3D9UM49l2trH+uDrJ+2xlaP6UPJi0NLm3Lxt0UBYawosXvt0egq9/2hl6gxyM2hmaVU9h6+YIQ8NoYdoAgvOjSMd8mOFYfbsrNJxOaJhrb7UPZlMM5QP6oEPDPF1FjZ2rXZP9TM/YUx3STTIwjfgsLZa9VqfP/3dvyF3gX5zPScbQ8GCUY7Xgqo/fw1+nOjsboSH9vfQ99BQc3g1xCWQdcPVUrx8Zz3Grh/M76eFabs7ydhH/OrDrs2eMFBpOITRsziXHLLTP9ocaQLOCQYeG+Y7d2FlisJF5pteH513PMsyzN1m04OU/HVzrJlKFd6b+PLjNsGVBjwFYtuX/aYL+G2Pkh8/jU+LxsvjPIwgNyfvSN+8pBPunnbE3mG8qM70g9zrLsMfQMPt+lxn3Mvw9u/ZLxGPPsTz+whgpNJxMaJjzJf6j2VCFv4QPOTRsX6LvMxy/sXO1l9+HTM/XQymz42YDa77xyYzQ2yEGDJm70d72dI6Hgw0NCwr5Jxl4m+UpNIRXL31fepxt+E+7P9tG4TWq9zC86LFGPzKf71aP53qU6RwPM5/XeeTjz7WsenAvAkJDoeH6L1i9dVJ+9bL0fByFLv0bQWh4nun4jZ3lBWhFBBFZt0OItMfmUs9NPVu7sCXg757PYj/Jx3zjej/7SmbqGi+06mfcTDuLdATbDyA0JP6L30mPIdE/7aypIpuj1MfVc6iadZZhAaFh8kCkPb+nzOe0G/kccoWeD0PcRkBoKDQMeMF6KOKFabH076i4+3vooWGumTf19dNJedlrsp/12eox0K3viaxLYiOFeEt3GG4CqusfJd/77Th/mfme2+8xXDofXGBoafIy1zXlOPKoxkJDeOvFb7OHEOUtZyUtV+5xOfJL85GGHlln5r04t3rW6G3mc3lIcB5fMj6zv4yTQsMphIZZQ6Xl3Ze0ZHnwoWHOmTffroydy4dSOcP6x172l6vvvZx7zH27uox4fZ5WrnHdZKSw8LCXwLDnZfGL5jUDCw1LmRFccDOQpDNlNfQSGsIHL38/CgiKfjdI+dpzLeqX/psCanHb4/mXcC9cRt4HsI/AsAnDE12n84znMDdOCg2nEBq2L8h9N0V580N8n7NV2pfPvWwv3AmD0qwzbyItDx39uJmzi+3vffgyBznZx5VI40XQFwXNPq3XP0tYatrM8swdGPa4NPnVmHc/qNCwx2XubbB8U+oXd5mC4MrvJaEhvPfyt9nzMtzX6i7F25lr8LXdY7GUGuxMPDSMFiK353Q3putYL3ke896aEa73fLZotvTYOh/R8yM0TP1yWepLVhMeXp/l6tDbvEDVsy+XXZ44hNAw+8yb/kODlQKi42rezoh7bJ1neObyNyKqg4Hv1XaWZyj/EtHbqPdEpJmPz/fWQR+z7tox7GFqAdi/oWF1OqDQsNfA6t0vMAr44u7fY7z+mfganA/ldxZCQ/p5ATwsKCz6rZ7xd5BqX786mGzDhRJmFhYR0BQWGr5cur6xxrl8aa9vX4H4XeJrdZ/5fC5K2kLgjS8+jj55lk9G8vwIDdO/OJS/pGsRslw0XUAjzuRpApzFMu2LHpenHSa9vvlD4Ytim9ssrvfRq9k12ZcL9nS/PSYOqPd6+QIiYrgRvXFLPW7Vsy4zBIjNz8i5JPyN576YLwOGszS51718PxwHf3/Z0Oe+qPXzmPwLlqqIz5gIDSk7OKwKDIxezjr72QYDu6u+TLahwm773//ocebZMuf5pcd7oMTQ8J92D78fy8xALSAszBLutOfYx/25X8h49TsoXHbcuh3J8yM0zHF+w5qdsVhG1Sy/q4+7VYefrTeCot//7qD5s/V/m3tGYV+h4WI23eSWmP8nKFw+UEk/dvYb1J9HDt43Zv01oLgdzHVZhB9VO1btx5iV93wdvzazsftualWfW4aZrCuE8tUgfo/1PDNztvxegVWuGf+dcSXH7+iC7luEhpQbGm7PymiKstLef21o8J7HAZ3LY58vyYWHhi89tEvJT184a693SWHwVuJr9aXHYPS2j/Dw+Z+99lqvuz/llxE8P0LDHOfXx6b55AkN++2UfdtHeNjMeqvDlHW7bmaYKTnL2Vn4/SDgaJ3gonmh73uG7uLZ2Yt83+QNc+vn8r9ffsxffgHywt7i3zf3dWGhWFmztQps8lXe0uT6y5R1vqzL0PSnmS2aIzDUuVoWpAgU9iLI2/YLuP5DCA2H4Eem6zXv+Tzv2jFjL3Yg184k3H8RCA8hyBUajiQ0/BMu9R9kCA3TvETPez3H5gWwCUT2Ygdy7UzC/ahhSoZmFr1fk/8EiHX9mmt09EZgtf9vsNVzUPhH/P0zB9l9t2/1s13YdgRZO8cPdWlyaNOfBMvu61l/zR6Gua5dAZ2rERoyrODwQuiTXSlBgNAwzkzIzUzXq8/Zhm+5aZeR14FSvR/pbmvnxTF/ffG//3bS/jeXs7TbJAgNhYarBjDlNkYRGoa+RD8UFDTc1KFPG1Id/Amnvlc7L14gv74RXp38WV6ecrZVjtBwcU3KWCIvqBIarj9+7Rb5btffkvnlQre+lybX42usc1kEiPN174XmS4k+VjoU0OkcoSHDCg03ClvmOXZFbJYsNIxmnvmaHaj50nZH8PwIDfuYgVBSwCQ0jBWIHKh3WUFIM7Om9BlRJQYuifYha66HGq/wnFz/LPbdruQA+NvV5ajrs2iwUnWW3f/H1fmnjVgSz66WgaAIrPvy+yhkSO5yVlBHWqFhsNuertsvtRcaCg2Thhl1F9FHL8bjCQ2b61rPCFHzomZPFbZMeQjXZp5w3BtO593+g6+H1HvbBV/PUmfNZxzzBxmqTuQaIDRkuMHhrpBhOoGh0LD8UOqD67ZZ2DJloaHQcFSh4YsX6DJfury8rHdN6/3/zCItKjRsw40LdV8qqLpM/HxsqPPSz8h+8e91ixlt5c2ULSBsnXhoaJYhQkOCXwzrBgdD66gsMBQa9uFXz9du3zUQGgoNk4dMGz0vIRIaxn9Z3Ff3wkLDxXMmzP18ZlvyfeBmZliPZpZWkTNHC1ia3NZmZ6LjyJO9DBEaIjgUGAoN86ifj80Crt/ctRAaCg0Tn/uiq7KZUOlewPeyX1NLYosKDdtrognHRy/6ufaZ9CXJKALDFyFwpX4DCVTz1F/HZISGRH1BrDuePgodgv0o/DoLDddzUNA1/Ol6vOur0FBoGC/UuD7z0pzAi47Bma/nT/V/95p87eWalLicsozAcM84J/BaM4w/LOpeLmQfyEmGht+ubmQcCA1J8ZK4M7Nv2rrqwHV/ANc4R+hxP7Jre1bgdbz0zL1pcwTPj9CwpGfte/XVMsqoLzEPvV7Pb1eXrsMbMiyDFVqVFxj+GePUfhSBYXM965ny5XQovyiqNtPq3P6Yqus6QkOoXxQ3BBIru3s2iIE5V+jx7HAkS96L+sDz6lqeefY6nsby/AgNC3vWFs0Czr1Ix3gRv/7Z+/UUUv0nqOr/mlQHE3uhLyIwfBGmWKI8gsDwz/Us5ffV83NdWF1uJxQY7sg1EBqS44VxPpLQJ7V6qejGgK5rttBjtujOPeQl78XuTfnieu4PvMYx3QoNhYaJg43DmaYBQbMMS9mQvW2O4lou3BZxTRYz3qZ5Tb5d3fX5kl/PSPI8tGFLYUHXmuNb//uF1iH4//1d1GfoSWxRsfg9KzBEaEjWl8addhadQOKNgCL1/mljCD3qGZgDXa5cfGD4qsae0wxhm9Bw2qFhU5d61mE9U23qs6LWm0W1W9S1rIOSOqyZ/PWpihhz2muyNaHZQH+WcJaw75vnobkOvTe8i3Y951f3PY/5v4qryfiX4t/rlIzQkD5fHutZh/Y6/HfvwsF2ouoj9Kj3mWtD1sEsSR5KYPiixhvtcU/5udwc4/MjNCw2PNwsZglY+YHhZal7K7VLzy8mPbOqsKCkuSZT2HtyEaTPC/xSZFoNg+qZWT0tC08bGvbcZKjQGZuz0rpLx52tLDBEaEgRocTZhJcs1+d9niOYGGPoUYdwAwi16uDpcODX93SCz2Z93fbG/PwIDQt+5r5XO5prvD/TpLTZhcW+YPe3FLPYsKS9JmNdMlv0rLZ2RtYUZnyel9LdN8E13Opx7C9uafKrL/zGc283Xz5cn5Vab4SGTDc83Gz38ZtKePjYhqWjWLLQd+jx/M9BobNWL0Z2jacy6zBrkC80FBp+/IJ2/WPEIceybmfH1dEQl/m1L9lTmXV4PoRrtJgJOqLwsJ4NNKBZbe3en+ObmbW4DoP4QiMseO9pNnyBS5M7dak7TI+jIValQzJCQ0oPJuqZh0cj3kvtoV2WvTGy69Z76NHeO+cFXefRLUtp67zzrBrhs3n/7EdfgVOG0LuI+/H5n6+pn71RPnfNy0h1OKm9werOq8fVfCxLo5rZo+NcwnbfBNsDvE4jCA9vh9xgo515eD7ovVwXx34+hbDwxXOz2cs1G0gw3jYAGt4XRYsl9YfyCISGDC2c2Hv2azb8Lq4P7SzK0X6gKGmmVBtK9Bk6n48tFH6nztvtfT3k5/O+lGcz8Qy8+jyLWWIyS7sX6enon70mfGqCjrEt87ttwqfj6mBMjQPefKFc7O825Flu9805jCQoGWB4WI1pv7zFnofV0aC+FFnMKjwa6zLkT69ZPU7n3vZgYEtlm1nm9ThZeii+uJfnU72XERoyroBit50FNJQZiKMPCksNDdvjqfc6PJnl67D82F7vyU3nb2t90Ab8pTc2emyXWB+VuIS1rWPMe7bIPVPbWcE/Im9HUddtPrnnr57xUb+4DnP5a9UENcfV3hT3TWpnjx40ezTWMzzKvlaPzT1W32sj3hT/RXhV3n6iixf7k7E3JViE6tW80C9FHpsQ6Hv11btZO1M039h1PpKxvowAsbluzQxxy5ARGjLqgOqwfeksZZnkXftyfjjR4KjYPdnamYfnCWbEPbUB1GFJs7gKuRcO2hD1tufn8qY9jkE9l+0szsN29uFFO84t66wNzAfxhUX7hdBJe9yrnOdFW59JjrkfvsQtQsTzYl66Fy+VVTsz7WRKy/hWvHZbzYvlYhbibc/X7Ka9XodTfalsA8TDJiztI9RdhAtVGxRO9RpsNtegr2diMfPxvA3Ld4xTb1yjXDMOxzSztv6ibLH3Ye57+ra9n0e5fRJCQ1j2JfugfYn81b5UppiV+PuF9bS1N4WlqEMODV8c4+8ZcRdrhlkP7bX/Mdb9ChPdGxvtc3L64vmJGeA+tn/nZfsz9ut9F9Ue/sxw2F0Eic3yy/MmiIi9DPB3wFHPpFjMHtwXDkYJrfba61a1Yi6ffWyv2eWfayYY+SzAWlyPxUzEKnI4VQmolhzPmvu1CRKrxZ6aEcew37OfLdVc7prkmeU+uKXJK4aIu+3s2vP2nr6NGBDO/S5GaAjLh0W7rxy+CP5ee/1nLUP4vMbFh4YfHPvXN+6JE9e/t+fzvWfz8NWfG+2eZ5A5mNpd2QSXExcboPzrcBEwvtLM1Hrx50a8X2TPL/0fX4ffYZTnKM01qPd67V6Hozfqf6D+UX9/5GmIcnz90zjf6t7P8//8bhZ208d9qgjAUoPFgENDAABghc/+9dL5PEuTd9UbCh4LFAFYarAQGgIAwDQ++8daGv7JnrhqDYWPBYoALDVYCA0BAGD8n/ub/T1zzDKc7tJkGMx4oAjAUoOF0BAAAMb/uT9PAxRLk2EI44EiAEsNFkJDAAAY92f+XA1Q5lf36g0DGBMUAVhqsBAaAgDAuD/zZ2uAcn2m3jCAMUERgKUGC6EhAACM+zN/jgYote/VV/WGAYwJigAsNVgIDQEAYLyf93M1QLE0GYYzLigCsNRgITQEAIDxft7P1gDF0mQYzLigCMBSg4XQEAAAxvlZP18DlHpp8raaw0DGBkUAlhoshIYAADDOz/q5GqB8u7pTbxjQ2KAIwFKDhdAQAADG+Vk/VwOUeeXzPgxpbFAEYKnBQmgIsM7Yuflsrx7fXqhe+NH+b/NnO2rW6zU4fGbJHDC9cfK42s3UAMXSZBja+KAIwJIvXUJDgOXHzDqoulhjHKxDrA017PUaXLoGwKTGy29XvzKFhrfqDQMbHxQBBvCg1hsTH1d79XT+N9X/7vnPJH75EhoCLD9mXgSMhftq6BoAZHvPyNUA5bg6UXMY2BihCFDwA/q92n7+BXu+wsbCv1KFh0JDgKXHy43AsXBPHXu/BrvqCExivDyu5hmXJm+pOQxsjFAEKPUX+PVZSR3JhIYAS4+Xu4Fj4YY6ugYAWcbL53eHTF2Tb9QbBjhGKAIU+GAeV/uBU/93E7yACQ0Blhsv5wHj4L0a9n4N7tQQmMg7R74GKJYmwzDHCUWAwh7K79VW8L4ix9c/E7yACQ0BlhsvfwWMgxdq2Ps1+KWGwCTGynwNUP5Jvf86kGicUAQo7KFsmpsE/2I+T/ACJjQEWG68vAsYB83E6P8azNUQGP04mbMByvyqUnMY6FihCFDYQxllX5HqNMELmNAQ4POxUhOU4V+DXXUERj9W5myAclwdqTkMdKxQBCjogayXJhf6i1loCLDUWKkBh2sAUP5YmasBiqXJMOyxQhGgoAcy1jd+x1X0mSpCQ4ClxkpNUIZ9DTRBAabwzpGvAYqlyTDs8UIRoKAHMsZ+hvXeJP/395cEL2FCQ4DPx0pNUIZ9DTRBAcY/TuZsgHJcHao5DHi8UAQo6IGcX51HCA2TvPAIDZnMc/h/f39pZv0uPlDffvLM3TbPrQ/E/DtWaoIy7GswV0Ng1GNkzgYo9c/56+8NdYcBjxmKAAU9kPOri1K/zRMaMqEP0ndrPnsCn6nfP5qgjOEa7KojMOpxMmcDlG9Xl2oOAx8zFAEKeiDDlwrcp1ia3L6ICQ0Z/zM4v/4R8Pydq+HE7x8NOFwDgPLfN+4sTQaWHjMUAQp6II+rvaDp/9+rnYQvYkJDxv8Mfr4c2UxDPhonNUEZ9jXQBAUY9xj5vdrOOMvQ0mQYw7ihCFBcaHG/1i/l4+og8YvYZobQUOhCf89evZdhYV3LGdg9pAnK0K+BJijAuMfI4+owX9fk6x9qDiMYNxQBCgwujquT51+2j0uEhQ/Nn/3r781ML2M3iUPDbfcAvT1736uvQR+OfZvuHtIEZejXYK6GwKjHyOPqIFNoeOtzEYxk3FAEKPoX++7iG8HqtOO42k+5FPmDl7F6r6inRIGh/eDo+3kL2Rjc0tKp3z+aoIzhGuyqIzDqcfKvvzeCtmJZco9ngSGMaNxQBGDFl7J6mfJZvf9WhKCwDiAvnu2rLb3f2/WH3PU/IFtaamzUgMM1ACh/rGyCw+sfzfZGcULCenVU1Xz5+r3aUmMY2ZihCEDAC1rdHGWvbmDSqsPE6h3nL/7coRkdFHc/a4JC2HioCcqwr4EmKMD0xs3v1c7zZ5ijP6uZvl1dNgHgW46vz16sejpqVkRl2iIJ6HGcUAQAJv/LUBMUQu8hTVCGfg00QQEAeP35ShEAmPwvQ01QCL2HNEEZ+jWYqyEAwKvPV4oAwOR/GWqCQsj9M5t90QRl8NdgVx0BAF59xlIEACb/y1ATFELun9nsqwYcrgEAwOg+YykCQOaB93u1NTuuDhcbSdebSv/ZZPr8xQbTu/U+e+qV6ZpMoAlKc0/Vft9jjRcbnH+7+tX5dyO+BxfL0Zvn7+d/nr/n81757xtgE5R2Zt5u25jqd5Oqy1eNq07aP7NV/DUNuwa3K/6szWcHbY1+vGj4dfmilnXdjOEAwLA/NysCFPRA/vX3Zt1Q4d2X+vn1j/al9vD5pXc74gv0VidM6AZZ3UDhuDqI+bPffaH/HaotOrn9ereD26Jz20ak+m+8qsPPF+d/+SJM2Vv1Zy5q/Pz3fbu6WzqM+nb11Jx7XfPI4U3bLe/jGr+85xJf8/9cg2a58HvH1R7T4lrsL3MtXvy9bwtrgnLy4d/92vO9kKWWi+fopK3bPwFum5pnbvbyn/Gwey/8DvhO6vNc+u9b/D0PS5zzY30PrhAinQ+hCUobdh3VP3ON46wDsYMVf95uGzzWIdrZi3CtehGu1eHbZoRzC7kGP5f4+zfaYHLVfRPrc/2a89kBAIj2+VERoJCwcH51scbLfBUamDVhQB1OFTC7qm1GsfqMrzoECAy1mvBpuTDh98+8q6/bUuHcIuz4J0j985YMRz4Nz9a51+p7ZI0ZWCteg5MmrFnn+n9wLZprsM7fm9JxdZgoKNxa+5le9j48ro5Sz0BcBNorXLPj63dDn6Yzdv0lw6o1qf/8kuNrPVOt5CYoz//stKHfPxE8fHbMz/9srxiuPdbHGHiOIdfg8IO/d7MNPB8D63YxhBmbAABCQyjtQWxmsa0dPhys/XProKyg5g/PdbjpI4RpQoVVAsN/f+buJwHoRfyw6fosJLBpZ68Wt3dfHfoGhnDz9885YOlxKt+ufsUPC4P2ZVzVY7Lgsw62Vw8979+py85Ks3tfW2JWaMlNUNrA61eksPC1mzocfOfn3uQMTyNcg613gtbzBHU787kHABAaAsu/IIeFJbvrh4bN8tt1Q4+bqHWog7uwOqz94t0GfOsECtvvhDdV4sBm7fAuKMj8dvWQ7DkICzPr63/09j0euPQ4nSri+HHeY/h5ucyM29XGpXpJ8urH8c61D5thukRAX2oDjnbp71OiwPC3+u8/evVzN9b8u/YDzjXkGty/8fcdJa7dD59/AAChIZDmBfmlgOXJYTOwrqO+9Kwd3MWoQ7N/3urLdd/4ew6yLYP9YDlmstAwsM4fHlfIbLAPwvMoS8MThW1Rxo51ZsimmHX4vdqJOCaerDMD943aPAVeo6VC8tKaoLSh3UXisPC10xc/fzfWbL8VznkeI8BrZyz+ylSzuc9AAIDQEIj/gvxif7G1f+5iSe5TQEhzGLkO84AX/KAX77Vmar2aaRk0azPy7LqP6xx4nBH2VXwjMN4KDHee3psRVtxehoGhb3NOfc8ufH9f0a0o98M6wfZxtf9iLDnMORu0pCYoa+wlGNNhQID3GHjeIddgr8fa7fscBAAIDYG4L8gR9kVrm0OEhEc7kesQEoJcBP7sNWZc/jvTsrcAZ9EAZCNjOBu0h+YHx3QSWIdf79zj20UGhp/swbhEwHpf6Hndx2iQstb5tYFlxMBw6dnUpTRBaffge+wpMPy9VHlv3c7Mged+G3DMX9rjfuypZps+CwEAQkMg3gtyYPDQvlwfhczsSlCH24A6rL+B/rp7Kbaz/HqZYRiwTDx4Ofy8Oo1+7UMa4LyaZfbqXA+LDQ3XCN3bELTUwDDKtgVr7vH6+Kc+MbtGLzGTt5QmKPXS3p4Dw5ddkB9yNgcJvAYXBYStJz4LAQBCQyDWC/Kn+7gtFx5pgtIGDes2QdkJnrUXZ7bh0yqNKOo/m2JWX0D9Q5cmP7wfRAc2V0m4B+Aaddopdql1xFnI6zZBaceQ28jn8elS/BKaoLQdku8LCAx7WaYbeA3mBYStdcj6xWciAEBoCIS/IEdqShE2u0sTlOaYY85qyri3YWD4dBv5GTgKO/frsyTBeFrnKz4f24MJDEO3TVhrqfr1j/o+iB3sLtk5udcmKO0su9uBB4Z9NUGplwbfFHL+hz4TAQBCQyDCC/KfF/O7oJ9dVhOUkOAofxOUekZT6JLauKoVr33IsT9GvfbPxx44G2z7g/Mss3PyCjPxmtnIoZ2lX4bdTU2q03qW8guH7azMi2g/Z829Ddc6hhTh8JLBZ99NUJ7/+TmCwLCvJihPBdXgxmciAEBoCIS/IMeYzVNaE5Swl/4emqAUOOtrhdmWwSHLCsuhPw3EwoKdm6TnWcIswzhB3mMTzC8zc64eG+rlvuGzXw/WfB7v1wop48/eXeqLkT6boDz/szuCwLDPJiilsUQZABAaAoEvyP++1M7XD+kCGkRoghI6w2u/DlSaJZUxZyu+0wzknesftiR4ib3ekt+HSwQ7zSy6usN0OYHhxWrhboQ9M+t77IPZmOnCytUb5gSHyPGe14cllyb31gSl/dmx9zGsnp22y3332///Z4ZQrq8mKOvW56heTtz+/zGvwVefiwAAoSEQ5wU5pAlKHVglmt21VnAXtlQ6fxOUFZZON8HPB0FRtGYqK4Snwee94h6KHwRTVcB9uPYS2D8/P2z58kXS8WEx4+8pd3DXGZ/CAtfzlX9mcGfvpd229171Zo0/2CfzVWDVWxOUyMuSfz3b/uTn7bR/bkxNUJZtUlKHg5sfHMNppJ8199kIABAaAnFekAMCk6AZbsfXP6OHI5mW5SYL7N7bZ23JaxQ8427F5erhgfX1WYTrvhV4zufBxxCyV2DADNfkz2hgmP7nGNbb73OtfTbb5+Ak2fNYB6D18/7GM9mOxefN/VAHycs+tz01QWkDvCh7Ca4a2LUhXexuw1sBtZgnDAx/LRvsBu5t+efn+WwEAAgNgRgvyEEdbANn9h1FrsPQmqAss3/c/hrHchv4c29XvAdCZpGFN3AI7XYbMNM2TnAaHsp98Ezs9x3qBj+ba+y5OovViCUgwF8xsOqlCUr930YIqO7XDevq/+7Z3cCboHwWph6ueBxfItTkro/PIgAAQkMo8eELe0Fee5ZVcU1QQpZK99ME5ePAcM36RJhltVKAGtjs4jb4uoeFlvfBP7/e77CnGa4fHtdiuf5DQFh3EysgCwoN1wgug/Z4fW8Je+RO76+CouxNUCLNMnwImd3XHsdGpP38SmuCUgeGO2sey1lokOuzEQAgNATCX5B7bIISe8ZO2FLpopqg3NdLbtc+nr/+3gwNSFa7/4L2tQxqhhM+ky58aXDg0vRkL/dBx1U/n2s0PXn/WAK6T684Izl6E5RFYJhuNmh4A479NX9ujFmG+5FqUC9VfhpRE5S1A8MX9QgKc302AgCEhkD4C3JA99rAmX230WsxjiYojyGBYZRaPFsxnArtoLx+QBoyy7au0V9/b0aodXFNUIJnGUbcOiB3g6LoTVASBoaRAqLNNX7mZoSQ7jJyHQ5H0gTlKSQwfHFMQfs9+nwEAAgNAU1Q/g3uxtEEZY09DN8J0+6zhYbhHZTXCmWaGZVB4ej1jyi1LrAJSmCQ+xhzFnDQfp9rzEiO/DweJR/DwxpwPKz5M08iBGPbkeuwOZImKIeR6nEvNAQAhIZAaDigCcr/grsGF9IEpTqNGNRU2ULD8EYg8zWveVg4FGFPzVKboITtsRkpTG2W7QfNRl7ry4XAmZ/JZ4G+EQ5lb4ISodHGWYI67I+gCcp5xHpUQkMAQGgIhIYDmqD8bxRNUG4j3xfZQsM2qHnIHVIFzfCLFAiV2AQlePZxjDC1rkvY9Vn7WKL83Pp+TtSg5o1wKKQBx+kaP283QvOTFB2kTwfeBOU+Zl2EhgCA0BCIEQ5pgvK/ETRBiR6iZg4Nw/YWXPmFP8KS6N1IAV1xTVACZ9qtHV7Xe1O248JtlGfi+Zle+RhiNUGJtE3AEsFQ9iYoEWbVnSWqxeXAm6DsR66H0BAAEBoCAQ+dJigvQ6vHwTZB+Xb1K0E97rOGhsfXZwHnf7fG+V2UcP+V1gQlQtORo6XGnXom4XF1UC+pb5fn30cJ6wIbAgXP/FwzxA4IhvpognIf+DO3E9XisY/ALlITlMsE9bgTGgIAQkNg/YeuzyYoITPZ4jdB2R5sE5RIHXzfuD6PAcd1v0YdDnLdi8FL4yPOIiutCUrwdSjFmtcoShOUyLN+PwmGsjZBqZuFBIZjN4nqEHpcWz1dg2jdkt84rofcDXIAAISGMK7QsL8mKGEz+2I3QTnIGZC9qkNoE5TzJPdG2DFVK/+88D0ut1c4t5Cl0A+xlsaX2AQlYhOQHq3fiCXC+V9kHcMzN0GpO/sGBmRHieow5CYoF4lqUly4CwAgNIQhPXT9NUHZLmkmz6CboCSY1VTPXMwdZAbv7bjkzLIIswzn0epcWBOU4KXJfauPPfD6BDdByTjLsA2GcjdBOQ+cUbeRqA5nA26CspugHhuBx/Qr530MACA0hDJDw76aoBwEBQOxm6CE1aG/JihrNHrIEmbNq9O1fm7YUt35cud2/TNon7yI915pTVCi7K/Zn/uQPVab8w9vgnKbdfzupwnKfV/h3CfHFdL0o88mKLeJ6rFTYrMaAAChIQzlgdME5XdQstXXEtEIHXwPEoWG8z6OK2jm6xL7XDYzKINm0q2/7PXN4ymtCUqM/fx62b/w+meMWZfBYXnkbROWCIayNkGJsG9gqq7JG+0sxiE2QZknqkmRy8gBAISGMJQHThOUOEFJn01QIi9PjRRm/bNO19pFPQI6KC+xj2JQWJ2g4UxxTVCGtJ/hIvw9X/deG9Lz+EEwlLsJymHumY1LHtd+4HFt9XQNgn72J8cVus/ijs9JAIDQEKYdGmqCsghKLgOCi7vAOpz3dQ2ShVkBy2YDl61/GIIUN8uwxCYoIcv083lswt8UHcPDQtPb7GN4/iYo5zlnNq5wXD8H2gTlPuG9cdtXTQAAhIYwhgdOE5Twpcnfrn4FXoPbUkKsaGFW2L0R1qTkg9mvxc0yLK0JSnjzm5Rum1mogXsWLhEa3pX2PCYMhk7X+Hl3hQZkDwNtgvIzUT1Cm6Bc5L6XAQCEhlDaA/ft6mHqTVAC91YMq0OmbsHZw6zj6rC3mrwTKJc2y7Ctc2lNULYLCQjvm+0L6pCwHisyLfmNMPNzP+v43U8TlOKCqAhLk8+GdA2WPK69EvdZBAAQGsJQHrbQWUUjaILSBFQhwekiKNgNCGnCmqCk2s8wNEgN3GMuaPblOw1YgmbVJphl2BxTeU1QAjtmN/uUrqC+z6rTRr2/auJZhBmex82sxzu8Jiiniepw2Vdwl/sarHBcZ/YzBACEhkBIQLA/9SYoUTrF9tUEJXAvxXePKTxIDZ4BFxSmvdEcJLjhT6pl4KU1QUm4n+RAxsT5kM6/hyYou4FB1EmCGmwHHlOfTVBKXa5tP0MAoNzP7IoAmR62ZoZPP7P9SmiCEmWWYZ9NUAL3UvwgONkPDFLPIxzDSay6tNf5LuDvSzPLsMQmKM/PVkCdboY/JgY1Jcq+B1wPTVCOAsO5w8Jq0HcTlFTLtUOXJtvPEAAo9zO7IkC2F+RhNkGJFJYEhqb9N0EJ2Evxw2MK6SQduJ/hi/BqP1Z4FT6bNNEsw8KaoAQ/E4lC7MxjYkBTouo0+/Hmb4JyWtL+ffUS2gizDPtsgnKS6L74ZT9DAGC0OYYiQKaHLawJytqz/YKXikaY9RXeFCM8uIvQ8CP6/m9RuucG7mfYHEdYsPwY8To/ptqnrrQmKItjuv7Z5wzTXsfDQpsSfRAM9dEEJTQ03I1cg9sIoWGfTVB2E9wXddfkp8Dj2vYZCQAo9nO7IkCGBy28CcpOT2HJP1HOP2yWZf9NUCJ3kP63LoGzLyMuUQ0K+9pZeIGNRpLN5oxwbBflHVP+mXZRz314TVBCZ9ltrvEzf5QSktVLnSMEhqFNUEKvwUaC+2IeeEy3Q36OAQChIRDjQQtb/hkUWAXOZrqPcO7zKIFhn01QInaQ/nM8MWZfRlia/Od4QpaK1gFQ6PLfeh/EBMHsi4DurrSALqhBUaJl3BnHxKE1QQnZX/BhzZ95XkJo2DY/eYoUGm71dA1uE9wTm4ENUGpHQ36OAQChIRAlHOixCcq3q5u+QsMmTIqxLLnvJigRO0i/OJ6LwHpEnf0Y2EH5KGx/urTLTSMshT1Iclxhz+agmycMsAnKz9yNLiKEhgcRzrteEnwXKTB87PEanCe4J0Kvz1OK2Y8AAFE/8ygCZHlBvugrsKqDvz5Cw3qvvcCfHbcJStBMs+fziNgII0LH5OhBZlAH5ZDu3AtV0ufve7XT17L4T+7Jy75C9ALGRE1QPv+ZZ313T47Q5KOUJig3dQAa8X7Yi1CPQe9LCgAIDYF44UAvTVDal/OgQGfNkCZuYBjaBOWvvzciHEOU2U3trLfw2gTsc/nmccUIMgs5lzfO7SDw+JI0Kgiabfft6mmw4+Ewm6CELM/dX/PnnvYZSkWYSRe7CcpTCSFdeyz3EeqxM9RnGAAQGgKxHrIem6C0wUTW7rxJAsPQJiih++1FbNQRuCwz2fLM9rrlDwwTLP3+b80DG84k6Jy9OK7rHyUeV/LroQlKrtDwcd3ZdQkCw76boJQ2+3LQ2wsAAEJDINZD1mcTlNDAchGUnawUzoXMqiyzCcra9UgSGNb3xIpB7grH95g1MKzvlQwBUHDdU+1pGBpmJl7WnXBM1ARluZ97ECGc2lvxZ9bNPS4TBIZ9NkF5vY/gYcBxxApTd3w+AgCEhkC/TVBizB5bBDsfBnaL5YbBAUi6JighTT7edr5KmNt2Sr6M87PTdc0NbMzRa/fnT87rssRwLnjZdM4a1mNJE/Y9P+eBy7U1QVn658aYXXe/bFjX7tP3kCgw7LMJylt+rPjzY4apZhkCAEJD4M8Lcm9NUNqfHyW0e2t222ImY3WabHZhGU1Q3m+OUocoHwSqzV6Ki87CsWbwPaacmVffbxlDwyrjM1iFH+/HzTfaJbfnzc9actlwlL0262cv1czTRVB49J/6hT6PmqAs+3O/RAqpPgwO27Aw1ezCEpqgvOeuncH45YOfu/Fs3i71jjXTcWuon2cAAKEhEPsh67EJSvuCHnPJ6W0TjNQzt9IEcaU2Qfmsc3A1O74+W8y2/O3qollKPKCZeVGXcX++xHo74zN4GS0obvYh7Fzn8zee8fMVns/bOMf2fFyBHb7bkPCgCY8/fr7XDnw1QVn5599FDMl+tfsk/vYrUkOPITRB+Wzvx6rtVv2yPpcJfu6Rz0UAgNAQ+B1Y9doEJW4oESksWifELKEJSv+1u0x+vx5Xe5man5xlfQ7jNJ9JMotyETbHDLCvfyyCv+f7/oNZqc3MyPrPHFcnzfL9Vb7cCLgXNUFZ+eefZQr1VpkpN9QmKH0b5B6kAIDQEEgXwvTWBKXHwOTDpYVrhYalNEHpz2OWhiExGucsM1svwn292jMQ2KU44b579RcDA7wfzwOexyNNUFb6+buFBV/rNgLZ6ukalBS2bvlcBAAIDYEXYUV/TVCKC83qJhvrBCTlNUHJL1H33ndC5sfE57Kb/TnM/QysuBdp0L6n/dyP8/WvRdC+mZNpgvLqGB5KmSn37GQETVD6cOAzEQAgNATihQERmqA0x/C92i4gaLhv90vbXyM0LLEJyiACmjXrdVNKmBbtnHLP5ltx78nBzTYM2S5AE5R1juFHAaHXfdsU5GIkTVByso8hACA0BN4MYHptgpIlCFpmaW3b9GKtmZdlN0FJ7PpH9ns2XQfl7MuSX4VVj9mu2xrdjAcz2zBg24SmCUpIc6CJNUF5dRx9zjasG4Vst8eyTuOUPpqg3BcSGJ76LAQACA2BtwKr3pug/BsEVYc9hQy3LwOUtYKR6TZBOe/lvk21lLeHZcndUC7bvob3a44XG02wWn5ouPbM3+AZlRNrgvLqWOY9hV63v/fia2caDqUJyk4BMxTPfRYCAISGwHvhS+9NUP4cSz3DJ3cX5XovwVfnsFYo0lcTlMXszNt+wpn8yzBfBDtf4weGebslvxPKbQbNcsswO7RZwl96cBgw208TlKBj+dJDCFYvRd54cQx7A2mC8tj+t1s9zjic+xwEAAgNgfcfrgKaoLwKJLbzhCZvLylea6lwr01Qrn+0x5wzOHzMvQQzynX65Br2uSy5G1pdn2VYmvw18DndCdrWIO0XAZeB9dcEJex4NjOGYGdv/PyhNEG5fPHf77TLq3Mu5d71GQgAEBoCn4SG/TdB+e9Le3WYNDisw453Qq+1ZrD12QSl7VjchGh59oS8/b33Y+/3bszQKjBES3Bulwnv/5sox7gIbs8LCw1vQ5cHa4IS5ZhSh2AP73X6HVATlLNXf8fXZ3eZuktv+fwDAAgNgbTBS8QmKG8Eh7sJmkLcN0sPP5hRtta+igF1CJ4x92pPyedjOUkUuN73PbswXbDW3zLrd89tsVT/PEFg+BRzH9L2ntsrYNbhp8/2CqHhfaznMVNoGDKrbyfhcW0/u0nQHfmwXgb9wc+tci/TXfMa7Lzx99TLu88S7vtY1BgOACA0hJIfrIKaoLx5fIu9084jBCV3dRi45M/cXikAqZcyBgQVQU1Qvl09fXBdLyKGMfMi798YTUPq2ZmFLEv+IJC7i3QtL9bpmLzCsR4knSH5xjLkZj/Q2CHoOls2LGYvH/RyjyyajzzFmqmX4PhOArs7/9POwDtc8uedrvh3//gohFzhHJf9eXUtTj75+76uGX4G1Q4AYLDZhiJAkpf8YpqgfHicdXi42OdtlWWDj01IssbsuGaWVx3m1bP2mgDhDfW/j9AlNbQJyod/dx0eLv7+VZdb3jbnWMgy5HfPL2gvyD+z7oo+xxfPwNcmQF8lQKyDrEWwdpSzo++f+25xfW4jhIP1zNmqCejrZ7J+9hKPPUuNAYtxYL+Ee6jtGLzbBmYf2a9nAPZ0fIcrBmF1sPlrndlx7fLow09qsfuygUqEc9z85GeetD/zywp/51YbCt+tMavwpI9rDQDQy+dhRYAED1ZYw4Wql2NeBBLvv8zXAUkPSwTXPp/AJigrhE5bbd0O3qnbXvPvE85Ei3vvrrGMPOPy+sTPwEZ7LfdeXcOT9n/fDenmneSY69BzcWzzD0O4l/digedBhHthEa7ttgHmW+HaPOWy6QHXbfuTuu21/34QYzgAQNTPSooACR6soKWE5e0DN9BrENIE5XCSNYvRYTuwwy4AAACFvCMqAkR+qOrldyHBy3G1p46B1yByE5QJ3bd3gYHhQ87lugAAACR8T1QEiPxQNctUA4IXywZjXIPdIewpWVTNQvcxFHgDAACM6z1RESDyQxXWoOBWDSNcg5AmKBO8BoH1+tPt2r0HAAAwondFRYCID1Rw+GI/wyjXIWTW3MTCr6aRRvg+hndTnJ0JAAAw6vdFRYBID1Pdfbje0y1sL70ttYxwLcKaoAyy8++a92y99+N9YGD4VAeP7jsAAICRvTMqAkR4kOo99EIDw29XN2oZ4VpogrJ8reZXF8HLks2OBQAAGOc7oyLAEg9K3Vn2uDp6dtgEhL/Nq9N6OWt48NLMcDtR6wjXKqQJSm0iy2yj7GM4v6rccwAAACN9b1QEWOJBiTIj65Mlnn/9vanWEa6VJiif1yjGPobzq0fL6QEAAEb87qgI8MlDUgcsKQPDxvUPtY50vUKaoMyvzkdfnxj7GC5mxh663wAAAEb8/qgI8MlDEmv58UcztswyjHe9wpqgzEdfnxizZr9d/XKvAQAAjPz9URHgk4cktMGJoCrftQpvgjLqLsBR9jGsn4fnOrvfAAAARv6OrQjwwQPyvdpKvJfhzVQab2S5XqFNUEYchgXX5t+Qe9e9BgAAMIF3bEWADx6QultyymXJGknEvl7zgAD3brR1+evvzSgzZo+vz9xnAAAAE3nHVgT44AE5vj5LuCx5X40jX6+QJigj3qfv+fyqCPfsrVmxAAAAE3rHVgT44AGJ0TTiv+HU0+y4OlDfBNdLE5Q37uHqNMo9+73ado8BAABM6B1bEeCDByTODK1uE4mRN9vo7VqFNkEZ4V59EfcxnLvHAAAAJvaerQjwwQNSNyqJFxpe1HvLqWuia/W9+qoJyqua1EuKw4PuS/cXAADABN+zFQE+eECOr39GCAur2fdqRz1TX6uAWXUj3M8weObl76DbPoYAAADTfM9WBPjgAWmCl+sfs/nV/aqzs2bH1YmwMPP1Oq6OVtzX8HbMS2/bejys3PCkDsvtYQgAADDtd2xFgCUfljpArGezHVeHi+YSLxxXe+2/21WrIq7VZns99l9dp5Pf12lKM+hm36utd+px+Kce9toEAADg5bukIgAAAAAALykCAAAAANChCAAAAABAhyIAAAAAAB2KAAAAAAB0KAIAAAAA0KEIAAAAAECHIgAAAAAAHYoAAAAAAHQoAgAAAADQoQgAAAAAQIciAAAAAAAdigAAAAAAdCgCAAAAANChCAAAAABAhyIAAAAAAB2KAAAAAAB0KAIAAAAA0KEIAAAAAECHIgAAAAAAHYoAAAAAAHQoAgAAAADQoQgAAAAAQIciAAAAAAAdigAAAAAAdCgCAAAAANChCAAAAABAhyIAAAAAAB2KAAAAAAB0KAIAAAAA0KEIAAAAAECHIgAAAAAAHYoAAAAAAHQoAgAAAADQoQgAAAAAQIciAAAAAAAdigAAAAAAdCgCAAAAANChCAAAAABAhyIAAAAAAB2KAAAAAAB0KAIAAAAA0KEIAAAAAECHIgAAAAAAHYoAAAAAAHQoAgAAAADQoQgAAAAAQIciAAAAAAAdigAAAAAAdCgCAAAAANChCAAAAABAhyIAAAAAAB2KAAAAAAB0KAIAAAAA0KEIAAAAAECHIgAAAAAAHYoAAAAAAHQoAgAAAADQoQgAAAAAQIciAAAAAAAdigAAAAAAdCgCAAAAANChCAAAAABAhyIAAAAAAB2KAAAAAAB0KAIAAAAA0KEIAAAAAECHIgAAAAAAHYoAAAAAAHQoAgAAAADQoQgAAAAAQIciAAAAAAAdigAAAAAAdCgCAAAAANChCAAAAABAhyIAAAAAAB2KAAAAAAB0KAIAAAAA0KEIAAAAAECHIgAAAAAAHYoAAAAAAHQoAgAAAADQoQgAAAAAQIciAAAAAAAdigAAAAAAdCgCAAAAANDx/6Dfdq4bkpSWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 4,
     "metadata": {
      "image/png": {
       "height": 100,
       "width": 150
      }
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Image(filename = fig_dir + \"ogs-jupyter-lab.png\", width=150, height=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d2ff19ce",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtsAAAJMCAIAAACkcf0EAAAABGdBTUEAALGPC/xhBQAACjFpQ0NQSUNDIFByb2ZpbGUAAEiJnZZ3VFPZFofPvTe9UJIQipTQa2hSAkgNvUiRLioxCRBKwJAAIjZEVHBEUZGmCDIo4ICjQ5GxIoqFAVGx6wQZRNRxcBQblklkrRnfvHnvzZvfH/d+a5+9z91n733WugCQ/IMFwkxYCYAMoVgU4efFiI2LZ2AHAQzwAANsAOBws7NCFvhGApkCfNiMbJkT+Be9ug4g+fsq0z+MwQD/n5S5WSIxAFCYjOfy+NlcGRfJOD1XnCW3T8mYtjRNzjBKziJZgjJWk3PyLFt89pllDznzMoQ8GctzzuJl8OTcJ+ONORK+jJFgGRfnCPi5Mr4mY4N0SYZAxm/ksRl8TjYAKJLcLuZzU2RsLWOSKDKCLeN5AOBIyV/w0i9YzM8Tyw/FzsxaLhIkp4gZJlxTho2TE4vhz89N54vFzDAON40j4jHYmRlZHOFyAGbP/FkUeW0ZsiI72Dg5ODBtLW2+KNR/Xfybkvd2ll6Ef+4ZRB/4w/ZXfpkNALCmZbXZ+odtaRUAXesBULv9h81gLwCKsr51Dn1xHrp8XlLE4ixnK6vc3FxLAZ9rKS/o7/qfDn9DX3zPUr7d7+VhePOTOJJ0MUNeN25meqZExMjO4nD5DOafh/gfB/51HhYR/CS+iC+URUTLpkwgTJa1W8gTiAWZQoZA+J+a+A/D/qTZuZaJ2vgR0JZYAqUhGkB+HgAoKhEgCXtkK9DvfQvGRwP5zYvRmZid+8+C/n1XuEz+yBYkf45jR0QyuBJRzuya/FoCNCAARUAD6kAb6AMTwAS2wBG4AA/gAwJBKIgEcWAx4IIUkAFEIBcUgLWgGJSCrWAnqAZ1oBE0gzZwGHSBY+A0OAcugctgBNwBUjAOnoAp8ArMQBCEhcgQFVKHdCBDyByyhViQG+QDBUMRUByUCCVDQkgCFUDroFKoHKqG6qFm6FvoKHQaugANQ7egUWgS+hV6ByMwCabBWrARbAWzYE84CI6EF8HJ8DI4Hy6Ct8CVcAN8EO6ET8OX4BFYCj+BpxGAEBE6ooswERbCRkKReCQJESGrkBKkAmlA2pAepB+5ikiRp8hbFAZFRTFQTJQLyh8VheKilqFWoTajqlEHUJ2oPtRV1ChqCvURTUZros3RzugAdCw6GZ2LLkZXoJvQHeiz6BH0OPoVBoOhY4wxjhh/TBwmFbMCsxmzG9OOOYUZxoxhprFYrDrWHOuKDcVysGJsMbYKexB7EnsFO459gyPidHC2OF9cPE6IK8RV4FpwJ3BXcBO4GbwS3hDvjA/F8/DL8WX4RnwPfgg/jp8hKBOMCa6ESEIqYS2hktBGOEu4S3hBJBL1iE7EcKKAuIZYSTxEPE8cJb4lUUhmJDYpgSQhbSHtJ50i3SK9IJPJRmQPcjxZTN5CbiafId8nv1GgKlgqBCjwFFYr1Ch0KlxReKaIVzRU9FRcrJivWKF4RHFI8akSXslIia3EUVqlVKN0VOmG0rQyVdlGOVQ5Q3mzcovyBeVHFCzFiOJD4VGKKPsoZyhjVISqT2VTudR11EbqWeo4DUMzpgXQUmmltG9og7QpFYqKnUq0Sp5KjcpxFSkdoRvRA+jp9DL6Yfp1+jtVLVVPVb7qJtU21Suqr9XmqHmo8dVK1NrVRtTeqTPUfdTT1Lepd6nf00BpmGmEa+Rq7NE4q/F0Dm2OyxzunJI5h+fc1oQ1zTQjNFdo7tMc0JzW0tby08rSqtI6o/VUm67toZ2qvUP7hPakDlXHTUegs0PnpM5jhgrDk5HOqGT0MaZ0NXX9dSW69bqDujN6xnpReoV67Xr39An6LP0k/R36vfpTBjoGIQYFBq0Gtw3xhizDFMNdhv2Gr42MjWKMNhh1GT0yVjMOMM43bjW+a0I2cTdZZtJgcs0UY8oyTTPdbXrZDDazN0sxqzEbMofNHcwF5rvNhy3QFk4WQosGixtMEtOTmcNsZY5a0i2DLQstuyyfWRlYxVtts+q3+mhtb51u3Wh9x4ZiE2hTaNNj86utmS3Xtsb22lzyXN+5q+d2z31uZ27Ht9tjd9Oeah9iv8G+1/6Dg6ODyKHNYdLRwDHRsdbxBovGCmNtZp13Qjt5Oa12Oub01tnBWex82PkXF6ZLmkuLy6N5xvP48xrnjbnquXJc612lbgy3RLe9blJ3XXeOe4P7Aw99D55Hk8eEp6lnqudBz2de1l4irw6v12xn9kr2KW/E28+7xHvQh+IT5VPtc99XzzfZt9V3ys/eb4XfKX+0f5D/Nv8bAVoB3IDmgKlAx8CVgX1BpKAFQdVBD4LNgkXBPSFwSGDI9pC78w3nC+d3hYLQgNDtoffCjMOWhX0fjgkPC68JfxhhE1EQ0b+AumDJgpYFryK9Issi70SZREmieqMVoxOim6Nfx3jHlMdIY61iV8ZeitOIE8R1x2Pjo+Ob4qcX+izcuXA8wT6hOOH6IuNFeYsuLNZYnL74+BLFJZwlRxLRiTGJLYnvOaGcBs700oCltUunuGzuLu4TngdvB2+S78ov508kuSaVJz1Kdk3enjyZ4p5SkfJUwBZUC56n+qfWpb5OC03bn/YpPSa9PQOXkZhxVEgRpgn7MrUz8zKHs8yzirOky5yX7Vw2JQoSNWVD2Yuyu8U02c/UgMREsl4ymuOWU5PzJjc690iecp4wb2C52fJNyyfyffO/XoFawV3RW6BbsLZgdKXnyvpV0Kqlq3pX668uWj2+xm/NgbWEtWlrfyi0LiwvfLkuZl1PkVbRmqKx9X7rW4sVikXFNza4bKjbiNoo2Di4ae6mqk0fS3glF0utSytK32/mbr74lc1XlV992pK0ZbDMoWzPVsxW4dbr29y3HShXLs8vH9sesr1zB2NHyY6XO5fsvFBhV1G3i7BLsktaGVzZXWVQtbXqfXVK9UiNV017rWbtptrXu3m7r+zx2NNWp1VXWvdur2DvzXq/+s4Go4aKfZh9OfseNkY39n/N+rq5SaOptOnDfuF+6YGIA33Njs3NLZotZa1wq6R18mDCwcvfeH/T3cZsq2+nt5ceAockhx5/m/jt9cNBh3uPsI60fWf4XW0HtaOkE+pc3jnVldIl7Y7rHj4aeLS3x6Wn43vL7/cf0z1Wc1zleNkJwomiE59O5p+cPpV16unp5NNjvUt675yJPXOtL7xv8GzQ2fPnfM+d6ffsP3ne9fyxC84Xjl5kXey65HCpc8B+oOMH+x86Bh0GO4cch7ovO13uGZ43fOKK+5XTV72vnrsWcO3SyPyR4etR12/eSLghvcm7+ehW+q3nt3Nuz9xZcxd9t+Se0r2K+5r3G340/bFd6iA9Puo9OvBgwYM7Y9yxJz9l//R+vOgh+WHFhM5E8yPbR8cmfScvP174ePxJ1pOZp8U/K/9c+8zk2Xe/ePwyMBU7Nf5c9PzTr5tfqL/Y/9LuZe902PT9VxmvZl6XvFF/c+At623/u5h3EzO577HvKz+Yfuj5GPTx7qeMT59+A/eE8/vsbQFrAAAAIGNIUk0AAHomAACAhAAA+gAAAIDoAAB1MAAA6mAAADqYAAAXcJy6UTwAAAAJcEhZcwAAXEcAAFxHAcbToaUAACAASURBVHic7d3dk2RnfdjxM9Ls7sySfyFXTtlxUtxuqbiIShXy/makOEYGgQWOeNuVhIgNwtgubEArQNJKOwpGicybzYJdwtgxDrYFKt24VHubK6cq/0Ri7czsLnQuWmqGeenzTPc55/e8fD41F7A6c/pMd0+f7zzPc7o3ZrNZBwAQ6o7oAwAAUCQAQAYUCQAQT5EAAPEUCQAQT5EAAPEUCQAQT5EAAPEUCQAQT5EAAPEUCQAQT5EAAPEUCQAQT5EAAPEUCQAQT5EAAPEUCQAQT5EAAPEUCQAQT5EAAPEUCQAQT5EAAPEUCQAQT5EAAPEUCQAQT5EAAPE2ow8AGMtdFy5EH8LxXrt+PfoQgOwYIwEA4ikSACCeWRtoTvqkybHzPmt+O8CxFAk0ZOIFHPOb0yVACrM20Iqo9aTWsQIpFAk0ITYLRAnQS5FA/XIIghyOAciZIgEA4ikSqFw+gxP5HAmQIUUCAMRTJABAPEUCAMRTJABAPEUCdZq/U2pui0m9iytwEkUCAMRTJABAPEUCAMRTJABAPEUCAMTbjD4AYBS5XWWzkO2BAbGMkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQDHu+vChVP9O8A6NmazWfQxAAFGDYvXrl8fb+dAlYyRAADxFAkAEM+sDQAQzxgJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8RQJABBPkQAA8TajDwDI1N88cP/yDd72jWvTHAnQgo3ZbBZ9DEAueivkJOoEWJMiAbpujRY5SJcAK1Mk0LpBWuQgXQKswMpWaNrgOTLSPoHqGSOBRk3QDQZLgHTGSKBF0wxjGCwB0ikSaM6UoSBKgESKBNoyfSKIEiCFIoGGiAMgW4oEWhGYI0oI6KVIgCmIEmA5RQJNEARA5hQJMBFVBCyhSKB++aRAPkcC5EaRAADxFAlUzrAEUARFAkxKIQHHUiQAQDxFAgDEUyRQM1MkQCkUCTA1nQQcpUgAgHiKBACIp0gAgHiKBACIp0gAgHiKBACIp0gAgHiKBACIp0iAqb3tG9eiDwHIjiKBmjn3A6VQJABAPEUCAMRTJMCkTCQBx1IkUDkFABRBkQDTkUfASRQJ1E8HAPlTJMBEhBGwhCKBJqgBIHOKBFoRGyWSCFhOkQCje9s3rt114UL0UQBZUyTQkJCBikWOiBJgCUUCbZk4Sg6NjogS4CQbs9ks+hiAiSyC4Omf+5kJbu6kyZrXrl+f4NaBshgjgVYcjIPH/vb/jH1zS9aOGCkBjjJGAk04aazibx64f4ybS1nKaqQEOEiRQP2WxMHgUXJwnYooAdIpEqjf8jKYZ8H6XXLsmllRAiRSJFC5lBw5aIU0WX79jigBUigSqNlpc+SQJXVyqquIRQnQS5FAtcbogKP7TNyJKAGWc/Uv1GmaHEnXe3MuCYbGKRJo0YADEuklIUqAJRQJVGjN5SPjESXASRQJ1CbbHEk8AFECbVIkUJUiTueiBDhKkUA9CrqeRZQAhygSaMX6OTJs0IgS4CBFApUIXD6ycjqIEmBBkUANMl/NuoQoAeYUCRSv9HO2KAE6RQKlK2g16xKiBFAkULMicmROlEDjFAkUrNzlI8cSJdAyRQKlCsmRsStHlECzFAkUKbcT84DHI0qgTYoEylPHatYlRAk0SJFAbUrPkTlRAq1RJFCYylazLiFKoCmKBErSTo7MiRJohyKBYmRy9p24e0QJNEKRQBnyX806XhmIEmiBIoEahOfI2EQJVE+RQAFaWz5yLFECdVMkkDs5siBKoGKKBLLmFHuIKIFaKRLIV/6rWUOIEqiSIoFSBeZIeAmJEqiPIoFMlbh8ZMoOECVQGUUCOSoxR6YnSqAmigSy4zyaTpRANRQJ5MVq1tMSJVAHRQIlkSPHEiVQAUUCGSlo+UhWB9OJEiifIoFcFJQjSwSe+EUJFE2RQBacLAchSqBcigTiWc06IFEChVIkkDs5clqiBEqkSCBYHctHciNKoDiKBCLJkfGIEiiLIoEwpZ8R8w8mUQIFUSQQo+LVrFmd5kUJlEKRQI7KzZEMiRIogiKBAJaPTEyUQP4UCUxNjoQQJZA5RQKTctoLJEogZ4oEplPfatb6DliUQBRFArko7uy+RM7ndVECeVIkMBHLR/IhSiBDigSmIEdyI0ogN4oERufclidRAllRJDCu+laz1kSUQD4UCUSSI+FECWRCkcCIWlg+UsFPIUogB4oExtJCjixR1llclEA4RQKjcAIrjiiBWIoEhmc1a6FECQRSJDA1OZIzUQJRFAkMrPHlIxUQJRBCkcCQ2syR+n4uUQLTUyQwGGepQ4q+Q0QJTEyRwDCsZq2PKIEpKRKYghwplCiBySgSGECby0caIUpgGooE1iVHqidKYAKKBNbiVNQIUQJjUySwOqtZF1r4SUUJjEqRwFhaOEn3quwkLUpgPIoEVmT5SJtECYxEkcAq5EjLRAmMQZHAqTnfIEpgcIoETsdqVuZECQxLkcCQWs6RBn92UQIDUiRwCpaPrKDus7IogaEoEkglRziWKIFBKBJI4qTCEqIE1qdIoJ/VrPQSJbAmRQLrkiPMiRJYhyKBHpaPkE6UwMoUCSwjR07FHdKJEliVIoETOXMMpbV7UpTAChQJHM9qVtYhSuC0FAmsQo7QS5TAqWxGHwDkyPIRBvHa9evLn0t3Xbjg6VScoVLypId+7P1nyxgJHCZHGJCREkikSOCnOD2sSbEdJUogxcZsNos+BsiF1ayDOOlubPze8+xauP+JV5dvcO3xu6c5kgEteXyHemQnuIlYigR+wnzNIBTJSVqOkt4KOVZZaXLs4zvsYzrBTQSyshXeIEcYW5sLXVdrkUPfW1aasBrrSKDr5AhTaWpNyf1PvLpOjoy0K7KlSKCqcwD5ayRKxggIUVI3RULrWp7aJ0rdUTLqeIYoqZgigWXkyArcaSlqjZIJisEMTq0UCU2zfGRKhZ5ix1NflEwZCqKkPoqEdskRwtUUJdMngiipjCKhUQW90FO3OqJEHLA+RUKLrGYlK3VESQglVBNFAofJEaZXdJTEZoEoqYYioTmWj0zA3biCQqNEEDAURUJb5EisPM+p+Sg0SsKpojooEhri1Zz8lRUlUoABKRJaYTUrpSgrSjKhjSqgSKDr5AiZKSJKRADDUiQ0wfIRilNElGRFIZVOkVA/OUKhRAlNUSRUzkt2FKk3iGyjxIAEg1Mk1Mxq1gxpxNPKNkpgWIqEdskRSiFKEhm5KZoioVqWj1ATUUL1NqMPAEYhR6jPa9evL39i33Xhgud2bpRiOmMkVMhLALUyUkLFFAm1sZo1H+7qMYgSamXWhrY4R+bA5MKaTN8UZNgHou7cNEZCVSwfoRFGSqiPIqEecoSmiBIqo0iohBdfGiRKqIkioQZWs9IsUXLQtcfvjj4EVqdIqJ8coW4hUeLcz+AUCcWzfCRn7v9pGCmhAoqEssmRQjlBDk6UUDpFQsG8wsJBjUeJiaTSKRJKZTUrHDVllCgAhqVIqJMcoVltjpTIowooEopk+QgsMVmU6AAGpEgojxwpi0ckRFMjJcKoDoqEwtT0Mto4D+XYpokSNcBQFAklsZoVTqWFKJFE1VAk1EOOwFF1T9/IkZooEoph+QisZoIoCSkDOVIZRUIZ5Aiso9YooSaKhAIUPaoMmagpSq49frcAqo8iIXdWs1bAY5SJOqJEi9RKkVA2p7qiGf2a3jRRMl40yJGKKRKyZvkIDK7QS4LN1FRPkZAvOQIjmSxKBmkILdKIzegDgOMZz4dRvXb9+vLfsrsuXBik+xcxcf8Tr678vbRgYzabRR8DHGY1a31Oekw9lLFCfteWp0mhFbLknhzqPpzgJmIpEnJkvqZKxz6sHs1w/gA4raFGcE+6Y8fef7bM2pAdOdKUoaYGWNmo0zcrzNScpNCxE9JZ2Upe5AhMr+7PvqEUxkjIiFc9iDLZQtcKjH0/NHs/GyMhFyazIZaREmIpEsogR2ACooRAioQsWD4CmRAlRFEkxJMjjfBQlkKUEEKREMxLG54DGRIlTE+REMlqVsiWKGFiioR8yRGIJUqYkiIhjOUjkD9RwmQUCTHkCJRClDANRUIAr1/N0pqFEiVMQJEwNatZOcr5LH+ihLEpEvIiRyBbooRRKRImZfkIFE2UMB5FwnTkCFRAlDASRcJEvEhBNUQJY1AkTMFqVqiMKGFwioR4cqQpHu5qeCgZliJhdJaPkMKf1GV51+VX33X51Z+59wsnbeBXm9NSJIxLjkB93n351Y2um3/9g+Oi5LXr1++6cEFlciqKhBF5PYL6PPDkqxsb3cGvQ39azHMk6vAo12b0AVAtq1mhMu/9/Ktd123Mfnzo399z+ZXF/5YjrEyREEOOQOYOhcXP3vvkko1/9t4n//d3Pj7yEVE5RcIoLB9hCX9GZyXlsfi5ey933ax3m7/9zicGOihapEgYnhxhNXdduODpMZ6VK/Af3nt58C3hKEXCwPztC7GG/R38+Xc80c16RkcOevCJHwx46zRFkTAkq1lhMmPX/8+/44nVvtGsHKtRJExHjsAKQs7u/+jez/UuHDnJ+y6/LEpYgSJhMJaPwDryOYX/43s/t+Ye3n/55fmv/P1PvDrEEdEERcIw5Agkyqc8jrUkRzbu/MkpY/aj28v38/7LL7/4ibcPdlg0QJEwgMxfYclQC6P6xf2Ab733s13Xdd3sjs2TTg1vzOP8+Pbtjb693bG5+YEvvuKdwUmnSFiX1awMqMQLgIsrj4WDd/UHv/jKki0P+vHtWwk5cmb+P7Y3bu/OnGhI4onCuIo7u8BJ6iiPY33oiz/sjYyu6350u2emZu7Ozc3FaMrGqstjaZAiYS2Wj1CfisvjWB9OzpGUze48MuNzfuPWjdmZFQ6M1igSVidHKFpr5XGsDz/1SrfRXxo/unUzZW93njl76F/mu37Lxq3XRQl9FAkrKvfVnNaU+1wdNes/8tQPuzeLYbnbt271brN5Zh4cJ87RiBJ6KRJWYTUr6xv8chvlke7SUz9MvAbm1q2bvdVy5sjQyIJ1JKRTJAxPjrCO3sttlMeaLj31w5TNbqXN1CzJkYP+bpa0GS1TJJya5SNMQ3mMITFHbicMjXRdt9mXI/Od/D85QgJFwunIESZQSosU94R/+OkfJCxj7W7fTMuRs2d7P/tmo5v939m5pIOjeYqEUyjlPEGeyn3+FFceRz3y9A+6xHWsN/d7N9s8mxoZcoR07RbJgH/rz3dVwWvWclazkkh55OaRp19O3PL2zf61IylDI13XPfvY2zuftMdptFUk6S+Uiy17X6HKffEdVq0v5SxR7pO/qafro0//YCNhcOTWzb2UvZ05u5Wy2ZXH/mnKZnBQK0Wy8ktnI+MfvSwfaZnyKNdHn0laOHJrfy9lQufMuaQceeajcoRV1F8kS15M068wPOlyxHJfqU9FjjSi3OezJ+GxHnvmBymdMWyOPC1HWFXNRXLSy2vii9dis/l+jg6WlPvyfSqN/JhNKfcxVR6JPvbMy13aOtab+/2TNWffaJH+tSNPffTtCbcJx6u2SI59zV3t5ezgO0uW+FHp67CatWjllkfnqbWGj115OSlGuu7mXkKObCUNjTz1qBZhXRUWyYAtcujbF4Mlg7/7daGcMzJR7rPRr9LgPnYl6bKaW3u7XcIgypmt7ZS9tZYjQz1pT3oJHXv/2aqtSMbIkYP7OTiDUz3LR3JT7hPPs2Uav3blr1M+rWY/YWik67pzW1spMzVfePSfpewNelVVJKPmyGJv5Z4Vuq77yivvX/zvB+95ccmWciRQuc+xFZ4Ypf9O5ePXr/x10sKRvRsp1XJ263zKjX6+L0euPX53yn6g67qN2aySD2acIEeW3NBIt7W+gwmy3MFAsXxkGuWeiYd9ApT1O5WhT1z5q8Qt9/d2e7c5lzZT03Xd5Uf/eeKWh/S+bVqJHbPydZ1Z3USsSopkshwJubnVpLfIIQ/e86IBkmEpj+UUyTpCcmSFFlntzVvLSpMJTg1FnH1WVsOszfQv95mPM6/cIotv/8CTb/3yx//Xsf+1mqf+GHJ+ViznYS3U41f+MmWmZr5wpHfLxIUjTzz6LxJu8yfWeSP5xfeWlSaspoYiOVabr7BrtshBH3jyrUcHS9q8V49SHuTgk1f+KuXt4ff3bqRUy7m0hSOfO83oyIAfajPflS6pW/FFEjWEleEwyYA5stjhwR+zwZNZbg9xugYfrNZ8Mm2yZn/v9ZTNzm29JWWzqBw5uE9RUrGyi6TcE8bgBs+RxW4zbK/BlfsDKo8G/caVv+zS3o81ZXTkzaGR/smazyZP1oz6eb+ipGJlF8mxJnuNzudUPVKOLHZezWkvk8drBdU8BMfK51cpf59KWzjSnS5HenzmNAtHRs2RgzehS+pTcJF4CZsbNUcWN7H8zUtyU+5zo+7yOK3WPrSh16eufD9ls5v7N7qu6/3I37PnzqcMjXzm0X+ZcqNzE+TIwdsSJZUpuEiYUp5Rojxox29e+X7STM3+jZS9nTuXNDryu7nmyOIWRUlNSi2STN7AIHy0eYIBkhwoDxr3W89+v3fMo+u6/b3X0y6rSVrH+juPZJ0j1KfUIqGbPEcmGCZRHnDIbz/3P7uEKZiu6/Z2+0dHtrZT17F++uF/1X+TGTBMUhNFwikMFSXKA1LMcyTF/m7/6Mi57aShkRVaJHaARJRUo8giKfd8NqAi5mvKfaSUx8TCJ0Az9Onn/iLlI/HmQyO9ObK1nbSO9bcf/tcJt/lTzNcwlCKL5CTOIlHKPZd4zmSu2cttPv3cX6Rstr+b9n6s20nrWFfIkUwYJqlDVUUSorWXy/nETXEJ0trDRNF+57nvpXRGysKRLnl05Lce/jcpezvEAAkDUiRFKmLKJoTyoGi/+9z3uuT3Y72jb7vEN0Druu43V8qRrBgmqUB5RVLcX+eMQXlQn3mOpNjf6x8dScyRdVrEAAnDKq9IaIryoBGfee7Pp5+p+dTD/zZlb6UwTFI6RUIWlAct+8zV76W85cj+jbQ3QDufdJXvpy4VP1NDZRQJk1IenKTZC4A/ezVpdOQ0OdI/OvIbl9YdHTFlw+AUCaNQHgyl4guAP3v1zxO33L/xeu82iUMj3RA5AmNQJKyl1lMFjO1zaUMjXdft3/i7/vdjPf/3Unb1ydpbxFKSoikSViFEYB1PXP0fiTM1XcLFwIkzNY9f+ncJtwlhFAnApC5f/bOky2pupF1Wcz7psppPXPr3KXtjcG2ujlqNImEVh37HDJlAoifTcmQ3YeFI13XbaWtHPi5HKEF5RdLsgvycLR4RacI66v7t/vzVP+vS3o91N+GymsQW6bru1+UIhSivSOi67sF7Xgx8I/kl5wxpwhgquNzmC1f/NHEda3KO9M/U/Nql/5B2m4xo2KduxcneKRJGcvDXpvRzCazpC1f/NGWzYWdq5AjFuSP6AIYUFY93XbhQd7ce9OA9L572W+b3Tzt3ERz0xavf3ehmvV97u69vbHS9X+fPn0/ZmxyhREUWib+5u5XKYBDrhIU0oTVPXf3uRtf1fu3uvt51s96v89vnU/b2Xy79QvTPDaswa8MpPHjPi1/uBugJczpU7+mdP+m6lA+r6W7cuNG71fnz8w/y7V878tjFd/TfZL28PVrRihwjYS5kmOS169cHX6hl4ITKzHMkxe6NG71jHm/mSI/HLr5jyhxx7mdwpY6RnHSVYAVr8rN1MIAO3slDxYTrdOiquAA4MUeGXcfa+NAIdSi1SPIxf/WMOolOdhnwkvGYxc8+eJp06oQ3lfLHxpWd76SMPN9I+yDf82lX+T568d6EnUHuCp61KeLlaQJRS1yPmk/omNOhWc/svDTrZr1fu7s3NhJsnz+fsjc5smAiqXQVjpGU8rfUgMYeKTlt9JjToTVXdl7qUt+Ptf/TarbTFo50XffoxfsStxzDtcfvvv+JVwMPgMqUXSThU86xUzYHjRcla47BmNOhes/uvDT0+7H2eyS0RTJkgKQCBc/azB17TmpzkH+M6ZsB92lOhypd3fnjO7of937t3Xh978brvZfVvOX8dsre8skRHcCAyh4jWWKCuZt8BkgWBhwpGW95ijkdlgsf+0x3deePUzbb3b2R8sYk29tJkzWXLv5iymZNEUZ1KH6MpHMS+mkP3vPi+jEx2WrZwQdOFqMmpZzSSJfbY3p1549S3mh1d7d/4UjXddvb2yl7yzBH1ABD2ZjN+i8tK8Kxr1bjxUqGAyRHnXa8JIfLdsY462T+MHHUSU+DTB7KnZ1vJ265t7vXu83W9lbi3i5e/KXELacXuMQ1nySa4DQ08ZluYvXM2hw70jvS3E1uf6udZFEYy9MkhxBZMKdD5p7f+XbqOtbd3d5ttre3U3b1kYxbJFY+OcL66hkjmZsyUZ3eJmPgpDXZjpE8nzY6srfX3yJd121tVZUj0w+T5JYjxkjWVMM6koPGvvRGjoRwnQ45+K8730756N393Rsbs1nv1/bWdsreSsmRLr8+oDj1zNoszM9bh8406y/78DYYOTCnQ5Qv7XwrZbJmb2835bKara3tlLeH/9DFdybcZkYme8809VOl2mZtDhpq4PfQfpy3cmNOp0r5jE7/3s63Erfc2+u/rGZrK/X9WD9YWo4sjB0l2eaIWZs11Vwkc0tOV72PohYpizSpSSZLSUJypNwWOWikLsk2RzpFsrYKZ20OOXYSZ+5UJ7BqHvKKmdNhWF/euZb0UTUJLdJ13fbW+ZSZmg9cvD9lb/kbfAYn5xZhEPUXydySLkn5Rorj83RY05ef/2bKR+ft7+3dkbBy5NzWVlKOfOSXU46tFPOGWL9LtEgjWimSufS/oZ1yqjF4mnQGThrwwvPfTBkd2d/bTdnsXNpVvg/VlSMLi55YIU20SFPqX0cCRw2+4kSaDC5qHckLz/9h4pb7e/3vx3puK/X9WB/6yLsSt6zA8jQptELWWbOYz03EUiQ0zWLYnE2/iO+/7/xB4pZ7+/u922ydO5eyq1+9+O7EGyUTQ71unPRkHnv/2Wpr1gYOMadTnPE+1jsxR/b397quf4XJuXNJoyNyBBYUCXSd63Sa9+LON5IWjuzvJy0cOXcuZR3r+y8+kLAzaIVZGziROZ1Yky0l+f2dr6dsljJT0yVP1rzv4ntSNoN2KBLoJ01CTFAkX9n5WuKW+zdv9m5z7uzZxL09ePG9iVtCO8zaQD9zOlX66s7XUqZguq7bS8iRrbQc+RUtAidQJHA63nutDl9NGx3Zv7nfpaxjPZs0UyNHYAmzNrAuczrjGekC4K/tfDVls3mO9ErMkfde/JWUzaBZigSG5L3XhjXGUpKvR+TIe+QI9DFrA0Myp5Ozr+98JXHLlBx5s0X6/6h7z8UHE28XWqZIYBTeey03X9/5/cQtb9682btw5OzZs2kt8r7EGwXM2sB0zOmc1lCzNt/YeTFls5sJ19R0b+RIvwcuvj9lM2DOGAlMx5xOiD/YeTHlKt+UoZEuOUfeLUfglIyRQCTX6fRa83KbP7z631I2u3n7dspmZzeT/op716X/nLIZcJAigVyY0znW0btl+By5dat3m7NnziTeqByB1SgSyI6Bk2PtvnDfSf9p+6GXjv33b159oXe3KaMjiUMjv3zpoZTNgGMpEsiXNFlSISc5VCdLouTW7f6hka7rzmwmjY7IEViTIoEytDans0KLHLJIk2vHRcmwOXK/HIG1KRIoTPUDJ+u3yEHzLrl29fcO/uOttHWsZ9Ima+6/9MEVDgw4RJFAqepLk2FbZGEeJd+6+qX5/711+0e933Jm887Enb/z0odWPjDgIEUCNahgTmekHFnYfuilbz33pds/6h8d2bwzaWjknQ9rERiSIoGqFDpwMnaOzL30+t292yQuHPklQyMwNEUCdSooTabJkbk/ef2fLPmvd6bmyIcHOhzgJxQJ1C/nOZ0pc2TuuzeOj5I70yZr/tOljwx6OMAbFAk0JLeBk+lzZO67Pz1Sspl2TU3Xdb8oR2A0igRalEOaROVIN1/l+syz8/+9eecdKd/yHx++NOYRAYoEmhcypxOYI3PbD730R1eu3HFHYo48PPbxAIoEeMOUAyfhRdJ13ff2+i+96eQITEWRAIeNnSY55EjXddsPvfTSc88t3+Y+OQJTUSTAMmPM6WRSJN08Sp69ctJ/ve+RR6c8GGicIgGSDJUmr/zq3x9kP4OYv8H8d547HCX3PqxFYGqKBDidddIkqwGSuXmUfPfZZxb/8guPfDTucKBdigRY3WnrJKsBkrl5kXRvRokcgSiKBBhAYprkXCRArKRr8QGWe+369fnX8m0mOx6gOKnvnQyQ4mB2jHEV8eB2X7jPMAnkQJEAY1nUyTxNclvTCmTFrA0wut4JHQBFAgDEUyQAQDxFAgDEUyQAQDxFAgDEUyQAQDxFAgDEUyTARLw1KrCEIgEA4ikSoGlGbiATigQAiKdIgOkYkABOokiAdikkyIciAQDiKRJgUvkMS+RzJECnSIDpSQHgKEUCtEgVQW4UCRBAEACHKBIgRmCU6CHIkCIB2iJHIE+KBAgzfRzIEciWIgEiTZkIcgRypkiAYNOEghyBzG3MZrPoYwDouq7bfeG+kfYsRyB/igTIyOBRokWgFIoEyM4gXaJFoCyKBMjUyl2iRaBENq6TjQAAA0BJREFUigQoQG+dqBAonSIBAOK5+hcAiKdIAIB4igQAiKdIAIB4igQAiKdIAIB4igQAiKdIAIB4igQAiKdIAIB4igQAiKdIAIB4igQAiKdIAIB4igQAiKdIAIB4igQAiKdIAIB4igQAiKdIAIB4igQAiKdIAIB4igQAiKdIAIB4igQAiKdIAIB4igQAiKdIAIB4igQAiLcZfQBAbe66cCF949euXx/vSI51qsPrIo4Q2mSMBIh02j4AaqVIgIYIIMjWxmw2iz4GoH7LU2CymZFMDgM4yhgJMIUcTvbzHDnpSHI4QmiZIgEm5cQPHEuRAPEs7wAUCTC1kGGS5VM2QDhFAmTBMAk0TpEArTBAAjlTJECAiePAAAzkT5EAudAN0DJFAsSY+F3RTNlA5hQJkBHDJNAsRQLUzwAJ5E+RAGEmCAWDLlAKRQLkRUNAmxQJEGnUYRJrWqEgigTIjmESaJAiAWpmgARKoUiAYCNFg4EWKIsiAXKkJ6A1igSIN/gwiTWtUBxFAgDEUyRApkzcQFMUCZCFAWdYTNlAiRQJkC/DJNAORQJUxQAJFEqRALmQEdAyRQJkzcQNNEKRABlZc5jElA2US5EAuTNMAi1QJEBVDJBAoRQJkJeVk8JQChRNkQAFUBtQPUUCZGeFYRJrWqF0igQog2ESqJsiASphgASKpkiAHJ0qLwyfQAUUCVAM5QEVUyRAphKHSaxphTooEqAkhkmgVooEKJ4BEqiAIgHy1ZsahkygGooEKIwKgSopEqBU1rRCTRQJkDXBAY1QJEB5TNxAfRQJkLtjh0lM2UBlFAkAEE+RAKUyQAI1USTAFBYrP1ZbAiI+oHqKBACItzGbzaKPAajZSYMiKwx7HNzVyqMmSwZpjMRAoM3oAwBqkzgvc3SzYYNghemhY79FpsA0zNoA5VEJUB9FAhRDiEDFrCMBAOIZIwEA4ikSACCeIgEA4ikSACCeIgEA4ikSACCeIgEA4ikSACCeIgEA4ikSACCeIgEA4ikSACCeIgEA4ikSACCeIgEA4ikSACCeIgEA4ikSACCeIgEA4ikSACCeIgEA4ikSACCeIgEA4ikSACCeIgEA4ikSACCeIgEA4v1/egBDH7LdhREAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 5,
     "metadata": {
      "image/png": {
       "height": 100,
       "width": 150
      }
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Image(filename = fig_dir + \"h2m-tet.png\", width=150, height=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e520ca69",
   "metadata": {},
   "outputs": [],
   "source": [
    "#3-description (markdown)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77ae2401",
   "metadata": {
    "tags": []
   },
   "source": [
    "## H2M process: Liakopoulos benchmark"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f2dd89e",
   "metadata": {},
   "source": [
    "**Problem description**\n",
    "\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0785784b",
   "metadata": {},
   "source": [
    "**Analytical solution**\n",
    "\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "271ac1cc",
   "metadata": {},
   "source": [
    "**Numerical solution**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f66a6aec",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inspecting vtu_mesh_file\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table><tr><th>Header</th><th>Data Arrays</th></tr><tr><td>\n",
       "<table>\n",
       "<tr><th>UnstructuredGrid</th><th>Information</th></tr>\n",
       "<tr><td>N Cells</td><td>100</td></tr>\n",
       "<tr><td>N Points</td><td>202</td></tr>\n",
       "<tr><td>X Bounds</td><td>0.000e+00, 1.000e-01</td></tr>\n",
       "<tr><td>Y Bounds</td><td>0.000e+00, 1.000e+00</td></tr>\n",
       "<tr><td>Z Bounds</td><td>0.000e+00, 0.000e+00</td></tr>\n",
       "<tr><td>N Arrays</td><td>2</td></tr>\n",
       "</table>\n",
       "\n",
       "</td><td>\n",
       "<table>\n",
       "<tr><th>Name</th><th>Field</th><th>Type</th><th>N Comp</th><th>Min</th><th>Max</th></tr>\n",
       "<tr><td>bulk_node_ids</td><td>Points</td><td>uint64</td><td>1</td><td>0.000e+00</td><td>2.010e+02</td></tr>\n",
       "<tr><td>bulk_element_ids</td><td>Cells</td><td>uint64</td><td>1</td><td>0.000e+00</td><td>9.900e+01</td></tr>\n",
       "</table>\n",
       "\n",
       "</td></tr> </table>"
      ],
      "text/plain": [
       "UnstructuredGrid (0x7f651d4f4640)\n",
       "  N Cells:\t100\n",
       "  N Points:\t202\n",
       "  X Bounds:\t0.000e+00, 1.000e-01\n",
       "  Y Bounds:\t0.000e+00, 1.000e+00\n",
       "  Z Bounds:\t0.000e+00, 0.000e+00\n",
       "  N Arrays:\t2"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mesh = pv.read(vtu_mesh_file)\n",
    "print(\"inspecting vtu_mesh_file\")\n",
    "mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "9e5ae5ee",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inspecting mesh and initial conditions\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f651393b220>,\n",
       " <matplotlib.lines.Line2D at 0x7f651393b670>]"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2dffAkx1nfP8/s7/R+p7Ott/OdDp2CbBlxAswhERkbDLGRDYkqlT9wnJfCFXzIxlReigqpSlVCKpWyU1AVm8IgzooAB4IqZSlERQlk/gglpUBGNgLpTrJU4mSke5UlS3dnWbLufvvkj5ne7emd2Z3dnd552edT9bubnZmd7Z6Zfvrpbz/dLaqKYRiG0X2SphNgGIZh1IMZdMMwjJ5gBt0wDKMnmEE3DMPoCWbQDcMwesJGUz982WWX6TXXXNPUzxuGYXSSr3zlKy+q6uVFxxoz6Ndccw1f/vKXm/p5wzCMTiIif1t2zCQXwzCMnmAG3TAMoyeYQTcMw+gJZtANwzB6ghl0wzCMnjDToIvIXSLygogcLDkuIvKrIvKMiDwmIu+sP5mGYRjGLKp46L8N3Drl+AeA67K//cBvLJ8swzAMY15mGnRVfRD4xpRTbgM+rykPA9tFZEddCTQMwzCqUcfAop3A897nI9m+4+GJIrKf1Itn9+7dC//g5x48zMXnNzYmai1Q0nnyBal8fiJCIul3RCCR/P+SHU8kvar77O8n+99dxx0Xf//oO+461dNpGE1z+MVvctv37oxy7TqsYlFJKlw1Q1UPAAcA9u3bt/DKGjdf+2Zu3LV90a8bkVBVhgpDVdT7X8nvd+f5/yvp8dHn4DpDBSi5vq3RsjY4x8FV4kkyw0EAkmSGg4Ag3nXGDkjRvnY7DnUY9CPA1d7nXcCxGq5bylePnzGPzGg1SkllNBxXXq4ySs9J96sqwyHeOe465M8J9uP9lqK9Lh/OQdCgws85Akw6EBMOhecgoLA5dPc324f3bLJ7XQdPnzzDPR+7pZ6LBdRh0O8DPiEidwM3A6dUdUJuqZPrd2xl765LF/pu6B3mH1zgNRa+HPmHPHlOUAADT7UrBW1eycXI47zHLZJMlYtCb7DMOxzLUQVSFpLb33Yvct15/MipaNeeadBF5PeBHwEuE5EjwH8EtgCo6h3A/cAHgWeAbwEfiZVYx+ceepZLFtTQs3LhFTDJF54kX0jGBUeC73kF0ytoBNf2f2NdjeOomUyJnu7fq4q6uds/1tzHz0HC32rovqvC5kgPMl3ISHn0+ZcXdkhnMdMqquo/nnFcgZ+rLUUV+Oi795iG3jEWaRm5ZnGop5dKFSUtKGM9cI4DZJV8MkM3l8yBY9JB6GrLqJOhIqahG23A13JznbvDGZURqU4+rTKCbLtETzeKcc+EuTvovXuM/71qz3WeZ9J2DX3lzNLQ6462WEdv0DT02TivLpFkKW8wJ0cl5LzDrkZbGOU0qqG3kVka+khLLdVjZaKQJFlAsx8nnSShJh7q6eNCuu46+aK4juLCOPQq/ReQGtOgmTxhCLPQNcJm8hLPy1Xgm6psDtPcGMYsGtXQ24hp6P2iim6eNqPzLaU6Io+M9cA5Dr4GXjXyaJZuXhj7XhAYsQo6adC/eOgkB4+ebjoZvaYpycVGnBoxEGTUd3FWh5XGBkxzCMJY9lkx8b7zcPDYKfbu2hsln5006O+/4Urz0HuMjTg1YuA7CKnjkJSOOI3ZB3Lw2LZoeeykQbcoFyMWYeSKjfDsL7O863kiYtz+wmeM964oPHXiDF+wKJcx11x2MXsuvzhq7HKo01pQwXrgPLRBImwsqLkuOtrT3rH+op5desyiXPJ86fBLHD/1WrTRhUlWoNOCllhB6zF+QfMr9nMKqsOZ2umElz2H/KNK8dR2RjM4yS3SM3H25Ynjp7lpz5uj/EYnDfr1O7bx9iu3Fh7zvfNZemnorW96TaZ1L2eR3+1W4Yeeusp94FXmMT10Y/3Y+aYLo127kwb94cMv8cwL3yw9XjxHiJunpUJMs2fG/GiPrmzXQRvyUWU7/4wXa6k5T1kyFRvN/hcd/54KKu54eoc2R52o1ptqVOfIy69Fu3YnDfp7334Fe3fGCcw3ukOur4T6ZsQcajpYSIep7FImpRj9x/WfTcShF7TqXOjrLM6lo9Ci0EmDfub1s5x67WzTyTA6gvPonSFP93n7PQ/b6eTjs8KtcXhjVQ0d2tOiWVULz6cNeVrmHoRBFkr5XDxVePJ4vDE0nTTob2wOOf16atBLC1KFCZKM9UCYf6BRAiRJUqiNF3Wi5zTzwKMz+ssiEXVbBvHeiW4a9HNDzm4OR4VpI/E0Uq+A4ppHYSGbo3lkdI/MaV5oYJErhEPg3HCxKBfnKgh0crtu2pCnWPdACpyBWRF1r53dnOMX5qOTBn2oytlNzembrnk7a4RgnFfWaBvLTgWwqIfu9hv9ZFlHIf1+vPR10qBfdN4G2y/aMnfT1zzy/qJzVux1e+gj2uJimose7R7MO0eQG9eyke27eMHV1qrQSYP+1IkzvPKtNyqdGy4Tl6sAkrGhN8+q2/gx33VM1jWPh26OgjEPpqEHXHflJXznFZfMPG8cjVCt43TaNUS6s10HbchHzG1In7urxn0xTlURkfxxBXVx6f71su+Npk1NZlQaWQetsb68+eLzol27kwYd4PyNpOkkGB0mlFIWiUNfSNaZ4jwY7WD0vPCfezAh18gZdOfkxzaQnVPEV0+e4d3XXR4l7Z006Fdtu4C3bo83fNZoJ06uHjpjGxjNCX07K1BlLTRnXP14dPWMtDviznerWqWYm12VmPJ8DMatqmp9dGURdem1Jq8fcy2HThr0hw9/g+e+8a1RYTP1e2yU+n4vCju8Pa28Lg19Q2A0X7b0/76uM+OKfNJJGHniC0TUlbXQHj9ymlu/+6ooeemkQb/lO9/C3p2XWodUz/ElCr8ZXJeHPlQ3h/lkodxUEDFtZF3wI1L8cS2LDiabtkzdzu0XRctHJw367/zZ17jk/I1cUP9oNjvyHld6U8cPJe/FBQ8r8OIKm1aj7zd5Byapu1O0LfiFI+dl+8+WGZOtTSuIULiA9MC6aNaSUR+HwiZxKvQnbfrcPPvfcy037to+8rhsqbL1kVyaYB7dtIoEVNV5WNXCwsZqiVlGO2nQY05u03Vsjpp6GWmroTNQYa6gsgmc3OAl8GWfyZkg7VnWh+tQT5cKdA5d3pnLPc9h3unDdx6ZdBbLZD93fb9ufurEGe79+Lui5LOTBv0dO7bZItE9IVwdvXS49DB/zngOFW//cEYLzTOyZNeZ5S2Ng1rMW+46xa2gyVGdfotsnsVKpmnoPo/bEnR5Hjh0gsePxrspXaTLGvpkYSrpVJpjqHWVpQZFzAvuKzmPfMUtK/fdsqJ46NgpPrnrxij57qRB//EbrjIPvcM02ffh9zXM2ja6iwgMJK3V/Wk+pvZhJMVeeh1euY956AEHHjzM1guqJT3Vr2ZHuZTFKhN0cCVJ+QMsWurOnTuPIVlku25iprdJj9y/V1W2je7iWq2bqqSLBLWjon70+ZfZuyvOimudNOguyqUqoaeXazaFHmBBB0lOz/U+lza3Aq/TaA/zPo95Kivn6YX6a50xzJOhuRYRY4ypZNBF5FbgM8AAuFNVPxUcvxT4XWB3ds1fUdXfqjmtI5bV0JuKbQ49RKMbVPXqff01J/0M/cq/oMO3gh5bOF9ITmaqr1VVN7Fbp23NdxkHj51i7669Ua4906CLyAD4LPA+4AjwiIjcp6pPeKf9HPCEqv59EbkceEpEfk9Vq81xOyfLaOg68pxneN3eSMQYGq7RDRaSkyT1fBJJRvvm1XAnJCR37YTpjkQFDddolqY19JuAZ1T1MICI3A3cBvgGXYGtkr5FlwDfAM7VnNYRXz1+prOeblfTvc7Mq7u7bVfBD3UYp5PXaxGMvfj5O3+75KHnolcix5Pnn002b1QNt6fpOPSdwPPe5yPAzcE5vwbcBxwDtgI/parD8EIish/YD7B79+5F0gvA9Tu2RutUMBZj2X6KaasHhcYqbFk5/bruZngd+Lo3Vpkvzag1U6GjvKg/o2rkCsTrp2jaQy/KRfjG/zjwV8CPAn8H+BMReUhVc0M6VfUAcABg3759C5eaLx46GXUKSp82aHtNeFiL/GbT/RKLetLTto1m8Sv8tLUyXFm/RCx59ODR03yyKQ2d1CO/2vu8i9QT9/kI8ClVVeAZEXkWuB74i1pSGfD+G660OPSGaXqEJ3jN8BZ76MZyCGn/w4ZXya9iLp30e3H6Ipr20B8BrhORPcBR4EPAh4NzngN+DHhIRK4E3g4crjOhPp976FkuibjQqk/mcDJI/JdDcvvJXopBMvlyjPZnL8b4OnnPdZCMr5kkk9ET82zXQZXfESkIySu4BzEKoH9/Ym4bzTNyFEYVbbcr3Ebj0FX1nIh8AniAtPP+LlU9JCK3Z8fvAP4z8Nsi8jipTfpFVX0xSoqBj757z0o89EUiYuqam1sbfmebCOdaZ3Ly0oq8UItf7x+V3FxVvR+4P9h3h7d9DHh/vUkrp8tRLl3B7u/qcJJSWrEPCweqTTgIod47R+SML2817TiA78SMJbgq63UWzWo4LcplIpolu/6qO6ufPnmGez52S5Rrd3KkqEW5rIamdfJ1qVSc9wxu/pH1oygipa8jbJvW0FvHKjX0crTCnCP5TpbxgJF65oeJmzv18jUZkZIvIPPp6bm85vomyN1TY33I6+QtaDZExOZyCViVhj6LPq+YZBp6O1A0Xwkm5br5qJJNZjgaCJLkHY2c50veqTC6QycN+irj0OehT96sGfJ2IMiow3yVI06dZm3UT6NzubSRNsahN6k3x/DYzUNvD4IwkLG+HltvtoiYuJiGHmBRLvGx+9tOirzpiRV1QgdBC0ZKFkTO1BkRM4pccWlc8Tqei0WuuAzGffctyiXAolwWZ9GWRFHBmTe2fp0iV2Lhe89tp6iFYLH15qFP0FYNvStU0fpjFUCjv/ge8qrX8ewSpqEHtFFDbzN1ROOEUTZVZ0I01odkVMEnNv/7FMxDD2hHHPqipPre1IVmJ6JZKsxYKN5sh+Sbn0kibCTCQCabqGWedr4JG3jzRQWwoGPO5clYD1wF3rY1PNuGxaEHtCUOfVGqzAmedgrlPeFCD9n7PK0DzP1vrAduYFixxjxZQS8aEdOmEZhGRw26RbkYxnSE8RqnZ3WYj0OvMJd4WUTMUJXNYfbnbQ9VObep4+MVI1Im56TxnBLy++NjUS6NYFEueWLNCplrLTQ4atVYjBgRMc6TX9Wo1dhzkzeBaegBFuUySZGGvczc5KUFMGium07eb3w5cF4vv44Vg/r4blmUS4BFuUxGpORCu8JIlRl6fZXIldy5PS1oxiTOE98iyVw6fFlESpJ6ANODAjoYuTIP5qEHdDvKZR60ILIlfdkHIgyScQRLkn0e/Xmf/SiX0fnBOYMEBkmSncPo3CJv3oz5eqGaRq5knxpNSx+wKJeArke5zMOqPXH/PPPE1wPNZnSE+SNizBNvF5006BblshrsHq8PZzeHWYQK46iVYToHy+Yw2/YjXGZEuYy+G0bFbKb/D71rplElfoRJU9urwaJcAtYlyqXJ+daHqpjjNA7/6/N2IkIykFxHep0jPM07z2MaesCBBw+z9YJOJn1O8iNFw2bxINPOfU18Y+D0dKeJk2nkCYMsPGx8Trn+Pgo/sxWFcnnt47aTx22E52owDT1g/3uuXRsNfRrT9PWcJx7o5nWNOG3as4y9baTvRJluXucc7DbitB46adBNQ28HTXuWq9juA66DfFpFPW1u9KHC5nCY6uuqo22ng+c08VA3H+qSo0P7V7k+deIM9378XVGu3UmD3hcNPdYIz6pzk6tiOvkM2tBSWHYbUtkt9ZQTG+HZMKahB/RLQ4838+IsDX3WTIk273lKG1oKy247RwFMJ28a09ADTEOfTZNrnLr/2+KhmoYehyJ93a/YS52VBRyEouub5z9JJw26aejtpi1eaR3bXSEnpVE8YGyiop6ip+e17MmwV1fBO10dLQmBzc4Z6fgVdPPCOV7wHJFsf1dbGaahB3RJQ2+Dp2wsThtaClW3kwQ2cjKaecptxDT0gN988DBbOzSXS6h/T4Zoebq4UFDIAp27oJA5/XvU+eV1ShWmySvkSElagjSGa5BODVmD0rlguoQvv7R9WxVQ2Oyo57ouPPrcK6ah+/ysaehL0WR0jTUa1gNFPWfD4td9Ykp6nTToTx63udCN9cOvhDWoPGuZ8sGT7Ip0+Hlt5Ci9Q18rnyZD1rm6UXs9h8Y1dBG5FfgMMADuVNVPFZzzI8CngS3Ai6r6wzWmM8c7dmxrlYfedDx5qKUa/UQC6cs83m7SqIYuIgPgs8D7gCPAIyJyn6o+4Z2zHfh14FZVfU5EroiVYGhrHHq8ePJphcyVMydlLFvuVKfr5vm0lBiUkpWSpt+P6Zq/kWI6efdpWkO/CXhGVQ8DiMjdwG3AE945HwbuVdXnAFT1hboT6mNx6HFpcpZHY30oakWsQ1RO0xr6TuB57/MR4ObgnLcBW0TkT4GtwGdU9fPhhURkP7AfYPfu3YukFzAN3egveZ28OEZ8uq5cvXItqqRXZetcPqfFrxc7D92PX29aQy96xOEd2QC+H/gx4ELgz0XkYVV9Ovcl1QPAAYB9+/YtfFeb0NCbjiefVaubht4Pxl6rreHZV5qOQz8CXO193gUcKzjnRVV9FXhVRB4Evgd4mgg0paH7c4EXNvH8AoO/AEC+qZgk5QUqvb47p3268qRBKctD2T2aNEwTiyQU6PLufq0DqraGZ59pei6XR4DrRGQPcBT4EKlm7vN/gF8TkQ3gPFJJ5r/VmVAf09CboeloHm3AuLVlFOgqt+tmkbTkWhE2O2RlZhp0VT0nIp8AHiANW7xLVQ+JyO3Z8TtU9UkR+WPgMWBIGtp4MFaibS6X9aPJ5+3/dh3brlIba75jKW7qAt8z4sRzmvGwWItvolJchKJO97wmHujpgawZ6ulV7rE7L7Y9b1pDR1XvB+4P9t0RfP5l4JfrS1o5seZyCSf1L3voeAWoVE8v8lzDFy8orI6uelV9246FIGwkY+nJ4snXi6Y19NbxxUMnOXg0TqRLOL9JqmfHjb8umkN8Xu+vbur2Sru6XQd52QhUh5OedeBhVo1ymXBCCjrvjXZx8OhpPrlrb5Rrd9Kgv/+GK2vR0JvWhENvsE0aaNPecRu262SQCBsimCZsmIceUG+US7MjPPvMtAEeExEylBi3wIhViZxB8mGeTXvowMhRAFsxaN1pOsqldViUS/tpOm6/Li+7DS2FVW/HJEZ6i/oZ1mHEaRGdNOgW5WJMo+53o6onrwpDvMiVGfp3naM817U8OAlzMoplysjRYXk/RN75mC6zLto98fTJM9zzsVvqugU5OmnQ54lyqcNTXFYn7yKr8tj6gtPGUw8uWXtP0SjHNPSAzz30LJfMsWLRMqsFhQVHvONVRkt2ldgGfZ4oobqNYixGjsJI7jGd3JjENPSAj757j2noHcZWTDJWxaL6elejizpp0E1DN5rAl9LGGu3q4sntnZ+PUSVOdr+H7l5OGUQ49J9blX4P7xmT31+GaegBoYYef4Tn9IiKPhY009AncTr5htgIT2NxTEMPKNLQY8WSl2q6PdDJp9GUQVe0sfh1qfA88zq56TfG/JiGHmAaer+xFZOM2DjHYeT4BaN3J3RzzzEsnPI52y8J5S21bH+MkciOThp009CNWPiads74T2irQdjrqPO1+gyAE3os/ZTv2sqs/otyB6JcN6+yatJTJ87wBdPQx1xz2cXsufzivGc2Q0MvLYQVNXSTMNcD56ENEmGjNO47L9UURjoUeGZjKahYHrR3rL+oZ5ceMw09z5cOv8TxU6/NDjEqKHhVVsdJsgKdFrTEClqP8QuaX7GfU1AdRh/taQ55i3BKSKRn4uzLE8dPc9OeN0f5jU4a9Ot3bOPtV24tPOZ757P00tBb3/SaTOteziK/262iaErjgVeZx/TQjfVj55sujHbtThr0hw+/xDMvfLP0eJV1OksLZ9LfyJU+kn/Gi7XUnKcsmYqNZv+LN8WxCiruOEDqAKRYb6pRnSMvvxbt2p006O99+xXs3Rkn7MfoDv7AkTpnanTnDYcAw9JOUaNflIXqTkTEFIZCj0NfZ3EunT85Cp006GdeP8up1842nQyjI4w18ky3zoruuEJIC7Iz1E7/xtPWNTvO6HN1Dd1afN1g2tiLqiN6q/Dk8TirrUFHDfqLr77B1156de7vOa3cFUgraOX0baRo4QClKd5WXRp6X+6fUR9XbD0/2rU7adBff2OTV799bqHvyqgwg2mfVejmPVLyBjhVyNMDLrjEaeKLaehZpScycg40PW3cqa6jfwxjxFmTXPJceuEWrth2Qflw/4prNprv1G9ia+jTBhYZ64Eqk/bFb80VdNSfOPV6tPR00qBfeN6Ai84bTA7rHsIwix2uNCy86Yy0mR7GLc47HmGQjBermMdRKCrIRn/RMvsS2qesxbb1gi3R0tJJgz5IhIvOGzSdDKND+NEwq/DQ/aH/5jp0A9ew8ivgcYd4yejzgpbfLJ48cYZ/9P27ouShkwb90gu3sP2i85pOhtEkXmGqMjoT6o1yEYGBCIM+NWHWnCKDDhS2yMom8KrSEX7+lnjOaCcN+p88cZK//NuXC48liaepU22ZuELty8pp68nPHV4yjUPJsnbLRrkYxqLEfH06adDf87bLuP6qbRP7XScYVF+SrGxqVVdbuyaUiLR6u27akKdVbS8S5TKvhl62hJmFNq4fL33zjWjX7qRBf/nVs7zyreUGFvlem7HeaE5aSSv0TdV0TVMdzj3p1tQ5hDzN3mmzRvMM/WcXOINh/8r4eeanyR31pfj9J9n+zeH4OT998gw3X/uWKPnopEHfc/nFvO3KS5pOhhEZp1vHXuDCt6l+x9bIs04/pbHr2Xl93q6bNuRp2rYvzS4S/jxPiwzg8aM2fW6Ou/5fugRduKTc5HqM3lJz4k3ahXswge5a0MlRtPScayobcRk/g2q6eJXV3MsKqP/8B4mYDLKmjCt4JR3/U38L6umTZ/hB89DH7H/Ptdy4a3uuSbQKT8736Pyh8W3broM25KPJbT/UMOcoJDO8Ma+TddmlylxFZBhVqWTQReRW4DPAALhTVT9Vct4PAA8DP6WqX6gtlQFtWILO//22bddBG/LRhm1F2RzFoA8ndVPfUfAcADTvMMyz6EXuOkYtZI8kGy/gnkP+GeT08mH+WeEdd/v9Z1/2Lrjr+/XyUyfOcO/H3xUlnzMNuogMgM8C7wOOAI+IyH2q+kTBef8VeCBGQn2u37E12qrZxmoJW0+THVDxlhmsMuBHENJoFxhI0poWRMztmDSZp7IQ1BjLCYYh0j6PN7wE3U3AM6p6GEBE7gZuA54Izvt54B7gB2pNYQF/fPBE1HX5jNVRJjNMGzNQ11KDVdQM3+iv23ZMYqd9Muok86yHMxwCJqOQQofAXX/aouDT3q2DR0/zyV1757hb1ali0HcCz3ufjwA3+yeIyE7gHwI/yhSDLiL7gf0Au3fvnjetI2797qu4cdf2hb9vNEuTfR/G6lmV9+8jYUVfc8f5tMCJWf0eTXvoRakLi8angV9U1c1pmVHVA8ABgH379i1cvA48eJitF8zbn1vyEN3+idGCXoQMJd6g1+SC/AO1iZnKGXvagWedxB39mX5uMONrzqq8/9HvKaCwueLfncWjz70STTKuYhWPAFd7n3cBx4Jz9gF3Z8b8MuCDInJOVf+gllQGuCiXqowXtZjUUXOdHMMpXqLbP+dI1FV5hU14QcYkRfrrRIUU6rIz9NfC+OcC58HoBjHLaBWD/ghwnYjsAY4CHwI+7J+gqnvctoj8NvCHsYw5xF3Cqeus2gsyxuSchBiS0RR9t83P3Tk1ddc5fuRKrsPcv29Uj0iZjE7y9HHy+3WGTj6NRqNcVPWciHyCNHplANylqodE5Pbs+B1RUjaFd+zYZhp6i2myRZSuFtSMt+oK+MC85RwxW49dbBE1raGjqvcD9wf7Cg25qv708smazmIaejGqk7p5F0ecxvKCFsUvKEjJvWRKpEoyo1AyHnYd9o+02VtdZ2I8l7bq5NNoWkNvHfNq6LOwqAtjFRR5iW2IujBWS9MaeuswDd3oGq4Cn3AGVhAX3eUWi6+TTxvhObp/NY/wjMFTJ8/wv5vS0NuIaejNU97xFIzUHObPGXc8efuHM1pEnlGD8f+xRxjWSaY8gQiD7JeM2bhWx6pbNk4OjdG6aVxDbxtfPHSSg0dX46XHNBgxhzmvIt+FHVCRZ0J0BdARe9tYLaNwX1zFPsyNwKwWndLuuXOaHinaOt5/w5XmoTfAus1uaTTDIHFBCUkvZ7c0Dz2gziiXqqTRMJ7XWPaSBc02P6rF90Ld2qelYVOjc7zmX8MeepNzrqzSKzcPvVmco7yp8eYkb5JHn3/Zolx86o5yqUqVGf/Sjpi8jlyoL3ufizqyirzgdaap/DchnYW6cVfiq43m6aRBb8N86Ktm3fLbJlbdOvA7gXNS1HCKg4CLBCmOdJmcfXC839efY1BUeRVJc5NTJ1fPR5WoIXcvm6bRkaJtxOZDj0MbRni2oVO5Da0hkXTEaV8jYnItiQhrd4atl9iRK/NgGnpAExp6FSZHPwZ6+kzdfDxTIOS15/zoyvxsgqHGvFweKA/lCnVzSgpgUNCq6OmujDWllZuGvlqcowD91MmnYRp6QFMa+jTaEJfdddrgGYe0oaUQbhfp5jbi1ICOGvRVxqEvgqKtiQgxlqcNLQW37VfmYy05GDkaOgglHe5loyPLwlLr8BvaUCHGjA6rwsFjp9hrcehjuhCH3uRambWkv4EX3ZiN7z2njkMC5Ad1LRuzjUChp8+kU2LMj2noAesY5bJq7P62H0E8B2FYOdpj5HVP8fSLIklmRcQ4Tz43f0qBpz85AM1rRZCPcilMA3knJlZ0TiyePnmGez52S5Rrd9KgW5TLfNQxwrPJ1oJRjiAM/F7lJtOSJaHO1kKbRnjWhXnoAZ976FkuOb+TSfMAAx0AAA/ySURBVCftzS/pvHL7617fNBE2EmEQ/N60TjPT8415cRV4X0d41oVFuQR89N17Wq+hT6PJEafGeqCMR5xORqtMVtBFkTNlurmNOG0vnTTopqEbxnRcBMdwCGczfX2eEadlETFDVTaH2Z+3PVTl3KaOjxd8LyfreXKfr8tP09Pj434krm0xDT3ANPQ8sUZ45loLlOjw5vW3Ft9Dr+2azsNn/UZ41oVp6AHd1tBns8j6pokIg8T78z77GrrbToJzBsn4uvOO8LTW0vowrsRNJ18U09ADuq6hT6PJEae+t21x6P1nwhPOHAQbcdpdOmnQTUOPj93fflMWh+7i14eqnBum63huDrNtXz+foaGPvhtq7puT+ro/UjW2fj3rrqSYhr5S1lFDb8Jz96NvzMD3i2n6elHEi3nu9WEaekDfNfQyEkmX50ok08Q9fbwODX2QJNk55L8vAmJGfV0YOQojj9V08joxDT2gzxr6NJqcr9yiWdYHdZPLWTRL5+ikQTcN3TAWx5fSpsWnx9LQU82+yXjz8juTYhr6Sumrht70mqVWSa4HzuveIklhSOp4egcbQRoD09AD+qyhLzPHyzIa+qgQSr5ZbHO29BPVdM6V7FOjaVk3TEMPWFcNfRZ1zKpYutAB/vk2q+K64vR1sFkV20gnDbpp6IYxHT+2e5F50ktluuwaszT0vBS4nvOel9G4hi4itwKfAQbAnar6qeD4PwF+Mfv4TeBjqvrXdSbUp8saeldWMrKRot3GecOJJHN70UXesa1kVB+NaugiMgA+C7wPOAI8IiL3qeoT3mnPAj+sqi+LyAeAA8DNMRIMq9bQ/Z7v5bed9ixCNj95+vIPkvH+tHCJd07WHPX3pyM9gnPEK5jjQopk+nnZYBC/04p8B9moUCczOswYh66FA1CsUlg9rgK3ucnbR9Ma+k3AM6p6GEBE7gZuA0YGXVX/zDv/YWBXnYkMMQ19fppsGbgIGt/r78K2MT/lo0ADB6EoDt1aBktTxaDvBJ73Ph9huvf9L4A/KjogIvuB/QC7d++umMRJTEPvFu5Z+c+sK9ttpo7O7rKpkXHfK9HTy9KTVuBZGkZa+hRHYZhPI37HPpN5Kxsopx3qmX/qxBnu/fi7oly7ikEversL756IvJfUoP9Q0XFVPUAqx7Bv376Fn0DbNfSY0Saz5iYfqmKOSn20oaVQtr2MTm7RJs3RdBz6EeBq7/Mu4Fh4kojcCNwJfEBVX6onecU8cOgEjx+Nd1OWxW8KIvlmpx83XtocXUCvDgd0GPXRhpaCM+ROzkor9nQloiorEM2KXAmjTIomg1u05dJ0xVd1e1UcPHaKvbv2Rrl2FYP+CHCdiOwBjgIfAj7snyAiu4F7gX+mqk/XnsqAH7/hqlZp6E2P8MylxfTfWmmLsXH9EEkCG1klXzbQy0Z4tptGPXRVPScinwAeIA1bvEtVD4nI7dnxO4D/ALwF+PXswZ9T1X2xEn3gwcNsvaBdIfSzV/qRiaZr2uzNTz2aJDMK3uic9Hx3bl2ehqKFzevpETGzO8C6POK0LR46hDMhWuXdRZqOckFV7wfuD/bd4W3/DPAz9SatnP3vubZVHnrfaLIPoC3esBEfG3FaP+1ycytiUS79pWkPONxeNX7UxiIjPKeFkcL0yBUn68TO32i7MH9TRo1696Bq5MxEhAzNR8Q0HeXSOpqIcmk6jttohlW3DgRhI0kjV2bJWL73aXHc3aHpKJfW8ZsPHmbrimdbdC+83+Sb0MVn6t/jWRMhX1jz0S/5mRWrlDPndMQsk6q+UQjuR+I1jYMmM8H9KNTmmeygC+9vE/Ym51GuaNtV6kY/efS5V2y2RZ+fNQ29MboyF43RXUYd884J6tmKSTFlrU4a9CePn246CYaxMH5YaqHW7WnHRVr3tNGSuUrUrzzpltdf1Ok+6TzMGVfP9L4I12EfG9PQA96xY1sUD73peHLr6F0PkgQGiEV3rCmmoQfEjENvMp58GVahoc/6/UJjNKduPjG6tvQ55DsJu2Sn3LMaqsKw2bQYq8c09ACLQ28nTbZwjPXA6evQ3fh109ADTEM32k6oy8aac2WsuY8r0b7j9PW8Jt6d+HXT0AOW0dD7GqXhxzIbzePmXOljlIaxHKahBywz2+JY2+1nHPU6eGhtJHOUK0ekTHjWviPB9OkTNLiOSU7d4tCxU3xy141Rrt1Jgz7PbItFoVuFnvUwf05h+NiwuDDlCusw7/XD6kcb1sEq09vW7XnIV/L1OwRlHc02E2L3MA89YN4olzaP8FxXmpnRkdFzhnzlV7Y9DyNHYVQhmOtsTNL4bIttw6Jcuo2OZIMZUS3DGa2hrCMLZkgXw6APZA5D24aWQhdbeGWsOk99HnFaRCcNus22aCzKIu+N+44fJTHSylndTIjG/OScgNJAiPH9zn0mzohTi3IJCGdbbHqEp6NPHlYb8tH0doiQzoTopJ7SicrCOGZfSsJWDFp3TEMPKNLQuzrC05jEv/cEunm+z6JET597TdZx09z9fhl5ndy8ZmN+TEMPMA2937RhxSRHG1oKXWzh+bQhT0XbTY04jSmfddKgm4ZuxCRfaQzrrTSYrqdbi291qE5q4r6+7suseOf4fSmFUu+MEadPnTjDPR+7JUqeOmnQr37zRex+y0XpB/9GTtHQJ7Rypo/qDHV2kfwEWG3brps25KmJbeehbSQy9swm5LtA4imKdCjwzMYSUrG2bsZ8PXjMNPQ8f/DoUf706RdGn8v17sXix3OdXWTC+brhWoVtyruSe74xYtaV1LorAqLjJroKQ5z3nSXGMBbg6ZNn+P7veFOUa3fSoP/k9+zgu3ZsK/V6jP7it55ije515w2HAMNcvLvRX/wnPNlxXhKxtEBs+1XbLoiWh04a9CePn+bVb5/LhxIy1rSM9aFuD92Pjhl4izVPFl5oV/PFqJfyUGVX2Q91uFDn/GNHT/He66+IkupOGvTv2nEpe3fGCfsxuoPrK1mVh57v6Gou30YcyiJ9lGDEaSDRDhLJZtas1p/11m+ah57jzOtnOfXa2aaTYXSEkeEfteCyFh1jYw/jGHPnYcHYeGt2nNHn6rMmWkRWN5gWulka2eJ73xV/J+Z6Dp006C+++gZfe+nVub/nyzJpt5cVtDJixiU3QRitUuZtOdllICAkS0e59OX+GfVxxdbzo127kwb99Tc2efXb5xb6rowKM1ikQhW6eY+UvAEGIXPOUdLnr5JW65rtkKyaL4tycee7W6JZ2I1zDjQ9LfPKs0R09P4Z8Ti7GW8h2U4a9Gsuu5gb3rqt6WQYLSe2hj4xpsGbWMvoJ06yG70jCzzqJ4+f5tbv3lF30oCOGvQtA+G8jaTpZBgtwte+fX0T4mno2ccsPl4YDFxqTGbpJM44z3h8/pQALlpqHt5y8XmLpK4SnTToX/7ayzz/jW/NPtEfgML0xRSq9lCvCzFHocZgPKozH3pYtrRgnRp6V+6R0Q6+fubb0a7dSYP+tqu28p2XXzL1HOd9wZQIhGBOY2M9UK1XQ3eVB1RbIMEc+PXmzU176CJyK/AZYADcqaqfCo5LdvyDwLeAn1bVv6w5rSNUlfM2ZpUKKzVGOW7QSCwNfWIOfXeOeQ6tZ/RuBCGJ43UVqs0GWja6+Ksnz/Du6y6PkvaZBl1EBsBngfcBR4BHROQ+VX3CO+0DwHXZ383Ab2T/R2HHpReyc/tFsS5vtJTSjs2KI/T8mRDn1dAl8MSN/uH6XIokuInFR+aYbjfk8aPNTs51E/CMqh4GEJG7gdsA36DfBnxe0zvysIhsF5Edqnq89hQDjz73chqHbs5ONFoZh+73hyTh6j7xNHSLJzdC0tabskn6/zwsGnJdhSoGfSfwvPf5CJPed9E5O4GcQReR/cB+gN27d8+b1hEfedee7HomrBiG0S2uu3J6/98yVDHoRTYzrJKqnIOqHgAOAOzbt29h/9pCFg3DMCapYhmPAFd7n3cBxxY4xzAMw4hIFYP+CHCdiOwRkfOADwH3BefcB/xzSflB4FQs/dwwDMMoZqbkoqrnROQTwAOkYYt3qeohEbk9O34HcD9pyOIzpGGLH4mXZMMwDKOISnHoqno/qdH2993hbSvwc/UmzTAMw5gH6100DMPoCWbQDcMweoIZdMMwjJ5gBt0wDKMniJu/YOU/LPJ14G8X/PplwIs1JqcLWJ7XA8vzerBMnr9DVQtn92rMoC+DiHxZVfc1nY5VYnleDyzP60GsPJvkYhiG0RPMoBuGYfSErhr0A00noAEsz+uB5Xk9iJLnTmrohmEYxiRd9dANwzCMADPohmEYPaF1Bl1EbhWRp0TkGRH5dwXHRUR+NTv+mIi8s+p328qieRaRq0Xk/4rIkyJySET+5epTvxjLPOfs+EBEHhWRP1xdqhdnyfd6u4h8QUS+mj3rv7va1C/Gknn+19k7fVBEfl9ELlht6hejQp6vF5E/F5Fvi8gvzPPdSuhohfLm/0in5/0b4FrgPOCvge8Kzvkg8EekqyT9IPClqt9t49+Sed4BvDPb3go83fc8e8f/DfA/gT9sOj+x8wv8DvAz2fZ5wPam8xQzz6TLVz4LXJh9/l/ATzedp5ryfAXwA8B/AX5hnu9W+Wubhz5akFpV3wDcgtQ+owWpVfVhYLuI7Kj43TaycJ5V9biq/iWAqp4BniQtDG1nmeeMiOwCfgK4c5WJXoKF8ysi24D3AP8dQFXfUNVXVpn4BVnqGZNO7X2hiGwAF9GNFdBm5llVX1DVR4Cz8363Cm0z6GWLTVc5p8p328gyeR4hItcA3wd8qfYU1s+yef408G+BYawE1swy+b0W+DrwW5nEdKeIXBwzsTWxcJ5V9SjwK8BzpAvNn1LVL0ZMa10sY4NqsV9tM+jLLEhdaaHqFrL0ItwicglwD/CvVPV0jWmLxcJ5FpGfBF5Q1a/Un6xoLPOMN4B3Ar+hqt8HvAp0oX9omWf8JlLvdA/wVuBiEfmnNacvBsvYoFrsV9sM+jILUnd1oeqlFuEWkS2kxvz3VPXeiOmsk2Xy/C7gH4jI10ibpT8qIr8bL6m1sOx7fURVXcvrC6QGvu0sk+e/Bzyrql9X1bPAvcAtEdNaF8vYoHrsV9MdCUGHwQZwmLRmdh0DNwTn/AT5jpS/qPrdNv4tmWcBPg98uul8rCrPwTk/Qjc6RZfKL/AQ8PZs+5eAX246TzHzDNwMHCLVzoW0U/jnm85THXn2zv0l8p2itdivxm9CQUY/SBqt8TfAv8/23Q7cnm0L8Nns+OPAvmnf7cLfonkGfoi0WfYY8FfZ3webzk/s5+xdoxMGfdn8At8LfDl7zn8AvKnp/Kwgz/8J+CpwEPgfwPlN56emPF9F6o2fBl7JtreVfXfePxv6bxiG0RPapqEbhmEYC2IG3TAMoyeYQTcMw+gJZtANwzB6ghl0wzCMnmAG3TAMoyeYQTcMw+gJ/x8xf75vXdtEmQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(\"inspecting mesh and initial conditions\")\n",
    "#file\n",
    "reader = vtk.vtkXMLUnstructuredGridReader()\n",
    "reader.SetFileName(vtu_mesh_file)\n",
    "reader.Update()  # Needed because of GetScalarRange\n",
    "data = reader.GetOutput()\n",
    "#ic-values: pressure = data.GetPointData().GetArray(\"OGS5_pressure\")\n",
    "#points\n",
    "points = data.GetPoints()\n",
    "npts = points.GetNumberOfPoints()\n",
    "x = vtk_to_numpy(points.GetData())\n",
    "triang = tri.Triangulation(x[:,0], x[:,1])\n",
    "#plt.triplot(triang, 'go-', lw=1.0)\n",
    "plt.triplot(triang,lw=0.2)\n",
    "#plt.tricontour(triang, pressure, 16)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b327862b",
   "metadata": {},
   "source": [
    "**Running OGS**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "66ef056a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/ok/ogs-devs/build/release/bin/ogs liakopoulos_TH2M.prj > log.txt\n",
      "computation time:  33.25  s.\n"
     ]
    }
   ],
   "source": [
    "#run ogs\n",
    "import time\n",
    "t0 = time.time()\n",
    "#print(f\"ogs {prj_file} > out.txt\")\n",
    "print(f\"{ogs_bin_dir}ogs {prj_file_test} > log.txt\")\n",
    "! {ogs_bin_dir}ogs {prj_file_test} > log.txt\n",
    "tf = time.time()\n",
    "print(\"computation time: \", round(tf - t0, 2), \" s.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50bebe96",
   "metadata": {},
   "source": [
    "**Spatial Profiles**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "id": "c2f0bf2c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[   0.  120.  300. 4800. 7200.]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEaCAYAAAAlqOH8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZhUxb3/8fdHQXFhEwSX0QsKGGVViGASEhQBl6jBoOIS0OtPE+MSTTSYxMSYYCKJ0VwTl3iFC6hREBPxXhVjMIBxB0UW0QCCcRA3dlRk+/7+OGfGZpiBbqa7h5n5vJ6nn+muc6pO1TT0d6pOdZUiAjMzs0LbpaYrYGZm9YMDjpmZFYUDjpmZFYUDjpmZFYUDjpmZFYUDjpmZFYUDjlk1Seot6c0szz1f0j8LXadcSBotaXhN16MQcnlvrPAccAxJazMemyV9mvH6XEk/l3RfJflCUrv0+VBJMyStllQq6TeSGmScu1jSekktK5QxMy2nTRV12+rDUFKbNE+D9PU1kuZIWiNpkaRrKqnn+xXq00DSB5Ky/iJa2objK6ZHxDMRcVi25Vjx1NR7U/HfqCUccIyI2LvsAfwbOCUj7f4si9kTuBJoCfQE+gJXVzhnEXB22QtJnYE9qt0AEDAEaA6cAFwmaXCFc1YCJ2a8PglYkYdr207KH/Y7Hwccy4uIuDP9a3J9RCwB7ge+XOG0e0kCQ5mhwNg8XPs3EfFKRGyMiDeBiVlce0g+rg0gqY+k0ozX10pamPa4Xpc0cBt5fyvpn5KaStpF0nWS3k57X2MlNU3PK/uL+WJJ70paKukHGeXsLun36bF30+e7p8e2GsbL7J1WUqeLJC2QtFzSo5IOSNMl6da0bqskzZLUqYoymkoamdZziaThknZNj90paULGuSMkTU7L75P2kH8s6aO0V3luhXbeLOnfaa/1Lkl7ZL4PkoZJeg/4n0rem8Vpj3iWpI/TOraW9ET6fv1dUvOM83tJek7SSkmvSeqTcWyKpF9KejbN+7eMHvy09OdKJSMFx1T1b6A+ccCxQvkqMLdC2gtAE0mHpx8+ZwFbDNVJOkfSrB29qCQBvSu59iPAVyU1k9QsPWdihbzXSvq/Hb12hoVp+U2BG4D7JO1f4Vq7SPpvoAvQPyJWAeenj2OBQ4C9gT9WKPtYoD3QH7g2Y4jvJ0AvoBvQFTgauC7Xiks6Dvg1cCawP/A28GB6uD/J+9oBaEby/i2roqgxwEagHXBkmvf/pcd+AHRJA2Fv4EJgaHy+ztZ+JD3lA0n+KLlbUtmw2Ij0+t3Ssg8EfpZx3f2AfYD/AC6uom7fBPql5ZwCPAH8OL3mLsAV6e/iQOAxYHha5tXAw5L2zSjrHOACoBWwG5/36r+a/myWjhQ8X0Vd6peI8MOP8gewGDi+QtrPgfUkw1KZjwDaVVLGBUAp0LJiuSQfgr8mGfp6CmiQltOmivqMBtZVuO7qNE+DSs6/AXgN2D0jLUg+nO4Bvg18B/jvNC2q87tJ0/sApdvINxM4LX1+PvAiMA54GNgt47zJwHczXh8GbEh/R23Sdnwh4/hvgJHp84XASRnHBgCLM675zwp1Kn/v0t/x8PT5SOA3GeftndahDXAc8C+SwLbLNtrbGvgM2CMj7WzgHxmvjwaWkwS0syv8LjcCe2WkjQd+SjJ0+jFwaMaxY4BFGXnXA42qem/S9/DcjNcPA3dmvL4ceCR9Pgy4t0LbniQJjgBTgOsyjn0XmJQ+L3u/tvo3Wp8fHuO0bI2PiPMyE1TJDXdJ3wBuIvlg/qiScu4lGW5oS/ZDWjdHRPlf60omGCyq5NqXkQyV9Y6IzyopZyxJsBPJh0lBSBoCfJ/kQweSD+3MyRLtSHshEbE+I/0Akg/gMm+TBJvWGWnvVDjeeRt5D9iB6h8AvFL2IiLWSloGHBgRT0v6I3A7cLCkvwJXR8TqCmX8B9AQWJp0OIGk5/BORrkvSXqLpGcwvkL+FRHxcSVt2ZfkXuGMjHIF7Jpx7ocRsW47bXw/4/mnlbzeO6MdZ0g6JeN4Q+AfGa/fy3j+SUZeq4SH1CxvJJ1A0nM4JSJmV3ZORLxNEixOAv6Sx2v/J3At0DciSqs47RmSYaLWQEGmJkv6D5LfwWVAi4hoBswh+WAsM4+kF/hExlARwLskH3JlDib5az/zA/GgCsff3UbesmMfk3xQl9Vxv200YYtyJO0FtACWAETEbRHRHehIMiR1TSVlvEPSw2kZEc3SR5OI6JhR7qXA7un1flghf/P0uhXb8hFJQOiYUW7TSCa7lMnn8vfvkPRwmmU89oqIm7LI62X4K+GAY3mRjv3fD3wzIl7azukXAsdV+Cu2Otc+F/gV0C8i3qrqvEjGOk4BTk2f74iGkhplPCqOEuxF8mHzYVq3C4CtbqxHxAMk9w3+LunQNPkB4CpJbSXtnbZpXERszMj6U0l7SupIErTGZeS9TtK+6Y3rn/H5/bHXgI6SuklqRDJEWpU/Axek5+6e1uHFiFgs6YuSekpqSBLE1gGbKmnbUuBvwO8kNUnvVx0q6Wvp76QDyX2R84BvAT+U1K1CMTdI2i29x/N14KGI2EwSzG+V1Cot60BJA7bRnuq4DzhF0gBJu6bvdx9JJVnk/RDYTHIvzlIOOJYvPyW5Sf64Pv8OzxOVnRgRCyNiemXHlHzvp+IN/+0ZTvJX+MsZ176rimvPjYhKy09nRlVa5wyPk/yVXfb4eYXyXwd+BzxP0jPpDDxbRV3GAL8Ank6HCUfx+ZDjIpIP9MsrZJsKLCC533NzRPwtTR8OTAdmAbNJhsWGp9f5V3qdvwPz2UbvLiImk7yXDwNLgUOBsinmTUg+8FeQDHMtA26uoqghJDfRX0/PnwDsnwbo+4AREfFaRMwnCbz3pgEOkmGqFSS9mvuB70TEG+mxYWn7X5C0Om1TQb5nExHvAKel9fuQpMdzDVl8bkbEJ8CNwLPpDLdehahjbaMd/0PPzIol475Vwwo9njolnXZ8X0Rk04uwWsY9HDMzKwoHHDMzKwoPqZmZWVG4h2NmZkXhL35WoWXLltGmTZsdzv/xxx+z1157bf/EOqS+tbm+tRfc5vqiOm2eMWPGRxGxb2XHHHCq0KZNG6ZPr3TmblamTJlCnz598lehWqC+tbm+tRfc5vqiOm2W9HZVxzykZmZmReGAY2ZmReGAY2ZmReF7OGZWa2zYsIHS0lLWrdvegtD507RpU+bNm1e06+0Msmlzo0aNKCkpoWHDhlmX64BjZrVGaWkpjRs3pk2bNmRsUVBQa9asoXHjxkW51s5ie22OCJYtW0ZpaSlt27bNulwPqZlZrbFu3TpatGhRtGBjlZNEixYtcu5pOuCYWa3iYLNz2JH3wQHHzMyKwgHHzCxLK1eu5I477shLWZ999hlnnXUW7dq1o2fPnixevDgv5e7MHHDMzLKUz4AzcuRImjdvzoIFC7jqqqsYNmxYXsrdmTngmJll6dprr2XhwoV069aNa665plplTZw4kaFDhwIwaNAgJk+eTMXV+z/++GNOPvlkunbtSqdOnRg3LtlRfMaMGXzta1+je/fuDBgwgKVLlwKwYMECjj/+eLp27cpRRx3FwoULq1XHfPO0aDOrlW7437m8/u7qvJZ5xAFNuP6UjlUev+mmm5gzZw4zZ86s9Hjv3r1Zs2bNVuk333wzxx9//BZpS5Ys4aCDDgKgQYMGNG3alGXLltGyZcvycyZNmsQBBxzAY489BsCqVavYsGEDl19+ORMnTmTfffdl3Lhx/OQnP2HUqFGce+65XHvttQwcOJB169axefPmnH8HheSAY2aWJ88880zW51a2F1nFmV+dO3fm6quvZtiwYXz961+nd+/ezJkzhzlz5tCvXz8ANm3axP7778+aNWtYsmQJAwcOBJIvZu5sihZwJI0Cvg58EBGd0rR9gHFAG2AxcGZErEiPdQH+BDQBNgNfjIh1kroDo4E9gMeB70VESNodGAt0B5YBZ0XE4rSsocB1aVWGR8SYQrfXzAprWz2RmpJLD6ekpIR33nmHkpISNm7cyKpVq9hnn322OKdDhw7MmDGDxx9/nB/96Ef079+fgQMH0rFjR55//vktzl29Or+9vUIo5j2c0cAJFdKuBSZHRHtgcvoaSQ2A+4DvRERHoA+wIc1zJ3Ax0D59lJV5IbAiItoBtwIj0rL2Aa4HegJHA9dLap7/5plZXde4ceNKA0qZZ555hpkzZ271qBhsAE499VTGjEn+9p0wYQLHHXfcVj2cd999lz333JPzzjuPq6++mldeeYXDDjuMDz/8sDzgbNiwgblz59KkSRNKSkp45JFHgGQW3CeffJKvpudF0QJOREwDlldIPg0o622MAb6RPu8PzIqI19K8yyJik6T9gSYR8Xwk/dGxGXkyy5oA9FXy7g0AnoqI5Wnv6Sm2DnxmZtvVokULvvzlL9OpU6dqTxq48MILWbZsGe3ateOWW27hpptu2uqc2bNnc/TRR9OtWzduvPFGrrvuOnbbbTcmTJjAsGHD6Nq1K926deO5554D4N577+W2226jS5cufOlLX+K9996rVh3zTZWNIxbsYlIb4P8yhtRWRkSzjOMrIqK5pCtJhsZaAfsCD0bEbyT1AG6KiOPT83sDwyLi65LmACdERGl6bCFJr+Z8oFFEDE/Tfwp8GhE3V1K/i0l6T7Ru3br7gw8+uMNtXbt2LXvvvfcO56+N6lub61t7oebb3LRpU9q1a1fUa27atIldd921qNesadm2ecGCBaxatWqLtGOPPXZGRPSo7PydddJAA+ArwBeBT4DJkmYAlQ1SlkXMytZZiG2kb50YcTdwN0CPHj2iOrv8eZfAuq++tRdqvs3z5s0r+kKaXryzao0aNeLII4/Mutya/h7O++kwGenPD9L0UmBqRHwUEZ+QTA44Kk0vychfArybkeegtKwGQFOSIbzy9ErymJlZkdR0wHkUGJo+HwpMTJ8/CXSRtGcaPL4GvB4RS4E1knql92eGZOTJLGsQ8HR6n+dJoL+k5ulkgf5pmpmZFVExp0U/QDLbrKWkUpKZYzcB4yVdCPwbOAMgIlZIugV4mWT46/GIeCwt6hI+nxb9RPoAGAncK2kBSc9mcFrWckm/TMsC+EVEVJy8YGZmBVa0gBMRZ1dxqG8V599HMjW6Yvp0oFMl6etIA1Ylx0YBo7KurJmZ5V1ND6mZmVk94YBjZpalfK4WPW3aNI466igaNGjAhAkTytNnzpzJMcccQ8eOHenSpUv5gp0AixYtomfPnrRv356zzjqL9evX56UuxeKAY2aWpXwGnIMPPpjRo0dzzjnnbJG+5557MnbsWObOncukSZO48sorWblyJQDDhg3jqquuYv78+TRv3pyRI0fmpS7F4oBjZpalfG5P0KZNG7p06cIuu2z5MdyhQwfat28PwAEHHECrVq348MMPiQiefvppBg0aBMDQoUPLl7HJNHfu3PLVCbp06cL8+fMBuO+++8rTv/3tb7Np0yYgWZH6qKOOomvXrvTtW+kt9bzZWb/4aWa2bU9cC+/Nzm+Z+3WGE7deYqZMPrcnyMZLL73E+vXrOfTQQ1m2bBnNmjWjQYPkY7ukpIQlS5Zsleeuu+7ie9/7Hueeey7r169n06ZNzJs3j3HjxvHss8/SsGFDvvvd73L//fdz4oknctFFFzFt2jTatm3L8uWFncDrgGNmlie5bE+wPUuXLuVb3/oWY8aMYZdddslqOwOAY445hhtvvJHS0lJOP/102rdvz+TJk5kxYwZf/OIXAfj0009p1aoVL7zwAl/96ldp27YtwFarVeebA46Z1U7b6InUlHz1cFavXs3JJ5/M8OHD6dWrFwAtW7Zk5cqVbNy4kQYNGlBaWsoBBxywVd5zzjmHnj178thjjzFgwADuueceIoKhQ4fy61//eotzH3300UqDVqE44JiZZSmb7Qmqa/369QwcOJAhQ4Zwxhmff7VQEsceeywTJkxg8ODBjBkzhtNOO22r/G+99RaHHHIIV1xxBW+99RazZs2if//+nHbaaVx11VW0atWK5cuXs2bNGo455hguvfRSFi1aVD6kVshejicNmJllKZ/bE7z88suUlJTw0EMP8e1vf5uOHZMN5caPH8+0adMYPXo03bp1o1u3buX3jEaMGMEtt9xCu3btWLZsGRdeeOFW5Y4bN45OnTrRrVs33njjDYYMGcIRRxzB8OHD6d+/P126dKFfv34sXbqUfffdl7vvvpvTTz+drl27ctZZZ1WrTdtT1O0JapMePXrE9OnTdzh/Ta+qWxPqW5vrW3uh5ts8b948Dj/88KJe06tFV62y90NSldsTuIdjZmZF4YBjZmZF4YBjZmZF4YBjZmZF4YBjZmZF4YBjZmZF4YBjZpalfK4Wfdddd9G5c2e6devGV77yFV5//fXyY2PGjKF9+/a0b9+eMWPGlKd7ewIzs3oinwHnnHPOYfbs2cycOZMf/vCHfP/73wdg+fLl3HDDDbz44ou89NJL3HDDDaxYsQLw9gRmZvVGPrcnaNKkSfnzjz/+uHxNsyeffJJ+/fqxzz770Lx5c/r168ekSZO8PYGZWU0Z8dII3lj+Rl7L/MI+X2DY0cOqPJ7v7Qluv/12brnlFtavX8/TTz8NwJIlSzjooIPKzynbhsDbE5iZWblcF++89NJLufTSS/nzn//M8OHDGTNmTJXbEHh7AjOzGrKtnkhN2dHtCQYPHswll1wCJD2XKVOmlB8rLS2lT58+dWJ7At/DMTPLUjbbE8ycOXOrR2XBpuzeCsBjjz1Wvq30gAED+Nvf/saKFStYsWIFf/vb3xgwYMAW2xMAWW1PcOqppzJr1iz69u3LhAkT+OCDD4BkYsLbb7/NMcccw9SpU1m0aFF5eiG5h2NmlqXM7QlOPPFEfvvb3+5wWX/84x/5+9//TsOGDWnevHn59Od99tmHn/70p+XDXz/72c/Kh7pGjBjB4MGDue666zjyyCOr3J7gvvvuo2HDhuy3337l+cu2J9i8eTMNGzbk9ttvp1evXuXbE2zevJlWrVrx1FNP7XCbtsfbE1TB2xPkrr61ub61F2q+zd6eoDi8PYGZmdVqDjhmZlYUDjhmZlYUDjhmZlYUDjhmZlYURQs4kkZJ+kDSnIy0fSQ9JWl++rN5hTwHS1or6eqMtO6SZktaIOk2pd9akrS7pHFp+ouS2mTkGZpeY76koYVvrZmZVVTMHs5o4IQKadcCkyOiPTA5fZ3pVuCJCml3AhcD7dNHWZkXAisiol2abwQkQQ24HugJHA1cXzGwmZllI5+rRZeZMGECksj8GsYPf/hDOnbsyOGHH84VV1xRvqxNVdsTRARXXHEF7dq1o0uXLrzyyit5rWO+FC3gRMQ0oOLXWE8DyjZ7GAN8o+yApG8AbwFzM9L2B5pExPORvANjM/JkljUB6Jv2fgYAT0XE8ohYATzF1oHPzGy78h1w1qxZw2233UbPnj3L05577jmeffZZZs2axZw5c3j55ZeZOnUqUPX2BE888QTz589n/vz53H333eXL5OxsanqlgdYRsRQgIpZKagUgaS9gGNAPuDrj/AOB0ozXpWla2bF30rI2SloFtMhMryTPFiRdTNJ7onXr1lusZ5SrtWvXVit/bVTf2lzf2gs13+amTZtuc2mZQti0aVP5NX/wgx+wcOFCunTpwrHHHsvw4cOrVfawYcO47LLL+MMf/sDHH3/MmjVr+OSTT/jkk09YtmwZEcFnn33G3nvvzerVq5k8eTJ/+tOfWLNmDYMGDeLXv/415513HhMmTOCMM85g7dq1dOzYkeXLlzN//nz222+/Ldpx6aWX8uqrryKJ8847j8suu4y33nqLH/zgByxbtow99tiDP/zhDxx66KEsXLiQK6+8ksWLFwNw6623bhEYAdatW5fTv4eaDjhVuQG4NSLWVlhYrrJV5mI7x7aVZ8vEiLuBuyFZaaA636iu6W9k14T61ub61l6o+TbPmzev/Bvw7/3qV3w2L7/bE+x++BfY78c/3iIt81v3v/vd73jzzTeZNWtWpflzWbzz1Vdf5f333+fMM8/kjjvuYK+99qJx48Ycf/zx9O3blw4dOhARXHbZZfTo0YOPPvqI5s2b07x5ckfgsMMO4/3336dx48Z88MEHdOjQobyeBx98MKtWrSpfnw1gxowZfPDBB+U7i65cuZLGjRvz/e9/n7vuuov27dvz4osvcs011zBx4kR+/OMf07dvX6688ko2bdrE2rVrt1p9oFGjRhx55JFZ/35rOuC8L2n/tHezP/BBmt4TGCTpN0AzYLOkdcDDQElG/hLg3fR5KXAQUCqpAdCUZAivFOhTIc+UwjTHzOqzbLcn2Lx5M1dddRWjR4/e6tiCBQuYN28epaXJYE6/fv2YNm1apUv6lP1Bns3WBYcccghvvfUWl19+OSeffDL9+/dn7dq1PPfcc5xxxhnl53322WcAPP3004wdOxaAXXfdlaZNm2bVtm2p6YDzKDAUuCn9OREgInqXnSDp58DaiPhj+nqNpF7Ai8AQ4A8VynoeGAQ8HREh6UngVxkTBfoDPypwu8yswCr2RHYG2fZw1qxZw5w5c8p7i++99x6nnnoqjz76KP/4xz/o1asXe++9NwAnnngiL7zwAr17965ye4KSkhLeeefzOweVbV3QvHlzXnvtNZ588kluv/12xo8fz+9//3uaNWu21YZyhRq2LOa06AdIgsFhkkolXUgSaPpJmk9yv+amLIq6BLgHWAAs5PNZbCOBFpIWAN8nnfEWEcuBXwIvp49fpGlmZjnJ1/YETZs25aOPPmLx4sUsXryYXr168eijj9KjRw8OPvhgpk6dysaNG9mwYQNTp07l8MMP3+b2BKeeeipjx44lInjhhRdo2rQp+++//xbX/Oijj9i8eTPf/OY3+eUvf8krr7xCkyZNaNu2LQ899BCQ9JRee+01APr27cudd94JJPd/Vq9eXe3fXzFnqZ0dEftHRMOIKImIkRGxLCL6RkT79OdWgSAifh4RN2e8nh4RnSLi0Ii4LJ2tRkSsi4gzIqJdRBwdEW9l5BmVpreLiP8pTovNrK7J3J7gmmuuKcg1Bg0axKGHHkrnzp3p2rUrXbt25ZRTTgGS7QluueUW2rVrx7Jly8q3JzjppJM45JBDaNeuHRdddFGlM+mWLFlCnz596NatG+eff375Zmz3338/I0eOpGvXrnTs2JGJEycC8F//9V/84x//oHPnznTv3p25c+duVWauanpIzcysVvnzn/+c9zIzZ3rtuuuu/OlPf6r0vEMOOYSXXnppq3RJ3H777du8RteuXSv9fk7btm2ZNGnSFmlr1qyhdevW5cEnX7y0jZmZFYUDjpmZFUVOQ2rpTK/2QKOytHQFATMzs23KOuBI+n/A90i+xzIT6EUy6+y4wlTNzMzqklyG1L4HfBF4OyKOBY4EPixIrczMrM7JJeCsi4h1kGwFEBFvAIcVplpmZlbX5BJwSiU1Ax4BnpI0kc+XlTEzq/PyuVr0VVddRbdu3ejWrRsdOnSgWbNmAMycOZNjjjmGjh070qVLF8aNG1eep95sTxARAyNiZUT8HPgpyTf7v7HtXGZmdUc+A86tt95avhLB5Zdfzumnnw7AnnvuydixY5k7dy6TJk3iyiuvZOXKlUDt355guwFHUiNJV0r6o6RvS2oQEVMj4tGIWF+MSpqZ7QyuvfZaFi5cSLdu3fK60sADDzzA2WefDUCHDh3KV3k+4IADaNWqFR9++CERwdNPP82gQYMAGDp0KI888ggAEydOZMiQIUiiV69erFy5kqVLl25xjU2bNnH++efTqVMnOnfuzK233grAwoULOeGEE+jevTu9e/fmjTeSFbjff/99Bg4cWL7awXPPPVftdmYzS20MsAF4BjgROIJkAoGZWY15Zvy/+OidtXkts+VBe9P7zA5VHr/pppuYM2fOVotdlslle4Iyb7/9NosWLeK447ae8PvSSy+xfv16Dj30UJYtW0azZs1o0CD52C4pKWHJkiVAsmzNQQcdVJ6v7FjmemozZ85kyZIlzJkzB6C813TxxRdvsT3Bd7/7XSZOnMgVV1zB1772Nf7617+Wb09QXdkEnCMiojOApJHA1usqmJlZ1tsTZHrwwQcZNGgQu+666xbpS5cu5Vvf+hZjxoxhl1122eYWBHVpe4INZU/SnTSrfVEzs+raVk+kpuxID+fBBx/cah201atXc/LJJzN8+HB69eoFQMuWLevF9gRdJa1O96FZA3TJeF399arNzGqJfG1PUObNN99kxYoVHHPMMeVp69evZ+DAgQwZMmSLnke92J4gInaNiCYR0Th9NMh43aTaNTAzqyXyvT3BAw88wODBg7cY/ho/fjzTpk1j9OjR5dOmy3ogtX17AlU29lfpiVL3iJhRIe2UiPjfatdiJ9SjR4+YPn36Duev6b3fa0J9a3N9ay/UfJvnzZtX6VbLhbRmzRoaN25c1GvWtGzbXNn7IWlGRPSo7Pxcvvj535I6ZxR6NnBdDvnNzKwey2W16EHABEnnAl8BhgD9C1IrMzOrc7IOOBHxlqTBJEvbvAP0j4hPC1YzM7NKRMRWU36t+LK9HZNpuwFH0mwgs+R9gF2BFyUREV1yvqqZ2Q5o1KgRy5Yto0WLFg46NSgiWLZsGY0aNdr+yRmy6eF8fceqZGaWXyUlJZSWlvLhh8XbGWXdunU5f7DWdtm0uVGjRpSUlORU7nYDTkS8nVOJZmYF0rBhQ9q2bVvUa06ZMoUjjzyyqNesaYVqcy6z1MzMzHZYNqtFHyMPlpqZWTVl08MZCsyQ9KCk8yXtV+hKmZlZ3ZPNPZzvAEj6Asn2BKMlNQX+AUwCno2ITQWtpZmZ1Xq57Pj5RkTcGhEnAMcB/wTOAF4sVOXMzKzuyGWlgXLpFz4fTx9mZmbb5VlqZmZWFA44ZmZWFFkHHEntJY2SdPv2z640/yhJH0iak5G2j6SnJM1PfzZP0/tJmiFpdvrzuIw83dP0BZJuK5uyLWl3SePS9BcltcnIMzS9xnxJQ3ek/mZmVj259HDuBR4CegNI6iRpbA75RwMnVEi7FpgcEe2ByelrgI+AUyKiM8m07Hsz8twJXB79oD4AABRmSURBVAy0Tx9lZV4IrIiIdsCtwIi0nvsA1wM9gaOB68sCm5mZFU8uAWeXiHgC2AQQEXOATtlmjohpwPIKyacBY9LnY4BvpOe+GhHvpulzgUZpD2Z/oElEPB/JUqVjy/JUKGsC0Dft/QwAnoqI5RGxAniKrQOfmZkVWC6z1N6V1JZ05ej0w3yPal6/dUQsBYiIpZJaVXLON4FXI+IzSQcCpRnHSoED0+cHkmybQERslLQKaJGZXkmeLUi6mKT3ROvWrZkyZcqOtou1a9dWK39tVN/aXN/aC25zfVGoNucScK4E7gH2k3QBSS9hzrazVI+kjiRDY2UbvVW2xE5s59i28myZGHE3cDckW0xXZyvdmt6KtybUtzbXt/aC21xfFKrNWQ2pSdoFOIckyFwBHAJMBb5Vzeu/nw6Tkf78IOOaJcBfgSERsTBNLgUy18MuAd7NOHZQmrcB0JRkCK88vZI8ZmZWJFkFnIjYDBwfERsjYkJE/DQi7oiIddW8/qMkkwJIf04EkNQMeAz4UUQ8m1GPpcAaSb3SIb0hZXkqlDUIeDq9z/Mk0F9S83SyQP80zczMiiiXSQOvSrp+R1eOlvQA8DxwmKRSSRcCNwH9JM0H+qWvAS4D2gE/lTQzfZTd37mEZGhvAbAQeCJNHwm0kLQA+D7pjLeIWA78Eng5ffwiTTMzsyLK5R7OQUBn4BJJLwKzgFkR8VA2mSPi7CoO9a3k3OHA8CrKmU4ls+PS3tYZVeQZBYzKpp5mZlYYWQeciDgTki9YAh1Jgk9Pku/mmJmZbVPWAUfSyyS9mtnpz8cjYsy2c5mZmSVyuYdzGklvZjfgO8BiSYsLUSkzM6t7chlSe5dkOvEkAEmHk8wGMzMz265cFu88OPN1RMwjuZdjZma2XbnMUhsn6SBgEcl9nHXA4QWplZmZ1TlKvhuZQwapHckMtX2AJyOidDtZaqUePXrE9OnTdyjvqAt+QWxunecamZkVR+zyHv/vf67fobySZkREj8qO5TKkNlVSk4hYALQGGpOxFI2Zmdm25DKk1iwiVkvqDlwE/B/w33y+nIyl/vN/fuYF/+qB+tZecJvri0Ktjp1LwNmQLoo5BBgREeMl7diYk5mZ1Tu5BJw/AK8Bjfh8Z869814jMzOrk3L5Hs4YSQ8DmyLi03TywPOFq5qZmdUlOU0aINlm+lNJ3wG+Dny7YDUzM7M6JZelbSpOGmhOMmnAzMxsuzxpwMzMiiKXgHMbnjRgZmY7KJdJA2Ml/QVPGjAzsx2QSw+HiFib8XwBcEHea2RmZnVSLrPU2ksaJen2QlbIzMzqplxmqd1LsgFbbwBJnSSNLUitzMyszskl4OwSEU8AmwAiYg7QqSC1MjOzOieXgPOupLZAAEgSsEdBamVmZnVOLpMGrgTuAfaTdAFwAjCnILUyM7M6J6uAI2kX4BySIPMNoCswFRhVuKqZmVldklXAiYjNko6PiF8BE9KHmZlZ1nK5h/OqpOvTezdmZmY5yeUezkFAZ+ASSS8Cs4BZEfFQQWpmZmZ1Si5L25wJIGl3oCNJ8OlJ8t0cMzOzbcplpYETJb1A0rO5FvhXRFxdsJqZmVmdkss9nDuAHwC9gLuB30o6uyC1MjOzOieXgPN+RDwbESsi4u/AAOAn2WZO12H7QNKcjLR9JD0laX76s3nGsR9JWiDpTUkDMtK7S5qdHrutbBKDpN0ljUvTX5TUJiPP0PQa8yUNzaHNZmaWJ7kEnMWShkvaLX29AViTQ/7RJN/jyXQtMDki2gOT09dIOgIYTHKv6ATgDkm7pnnuBC4G2qePsjIvBFZERDvgVmBEWtY+wPUk95uOBq7PDGxmZlYcuQScAE4H3pH0T2ABMEVS+6wyR0wDlldIPg0Ykz4fQ/Kl0rL0ByPis4hYlF7raEn7A00i4vmICGBshTxlZU0A+qa9nwHAUxGxPCJWAE+xdeAzM7MCy2WW2tkAkhqRLNrZNX3cI+mQiDhoB67fOiKWpuUvldQqTT8QeCHjvNI0bUP6vGJ6WZ530rI2SloFtMhMryTPFiRdTNJ7onXr1kyZMmUHmpRYu3ZttfLXRvWtzfWtveA21xeFanNOG7ABRMQ6YHr6KJTKvlwa20jf0TxbJkbcTTIhgh49ekSfPn22W9GqTJkyherkr43qW5vrW3vBba4vCtXmXIbUCuH9dJiM9OcHaXopyRdNy5QA76bpJZWkb5FHUgOgKckQXlVlmZlZEdV0wHkUKJs1NhSYmJE+OJ151pZkcsBL6fDbGkm90vszQyrkKStrEPB0ep/nSaC/pObpZIH+aZqZmRVR1kNqkqYCp0TEaknfARoBd0TE+izzPwD0AVpKKiWZOXYTMF7ShcC/gTMAImKupPHA68BG4NKI2JQWdQnJjLc9gCfSB8BI4F5JC0h6NoPTspZL+iXwcnreLyKi4uQFMzMrsFzu4TRLg0134CLg/4D/5vNexTaVTTqoRN8qzr8RuLGS9OlUstNoem/pjCrKGoW3UjAzq1G5BJwN6b2RIcCIiBgvqZATB8zMrA7JJeDcBrxGMpR2bZq2d95rZGZmdVIu38MZK+kvwKaI+FRSO+D5wlXNzMzqklwmDbQAzgTWSZoLzI6ICwpWMzMzq1NymRb9V2Bf4FfAb4FVkuYVpFZmZlbn5HIPp3FE/ELS6RHxNUnfBNoVqmJmZla35NLDWZf+/EzSHhHxMHBSAepkZmZ1UC49nJvTpf7HA6MkPUcVi2CamZlVlHUPJyIeTpf4/x3wOMn6ZKcVrGZmZlan5DJLbT4wm+S7ODNJlrVZXKB6mZlZHZPLPZw/Ae8By4ATgTnpVs+/kNSwILUzM7M6I5d7OOdFRLeyF5LuAi4AVgO3AJfnuW5mZlaH5NLDWSWpS9mLiJgJ9IqIm4Ev571mZmZWp+TSw/k2cL+kmST3cA4DNqfHdst3xczMrG7JZZbaG8DRwCSgFbAA+LqkvYAHC1M9MzOrK3Lp4ZBugvZQ+sg0PG81MjOzOqmmt5g2M7N6wgHHzMyKwgHHzMyKwgHHzMyKwgHHzMyKwgHHzMyKwgHHzMyKwgHHzMyKwgHHzMyKwgHHzMyKwgHHzMyKwgHHzMyKwgHHzMyKwgHHzMyKYqcIOJK+J2mOpLmSrkzTukl6QdJMSdMlHZ1x/o8kLZD0pqQBGendJc1Oj90mSWn67pLGpekvSmpT7DaamdV3NR5wJHUCLiLZ3K0ryaZu7YHfADdERDfgZ+lrJB0BDAY6AicAd0jaNS3uTuBioH36OCFNvxBYERHtgFuBEUVompmZZajxgAMcDrwQEZ9ExEZgKjAQCKBJek5T4N30+WnAgxHxWUQsItl59GhJ+wNNIuL5iAhgLPCNjDxj0ucTgL5lvR8zMyuOnHb8LJA5wI2SWgCfAicB04ErgScl3UwSGL+Unn8g8EJG/tI0bUP6vGJ6WZ53ACJio6RVQAvgo8yKSLqYpIdE69atmTJlyg43au3atdXKXxvVtzbXt/aC21xfFKrNNR5wImKepBHAU8Ba4DVgI3AJcFVEPCzpTGAkcDxQWc8ktpHOdo5l1uVu4G6AHj16RJ8+fXJrTIYpU6ZQnfy1UX1rc31rL7jN9UWh2rwzDKkRESMj4qiI+CqwHJgPDAX+kp7yEMk9Hkh6LgdlZC8hGW4rTZ9XTN8ij6QGJEN0y/PfEjMzq8pOEXAktUp/HgycDjxAEiy+lp5yHEkQAngUGJzOPGtLMjngpYhYCqyR1Cu9PzMEmJiRZ2j6fBDwdHqfx8zMiqTGh9RSD6f3cDYAl0bECkkXAf+V9kjWkd5biYi5ksYDr5MMvV0aEZvSci4BRgN7AE+kD0iG4+6VtICkZzO4OM0yM7MyO0XAiYjelaT9E+hexfk3AjdWkj4d6FRJ+jrgjOrX1MzMdtROMaRmZmZ1nwOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVhQOOmZkVxU4RcCR9T9IcSXMlXZmRfrmkN9P032Sk/0jSgvTYgIz07pJmp8duk6Q0fXdJ49L0FyW1KWb7zMwMGtR0BSR1Ai4CjgbWA5MkPQaUAKcBXSLiM0mt0vOPAAYDHYEDgL9L6hARm4A7gYuBF4DHgROAJ4ALgRUR0U7SYGAEcFYRm2lmVu/tDD2cw4EXIuKTiNgITAUGApcAN0XEZwAR8UF6/mnAgxHxWUQsAhYAR0vaH2gSEc9HRABjgW9k5BmTPp8A9C3r/ZiZWXHUeA8HmAPcKKkF8ClwEjAd6AD0lnQjsA64OiJeBg4k6cGUKU3TNqTPK6aT/nwHICI2SloFtAA+yqyIpItJeki0bt2aKVOm7HCj1q5dW638tVF9a3N9ay+4zfVFodpc4wEnIuZJGgE8BawFXgM2ktStOdAL+CIwXtIhQGU9k9hGOts5llmXu4G7AXr06BF9+vTJqS2ZpkyZQnXy10b1rc31rb3gNtcXhWrzzjCkRkSMjIijIuKrwHJgPkkP5S+ReAnYDLRM0w/KyF4CvJuml1SSTmYeSQ2Apul1zMysSHaKgJMxIeBg4HTgAeAR4Lg0vQOwG8kQ2KPA4HTmWVugPfBSRCwF1kjqld6fGQJMTC/xKDA0fT4IeDq9z2NmZkVS40NqqYfTezgbgEsjYoWkUcAoSXNIZq8NTYPEXEnjgddJht4uTWeoQTLRYDSwB8nstCfS9JHAvZIWkPRsBhepXWZmlpL/0K+cpA+Bt6tRREsqTEqoB+pbm+tbe8Ftri+q0+b/iIh9KzvggFMgkqZHRI+arkcx1bc217f2gttcXxSqzTvFPRwzM6v7HHDMzKwoHHAK5+6arkANqG9trm/tBbe5vihIm30Px8zMisI9HDMzKwoHHDMzKwoHnGqQdEK6J88CSddWclzpvjwLJM2SdFRN1DOfsmjzuWlbZ0l6TlLXmqhnPm2vzRnnfVHSJkmDilm/QsimzZL6SJqZ7lc1tdh1zLcs/m03lfS/kl5L23xBTdQzXySNkvRB+uX6yo7n//MrIvzYgQewK7AQOIRk2Z3XgCMqnHMSyWoHIlmE9MWarncR2vwloHn6/MT60OaM854m2YdpUE3XuwjvczOS1T4OTl+3qul6F6HNPwZGpM/3JVm1ZLearns12vxV4ChgThXH8/755R7OjjsaWBARb0XEeuBBkn13Mp0GjI3EC0CzdN+e2mq7bY6I5yJiRfryBbZcULU2yuZ9BrgceBj4oJJjtU02bT6HZHHdf8MW+1XVVtm0OYDG6VqNe5MEnI3FrWb+RMQ0tr2Icd4/vxxwdlz5HjupzP13cjmnNsm1PRfy+Xp2tdV22yzpQJJNA+8qYr0KKZv3uQPQXNIUSTMkDSla7Qojmzb/kWTDyHeB2cD3ImJzcapXI/L++bWzLN5ZG2Wzx05W+/DUIlm3R9KxJAHnKwWtUeFl0+bfA8MiYlMd2Ug2mzY3ALoDfUkWy31e0gsR8a9CV65AsmnzAGAmySr2hwJPSXomIlYXunI1JO+fXw44O66qfXlyPac2yao9kroA9wAnRsSyItWtULJpcw/gwTTYtAROkrQxIh4pThXzLtt/2x9FxMfAx5KmAV2B2hpwsmnzBSTb3gewQNIi4AvAS8WpYtHl/fPLQ2o77mWgvaS2knYj2fLg0QrnPAoMSWd79AJWRbJvT2213Tanexr9BfhWLf5rN9N22xwRbSOiTUS0ASYA363FwQay+7c9kWQL+AaS9gR6AvOKXM98yqbN/ybp0SGpNXAY8FZRa1lcef/8cg9nB0XERkmXAU+SzHAZFRFzJX0nPX4XyYylk4AFwCckfyHVWlm2+WdAC+CO9C/+jVGLV9rNss11SjZtjmRr+EnALJLdeO+JiEqn19YGWb7PvwRGS5pNMtw0LCJq7bYFkh4A+gAtJZUC1wMNoXCfX17axszMisJDamZmVhQOOGZmVhQOOGZmVhQOOGZmVhQOOGZmVhQOOGZmVhQOOGZmVhQOOGYFIOl4SffmcP6XJN2wI2VKKpF01jby/UnSl7Oti1mhOOCYFUZX4NVsT063dbh+B8vsS7KvSVV6kmwVYVajHHDMCqMrsJ+kZyS9J+l4gHStromSpkt6SdJhafpDkr6SPj9c0rR0l8VrJC2oqsw0zy3AoHT3zbaZlZB0OPCviNiUkdZG0huS7pE0R9L9aVnPSpov6ejC/3qsPnLAMSuMriSrKfcGvgucK6khySra30/Xl/s5ULaVcSdgtqQGwP0ke610IdmBck5VZUbEP0kWnjwtIrpFxKIK9TgRmFRJ/doB/wV0IVnx+BySrSSuJtnZ0izvvHinWZ6lgWUf4OY0qQGwEvgG0BF4OF3YtAHwjKRGQMOIWCXpTOC1iCgbOnsd+GAbZUKyavGbVVRnAJUvurgoIman9Z0LTI6ISBembJN7q822zwHHLP+OIAkaZbtBdiHppXQFfhIRIzNPltSdJLCUnTsz43Ankh5KpWVKakGybPyGipVItw1oFhGV7WHyWcbzzRmvN+PPBSsQD6mZ5V9X4LWM111IlvFfCgyQtAuApM5Kujqd0+MAy0i2b0ZSN+C8tKyqymxL1ZtiHQv8Iw/tMcsLBxyz/OvK5wEEkl7KHGAUyf+5eZJmkuynEmwZcO4Fekh6GfhPYHFEvLWNMt8g2c9kjqQvVahHVfdvzGqE98Mx24lI2jsi1qbPrwGaRsR1O1jWK0DPyobbzGqCA47ZTkTST0m2N94APEsyo+2zbecyqx0ccMzMrCh8D8fMzIrCAcfMzIrCAcfMzIrCAcfMzIrCAcfMzIrCAcfMzIrCAcfMzIri/wNzcoRR8ZIhjAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import vtuIO\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "#Read simulation results\n",
    "#pvdfile = vtuIO.PVDIO(f\"{pvd_name}.pvd\", dim=1)\n",
    "pvdfile = vtuIO.PVDIO(\"result_liakopoulos.pvd\", dim=1)\n",
    "print(pvdfile.timesteps)\n",
    "yaxis =  [(0,i,0) for i in np.linspace(start=0.0, stop=1.0, num=100)]\n",
    "r_y = np.array(yaxis)[:,1]\n",
    "time = [0,120,300,4800,7200]\n",
    "variable = \"gas_pressure\"\n",
    "\n",
    "for t in time:\n",
    "    pressure_yaxis_t = pvdfile.read_set_data(t, variable, data_type=\"point\", pointsetarray=yaxis)\n",
    "    plt.plot(r_y, pressure_yaxis_t, label='t = '+str(t)+' sec')\n",
    "titlestring = \"TH2M:H2M: Liakopoulos experiment\"\n",
    "plt.title(titlestring)\n",
    "plt.xlabel(r'$height\\;/\\;\\mathrm{m}$')\n",
    "plt.ylabel(r'$gas\\;pressure\\;/\\;\\mathrm{kPa}$')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.savefig(f\"{prj_name}.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb6272f0",
   "metadata": {},
   "source": [
    "**Temporal Profiles**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d49ef766",
   "metadata": {},
   "source": [
    "**Contour plots**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "d5223cfe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Vertical cross-section\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANkAAAJDCAYAAACG3LF6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZhcZZn38e9dXb2nt3RC9g3IHvawyuaCLCqg4sLMCDgqA4qivI6ijDMi6Awjo4wiOLzooL4jyqphQEBAQRGQnSQkhBCS0NnX7qQ73dVVfb9/nOpOp9NdXd2pp6pJ/z7XVRfdp87yPCG/PKdOnec+5u6ISDixQjdAZH+nkIkEppCJBKaQiQSmkIkEppCJBNZvyMzsp2a20cwW9fG+mdkPzGy5mb1iZkfmvpn5NRz7LOFkM5LdBpyR4f0zgenp18XAzfverIK7jeHXZwmk35C5+xPA1gyrnAP83CNPA7VmNi5XDSyE4dhnCScXn8kmAG91+70hvWx/Nhz7LIMUz8E+rJdlvd6rZWYXE51eUVFhR007KBeHD2P6zDivv5ZM9fH2oPpcRPyoynhdjloYRlNy02Z3H13oduxPcvG3vAGY1O33icDa3lZ091uAWwDmHVrid98/KgeHD6PhrSTvecem9r7eZhB9rik+wE+oPy+n7cy1BzfcvKrQbdjf5OJ0cQFwQfqK23FAo7uvy8F+h7Lh2GcZpH5HMjO7HTgVGGVmDcC/AMUA7v5j4AHgLGA50AJ8MlRj8+WKy7bx7FMJgNLh0mcJp9+Qufv5/bzvwOdy1qIh4Hs3Rp+bZk1e94K7z+/5/v7YZwlHd3yIBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkElosHsw/KtlQFdzUdWajDZ+n+QjdA9gMayUQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAitYIZ2mZBm/Xz+rUIfPkgrpyL7TSCYSmEImEphCJhKYQiYSmEImEphCJhKYQiYSmEImEphCJhKYQiYSWMFuq0ok4qxsGF2ow4vkjUYykcAUMpHACna6OBx5WTGJWRML3YzMNhS6AfsfjWQigSlkIoEpZCKBFewzmSWM0tUlhTq8SN5oJBMJTCETCUwhEwlMIRMJTCETCUwhEwlMIRMJTCETCaxgX0YXJaBqlRfq8CJ5o5FMJDCFTCQwhUwkMIVMJDCFTCQwhUwkMIVMJDCFTCQwhUwkMIVMJDCFTCQwhUwksMLdINzq1LzRVqjDi+SNRjKRwBQykcAUMpHAFDKRwBQykcCyCpmZnWFmr5nZcjO7spf3a8zsPjN72cwWm9knc9/U/NqyZRnAvOHUZwmj35CZWRHwI+BMYA5wvpnN6bHa54BX3f0w4FTgP8zsbfs0CfcOXnt9AcAyhkmfJZxsRrJjgOXuvsLdE8CvgHN6rONAlZkZMALYCiRz2tI8ampqoKK8HiAxXPos4WQTsgnAW91+b0gv6+5GYDawFlgIXO7uHTlpYQG0tTVSWlrTfdF+32cJJ5uQWS/LetZyOx14CRgPHA7caGbVe+3I7GIze87Mnku0Nw+4sQU2HPssOZBNyBqASd1+n0j0r3d3nwTu8chy4E1gVs8dufst7j7f3eeXFFcOts3BlZbW0NbW2H3Rft9nCSebkD0LTDezaekP9h8HFvRYZzXwbgAzGwPMBFbksqH5VFU1gZZdmwFKhkufJZx+Q+buSeAy4CFgCXCHuy82s0vM7JL0atcAJ5jZQuBR4KvuvjlUo0OLxYqYMf1sgBkMkz5LOFndhe/uDwAP9Fj2424/rwXem9umFdao+pkAi9x9fuey/b3PEobu+BAJTCETCUwhEwlMIRMJTCETCUwhEwlMIRMJTCETCaxgJeGsNUHJ0oZCHV4kbzSSiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQRWsMfZDkepshiNB5cVuhmZPV7oBux/NJKJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKbH2eaBmf0UeL8V7f7jTra1sOKxn5PYuY2SEXUc+K4LiJdWkGxt5o3HfkbLpreon340k0/40F77W/77n9DWtJW5H/5HADpSSVY+/ktaNjdQVFbJge/8BKVVIwHY8vqzrHvpEQDGHf4e6qcfHb7DgJmdAfwnUATc6u7/1uN9S79/FtACXOTuL2Ta1sxGAr8GpgIrgY+6+7b0e18DPgWkgC+4+0OBu5g1jWT5cRtwRvcF619+lOrx05n3ka9RPX46619+DAArijPhyDOYeMwHet3RtpWvEIuX7rFs82vPUFRawbyPfp0xc09mzbP/C0RBXvviw8w6+3JmnX05a198mGRbS+5714OZFQE/As4E5gDnm9mcHqudCUxPvy4Gbs5i2yuBR919OvBo+nfS738cmEv053xTej9DgkKWB+7+BLC1+7Ltqxd3jSr1049m++pFABQVlzJi7IF0H/U6pdrb2LDoccYd/p49ljeuXkT9wfMBqJt2KE1rX8fdaWpYSvX4GcRLK4iXVlA9fgZNDUtDdLGnY4Dl7r7C3RPAr4BzeqxzDvBzjzwN1JrZuH62PQf4WfrnnwHndlv+K3dvc/c3geXp/QwJClmBJHftoLiiGoDiimqSu3b2u83a5x9kzLxTicVL9lieaG6iZEQtABYroqiknFRbM4mWRkoqa7vWK6msJdHSmMNe9GkC8Fa33xvSy7JZJ9O2Y9x9HUD6vwcM4HgFo89kWTjp1DLftrWjz/cXL2xfDLR2W3SLu9+Syza0bFlDW9NmJh13Dm07tvZ413vZwnpdbFhWx9vHPvd2kJ6t6WudbLbtaTDb5I1CloVtWzu4+/5Rfb4/a/K6VnefP5B9xsuraG9poriimvaWJuLlIzKu37xxFS1bGlj462vxjg6SrTt57f6bmPm+z1JSWUNi53ZKKmvxjhSpxC6KSisoqaxhx7o3uvaRaN5O1biDsmrfPva5AZjU7feJwNos1ynJsO0GMxvn7uvSp5YbB3C8gtHpYoHUTp7LltefBaIrgLWT52Zcf/TsEzj0/H/hkI/9EzPffxml1aOZ+b7PAlAzeS5blj8HwLY3X6F6/HTMjOqJs2has4xkWwvJthaa1iyjeuKssB2LPAtMN7NpZlZCdFFiQY91FgAXWOQ4oDF9Cphp2wXAhemfLwR+2235x82s1MymEV1M+Wuozg2URrI8MLPbgVM9leSV27/F+CNPZ+yh72LFYz9n87K/UlJZy4HvvrBr/YW/vpZUohXvSLF91SKmn3Ex5XVj+9z/qBnH8ubjv2TRHd+hqLSCA9/5CQDipRWMO+I9LP3tDQCMO+I04qUVQfsK4O5JM7sMeIjoMvxP3X2xmV2Sfv/HwANEl++XE13C/2SmbdO7/jfgDjP7FLAa+Eh6m8VmdgfwKpAEPufuqeAdzZK5F+bUtab4AD9h1EcKcuxsPbj+pufdff68Q0u8n1On57M5XawcPclnn/OlnLYx157/yf/JaZ9Fp4siwSlkIoEpZCKBKWQigSlkIoEpZCKBKWQigSlkIoEpZCKBKWQigenexTwYSPkBgHUvP8qW156BWIxJx51LTfqm3q1vvMC6lx/FMIorqpl26t8QLxsxJMsPZJKpjECP9XJWhiB9s/GNwKlAB3CVu98dqIt70EiWH7eRZfmBXdvWs23Fi8z58FeYfvpnWP2Xe/CODrwjxVtP/5aZZ13KnA99mfKR49j46pPA0Cs/kIVeywh0F6AMwVXARnefkd7f44H6theFLA8GUn5g++rF1B14BLGiOKVV9ZRV19O8aXXnnki1J3B3Uu1tXTOrh2D5gf70VUagu1yXIfh74F8B3L3D3TfnsD8ZZRUyMzvDzF4zs+Vmtte/Oul1TjWzl8xssZnl7V+JUDa1rQaYZ2bLs1h97kD73Ff5gfbmPUsGFFfW0t7SiMWKmHzCh3n13ut55farad22gVEzjgWGZPmB/vRVRqC7nJUhMLPOP4RrzOwFM7vTzMbkpiv96/czWbdh+zSiRj9rZgvc/dVu69QCNwFnuPtqM+vtD+1tw72DV5ueAFgGzN+Wqmi7q+nIDFvcvwGoMbPn0gv2ofxA76UEvCPFpiV/Yc65V1BSVc9bT93L+pcfZdwRp/W5zb6UH9iWqqCfPo/q1l/o0WczewTobRLcVVk1ILdlCOJEs6WfdPcrzOwK4HrgE1m2ZZ9kc+Gja9gGMLPOYfvVbuv8DXCPu68GcPeNe+3lbWR7+0YqimrYlWpKuHtiwtza/jbZnKvyA8WVtSSat3et1968neKKalq2rAGgtDqa41U37TDWvxJ9jgtRfiALGfvs7u/p6z0z66uMQHeZSgoMtAzBFqKJofeml99JdHEkL7I5XcymEtAMoM7M/mhmz5vZBblqYCG0dTRTXpS55kYPswfa577KD9ROnsu2FS/SkUrStmMLrU2bqRw9meKKGnZt30B7+rSyae0yymqjM54hWH6gP32VEeguZ2UIPJqZfB/RlUWAd7PnIBFUNiNZNsN2HDiKqPHlwFNm9rS7L9tjR2YXExWypCw2oL/EeTbg2eKvA9/I0Od/BaqyKT9QXjeWummHs/juf8diMSYf/yEsFqOksobxR7yXZff/CIsVUTKijqknfxwYeuUHstBrGQEzG090qf6sAGUIvgr8wsxuADaRLneQD9mELNvKQ5vdvRloNrMngMOIPtN0SZ+z3wJR+YHBNjq00tgIdqX6r4PYTRLot889yw/MOOvSXnc27vD37FXAFKJiOqNnn7DX8li8mIO61QjpbtSMY7sukAwV7r6F6B/knsvXEtX96Pz9AaJaIFltn37v28C3e1m+Cjh58K0evGxOF7OpPPRb4CQzi5tZBXAssCS3Tc2fmuIDaEk1ApSk+9yfGG/zPks4/YbM3ZNA57C9BLijs/JQt+pDS4AHgVeISnHd6u6LwjU7rJjFmFN9EkSfNbMJzmze5n2WcLK6raq3YTtd1qv7798Fvpu7phXW6NIpAIvcff6EubX9ndoudvcb8tAseRsq2L2LXlZCYtbEQh0+O+sL3QDZH+i2KpHAFDKRwBQykcAUMpHAFDKRwBQykcBUfiAPclF+oCOZ4I1Hf07bjs2YxaiZPIeJR78fQOUH0tubWT1wF3A0cJu7X5Zev4LozvuDiMoS3Ofuvc6LDEEjWX7cxj6WHwAYe8ipzDvvSmafewXNG1bS+FZ0M4rKD3Rt3wp8A/hyL8e83t1nAUcA7zCzM3PYn4wUsjzIRfmBWLyEqvEHAxArilNRP5FEczTLWeUHou3dvdnd/8yez7LG3Vvc/Q/pnxPAC0Q3uueFThfzyWDHlGjmUHvbDlpn16T/NtTQfu9Odkwxmhc2Uj5pStd6jKmlqbIJn7J7xlFq1y62rVvMhPeezI6RRmt7E4mZdeyoMSBOrKKcxtEt7GhowifUde3L36xlR7yJ4inZzY4OaI/yAX3MpO9tHmPndIJstu9Vehb/B4hOQ/NCIctCU7KM36/PNNkx81T8XPJUinV3/oLaY0+iZGR9emFvdQaM3usPZBewpmQZD62fnWGNIVl+IPMOzeLA7cAPOmf654NClhsDLj9QVFlFckcT8apqkjuaKKpMlx+oriHZuLv8QLKxkXhVTdfvGxbcSUn9KOpOOKVrWby6lmTjdopravFUilRbK7HyCuLVtexaubsOULJpO+VTDx50J3sYiuUH+nML8Hq+b+bWZ7ICGTFrLk0vRuUHml58lhGz5gFQOWseTQtfpCOZpH3bFtq3bqJs4mQANj/yAB1tuxh95rl77+ulaF87Xn2FimkHY2ZUHjyT5uXLSO1qIbWrhebly6g8eGYee9mnUOUH+mRm1wI1wBf3od2DopEsD8zsduBUTyZZcf3V1L/zdEae9G7W/vrnNL7wDPGaOsZ/LCoRUnrAWKrmHc6qH14HsRgHvO/DWCxGe+N2tj7xCCWjDmD1j78HQO2xJ1Jz1HFUH3ks6+/5JW/e8G1i5RWM+0i0r6KKSupPPY3V//V9AOpPfS9FFZUF+TPoIUj5gfQ+VgLVRBNuzwXeCzQRnaYuBV6w6JT5Rne/NXhPUcjywt3PByibMMmnXHJF1/JJn+y9/ED9KadRf8ppeywrrqllxre+1+v6seJixn+s9/IDNUceS82Rw6r8wNQ+Dluwqz06XRQJTCETCUwhEwlMIRMJTCETCUwhEwlMIRMJTCETCUwhEwlMIRMJTLdV5YmZXW7xOCt/eB01Rx1H3Qmn0LZ+DRsW3EVHoo3i2pGMPe/vKCorA6Bt/Vo2LLiTjrZWMGPyP3yJWHExrWvfYv09t+PJdiqnz2b0WR/EzOhIJll/zy9pW/sWReWVjPvoBRTXjSxwr3uX7/ID6W3OB75ONF1mLfB3+XputEayPDCzecBnikcdwJTPfpnmZa+S2LKJ9b+5g1GnvY+pl32FEXMOYduTfwDSc8bu/h/GnH0eUz//VSb9/eewoiIANtx3F2PO/ihTL/86iS2baXk9munc9MIzFJWVM+2LV1F7wils+v3/Fqq72chr+YH0PLL/BN7p7ocSPRjlMvJEIcuP2cDTFothRUWUTz2Ina8upH3LRsqnRo+XrThoBjtffQWA5jdeo3TMOErHRg80LaqoxGIxkjua6Ghro3zy1OhJmofPZ+fShQDsXLKI6sOjcgZVcw6lZUVUgmCIymv5AaKbgw2otOgW/Gr2fsbe7pXNRmbx6vcZx510upgfi4Bve0eKjkSC5mVLKJswiZIDxtG8dDEjZs9j56KXaU9P1mzfvAnMaPjZf5Fq2UnVvCMYedK7SDY1Uly9ewJnvLqWZFMTAMkdjcRrov/vVlREUWkZHS3NXZNBh5i8lh9w93YzuxRYCDQTPRn1cxk2WZt+ZbpzvwiYnOm4nRSyLCQScVY2jM60Ssap+O6+xMyua9+y+dY1v7iF0rHjIRZj7LkfY+MD97Lljw8zYtbcrlNC7+hg16o3mfIPX8SKS2i47WbKxk8klv681qs+SxAMTiIRZ3XDqEyrvG3KD5hZMXApUaWqFcAPga8B1/axyRJ3P6Kffb6Y7fEVstzot/yAu/+kbMKkWyd96jI2//5+4jW1lIwew8QLLwEgsXkjO5dFzwqP19RSMfWgrlGocsZsWtc1UH3YfNqbGrv2mWzaTry6OtqmjxIEAb2dyg8cnm7TG+nj30EvnwO7Ob6f/WW7DqDPZHnTeUrTvn0bO5YspOqQI0ju3AFEI9eWxx+h9ujoedCVB8+kbcNaOhIJPJVi18o3KBk9lnhVNbGSUna9tTIq+fbSc1Smyxb0VYJgiMp3+YE1wBwz6zwdOY0MT1B1967PdGZWZ2bHmNnJna+e6/RHI1n+3J3YuJ41/3MrY973IYrKK9j21BNs/+uTAIyYfQjVRxwDQFF5BXUnnBKVDTCjcvpsRsyMLqyN+cB5rL/3dry9nYrps6icHlWU6qsEwRCV1/ID7v6qmV0NPGFm7cAq4KL+GmlmnwYuJxpFXwKOA54C3jWQzipkeeLuJ/UsP1B3/MnUHX9yr+tXHzaf6sP2PhsrmzCJqZd9Za/lmUoQDDWFKD+Qfvzyj3t7L4PLib5ze9rd32lms4CrB7gPnS7KbqmWZoCi9Je9Aq2dp4VmVuruS4EBl/vSSCZdVlx/NURf/D5f4KYMFQ3p78N+A/zezLaR4fu1vihk0qVk1Bja1q9Z6O7zS6dNHLLfZOeLu38w/eM3zewPRHUbHxzofnS6KF0mfeYLhW7CkGBmZWb2RTO70cz+wczi7v64uy9I330yIAqZdEl/JpPoVq35RHeInAn8x77sTKeL0qXhtpsBSgrdjiFgjrsfAmBmPwH+ui8700gmXUafeQ7ADDObXui2FFh75w/untzXnWkkky4jZsyB6Iva3xW4KYV2mJk1sfv+yXIz25H+2d29eiA7U8ikpx1ETz75U6EbUijuXpTL/Slk0uX1a78G0Z3qw3okM7MFmd5397MHsj+FLE8GWn5g6xOP0PjCM2AxDjjrg1ROj570GbL8wPR/+leW/fMVL4b+nixg+YHTiO5rLAESwD+6+2M99rkAONDd52Vo4vFEswB+CTzDPj4RRhc+8mCg5QfaNq6naeGLTLnsq0y84GI2/u/deEcHoPID/Wy/GfhA+srghcAveuzzQ8DOLNo3lmi+2TyikJ9GNLXncXd/fCAdBYUsXwZWfmDpIqoPOYJYPE5xXT3FI0fR2rA6ePmBVTfv09dBAxGq/MCL6ZuMARYDZWZWCmBmI4Ar6HuiZhd3T7n7g+5+IdGd98uBP5rZ5wfWzYhClh+LgJO7lx9INm3vKj8A7Fl+oGl3KQGAeE0NyR2Ngyo/MBCJTRsgmnf1yuC7mpU9ygcA2ZYfmDCA7T8MvOjubenfryH6UrklmwaaWWl65Pt/RKUKfgDck822PekzWRYsYZSuzvgdbU7LD/Qtw8iUg/IDU79wJW9+79rlwAcsYSv3pc+FLD9gZnOB64geZYuZHQ4c7O5fMrOpWWz/M6JTxd8BV7v7oizb3CuFLDdyWn6guLqGZHpUA0g2NhKvqiFeXRu0/EBx7UiAhLuvKpswqb/Vh2T5ATObCNwLXNBZboDoQsZR6QmdceAAM/uju5/aRxM/QVRwZwbwhW4zzI1BfE+m08U8GVD5gVnzaFr4Ih3JJO3bttC+dRNlEyer/EA/5QfS01LuB77m7k927sjdb3b38ekJnScCyzIEDHePuXtV+lXd7VU10ICBRrJ8yrr8QOkBY6madzirfngdxGIc8L4PY7Ho38OQ5Qd2rV6Zm572L1T5gcuAg4FvmNk30sve6+79FdoJSiHLk4GWH6g/5TTqTzltr+Uhyw+kR8LZZvar0vET92lfmYQqP+Du19LP1UN3X0n0eatPZvaCux+5r+t0Usiky5izP0Ljc08tAb4JfKzAzSmk2f1cYTWiCZxZUchkL+6+NIsLH/uzWVmsk8p2ZwqZSA/uviqX+9PVRZHACjaSpcqMxoNKC3X47Pyx0A3Ir00PZrz5XAZJI5l0aXnz9UI3YUgws+Mth18yKmQie7sQeN7MfmVmF5lZb7eHZU0XPqRL2/q1AIeY2YKQ35MNde5+CUC6LPeZwG1mVgP8gaju4pPunvXVRY1k0qV0zHiApcCNBW7KkODuS939++5+BtFDJv5MdHfJMwPZj0Yy6and3R8e5t+T7cXddxHdfbLXHSj9UcjypLfyA63r1rDxvjvxZDK6R/H9H6Z84pSubdq3b2PljddRf+rpjDzxnUDY8gO1x76DDb+9I6f97k3A8gPHAJ3TbQz4prvfm97mKOA2oJwoKJf7QGe1DpJOF/Ogr/IDmx++j/pTT2fKZ79M/bvOYPPDe5YM2PTgb7puAO4UsvxAzVHHDbqPAxSq/MAiYL67Hw6cAfyXmXUOJDcDFwPT068zQnSsNwpZfvRafgCMjrbogY0dra3Eq3bPoti5ZCHFdfWUjB7TtSx0+YE8ClV+oKVbMdIy0pM803POqt39qfTo9fM+jrkHMzvTzJ4xs9fM7A4zy/oRtt0pZPnRa/mB0Wedy6aH72PF9d9i00MLGHXa+wDoSLSx9U+PUX/q6XvsJHT5gTwKVn7AzI41s8VEdewvSYduQnr73vaVyU1EdUGOIzoN/a6ZnZ/FdnvQZ7LcGFT5gca/PsnoM86hau5h7Fj0Eht+82smXnQpWx57iLoTTiFW2vOOmLDlBwZoSJYfcPdngLlmNhv4mZn9brD7AjZ0m/z5iJk9RXRl8fYstu2ikGWhKAFVqzL+PxlU+YHNv7+f0WdFj8AaMfcwNvz21wC0Nqxix6svs+nh++ho3QVmWHGcqjmHBS0/kMs+F6r8QLfjLzGzZqK5Yw3p7XvbVyYrzexa4FvpU9Z2ogrLA6LTxTzprfxAvKqaXSujMhS7VrxO8cjRAEz69Oc58IpvcOAV36D2uJMZedJ7qDv2JJUf6L/8wLTOCx1mNoXo0bMr06eUO8zsuPTtUhf0ccyeHPgQ8JaZ/ZndpeEG9EAOjWT5s1f5gTHnfJSND/wG70gRixcz5pyP9LuTkOUH8ihU+YETgSvNrB3oAD7r7pvT713K7kv4vyOLUuTufn66XWVEI+Jh6detZnagu2f1ZaJClie9lR8on3IgUy69IsNWMOpde15pDll+IF8Clh/4BT2qBnd77zn6KTuQob2twHPp14DpdFEkMIVMJDCFTCQwhUwkMIVMJDCFTCQwhUwkMH1Plkde4rRNThS6GZJnGslEAlPIRALT6WKemNnlVhxn7de/x4hTjqH69BPZdNMvaV+3CYCOll3EKsoZf83leDLFlv++m8SqNXiqgxHvOJKa90flB9pWNrDl1jvxRJLyQ2dS97cfwMzw9iSb/+8dJFauITaigtGXnk989MDKD+RLgcoPlBAVCDqV6L7Gq9z97iAd7EEjWR50lh+IjxvNuGsuZ9fLS2lfv5nRn/0bxl9zOeOvuZyK+fOomD8XgJZnF+LtScZf+yXGffPz7PjDMyQ3bQVg689+Q/1FH2L8dV+mfcNmWhcuA2DnE88Sqyhnwr//I9XvPZFtdz5YqO5moxDlB64CNrr7jPT+Hg/Ss14oZPmxR/mB0pnTaHlhcdeb7k7LswupPPbwaIGBtyXwVApvb8ficX8VZ0cAAA/8SURBVKy8jOT2Jjp2tVF68BTMjBHvOLJrPy0vvsqIE6PHZVUcPY/WV5er/MCeEzP/HvjX9Hod3e7OD04hy4+o/EAqRUdbgl2vvEZqy+5nQrcte5Oi6hEUjx0FQMX8Q7DSEhq++B3WXPFvVJ95EkUjKkhtayI+cnf5gaK6GlLbovIDqW1NFI3cXX4gVl5Gx86W/PVwYPJafiD9mFuAa8zsBTO708zGkCcKWW6MMrPnur0u7v6muy8Brktu2MLG//gpJZPGQdHuP/rmp1+m8tjDun5ve/MtiMWY+P2vM+H6r9L04J9o37il9xIDuw+y97KwczYz9tnMHjGzRb28zulrhz0MuvyAu88Fjga+lp4LFieaDf1k+umYTwHXZ9mOfaYLH1koanVq3mjLtEpW5QdKp028dezXL2HbXQ8Sr4tGJE+laHl+MeO++fmudZufeonyQ2Zg8SKKqkdQOn0KiZVrKJ0xleTW3eUHUtsaKaqLyg8UjawhtXU78ZE1eCpFx65WYpX7UH6gtYOa5a2D7vMQKz/wPNAC3Jt++07gUxn6llMayfKks/xAcst2Wp5bTMVx0cjVung5xeNG73EaGK+vpXXJG7g7HW0JEm+8Fa1TW02svJS25atxd3Y++QIVR0TXAioOn8POP78AQMuziyibfZDKD+wuP+DAfURXFiGa8PlqLjuUiUay/Lk7sWYDG2+4jZEXnENRepRpfmbPU0WAqncfz5Zb72LdVd8HoPLEo6JTTGDkBeemL+G3U37oTMoOnQnAiJPns/mWO1jzle8Sqyxn1KUDrlyWT4UoP/BV4BdmdgOwCfhk8F6mKWR54u4nlU6b6N1PCwFGfeaje60bKytl9GV/2+t+SqdNZPy3v7TXcisp7nOboaZA5QdWAScPvtWDp9NFkcAUMpHAFDKRwBQykcAUMpHAFDKRwBQykcCyCpmZnZF+ENpyM9trWkK39Y42s5SZnZe7JhZGY8NSgHlmtjyL1Y/aH/osYfQbsn7m9fRc7zqib+jf1ryjg9V/uQdgGVGf+9PY/yoyXGUzkmWa19Pd54G76f1mz7eV5k2rKauuB0ik+9yfZP+ryHCVzW1Vvc3rObb7CmY2Afgg8C6iKQZva+0tjRRX1va/4gAMpPzArkWvs/3OB/FUEiuKU/uxMymfczCg8gPZbG9mk4luAP6mu19vZhVEd94fBKSA+9y9z489uZbNSJbNvJ4bgK+6eyrjjswu7px/lNw1ZJ9nnFPpPi8C/s2qKrMqP1BUVcHoL17I+Gu/RP1nPsKWW+7o2p/KD2S1/ffZ+/lj17v7LOAI4B1mdmauOtOfbEKWaV5Pp/nAr8xsJXAecJOZ7TWl3N1vcff57j4/Xl45yCaHV1xRQ3vz9v5X3K2ODH0Grgb+J149IqvyAyVTJhBPzxMrnjAGb0/i7UmVH+in/ABA+v/BCqDrDzhdluAP6Z8TwAvs+XjboLIJWaZ5PQC4+zR3n+ruU4G7iKYY/Cbnrc2TytGTaG3aDFCS7nN/tpG5zwMqP9Bdy3OLKJkyHiuOq/xAP+UHzKySaErL1X0dOF2K4ANEI2Be9BuydGGSznk9S4A73H2xmV1iZpeEbmAhWKyIycd/CGAGUZ/7UwZcl6vyA50Sazaw/Y7fMfKiD3buqO8WqPwAROH6vrvv7HWH0YTO24EfuPuKLNuxz7KaT9bbvB53/3Ef616UzT5TJbBjypCduUtsyhx4mEXuPr+m+AAvWdqQafWluSw/AJDc2simH/yC+os/SvEB9UBUYiBf5QestZ1++jwUyw8cC5xnZv8O1AIdZtbq7jem378FeN3db8jUsVzTHR95MpDyAx3Nu9j4/f+m9rzTKZs+tWu5yg9kLj/g7id1+9hyA/CdzoCZ2bVADfDF3HcnM82Mzp+syw80PfoXkhu20LjgMRoXPAbAmH/8FEXVI1R+IMP2fTGziUTFTZcCL6T/8bnR3W/Nfdf2ppDlyUDKD9Se/W5qz95rdj2g8gOZtu+xzje7/dxA6E+oGeh0USQwhUwkMIVMJDCFTCQwhUwkMIVMJDCFTCQwhUwkMIVMJDCFTCQw3VaVJ72VHwBo+v2T7Hj0KSwWo/ywWdR9LLqrKPHWOrbcdi++qxVixrh/vgwrKVb5gSy271l+IL3sfODrRNNl1gJ/l6/nRmskywMzmwd8Jj5u9B7lB1qXvMGuF5cw/povMv47V1B9ZvRkH0+l2Pxfv6b+wnMZ/50rGHPlxRAvAlR+IMvt9yg/kJ5H9p/AO939UOAVojmSeaGQ5cds4GmLxfYoP7Djsaepft8pWHF0QlFUPQKA1kWvUzJpLCWTx0fLR1RisZjKDwyy/ADRzcEGVFp0C341e5fQCEani/mxCPh29/IDpVMn0L5+M23LVrL97oex4jh1HzuL0gMn0b5+M2BsuP4ndOxopuLYw6g565RBlR8oqhqStVT2KB/QOdeuh0xV0nrdvlv5gdOAL3du6O7tZnYpsBBoBl4HPpfbLvVNIctCU3LTQw+uv2nvAhy7lZnZc91+vyVdQAfoekj4dckNW27do/xARwcdzbsY+43PknizgU03/ZIJ3/0K3tFB6+srGfcv0eewDf9+KyVTJxArK+27BTkuP9CU3PTQgxtuHnSfzewRYGwv212VZRP2qfxA9wmrZlYMXEpUqWoF8EPga8C1WbZlnyhkWXD3M3Kwj73KD7Sv3UTFUXMxM0oPnISZ0bGjmXhdDWUzp3WNQuWHziSxcg2VJxyRt/ID+9rnoVR+AHgm3aY30se/g14+B4aiz2R50lv5gYoj59C65A0A2tdvwlMpYlWVlB0ynfaG9XS0JfBUirbX3qR4/BiVHxh8+YE1wBwzG53e5jSyK5CUExrJ8mev8gMjTp7Plp/cxdqrvo/Fi6j/9EcwM4oqK6g6/STWX30jmFF+6EwqDp8FoPIDGbbvi7uvNbOrgSfMrB1YBVwUoF+9skJdgSqbMMmnXHJFQY6drWX/fMXz/VWhGojeyg8MNasuujKnfRadLooEp5CJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppnReVRSkmTyxLzU0xy0VYVuwH5II5lIYAqZSGAKmUhgCplIYAqZSGAFu7roJU7b5EShDi+SNxrJRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAKLF+rAJSVJpk7cVKjDZ2VVoRsg+wWNZCKBKWQigSlkIoEpZCKBKWQigSlkIoEpZCKBKWQigSlkIoEpZCKBFey2quGoOt7K6WOXFLoZGf2p0A3YD2U1kpnZGWb2mpktN7Mre3n/b83slfTrL2Z2WO6bml9b//omwLzh1GcJo9+QmVkR8CPgTGAOcL6Zzemx2pvAKe5+KHANcEuuG5pPnurg9R8+CrCMYdJnCSebkewYYLm7r3D3BPAr4JzuK7j7X9x9W/rXp4GJuW1mfjW9tp7y8bUAieHSZwknm5BNAN7q9ntDellfPgX8bl8aVWiJzTspPaCq+6L9vs8STjYXPqyXZd7rimbvJPoLd2If718MXAz0/Es8tHiv3dvnPteMK89RA+XtJJuRrAGY1O33icDaniuZ2aHArcA57r6ltx25+y3uPt/d5xfXVgymvXlRMrqKto07ui/KSZ8r60qCtFeGtmxC9iww3cymmVkJ8HFgQfcVzGwycA/wCXdflvtm5lf1zLHsWrMdoGS49FnC6Tdk7p4ELgMeApYAd7j7YjO7xMwuSa/2z0A9cJOZvWRmzwVrcR5YUYyDP/8ugBkMkz5LOFl9Ge3uDwAP9Fj2424/fxr4dG6bVlj1xx4IsMjd53cu29/7LGHotiqRwBQykcAUMpHAFDKRwBQykcAUMpHAFDKRwBQykcAUMpHAFDKRwBQykcAUMpHAFDKRwBQykcAUMpHAFDKRwBQykcAUMpHAFDKRwBQykcAUMpHAFDKRwAr2fLLqeCunjV1aqMNn5fFCN0D2CxrJRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAJTyEQCU8hEAlPIRAIrWAXh4aiuqIXzql8odDMy+qdCN2A/pJFMJDCFTCQwhUwkMIVMJDCFTCQwhUwkMIVMJDCFTCQwhUwkMIVMJDCFTCQwhUwkMIVMJDCFTCQwhUwkMIVMJDCFTCSwgs2M1ixhGS40kokEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSmkIkEppCJBKaQiQSWVcjM7Awze83MlpvZlb28b2b2g/T7r5jZkblvan796Y+tAPOGU58ljH5DZmZFwI+AM4E5wPlmNqfHamcC09Ovi4Gbc9zOvEqlnG/9UxPAMoZJnyWcbEayY4Dl7r7C3RPAr4BzeqxzDvBzjzwN1JrZuBy3NW9eeamdyVOLABLDpc8STjYhmwC81e33hvSyga7ztrFhfYpx44u6L9rv+yzhxLNYx3pZ5oNYBzO7mOjUCqBt1uR1i7I4fiHUAdXAzG7L9vc+d5rZ/yoyENmErAGY1O33icDaQayDu98C3AJgZs+5+/wBtTZPzOx44JtAfXrRft/nTmb2XKHbsL/J5nTxWWC6mU0zsxLg48CCHussAC5IX3E7Dmh093U5bms+PUt0QaNkGPVZAul3JHP3pJldBjwEFAE/dffFZnZJ+v0fAw8AZwHLgRbgk+GaHF63Pt8DLGEY9FnCMfe9Pkbk58BmF6dPpYasXLdxOPZZChgykeFCt1WJBBY8ZEP9lqws2jfLzJ4yszYz+3KO9rnf9VkycPdgL6ILJW8ABwIlwMvAnB7rnAX8juh7p+OAZ0K2aRDtOwA4Gvg28GX1Wa+BvkKPZEP9lqx+2+fuG939WaA9V/tk/+uzZBA6ZEP9lqwQxx6OfZYMQocsZ7dkBRLi2MOxz5JB6JDl7JasQEIcezj2WTIIHbKhfktWNu0Lsc/9rc+SSegrK0RX0pYRXdG6Kr3sEuCS9M9GNCn0DWAhMD+fV36yaN9Yon/9m4Dt6Z+r1We9sn3pjg+RwHTHh0hgCplIYAqZSGAKmUhgCplIYAqZSGAKmUhgCplIYP8fbt8BmNvvWCQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 144x720 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "theme = \"Vertical cross-section\"\n",
    "print(theme)\n",
    "import matplotlib.tri as tri\n",
    "file_vtu = \"result_liakopoulos_t_7200.vtu\"\n",
    "m_plot=vtuIO.VTUIO(file_vtu, dim=2)\n",
    "triang=tri.Triangulation(m_plot.points[:,0],m_plot.points[:,1])\n",
    "u_plot = m_plot.get_point_field(\"gas_pressure\")\n",
    "p_plot = m_plot.get_point_field(\"displacement\").T[1]\n",
    "fig, ax = plt.subplots(ncols=2, figsize=(2,10))\n",
    "contour_left = ax[0].tricontourf(triang, u_plot)\n",
    "contour_right = ax[1].tricontourf(triang, p_plot)\n",
    "fig.colorbar(contour_left,ax=ax[0],label='$T$ / [K]')\n",
    "fig.colorbar(contour_right,ax=ax[1],label='$p$ / [MPa]')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "78afcf25",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Thu Aug  4 15:15:01 2022\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "print(time.ctime())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4beb4f2a",
   "metadata": {
    "tags": []
   },
   "source": [
    "**OGS links**\n",
    "- description: https://www.opengeosys.org/docs/benchmarks/\n",
    "- project file: https://gitlab.opengeosys.org/ogs/ogs/-/blob/master/Tests/Data/TH2M/H2M/Liakopoulos/liakopoulos_TH2M.prj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf4def87",
   "metadata": {},
   "outputs": [],
   "source": [
    "**Credits**\n",
    "- Wenqing Wang for set-up the OGS benchmark, https://www.opengeosys.org/docs/benchmarks/liquid-flow/liquid-flow-theis-problem/\n",
    "- Jörg Buchwald for ogs6py and VTUInterface, https://joss.theoj.org/papers/6ddcac13634c6373e6e808952fdf4abc"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}