{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# A simple heat problem for Python-Post-Processing\n", "\n", "This tutorial starts with a short description of the problem and how to simulate it with cfs. All needed files are provided. Afterwords some simple python postprocessing is done." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Short description of problem\n", "\n", "The geometry of the domain is a cylinder with the radius $r$ and the height $h$, due to symmetry reasons only on quarter is used for computation. The top volume is called V_all with the mantle surfaces S_mantle and the symmetry surfaces S_x and S_y. S_top is located on top of the cylinder and S_bottom on the bottom.\n", "\n", "|Sketch of problem|\n", "|:-:|\n", "|![](mesh.png)|\n", "\n", "S_bottom does have a prescribed temperature $T_{bot}$. We assume an ambient temperature of $T_{air}$ and a heat transfer coefficient $\\alpha$. Furthermore the domain moves with a certain velocity $v$ forward. On S_x and S_y the heat flux is zero, due to symmetry reasons. Because the temperature of the surrounding air is lower that the initial temperature of the cylinder, the cylinder cools down, while moving forward.\n", "\n", "| Description | Variable | Unit | Value |\n", "|:-----------:|:--------:|:----:|:-----:|\n", "|Density of iron|$\\rho$|kg/m³|7874|\n", "|Heat capacity of iron|$c$|J/(kg·K)|444|\n", "|Heat conductivity of iron|$k$|W/(m·K)|79.5|\n", "|Heat transfer coefficient|$\\alpha$|W/(m²·K)|20|\n", "|Heat source density|$\\dot{q}_v$|J/m³|50|\n", "|Bottom temperate|$T_{bot}$|K|293|\n", "|Air temperature|$T_{air}$|K|273|\n", "|Radius of cylinder|$r$|m|0.2|\n", "|Height of cylinder|$l$|m|1|\n", "|Velocity of cylinder|$v$|m/s|0.001|" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulation\n", "Download these files:\n", "\n", "- Journal file: [mesh.jou](mesh.jou)\n", "- Material file: [mat.xml](mat.xml)\n", "- Simulation file: [simulation.xml](simulation.xml)\n", " \n", "\n", "Create the mesh using mesh.jou and Trelis and run the simulation with cfs:\n", "\n", "- Terminal command for meshing: trelis -batch -nographics -nojournal UnitCube.jou\n", "- Terminal command for simulation: cfs simulation\n", "\n", "The simulation results should be in the ./results_hdf5/." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Postprocessing with pyhton\n", "For Python-postprocessing we are going to use the python-library [hdf5_tools.py](https://gitlab.com/openCFS/cfs/-/blob/master/share/python/hdf5_tools.py). This library is also enrolled automatic with every openCFS-version and can be found under /path/to/install/dir/CFS/share/python/hdf5_tools.py.\n", "\n", "In this tutorial we are going to use two functions from hdf5_tools.py:\n", "* get_result()\n", "* get_coordinates()\n", "\n", "How each function works, is described in the according docstring.\n", "\n", "Lets start with reading the nodal result (temperature) out of the CFS-file with get_result():" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import sys\n", "sys.path.insert(0, \"./Devel/CFS_SRC/CFS/share/python/hdf\")\n", "from hdf5_tools import get_result\n", "\n", "#Reading in the cfs-file, instert here the path to your cfs-file\n", "hdf5=f'./results_hdf5/simulation.cfs'\n", "\n", "#Reading the temperatue out of the cfs-file\n", "T=get_result(hdf5,\"heatTemperature\", region=\"V_all\", multistep=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After reading the CFS-file, the array T contains now the nodal temperature of all the nodes in the region V_all. Now we could simply search for the maximum and minimum temperature:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The maximal temperature is 20.0°C.\n", "The minimal temperature is 19.189°C.\n" ] } ], "source": [ "T_max=T.max()\n", "print(f'The maximal temperature is {np.round(T_max,3)}°C.')\n", "T_min=T.min()\n", "print(f'The minimal temperature is {np.round(T_min,3)}°C.')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The maximal temperature is 20°C and its clear, that it occurs on the S_bot since we prescribed the temperature there\n", "\n", "And where does the minimal temperature occurs? For this we going to use get_coordinates()." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "from hdf5_tools import get_coordinates\n", "# Getting the node, where the minimal temperature occurs\n", "Idx_min=np.argwhere(T==T.min())\n", "\n", "#Get the coordinates for each node\n", "X=get_coordinates(hdf5, region=\"V_all\")\n", "#Pluggin in the indices for maximal and minimal temperatures:\n", "X_min=X[Idx_min]\n", "# Reshaping X_min\n", "X_min=X_min.reshape(3)\n", "\n", "print(f'The coordinates for the minimal temperature are {X_min} ([x,y,z])')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Great, but actually i want them in polar coordinates. Well thats quite simple to achieve:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "In radial coordinates it is [ 2.00e-02 -1.35e+02 5.00e-02] ([r,phi,z]; phi in grad)\n" ] } ], "source": [ "#Quickly write a funtcion which convertes carthesian coord into clyindirc coords\n", "def cart2pol(X):\n", " r=np.sqrt(X[0]**2 + X[1]**2)\n", " phi=np.arctan2(X[1],X[0])\n", " phi=phi*360/(2*np.pi) #for degrees\n", " z=X[2]\n", " return [r,phi,z]\n", "\n", "print(f'In radial coordinates it is {np.round(cart2pol(X_min),2)} ([r,phi,z]; phi in grad)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Cool, that seems plausible.\n", "\n", "But how to i get the temperature distribution along the z-direction in the middle of the cylinder?\n", "\n", "For this we read out all the indices where x==0 and y==0 is zero, and then we could simply use the indices and get the according temperatures, right?\n", "(Because we are using symmetry for our mesh, there are actually nodes in the middle of the cylinder along the z-direction at x==0 and y==0, otherwise you have to interpolate or take the nearest node.)\n", "\n", " Lets give it a try:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'Temperature in °C')" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3gVZfr/8fedAiFUIQEpoReVDkEpQWBtiFgR1F2wrigqoq7u19+ufrfv2nZVwAVRAdvi2hVUrCBFEEKTIL0jaEBAitLv3x9n8BvxJDmEnJyUz+u65jpzZp6Zuceo95mZZ+7H3B0REZFjxcU6ABERKZ6UIEREJCwlCBERCUsJQkREwlKCEBGRsBJiHUBhSklJ8YYNG8Y6DBGREmPevHnb3D013LpSlSAaNmxIZmZmrMMQESkxzGx9but0i0lERMJSghARkbCUIEREJCwlCBERCUsJQkREwopagjCzNDObYmZLzWyJmQ0Lllc3sw/NbGXweVIu2/c2s+VmtsrM7o1WnCIiEl40ryAOAb9x91OBzsCtZnYacC/wsbs3Az4Ovv+EmcUDTwDnA6cBVwXbiohIEYnaexDuvgXYEszvNrOlQF3gYqBn0OxZYCrwP8dsfjqwyt3XAJjZS8F2X0Yj1lfnbWLU1FV0b5ZKRtMU2qRVJaVieeLiLBqHExEpEYrkRTkzawi0Bz4HagXJA3ffYmY1w2xSF9iY4/sm4Ixc9j0YGAxQv379AsX3u9cXc+DwEVZv3cv4z9b9bH2FxHi6NU2hd6uTubxjvQIdQ0SkpIl6gjCzSsBrwB3uvsssol/l4RqFHdnI3ccAYwDS09MLNPrRsr/0ZunXu5ixchszVm1j+sptP1n/w8HDfLT0Gz5a+g0jPllJk9RKNE6pSJOalULzqRWpUbEcEZ6biEiJENUEYWaJhJLDi+7+erD4GzOrHVw91Aayw2y6CUjL8b0esDlaccbFGS3rVKVlnarc1KPJT9bt3X+Itdv2sjJ7N4s2fsfWPftZs3UvM1dtY/+hIz+2q1ohkcapFTnl5Cr0aJ5Ct6YpVE5KjFbIIiJRZ9EactRCP6efBba7+x05lj8MfOvuDwS9k6q7+2+P2TYBWAGcBXwFzAV+6e5L8jpmenq6F1UtpiNHnM3f/cDqrXtZnb2HNdv2sDp7L1lffcfu/YdIjDc6NaxOrxY16XVKKk1SK+kKQ0SKHTOb5+7pYddFMUFkANOBxcDRn9q/I/Qc4mWgPrAB6O/u282sDvC0u/cJtu8DPAbEA2Pd/W/5HbMoE0RuDh4+wvz1O5iyfCtTl2ez7OvdANQ7qcKPyaJL4xQqlIuPaZwiIhCjBBELxSFBHOurnT8wdXk2U5ZtZeaqbfxw8DBJiXGc36o2/dPr0blRDfWWEpGYUYIoJvYfOsyctduZnPU1by/azO59h6hfPZn+HevRr2M96lSrEOsQRaSMUYIohvYdPMzkrK95OXMjn63+FjPo3iyVK9LTOPu0mpRP0C0oEYk+JYhibsO33/PqvI28Mm8TW77bR7XkRC5pV5eruzSgcWqlWIcnIqWYEkQJcfiIM2PVNl7O3MiHS77h4JEj9GlVmyE9m9CqbtVYhycipVBeCaJUDTla0sXHGT2ap9KjeSpbd+9n7My1vDBrPe8s3sKZzVMZ0qMJnRtXV3dZESkSuoIo5nbtO8jzs9YzbuZatu05QPv61bilZ1POOqWmej+JyAnTLaZSYN/Bw7ySuZEnp61h044faF6rEkN6NqFvmzokxmtYDxEpGCWIUuTQ4SNM+mILo6auZvk3u2mUUpHfnteC3q1O1q0nETlueSUI/fQsYRLi47ikfV3eG9adMYM6khBnDHlxPv1GfcbcddtjHZ6IlCJKECVUXJxxbsuTeW9Ydx7s15qvdv5A/9GzGPxcJquy98Q6PBEpBXSLqZT44cBhxs5cy6ipq/nh4GGu6JTGHWc3o2blpFiHJiLFmJ5BlCHf7tnPiE9W8cLs9ZRLiOPG7o258czGVCqvHs0i8nN6BlGG1KhUnj9e1JKP7upBrxY1efzjlfR6ZCpvLfyK0vRjQESiTwmilGqYUpEnftWB12/pSu2qSQx7aSG/fOpzVmXvjnVoIlJCKEGUch3qn8Qbt3TjL5e0Ysnm7zj/8ek88N4yvj9wKNahiUgxpwRRBsTHGYM6N+CTu3tycbu6jP50Nef8axqTs77WbScRyZUSRBmSUqk8j/Rvyys3d6FyUgI3vzCP68bPZf23e2MdmogUQ0oQZVCnhtWZNDSD+y44lblrt3POo9N47KMV7D90ONahiUgxogRRRiXEx/Hr7o355O6enNfyZB77aCV9h89gwYYdsQ5NRIoJJYgyrlaVJEZc1Z5x13Ziz/5D9Bv1GX9750t+OKCrCZGyTglCAOh1Sk0+uPNMrjy9Pk9NX8v5j09j9ppvYx2WiMSQEoT8qHJSIn+/tDX/ufEMjjhcOWY29725mN37DsY6NBGJASUI+ZmuTVKYfEd3bshoxIufb+C8R6cxdXl2rMMSkSIWtQRhZmPNLNvMsnIsa2tms8xssZlNNLMquWw7zMyyzGyJmd0RrRgld8nlEri/72m8NqQryeUTuHbcXO56eSHffa+rCZGyIppXEOOB3scsexq4191bA28A9xy7kZm1Am4ETgfaAn3NrFkU45Q8dKh/Eu/cnsFtvZry1sLNnPfYNKav3BrrsESkCEQtQbj7NODYEWxaANOC+Q+BfmE2PRWY7e7fu/sh4FPg0mjFKfkrnxDP3ee14I1bulIpKYFBz8zhf9/KUrkOkVKuqJ9BZAEXBfP9gbRc2pxpZjXMLBnok0s7AMxssJllmlnm1q36ZRtNbepVY9LQDG7IaMRzs9ZzwfAZzNd7EyKlVlEniOuBW81sHlAZOHBsA3dfCjxI6ApjMrAIyPWnqruPcfd0d09PTU2NTtTyo6TEeO7vexr/ufEMDhw6wuWjPuOfHyznwKEjsQ5NRApZkSYId1/m7ue6e0dgArA6l3bPuHsHdz+T0G2qlUUZp+Sva5MU3rujO5d1qMeIT1Zx2aiZrPhGpcRFSpMiTRBmVjP4jAPuA0bn064+cBmhZCLFTJWkRB7p35YnB3Vky8599B0xg6enr+HIEVWIFSkNotnNdQIwC2hhZpvM7AbgKjNbASwDNgPjgrZ1zOzdHJu/ZmZfAhOBW91dN7qLsfNansz7d55Jj+ap/PWdpfzq6c/Z8t0PsQ5LRE6QxqSWQuPuvJK5iT9OXEJifBwPXNaa81vXjnVYIpIHjUktRcLMGNApjXdu707DGskMeXE+9772hbrDipRQShBS6BqlVOTVIV25pWcT/pu5kb7DZ/DFpp2xDktEjpMShERFYnwcv+19Cv/5dWd+OHiYy/79GaOmrtYDbJESRAlCoqpLkxq8N6w757asxYOTl+kBtkgJogQhUVctuRxP/LIDD/Vrw6JNO+n92HQmZ22JdVgikg8lCCkSOR9gN6iRzM0vzOe+Nxez76BGrhMprpQgpEg1SqnIqzd3ZfCZjXlh9gYueWImq7L1BrZIcaQEIUWuXEIcv+tzKuOu68TW3fu5cMRMXp67kdL0To5IaaAEITHTq0VN3h3Wnfb1q/Hb175g2EsLNbypSDGiBCExVatKEs/fcAZ3n9ucdxZvoe8IvTMhUlwoQUjMxccZt/2iGf8d3JmDh47Qb9RnKvonUgwoQUixkd6wOu8O606vFjX56ztLueHZuWzf+7MhQ0SkiChBSLFSLbkcTw7qyJ8vbsnMVd/S5/HpzF137Mi1IlIUlCCk2DEzru7SkNdv6UpSYhxXjpnNv6eu0i0nkSKmBCHFVqu6VZk4NIPzW53MQ5OXc934uXy7Z3+swxIpM5QgpFirnJTIiKva89dLWjFrzbf0GT6dOWt1y0mkKChBSLFnZgzs3IA3bulKcrkErhwziyem6JaTSLQpQUiJ0bJOVd6+rRsXtKnDw+8v55pxc9imW04iUZNrgjCzu4JxpI9dPtTM7ohuWCLhVU5KZPiV7fj7pa35fO12+jw+nc/XfBvrsERKpbyuIK4Hng+zfEywTiQmzIxfnlGfN2/pRsXyCVz1lHo5iURDXgnC3f1nbym5+37AoheSSGROq1OFt2/rxvmta/PQ5OXc8OxcdujFOpFCk+czCDOrFckykVipnJTIyKva85fgxboLhk9n/oYdsQ5LpFTIK0E8DLxjZj3MrHIw9QQmAo8USXQiETAzBnVpyKtDuhAXZwwYPYtnZqxV+XCRE5RrgnD354D7gT8D64C1wJ+AP7j7s/nt2MzGmlm2mWXlWNbWzGaZ2WIzm2hmVXLZ9k4zW2JmWWY2wcySjvO8pAxqU68a7wztTq9TavKXSV8y5IX57FL5cJECy/MWk7u/5+493L2Gu6cE8+9FuO/xQO9jlj0N3OvurYE3gHuO3cjM6gK3A+nu3gqIB66M8JhSxlVNTmTMoI78vs+pfLj0G/oOn0HWV9/FOiyREimvbq7JZnZTMFU83h27+zTg2FdeWwDTgvkPgX65bJ4AVDCzBCAZ2Hy8x5eyy8y48czGvHxTZw4ePsJloz7jxc/X65aTyHHK6wriRWATsBH4TyEdLwu4KJjvD6Qd28DdvyL0jGMDsAX4zt0/yG2HZjbYzDLNLHPr1q2FFKaUBh0bVOed27vTuXENfv9GFne9vIjvDxyKdVgiJUZeCaIKsCaYqhXS8a4HbjWzeUBl4Gd9Es3sJOBioBFQB6hoZgNz26G7j3H3dHdPT01NLaQwpbSoXrEc46/txF3nNOfNhV9xyRMzWZW9J9ZhiZQIeSWIgcDVwHXAoMI4mLsvc/dz3b0jMAFYHabZ2cBad9/q7geB14GuhXF8KZvi4ozbz2rG89efwbY9B7h45AzeXqS7liL5yasX0xZ3/3/u/j/uvqEwDmZmNYPPOOA+YHSYZhuAzsEzEAPOApYWxvGlbMtolsI7t2dwSu0q3D5hAX94K4v9hw7HOiyRYitqxfrMbAIwC2hhZpuCuk5XmdkKYBmhB8/jgrZ1zOxdAHf/HHgVmA8sDmIcE604pWypXbUCLw3uzI3dG/HsrPUMeHI2m3Z8H+uwRIolK009O9LT0z0zMzPWYUgJMTlrC/e88gXx8cajV7SjV4uasQ5JpMiZ2Tx3Tw+3TuW+pczq3ao2E4dmULtqBa4bN5dH3l/OYRX8E/lRQiSNzKwr0DBn++BNa5ESrWFKRd64pSt/eGsJI6esYsHGHQy/sj01KpWPdWgiMZfvFYSZPU/ovYQMoFMwhb0cESmJkhLjefDyNjx0eRsy1+3gguEzmLdew5qKRHIFkQ6c5qXpYYVIGAPS02hZpwpDXpjPFU/O5nd9TuW6bg0JdaYTKXsieQaRBZwc7UBEioOWdaoycWgGvU6pyZ8nfcltExawZ7/evpayKZIriBTgSzObA/w4ALC7X5T7JiIlV9UKoYJ/T05bw0OTl7F0yy5GD+xI81qVYx2aSJGKJEH8MdpBiBQ3ZsbNPZrQtl41hk5YwMUjZ/JAv9Zc3K5urEMTKTJ6D0IkH9m79nHbfxYwZ912BnVuwH19T6V8QnyswxIpFAV6D8LMZgSfu81sV45pt5ntilawIsVNzSpJvHjjGQw+szHPz17PFU/OZvPOH2IdlkjU5VWLKSP4rOzuVXJMld097EhwIqVVYnwcv+tzKqMHdmBV9h4uGD6d6StVXl5KN71JLXIcereqzdu3daNm5SSuHjuHER+v5IjevpZSSglC5Dg1Tq3EG7d25ZJ2dfnnhyu44dm57Pz+Z0ObiJR4ShAiBZBcLoF/DWjLXy5pxYxV2+g7YgaLN2nsayldIkoQZtbAzM4O5iuYmTqES5lnZgzq3IBXbu7KkSNOv9Gf8dKcDRr7WkqNSGox3UhofIYng0X1gDejGZRISdIurRqTbu/OGY2qc+/ri7nn1S/Yd1ADEUnJF8kVxK1AN2AXgLuvBFQ4XySH6hXLMf6607n9F015dd4mLv33Z6z/dm+swxI5IZEkiP3u/uMTODNLAHQNLXKM+DjjrnNbMO7aTmze+QN9R8zg46XfxDoskQKLJEF8ama/AyqY2TnAK8DE6IYlUnL1OqUmk4ZmUL96Mjc8m6mBiKTEiiRB/A+wldD40DcB7wL3RTMokZIurXoyrw3pyhXpaYycsoprxs7h2z37899QpBjJsxaTmcUBX7h7q6ILqeBUi0mKo//O3cD9by0hpWI5nvhVB9rXPynWIYn8qMBjUrv7EWCRmdWPSmQiZcAVnerz+pCuxMUZA56cxfOz1qkrrJQIkdxiqg0sMbOPzezto1O0AxMpTVrVrcqkoRlkNE3h/reWcNfLi/j+gAYikuItkvEg/hT1KETKgGrJ5Xjmmk6MnLKKRz9awdItuxg1sCONUirGOjSRsPK9gnD3T8NN+W1nZmPNLNvMsnIsa2tms8xssZlNNLOfVYU1sxZmtjDHtMvM7jj+UxMpfuLijNvPasb4607n6137uGjEDD5Y8nWswxIJK5I3qXOOB7HPzA5HOB7EeKD3McueBu5199bAG8A9x27k7svdvZ27twM6At8HbUVKjR7NU5k0NINGqRUZ/Pw8HnhvGYcOH4l1WCI/EckVRM7xIJKAfsDICLabBmw/ZnELYFow/2Gwr7ycBax29/X5HU+kpKl3UjIv39SFq06vz+hPV3P12DlsU1dYKUaOu5qru78J/KKAx8sCLgrm+wNp+bS/EpiQVwMzG2xmmWaWuXWrBnCRkiUpMZ5/XNaahy9vw7z1O+g7fAbz1u+IdVgiQGS3mC7LMV1uZg9Q8FIb1wO3mtk8oDKQaxF9MytHKJm8ktcO3X2Mu6e7e3pqamoBwxKJrf7pabx+S1cSE4wrx8zi2c/UFVZiL5JeTBfmmD8ErAMuLsjB3H0ZcC6AmTUHLsij+fnAfHdXMRspE1rWqcqk27pz58sL+cPbS1iwYQd/v6w1yeUi+c9UpPBF8m/e0+4+M+cCM+sGZB/vwcysprtnB29o3weMzqP5VeRze0mktKmanMjTV6fzxJRV/OujFSzdsptRAzvQOLVSrEOTMiiSZxAjIlz2E2Y2AZgFtDCzTWZ2A3CVma0AlgGbgXFB2zpm9m6ObZOBc4DXI4hPpFSJizOGntWMZ687nezd+7h45EwmZ6krrBS9XGsxmVkXoCtwB/BojlVVgEvdvW30wzs+qsUkpc2mHd9zy4vz+WLTd9zcowl3n9uchHiNFCyFp6C1mMoBlQjdhqqcY9oFXF7YQYrIz9U7KZlXbu7CL88IdYUd9Iy6wkrRybOaK4TGoy4p7yHoCkJKs1cyN3Lfm1mclByqCtuxgarCyokrcDXXwPdm9rCZvWtmnxydCjlGEcmHusJKUYskQbxI6KFyI0KF+9YBc6MYk4jk4mhX2O7NUvnD20u4878LVRVWoiaSBFHD3Z8BDgaF+q4HOkc5LhHJxdGusL85pzlvLdrMpU98xtpte2MdlpRCkSSIg8HnFjO7wMzaA/WiGJOI5OPYrrAXjZjB+6oKK4UskgTxVzOrCvwGuJtQRdY7oxqViETkzOapTAyqwt70/DwenKyqsFJ48kwQZhYPNHP379w9y917uXtHd9eIciLFRM6qsKOmqiqsFJ78xqQ+zP9VXxWRYupoVdiHgqqwF46YwYINqgorJyaSW0yfmdlIM+tuZh2OTlGPTESO24D0NF4b0pWEeGPAk7N4fpa6wkrBRfKi3JQwi93dCzomRNToRTmRkO++P8gd/13AlOVbubR9Xf5+aWsqlIuPdVhSDOX1oly+1VzdvVfhhyQi0VQ1OZFnrunEyCmrePSjFSzdsovRAzvSMKVirEOTEiSSAYNqmdkzZvZe8P20oDKriBRjcXHG7Wc1Y/x1p/P1rn1cOHIGH36p4VUkcpE8gxgPvA/UCb6vIFThVURKgB7NU5l4WwYNa1TkxucyeWjyMg4f0XMJyV8kCSLF3V8GjgC4+yHgcFSjEpFClVY9VBX2qtPT+PfU1Vw99nO+VVdYyUckCWKvmdUgGIfazDoD30U1KhEpdKGusG146PI2zF23g77qCiv5iCRB3AW8DTQxs5nAc8DQqEYlIlEzID2N14d0JT4u6Ao7e726wkpY+SYId58P9CA0utxNQEt3/yLagYlI9LSqW5VJQzPo1jSF+9/M4jcvL+KHA7pzLD8VSS+mJOB24C+Eyn3fGiwTkRKsWnI5xl7TiTvObsYbC7/i0n/PZJ2qwkoOkdxieg5oCYwARgKnAc9HMygRKRpxccYdZzdn3LWd2PKdusLKT0WSIFq4+w3uPiWYBgPNox2YiBSdni1qMmloBg1qJHPjc5k8/L66wkpkCWJB0HMJADM7A5gZvZBEJBbSqifz6s1duSI9jSemrOaasXPUFbaMiyRBnEGoYN86M1sHzAJ6mNliM9PDapFSJCkxngcvb8OD/VozZ912Lhwxg4Ubd8Y6LImRSBJEb0LjUfcIpkZAH6AvcGFuG5nZWDPLNrOsHMvamtmsILlMNLMquWxbzcxeNbNlZrbUzLocz0mJyIm5olN9Xru5K3FxRv/Rn/GCusKWSZF0c10P7AKqAjWOTu6+PliXm/GEkktOTwP3untr4A3gnly2fRyY7O6nAG2BpfnFKSKFq3W9UFfYrk1SuO/NLH7zirrCljWRlPv+C3AtsJrgbWoiLPdtZg2BSe7eKvi+C6jq7m5macD77n7aMdtUARYBjf04f7Ko3LdI4TtyxHn845UM/2QlLWpV5slBHWlQQ1VhS4u8yn1HcotpANDE3XsGQ472OoGxILL4vxHq+gNpYdo0BrYC48xsgZk9bWa5/ttoZoPNLNPMMrdu3VrAsEQkN3Fxxp3nNGds0BW274gZfKSusGVCJAkiC6hWSMe7ntCLdvOAysCBMG0SgA7AKHdvD+wF7s1th+4+xt3T3T09NTW1kMIUkWP1ytEV9tfPZfLI+8vVFbaUy3fAIOAfhLq6ZgE/9nlz9+Meq9rdlwHnAphZc+CCMM02AZvc/fPg+6vkkSBEpOgc7Qr7h7eWMHLKKhZt2snjV7anesVysQ5NoiCSBPEs8CCwmKDkd0GZWU13zzazOOA+YPSxbdz9azPbaGYt3H05cBbw5YkcV0QKz9GusB0aVOP+t5bQd/h0/j2wI+3SCutGgxQXkdxi2ubuw4O3qD89OuW3kZlNIPTORAsz2xSMQneVma0AlgGbgXFB2zpm9m6OzYcCLwbvWbQD/n6c5yUiUaausKVfJL2Y/kXo1tLb/PQW0/zohnb81ItJpOjt/P4Aw15ayKcrtnJZh7r87ZLWVCgXH+uwJEJ59WKK5BZT++Czc45lDhS0J5OIlCLVkssx7tpODP9kJY9/vJIvN+9SV9hSIt8riJJEVxAisTV1eTbDXlrIEXceHdCOs0+rFeuQJB8n9B6EmdUys2fM7L3g+2nB8wQRkZ/oeUxXWFWFLdkieUg9HngfqBN8XwHcEa2ARKRkO9oV9spOqgpb0uWaIMzs6POJFHd/maCLq7sfAlSQRURylZQYzwP92vBQvzY/VoVdsGFHrMOS45TXFcSc4HOvmdUgqMMUjA3xXbQDE5GSb0CnNF4fEuoKO+DJWTyvrrAlSl4JwoLPuwh1cW1iZjMJDUE6NNqBiUjp0KpuqCpsRtMU7n8zi9+8rKqwJUVe3VxTzeyuYP4N4F1CSWM/cDagwYJEJCLVksvxzDWdGPHJKh77eAVfbtnF6IEdaZiirrDFWV5XEPFAJUJF9SoSSibxQHKwTEQkYnFxxrCzmzHu2k58vWsfF46YwQdLvo51WJKHXN+DMLP57t6hiOM5IXoPQqRk2Lj9e255cT6Lv/qOIT2b8JtzmpMQH0mnSilsBX0PwvJYJyJSYGnVk3nl5i5cdXp9Rk1dzdVj57BNXWGLnbwSxFlFFoWIlDlJifH847LWPHx5G+at30Hf4TOYt15dYYuTXBOEu28vykBEpGzqn57G67d0pVxCHFc8OYvxM9eqK2wxoZt+IhJzLetUZeLQDHq2SOWPE79k2EsL+f7AoViHVeYpQYhIsVC1QiJjBqVzz3ktmPTFZi55YiZrtu6JdVhlmhKEiBQbcXHGrb2a8tz1Z7BtzwEuGjmT9xZviXVYZZYShIgUOxnNUpg0NIOmNSsx5MX5/O2dLzl4+IRGPJYCUIIQkWKpTrUKvHxTF67p0oCnpq/lV099TvaufbEOq0xRghCRYqtcQhx/urgVj13RjsVffccFI2bw+ZpvYx1WmaEEISLF3iXt6/Lmrd2oXD6BXz79OU9NW6OusEVACUJESoQWJ1fmrdu6cc6ptfjbu0u55cX57N53MNZhlWpKECJSYlROSmTUwA78vs+pfPDlN1w8cibLv94d67BKLSUIESlRzIwbz2zMi78+g137DnHJEzN5a+FXsQ6rVIpagjCzsWaWbWZZOZa1NbNZZrbYzCaaWZVctl0XtFloZirPKiI/07lxDd69PYPWdasy7KWF3P9mFvsPaSCiwhTNK4jxQO9jlj0N3OvurQkNQnRPHtv3cvd2uZWhFRGpWSWJF288g8FnNub52esZ8ORsvtr5Q6zDKjWiliDcfRpwbMG/FsC0YP5DoF+0ji8iZUNifBy/63Mqowd2YHX2HvoOn860FVtjHVapUNTPILKAi4L5/kBaLu0c+MDM5pnZ4Lx2aGaDzSzTzDK3btW/FCJlVe9WtXn7tm7UrJzENePm8PhHKzlyRF1hT0RRJ4jrgVvNbB6hYUsP5NKuWzCa3flB+zNz26G7j3H3dHdPT01NLfyIRaTEaJxaiTdu7cql7ery6EcruP7ZuezYm9v/ZiQ/RZog3H2Zu5/r7h2BCcDqXNptDj6zCT2rOL3oohSRkiy5XAL/HNCWv13ais9WfUvfETNYtHFnrMMqkYo0QZhZzeAzDrgPGB2mTUUzq3x0HjiX0K0pEZGImBm/OqMBrw7pAkD/0bN4YfZ6vX19nKLZzXUCMAtoYWabzOwG4CozWwEsAzYD44K2dczs3WDTWsAMM1sEzAHecffJ0YpTREqvNvWqMWloBl2b1uC+N7O4878aiOh4WGnKqOnp6Z6ZqdcmROSnjhxxnpiyin99tIKmqZUYNbAjTWtWinVYxbDJgP0AAAwJSURBVIKZzcvtdQK9SS0ipV5cnDH0rGY8f/0ZbN97gItGzuDtRZtjHVaxpwQhImVGRrMU3rm9O6fVrsLtExbwv2/p7eu8KEGISJlyctUkJgzuzK8zGvHcrNDb15t2fB/rsIolJQgRKXMS4+O4r+9pjB7YgTXZe+g7YgZTlmfHOqxiRwlCRMqs3q1q8/bQDE6uksR14+byyPvLOaSxr3+kBCEiZVqjlIq8eWs3BqTXY+SUVQx8RmNfH6UEISJlXlJiPA9d3pZH+rdl4cad9Bk+g5mrtsU6rJhTghARCVzesR5v3ZpB1QoJDHzmcx77aAWHy3DBPyUIEZEcWpxcmbdvy+CSdnV57KOVXDN2Dlt37491WDGhBCEicoyK5RP414C2PHBZa+as284Fw6cze823sQ6ryClBiIiEYWZceXp93rylGxXLJ/DLp2bzxJRVZWqMCSUIEZE8nFanChOHZnBBmzo8/P5yrhs/l217ysYtJyUIEZF8VCqfwPAr2/GXS1oxa823nP/4dGasLP29nJQgREQiYGYM6tyAN2/pRpWkBAaN/ZwHJy/jYCl+sU4JQkTkOBy95XRFehqjpq6m/+hZbNxeOms5KUGIiByn5HIJPNCvDSOuas/q7D30eXw6E0th+XAlCBGRArqwbR3eHdadprUqMXTCAv7n1S9K1Yh1ShAiIicgrXoyL9/UhVt7NeHleRu5cMQMvty8K9ZhFQolCBGRE5QYH8c9553CCzecwe59h7jk3zN58tPVJb5MhxKEiEgh6dY0hfeGdadn81T+8d4y+o/+jNVb98Q6rAJTghARKUQ1KpXnyUEdefzKdqzeupc+j0/n6elrSuTVhBKEiEghMzMubleXD+88k+7NUvjrO0u54slZrN22N9ahHRclCBGRKKlZJYmnrk7nn/3bsuKb3Zz/+DTGzlhbYuo5RS1BmNlYM8s2s6wcy9qa2SwzW2xmE82sSh7bx5vZAjObFK0YRUSizczo17EeH9zZgy6Na/DnSV9y5ZjZrCsBVxPRvIIYD/Q+ZtnTwL3u3hp4A7gnj+2HAUujE5qISNE6uWoSY6/txMOXt2Hpll2c//h0Rn6ysli/NxG1BOHu04DtxyxuAUwL5j8E+oXb1szqARcQSigiIqWCmdE/PY0P7jqTjGYpPPLBCno8PJXnZ68vljWdivoZRBZwUTDfH0jLpd1jwG+BfP+JmdlgM8s0s8ytW7cWTpQiIlFUu2oFnro6nVdu7kKD6snc/2YW5/zrUyYu2lysnk8UdYK4HrjVzOYBlYEDxzYws75AtrvPi2SH7j7G3dPdPT01NbVwoxURiaJODavzys1deOaadMonxDN0wgIufmIm01cWjx+7RZog3H2Zu5/r7h2BCcDqMM26AReZ2TrgJeAXZvZCEYYpIlJkzIyzTq3Fu8O688/+bdm+9wCDnpnDr56ezRebdsY2NvfoXc6YWUNgkru3Cr7XdPdsM4sj9BB7qruPzWP7nsDd7t43kuOlp6d7ZmbmiYYtIhIz+w8d5oXZG3hiyiq27z1AeoOT6NjgJNrXr0b7+idRq0pSoR7PzOa5e3q4dQmFeqSfHnQC0BNIMbNNwB+ASmZ2a9DkdWBc0LYO8LS794lWPCIiJUH5hHhuyGjEgPR6jJu5jk+WZTN25loOTgv9mK9brQLt6lejQ/1Q0mhZpwrlE+KjEktUryCKmq4gRKQ02nfwMF9u2cWCDTuZv2EHCzfs5KudPwBQLj6OtmlV+e/gLsTF2XHvOyZXECIiUjiSEuPpUP8kOtQ/iRtoBMA3u/axYMNOFmzYwa59BwuUHPKjBCEiUgLVqpJE71Yn07vVyVE7hmoxiYhIWEoQIiISlhKEiIiEpQQhIiJhKUGIiEhYShAiIhKWEoSIiISlBCEiImGVqlIbZrYVWB/rOI5TCrAt1kEUMZ1z2aBzLhkauHvYsRJKVYIoicwsM7c6KKWVzrls0DmXfLrFJCIiYSlBiIhIWEoQsTcm1gHEgM65bNA5l3B6BiEiImHpCkJERMJSghARkbCUIIqAmVU3sw/NbGXweVIu7Xqb2XIzW2Vm94ZZf7eZuZmlRD/qE3Oi52xmD5vZMjP7wszeMLNqRRd95CL4m5mZDQ/Wf2FmHSLdtrgq6DmbWZqZTTGzpWa2xMyGFX30BXMif+dgfbyZLTCzSUUXdSFwd01RnoCHgHuD+XuBB8O0iQdWA42BcsAi4LQc69OA9wm9CJgS63OK9jkD5wIJwfyD4baP9ZTf3yxo0wd4DzCgM/B5pNsWx+kEz7k20CGYrwysKO3nnGP9XcB/gEmxPp/jmXQFUTQuBp4N5p8FLgnT5nRglbuvcfcDwEvBdkc9CvwWKCm9Ck7onN39A3c/FLSbDdSLcrwFkd/fjOD7cx4yG6hmZrUj3LY4KvA5u/sWd58P4O67gaVA3aIMvoBO5O+MmdUDLgCeLsqgC4MSRNGo5e5bAILPmmHa1AU25vi+KViGmV0EfOXui6IdaCE6oXM+xvWEfp0VN5HEn1ubSM+9uDmRc/6RmTUE2gOfF3qEhe9Ez/kxQj/ujkQrwGhJiHUApYWZfQSEGz3895HuIswyN7PkYB/nFjS2aInWOR9zjN8Dh4AXjy+6IpFv/Hm0iWTb4uhEzjm00qwS8Bpwh7vvKsTYoqXA52xmfYFsd59nZj0LPbIoU4IoJO5+dm7rzOybo5fYwWVndphmmwg9ZziqHrAZaAI0AhaZ2dHl883sdHf/utBOoACieM5H93EN0Bc4y4MbucVMnvHn06ZcBNsWRydyzphZIqHk8KK7vx7FOAvTiZzz5cBFZtYHSAKqmNkL7j4wivEWnlg/BCkLE/AwP31g+1CYNgnAGkLJ4OiDsJZh2q2jZDykPqFzBnoDXwKpsT6XPM4x378ZoXvPOR9ezjmev3dxm07wnA14Dngs1udRVOd8TJuelLCH1DEPoCxMQA3gY2Bl8Fk9WF4HeDdHuz6EenasBn6fy75KSoI4oXMGVhG6p7swmEbH+pxyOc+fxQ/cDNwczBvwRLB+MZB+PH/v4jgV9JyBDEK3Zr7I8XftE+vzifbfOcc+SlyCUKkNEREJS72YREQkLCUIEREJSwlCRETCUoIQEZGwlCBERCQsJQgp08xsqpmlB/PvFkbVWDO7xMxOy/H9z2aW60uFx7HfdDMbfqL7EYmUurlKmWZmU4G73T0zl/VG6L+TiOvomNl4Qv3dXy2UIEViRFcQUmqY2c1mtjCY1prZlDBtKpjZS0HN/v8CFXKsW2dmKWbWMBiz4N/AfCDNzO4xs7nBdn/Ksc3VwbJFZva8mXUFLgIeDuJoYmbjzezyoP1ZwbgAi81srJmVz3HsP5nZ/GDdKWFi73l0PAEz+2Ow/VQzW2Nmt+fyz2SPmT1oZvPM7CMzOz3HNhed0D9wKfWUIKTUcPfR7t4O6ESoNs6/wjQbAnzv7m2AvwEdc9ldC0Llm9sH880IlX1uB3Q0szPNrCWhwoS/cPe2wDB3/wx4G7jH3du5++qjOzSzJGA8cIW7tyZUwmFIjmNuc/cOwCjg7ghO+RTgvCCuPwR1jo5VEZjq7h2B3cBfgXOAS4E/R3AMKcOUIKQ0ehz4xN0nhll3JvACgLt/QajsQzjrPVTXH0KVdM8FFhC6ojiFUML4BfCqu28L9rc9n7haAGvdfUXw/dkgnqOOFq+bBzTMZ18A77j7/uD42UCtMG0OAJOD+cXAp+5+MJiP5BhShilBSKliZtcCDYA/Bd8vzXHbKT1oFsmDt705dwv8I7giaOfuTd39mWD58TzEC1cSOqf9wedhIqu0vD/HfG7bHPT/e9B45Og2wTMVVXOWPClBSKlhZh0J3ZoZePShsru/keN/7JnANOBXQftWQJsIdv0+cH0wjgFmVtfMahIqQjjAzGoEy6sH7XcTGlLzWMuAhmbWNPg+CPi0AKcqUiSUIKQ0uQ2oDkwJrhjCDfE4CqhkZl8QGuVrTn47dfcPCI0nPMvMFgOvApXdfQmh5xifmtki/u+Zx0vAPcHD6CY59rMPuA54JdjPEWB0Ac9VJOrUzVVERMLSFYSIiISlBCEiImEpQYiISFhKECIiEpYShIiIhKUEISIiYSlBiIhIWP8feTUBAATsBNAAAAAASUVORK5CYII=\n", "text/plain": [ "