diff --git a/docs/links.md b/docs/links.md new file mode 100644 index 0000000..c8f0ceb --- /dev/null +++ b/docs/links.md @@ -0,0 +1,5 @@ +Delaynay triangulation + +http://code.activestate.com/recipes/579021-delaunay-triangulation/ + + diff --git a/notebooks/2015-05-31-experimenting-dual-numbers.ipynb b/notebooks/2015-05-31-experimenting-dual-numbers.ipynb index 5a0e894..3a026f7 100644 --- a/notebooks/2015-05-31-experimenting-dual-numbers.ipynb +++ b/notebooks/2015-05-31-experimenting-dual-numbers.ipynb @@ -1,5 +1,29 @@ { "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#Automatic differentiation\n", + "\n", + "Author(s): [Jukka Aho]()\n", + "\n", + "https://github.com/ovainola/JuliaFEM/issues/14\n", + "\n", + "Some tests considering automatic differentiation. Basic question is: *how much do we lose on computation time using automatic differentiation compared to the analytical tangent stiffness matrix*?\n", + "\n", + "It is already known that evaluation of Jacobian takes $n$ function calls, where $n$ is number of unknown parameters. Analytical solution takes only one call, but the function is more cumbersome and needs more cpu time. In linear 2d quadrangle case I assume that we need 8 calls per element. We can expect that autodiffed stiffness matrix evaluation is a bit slower, but interesting question is that are these two strategies even in same decade. FENiCS uses automatic differentiation succesfully so in principle it should work.\n", + "\n", + "**TODO**\n", + "\n", + "- Evaluate 10000 tangent stiffness matrices analytically (nonlinear_stiffness) and 10000 stiffness matrices using automatic differentiation (Kt).\n", + "- @time it.\n", + "- Check that analytical and automatically differentiated matrices match, they should be 1:1.\n", + "- Optimize tangent stiffness matrix code for both cases\n", + "\n", + "ForwardDiff is not behaving nicely (or at least I cannot debug it), so we need to try DualNumbers directly. " + ] + }, { "cell_type": "code", "execution_count": 1, @@ -24,7 +48,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Example: solving scalar equation" + "## Example: solving scalar equation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Find root for $f(x) = -2 + x + x^2$. Roots are $\\{-2, 1\\}$." ] }, { @@ -40,7 +71,7 @@ "iVBORw0KGgoAAAANSUhEUgAAAuwAAAGcCAYAAACC6vWgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XlUVPX/P/DngIBL5r6vbRqS5IZJmn3cUFFHzQUXUEHTEtxF3ApySzSxBFMz3IJQ0UQwZVFzoVzBNAvM+rgrKfpBCBQV7u+PfvmNcMGZe+c9857n4xzO6VzuzH3e82zq5fXO++oURVFARERERERmyUZ0ACIiIiIiejwO7EREREREZowDOxERERGRGePATkRERERkxjiwExERERGZMQ7sRERERERmjAM7EREREZEZ48BORERERGTGOLATEREREZkxDuxERERERGbsmQf2n3/+GQMGDMBLL72EcuXKoUqVKnjzzTcRGRlZbN+0tDR069YN5cuXR5UqVTBs2DBkZmaqEpyIiIiIyBqUetYXXLx4EX/++SdGjBiB2rVrIy8vD1u2bIGXlxfOnz+PWbNmAQAuX76M9u3bo1KlSvj444+Rk5ODTz75BD/99BOOHj0KOzs71U+GiIiIiEg2OkVRFGPfpLCwEC1btsStW7dw4cIFAMDYsWOxYcMGpKeno27dugCAPXv2oEuXLli1ahXeffddYw9LRERERCQ9Ve5ht7GxQd26dYtcNd+6dSt69uz5cFgHgE6dOqFRo0bYvHmzGoclIiIiIpLeM98S87e8vDzk5eXh9u3biI2NRUJCAsLCwgAAV65cwY0bN9CqVatir3NxccGuXbsMT0xEREREZEUMHtgnT56ML7744q83KVUKy5Ytw+jRowEA165dAwDUqlWr2Otq1aqFW7du4f79+7yPnYiIiIjoKQwe2CdNmoSBAwfi6tWriIyMhJ+fH8qUKYPhw4fjzp07AAAHB4dirytdujQA4M6dO48c2DMzM5GQkICGDRuiTJkyhsYjIiIiIjIbd+7cwfnz59G1a1dUrVr1mV5r8MDeuHFjNG7cGADg6emJrl27YuLEifDw8Hg4aOfn5xd73d27dwHgscN4QkICPD09DY1FRERERGS2IiIiMHTo0Gd6jcED+7/169cPSUlJSE9Pf3grzN+3xvzTtWvXUKVKlcfeDtOwYUMAf52Mo6OjWvFIoEmTJmHp0qWiY5BK2Kd82Klc2Kdc2Kc80tLS4Onp+XDWfRaqDex/3wZjY2ODOnXqoFq1ajh27Fix/Y4ePYpmzZo99n3+vvLu6OiIFi1aqBWPBKpQoQK7lAj7lA87lQv7lAv7lI8ht3w/87KON27cKLbt/v372LBhA6pUqQInJycAf11x37FjBy5fvvxwvz179uDs2bMYMGDAMwcly5WRkSE6AqmIfcqHncqFfcqFfRJgwBX20aNHIycnB+3bt0ft2rWRkZGByMhI/Prrr1i7di1sbW0BADNnzkR0dDQ6dOiACRMmICcnB4sXL4azszO8vb1VPxEyX1euXBEdgVTEPuXDTuXCPuXCPgkwYGAfNGgQwsPDsWLFCty8eRPPP/883njjDYSFhaFTp04P96tbty7279+PyZMnY/r06XBwcEDPnj2xZMkSLudoZVq2bCk6AqmIfcqHncqFfcqFfRJgwMDu4eEBDw+PEu3bpEkTxMfHP3MoksvgwYNFRyAVsU/5sFO5sE+5sE8CAJ2iKIroEP+UmpqKli1bIiUlhV+yICIiIiIpGDPjPvOXTomIiIiIyHQ4sJPm+CVjubBP+bBTubBPubBPAjiwkwm4ubmJjkAqYp/yYadyYZ9yYZ8E8B52IiIiIiLN8R52IiIiIiJJcWAnIiIiIjJjHNhJc8nJyaIjkIrYp3zYqVzYp1zYJwEc2MkEFi1aJDoCqYh9yoedyoV9yoV9EsCBnUxg48aNoiOQitinfNipXNinXNgnARzYyQTKli0rOgKpiH3Kh53KhX3KhX0SwIGdiIiIiMiscWAnIiIiIjJjHNhJc/7+/qIjkIrYp3zYqVzYp1ysoc8NGzbAz88PhYWFoqOYLQ7spLn69euLjkAqYp/yYadyYZ9ykb3P/Px8fPDBB7h+/TpsbDiWPo5OURRFdIh/MuaxrURERERkOcLCwjBhwgScPn0ajo6OouNoypgZl3+UISIiIiKTy8vLw/z58+Hp6Sn9sG4sDuxEREREZHJhYWHIzMxEYGCg6ChmjwM7aS49PV10BFIR+5QPO5UL+5SLrH1mZ2cjODgYo0aNwosvvig6jtnjwE6amzZtmugIpCL2KR92Khf2KRdZ+1y6dClyc3Mxa9Ys0VEsAgd20lxYWJjoCKQi9ikfdioX9ikXGfu8efMmlixZgrFjx6Ju3bqi41gEDuykOdmXpLI27FM+7FQu7FMuMva5ePFiFBYWYvr06aKjWAwO7ERERERkEhkZGVi2bBkmTpyI6tWri45jMTiwExEREZFJLFiwAPb29pgyZYroKBaFAztpLjg4WHQEUhH7lA87lQv7lItMfV68eBGrVq2Cv78/KlWqJDqOReHATprLy8sTHYFUxD7lw07lwj7lIlOfc+bMQYUKFTBhwgTRUSyOTlEURXSIfzLmsa1EREREZH7S09Ph5OSEkJAQqx3YjZlxeYWdiIiIiDT1wQcfoG7dunjvvfdER7FIpUQHICIiIiJ5HT9+HFu2bMGaNWvg4OAgOo5F4hV20lxmZqboCKQi9ikfdioX9ikXGfqcOXMmXn31VXh5eYmOYrE4sJPmfHx8REcgFbFP+bBTubBPuVh6n9999x2SkpIwf/58lCrFGzsMxYGdNBcUFCQ6AqmIfcqHncqFfcrFkvtUFAUzZsyAi4sL+vbtKzqOReMfdUhzXO1HLuxTPuxULuxTLpbcZ2xsLI4cOYKkpCTodDrRcSwar7ATERERkaoKCgowc+ZMdOzYEZ07dxYdx+LxCjsRERERqSoyMhK//PIL1q5dKzqKFHiFnTQXHh4uOgKpiH3Kh53KhX3KxRL7zM/PR2BgIPr27YvWrVuLjiMFDuykudTUVNERSEXsUz7sVC7sUy6W2OcXX3yBixcvYt68eaKjSMOggf3YsWPw8/ODk5MTnnvuOTRo0AAeHh44e/Zskf1GjBgBGxubYj+Ojo6qhCfLsHz5ctERSEXsUz7sVC7sUy6W1md2djbmzJmDESNGoEmTJqLjSMOge9iDg4Nx6NAhDBgwAM7Ozrh27RrCwsLQokULHD58GE5OTg/3dXBwKPbXORUqVDAuNRERERGZnSVLliAnJ8eil6M0RwYN7FOmTIGLi0uRBfA9PDzQtGlTLFy4EF999dXD7XZ2dhgyZIjxSYmIiIjIbGVkZGDJkiUYP3486tWrJzqOVAy6JcbV1bXY06pefvllNGnSBOnp6UW2K4qCwsJCZGdnG56SiIiIiMza3LlzYWdnh+nTp4uOIh3VvnSqKAr++OMPVK1atcj2vLw8PP/886hYsSKqVKkCPz8/5ObmqnVYsgB6vV50BFIR+5QPO5UL+5SLpfT522+/4YsvvsCMGTNQuXJl0XGko9o67JGRkbh69WqRbwTXrl0bAQEBaNGiBQoLC7Fr1y58/vnnOHnyJPbt2wdbW1u1Dk9mzM/PT3QEUhH7lA87lQv7lIul9Dl79mzUqFED48aNEx1FSjpFURRj3yQ9PR1vvPEGmjZtioMHDz7x8bMff/wxZs2ahaioKHh4eBT7fWpqKlq2bImUlBSLfhwvERERkTU4fvw4XFxcEB4eDh8fH9FxzJYxM67Rt8RkZGSgR48eqFSpErZs2fLEYR0AJk2aBBsbG+zZs8fYQxMRERGRQIqiICAgAE2aNMGwYcNEx5GWUQP77du30b17d2RnZyM+Ph41a9Z86mtKly6NypUr49atW0/cz93dHXq9vsiPq6srYmJiiuyXmJj4yPu7fH19iy0nmZqaCr1ej8zMzCLbAwMDERwcXGTbxYsXodfri32JNjQ0FP7+/kW25eXlQa/XIzk5ucj2qKgoeHt7F8vm4eHB8+B58Dx4HjwPngfPg+dh8eeRlJSEvXv3oly5csjKyrLY8wDU7SMqKurh/NqjRw/UrFnTqL99MPiWmLt378LNzQ0nTpzA7t278cYbb5TodTk5OahQoQLGjBmDFStWFPs9b4mRT0xMDPr06SM6BqmEfcqHncqFfcrFnPssLCxEq1atULZs2afeEk0CbokpKCiAh4cHjhw5gujo6EcO6/n5+cjJySm2fe7cuQCAbt26GXJoskBRUVGiI5CK2Kd82Klc2KdczLnPqKgonDhxAsHBwRzWNWbQFfaJEydi2bJl6NWrFwYMGFDs956enjh//jyaN2+OIUOGoHHjxgCAhIQE7Nq1C927d8e33377yPfmFXYiIiIi83b37l28+uqraN68ObZt2yY6jkUwZsY1aFnHkydPQqfTIS4uDnFxcUV+p9Pp4OnpiUqVKqFXr15ISkrC+vXrUVBQgFdeeQUff/wxpk6dashhiYiIiMgMhIWF4fLly0hISBAdxSoYNLB/9913T92nQoUK2LBhgyFvT0RERERm6tatW5g/fz5Gjx798C4K0pZqTzolIiIiIvnNnz8fDx48QGBgoOgoVoMDO2nuUUsvkeVin/Jhp3Jhn3Ixtz7PnTuHsLAwTJs2DTVq1BAdx2pwYCfNubm5iY5AKmKf8mGncmGfcjG3PmfNmoUqVapg8uTJoqNYFYPuYSd6FoMHDxYdgVTEPuXDTuXCPuViTn0eP34cUVFRWL16NcqVKyc6jlXhFXYiIiIieiJFUeDv7w8nJyeMGDFCdByrwyvsRERERPREO3fuxL59+7Bjxw6UKsXx0dR4hZ00l5ycLDoCqYh9yoedyoV9ysUc+nzw4AECAgLQoUMHuLu7i45jlTiwk+YWLVokOgKpiH3Kh53KhX3KxRz6XLNmDX7++WcsWrQIOp1OdByrpFMURREd4p+MeWwrmae8vDyULVtWdAxSCfuUDzuVC/uUi+g+s7Oz8corr6Br1658IKaRjJlxeYWdNMf/cciFfcqHncqFfcpFdJ/BwcHIzs7G/PnzheawdhzYiYiIiKiYixcvIiQkBFOmTEG9evVEx7FqHNiJiIiIqJiZM2eiQoUKCAgIEB3F6nFgJ835+/uLjkAqYp/yYadyYZ9yEdXn0aNHERkZiblz56J8+fJCMtD/4cBOmqtfv77oCKQi9ikfdioX9ikXEX0qioIpU6bgtddeg4+Pj8mPT8VxlRgiIiIiemjr1q3o378/EhIS4ObmJjqONLhKDBEREREZLT8/HwEBAejWrRuHdTPCZ8sSEREREQBg+fLlOHfuHLZv3y46Cv0Dr7CT5tLT00VHIBWxT/mwU7mwT7mYss/MzEzMnTsX7777LpycnEx2XHo6DuykuWnTpomOQCpin/Jhp3Jhn3IxZZ9BQUEoLCzEnDlzTHZMKhneEkOaCwsLEx2BVMQ+5cNO5cI+5WKqPn/++WesXLkSCxcuRPXq1U1yTCo5XmEnzXGJMbmwT/mwU7mwT7mYok9FUTB58mQ0bNgQ48aN0/x49Ox4hZ2IiIjIiu3cuROJiYnYtm0bHBwcRMehR+AVdiIiIiIrdf/+fUyePBkdO3ZE7969Rcehx+DATpoLDg4WHYFUxD7lw07lwj7lonWfn3/+OX777TcsXboUOp1O02OR4Tiwk+by8vJERyAVsU/5sFO5sE+5aNlnZmYmgoKCMGrUKDg7O2t2HDKeTlEURXSIfzLmsa1EREREVDJ+fn746quvcPbsWa4MYwLGzLj80ikRERGRleEyjpaFt8QQERERWREu42h5OLCT5jIzM0VHIBWxT/mwU7mwT7lo0WdcXBwSExOxZMkSLuNoITiwk+Z8fHxERyAVsU/5sFO5sE+5qN3n3bt3MWnSJLi5uUGv16v63qQd3sNOmgsKChIdgVTEPuXDTuXCPuWidp9Lly7FxYsXsWPHDi7jaEF4hZ00x9V+5MI+5cNO5cI+5aJmn1euXMH8+fMxbtw4ODo6qva+pD0O7ERERERWYPr06Shbtiw+/PBD0VHoGfGWGCIiIiLJ/fDDD4iIiMDq1atRsWJF0XHoGfEKO2kuPDxcdARSEfuUDzuVC/uUixp9FhYWYvz48WjZsiW8vb1VSEWmxoGdNJeamio6AqmIfcqHncqFfcpFjT7Xrl2LlJQULFu2DLa2tiqkIlPTKYqiiA7xT8Y8tpWIiIiI/k9WVhYaNWoENzc3REREiI5j1YyZcXmFnYiIiEhSc+bMQV5eHoKDg0VHISM888B+7Ngx+Pn5wcnJCc899xwaNGgADw8PnD17tti+aWlp6NatG8qXL48qVapg2LBhfAIbERERkQn8/PPPWLZsGWbNmoU6deqIjkNGeOZVYoKDg3Ho0CEMGDAAzs7OuHbtGsLCwtCiRQscPnwYTk5OAIDLly+jffv2qFSpEj7++GPk5OTgk08+wU8//YSjR4/Czs5O9ZMhIiIiIkBRFPj5+eHFF1/E5MmTRcchIz3zFfYpU6bgwoUL+PTTT+Hj44NZs2bh4MGDePDgARYuXPhwvwULFuDOnTvYu3cv/Pz8MGPGDGzevBknT57EunXr1DwHMnN89LFc2Kd82Klc2KdcDO1z8+bN2LdvH0JDQ+Hg4KByKjK1Zx7YXV1dUapU0QvzL7/8Mpo0aYL09PSH27Zu3YqePXuibt26D7d16tQJjRo1wubNm42ITJbGz89PdARSEfuUDzuVC/uUiyF9/vnnn5gyZQr69u2Lrl27apCKTE2VL50qioI//vgDVatWBfDXo29v3LiBVq1aFdvXxcUFJ06cUOOwZCHc3NxERyAVsU/5sFO5sE+5GNLnvHnzcPPmTYSEhGiQiERQZWCPjIzE1atX4eHhAQC4du0aAKBWrVrF9q1VqxZu3bqF+/fvq3FoIiIiIvr/zpw5g5CQEMycORMNGzYUHYdUYvTAnp6eDl9fX7z55psYPnw4AODOnTsA8Mh7pkqXLl1kHyIiIiIynqIoGDduHOrVqwd/f3/RcUhFRg3sGRkZ6NGjBypVqoQtW7ZAp9MBAMqUKQMAyM/PL/aau3fvFtnncdzd3aHX64v8uLq6IiYmpsh+iYmJj/xChq+vb7HH+aampkKv1xdbWjIwMLDY+qQXL16EXq8vcl8+AISGhhb7EOTl5UGv1yM5ObnI9qioqEc+AtjDw8OqziMmJkaK8wDk6MPY84iJiZHiPAA5+lDjPGJiYqQ4D0COPow9j79/b+nn8TdrP4+/378k5xETE4OkpCQEBwdj4MCBZnUef7P0Pkp6HlFRUQ/n1x49eqBmzZrw8fEplqPEFANlZWUpzZo1U6pWraqkpaUV+d3ly5cVnU6nLF68uNjrPD09lapVqz72fVNSUhQASkpKiqHRyMwMHDhQdARSEfuUDzuVC/uUS0n7zM3NVerXr6/06NFD40RkKGNm3Gdehx346yp5r1698Ntvv2H37t149dVXi/y+Tp06qFatGo4dO1bstUePHkWzZs0MOSxZqE2bNomOQCpin/Jhp3Jhn3IpaZ8ff/wxMjIysHfvXo0TkQjPfEtMQUEBPDw8cOTIEURHR+ONN9545H79+vXDjh07cPny5Yfb9uzZg7Nnz2LAgAGGJyYiIiKih3799VcsWrQIAQEBeOmll0THIQ088xX2KVOmIC4uDr169UJmZiYiIiKK/N7T0xMAMHPmTERHR6NDhw6YMGECcnJysHjxYjg7Oz/yviMiIiIiejaKosDX1xd169bFjBkzRMchjTzzwH7y5EnodDrExcUhLi6uyO90Ot3Dgb1u3brYv38/Jk+ejOnTp8PBwQE9e/bEkiVLYGdnp056IiIiIisWHR2N3bt3Y+fOnU9d0IMs1zPfEvPdd9+hoKAAhYWFxX4KCgqK7NukSRPEx8fjzz//xM2bN7FhwwZUq1ZNtfBkGfg3KnJhn/Jhp3Jhn3J5Up/Z2dmYOHEi+vbti+7du5swFZmaKg9OInoSPnVPLuxTPuxULuxTLk/qMygoCLdv38ann35qwkQkgk5RFEV0iH9KTU1Fy5YtkZKSghYtWoiOQ0RERGR2Tp48iZYtW2LBggWYNm2a6DhUAsbMuLzCTkRERGRBCgsLMXbsWDRu3BgTJ04UHYdMwKB12ImIiIhIjHXr1uGHH37Avn37YG9vLzoOmQCvsJPm/v1IYbJs7FM+7FQu7FMu/+7z5s2bmDZtGjw9PfH2228LSkWmxoGdNLdo0SLREUhF7FM+7FQu7FMu/+5z+vTpePDgARYvXiwoEYnAW2JIcxs3bhQdgVTEPuXDTuXCPuXyzz4PHjyIL7/8EitWrEDNmjUFpiJT4xV20lzZsmVFRyAVsU/5sFO5sE+5/N3nvXv38N5776FNmzYYPXq04FRkarzCTkRERGTmlixZgjNnziA1NRU2Nrzeam3YOBEREZEZ+/333zFnzhxMnjwZzs7OouOQABzYSXP+/v6iI5CK2Kd82Klc2Kdcpk6dCl9fX9SoUQOBgYGi45AgvCWGNFe/fn3REUhF7FM+7FQu7FMumZmZSEhIQFxcHMqVKyc6DgmiUxRFER3in4x5bKslURQFiqLwPjQiIiJ6pKysLLz66qto164dtmzZIjoOGcmYGZfTogCFhYXo3bs318olIiKix5o5cyby8vLw2WefiY5CgnFgF8DGxgaNGzfGRx99hN9++010HCIiIjIzhw8fxsqVKzFv3jzUqVNHdBwSjAO7IEFBQahZsybGjBkDM7srSXXp6emiI5CK2Kd82Klc2Kflu3fvHkaNGoVWrVqhU6dOouOQGeDALki5cuWwatUq7N27F+vXrxcdR1PTpk0THYFUxD7lw07lwj4t36JFi5Ceno7Vq1djxowZouOQGeDALpCbmxs8PT0xZcoUXL9+XXQczYSFhYmOQCpin/Jhp3Jhn5btzJkzmDt3Lvz9/fH666+zTwLAgV24kJAQ6HQ6TJo0SXQUzXCJMbmwT/mwU7mwT8tVWFiI0aNHo379+vjwww8BsE/6Cwd2wapVq4aQkBB8/fXX2LVrl+g4REREJEh4eDgOHDiAVatWoUyZMqLjkBnhwG4GvLy80LlzZ7z//vvIzc0VHYeIiIhM7Nq1a/D394e3tzc6duwoOg6ZGQ7sZkCn02HlypW4fv36w78Ck0lwcLDoCKQi9ikfdioX9mmZxo8fD3t7eyxevLjIdvZJAAd2s/HSSy8hKCgIn376KY4fPy46jqry8vJERyAVsU/5sFO5sE/LExsbiy1btuCzzz5DlSpVivyOfRIA6BQzWwTcmMe2Wrr79++jdevWUBQFx44dg52dnehIREREpKHs7Gw0adIEzs7O+Pbbb6HT6URHIo0YM+PyCrsZsbOzQ3h4OE6fPo1FixaJjkNEREQaCwgIQFZWFlasWMFhnR6LA7uZadGiBfz9/TFnzhykpaWJjkNEREQa2b9/P1auXImFCxeiQYMGouOQGePAboYCAwPRsGFDjBw5EgUFBaLjGC0zM1N0BFIR+5QPO5UL+7QMeXl5GDVqFNq1a4exY8c+dj/2SQAHdrNUunRphIeH49ChQ1i+fLnoOEbz8fERHYFUxD7lw07lwj4tQ2BgIC5duoQvv/wSNjaPH8fYJwEc2M1Wu3bt4OvrixkzZuD8+fOi4xglKChIdARSEfuUDzuVC/s0f8eOHUNISAiCgoLQuHHjJ+7LPgngKjFmLScnB6+99hoaNWqExMREfhmFiIjIwt27dw8tW7aEvb09jhw5glKlSomORCbCVWIkVb58eaxatQq7d+/GunXrRMchIiIiI3388cdIT0/HmjVrOKxTiXFgN3PdunXDsGHDMHnyZFy7dk10HCIiIjLQ6dOnMX/+fEyfPh2vv/666DhkQTiwW4ClS5fCwcEB7733HszsDqYSCQ8PFx2BVMQ+5cNO5cI+zdODBw/g4+ODl19+GbNnzy7x69gnARzYLULlypWxcuVKxMbG4uuvvxYd55mlpqaKjkAqYp/yYadyYZ/mKSQkBCkpKQgPD4eDg0OJX8c+CeCXTi3KkCFDEB8fj19++QU1a9YUHYeIiIhKIC0tDc2bN8e4ceOwePFi0XFIEH7p1EosW7YMdnZ2FntrDBERkbV58OABRowYgYYNG2LOnDmi45CF4sBuQapWrYqVK1di+/btiIqKEh2HiIiIniIkJATHjx/H2rVrUaZMGdFxyEIZNLDn5uYiMDAQ3bp1Q+XKlWFjY4P169cX22/EiBGwsbEp9uPo6Gh0cGvVt29fDBo0COPGjUNGRoboOERERPQYaWlp+PDDDzF58mS4urqKjkMWzKCB/caNG5g7dy7OnDmDZs2aAcBjH+rj4OCAiIiIIj+ffPKJ4YkJoaGhsLW1xfvvv28Rt8bo9XrREUhF7FM+7FQu7NM8qHUrDPskADBoxf7atWsjIyMD1atXR0pKClxcXB67r52dHYYMGWJwQCquatWqWLFiBfr3749NmzZh0KBBoiM9kZ+fn+gIpCL2KR92Khf2aR7+vhXm+++/N+pWGPZJgIFX2O3t7VG9enUAeOoVXkVRUFhYiOzsbEMORY/Rr18/eHh4wNfX1+xvjXFzcxMdgVTEPuXDTuXCPsX7560wbdq0Meq92CcBJvjSaV5eHp5//nlUrFgRVapUgZ+fH3Jzc7U+rFUIDQ1FqVKlMGbMGIu4NYaIiEh2XBWGtGDQLTElVbt2bQQEBKBFixYoLCzErl278Pnnn+PkyZPYt28fbG1ttTy89KpVq4ZVq1ahb9++2LBhA4YPHy46EhERkVULDg5W5VYYon/S9Ar7ggULsGDBAvTv3x8DBw7E2rVrMX/+fHz//ffYsmWLloe2Gn369MGwYcMwfvx4XLx4UXScR4qJiREdgVTEPuXDTuXCPsX58ccf8dFHH2H69OlG3wrzN/ZJgIB12CdNmgQbGxvs2bPH1IeW1meffYbnn38eI0eORGFhoeg4xXDNeLmwT/mwU7mwTzHy8/Ph5eUFR0dHBAYGqva+7JMAAQN76dKlUblyZdy6deuJ+7m7u0Ov1xf5cXV1LfYnzcTExEcueeTr64vw8PAi21JTU6HX65GZmVk3YVBXAAAgAElEQVRke2BgIIKDg4tsu3jxIvR6PdLT04tsDw0Nhb+/f5FteXl50Ov1SE5OLrI9KioK3t7exbJ5eHioeh4VK1ZEeHg4du/ejXfeecfszmPTpk1W1Yfs57Fp0yYpzgOQow81zmPTpk1SnAcgRx/GnsemTZukOI+/Wcp5BAUF4cyZM3j//ffRv39/1c7j7z7Zh2WdR1RU1MP5tUePHqhZsyZ8fHyK5SgpnWLktxWPHz+O1q1bY926dRg2bNhT98/JyUGFChUwZswYrFixotjvU1NT0bJlS6SkpKBFixbGRLM6Y8eOxfr16/Hjjz/ilVdeER2HiIjIKvzwww946623MG/ePMyYMUN0HDJTxsy4ml1hz8/PR05OTrHtc+fOBQB069ZNq0NbrUWLFqFWrVoYMWIECgoKRMchIiKSXm5uLoYPH47WrVsXu0JMpBaDV4kJCwtDVlYWrl69CgCIjY19+KXH8ePH49atW2jevDmGDBmCxo0bAwASEhKwa9cudO/eHb1791YhPv3Tc889h3Xr1qF9+/ZYsmQJpk2bJjoSERGR1AICAnDlyhXs3LkTpUppuvgeWTPFQA0bNlR0Op2i0+kUGxsbxcbG5uE/X7hwQcnKylK8vLyUV155RSlXrpxSunRppWnTpsrChQuVBw8ePPZ9U1JSFABKSkqKodGs3tSpUxV7e3vl5MmToqMoiqIoI0aMEB2BVMQ+5cNO5cI+TScpKUkBoISGhmp2DPYpD2NmXIP/KHju3Lmn7rNhwwZD356MMHfuXCQkJMDT0xNHjx5F6dKlhebhU9rkwj7lw07lwj5N43//+x+8vb3RqVMnjB07VrPjsE8CBKwSQ9orXbo0IiIicObMGcyePVt0HAwePFh0BFIR+5QPO5UL+zSNsWPH4s8//8TatWthY6PdOMU+CeDALi1nZ2fMnz8fISEh+O6770THISIiksbXX3+NjRs3YsWKFahXr57oOGQFOLBLbNKkSWjfvj2GDx+OrKws0XGIiIgs3sWLFzF27FgMGTIEgwYNEh2HrAQHdonZ2tpi/fr1uH37Nvz8/ITl+PcDD8iysU/5sFO5sE/tFBYWYvjw4Xj++eexfPlykxyTfRLAgV16DRo0wPLlyxEZGfnwaWmmtmjRIiHHJW2wT/mwU7mwT+2EhIRg//79WL9+PSpWrGiSY7JPAlR40qna+KRT9SmKgkGDBiExMRE//fQT6tata9Lj5+XloWzZsiY9JmmHfcqHncqFfWrj1KlTcHFxwbhx4/DJJ5+Y7LjsUx5m+aRTMh86nQ4rVqxAuXLlMGLECBQWFpr0+PwPjVzYp3zYqVzYp/ru3r2LoUOHonHjxpg/f75Jj80+CeDAbjUqV66MdevWYc+ePQgJCREdh4iIyGLMnDkTv/76KyIjI+Hg4CA6DlkhDuxWpHPnzpg6dSpmzpyJ1NRU0XGIiIjMXkJCApYuXYqFCxeiadOmouOQleLAbmXmz5+P1157DUOGDEFubq5Jjunv72+S45BpsE/5sFO5sE/1XL9+HcOHD0fXrl0xYcIEIRnYJwEc2K2Ovb09oqKicOnSJUyaNMkkx6xfv75JjkOmwT7lw07lwj7VoSgKvL29UVhYiHXr1mn6NNMnYZ8EcJUYq7V69WqMHj0a33zzDfr27Ss6DhERkVkJDQ3F+PHj8e2338Ld3V10HJIAV4mhZzZq1Cj07dsXo0aNwuXLl0XHISIiMhunTp2Cv78/xo8fz2GdzAIHdiul0+mwevVqlClTBsOGDUNBQYHoSERERMLduXMHgwcPRqNGjRAcHCw6DhEADuxWrUqVKvjqq6+wb98+TR8CkZ6ertl7k+mxT/mwU7mwT+NMnToV//3vfxEVFYXSpUuLjsM+CQAHdqvXoUMHBAQEYPbs2Thy5Igmx5g2bZom70tisE/5sFO5sE/DxcbG4vPPP0dISAicnJxExwHAPukv/NIp4f79+2jfvj0yMjJw4sQJVKxYUdX3v3jxIr/lLhH2KR92Khf2aZhLly6hWbNmeOutt7Bt2zbodDrRkQCwT5nwS6dkFDs7O0RFReF///sfRo8eDbX/DMf/0MiFfcqHncqFfT67Bw8eYMiQIShXrhzWrFljNsM6wD7pLxzYCQDQsGFDfPnll4iOjsbq1atFxyEiIjKZOXPm4NChQ4iKikLlypVFxyEqhgM7PdS/f3+89957mDBhAk6fPi06DhERkeb27t2LefPm4aOPPkLbtm1FxyF6JA7sVERISAhefvlleHh4IC8vT5X35LJYcmGf8mGncmGfJXf9+nUMHToUHTt2xPTp00XHeST2SQAHdvqXMmXKYNOmTTh37hwmTJigynuqNfiTeWCf8mGncmGfJVNYWIjhw4ejsLAQERERsLW1FR3pkdgnAVwlhh4jPDwco0aNQlRUFAYNGiQ6DhERkaoWL16MadOmIT4+Hl27dhUdh6wAV4kh1fn4+GDw4MEYPXo0zp49KzoOERGRag4fPoyZM2ciICCAwzpZBA7s9Eg6nQ6rVq1CrVq1MGDAANy5c0d0JCIiIqPdvHkTHh4ecHFxwdy5c0XHISoRDuz0WOXLl0d0dDTOnDlj1P3smZmZKqYi0dinfNipXNjn4xUWFmLYsGHIzc3Fpk2bYGdnJzrSU7FPAjiw01M4OzsjNDQUq1evRmRkpEHv4ePjo3IqEol9yoedyoV9Pt7ixYuxc+dOfPXVV6hXr57oOCXCPgngwE4lMHLkSHh5eWHMmDFIS0t75tcHBQWpH4qEYZ/yYadyYZ+PdvDgQcyaNQszZsxA9+7dRccpMfZJAFeJoRLKzc2Fi4sLbGxscOTIEZQrV050JCIiohK5fv06mjdvjpdffhl79uxBqVKlREciK8RVYkhz5cqVQ3R0NM6dOwc/Pz/RcYiIiEqkoKAAQ4cOxYMHDxAVFcVhnSwSB3YqMScnJ3z++edYt24d1q5dKzoOERHRU82fPx979uzB119/jdq1a4uOQ2QQDuz0TIYPH46RI0di7Nix+PHHH0v0mvDwcI1TkSmxT/mwU7mwz/+ze/duBAUFITAwEJ06dRIdxyDskwAO7GSA0NBQODo6ol+/fvjf//731P1TU1NNkIpMhX3Kh53KhX3+5dKlSxg8eDA6d+6M2bNni45jMPZJAL90SgY6d+4cWrZsibZt22L79u2wseGf/YiIyDzk5+ejffv2yMjIQEpKCqpWrSo6EhG/dEqm98ILLyAiIgI7duzAwoULRcchIiJ6aNKkSfjxxx+xZcsWDuskBQ7sZDB3d3d88MEH+OCDD7B7927RcYiIiLBhwwasWLECoaGhcHFxER2HSBUc2MkogYGB6Ny5MwYPHoxLly6JjkNERFbs5MmTGDNmDLy9vfHuu++KjkOkGg7sZBRbW1tERkaibNmy6N+/P/Lz84vto9frBSQjrbBP+bBTuVhrn1lZWejXrx8cHR2xfPly6HQ60ZFUYa19UlEGDey5ubkIDAxEt27dULlyZdjY2GD9+vWP3DctLQ3dunVD+fLlUaVKFQwbNgyZmZlGhSbzUrVqVWzZsgU//vgjJkyYUOz3fNCSXNinfNipXKyxz8LCQgwbNgw3b97Eli1bUKZMGdGRVGONfVJxBg3sN27cwNy5c3HmzBk0a9YMAB75J9nLly+jffv2+O9//4uPP/4YU6dOxbfffosuXbrg/v37xiUns+Li4oLly5dj1apV+PLLL4v8zs3NTVAq0gL7lA87lYs19jlv3jzs2LEDERERePHFF0XHUZU19knFGfR83tq1ayMjIwPVq1dHSkrKY7/UsWDBAty5cwcnTpxA3bp1AQCtW7dGly5dsG7dOt5fJplRo0YhJSUFvr6+aNq0Kd544w3RkYiISHJxcXEIDAzEnDlz0KNHD9FxiDRh0BV2e3t7VK9eHQDwpGXct27dip49ez4c1gGgU6dOaNSoETZv3mzIocnMffbZZ2jVqhXeeecdZGRkiI5DREQSO3PmDDw9PdGnTx/MmjVLdBwizWj2pdMrV67gxo0baNWqVbHfubi44MSJE1odmgSyt7fHli1boCgKBgwYgHv37iEmJkZ0LFIR+5QPO5WLtfSZnZ2NPn36oE6dOtiwYYO0D/Czlj7pyTT7t/vatWsAgFq1ahX7Xa1atXDr1i3exy6pWrVqYcuWLThy5AgmT56MqKgo0ZFIRexTPuxULtbQZ2FhIby8vHD16lXExMSgfPnyoiNpxhr6pKfTbGC/c+cOAMDBwaHY70qXLl1kH5LPm2++ibCwMCxfvhzdunUTHYdUtGnTJtERSGXsVC7W0Oe8efMQFxeHyMhINGrUSHQcTVlDn/R0mg3sfy+p9Kh1ue/evVtkn0dxd3eHXq8v8uPq6lrsr4YSExMfuUapr68vwsPDi2xLTU2FXq8vtqxkYGAggoODi2y7ePEi9Ho90tPTi2wPDQ2Fv79/kW15eXnQ6/VITk4usj0qKgre3t7Fsnl4eFjFeYwePRrvvvsu3nvvPfTv399iz+OfLLkPngfPg+fB85DhPP7+kmmjRo3Qpk0biz0PQI4+eB6PPo+oqKiH82uPHj1Qs2ZN+Pj4FMtRUjrlSd8aLYHjx4+jdevWWLduHYYNG/Zw+5UrV1CvXj0sWrQIU6dOLfIaLy8vxMfH48aNG8XeLzU1FS1btkRKSgpatGhhTDQyA/n5+ejQoQPOnz+PY8eOoU6dOqIjERGRhUpLS0ObNm3QsWNHbN26Vdr71klOxsy4mv2bXqdOHVSrVg3Hjh0r9rujR48+XL+d5Obg4IBvvvkGtra26Nu3L2+DIiIig9y6dQu9evVC/fr1pf6SKdGjaPpve79+/bBjxw5cvnz54bY9e/bg7NmzGDBggJaHJjMyY8YMbN++HadPn8a77777xKVAyfw96q8hybKxU7nI2Of9+/cxcOBAZGVlITY2Vuovmf6bjH3SszPowUkAEBYWhqysLFy9ehUAEBsbi4sXLwIAxo8fj+effx4zZ85EdHQ0OnTogAkTJiAnJweLFy+Gs7Mz/wW0Im5ubmjRogXWrl2LQYMGoWnTpggICBAdiwzEp+7Jh53KRcY+J0+ejP3792P37t144YUXRMcxKRn7pGdn8D3sL7zwAi5cuPDXm+h0AP56iJJOp8O5c+dQv359AMAvv/yCyZMnIzk5GQ4ODujRoweWLFmCatWqPfJ9eQ+73GbPno0FCxYgNjYWPXv2FB2HiIjM3BdffIExY8Zg5cqVGDNmjOg4RAYzZsY1+Ar7uXPnSrRfkyZNEB8fb+hhSDJz5szB6dOnMWTIEBw+fBhNmjQRHYmIiMzU/v374evrC19fXw7rZNX4jQ0yKRsbG3z11Vdo0KAB9Ho9bt68KToSERGZoXPnzqFfv35o3749li5dKjoOkVAc2Elz/14/tXz58oiNjcXt27fRr18/3Lt3T1AyMsS/+yTLx07lIkOf2dnZ0Ov1qFixIqKjo2FnZyc6kjAy9EnG48BOmlu0aFGxbS+88AJiYmJw6NAhvPfee1w5xoI8qk+ybOxULpbe54MHD+Dh4YFLly4hLi4OlStXFh1JKEvvk9TBgZ00t3Hjxkdub9u2LcLDw7F27Vr+B8mCPK5PslzsVC6W3uekSZOQlJSELVu2wNHRUXQc4Sy9T1KHwV86JSqpsmXLPvZ3np6eOHPmDKZPn45XXnkF77zzjgmTkSGe1CdZJnYqF0vuMzQ0FGFhYVi5ciU6d+4sOo5ZsOQ+ST28wk7CffTRRxg4cCA8PT2RkpIiOg4REQmwc+dOTJw4EZMnT+aKMET/woGdhLOxscG6devQtGlT9OrVq8iTcYmISH4//fQTBg0ahJ49e/IWSaJH4MBOmvP393/qPmXKlMH27dthZ2eHXr164c8//zRBMjJESfoky8JO5WJpfWZkZKBnz5546aWXEBkZCVtbW9GRzIql9Una4MBOmvv7qbdPU7NmTezYsQO///47PDw88ODBA42TkSFK2idZDnYqF0vqMzc3F71798b9+/cRFxeH5557TnQks2NJfZJ2dIqZradnzGNbSQ5JSUlwd3eHj48PVq5cCZ1OJzoSERGprKCgAO+88w727NmDAwcO8P/5JD1jZlxeYSez06VLF3zxxRf44osvEBwcLDoOERGpTFEUTJgwAd9++y2io6M5rBM9BZd1JLPk7e2N8+fPY8aMGWjQoAEGDx4sOhIREalkyZIlWL58OVatWoXu3buLjkNk9niFnTSXnp5u0OuCgoIwbNgwjBgxAvv371c5FRnK0D7JfLFTuZh7n5s3b4a/vz9mzpyJ0aNHi45j9sy9TzINDuykuWnTphn0Op1Oh9WrV6Ndu3bo06cP0tLSVE5GhjC0TzJf7FQu5tznwYMH4eXlhaFDh2LevHmi41gEc+6TTIcDO2kuLCzM4Nfa29vjm2++Qd26ddG9e3dcvXpVxWRkCGP6JPPETuVirn2eOXMGvXv3Rtu2bbFmzRouKFBC5tonmRYHdtKcsUtSVahQATt37sSDBw/g7u6O27dvq5SMDMElxuTDTuVijn1evXoVXbt2Ra1atfDNN9/A3t5edCSLYY59kulxYCeLUK9ePSQkJODChQvo06cP8vPzRUciIqISyMrKQvfu3VFQUID4+HhUrFhRdCQii8OBnSyGk5MT4uLicPjwYXh5eaGgoEB0JCIieoK7d++id+/euHTpEuLj41GvXj3RkYgsEgd20pyaa6m3a9cOUVFR2Lp1KyZOnAgze+6XVeDa+PJhp3Ixlz4LCgowdOhQHD16FDt27ICTk5PoSBbJXPoksTiwk+by8vJUfb8+ffrg888/R1hYGBYuXKjqe9PTqd0nicdO5WIOfSqKAj8/P2zfvh2bN2/Gm2++KTqSxTKHPkk8nWJmlyiNeWwrWZegoCB89NFHWLNmDby9vUXHISKi/2/OnDkIDAzEl19+iZEjR4qOQ2QWjJlx+aRTsliBgYHIyMjAqFGjUKlSJfTp00d0JCIiq7dq1SoEBgZi3rx5HNaJVMJbYshi6XQ6LF++HP369YOHhwf27NkjOhIRkVXbuHEj3n//fYwbNw4zZ84UHYdIGhzYSXOZmZmavbetrS0iIiLQoUMH9O7dG0ePHtXsWPQXLfskMdipXET1uXPnTnh5ecHLywuffvopH4ykEn4+CeDATibg4+Oj6fvb29tj69ateP3119G9e3f8/PPPmh7P2mndJ5keO5WLiD4PHjyIfv36oUePHggPD4eNDccLtfDzSQAHdjKBoKAgzY9Rrlw57NixA3Xr1oWbmxvOnTun+TGtlSn6JNNip3IxdZ+pqano2bMn3nzzTWzcuBGlSvHrcWri55MADuxkAqZa7adSpUpITExE2bJl0aVLF2RkZJjkuNaGqzfJh53KxZR9pqeno2vXrnj11VcRExOD0qVLm+zY1oKfTwI4sJNkatSogaSkJNy9exddunTBzZs3RUciIpLShQsX0KVLF9SoUQO7du1C+fLlRUcikhYHdpJOw4YNkZSUhD/++ANubm7IysoSHYmISCpXrlxBx44dYW9vj8TERFSuXFl0JCKpcWAnzYWHh5v8mI6Ojti9ezfOnTuH7t27Iycnx+QZZCWiT9IWO5WL1n3+8ccf6NSpE+7fv4+9e/eidu3amh7P2vHzSQAHdjKB1NRUIcd1dnZGYmIifvnlF/Tq1YuPd1aJqD5JO+xULlr2efPmTXTp0gXZ2dnYu3cvGjRooNmx6C/8fBIA6BRFUUSH+CdjHttK9Cg//PAD3Nzc0LZtW8TGxsLBwUF0JCIii5OVlYVOnTrh0qVL2L9/PxwdHUVHIrIoxsy4vMJO0nvzzTcRFxeHAwcOYODAgbh//77oSEREFiUnJwfdu3fH+fPnsXv3bg7rRCbGgZ2sQocOHbBt2zbs2rULgwcP5tBORFRCubm56NmzJ3755RckJibC2dlZdCQiq8OBnaxGt27dsHXrVsTGxmLIkCEc2omInuLvYT01NRXx8fFo2bKl6EhEVokDO2lOr9eLjvBQr169sGXLFmzfvh1Dhw7l0G4Ac+qT1MFO5aJWn38P68ePH0d8fDxcXV1VeV96Nvx8EsCBnUzAz89PdIQi9Ho9oqOjERMTg6FDh+LBgweiI1kUc+uTjMdO5aJGn/8e1tu2batCMjIEP58EaDyw79u3DzY2No/8OXr0qJaHJjPi5uYmOkIxvXv3RnR0NLZt28ah/RmZY59kHHYqF2P7zM3NRY8ePTismwl+PgkASpniIBMmTICLi0uRbS+99JIpDk30WH8P7QMGDICnpyciIiJQqpRJPhJERGbp72E9JSWFwzqRGTHJdPLWW2/hnXfeMcWhiJ5Jnz59sHnzZgwcOBAFBQX4+uuvYWdnJzoWEZHJ/fnnn+jZsyeHdSIzZJJ72BVFQU5ODm87sFIxMTGiIzxR3759H64e079/f+Tn54uOZNbMvU96duxULob0efv2bbi5uT1cDYbDuvng55MAEw3s3t7eqFChAsqUKYOOHTsiJSXFFIclMxEVFSU6wlPp9Xps374diYmJ6N27N/Ly8kRHMluW0Cc9G3Yql2ft8+bNm+jUqRPS0tKwe/duDutmhp9PAgCdoiiKVm9+6NAhLF26FO7u7qhatSp+/vlnfPLJJ8jNzcUPP/yAZs2aFXuNMY9tJTLW3r170atXL7Ru3RpxcXF47rnnREciItLM9evX0blzZ1y7dg1JSUmP/P8yEanDmBlX03vYXV1di6zb2rNnT/Tv3x/Ozs6YMWMGdu3apeXhiZ5Zx44dkZCQAHd3d3Tt2hU7d+5EhQoVRMciIlLdlStX0LlzZ2RlZWH//v1o0qSJ6EhE9BgmX4f9pZdeQu/evfHdd99Bw4v7RAZr164ddu/ejV9++QWdO3fGrVu3REciIlLVhQsX8PbbbyM3NxcHDhzgsE5k5oQ8OKlu3bq4d+8ecnNzH7uPu7s79Hp9kR9XV9diX75ITEx85FPAfH19ER4eXmRbamoq9Ho9MjMzi2wPDAxEcHBwkW0XL16EXq9Henp6ke2hoaHw9/cvsi0vLw96vR7JyclFtkdFRcHb27tYNg8PD56HmZ9HVlYWnJ2dcf78ebRv3x5Xr161yPOQpQ+eB8+D56HeeXTs2BFvvvkmCgsLceDAAbzyyisWeR6y9MHzkPM8oqKiHs6vPXr0QM2aNeHj41MsR0lpeg/74/Tv3x+7du165MDOe9jl4+3tjbVr14qOYZC0tDS4ubnBzs4OSUlJfH4ALLtPejR2Kpcn9Zmamopu3bqhWrVqSExMRJ06dUycjp4VP5/yMGbG1fQK+40bN4ptO3nyJGJjY/nkLitiyV07OjoiOTkZpUqVQrt27fDTTz+JjiScJfdJj8ZO5fK4Pg8cOID//Oc/eOGFF3DgwAEO6xaCn08CNL7C3rFjR5QtWxaurq6oXr06fvnlF3zxxRdwcHDAoUOH0Lhx42Kv4RV2Mkd//PEHunXrhgsXLmDnzp1o06aN6EhERCW2Y8cODBgwAG3btsW2bdtQvnx50ZGIrI7ZXmHv27cvMjMzsXTpUvj6+iI6Ohr9+/fH8ePHHzmsE5mrGjVq4LvvvoOTkxM6deqEpKQk0ZGIiEokMjISffr0gbu7O7799lsO60QWSNOBfdy4cTh8+DAyMzNx7949XL58GevXr8eLL76o5WGJNFGxYkUkJCTg7bffRo8ePbBp0ybRkYiInig0NBSenp4YPnw4Nm3aBAcHB9GRiMgAQlaJIevy7293W7KyZcsiJiYGHh4eGDRoED799FPRkUxOpj7pL+xULsnJyVAUBTNmzMD48eMxdepUfPnllyhVStNHr5BG+PkkgAM7mcCiRYtER1CVvb091q9fj2nTpmHSpEnw9/dHYWGh6FgmI1ufxE5ls3DhQgwfPhwLFy7EkiVLsHjxYuh0OtGxyED8fBKg8ZNOiQBg48aNoiOozsbGBsHBwahTpw4mTpyIa9euYc2aNbC3txcdTXMy9mnt2Kk8cnJycPfuXWzatAlRUVEYNGiQ6EhkJH4+CeDATiZQtmxZ0RE0M378eNSsWRNeXl64fv06tm7dKv0XumTu01qxUzlkZGTA3d0dv//+O+Lj49GhQwfRkUgF/HwSwFtiiIw2cOBAJCQk4MiRI3j77bcfPhWViMhUfv31V7i6uuKPP/7AwYMHOawTSYYDO5EK/vOf/+DgwYO4ceMG2rRpg1OnTomORERW4sCBA3B1dUWZMmVw6NAhODs7i45ERCrjwE6a8/f3Fx3BJJydnXH48GFUrVoV7dq1Q0JCguhImrCWPq0JO7VcX331FTp37oxmzZrh+++/R/369dmnZNgnARzYyQTq168vOoLJ1KlTBwcOHHi4VvuqVatER1KdNfVpLdip5VEUBYGBgRg2bBi8vLywa9cuVKpUCQD7lA37JADQKYqiiA7xT8Y8tpXIXBQUFGDSpEkIDQ3F1KlTERwcDBsb/vmYiIyXn5+PkSNHIjIyEgsWLMD06dO5bCORBTBmxuUqMUQasLW1xbJly/Dyyy9j4sSJ+P333/HVV1+hXLlyoqMRkQXLzMxE3759cezYMWzatAkDBw4UHYmITICX/Ig0NH78eMTExCAxMRFt27bFhQsXREciIgt1+vRptG7dGmfOnMG+ffs4rBNZEQ7spLn09HTREYTS6/U4dOgQbt++DRcXF3z//feiIxnF2vuUETs1f7GxsXB1dUX58uVx7NgxtGnT5rH7sk+5sE8COLCTCUybNk10BOGaNm2Ko0ePwtHRER06dMCaNWtERzIY+5QPOzVfiqJg/vz56NOnD9zc3PD999+jQYMGT3wN+5QL+ySAAzuZQFhYmOgIZqFatWpISkqCt7c3Ro4ciUmTJuHBgweiYz0z9ikfdmqe8vLyMHjwYMyePRsffvghoqOj8dxzzz31dexTLuyTAH7plEyAS1L9H3t7e6xcuRJNmzbFxIkT8a20SzkAACAASURBVPPPPyMqKgpVqlQRHa3E2Kd82Kn5uXTpEvr06YP09HRER0ejf//+JX4t+5QL+ySAV9iJTE6n08HPzw+JiYkPl3hKTU0VHYuIzMTevXvRokULZGZm4vvvv3+mYZ2I5MSBnUiQjh07IiUlBVWrVkXbtm2xfv160ZGISCBFUbBo0SJ06dIFzZo1Q0pKCpo1ayY6FhGZAQ7spLng4GDREcxWgwYNkJycjCFDhmDEiBEYO3Ys7t27JzrWE7FP+bBT8bKzs9G/f38EBAQgICAA8fHxqFq1qkHvxT7lwj4J4D3sZAJ5eXmiI5i10qVL48svv8Qbb7wBPz8//Pjjj4iOjkadOnVER3sk9ikfdipWWloa3nnnHVy5cgXffPMN+vbta9T7sU+5sE8CAJ2iKIroEP9kzGNbiSzdkSNH0K9fP9y7dw+RkZHo0qWL6EhEpKFNmzZh1KhRqF+/Pr755hs0btxYdCQi0ogxMy5viSEyI2+88QZSU1PRvHlzdO3aFR9++CEKCgpExyIild29exfvv/8+Bg0ahJ49e+LIkSMc1onosTiwE5mZ6tWrY9euXZg7dy7mz5+Pzp0749q1a6JjEZFKzp49C1dXV6xduxarVq3C119/XaL11YnIenFgJ81lZmaKjmBxbGxsMGvWLOzduxdnzpxB8+bNsWfPHtGxALBPGbFT09m4cSNatGiB3NxcHDlyBKNHj4ZOp1P1GOxTLuyTAA7sZAI+Pj6iI1ist99+GydOnEDTpk3RpUsXfPjhh8Kfjso+5cNOtXfnzh28//77GDx4MHr16oWUlBS8/vrrmhyLfcqFfRLAgZ1MICgoSHQEi1ajRg3Ex8dj7ty5WLBgAdq3b49z584Jy8M+5cNOtXXq1Cm0atUK69atw6pVqxAZGYny5ctrdjz2KRf2SQAHdjIBrvZjPFtbW8yaNQvJycnIyMjA66+/joiICCFZ2Kd82Kk2CgsL8emnn8LFxQWlSpXC8ePHNbkF5t/Yp1zYJwEc2IksSps2bfDjjz+iT58+8PLygqenJ27fvi06FhH9y7Vr1+Du7o5JkybB19cXR44cgZOTk+hYRGShOLATWZjnn38eGzZsQEREBGJjY9GsWTMkJyeLjkVE/19cXBycnZ1x8uRJxMfHIyQkBKVLlxYdi4gsGAd20lx4eLjoCFIaOnQoTp48idq1a6N9+/bw9/fH3bt3NT8u+5QPO1VHdnY2Ro4cCb1eD1dXV5w6dQpdu3Y1eQ72KRf2SQAHdjKB1NRU0RGk9cILL+DAgQMIDg7GsmXLHj5BTUvsUz7s1Hh79uxB06ZNsXnzZqxevRrbt29HtWrVhGRhn3JhnwRwYCcTWL58uegIUrO1tYW/vz9SUlLg4OCANm3a4KOPPsL9+/c1OR77lA87NVxubi78/PzQuXNnvPTSS/jpp58watQozb9Y+iTsUy7skwAO7ETSeO2113DkyBHMmjULc+fORZs2bXDq1CnRsYiklZycjNdffx1r1qxBaGgodu/ejYYNG4qORUQS4sBOJBE7OzsEBQXh8OHDyM/PR8uWLfHBBx+Y5N52ImuRnZ0NPz8/tG/fHjVq1MDJkyfh5+cHGxv+L5WItMH/uhBJqFWrVkhNTcXs2bMRHByM5s2bcyUZIhXs2LEDTk5OWLduHZYuXYoDBw7glVdeER2LiCTHgZ00p9frRUewSvb29ggMDMSJEydQsWJFvPXWW/D19UV2drZR78s+5cNOn+769esYPHgwevXqBScnJ5w+fRoTJkyAra2t6GjFsE+5sE8COLCTCfj5+YmOYNWcnJyQnJyMzz77DOvXr4eTkxO2bdsGRVEMej/2KR92+niKomDdunVwdHREUlISIiIisGvXLrO+V519yoV9EsCBnUzAzc1NdASrZ2tri/Hjx+P06dNwdnbGO++8g169euG///3vM78X+5QPO320U6dOoX379vD29kb37t2RlpaGoUOHCl0BpiTYp1zYJwEc2ImsSsOGDbFjxw5s27YNp06dgpOTE+bNm4f8/HzR0YjMRnZ2NiZPnowWLVogMzMTe/bsQUREhLB11YmINB/Y8/PzERAQgNq1a6Ns2bJo06YNdu/erfVhiegxdDod+vy/9u48KorDjgP4dxXCjYoaFKiisqgoeICCRwRURKDiBUGNjYWo0Rqv5impV+LVJ+pTvFpjNdBKvCJGfQatRUE0csgpr4KAFyBEK3iALCjr9I88tyFQ5djdWcbv5715MbPsznfe9w38WGZnJk1CTk4OFi9ejLVr18LR0RH/+te/xI5GJCpBEHDkyBH06dMHX3/9NTZu3IisrCyMHj1a7GhE9I7T+MD++9//Htu3b8fvfvc77Ny5E23btoWvry9+/PFHTW+adMTJkyfFjkANMDExwaZNm5CVlQUrKyuMGzcOkyZNQkFBwRufxz6lh50CWVlZGDNmDKZPnw43Nzfk5OQgNDQU7733ntjRmox9Sgv7JEDDA3tKSgqOHj2KTZs2ISwsDLNnz8bFixfRvXt3LF++XJObJh1y+PBhsSPQGzg4OCAuLg5HjhxBeno6HBwcsHz58v97NRn2KT3vcqcPHjzA3LlzMWjQIJSUlCAmJgYnTpxAt27dxI7WbO9yn1LEPgnQ8MB+/Phx6OnpYe7cuap1BgYG+OSTT5CYmIj79+9rcvOkI44ePSp2BHoLmUyGoKAg5ObmYtWqVdizZw/kcjn2798PpVJZ52vZp/S8i53W1NRg8+bNkMvlOH78OMLDw5GdnQ0fHx+xo7XYu9inlLFPAjQ8sGdkZMDe3h6mpqZ11g8ZMgQAkJmZqcnNE1ETGRsbY82aNbh58ya8vLwwZ84cODs749y5c82+DCSRLhEEAd999x0cHBywYsUKBAcHo6CgAIsWLYK+vr7Y8YiIGqTRgb20tBRdu3att/71upKSEk1unoiaycbGBlFRUUhMTISZmRl8fHwwZswYXLt2TexoRM0WGxuLIUOG4MMPP0SfPn2QnZ2NHTt2wMLCQuxoRERvpNGBXaFQwMDAoN56Q0ND1eNEpLvc3NyQkJCA06dP4+HDhxg6dCgCAwORl5cndjSiRktLS4OXlxe8vLygr6+P+Ph4/PDDD+jbt6/Y0YiIGkWjA7uRkVGD13eurq5WPf7/+Pr6wt/fv84ybNiwep+WPn/+fIO37V2wYAEOHDhQZ116ejr8/f3x6NGjOuu//PJLhIWF1VlXWFgIf39/5Obm1lm/a9cuLFu2rM66qqoq+Pv748qVK3XWHz58GMHBwfWyBQUFvVP7ERwcLIn9AKTRR1P346uvvsKNGzeQlZWFiIgIxMTEoG/fvujevTsuXLjQavZDKn1oYj+Cg4MlsR9A3T5u3ryJoKAguLi4IDk5GV988QWuXr0Kd3f3VrUfv9SYPl4/r7Xvx2vv+n683m5r34/X3pX9OHz4sGp+9fPzQ5cuXRASElIvR6MJGjR27FjBwcGh3vrY2FhBJpMJZ86cqfdYWlqaAEBIS0vTZDTSokOHDokdgdQoMjJS2Lp1q9CpUydBX19f+PTTT4V79+6JHYtaQGrHaE5OjjBjxgxBJpMJNjY2woEDB4SXL1+KHUtrpNbnu459SkdLZlyNvsM+aNAg5OXloaKios765ORkAMDAgQM1uXnSEdOnTxc7AqnRrFmz8Pnnn+POnTvYsGEDoqOjYWdnh/nz56OoqEjseNQMUjlGc3JyMGPGDDg4OCAhIQF79uxBQUEBQkJCoKenJ3Y8rZFKn/Qz9kmAhk+JCQgIgFKpxL59+1TrampqEBERATc3N1hbW2ty80SkQaampli+fDnu3LmD9evX47vvvoOdnR3mzZv31psvEalTdnY2ZsyYgX79+uHy5cuqQX3+/PkNfo6KiKi10ejA/voDan/6058QGhqKffv2YfTo0SgsLMTmzZs1uWki0hJTU1OEhobizp07WLt2Lb7//nv07t0bQUFBSE9PFzseSZQgCLh06RJ8fX3h5OSEK1eu4C9/+QsHdSKSJI0O7ADwj3/8A0uWLMHBgwexePFiKJVKnDlzBiNHjtT0pklH/PrDItS6/b8+zczM8MUXX+Du3bvYs2cPUlNT4ezsjHHjxuHChQu8jrsOa03HqFKpxIkTJzBs2DB4eHiguLgYBw8exK1btzBv3jwO6mhdfdLbsU8CtDCwGxgYYPPmzSgpKYFCoUBSUhK8vLw0vVnSIfxrirS8rU8jIyPMmzcPN2/exJEjR/Do0SOMHTsWgwcPRkREhOoqUaQ7WsMxWlFRgd27d8PBwQFTp06FoaEhYmJikJWVhZkzZ/KmR7/QGvqkxmOfBGhhYCc6cuSI2BFIjRrbp56eHoKCgpCWlobz58/D2toaISEh+M1vfoOVK1eiuLhYw0mpsXT5GM3Pz8fixYthbW2NJUuWYMCAAUhKSkJ8fDx8fHwgk8nEjqhzdLlPajr2SQAHdtICY2NjsSOQGjW1T5lMBi8vL5w5cwZ5eXn46KOPsGvXLtja2mLatGmIj4/n6TIi07VjVKlU4uzZs/D19YW9vT0OHTqEhQsX4u7duzh27BhcXV3FjqjTdK1Pahn2SQAHdiLSIrlcjvDwcBQXF2P79u3IyMiAp6cnevfujS1btuDhw4diRyQRFRYWYu3atejZsyd8fX3x4MEDREZGoqioCBs3boSNjY3YEYmIRMGBnYi0ztzcHAsXLkRubi7i4uIwZMgQrF69GjY2NggMDMT58+ehVCrFjkla8OLFC0RHR8PHxwe2trbYunUrvL29kZSUhNTUVMyaNQuGhoZixyQiEhUHdtK4X99qmFo3dfYpk8ng4eGBb7/9FiUlJdiyZQtycnLg7e2Nbt26YdmyZcjMzOQpMxqm7WNUEARcvXoVCxYsgLW1NQICAvD06VPs378fpaWl2LdvH1xdXXl+ejPxe660sE8COLCTFnTr1k3sCKRGmurTwsICixcvRnZ2NpKSkjBlyhRERkZi0KBBcHR0xKZNm1BYWKiRbb/rtHWM5ubmYvXq1bCzs8OIESNw6tQpBAcHIzs7G1evXkVISAhMTU21kkXK+D1XWtgnAYBM0LG3rtLT0+Hs7Iy0tDQMHjxY7DhEJKKXL1/i/PnziIqKwsmTJ1FdXQ1XV1dMnToVU6dORc+ePcWOSG8gCAJu3LiBEydO4MSJE8jMzES7du0QEBCAjz76CKNGjULbtm3FjklEpBUtmXH1NJSJiKjF9PX14efnBz8/Pzx79gynT59GdHQ01qxZg+XLl2PgwIGYMmUKpkyZAgcHB55CoQMEQUBqaqpqSM/Ly4OpqSl++9vfYvXq1fD19eU56URETcSBnYhaBXNzc8ycORMzZ85EZWUlzp07h+joaGzevBlr1qxB9+7d4ePjA19fX4wePRomJiZiR35nPH78GLGxsTh79izOnTuH0tJSdOzYERMnTsS2bdswZswYDulERC3AgZ00Ljc3F3369BE7BqmJLvRpamqKgIAABAQEoLq6GnFxcTh79ixiYmKwd+9evPfee3B3d8f48ePh6ekJJycnnnrxBk3ttLa2FmlpaaohPSkpCUqlEv369cOMGTPg5+eHDz74AHp6/BEjBl04Rkl92CcBPIedtMDf3x+nT58WOwapia73mZ+fj5iYGMTExCAhIQHV1dVo37493N3d4eHhAQ8PDzg5OaFNG37m/rW3dfry5UukpaUhPj4ely5dwpUrV1BZWQkzMzOMHTsW48ePx/jx4/nhOB2h68coNQ37lI6WzLgc2EnjCgsL+YNcQlpTnzU1NUhOTkZ8fDzi4uKQmJiImpoamJubY+jQoXBzc4OrqytcXV3RuXNnseOK5tedFhcXIyUlBSkpKUhOTsa1a9fw/PlzmJqaYuTIkfDw8IC7uzucnZ2hr68vYnJqSGs6Runt2Kd0cGAnImqE6upqJCUl4erVq0hOTkZSUpLq7qo9e/aEs7MznJycVEv37t0l/UFWpVKJ27dv4/r168jOzkZWVhZSUlJQUlICALCxsVH9YuPu7o7BgwfzNBciombiVWKIiBrB0NBQdVoM8PMVTe7du4fk5GQkJycjIyMD27dvR3l5OYCfP+jq6OiIPn36wN7eHnK5HHK5HL169YKRkZGIe9I0FRUVKCgoUC35+fnIzs7Gv//9bygUCgBA586d4eTkhI8//hiurq4YOnQorKysRE5OREQAB3YieofJZDLY2trC1tYWQUFBAH4e4ktKSnD9+nVkZWXh+vXryMjIwLFjx1BRUaF6no2NDbp16wYbGxvVYm1tDWtra3Tq1AkdO3ZE+/btNfph19raWjx58gTl5eUoLS3F/fv36y23b9/GgwcPVM9p37495HI5+vfvj+nTp8PJyQmOjo6wtLTUWE4iImoZDuykcWFhYQgNDRU7BqmJ1PuUyWSqwdvHx0e1XhAEPHz4UPUOdUFBAYqKilBcXIyMjAwUFRWp3q3+5Wt16NABFhYWsLCwgLGxMYyNjWFkZKT6r4GBQZ3Tbl7/W6lUQqFQoLq6GgqFQvXvZ8+eoby8HOXl5Xj69Gm9/GZmZqr8crkc3t7esLOzg1wuh52dHSwsLOo9R+qdvmvYp7SwTwI4sJMWVFVViR2B1Ohd7VMmk8HS0hKWlpYYMWJEvccFQcCTJ09w//59lJWV1VseP34MhUKBqqoqVFZW4uHDh1AoFKipqVE9/5f/bdOmDYyMjFSLoaEh2rdvD3t7e9UvABYWFqpfCCwtLWFtbQ0zM7Mm79u72qlUsU9pYZ8E8EOnREREREQa15IZlxciJiIiIiLSYRzYiYiIiIh0GAd20rhHjx6JHYHUiH1KDzuVFvYpLeyTAA7spAUhISFiRyA1Yp/Sw06lhX1KC/skgAM7acFXX30ldgRSI/YpPexUWtintLBPAjiwkxbwaj/Swj6lh51KC/uUFvZJAAd2IiIiIiKdxoGdiIiIiEiHcWAnjTtw4IDYEUiN2Kf0sFNpYZ/Swj4J4MBOWpCeni52BFIj9ik97FRa2Ke0sE8CAJkgCILYIX6pJbdtJSIiIiLSRS2ZcfkOOxERERGRDuPATkRERESkwziwExERERHpMA7spHH+/v5iRyA1Yp/Sw06lhX1KC/skgAM7acFnn30mdgRSI/YpPexUWtintLBPAjiwkxaMGzdO7AikRuxTetiptLBPaWGfBHBgJyIiIiLSaRzYiYiIiIh0mEYH9sjISLRp06bB5eHDh5rcNOmQkydPih2B1Ih9Sg87lRb2KS3skwAtvcO+fv16REVF1VnatWunjU2TDggLCxM7AqkR+5Qediot7FNa2CcBgJ42NuLj49PkW7CSdHTu3FnsCKRG7FN62Km0sE9pYZ8EaOkddkEQUFFRAaVSqY3NERERERFJhlYGdk9PT7Rr1w4mJiaYOHEiCgoKtLFZIiIiIqJWT6OnxJiYmCA4OBienp4wNzdHamoqtm3bhuHDhyM9PR02Njaa3DwRERERUavX6IFdEATU1NQ06msNDQ0BAIGBgQgMDFSt9/f3h7e3N0aNGoWNGzfir3/9a73nKhQKAEBOTk5jo5GOS0lJQXp6utgxSE3Yp/SwU2lhn9LCPqXj9Wz7etZtEqGR4uLiBJlM1qjl5s2bb3ytYcOGCXZ2dg0+FhUVJQDgwoULFy5cuHDhwkVyS1RUVGPHb5VGv8Pet29fREZGNupru3Tp8sbHbWxskJeX1+Bj3t7eiIqKgq2tLYyMjBobj4iIiIhIZykUCty9exfe3t5Nfq5MEARBA5neyMXFBc+fP+dpL0REREREb6HRq8T85z//qbcuJiYG6enpGD9+vCY3TUREREQkCRp9h10ul2Pw4MFwdnZGu3btkJ6ejm+++QbW1ta4du0abwZARERERPQWGh3YV69ejR9++AF37txBVVUVrKys4Ofnhy+//JLDOhERERFRI4hyDjsRERERETWOVu50SkREREREzdMqBvbY2FiMHj0a7du3h7m5OVxcXHDs2DGxY1ELzZkzB23atMGECRPEjkLNdOHCBYSEhMDe3h4mJibo1asX5syZg59++knsaPQGNTU1CA0NhZWVFYyNjeHm5obY2FixY1EzXbt2DZ999hn69esHU1NTdO/eHUFBQcjPzxc7GqnJxo0b0aZNGzg6OoodhZopPT0d/v7+6NixI0xMTODo6Ihdu3Y1+vk6f0pMREQEZs+ejXHjxsHf3x9t27ZFbm4ubGxs8Mc//lHseNRMqampGD58OPT09DB27FicPn1a7EjUDC4uLnjy5AkCAwMhl8tx69Yt7N69G8bGxsjMzISlpaXYEakB06dPR3R0NJYuXQq5XI6IiAhcu3YNcXFxGDFihNjxqIkCAgKQmJiIwMBAODk5obS0FLt370ZlZSWSkpLQr18/sSNSCxQXF6N3795o06YNevTogevXr4sdiZro/PnzmDBhApydnREUFARTU1MUFBRAEARs2rSpUa+h0wP73bt34eDggE8//RTbt28XOw6piSAIGDFiBPr164fY2Fg4OjpyYG+lrly5gpEjR9ZZd/nyZbi7u2PlypVYv369SMno/0lJSYGbmxu2bt2qetOjpqYG/fv3x/vvv48ff/xR5ITUVImJiRgyZAj09P53L8SCggI4OjoiICAABw8eFDEdtdS0adNQVlaG2tpaPHr0CNnZ2WJHoiZ49uwZ7O3tMXLkSBw/frzZr6PTp8Ts3bsXgiBg3bp1AIDKykro8O8X1EgHDx7EjRs3sGHDBvbZyv16WAeADz74ABYWFsjNzRUhEb3N8ePHoaenh7lz56rWGRgY4JNPPkFiYiLu378vYjpqjmHDhtUZ1gHAzs4ODg4OPA5buYSEBERHRyM8PByCIEAmk4kdiZro0KFDePjwITZu3AgAeP78OV69etXk19HpgT02NhZ9+vTBmTNnYGNjA3Nzc3Tq1Alr1qzhoNdKVVRUIDQ0FCtWrODpEhJVWVmJiooKdOrUSewo1ICMjAzY29vD1NS0zvohQ4YAADIzM8WIRWomCAIePHjA47AVUyqVWLhwIebMmcPTmlqx2NhYmJubo6ioCL1794aZmRnatWuHP/zhD6ipqWn06+i9/UvEk5+fDz09PYSEhCA0NBQDBgxAdHQ0NmzYgNraWvz5z38WOyI10bp162BiYoKlS5eKHYU0JDw8HC9fvkRQUJDYUagBpaWl6Nq1a731r9eVlJRoOxJpwLfffouSkhJs2LBB7CjUTHv37kVhYSEuXrwodhRqgfz8fNTW1mLSpEmYPXs2wsLCEBcXh127duHJkyc4dOhQo15HawO7IAiN/k3C0NAQwP9OgQkLC8OyZcsAAJMnT0Z5eTl27NiBFStW1HuXiLSjOX3m5eVh586dOHLkCPT19TUZj5qhOZ3+WkJCAtauXYugoCB4eHioMR2pi0KhgIGBQb31rztVKBTajkRqlpubiwULFmD48OGYNWuW2HGoGcrKyrBmzRqsWbMGHTt2FDsOtUBlZSWqqqowf/58hIeHAwAmTZqEFy9e4Ouvv8a6detgZ2f31tfR2ikxly5dgrGxcaOWvLw8AICRkRFkMhmmT59e57WmTZsGhULBP92KqDl9Ll68GCNGjMDkyZNFTk8NaU6nv5Sbm4vJkyfDyckJ+/fvF2EPqDGMjIwa/MWsurpa9Ti1Xj/99BP8/PzQoUMHHD9+nOc8t1KrVq1Cp06dsHDhQrGjUAu9/p7661n29f8nJSU16nW09g573759ERkZ2aiv7dKlCwDAysoKt27dqneu8/vvvw8AePz4sVozUuM1tc+LFy/in//8J06cOIG7d++qHqutrUVVVRXu3bsHCwsLmJmZaSYwvVVzjtHXioqKMG7cOHTo0AExMTEwMTHRQEJSh65duzZ42ktpaSmAn7/vUuv09OlT+Pj44NmzZ7h8+XK945Rah/z8fPztb39DeHg4iouLVeurq6vx4sUL3Lt3D+bm5ujQoYOIKamxrKyscOPGjRbPslob2C0tLfHxxx836TkuLi4oKChAcXExevTooVr/+odN586d1ZqRGq+pfRYWFgIApkyZUu+xkpIS9OjRA+Hh4Vi0aJHaMlLTNOcYBX7+0+24cePw8uVLxMXF8cPEOm7QoEGIj49HRUVFnV+Qk5OTAQADBw4UKxq1QHV1NSZMmICCggLVBRuodbp//z5evXqFRYsWNfgzsUePHliyZAm2bdsmQjpqKhcXF8TGxqK4uBhyuVy1vqmzrE5fh/3UqVOYPHkyVqxYofrgzKtXr+Du7o7c3FyUlJTwXOhWoqioCBkZGXXWCYKAuXPnwtbWFitXrkT//v3Rs2dPkRJSczx//hyjR4/GzZs3ERcXh0GDBokdid7i9XXYt2zZgs8//xzA/67D3rlzZ1y9elXkhNRUSqUSU6ZMwblz53Dq1CmMHz9e7EjUAmVlZbhy5Uqd05kEQcCqVatQWVmJHTt2oFevXrxyTCuRmZmJwYMHY8aMGYiKilKtnzFjBqKjo3Hv3r1G/TVMpwd2APDy8sLFixcxZ84cODk54eTJk4iNjcW+ffswe/ZsseNRC9na2sLJyYk3TmqlJk2ahNOnTyMkJKTeh0zNzMwwceJEcYLRGwUFBeH777/H0qVL0atXL/z9739HamoqLly40OC19Um3LVmyBDt37sSECRMQGBhY7/GZM2eKkIrUzcPDA2VlZbxxUis0e/ZsfPPNN/jwww8xatQoxMfH4/jx43XekH4bnR/Ynz9/jlWrVuHo0aMoLy9Hnz59EBoaWu/kfWqdevTowTudtmI9evRAYWFhg/dFsLW1xe3bt0VIRW9TU1OD1atXIyoqCo8fP8aAAQOwfv16eHl5iR2NmsHT0xMJCQkNHocymQxKpVKEVKRunp6eKCsrw/Xr18WOQk30+lLkERERKCkpga2tLRYsWNCk04B1fmAnIiIiInqX6fSdTomIiIiI3nUc2ImIiIiIdBgHdiIiIiIiHcaBuF+LFAAAAFNJREFUnYiIiIhIh3FgJyIiIiLSYRzYiYiIiIh0GAd2IiIiIiIdxoGdiIiIiEiHcWAnIiIiItJhHNiJiIiIiHQYB3YiIiIiIh3GgZ2IiIiISIf9Fx6jYRCZ6zWcAAAAAElFTkSuQmCC" ], "text/plain": [ - "Figure(PyObject )" + "Figure(PyObject )" ] }, "metadata": {}, @@ -92,21 +123,21 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Works like a vessa of juna. We get derivative and function evaluated at same time. Beautiful solution." + "Works like a charm. We get derivative and function evaluated at same time. Beautiful solution." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Multidimensional functions" + "## Multidimensional functions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "From example:" + "From ForwardDiff example:" ] }, { @@ -193,6 +224,419 @@ "source": [ " That is, we need $n=2$ evaluations to get Jacobian. Each time." ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "# Comparisons" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "3x3 Array{Float64,2}:\n", + " 96.0 24.0 0.0\n", + " 24.0 96.0 0.0\n", + " 0.0 0.0 36.0" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Partial derivatives of bilinear Lagrange polynomials\n", + "dNdξ(ξ) = [[-(1-ξ[2])/4.0 -(1-ξ[1])/4.0],\n", + " [ (1-ξ[2])/4.0 -(1+ξ[1])/4.0],\n", + " [ (1+ξ[2])/4.0 (1+ξ[1])/4.0],\n", + " [-(1+ξ[2])/4.0 (1-ξ[1])/4.0]] \n", + "\n", + "E = 90\n", + "ν = 0.25\n", + "\n", + "μ = E/(2*(1+ν))\n", + "λ = E*ν/((1+ν)*(1-2*ν))\n", + "λ = 2*λ*μ/(λ + 2*μ)\n", + "\n", + "δ = eye(2)\n", + "C = zeros(2, 2, 2, 2)\n", + "C_sym = zeros(size(C))\n", + "\n", + "for i=1:2\n", + " for j=1:2\n", + " for k=1:2\n", + " for l=1:2\n", + " C[i,j,k,l] = λ*δ[i,j]*δ[k,l] + μ*(δ[i,k]*δ[j,l] + δ[i,l]*δ[j,k]) \n", + " end\n", + " end\n", + " end\n", + "end\n", + "\n", + "for i=1:2\n", + " for j=1:2\n", + " for k=1:2\n", + " for l=1:2\n", + " C_sym[i,k,j,l] = 1/4*(C[i,k,j,l] + C[i,k,l,j] + C[k,i,j,l] + C[k,i,l,j])\n", + " end\n", + " end\n", + " end\n", + "end\n", + "\n", + "D = E/(1-ν^2)*[[1 ν 0], [ν 1 0], [0 0 (1-ν)/2]]" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "1x4 Array{Int64,2}:\n", + " 1 1 1 1" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "a = 1/sqrt(3)\n", + "ipoints = [[-a -a], [a -a], [a a], [-a a]]\n", + "iweights = [1 1 1 1]" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "nonlinear_stiffness (generic function with 1 method)" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function nonlinear_stiffness(X, u)\n", + " n_integration_points = length(iweights)\n", + " n_problem_dimension = 2\n", + " n_solution_dimension = 2\n", + " n_shape_functions = 4\n", + " n_dofs = n_solution_dimension*n_shape_functions\n", + " K2_L = zeros((n_dofs, n_dofs))\n", + " K2_NL = zeros((n_dofs, n_dofs))\n", + " T = zeros((n_solution_dimension, n_shape_functions))\n", + " B_L = zeros(3, 8)\n", + " B_NL = zeros(4, 8)\n", + "\n", + " for m = 1:n_integration_points\n", + " ξ = ipoints[m, :]\n", + " w = iweights[m]\n", + "\n", + " Jᵀ = X*dNdξ(ξ)\n", + " ∇N = inv(Jᵀ)*dNdξ(ξ)'\n", + " ∇u = u*∇N'\n", + " F = I + ∇u # Deformation gradient\n", + "\n", + " detJ = det(Jᵀ)\n", + " #detF = det(F)\n", + " E = 1/2*(∇u + ∇u' + ∇u'*∇u) # Green-Lagrange strain tensor\n", + " P = λ*trace(E)*eye(2) + 2*μ*E # PK2 stress tensor\n", + " S = F*P\n", + "\n", + " B_L[:,:] = 0\n", + " for i = 0:3\n", + " B_L[1, i*2+1] = F[1, 1]*∇N[1, i+1]\n", + " B_L[1, i*2+2] = F[2, 1]*∇N[1, i+1]\n", + " B_L[2, i*2+1] = F[1, 2]*∇N[2, i+1]\n", + " B_L[2, i*2+2] = F[2, 2]*∇N[2, i+1]\n", + " B_L[3, i*2+1] = F[1, 1]*∇N[2, i+1] + F[1, 2]*∇N[1, i+1]\n", + " B_L[3, i*2+2] = F[2, 1]*∇N[2, i+1] + F[2, 2]*∇N[1, i+1]\n", + " end\n", + " K2_L += w*B_L'*D*B_L*detJ\n", + "\n", + " τ = zeros(4, 4)\n", + " τ[1:2,1:2] = S\n", + " τ[3:4,3:4] = S\n", + " B_NL[:,:] = 0.0\n", + " for i=0:3\n", + " B_NL[1, 2*i+1] = ∇N[1, i+1]\n", + " B_NL[2, 2*i+1] = ∇N[2, i+1]\n", + " B_NL[3, 2*i+2] = ∇N[1, i+1]\n", + " B_NL[4, 2*i+2] = ∇N[2, i+1]\n", + " end\n", + " K2_NL += w*B_NL'*τ*B_NL*detJ\n", + "\n", + " # Nodal forces\n", + " T += w*S*∇N*detJ\n", + "\n", + " end\n", + "\n", + " return T', K2_L+K2_NL\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "2x4 Array{Int64,2}:\n", + " 0 0 1 0\n", + " 0 0 0 0" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "X = [0 0; 1 0; 1 1; 0 1]'\n", + "u = [0 0; 0 0; 1 0; 0 0]'" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "elapsed time: 1.159613267 seconds (79896004 bytes allocated, 4.41% gc time)\n" + ] + }, + { + "data": { + "text/plain": [ + "(\n", + "4x2 Array{Float64,2}:\n", + " -60.6667 -24.0\n", + " -1.33333 -8.0\n", + " 133.333 35.0\n", + " -71.3333 -3.0,\n", + "\n", + "8x8 Array{Float64,2}:\n", + " 136.333 34.0 -29.3333 -2.0 … -34.0 23.0 2.0 \n", + " 34.0 81.6667 6.0 3.33333 -71.6667 -7.0 -13.3333\n", + " -42.3333 6.0 83.3333 10.0 -28.0 -45.0 12.0 \n", + " -2.0 -9.66667 10.0 52.6667 -11.3333 11.0 -31.6667\n", + " -127.333 -33.0 -20.6667 -19.0 55.0 -153.333 -3.0 \n", + " -34.0 -69.0 -28.0 -36.0 … 162.667 7.0 -57.6667\n", + " 33.3333 -7.0 -33.3333 11.0 7.0 175.333 -11.0 \n", + " 2.0 -3.0 12.0 -20.0 -79.6667 -11.0 102.667 )" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "@time nonlinear_stiffness(X, u)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "calculate_internal_energy! (generic function with 2 methods)" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function gauss2d(J, Wint)\n", + " a = 1/sqrt(3)\n", + " ipoints = [[-a -a], [a -a], [a a], [-a a]]\n", + " iweights = [1 1 1 1]\n", + " Wint[:,:] = 0\n", + "\n", + " for m = 1:length(iweights)\n", + " w = iweights[m]\n", + " ξ = ipoints[m, :]\n", + " Wint[:,:] += w*J(ξ)\n", + " end\n", + "end\n", + "\n", + "function calculate_internal_energy!(X, u, Wint, dNdξ, λ, μ, dim=2)\n", + " u = reshape(u, 2, 4)\n", + " X = reshape(X, 2, 4)\n", + " Wint = reshape(Wint, 2, 4)\n", + " I = eye(dim)\n", + " \n", + " function J(ξ)\n", + " Jᵀ = X*dNdξ(ξ)\n", + " ∇N = inv(Jᵀ)*dNdξ(ξ)'\n", + " ∇u = u*∇N'\n", + " F = I + ∇u # Deformation gradient\n", + " E = 1/2*(∇u' + ∇u + ∇u'*∇u) # Green-Lagrange strain tensor\n", + " P = λ*trace(E)*I + 2*μ*E # PK1 stress tensor\n", + " S = F*P # PK2 stress tensor\n", + " return S*∇N*det(Jᵀ)\n", + " end\n", + "\n", + " gauss2d(J, Wint)\n", + " u = reshape(u, 8)\n", + " X = reshape(X, 8)\n", + " Wint = reshape(Wint, 8)\n", + "\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "8-element Array{Float64,1}:\n", + " -60.6667 \n", + " -24.0 \n", + " -1.33333\n", + " -8.0 \n", + " 133.333 \n", + " 35.0 \n", + " -71.3333 \n", + " -3.0 " + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function Winte!(u, Wint)\n", + " calculate_internal_energy!(X, u, Wint, dNdξ, λ, μ)\n", + "end\n", + "Wint = zeros(8)\n", + "Winte!(reshape(u, 8), Wint)\n", + "Wint" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "g! (generic function with 1 method)" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using ForwardDiff\n", + "Kt = forwarddiff_jacobian!(Winte!, Float64, fadtype=:dual, n=8, m=8)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "ename": "LoadError", + "evalue": "`dual_fad` has no method matching dual_fad(::Function, ::Array{Int64,1}, ::Array{Float64,2}, ::Array{Dual{Float64},1}, ::Array{Dual{Float64},1})\nwhile loading In[14], in expression starting on line 2", + "output_type": "error", + "traceback": [ + "`dual_fad` has no method matching dual_fad(::Function, ::Array{Int64,1}, ::Array{Float64,2}, ::Array{Dual{Float64},1}, ::Array{Dual{Float64},1})\nwhile loading In[14], in expression starting on line 2", + "", + " in g! at /home/jukka/.julia/v0.3/ForwardDiff/src/dual_fad/multivariate_range.jl:18" + ] + } + ], + "source": [ + "J = zeros(8, 8)\n", + "Kt(reshape(u, 8), J)\n", + "Kt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This ForwardDiff package is just not working as expected. We have well defined function but for some reason it's just not working. Bad code." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "elapsed time: 1.198824455 seconds (532828332 bytes allocated, 27.60% gc time)\n" + ] + } + ], + "source": [ + "function calc_1()\n", + " for i=1:10000\n", + " K, T = nonlinear_stiffness(X, u)\n", + " end\n", + "end\n", + "@time calc_1()" + ] } ], "metadata": {