examples/notebooks/Poisson_non_periodic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving Poisson's equation on non-periodic topological domains\n",
    "In this notebook, we will show how to solve poisson's equation on single patch and multipatch domains. As we are running this in a notebook, everything will be run in serial and hence we are limiting ourselves to a fairly coarse discretization to avoid taking too much time. However, Psydac allows for hybrid MPI + OpenMP parallelization with barely any changed to the code. The lines that are impacted by that will be preceeded by their MPI equivalent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1 : Building the domain\n",
    "Psydac uses the powerful topological tools of SymPDE to build a large variety of domains. Here we show how to do a quarter annulus using a Square and a polar mapping and how to join two annulus to create a multipatch domain. Left commented is an example of a complex multipatch domain built using SymPDE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMEAAAD4CAYAAABVN4L5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOydd3hcV5n/P2eqRqPeu2R1yarultxt2U4jYYElIWxYFjYLCcuGBZ5QAhsgsIEfpMCySQgkhA2EBVLtOO523GVJlqxmWdao997rzNzfH/K9qzJzJceOI9v6Ps88ku49995zR+c95z1v+b5CkiQWsYhbGZqPugOLWMRHjUUhWMQtj0UhWMQtj0UhWMQtj0UhWMQtD91H3YEPAj8/PykqKuqj7sYibiAUFBR0SpLk7+jcDSkEUVFR5Ofnf9TdWMQNBCFEnbNzi+rQIm55LArBIm55LArBIm55LArBIm55LArBIm55XBMhEEK8JIRoF0KUOjkvhBC/FEJUCSGKhRDLppzbKYS4ePnct65FfxaxiCvBtVoJfg/sVDl/GxB3+fMg8ByAEEIL/Pry+WTgPiFE8jXq0yIWMS9cEz+BJEnHhBBRKk3uBv4gTcZtnxFCeAkhgoEooEqSpGoAIcSfL7ctv9I+bN++nfr6enJyctDr9ej1elxcXPD09MTHxwcfHx/8/f2JjIwkKCgIjWZRE1zEJK6XsywUaJjyd+PlY46Or3Z0AyHEg0yuIkRERMw6f+HCBRobG7l48eKcndHpdPj4+BAQEEBkZCTR0dHExcWRkZHB8uXLcXV1ne97LeImwPUSAuHgmKRyfPZBSfoN8BuAFStWzGrj6+uLj48PeXl5jI2NMT4+zuDgIN3d3XR3d9PT00NbWxtNTU20tLQov58+fZp3331XuY9WqyUsLIzIyEiWL19OTk4OGzZswGw2f6AXX8TCx/USgkYgfMrfYUAzYHBy/ANBCIHBYMBgMACTghEZGTnndf39/ZSUlFBUVERRURGlpaWUlZVx7Ngxnn76aXQ6HbGxsWRlZbFz50527tyJu7v7B+3mIhYYrpcQvAN85bLOvxrokySpRQjRAcQJIZYATcC9wGeuU58UeHh4kJ2dTXZ2NgDV1dXk5+eTlJTE2bNnOXnyJHl5ebzyyiu89NJL6HQ6UlNT2blzJ/feey9paWnXu8uLuJaQJOmqP8BrQAswweSs/wXgS8CXLp8XTFqBLEAJsGLKtbcDlZfPfXc+z1u+fLk0E+np6VJ6evqs4x8E+fn50uuvvy7Z7fZpx7u7u6VXX31V+od/+AcpMjJSYlJ1kyIiIqQHH3xQOn369DV5/iKuPYB8ycl4EtINmGi/YsUKaWYUaUZGBgBFRUVXff+DBw+i1WrZvHmzaruSkhL+/Oc/895771FUVIQkSURFRXHPPffw0EMPERcXd9V9WcS1gRCiQJKkFY7OLdoJZ8But9PX14e3t/ecbVNTU/nxj3/MuXPnsFgsfO9738NsNvPMM8+QkJBAVlYWv/nNbxgZGbkOPV/EB8WiEMxAf38/NpttXkIwFUuWLOGHP/whpaWlFBcX88UvfpGLFy/yL//yL4SGhvKv//qvNDQ0zH2jRVx33DTqUFpaGjabjZ///OdoNBo0Gg1arVZxnBkMBlxcXHBxccFkMuHq6orZbEanm24bqKmpIS8vj507d+Lh4XFV/RwdHeWHP/whe/fupbCwEL1ez2233cZjjz3GypUrr+rei7gyqKlDN2RmmTMIIfD09ESSJOx2O1arlZGREfr7+xkfH2diYmLWNSaTCTc3Nzw9PfH09KS1tRWtVoubm9tV96ejo4OMjAy+8pWv0NDQwM9+9jN2797Nrl272Lx5M48//jjr16+/6ucs4upw0wiBPPtnZWU5bWOz2RgdHWVkZITh4WEGBwcZHBxkYGCAurq6aUJy8OBBfH198fX1JSAg4AN5kauqqjCbzQQFBRESEsLrr79OU1MTP/jBD/jjH//Ihg0bWLduHU8++aRinl3E9cdNow5drXVIkiQGBwfZt28fXl5e6PV6uru7FcFwc3MjICCA4OBgAgMDZ6lRM9HX18e+fftITU0lKSlp1vn29nZ+9KMf8fLLLzM8PMzOnTv5+c9/TnLyYvzgh4FbRh26GgghsNvt2O124uLiiIyMRJIk+vr6aG9vp62tjfr6eqqrq9FoNAQGBhISEkJYWBhGo3HW/SwWCxqNhujoaIfPCwgI4Fe/+hXf+c53+Pa3v81rr71Geno6999/P7/4xS/w9fX9sF95EZexaB2agp6eHgDFMiSEwMvLi/j4eNavX88999zDxo0biYmJob+/n4KCAt555x2OHz9ObW0tVqsVgImJCWprawkPD3coIFMRHBzM73//e0pLS9mxYwevvPIKcXFxPPXUU9hstg/3hRcBLArBNHR3d6PT6ZxuiuUVIDMzk9tvv52cnBzi4+Pp6+vj7NmzvPPOO+Tn51NeXo7VaiUmJmbez46Li2P37t3s37+fwMBAvv71r5ORkUFeXt61er1FOMGiEExBT08PXl5e88o1EELg7e1Neno6d9xxB5s3byY0NJS6ujouXryIVqtlcHDwimfznJwcSktLefLJJ6mtrSUrK4uHH36YgYGBD/pai5gDN40QyHEgPT099PX1MTAwwMjICFarlfls/u12O729vVfsJINJgfD392f16tWKdUqv13P27FneffddysvLGR8fn/f9tFotjz76KGVlZWzevJn//u//ZunSpezfv/+K+7aIuXHTbIwlScJms3HgwIFZ5zQaDUajUXGWubq64ubmhru7O25ubri5uTEwMPCBPMUzUVdXpzjFOjs7qayspLS0lIsXLxIXF0dcXNyc+wQZERERvP322zzxxBM899xz3H777fzrv/4rP/vZz9Dr9VfVz0X8H24aIRBCoNVqyc7OVqw8ExMTymd0dJTR0VGGh4fp7Oyc5hPQarXKwBwcHKS3txcPD48rTsEcHR2lsbGRmJgY9Ho9wcHBBAcH09PTQ3l5OeXl5VRWVhIfH09CQsKcA1mSJAoKCsjMzOTMmTN8+ctf5plnnuHgwYP89a9/JTEx8cq/qEXMwk0lBEIIQkND52wrSZKSedbf309fXx+NjY0AymDV6XT4+/sTEBBAYGAgnp6eCOEoEe7/UF1djd1un7Uh9vb2Jjs7m76+PsrKyigvL8disZCUlERsbKxTYauvr6e5uZm0tDTi4+M5cOAAv/zlL/nOd77DypUr+fWvf80DDzwwz29oEc5w0wjBlUAIgdFoxGg0Kvb47u5uXFxcWL16Nd3d3XR0dNDe3k5LSwsArq6uBAcHExoaSkBAwKyBa7fbqa6uJiAgwGnMkaenJ1lZWXR1dSmZbBaLhczMTIKCgqa1HRkZobCwEF9fX+Lj44FJte6RRx5h48aNfPKTn+Rzn/scR48e5fnnn1ey6RZx5bglhWAm5E1xVFQU7u7uuLu7K2mZQ0NDtLW10dzcTG1tLRaLBYPBQFhYGBEREfj7+yOEoKWlheHhYcVzrQZfX182btxIS0sLRUVFHDt2jJCQEDIzMzGbzYoaZLPZWLly5SyBy8zMpKioiM997nO8/PLLlJSUsHv3bgIDAz+U7+dmxzURAiHETuBZQAv8VpKkJ2ec/yZw/5RnJgH+kiR1CyFqgQHABlidubY/TAwODmK1WvHx8Zl1zmw2Ex0dTXR0NFarlba2NhoaGqirq6O6uhqz2cySJUtoa2vDZDIREhIyr2cKIQgJCSEwMJDKykouXLighFnodDqam5tJT093uqq4u7vzxhtv8NOf/pTvfve7LFu2jLfeemsxOvUD4KqFYAqBVg6TqZV5Qoh3JElSuIMkSfp/wP+73P4u4GuSJHVPuc1mSZI6r7YvHxQzPcXOoNPpCA0NJTQ0FKvVSlNTEzU1NZSWThLveXh40NXVhZ+f35z7BxlarZakpCQiIiIoKCigsLBQ8VTPJzPt0UcfJSUlhfvvv59Nmzbx0ksv8elPf3pez17EJK6Fn2AVlwm0JEkaB2QCLWe4j8mc5AWD7u5utFrtFTFI6HQ6IiMj2bRpE3LVnJGREY4cOcLBgwdpaGiYl39ChtlsZt26dXh5eSFJEv39/Vgslnnd44477uDMmTMEBgZy//3389RTT837uYu4NurQlRBouTJJ1/iVKYclYL8QQgJekCb5hRxdq0q+JecQnD9/XjGXTk2oMZlMip9gpo59JZ7imbBarTQ3NxMWFsaqVauora3l0qVLnD59Gg8PD5KSkggPD5/Xvevr6+nt7SU5OZnu7m4KCwtpaWlh1apVuLi4qF6bmJhIbm4uO3bs4Otf/zqNjY1KgtEi1HEthGDeBFrAXcDJGapQtiRJzUKIAOCAEKJCkqRjs244B/mWLARVVVXK746g0WgUZ5mHhweenp709PTwQWugNTQ0MD4+TmxsrMJPFB0dTWNjI+Xl5eTm5lJRUUFqairBwcFO1STZGuTn50dycjJCCCwWC+fPn2f//v2sXbsWf3+HJbcU+Pv7c/LkSe655x6efvppurq6ePnllxcFYQ5cCyFwRqzlCPcyQxWSJKn58s92IcSbTKpXs4RgLshJNZ/4xCcApjnLxsfHGRkZYXR0lKGhISWZpqOjQ4ntqa6upq+vT0mi8fPzm5dX1mKx4OHhMW2AajQaIiIiCA8Pp6GhgdLSUk6cOIGfnx+ZmZmz9h6SJJGfn4/dbp9mDYqNjcXPz49Tp05x9OhRUlJSSExMVN1vmEwm9uzZw2c+8xn+8Ic/MDIywmuvvYZWq73Sr/SWwbUQgjzmQaAlhPAENgKfnXLMDGgkSRq4/Pt24IfXoE9KqIRaiILdbqeyspLi4mIiIiIYHBzk0qVLXLx4ESEEfn5+BAcHExIS4tBKI1M8ZmZmOhyYQggiIiIICwujurqasrIyDhw4QHR0NKmpqUrf6urqaGlpISMjY9a+xMvLi5ycHPLz8ykpKaG3t5eVK1eqJvVotVpee+01XFxc+MMf/kBfXx+7d+9eDLVwgqsWAkmSrEKIrwD7mDSRviRJUpkQ4kuXzz9/uenHgf2SJA1NuTwQePPyANIBf5Ikae/V9mm+0Gg0jI6OotVqlRnYarXS1dVFW1sbra2tFBcXU1xcjIeHB+Hh4URERCgDtaqqCp1ON6cqpdFoiI2NJSIigrKyMqqqqmhsbCQzMxM/Pz9FDXJmDdLr9axZswYvLy9KSkoYHBwkOztbNeVTo9Hw9NNPMzw8zN/+9jc+9alP8cYbbyyqRg5wy6dXHjlyBLvdztatWx2eHxoaorm5mYaGBjo7J624fn5+hIeHc/78eaKiolix4spcG319feTn59PV1YXRaGRiYoIdO3bMyzrV3NzMmTNn0Ov1rF+/Hi8vL4ftent7OXr0KBqNhrfeeouXX36Z++67j1dfffWWFIRF8i0nkEOvnQ0kmDRdxsXFsWXLFu68807S0tIYGxujsLAQu92OzWZjaGjI6fWO4OnpyebNm4mIiGBsbAyYVK3mMyGFhIQoAnv48GHa2tpmtZEFQKfTsWXLFn77299y33338dprr/GlL33pivp6K+CWFoKBgQGsVuu8w6ddXV1JTExkx44dmEwmDAYD9fX17Nmzh9OnT9Pb2zvvZ4+OjtLS0oK3tzc+Pj7k5uZy5syZeeUdeHp6snXrVsxmM8ePH59G6jVVADZt2oSbmxsajYZXX32Vu+++mxdffJEf/OAH8+7nrYBbOnZI9hQ7CpdQQ3t7OyMjI6xevRp/f38uXbqExWKhoaGBkJAQUlJSVFeXqdagNWvWYDabuXjxIqWlpXR3d5OVlTWnYLq6urJ582ZOnDjBmTNnsNlseHl5zRIAGRqNhr/85S9s3ryZH/zgB4SFhfGFL3zhit77ZsVNIwSyarJr1y6EEOh0umllm2TWOTmJxtXVlZ6eHjQazRUzzVVVVWE0GgkLC0Or1ZKenk5iYiJVVVVUVlayf/9+IiMjSUlJcVjco6amhtbW1mnWoKSkJPz9/Tl9+jSHDh0iMzNzzhxlg8HAhg0bOHHiBGfPnlXyImYKwNT2e/bsYc2aNTz00EOEh4ezffv2K3r3mxE3jRAIIdBoNAQFBSkCMTExwdjYGH19fYyOjk7TufV6PZIkodfraWxsxNfXd17VaIaGhmhpaSEhIWGa7d1oNLJ06VJiY2OpqKhQLEAJCQkkJiYqJs3h4WHOnz+Pv7//LGuQn58f27dvJzc3l4KCAvr6+sjIyFDdyMq1Eg4fPozNZiMmJkaVPc/T05P9+/ezcuVK7r33Xs6ePUtsbOyc730z46YSAtnU6Qh2u52RkRGGhoaURJrq6mpsNhtnzpwBJjfBAQEBBAUFERQU5NCuXl1dDeB0ljYajaSnpxMXF0dxcTHl5eXU1taybNkygoODpznFHPkWjEYj69ato7i4mMrKSvr7+8nKynKaL9Db28vx48dxcXHBbDZTUlKC2Wx2GFoiIzw8nDfffFPZ7Ofl5d3SlXduGiGYCxqNBrPZrAz0gYEBLBYLy5cvx8fHh87OTtrb22lsbKSmpgaNRkNAQADh4eGEhoZiMBiw2WxUV1cTHBw856rh6urKmjVriImJoaCggBMnTuDt7U1PTw+ZmZmqs7VGoyEjIwMvLy/y8/M5fPgwGzZsmOUXmLkJdnFx4fjx4+Tm5qLT6VTDuteuXcuvf/1rvvCFL/DpT3+a3bt335KmU7iFrUNTN8Xe3t7ExcWRnZ3N3XffzebNm4mNjWVgYIC8vDzeeecdTp8+TXl5OWNjY1ekPvj7+yv8RD09PQghMJlM87o2KiqK9evXMzIywqFDh+jr61POObIC6XQ6JRL19OnTdHV1qd7/n/7pn3jooYd47733+M///M95v9PNhltaCDQaDZ6entOOazQa/P39ycjI4Pbbb2fr1q1ER0fT1tbGhQsXEEIwMDDgkOHaGTQaDf39/Wg0Gtzd3Tl16hR5eXkKY50aAgMD2bx5M5IkceTIEXp6ehwKgAzZiebi4sKJEyfm5Ct65plnWLFiBT/4wQ84duyKQ7ZuCtzSQuDp6amqAggh8PX1ZdmyZQqFuslkorCwkN27d1NaWqo4u9QgW4PS09PJyckhMTGRmpoaDh06RH9//5zXe3l5sWXLFnQ6HUeOHOHw4cMOBUCGi4sLGzZsAOD48eOqvge9Xs/rr7+Ou7s79913H93d3U7b3qy4JYVA9hRfCcdQbW0tGo2GnJwctm7dSkBAAOXl5bz77ruUlpY6XRmGhoYoKirC39+f2NhYtFotaWlpbNiwgdHRUQ4ePKgwXajBzc2N5cuXY7VasVqtZGRkqO4r3N3dyc7OZnh4mNOnTzsNLYfJ/IyXX36ZlpYW/umf/mnuL+Mmw00jBHIOQXt7O52dnfT19TE8POyQBnFwcJCJiYl5C8HExAR1dXVEREQoDBXZ2dls376doKAgysvL2bNnD1VVVdMGm+wUA2ZZg4KCgsjJycHDw4NTp05RVlamGjbR29tLbm6u4vPIz8+f00Pt5+fHsmXLaGtr4/z586ptP/axj/H5z3+et99+m5dffnk+X8tNg5vGOiQz0B09enTWOYPBMI11Tp615ysEMuP0TLOol5cXWVlZdHd3c/78ec6dO0d1dTWZmZn4+/tTXV1NW1sby5Ytczhry17f/Px8ysrKGBgYYOXKlbNi/2fuAWAy8O/9999n69atqitCdHQ0fX19XLp0SQn8c4b/+q//4ujRozzyyCOsX7/+lvEf3DRRpOnp6UiSxIEDB5RSTePj44yNjc2qTCO/s16vV5JoAgMD8fLymmW7lySJffv2odVqycnJcdonSZJobGzk/PnzDA8PExERQVNTk0KvopYII0kSFRUVlJSUEBAQMM0v0NPTw/vvvz9rD9Df38/hw4cxGAxs2bJFNf3Sbrdz5MgR+vr6yMnJcegTkCSJuro63nrrLb75zW+yfv16Dh8+7PSeNxrUokhvGiGYbyi1zWbj8OHDjI2NERQURGdnp7I5NZlMBAcHEx4ejr+/PxqNhvb2do4ePcrKlStZsmTJnH2zWq2UlpZSWVkJwKpVq+adullbW0teXh6enp5s3LiR4eFhhwIgo6uri6NHj+Lp6cmmTZtUE22Gh4fZv38/JpOJbdu2TVttRkZGKCgooLm5GT8/P/72t7/x3HPP8fvf/57Pfe5z8+r7QseiEEyBJEm8/fbbhIWFKXkAIyMjCsFWa2srVqsVFxcXIiMj6evro7u7mzvvvHPOEk0yLBYLBQUFGI1GxsbGiI+PJzU1dV4pji0tLZw6dQqTycTY2Bh6vd6pFQigqamJkydPEh4ezpo1a1RXnJaWFo4fP058fDwZGRlIkkR9fT2FhYXYbDZSUlKIi4tjbGyMpKQkhoaGuHjx4hUHGC5EfOj5BEKInUKIi0KIKiHEtxyc3ySE6BNCFF3+fH++115rDA0NMT4+Pm0/YDKZiIqKIisri4997GOsXbsWX19fKisraW1tRafT0d7ePq94/6GhIc6fP09AQAC33347sbGxVFZWcuTIkXnlHQQHB5ORkaEQgq1du1ZV5w8NDSU1NZWGhgYuXLgw571jYmKorKykoaGBkydPkpubi7u7O9u3bychIQGNRoPJZOKFF16gq6uLf/u3f5uzzzc6rloIppBv3QYkA/cJIRxVnzsuSVLG5c8Pr/Daa4a5iLZ0Oh3h4eFkZ2crAW42m40TJ06wb98+6uvrnZobZ1qD9Ho9y5YtIysri4GBAQ4cOEB7e/uc/SspKcFoNCKEID8/f05fRGJiIhEREZSWlircqc6QlpaGi4sLp0+fVnwXmzdvnrVP2LFjB3fffTevvfbaTV8t56Mg37pW134gOPMUz4TdbqehoYHAwEDuuusuVq+epFI6c+YMBw4coLm5edbKIFuD0tLSpsUWhYWFsW3bNlxcXHj//fepqqpy+Mze3l5lD7B161bWrVvHwMAAx48fV/UuCyFYsWIFnp6e5ObmMjw87LDdyMgIubm5jI6OApOriDz7O8Ivf/lLjEYjDz/8sKqf4UbHtRACR+RbjvjR1wohzgsh3hNCLL3CaxFCPCiEyBdC5Hd0dHzgzsqe4rn08+bmZkZGRhTq9MjISHbs2MHatWuVleH9999X4nmmqkGOIkzd3d3ZsmULQUFBnDt3juLi4mlC5CgUIigoiDVr1tDd3c2ZM2dUB6JOpyMrKwu73T7LOSbr/vv27aO1tZW0tDTi4uJoaGhA7bsMDw/na1/7Gnl5efzhD39Q/b5uZFwLIZgP+dY5IFKSpHTgV8BbV3Dt5EFJ+o0kSSskSVrhiIRKkiSsViv5+fmcO3dOsdA0NDTQ1dXF2NgYdrt93p5ii8Wi0LErLyoE4eHh7Ny5k8zMTHp7e9m/fz/FxcWcPXsWmO0UmwqDwUB2djbR0dFUVFSQl5enMGI7iwUKCwsjMzOT5ubmOTf97u7urFixgq6uLmV/MDo6yqlTpzhz5oyi+ycmJpKamorZbCY/P1+1rtr3vvc9QkNDefzxx68oXupGwnUh35IkqX/K73uEEP8thPCbz7XzhVyzrLm5WUmomQmDwcD4+DjDw8OKOdBRnP7AwABtbW2kpKQ4VBU0Gg1xcXFERERw/vx5KioqgEndfK4Qa41Gw/LlyzGZTJSVlTEyMkJ3d7eqFSguLo6hoSEqKyvx9PRUzTiLiIigpaWF8vJJPuRLly5htVqVQh/y++h0OpYtW8bx48e5dOmS06o3RqOR7373uzz00EM888wzfPOb31R9vxsRV20iFULogEpgK5PkW3nAZyRJKpvSJghokyRJEkKsAv4GRDLJU6R6rSPMx0QqV6ORHWUDAwM0NzfT0dGBEAJJkhBC4OPjoxTlljeHhYWFWCwW7rzzzjk5QAcHB9m7d6/yzOTkZJKSkuYVm19UVERlZaXiiFNL87Tb7Zw4cYK2tjY2bdqkSsk4MDDA/v37lRpsq1evdnrvEydO0N7ezm233eY0xNtut5OUlER3dze1tbXzysBbaPhQTaSSJFmZJNjdB1wA/iKTb8kEXMAngVIhxHngl8C90iQcXnu1fYL/q0bj5eVFSEgICQkJ+Pr6IoTgYx/7GJs2bSIxMRG73U5JSQnvvfce+/fvp6KigpqaGkJDQ+cUANkapNFo2LZtm0KudezYMWXz6Qw9PT3U1tai1+ux2WyUl5er6vwajUZJyj9z5ozD+8u6/6FDh5R7qVXOgcnJw263U1xcrPrsJ554gs7OTn7+85+rvteNiFvKWfb+++8zNjY2K7l8eHiYxsZG6urqFBNqcHAwaWlpqlakqqoqzp07x/Lly4mJiUGSJGpqaigsLMRoNJKdne1w/zEzFKKhoYGSkhJiYmJYtmyZqsOrp6eHQ4cO4e/vz4YNG5S2o6OjnDt3jsbGRnx8fFi1ahUVFRXU1dWRk5Ojyn5RXFxMRUUF27dvd9rObreTmppKR0cHDQ0N867AuVCwSL6Fevi0q6sr8fHxbNu2DXd3dwwGA+3t7ezbt4+TJ086jNYcHBykuLiYwMBAoqOjgcnVJzo6mi1btgCT5FhNTU3TrptqBpX3AElJSSQkJGCxWJT9hTN4e3uTmZlJW1ubEprR0NDAvn37aG5uJjU1lS1btuDh4UF6ejoGg4H8/HxVR19iYiJ6vX7O1eDb3/42HR0dPPvss6p9vNFwywjB8PDwLE/xTHR3dzMwMEBKSgp33nknycnJtLe3s3///mn2d0mSFAfSihUrZs3c3t7ebNu2DU9PT06dOoXFYgEcm0FlpKWlERERQUlJyZz5BdHR0YSEhFBSUsL777/P6dOnMZvN5OTkTNuPGI1G0tLS6O7upq6uzun9DAYDSUlJtLa2qjrzPvOZzxAbG8uzzz6ralG60XDLCMF8SjLJBLuRkZEYjUZSUlK44447SEhIoKGhgb1791JRUcGlS5fo6OggPT3d6SbRxcWFTZs2ERQUREFBAUVFRU4FACZXkZUrV+Lj48PZs2dVcwXkemd2u12xYm3ZssWh6hYVFYWPjw/FxcWqJs7Y2FhcXFwUq5IjyNUzm5ubefXVV522u9FwSwmBXAvMEcbGxmhoaCAyMnIa1YrBYCA9PZ2dO3cSEGOHIf8AACAASURBVBBAcXExRUVF+Pj4KGqQM+h0OrKzs5XifJIkqQbDycXI9Xo9p06dcpgWKdv98/PzFQGUazM4ghCCzMxMRkdHuXjxompfExISlKQkZ/jiF7+In58fv/zlL9Ve/YbCTSMEclGON954g3feeYf33nuPI0eOkJubS3l5OS0tLbi7uzv1FNfU1GC3250mkri5uZGVlaVYWnp7e7l06dKcQXX9/f1KTbSJiQmam9XdICaTiTVr1jA0NDRLl29sbFR0/5SUFHbu3EloaChlZWUMDg46vaevry/h4eFUVlaqWq1iYmIwGo2qq4HRaOSzn/0s586d48CBA6rvcqPgphECmYFO1pc9PT2RJIn29nZKS0vp7e2lv7+fXbt2kZubS319vTLTSpKExWLB399f1RpksVjo7+8nPT2doKAgioqKOHXqlFM1Q94D6PV6cnJyCA0NpaioSFU/h0maltTUVIUDaWxsjNOnT3Pq1ClcXV3JyckhOTkZrVarFAiZy5uckpKimGKdQafTERcXR2tr6yx1TJIkOjo6OH36NGlpaej1en71q1+pPvNGwU2TXikz0Dkqpt3X18e+ffsICwsDJuPq6+rqFNpGT09PhoaGSE1NdXp/2RoUFBREfHw88fHxSpWbQ4cOsW7dumlqjqNN8Jo1azh27Bh5eXmYTCYCAgKcPi8hIYHW1lYlzshqtSrlmqaqPq6uriQlJVFSUkJraytBQUEO7+fu7s6SJUuorq4mKSnJqWMsJiaGCxcucOnSJVauXKnkV1ssFvr6+tDr9WRkZLBlyxaFLeNKuVwXGm6alUANMvdOQkKCkjOwZcsWYmNj6enpUeJsBgYGHOrhsjVIo9Eo1iAhBAkJCWzcuJHR0VEOHTqkkF05swLJOr+bmxunTp1SVWHGx8fRarXY7XaliEhycrJD3T8+Ph43NzeKiopUHW6JiYlIkqS6NzAajURFRVFXV8fZs2fZtWsX586dUyJV77rrLjIzM3nooYcYHh7mxRdfdHqvGwW3hBDIm2JZ1dFoNPj5+ZGRkcHGjRuBSR7SsrIyhU9oqooz1Ro0kwoxICCArVu3otPpOHr0KFVVVapWIDmIDuDkyZMOQ6QbGxvZu3cvbW1thISEYLVaFeuWI2i1WlJTU+nv71dVtdzc3IiMjMRisTjcG9hsNhoaGuju7sZut1NXV0doaChbt24lJyeH6OhoJbvuzjvvJDQ0lFdeecXp824U3DJC4OHh4TA9sra2FiEEmzdvnkWhUlNTw8DAACUlJQQFBTnNMXZ3d2fr1q2YTCZl1lSzArm7u7N69Wr6+vooLCxUjo+NjXHmzBklvXLbtm1kZ2fj7+/P+fPnVTe1YWFh+Pj4UFZWpmrDT0xMVDhVZQwPD1NaWsq7777L6dOnGR8fx2w2YzKZWLVqlRJuMhUajYZPfvKTlJSUzOngW+i46YVAzVNss9moqakhJCQEV1dXhUJFpjHJy8vj4MGDiiqgFs4wMjLC2NgYQggmJiYYGRlR7VdwcDBJSUnU1NTQ0NBAU1MT+/bto7GxkaVLl7Jt2zaF/UIm3SotLXV6PyEEKSkpDA8PU1tb67Sdh4cHQUFBVFVV0dLSwsmTJ3n33XcpLy/H29ub9evXc/vttyv3UnOeyUU+fvvb36q+60LHTS8E8uB0JASNjY2MjY3NCk329fVly5YthIeHMzExgd1up7W11ak5VI4F0uv1bNmyBbPZzIkTJ6YR6DrC0qVL8fLyIjc3l5MnT+Li4sK2bdtYunTpNN3fw8ODuLg4qqurVZ1ogYGB+Pj4UFFR4XRvMD4+jqurK6Ojoxw/fpzOzk4SEhK4/fbbWb9+vVJwPCwsDIPBQE1NjdPnpaamkpiYyNtvv636ngsdN40QjI2N0dPTwyuvvMI777xDfn4+/f39qp7iqqoq3NzcCAwMnHVucHCQ5uZmAgIC8PX1JT8/n7y8vFmqxsxgOF9fXzZs2IBOp+P48eOqKkxLSwvDw8PY7XZFpXLmzEtOTsZgMKjG9wghSE5OZmhoaFodM7mfeXl57Nq1i+rqaoUcWC5GOFN102q1hIeH09TUpOpp/vjHP05VVdU0te5Gw00jBGazmfr6ev7xH/+Ru+++m5UrV+Lp6UlGRgY/+9nPeOKJJzhw4IAyiHt6eujq6iI2NtYh4ZZsDVq1ahUbN24kOTmZ2tpajh49qgxsR8Fwcl/WrVvH2NgYp06dmiU4su5/8uRJXF1diYuLY2BgQDVJ3mAwkJiYSGtrq2pKZHBwMB4eHlRWVmK1WqmtreXgwYMcOHCA+vp6IiMjycnJYenSpQwMDKgyYERERGCz2VQdfJ/97GRt9j//+c9O2yx03DSh1HJUpbwiNDc3U11dTWlpKRUVFco/Uq4QL4cX33333bOyyyorKykqKppFuNXY2Ehubi4mk4nMzEylGIazTXB9fT1nzpxReH5gMndZZpBISkoiOXmSXOPgwYOMj4+zc+dOp/xGVquVPXv2KORczlBeXk5paSk6nQ6r1Yq7uzsxMTFERUUp7zoyMsLu3buJj48nPT3d4X0kSWL37t14eXkprNyOsGTJEjw9Pa+4hvT1hFoo9TVxlgkhdgLPMpkp9ltJkp6ccf5+4NHLfw4CX5Yk6fzlc7XAAGADrM46OhcCAwMdqjUy6uvreeutt3j77bd5++23+etf/0psbCwdHR388z//sxIvJFuDgoODZzHHhYWF4eLiwrFjx5QSSWpWoIiICDo7O6msrMTHx0dx0nl6erJ+/fppKlpmZiZHjhyhoqKClJQUh/eT43vOnz9PV1cXvr6+yjl532KxWJQVRa/Xk52dTUBAwKzVzmQyERQURH19PWlpaQ43/XJOdVVVFRMTEw7LVwHk5OTwu9/9jubmZtXqOAsV14t3qAbYKElSGvAj4Dczzm++zEf0gQRgPoiIiOCrX/0qhw4doq6ujq9//euMjIzw8MMPExcXxx//+EfsdruiBi1fvtzhwJg6S8uhGmqQfQtnzpyhrq6O5ORktm3bNmuP4u/vP6/4nujoaAwGg+LwGh0d5cKFC+zZs4cTJ07Q09NDcnIyUVFRjI2NOeRXnfqdjIyMqAbMTY1WdYZPfOIT2O123nzzTbWvYsHiuvAOSZJ0SpIk2dtzhsmE+o8MQUFB/PznP6empob/+q//YmJigs9+9rOsW7eOyspKMjMzZznF4P88wQaDgaysLCYmJjh+/LjTjeP4+Dj5+flKHoJc49hZEN984nv0ej3R0dE0NjZy4sQJdu/erRTrW7t2LXfeeScpKSnEx8djt9upr693eq+QkBC0Wq1qG5mMQG1fsHnzZlxdXTl48KDTNgsZ15N3SMYXgPem/C0B+4UQBUKIB51ddK14h6ZCr9fz8MMPU1VVxb//+79TUFDAN77xDfbt2zer7cxQiLCwMLKysujv73fICdTc3MzevXupr68nOTmZhIQEmpubVe3u7u7uREVFUV1d7dDPYLVasVgsyoBsbW0lOjqanTt3snnzZsLDw5WVycvLCy8vL9UBrtfrCQoKckgkJkOOr2ppaXHaxmAwsGzZMqUK6I2G68U7NNlQiM1MCsGjUw5nS5K0jEl16mEhxAZH187FO3Q1MJlM/OIXv+DIkSOEhobypS99iQceeECZ4Z3FAgUFBZGZmUlLS4sSfzQ+Ps7Zs2c5ceIERqORbdu2kZKSwtKlS3F1daWwsFA1vicpKQlJkpTUSZgMxy4sLGTXrl0UFBQoISByMJuzALbw8HC6urpUY5RCQ0MZGRlR9T8EBQUxNjam2mbDhg20traqrmILFddCCObFHSSESAN+C9wtSZJSVlGSpObLP9uBN5lUrz4SZGVlUVRUxN///d/zP//zP2zZsoXGxkbVWKCYmBgiIyMpKyvjwoUL7Nu3j7q6OpKSkqbp/jqdjvT0dPr6+uaM7wkLC6O6upq6ujqOHj3K3r17sVgsBAcHs2XLFrZv305qaipjY2Oqaopcy1gtXTMoKAghxKxc6KmQo13VVrGdO3cC3JA5BtdCCPKAOCHEEiGEAbgXeGdqAyFEBPAG8A+SJFVOOW4WQrjLvwPbAeexAdcBJpOJ//3f/+Wxxx7j5MmT7Ny5E6vV6tQKJIQgNTUVnU5HSUkJer2erVu3OqRin098z8jICDqdjomJCXJzcxkcHCQ1NZU777yTNWvW4OfnhxCCoKAgjEajqrpjNpvx9vZWFRQXFxe8vb1VN76urq64u7urtlm9erVC9Huj4apNpJIkWYUQMneQFnhJ5h26fP554PuAL/Dfly0Vsik0EHjz8jEd8CdJkvZebZ+uBX70ox9hNBr5/ve/z09/+lPuuOMOh+1aWlrIz89XokEDAgKc8vnL8T3Hjh2jtrZWCdeQE1aqqqpoampCkiR0Oh0Gg4HbbrvN4UZao9EQHh5OTU2NqvkyJCSEsrIyRkdHnfIoBQYGUlFRoXoff39/GhoaFNKymZCdeTei5/iaeIwlSdojSVK8JEkxkiT9+PKx5y8LAJIkfVGSJO8p1OwrLh+vliQp/fJnqXztQsFjjz3GT37yE/Lz8/nUpz41TZeXdf/jx4+j1+vZtm0b0dHRSvKJM8jxPRcvXmRsbIxLly6xd+9ejh49Snt7O/Hx8dx+++2kpaUxPDyseq/w8HBsNpuqp1nmUlWbxQMCAhRBdAZfX18mJiam1UW22Wx0dXVRVVVFXl4eMTExWCyWOQsKLjTcNJllHxa+9a1vUV9fz3PPPcf3v/99nnjiCWX2Hx0dVby+Wq0Ws9ms1C2T6wjPhBCC0NBQSkpK2LVrF3a7HR8fH1auXEl4eLjih4iIiKCoqIj6+nqnK4uvry9Go5Hm5mZF/58Jb29vDAYDbW1tREZGOmzj5+eHRqOhs7PTqbNLjmm6cOECGo2Gnp4e+vr6FIuRvBLYbDby8/PZtm2b8y91gWFRCOaBX/3qVxQXF/Pkk0+yZMkS3N3d8fDwIDs7e9oANRqNJCYmUlxcTEdHxzS+UJvNRmNjI1VVVUoGmsxS52iQGwwGgoKCaGhoID093aEKotFoCA4OpqmpCbvd7tBxJ4QgICBAqbTj6D5arRYvLy+lX1arld7eXnp6epSPvCLV1dVhNBrx9vYmODgYb29vvL29cXV1JSIigh//+McUFhYuCsHNBq1Wy1//+leWLl3K448/zptvvklmZqZDXT02NpaLFy9SXl7Oxo0bGRoawmKxKAnzbm5upKenMzw8TFVVlSrfaXh4OM3NzXR3d08LkZiKoKAgamtr6e3tdbpi+Pv709jYyMjIyCwnoDzgtVotnZ2d7N27l4GBAWWGlwd8SEgITU1NSkFzR8KUmpqKwWCYs2byQsOiEMwTwcHB/OQnP+HLX/4yr776qlL0byZkxobS0lKOHDmihCQEBwcTGxtLYGAgQggGBga4dOkSDQ0NJCQkOLyXbL5sbm52KgSy+bKtrU1VbYJJE6fZbJ42w08d8PB/bHTyDG8ymZQBPzExoZqwo9friYqKuuEyzRaF4ArwpS99iT/96U+88MILfO1rX5ulY4+NjVFTU6OkLnZ3d5OYmEhMTMysGdjd3R1vb29VITAajfj6+tLa2uqUCcPFxQV3d/dZ8T8TExOKStPd3Q2gFBKRr/P29iYsLAxvb2+0Wi3Hjh0jNjbW6f7C09MTq9XK0NCQ06DBqKioxZXgZsezzz7LypUrefTRR/nzn/+MJEl0d3dTVVVFQ0MDdrsdPz8/XF1d6e3tJSkpyWlodGhoKKWlpYyMjDilQAkICKC8vJzx8XGHBUVgcvPb0tJCRUWFMvCnWnFcXFzQ6/UYDAYyMzOVGX4qbDYbQgh6e3udCoHsmR4YGHAqBNHR0Rw4cED1nRYaFoXgCpGZmcltt93GG2+8wenTpxkcHKSnpwedTseSJUuIjY3F09OTtrY23n//fVpaWggPD3d4r5CQEEpLS2ltbXWaxC8LQVdXF8HBwYyPj8/atMoDvri4GJPJhLe3NxEREdNUmry8PJqampT0yZnQarW4u7srhc0dQR74amEYsbGxSJLEhQsXWLZsmdN2CwmLQvAB8J3vfIfdu3fz7LPP8o//+I8sW7ZsFoepv78/Li4uNDY2OhUCT09PjEYj7e3tDoVgfHxcccKdP3+ewsLCaQNQHvB+fn7U1NTMSgKaCi8vL2pqahgdHXU6Q7u5uakOcBcXF7RarWqb+Ph4ACoqKhaF4GbG2rVrSUtL4+jRo/zxj3906tENDg6msbFR1Xzp7+9PR0cH4+Pj02b3np6eaYNtcHBQSfSRZ3jZsmS1WqmpqXFauhVQSlENDg6qCkFbW5tTU6oQAjc3N9WUTFkI56KXX0hYFIIPiAceeIBvfOMb7N+/n9tuu81hm6CgIGpqauju7sbPz085LkdkyvUQhoeHeeutt5TzZrMZLy8vhVa9urqarq4uhbRrJnQ6HWazedo+YCZkBuvBwUGn9c7c3Nyw2Wyqq4XJZFKlk5GFYC7i4YWERSH4gHjggQd49NFH+ctf/uJUCOTBVlNTQ0dHhzLDT51J5dk8KipK0eNnlkLq6emhsbFRdXM8lypjNpsRQqi2kQe+2qbWxcVFNZTDbDbj5uamGsqx0LAoBB8QMnP08ePHlWNykn93d7cy4AGFu0eO6oyOjsbHx0dJfXzrrbcUUixHkFWZoaEhp0Igh2w4g0ajwWg0qqZuygI5V5vR0VGnKhNMWquuVeLT9cCiEFwF1qxZwwsvvMB7772HzWabppO7ubnh6+uLXq9nbGyMnTt3Oh3ARqNRVZWZapVxVmnHZDIxPj6OzWZzmr45lyozdSVwBqPRiCRJTExMqK5Kau+z0LAoBFeBzMxMJEmioKCArVu3Ehsbq8zw8gApLS3lwoULTgcmMOdmcz6DU24zOjqqWkJKbZaX+6xGtiVbwNSEYC5T60LDohBcBeQqlWNjY6xdu9ZhGzc3NyRJYnh4WFFrZsJkMqnq2QaDAY1GozqA5X3E2NiYUyHQ6/WqM7RWq1W4VNX6AuqC4u7ufkNtjG8aBrqPAtHR0bi4uCjVKR1hPrP4XDO0EAK9Xu+wdoKMqTO0Whu18/N5juz9VruPq6vrnITECwnXRAiEEDuFEBeFEFVCiG85OC+EEL+8fL5YCLFsvtcuZGg0Gnx9fVUTVqbO0M6g1+uxWq2q9c/mGsCyEDiqdyBDTttUg1wYRO08oNpGfp8bBdeLfOs2IO7y50HguSu4dkHDw8NDVZWZzwyt0+mQJEm1roBWq1U9Lzvj1AanRqOZs9CgRqNRvYdsEZpLCG6kOsfXYk+gkG8BCCFk8q2p3Bt3A3+QJv8DZ4QQXkKIYCBqHtcuaAwODtLa2jqL3n0q5EHhbOBoNBq0Wi3f/va3nd5jLkERQqDT6eb1nG99y/mCO9/nqK1cal7phYjrRb7lrM28ibs+DPKta4GhoaE59V95w+kMctFBNajVKl5oz5mYmFCy1G4EXIuVYD7kW87azJu4S5Kk33CZw3TFihULhko7MTGR9vZ2Ll265PD8xMQEb775JmlpaSQmJjpsY7FYKCgo4K677nI6gx45cgSYpDx0hOHhYXbv3s3y5cudrkoVFRUUFxfzd3/3d07Du/fv34+rqyvr1q1zeL6/v5+9e/eyZs0apyHXX/7yl3n++edVfRYLCdeLfMtZm3kRdy1kyKqBM8gzppouPh89Wwgxp74/1z3kNnPtPebaXIP6HkdOylcj61pIuC7kW5f/fuCylWgN0CdJUss8r13QmMsSMp/BOR+Li6zvq50HdevQfJ1h8/ETqFm7QkMnNdq5ipYvFFy1EEiSZAVk8q0LwF9k8i2ZgAvYA1QDVcCLwENq115tn64n5rLxy5vR+Xph1dqonddqtWg0mnk9R80PMNf76HQ69Hq96j5ITheV+VkXOq6Jx1iSpD1MDvSpx56f8rsEPDzfa28keHh4qIY8wPxt/Ffr6DIYDKoz9JUEyDnLgYDJYL2Z7zw1403GyZMn+fznP+/0WQsFi2ETVwlPT09GRkZUKQznKwRqM7TRaJwzQG6uWXw+3muz2YwkSYyMjDgNvzCbzXR1dVFRUeEwAcjV1ZXg4OBbayW4leHp6QlMxvzL9CczYTAY5lRBYO4ZGiZ1cUcFRGDuKFEXFxc0Go3qyiXHNw0MDGA2m1Uz3oqLi3F1dcXb23tWxtvy5cs5evSo6uSwULAoBFcJeeA3NjY6FQKTyaQaVWk0GhFCzDvM2ZkQmM1mOjs7ncb6azQazGaz08SasbExRZ06f/48BQUF0wRGznjz8/OjtraWtWvXOs2f3r59O7t37+bw4cPs2LHD6XstBCwKwVVCtpXX19c7TSx3cXFRjS8SQuDi4jKnmgKTzjlnRFxubm5MTEwwPj4+KztNhhzmIScATf1MHfByTnN0dLQyw8v3nJiYoK6ujt7eXqdC8OlPf5qvfe1r/OlPf1oUgpsd8iCYWTx7KkwmExMTE1itVqdOKkebzZnnAdVQaFmV6e/vn5ZHPDXjrb+/n4GBgWlV6GdmvFVWVtLf309WVpbD5+j1enx8fFRJwQICAli9ejW7du1SpYVfCFgUgquEXOZVLbVx6iwu7yFmws3NjdbWVqf30Ol0mEwm1Rxh+d719fXTcpqnZrzJgzEmJkZhn5uZHNPb20tLSwvDw8NOVa/g4GBKS0sdtunr68NisZCdnc2pU6d4/vnneeSRR5z2+6PGohBcJfz9/XF1dVXl6JzK9KAmBKOjo6qrhaen57SI1dHR0Wk5zTLdopzf4O7ujq+v77SMN5vNxq5du3Bzc3Na93lqeaaZtZxlhIWFUVpaSmNjo1Ips6mpCYvFQnt7OxqNho997GO88sorvPjii3z1q1+dMybpo8KiEFwlNBoNISEhqkIwH+Y2meKwv7/fIbHuyMgIWq2Wvr4+jh8/Tm9v77Q9hLu7OwEBAXR3dyNJEtu3b3dqlXFzc6Ojo8MpB6qXlxcuLi40Nzc7FQIPDw+F0Gt8fJzq6mpGR0dxdXUlNTWVJUuW4OLiwkMPPcTjjz/OX/7yF+69916n7/9RYlEIrgEiIiJUmZiNRiMGg0FVn5dXiN7eXkwm0yzWiqnm076+Pvz9/fH29lZmeHnAyznNaggICFB4U52RgoWEhFBfX+9wZZIkSal30NfXR19fH8HBwcTExBAUFDTtnt/85jd57rnneOyxx/jEJz6xIM2li0JwDbBkyRLef/99Vb6emaoMoDil5AEvhKCwsJD8/HyljYeHB4GBgYr9/cyZMyQlJTmNFPXz80OSJLq6upxSuAQGBiqEXs6IuOQKmq2trYSFTdZeHx8fp7a2FovFwsDAAHq9XhEYZ8Rgrq6ufP/73+fhhx/mscce46c//anDdh8lFoXgGiAtLQ2bzca5c+ecDgZPT09qampobGycRqgrz/BCCDQaDTqdjtTUVHx8fJRaxTIkSaKwsJCuri5VIdBoNLS1tTkVAnm2bm5udioEAQEBuLi4UFtbi9lsxmKxUFdXh81mw8fHh1WrVhEeHk5xcTFVVVUMDQ059TB/8pOf5MUXX+TZZ5/l7//+71m+fLnT7/KjwKIQXAOsXLkSgIKCArKzsxV2iak2+M7OTmw2G6dOnUIIoZBtyTZ4Ly8vSktLqaqqIiYmxmFohBACX19fVWIrnU43Z96zXq9XqtekpaU5dKxJkqSUf21ubkar1RIREUFsbOw07qOEhAQsFgvl5eXK9yDDarVSVlZGZWUlX/nKV/jqV7/K/fffT2Fh4YLKPFsUgmuAZcuWodVqOXXqFBkZGfT09CieV3nABwQE0NzcTHJyMomJiQ4tQL6+vlRWVtLb26tamaa5uVl15g0ODqa4uFi1TUREBHl5ebNKQc0sLwWT6tPatWsd8gy5uroSExNDVVUV8fHxyt6mq6uLs2fPMjAwQHR0NOnp6dhsNv7lX/6Fz3zmM7z++usLxlq0KATXAEajkbCwMC5cuMDIyAghISF4eXkpKo1Op8Nut/PWW28xPj7u1AQqD8bOzk6nQiCbNdva2oiOjnbYJiQkhOLiYpqbm4mLi3PYJjQ0lIKCAurq6pRi3lVVVbS0tCh6fmxsLJWVlXR3d6sO2OTkZGprayksLGTdunXK7G8ymdiwYYOilj344IOcPXuW3/3ud3zjG9/gqaeecnrP64lFIbhGyMjI4MCBA2zbts0pVbuPj49q7u3UyvHOzJceHh6YTCZaWlqcCoGHhwceHh40NjY6FQKDwUBISAg1NTU0NzczPDyMi4uLsumWHWAajYYjR45QXV2t1B6YCaPRSEpKCoWFhezZs4fR0VFl9p9pDXrhhReor6/n6aefRpIknn76aaffx/XCwliPbgJkZWUxPDzMuXPnnLbx9fWlt7dXNaw6ICBA2T84gjxLt7a2qmaahYeH09HRMatmgWw5ys3NpampCZvNhkajYe3atdxxxx2kpqZO8wD7+/vj7+9PRUWF06y1qTyso6OjrF69mhUrVjg0h2o0Gv7jP/6DFStW8Mwzz/DII4+oZtRdD1yVEAghfIQQB4QQly7/nMUWK4QIF0IcEUJcEEKUCSH+bcq5x4UQTUKIosuf26+mPx8lZErGw4cPO20TGBg4Z+X4oKAgrFbrrEJ8UxESEoLNZlPd/MqBfXKKo9Vqpbq6moMHD3Lo0CGamppYsmQJbm5uaLVawsLCnOYppKamMjo6SmVl5axzXV1dHDhwgIsXLyr3qKmpcTiwJUlSyka98MIL7Nixg2effZZ77rlHtcDIh42rXQm+BRySJCkOOHT575mwAl+XJCkJWAM8PINg62lJkjIuf27YDLPMzEzc3NymUbXPhK+vL1qtVnXwBgQEKOZLtTYGg4H6+nqnbdzd3fH398disXDu3Dl27dpFfn4+drudZcuWcdddd7FixQoSExPp6+tT7ZOfnx8hISFUVFQoXmqbzUZxcTGHDx9mYmKCDRs2kJWVxbJly2hv7rlqBQAAIABJREFUb6esbHqWrCwAtbW1JCcnk5mZybvvvsuDDz7Irl27yMjI4MyZM0778GHiaoXgbuCVy7+/Atwzs4EkSS2SJJ27/PsAk7nEDrmFbmRotVpWrFhBbm6u0+Vdq9XOy3wZEBBAS0uLU4YKeeZubm52qKLY7XYaGxsZGxtjeHgYi8VCUFAQmzdvZvv27cTGxiqqSmRkJCaTifJydb6zjIwM7HY7RUVFdHd3c+DAASoqKoiKimLHjh3K5nfJkiUsWbKECxcuKKvQTAFYunSpwoH0wgsv8Lvf/Y6WlhY2btzIv//7v6umiH4YuFohCLzMGsHln46zSi5DCBEFZAK5Uw5/5TI/6UuO1Kkp1y5I8q2p2LRpE52dnap1fIODg+nv71eNIwoNDWVwcFCV3jEiIgKr1UpTU5NybGRkhLKyMt59911OnTrFxMQEWq1WMXH6+/vP8glotVoSExPp7OxUjWJ1c3MjMTGRhoYGDh48yMTEBOvXr2flypWzTKfLli3D39+fvLw8WltbHQrAVNx333089dRTpKWl8fTTTxMfH8/vf//767ZXmFMIhBAHhRClDj53X8mDhBBuwOvAI5IkyWlWzwExQAbQAvzC2fWSJP1GkqQVkiStcObl/Khx1113AbBr1y6nbWQ6EjV1JzQ0FCGEqrrj7++P2Wymurqa9vZ2Tp8+ze7duykrK8PDw4Ps7GzuuOMOEhISaG1tVY1bio6Oxmw2U1JS4nT16e7uVnImtFotmzZtIjg42GFbrVZLdna2oh6qCYAkSeTn5+Pj48PBgwd56aWXsFqtfP7znychIYEnnnjiQy/4MacQSJK0TZKkFAeft4G2y5yiXP7pkG1JCKFnUgD+KEnSG1Pu3SZJkk2SJDuTVCyrrsVLfVTIyMggMDCQffv2OW3j5uaGh4eHqhC4uLgQGBhIfX2900FptVrx8PCgo6ODo0eP0tbWRlxcHLfddhsbN24kNDQUjUZDTEwMGo1GNcBPq9WSkpJCT0/PLK4gm81GSUkJhw4dwmq1KoVJzp8/PyeLtpeXl5Lq6WgVgsmNe0tLCykpKXh6evL5z3+e6upqnnzySUZGRvje975HaGgon/rUp/jrX//6oVC+X62f4B3gc8CTl3++PbOBmHzz3wEXJEl6asa5YFmdAj4OlF5lfz5SaDQatmzZwuuvv05/f78SHj0ToaGhVFRUqGZcRUREcPbsWTo6OqblLvf29ipxPPJ+ICAggHXr1jl0wplMJqKjo6muriY5OVnVg1xVVUVxcTEhISEYDAa6u7s5e/Ys/f39REVFkZGRgcFgwG63c/78eSorKx36M+TZvb6+nvj4eNra2jh+/DgrV64kMjJSaTcyMkJhYSF+fn7T/BlGo5FHH32UlStXcuzYMU6cOME777zD3/72N4xGI0lJSSQlJbFq1SoSEhKIjIxU8jrkWstXgqvdEzwJ5AghLgE5l/9GCBEihJAtPdnAPwBbHJhCfyaEKBFCFAObga9dZX8+cnz84x9nfHyc119/3WmbiIgIJElSTckMCwtDr9dTU1ODzWajvr6ew4cPs3//fmprawkLC2Pbtm1ERkbS3d2tqj/LHKhqm18hBJmZmYyNjVFcXKzM/rLuv2rVKkX3j4+PJzQ0lOLi4v/f3pkHR3We+fr5utXaV7QhIZCEJLShFYEEQiBAmN1xxubWTGYymZtUYlMzcc3kZmInuCbxeFKh7Liul6qxY7sysZM4iWNPLgyLAAkEYhNoQQsS2neEdrX2rfvcP6RzrKW7JcASCM5TRVnuPqf7O9X9nv7e7ffO8iNkA6itrSUiIoLY2Fi2bduGu7s7OTk5FBYWYjQaleOMRiPr16+flZHu6+ujo6ODgwcPkpGRQXt7O7/+9a955pln0Ov1fPbZZ/zLv/wLe/fuJTIyEi8vLxwdHc0mEC0h5tKrfxRJSEiQppYbP0oMDQ3h4eHB9u3bLfoGp0+fRqfTKfkFU1y7do3GxkZleoz8IQcGBipN711dXWRkZFgU/AUoKCigqqqKXbt2mf2Fkt9T9kWm3v1nMjY2xrlz5xgcHGT79u24uLjMMoC1a9cqxxsMBgoKCqipqcHT0xNfX18KCwuJjY01mYm+efMmlZWV7N+/f1qxnSRJ/OUvf8HX1xetVkt5eTnNzc10d3fzm9/8Rql8nYkQIk+SpART16yWTXzF2NnZsXXrVuULYq5Hd+XKlZSUlNDX1zdtlpkkSUodj+w32NrakpiYyPLly2ftq5ctW4a3tzcVFRUEBwebrUsKDw+ntraWwsJCUlJSZj1vMBgoLS2loaFBUbOLi4uzKCi2efNmMjMzuXjxItu2baOsrMykAcCXIWQPDw/y8vJob2/H0dGR4ODgWa89Pj5ObW0tfn5+s6pN+/r6GB8fx9vbm8DAwGmVq+np6SbXOhdq2cQC8Ld/+7cMDg7yxz/+0ewxAQEBCCGUGccjIyOUl5dz6tQpLl68SGdnJ2FhYbi5uWEwGPD29jY7EyA8PJzh4WHltUxha2tLREQELS0tswZtd3d3k5GRQVlZGQEBAWzevJnR0VHy8/MtOr8ODg5s2bKF8fFxzpw5Y9YApuLv76+0j/b393PhwoVp0o0wIRQwNjZmsmdCPtbcKNv7QTWCBeC5557DycmJTz/91Owx9vb2LF++nJqaGnJycjh+/DiFhYXKXX///v1ER0cTHh7OwMCAxWiSp6cnHh4elJWVWVSlDgkJwcnJiYKCAgwGAwaDgZKSEjIyMhgZGWHz5s1s2LABHx8fIiIiqK+vt2hYMNEs5O7uzvj4ODqdzmxPskxdXR3t7e3ExMSwbt069Ho9Z8+eJScnRwmFVldX4+zsbLLhp7u7G61Wa3FLd6+o26EFwMbGhl27dnHs2DHa2tpmKdONj4/T2NhIX18fo6OjNDY2EhAQQHBwsKLtL+Pr64uDgwO3b99W8gczEUIQHR3NuXPnqKioICLC9Ng3rVZLfHw8Fy5cID8/n66uLvR6Pf7+/sTFxU3b+4eHh9PR0UF+fr5SFj4TORN89+5dAgICaG5u5vz586SkpMy6DpgYJHLz5k08PDxYs2YNQghWrlxJWVkZVVVVNDQ04OnpSXd3N7GxsSavtbu7G1dX16+0F0H9JVggDh06xOjoKO+9957yWF9fHzdv3uT48ePcuHFD2Xu7ubmRkJBg8ouj0WgICwujq6vL4tALub6nvLzcoqaph4eH0uo5NDREcnIyiYmJs5xfjUZDYmIidnZ2XL582WQ16tRM8Pr165UpOufPn5+15ZoZDZK/4NbW1sTExCiJPblwsLq6mrKysmmZdUmS6O7u/kq3QqAawYKRmppKUFAQv/vd72hububixYucOnWKyspKvLy8SE1NZffu3YSHh9PZ2TlrXzyVgIAApb7H0h49Ojqa8fFxSkpMp1u6u7vJzMxEr9ej1WqxtrY2qz0EE35EcnIyY2NjZGdnKyXg5mqBXF1d2b59Ow4ODly6dImysjJlvXV1dYpinamh5ra2toSGhiJJEp6entjY2FBcXMzJkydJT0+noKCAyspKxsfHzWo33S/qdmiB0Gg0HDx4kCNHjvDb3/6WyMhIIiMjWb169bSIR2BgoNKJlZiYaPK15PqegoICiw30zs7OhISEUFFRoUgqwkTk5/bt25SWlmJjY0NycjJarZaLFy9SVFRkVkMVJjSINm3aRHZ2NpcvX2bz5s3k5+ebLYVwcHBg27Zt3Lhxg+LiYjo6Oli7di03b97E09PTbJMPTBiKJEnExcXh6urKwMAAzc3NtLS0UFNTo/RPFBUVUV9fj4ODA3Z2dlhbW6PT6TAajRYHCppDzRMsIE1NTQQHB5OamsrJkyfN7mPlmPiePXsUoa6ZGAwG0tPT0el07Ny50+yHPTo6Snp6Ora2tqSlpdHb28v169fp6elh1apVxMXFKTmGmzdvUlFRYVFdWqa+vp6cnBxF/t1cLZCMJElUVVVNKybcvXu32euTJIlTp05ha2trMndiMBi4fv06TU1NBAQE0N/fz8DAgDJQBODVV19Fo9FQXFw863xLeQJ1O7SA+Pn58cwzz5CZmWlxfldoaChCiDnreyIjI+np6bGYaZbj+z09PWRnZ3P27FmGh4dJTk4mKSlpmlq1LO1y48YNi9LxMJHldnd3Z2hoSKkonWtcbEhICOHh4RiNRqXUwpzo8N27d+nv7zeZN5Cvf3h4GDc3N8X/2L9/P88++yxf//rXOXDgAFZWVvflMKtGsMC89NJLjI+P88Ybb5g9xs7OjsDAQOrq6ix2WK1atQpXV1eKiooshkIdHR2xsrJStk67du1SqlenotVq2bhxI1qtlsuXL5sdJCI7tZ2dnfj4+NDf3092drbFwSMwEQ2qqKjAw8ODqKgo7t69S3p6OiUlJbNaTKurq7GxsTG5TnkNPT09s5xiIQQ6nQ47OzuEEPe1HVKNYIGJi4tj06ZN/P73v7fYHyCXPMzsyJqKRqMhLi6OwcFBk1KLRqORW7dukZmZqQh5DQ8PW5Q+dHBwYNOmTQwMDHD16tVZfcszSyFSUlJITEyko6OD8+fPmzXaqdGgDRs2EB4ezu7du/H19aW0tJSTJ09y+/ZtxsbGGBgYUIQDzBW/9ff3MzY29pVHhkA1gkXh8OHD9Pb2cuTIEbPHODg4EBQURF1dncWtiaenJ6tWraK8vHxanX1PTw8ZGRncunWLFStWsGfPHjZs2EB3d7dFw5Jfc926dbS2tpKbm6tEdMzVAvn7+5OSksLAwAAZGRkmFTRqa2u5e/cu0dHRih/g4ODAxo0bSUtLw83NjaKiIo4fP87Vq1eRJMli8dtCZIplVCNYBPbu3UtcXBy/+tWvLA7iCA8PR6vVmnTsphITE4NWqyU3N1ep+cnIyGBoaIhNmzaxceNGRQtJbnW01DUGE1GqyMhI6uvrlUpPc8VwMCEIsGPHDrRaLefPn6e6uloxnsHBQQoLC/H09DS5x1+2bBlbtmwhLS0Nb29vRVI+NzeXuro6k2oc3d3daDSarzw8CqoRLBqHDx+mu7ub119/3ewxtra2hIWF0dzcbDExZmdnR0xMDO3t7Zw6dYqSkhJWrFjB7t27FfFcmbi4OFxcXMjJyZlT0SEiIkIR3Dp79uyctUAuLi6kpaXh5eVFXl4eOTk5jIyMcOPGDSRJmpYUM8WyZcsUH8Df35++vj6uX7/O0aNHyc7Opqqqit7eXiVJ5uLisiCqdWqIdJEwGo3ExMTQ3NxMbW2t2Tva+Pj4tFCoqQ/daDRSVlambHPi4+PNRlVgYuZBRkYGTk5ObNu2zWylqfza6enp9Pf34+HhwbZt2+Z0NiVJUtZjZWXF2NjYnGuSOXfuHMPDw+zZsweYkHBpbGxUpCbhy/G1rq6uhIeH4+joqDj/U4mNjQUmQr8zWbBSaiHEMuBPQABQB/wvSZJmpT6FEHVAH2AAxuXFzPf8xwGNRsORI0fYv38/P/3pT3nrrbdMHmdlZUVsbCxXrlyhurp6VnJJr9dz/fp1uru78fHxob29ncbGRoKCgsx+WZ2dnUlKSuLSpUvk5uaSmJhoVoQ3NzeX/v5+XF1dFdGAmJiYOcOhERERuLi4cPnyZWBi+zI6OmqyF0Gmp6eHjo6Oaa/v4eGBh4cHsbGx9Pf309bWxt27d5WegStXrijn63Q6bGxs0Ol0WFlZMT4+/lCiQ/PRHZLZNqktNNUa7+X8Jc++fftITk7mww8/nLPR3tvbm5KSEqWn1mg0UlpaytmzZxkcHGTjxo2kpKQQFxdHe3s75eXlFt/b19eXqKgoGhoaTDrKM0sh0tLSlK3R9evXLardyedXVVWh1WpZvXo1dXV1nDx5kurqarNdb/LxpipPhRA4OTkRFBSkJPK2b99OWloaSUlJREVFKbOTzU3qnC8PWjbxNSB18u+PgSzgpUU8f8nxy1/+kk2bNvGDH/zAbL+BEIJ169Zx+vRp8vPzWbt2rXL39/PzIz4+XulNDggIoKWlheLiYuUuao6wsDD6+/spLS3Fzs5Oqdc3VwsUFxeHra2tYoyyw22KmpoaWltblW1QcHAwBQUF5OXlUVVVRVRUFD4+PsqdenR0lIaGBlauXDnnl1h2it3c3NBqtSYrWgGL2zxLLJbukAScEULkCSG+dx/nLwndofmQlJTEc889x5///Gdl62AKR0dHIiIiaG5u5syZM8rdf9OmTdOa84UQJCQkYG9vz9WrVy0KV8nG5ePjQ35+Po2NjWYNQD4+IiKCDRs20NHRoRTfzWRgYIDCwkK8vLwUw3J1dSU1NZWNGzdiMBi4dOkSGRkZNDc3I0mSIhQwH79BdorvtYF+viyW7lCyJEnxwB4mZBi33OtCl4Lu0Hx56623sLe35/vf/77ZrYJer1fGwgoh2Lp1q9n6HmtrazZt2sTIyAhXr1612HQvi++6u7tz7do1Lly4YFEXCCZ+bVJTUxkfHyczM3OaHpLsRwCzokFyv8Du3btJSEhgdHSUy5cvc+rUKcrKypSZa5ZYqPLpqSyK7pAkSXcm/9sG/IUv9YXmdf7jhq+vLz/84Q8pKCjgP//zP6c9J0d+5L2/7DRaEsYClJ6EtrY2iwp4MLFt2Lx5Mzqdjra2NlauXGmxGA4mHNa0tDRcXV25du0aubm5ishva2sr0dHRZuVcNBoNq1evZs+ePSQlJSGEYHh4WHHy29razBruwMDAgmWKlfU94Pmy7hCY1x1yEEI4yX8DT/GlvtCc5z+u/OQnP2HNmjW88sorSgNKb28v586do7i4GF9fX3bt2kVoaChRUVG0tLQo84nNERAQwJo1a6isrLR4rCyeNTo6iq2tLU1NTRaHkcvY29uTmppKaGgoNTU1nD59mps3b07bBllCo9GwatUqnJ2dsbKyYtWqVTQ1NZGVlcWJEyfIy8ubJTm/kJlimQd1jI8AnwkhvgM0AAdhQncI+EiSpL2AN/CXybuMFfCpJEnpls5/EtDpdHz00Uekpqbyve99jzfeeIOSkhKsrKxISkpi5cqVyp05JCSElpYWbt68ibu7u8UvRHR0NL29veTn52Nvbz9LKnFmKURoaCjZ2dlcu3aNkZGROffoGo2GmJgYli9fTnZ2NkajEUdHRwwGw7wc06GhIWWCTmxsLPHx8dy5c4fGxkZlMqZGo8HDwwNPT0/0ej1CiAXJFMuoybKHzLe//W3+67/+i8OHD7Nnzx7WrVtnUpVueHiYM2fOYGVlxc6dOy0WxY2NjZGVlUVvby/btm1T9t3maoHGx8e5evUqLS0thIeHs3bt2jnj7VVVVeTn5+Ph4UFHRwcODg7Exsbi6+tr8dxbt25x69Yt9uzZM6vDzGAw0NbWRmtrK21tbfT09CjP2dnZ4ezsjKOjo9JMY2dnh06nQ6fTKb6OEOKek2WqETxkenp6iIyMZGRkhJKSErNdY4CiO7pixQrlAzfH0NCQMjtg27ZtODs7W6wFMhqN5OfnU1NTw8qVK1m/fr3ZO/vAwACnT5/G3d2dLVu20N7eTn5+Pr29vXh7exMTE2OyX9poNHLixAlcXFzYsmXu2Mjo6CjHjx/H2dkZJycnRc3b3KSf+22qUdsrHzKurq58/PHH7N69m29961sWxXw9PT2JioqiqKiIsrIys6oS8KUI2Pnz58nKylJGtpqrBdJoNKxbtw5HR0eKioro7+9n06ZNs5xdOaQKkJCQgBACLy8vnnrqKaqqqrh16xZnzpzB39+fiIiIaXf7O3fuMDQ0ZLGdcyqjo6OMj48rShxTHx8eHmZ4eJixsTHGxsYwGAz3HUJVC+geAdLS0jh06BBnzpzh3XfftXisLEBbUlIypzPr6OjIli1bGBsbo6mpiaCgIIvCWEIIwsLC2Lx5M319fZw9e3aWakR1dTVtbW3ExMRMMxCNRsOaNWvYt28fYWFhNDU1kZ6eTk5OjpJbqKqqMumnmEN2imeGUa2trZWxuCtWrCAgIEBR31Y7y5Ywb775JhEREbz00ksUFBSYPU5Oji1btozr168rZcimkCSJ8vJyjEYjWq2W5uZmi409Mr6+vuzcuRN7e3uys7MpLCzEYDDQ399PUVER3t7eZmv/ra2tiY6OZu/evaxZs4ampiZOnz5NZmYmbW1tBAYGzvuL2t3dveBOMahG8MhgbW3NF198gVar5bnnnrM4mEIegmFtbU12drbJqTemaoGEEJw/f97iGFkZJycntm/fzurVqykvLyczM1OZKSZvgywhl3vv37+fqKgoxcmtrq5WNFjnYqEzxTKqETxChIWF8d5771FTU8M3v/lNi8fKg7IlSeLixYvTBLdMlUK4uLiwbds2dDodWVlZFgv4ZKysrEhISCA5OZn+/n66urrw8vIyO1PBFDY2NoSEhCCEwMPDA1dXV0pLSzl16pQy98xUJ91iZIplVCN4xPi7v/s7vvOd73D06FH+/d//3eKxzs7ObN68maGhIS5evMjo6KjFWiBHR0d27NiBs7Mzly9fpqqqal5rkmXXbWxsuHPnDhkZGRZHzM5ErhOKiopiy5Yt7N+/n5iYGGBCQyg9PZ2TJ0+Sl5dHU1MTw8PDDA4OMjo6uihGoIZIH0HkwRi5ubl8/vnnPPPMrKGg02hpaeHy5cu4urri6OhIQ0ODxVqgsbExrl27RktLC8HBwcTGxprdp0uSRFZWFt3d3ezatYuenh7y8/MZGhrC39+f6OjoWfLpM88/e/YskiTx1FNPzVqP3GTf0tJCe3u7oqJhY2PDyMgIQUFB+Pr64uTkhL29vUV/4n6balQjeERpb28nPj4evV7PlStXLEZ1YELoS244CQ0NJTo62uK+3Wg0UlxcTHl5OV5eXiQlJZnc5lRWVlJQUEBCQoLiDI+NjVFWVkZFRYUSFQoNDTWZwOvo6ODcuXPz6jQzGo10dXXR0dFBbW3tLL9BCIG9vT12dnbY2tpOa6iRFf+EEBQVFc16bTVPsATx9PTk2LFjbNmyhT179pCTk4Ovr6/JYyVJmhbKbG9vZ2xszGJXl1z+4OLiQl5eHmfPniUpKWmaHHp/fz/FxcUsX76cwMBA5XGdTkd0dDSBgYEUFxdTWlpKdXU1oaGhBAUFTTOG6upqrKysps0qs7QmuSeira1NmZKp1+vp6+ubpjqn1+sZHR1lbGxMKb4zGo33FSJVfwkecf7nf/6HZ599lpCQEK5duzar1GCmD+Dm5sbVq1dxdnZmy5Yt83Jie3p6uHLlCgMDA0RGRirqcllZWfT09LBr1y6zE3dgYmRUcXExra2tWFtbExISQnBwMJIkcfz4cVavXj3vBJl8TceOHcPX13faJBpzx8oKdwkJEzd6U1W0qgzjEubAgQO88847lJWVsW/fvmmqbzMNYO3ataxYsYLk5GT6+vo4d+6cxaHhMq6uruzcuRM/Pz9KSkrIysri1q1byjANSwYAE8msrVu3smPHDtzd3bl16xbHjx9XCuzmU2E6lcHBQUZGRublFAsh0Gq16HQ6VYHuceaFF17g8OHDZGdn87WvfQ2DwaAUw001ABkfHx+2bt3K6OgomZmZFhNqMjqdjqSkJBITE+np6aG0tBRnZ+c5J89Mxd3dnZSUFHbt2sWqVauUjO+NGzeorq6eU7ZRZjHKp6eiGsES4bXXXuPFF18kPT2d5557juvXr1vUBZInaFpZWXH+/Pl59QsIIZR6fyEEvb29SjXqveDi4qLoCQUEBDA+Pk5eXh7Hjh3j8uXLykwycyxWplhGdYyXEG+//TaDg4N89NFHDA0N8frrr1uMGjk7O7N9+3auXLnClStX5pRTh4loUFdXFwkJCWg0Gm7evMmZM2cICQkhIiLCYgn3VKqrq7G1tVWyy11dXTQ0NNDY2EhzczMajQZPT098fHzw9vZWDA8mjEBuvFkMVCNYYrz//vt0d3fzxRdfYGdnx2effWbxi2lnZ0dqaip5eXmUlpai1+tZv369ychRX1/ftGiQEILly5dTVFREeXk59fX1REVF4e/vbzEK09/fT0tLCxEREcpx7u7uuLu7ExMTQ2dnJ3fu3OHOnTtKTN/GxgZPT0+WLVumqF8vFg8UHZqPeJYQInTyGJnVwL9JkvSWEOJnwHcBWT7iJ5IknZzrfZ+k6JApDAYDhw4d4sMPP2T37t0cPXrUYjgUJpzoiooKioqKsLe3Z9OmTdP23JIkcf78efR6vcloUGdnJwUFBXR1deHs7ExUVJTZBprCwkIqKirYt2/fnE71wMCA0kTT2dk5TavVwcEBFxcXnJycpjXT2NraYm1tPeu9H4oCHV+KZx0RQrw8+f/TdIMkSSoHYicXogWamWi2l/m/kiT98gHX8USh1Wr54IMPsLOz45133mHHjh2cOHHC4lhTIQShoaG4u7tz9epVMjMziYmJITg4GCEElZWVdHR0sH79epNfXHd3d3bs2EFTUxPFxcVcvnyZZcuWERERMU1PyGAwUFtbi6+v75wGABNf9NWrVyuJuNraWm7cuEFgYCBjY2P09vZy9+7dWY348lwCnU6HVqtFo9HctwLdYotv7QCqJUkyP7ZFZd68/fbbuLq68tprr5GUlMTZs2fNDrmQ8fDw4KmnniInJ4eCggKlpbK4uBgfHx+L0SBZQmXFihXU1tZy+/ZtLl26hKurK6Ghofj5+dHY2Mjo6Oi89IRM0d/frwh/yT6BJEkMDQ0xMDDA0NAQw8PDjIyMKA014+PjFpU45uJBjWCaeJYQwqx41iR/DfxhxmP/JIT4eyAX+D/mtEgnRbu+BxMTW1QmePXVV/Hz8+Mf//EfSUxM5OTJk0RHR1s8x8bGhpSUFKqrqyksLKS1tRWNRkN8fPy87qQajYagoCACAwNpaGigrKyMnJwcpVzBwcFh1uzm+dLT0zPLKZbLJeb6ZVkwBbqvSHwLIYQ18DTw5ykPvwcEMbFdagHeNHf+4yS+9VXz3e9+ly+++AK9Xk9ycjL//d//PecmPQERAAALO0lEQVQ5QgiCg4MJCQlBkiQMBgP5+fkW5yfMRKPREBAQwO7du0lJScHe3l65Y8uN+5aEwGaymOXTU1kU8a1J9gD5kiS1TnntVkmSDJIkGYEP+VKUS+UeOXDgANnZ2bi5uXHw4EFeffXVOc/p6+ujsrKS5cuXExsbS3t7O6dPn1a60eaLEAIfHx+cnZ3RaDQEBwfT3t5OdnY2x48fJz8/n46Ojjm3LPJWZ7GN4EG3Q7J41hHmFs/6G2ZshYQQPvJ2Cvg6X4pyqdwHsbGxFBQU8PTTT/Ozn/2M3Nxcfv/735t0mI1GIzdu3ECj0bB+/Xrs7OxYsWIF+fn5FBYWUlNTQ1xcnEX1i6nIArv+/v7Ex8cTExNDS0sLDQ0N1NTUUFVVha2tLb6+vvj6+uLl5TVr+7LYmWKZxRDfQghhD+wEnp9x/utCiFgmBHvrTDyvco+4u7uTlZXF97//fT744ANiYmL4/PPPWbdu3bTj5GjQhg0blH4ABwcHUlJSlPj9xYsX8fHxISoqyqSEylTq6uowGAyKQ6zVavHz88PPz4+xsTGam5u5c+eOYhQajYZly5bh7e2Np6cnbm5uSqZ4rvf6qlGrSB9j/vCHP/DCCy8wOjrKz3/+c/75n/8ZjUZDX18fZ86cwdvbm+TkZJPOsMFgoLKykrKyMsbGxvD39ycyMtLkMG5JkkhPT8fa2podO3ZYXJPBYKC9vV3JDch3f7kQTgihtIO6uLhgY2Mz77Cn2lSjYpLKykqeffZZiouL2b59Ox9//LEy+XLXrl0Wu8JgYptz+/ZtKisrMRqNrFq1ivDw8GlbrNbWVi5cuEBiYuK8+gZmvn5nZyednZ3KMPOp/oiVlZWSKLO1tZ2WLJuqPqfRaEhJSQHuvZRaLZt4zAkJCSEvL4+XX36Zt99+m+joaJ5//nleeOGFOQ0AvpRQCQkJoby8nOrqaurr6/Hz82PNmjW4u7tTVVWlTMu8V6ytrfHx8VEa8GNjY1m1ahU9PT1KI01/fz99fX20t7dbrEQdHx+/r6Ya1QieAHQ6HW+++SZPP/003/zmNzly5AgFBQV8+OGHZmcezMTOzo7Y2FjCwsKoqKigurqapqYmXFxc0Ov1hISEPJA0ylShLVtbW5YvX27SKTcYDNOU5+TOMllbSe0nULHI1q1buXXrFs8//zwZGRlERkby5ptv3lM41NbWlujoaA4cOEB8fLwyU62uro6CgoJ5iXuZoqura15OsVarxcHBAVdXVzw9PfH19cXPz49Vq1ah0WhUI1CZGycnJ95//30uXbqEn58fP/zhD4mOjiYrK+ueXsfKykqpNF22bBnLly+nurqa06dPc+bMGW7fvn1Pibfu7m6cnJwWrXx6KqoRPKEkJSVRXFzML37xC5qamti+fTt/9Vd/RX39/Mu67ty5w8jICBEREWzcuJEDBw4o8i1FRUWcOHGCzMxMysrK0Ov1FpNlDyNTLKMawROMVqvl5ZdfpqKigm984xscPXqUsLAwDh06NC+pxqqqKhwcHJS9u42NDWvWrCEtLY29e/eydu1aJEmiuLiY06dPc/LkSW7cuEFDQ8M0xbyHlSmWUR1jFby8vPjd737HD37wA/71X/+V999/n08//ZRDhw7x4x//2GSbo16vp729naioKJMRGXn6ZkREhDKdpqWlhaamJmpra4GJrZm7u7viUC9WO+VM1DyByiwyMzN56aWXyMvLw9nZmW9/+9u88soruLu7K8fIAz0OHDhwT8O0jUYjPT09tLa2KvkBeeysPMBbbqRxcnLCwcFByRHM5fQ+rKYalceQHTt2kJuby4kTJ3jttdd46623+OCDDzh48CAvv/wyQUFB1NXVzWsQ90zkcompI6QuXLhAb28vgYGB6PV6urq6aGxsnHaeEAJbW9tZiTI5WSaEwGAwPJSmGpXHmH379rFv3z7OnTvHz3/+cz755BN++9vfkpyczI4dO9i8efMDv4cQgr6+Pry8vIiKilIeNxgMDAwMTFOdk32H0dFRBgcHlYYaOU9wvwp0SJK05P6tW7dOUll8SktLpW9961uSg4ODBEgrV66UfvSjH0kNDQ33/ZpDQ0PSn/70J+n27dsPvL6YmBgpJibG5HNArmTm+6RGh1TmTXh4OL/5zW+or6/nP/7jP3B2dub1118nMDCQzZs38+67795zsszcSKbFRDUClXvG3d2dw4cPU1JSwrVr1/jGN77B7du3efHFF/H29iYtLY13332XtjZLPVYTyOp4i10+PRXVCFQeiMTERD755BNaW1s5evQo+/fvJzc3lxdffBFfX1/i4+P50Y9+xIULF6ZNqpeRM8XzFfVaCFQjUPlK0Gq1PP3003z++ed0dnaSnp7OP/zDP9DV1cUbb7xBamoqbm5ubN++nVdeeYWMjAyGhoYeaqZY5kHFtw4CPwPCgQ2SJJkM3gshdgNvA1omOs6OTD4+p3iXKdQ8wdKiqqqKY8eOce7cOXJycpRRTzqdjsDAQMLCwli3bh0JCQkkJiZOy0fcCw+lqUYIEQ4YgV8BPzRlBJOCWxVMtFc2ATeAv5EkqVQI8TrQJX0p3uUmSZIl3SJANYKlTmVlJZmZmWRnZ1NQUEB9fT2Dg4PK8y4uLvj5+eHv74+/vz++vr6sWLFC0TySJR1nlm4/1M4yIUQW5o1gI/AzSZJ2Tf7/jwEkSfqFEKIcSJUmNIt8gCxJkkLnej/VCB4vjEYj5eXl5OTkKE3+DQ0NNDU10dnZabbwzs7ODhsbG6ysrNBqtcrIV1MO+cPOGK8Apqb/moDEyb/nLd6lim89vmg0GsLDwwkPD5/13MjICA0NDdTX19PY2Ehrays9PT3o9Xr0ej1DQ0NK0qy4uPi+RL/mNAIhRAZgSnfjsDShPTTnS5h47J5/fiRJ+gD4ACZ+Ce71fJWliTwHOSQkZMHeY04jkCQp7QHfowmY2sPnB8iTpFtl7aF5iHepqCwIixEivQGECCECJ6UY/5oJ0S74UrwL5hbvUlFZEB7ICIQQXxdCNAEbgRNCiNOTj/sKIU4CSJI0DvwTcBooAz6TJOnW5EscAXYKISqZiB4deZD1qKjcD2o/gcoTgTrCVUXFAqoRqDzxqEag8sSjGoHKE8+SdIyFEO2AKYEcD6BjkZezEDwO1/GoXYO/JEkmRxwtSSMwhxAi11wEYCnxOFzHUroGdTuk8sSjGoHKE8/jZgQfPOwFfEU8DtexZK7hsfIJVFTuh8ftl0BF5Z5RjUDliWdJGoEQYrcQolwIUTXZmzzzeSGEeGfy+SIhRPzDWKcl5nENqUIIvRDi5uS/f3sY67SEEOLXQog2IYTJ+dNL4XMAlp4MIxOKFdXAasAaKAQiZhyzFzjFRFdbEpDzsNd9H9eQChx/2Gud4zq2APFAiZnnH+nPQf63FH8JNgBVkiTVSJI0CvwR+NqMY74GfCJNcA1wnexce1SYzzU88kiSdBHosnDIo/45AEtzO2SqcX/FfRzzMJnv+jYKIQqFEKeEEJGLs7SvlEf9cwCWpjT7fBr3v5Lm/gVkPuvLZ6LepV8IsRf4f8DCdZsvDI/65wAszV8CS43793LMw2TO9UmS1CtJUv/k3ycBnRDCY/GW+JXwqH8OwNI0AkuN+zLHgL+fjE4kAXppUt/oEWHOaxBCLBeTY1eEEBuY+Kzmnqb3aPGofw7AEtwOSZI0LoSQG/e1wK8lSbolhHhh8vn3gZNMRCaqgEHgfz+s9ZpintfwHHBICDEODAF/LU2GXB4VhBB/YCKK5TEpuPBTQAdL43OQUcsmVJ54luJ2SEXlK0U1ApUnHtUIVJ54VCNQeeJRjUDliUc1ApUnHtUIVJ54/j8+PX0NDc0oFAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "from sympde.topology import Square, PolarMapping\n",
    "from sympde.utilities.utils import plot_domain\n",
    "\n",
    "from psydac.api.tests.build_domain import build_pretzel\n",
    "\n",
    "# Define the topological geometry for each patch\n",
    "rmin, rmax = 0.3, 1.\n",
    "\n",
    "# First quarter annulus\n",
    "domain_log_1 = Square('A_1', bounds1=(0., 1.), bounds2=(0., np.pi/2))\n",
    "F_1 = PolarMapping('F_1', dim=2, c1=0., c2=0., rmin=rmin, rmax=rmax)\n",
    "Omega_1 = F_1(domain_log_1)\n",
    "\n",
    "# Second quarter annulus\n",
    "domain_log_2 = Square('A_2', bounds1=(0., 1.), bounds2=(np.pi, 3/2 * np.pi))\n",
    "F_2 = PolarMapping('F_1', dim=2, c1=rmin+rmax, c2=0., rmin=rmin, rmax=rmax)\n",
    "Omega_2 = F_2(domain_log_2)\n",
    "\n",
    "# Join the patches\n",
    "\n",
    "from sympde import Domain\n",
    "connectivity = [((0,1,1),(1,1,-1))]\n",
    "patches = [Omega_1,Omega_2]\n",
    "Omega = Domain.join(patches, connectivity, 'domain')\n",
    "\n",
    "# Example of a complex multi-patch domain\n",
    "# Omega = build_pretzel()\n",
    "\n",
    "# Simple visualization of the topological domain\n",
    "plot_domain(Omega, draw=False, isolines=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Defining the Abstract PDE model using SymPDE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympde.calculus import grad, dot\n",
    "from sympde.calculus import minus, plus\n",
    "\n",
    "from sympde.expr.expr          import LinearForm, BilinearForm\n",
    "from sympde.expr.expr          import integral              \n",
    "from sympde.expr.expr          import Norm                       \n",
    "from sympde.expr               import find, EssentialBC\n",
    "\n",
    "from sympde.topology import ScalarFunctionSpace\n",
    "from sympde.topology import elements_of\n",
    "from sympde.topology import NormalVector\n",
    "\n",
    "# Define the abstract model to solve Poisson's equation using the manufactured solution method\n",
    "x,y       = Omega.coordinates\n",
    "solution  = x**2 + y**2\n",
    "f         = -4\n",
    "\n",
    "V   = ScalarFunctionSpace('V', Omega, kind=None)\n",
    "\n",
    "u, v = elements_of(V, names='u, v')\n",
    "nn   = NormalVector('nn')\n",
    "\n",
    "bc   = EssentialBC(u, solution, Omega.boundary)\n",
    "\n",
    "error  = u - solution\n",
    "\n",
    "I = Omega.interfaces\n",
    "\n",
    "kappa  = 10**3\n",
    "\n",
    "expr_I =- 0.5*dot(grad(plus(u)),nn)*minus(v)  + 0.5*dot(grad(minus(v)),nn)*plus(u)  - kappa*plus(u)*minus(v)\\\n",
    "        + 0.5*dot(grad(minus(u)),nn)*plus(v)  - 0.5*dot(grad(plus(v)),nn)*minus(u)  - kappa*plus(v)*minus(u)\\\n",
    "        - 0.5*dot(grad(minus(v)),nn)*minus(u) - 0.5*dot(grad(minus(u)),nn)*minus(v) + kappa*minus(u)*minus(v)\\\n",
    "        + 0.5*dot(grad(plus(v)),nn)*plus(u)   + 0.5*dot(grad(plus(u)),nn)*plus(v)   + kappa*plus(u)*plus(v)\n",
    "\n",
    "expr   = dot(grad(u),grad(v))\n",
    "\n",
    "a = BilinearForm((u,v),  integral(Omega, expr) + integral(I, expr_I))\n",
    "l = LinearForm(v, integral(Omega, f*v))\n",
    "\n",
    "equation = find(u, forall=v, lhs=a(u,v), rhs=l(v), bc=bc)\n",
    "\n",
    "l2norm = Norm(error, Omega, kind='l2')\n",
    "h1norm = Norm(error, Omega, kind='h1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Discretizing the domain, spaces and equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Discretize the geometry and equation\n",
    "from psydac.api.discretization import discretize\n",
    "from psydac.api.settings import PSYDAC_BACKENDS\n",
    "\n",
    "# The backends are as follows and can all be run with OpenMP by setting backend['omp'] = True.\n",
    "# For all but the Python backend, compilation flags can be accessed and changed using backend['flags'].\n",
    "# - Python backend: Generate and runs python files. Accessed via PSYDAC_BACKENDS['python']\n",
    "# - Pyccel GCC backend: Accessed via PSYDAC_BACKENDS['pyccel-gcc']\n",
    "# - Pyccel Intel backend: Accessed via PSYDAC_BACKENDS['pyccel-intel']\n",
    "# - Pyccel PGI backend: Accessed via PSYDAC_BACKENDS['pyccel-pgi']\n",
    "# - Pyccel Numba backend: Accessed via PSYDAC_BACKENDS['numba']\n",
    "\n",
    "backend = PSYDAC_BACKENDS['python']\n",
    "\n",
    "# Uncomment to use OpenMp\n",
    "# import os\n",
    "# os.environ['OMP_NUM_THREADS'] = \"4\"\n",
    "# backend['omp'] = True\n",
    "\n",
    "ncells = [10, 10]\n",
    "degree = [2, 2]\n",
    "periodic = [False, False]\n",
    "\n",
    "# MPI version\n",
    "# from mpi4py import MPI\n",
    "# comm = MPI.COMM_WORLD\n",
    "# Omega_h = discretize(Omega, ncells=ncells, comm=comm)\n",
    "Omega_h = discretize(Omega, ncells=ncells, periodic=periodic)\n",
    "\n",
    "Vh         = discretize(V, Omega_h, degree=degree)\n",
    "equation_h = discretize(equation, Omega_h, [Vh, Vh], backend=backend)\n",
    "l2norm_h   = discretize(l2norm, Omega_h, Vh, backend=backend)\n",
    "h1norm_h   = discretize(h1norm, Omega_h, Vh, backend=backend)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4: Solving the equation and computing the error norms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "> CG info       ::  {'niter': 109, 'success': True, 'res_norm': 5.704506834276634e-15}\n",
      "> L2 error      :: 9.14e-02\n",
      "> H1 error      :: 1.62e+00\n",
      "> Solution time :: 1.42e+01s\n",
      "> Evaluat. time :: 4.73e-01s \n"
     ]
    }
   ],
   "source": [
    "# Set the solver parameters\n",
    "# 'cg' -> Conjugate gradient method\n",
    "equation_h.set_solver('cg', info=True, tol=1e-14)\n",
    "\n",
    "import time\n",
    "\n",
    "t0_s     = time.time()\n",
    "uh, info = equation_h.solve()\n",
    "t1_s     = time.time()\n",
    "\n",
    "\n",
    "t0_d     = time.time()\n",
    "l2_error = l2norm_h.assemble(u=uh)\n",
    "h1_error = h1norm_h.assemble(u=uh)\n",
    "t1_d     = time.time()\n",
    "\n",
    "print( '> CG info       :: ',info )\n",
    "print( '> L2 error      :: {:.2e}'.format(l2_error))\n",
    "print( '> H1 error      :: {:.2e}'.format(h1_error))\n",
    "print( '> Solution time :: {:.2e}s'.format(t1_s - t0_s))\n",
    "print( '> Evaluat. time :: {:.2e}s '.format(t1_d - t0_d))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 5: Saving the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save the results using OutputManager\n",
    "from psydac.api.postprocessing import OutputManager\n",
    "import os\n",
    "\n",
    "os.makedirs('results_poisson', exist_ok=True)\n",
    "\n",
    "Om = OutputManager(\n",
    "    f'results_poisson/space_info_{Omega.name}',\n",
    "    f'results_poisson/field_info_{Omega.name}',\n",
    "# MPI version\n",
    "#   comm=comm,    \n",
    ")\n",
    "\n",
    "Om.add_spaces(V=Vh)\n",
    "Om.export_space_info()\n",
    "\n",
    "Om.set_static()\n",
    "Om.export_fields(u=uh)\n",
    "\n",
    "Om.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 6: Exporting the results to VTK and visualizing them with Paraview"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export the results to VTK using PostProcessManager\n",
    "from psydac.api.postprocessing import PostProcessManager\n",
    "from sympy import lambdify\n",
    "\n",
    "Pm = PostProcessManager(\n",
    "    domain=Omega,\n",
    "    space_file=f'results_poisson/space_info_{Omega.name}.yaml',\n",
    "    fields_file=f'results_poisson/field_info_{Omega.name}.h5',\n",
    "# MPI version\n",
    "#   comm=comm,    \n",
    ")\n",
    "\n",
    "Pm.export_to_vtk(\n",
    "    f'results_poisson/visu_{Omega.name}',\n",
    "    grid=None,\n",
    "    npts_per_cell=3,\n",
    "    fields='u',\n",
    "    additional_physical_functions={'exact_solution': lambdify(Omega.coordinates, solution, modules='numpy')}\n",
    ")\n",
    "\n",
    "Pm.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Results on the pretzel domain\n",
    "\n",
    "<p float=\"left\">\n",
    "  <img src=\"./paraview_images/poisson_pretzel_u.png\" width=\"32%\" />\n",
    "  <img src=\"./paraview_images/poisson_pretzel_exact.png\" width=\"32%\" />\n",
    "  <img src=\"./paraview_images/poisson_pretzel_error.png\" width=\"32%\" />\n",
    "</p>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "v_psydac",
   "language": "python",
   "name": "v_psydac"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  },
  "vscode": {
   "interpreter": {
    "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}