tp07-EF2D-sujet.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TP 7: Éléments finis en 2D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dans ce TP, on s'intéresse à la résolution approchée du problème aux limites suivant, avec conditions aux limites de Neumann homogène :\n",
    "\n",
    "$$\n",
    "(PN)\\left\\{\\begin{aligned}\n",
    "&-\\Delta u+u=f,\\ \\textrm{dans }\\Omega,\\\\\n",
    "&\\frac{\\partial u}{\\partial n}=0, \\textrm{dans }\\partial\\Omega\n",
    "\\end{aligned}\n",
    "\\right.\n",
    "$$\n",
    "où $f:\\overline{\\Omega}\\longrightarrow\\mathbb{R}$ est une fonction de classe $\\mathcal{C}^2$ donnée et où $\\Omega$ est un ouvert de $\\mathbb{R}^2$ donné. \n",
    "\n",
    "**Q1)** Écrire la formulation variationnelle du problème $(P)$ sous la forme \n",
    "\\begin{equation*}\n",
    "(PNv)\\ \\ \\ \n",
    "\\begin{cases}\n",
    "u\\in V  \\,\\textrm{tel que}\\\\\n",
    "a(v,u)=L(v),\\ \\,\\forall\\ v\\in V ,\n",
    "\\end{cases}\n",
    "\\end{equation*}\n",
    "où $V=H^1(\\Omega)$, $a$ est une forme bilinéaire dans $V$  et $L$ une forme linéaire de $V,$ que l'on explicitera. Montrer que le problème $(PNv)$ admet une unique solution.\n",
    "\n",
    "Soit $T_h$ une triangulation du domaine $\\Omega$ et $V_h$ l'approximation de $H^1(\\Omega)$ par des éléments finis $P^1$ associés à la triangulation $T_h.$ On note $\\{T_N\\}_{N=1,\\dots,Ntri}\\ $ les triangles de $T_h,$ $\\{S_I\\}_{I=1,\\dots,NSom}\\ $ les sommets de la triangulation et $\\{\\Phi_I\\}_{I=1,\\dots,NSom}\\ $ les fonctions de la base de $V_h$ définies par $\\Phi_I(S_J)=\\delta_{IJ},\\ I,J=1,\\dots,NSom.$ \n",
    "\n",
    "Le problème discret consiste alors à chercher $u_h\\in V_h$ tel que \n",
    "$$\n",
    "(PNv_h)\\ \\ \\ \\ \\ \\ \\ \\ \\ a(\\Phi_I,u_h)=L(\\Phi_I),\\ \\,\\forall\\ I=1,\\dots,NSom. \n",
    "$$\n",
    "\n",
    "**Q2)** Soit $u_h=\\displaystyle{\\sum_{I=1}^{{NSom}}u_I\\Phi_I}.$ Montrer que $u_h$ est solution de $(PNv_h)$ si et seulement si le vecteur $U=(u_1,\\dots,u_{NSom})^T$ est solution d'un système linéaire \n",
    "$$\n",
    "KU+MU=F,\n",
    "$$\n",
    "où $K$ et $M$ sont les matrices de $\\mathcal{M}^{NSom\\times NSom}(\\mathbb{R})$ définies par $K_{I,J}=\\int_{\\Omega}\\nabla\\Phi_I\\cdot\\nabla\\Phi_J,$ $M_{I,J}=\\int_{\\Omega}\\Phi_I\\Phi_J,$ et où $F$ est le vecteur de $\\mathbb{R}^{NSom}$ défini par $F_I=\\int_{\\Omega}f\\Phi_I.$ \n",
    "\n",
    "** Maillage de $\\Omega.$ **\n",
    "\n",
    "Dans un premier temps on considère $\\Omega$ le carré $[0,2]\\times[0,2].$ \n",
    "\n",
    "La triangulation $T_h$ de $\\Omega$ est représentée par deux matrices. La première, que l'on appellera dans le programme python $TabTri,$ contient la liste des triangles. Il s'agit d'une matrice de taille $NTri\\times3$ que, dans chaque ligne $N$, contient l'indice $I$ de chaque sommet $S_I$ appartenant au triangle $T_N.$ La deuxième, que l'on appellera $TabSom,$ contient les coordonnées de chaque sommet de la triangulation. C'est une matrice de taille $NSom\\times2$ que, dans chaque ligne $I$ contient les coordonnées du sommet $S_I.$ \n",
    "Le code suivant crée le maillage et ces matrices.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQgAAAEICAYAAACj9mr/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmYFNX1vz+HzRgXXHALoggqKAgIuOK+gkrAJQoKalQM\ni9EgqDAsI0a/CCoosoOGmCgRFRhlc4tBjSuKIiOLgIAaZBl2oyzD+f1x+mTatuveW1tPz/zqfZ55\nemb6dlXdrqq3bt0691xiZiQkJCRko0p5b0BCQkL+kggiISHBk0QQCQkJniSCSEhI8CQRREJCgieJ\nIBISEjxJBJHgBBGNJaIBqd/PI6Jv095bSUQXld/WJcRFtfLegISKATN3Le9tSMg9SQsiISHBk0QQ\nlZxU8/8eIlpARD8Q0VNEdBgRzSaibUT0BhEdmCr7AhF9T0RbiOhtImqUtpxJRPSgw/pOJaL3iWgz\nEa0hopFEVCPt/UuIaElqHaOJaC4R3Zb2/i1EtIiINhHRq0R0dNTfSYI7iSD+/+BqABcDOB5AWwCz\nARQAOARyDNyZKjcbwHEADgXwKYBnA6yrFEBPALUAnAHgQgDdAYCIagF4EUBfAAcDWALgTP0gEbVL\nbddVqW17B8DkANuQEBGJIP7/4ElmXsvM30FOug+ZeT4z/wRgGoCTAYCZn2bmbcy8A8D9AJoSUU0/\nK2LmT5j5A2bezcwrAYwDcG7q7csAFDPzVGbeDWAEgO/TPt4VwGBmXpR6//8ANEtaEeVHIoj/P1ib\n9vuPWf7el4iqEtHDRLSciLYCWJl6v5afFRHR8UQ0I3WrshVykusyfgPgGy3LMlLw27SPHw3gidTt\nyWYAGwEQgNp+tiEhOhJBJCjXA2gH4CIANQHUTf2ffC5nDIDFAI5j5v0htwy6jDUAjtSCRETpf0Pk\n8QdmPiDtZ29mfs9vZRKiIRFEgrIfgB0ASgD8GnLlD7qcrQC2E1FDAN3S3psJ4CQiak9E1QD0AHB4\n2vtjAfTVzlEiqklEvwu4HQkRkAgiQXkGwCoA3wH4EsAHAZfTG9Ia2QZgAoDn9Q1m3gDgdwCGQkR0\nIoB5EDGBmacBGALgH6nbk4UA2gTcjoQIoCRhTEJ5QURVIH0QNzDzW+W9PQm/JGlBJOQUIrqUiA4g\nor1Q1j8RtLWSEDOJIBJyzRkAlgPYAInJaM/MP5bvJiV4kdxiJCQkeJK0IBISEjzJy9GctWrV4rp1\n65b3ZiQkVFo++eSTDcx8iK1cXgqibt26mDdvXnlvRkJCpYWIVrmUS24xEhISPEkEkZCQ4EkiiISE\nBE8SQSQkJHiSCCIhIcETqyCIqA4RvUVEXxJRMRHdlaUMEdEIIlqWSm3WPO291qkUY8uIqE/UFUhI\nSIgPlxbEbgC9mPlEAKcD6EFEJ2aUaQNJVXYcgNshOQFARFUBjEq9fyKAjlk+m5CQkKdY4yCYeQ0k\n0QeYeRsRLYJk+PkyrVg7AM+kMgR9kBqMcwQk6cgyZl4BAET0j1TZ9M8GYudOoJpDFMeePQCR/Lgs\ns0YNe7kffwT23ttejll+qjho2HWZu3ZJvV3qs2dPtOv2Ux/X7/Knn4C99op2/7jWmxnYsQP41a+i\nW6af+sSxf4jclumKr0ApIqoLyV/4YcZbtZGWSgwyhLe2x/9P81j27ZDWB4466ijrtuy1F3DddcCJ\nlvZIYSHwu98BjRuby02eDCxeDAwaZC63Zg0wdizQsiXQtq257BtvAB98APTvby43bx7wyitu21lY\nKK+27dSytnJr1wKjRwPNmgFXXmku++67wOuv25c5fz4wfbq93M6dwEMPye/3328+qZilTI8ewKGH\nmpfrUm8tBwAFBXI8hV3mrl3Ag6m837b6uC6zpAQYMQJo2BDo2NFcdvVqEcmTT5rL+cFZEES0L4CX\nAPyJmbdGtwkCM48HMB4AWrZsaR1B1r+/HGDXXGOWRO3awCWXAHXqeJeZPFkOkuXLgYEDvcutWQP8\n3/8B99wD1KsnyzZJ4sorgaVLgauv9i4zbx6wezfQtKm8XnWVWRJ160q9N28Gevf2LgfIFcpUn7Vr\ngT//WZbToAFw8MFmSSxeDHz2GdChg3eZ+fNl+046ybzunTuBe++V77NWLWDjRvnbdFJt3w78979A\n165mSdjqDQCPPAJMnCj7dP16KW+ShG2Zu3bJ9g8eDBx0ELBhA9C3r7k+tmWWlMj7PXsCJ58MVK9u\n/u5Xr5aLUpQ4NUaIqDpEDs8y89QsRb4DkH4KHpn6n9f/Q1O1qhxcf/0r8GWIG5bJk+XE7NzZXE7l\nMGQI8Otfy0H63Xdy5Q/KvHnAyy/LVaRqVbmaPvsssHCh+XMXXAA0bw48+mjwdaschgwB9tkHuO02\nOSCnTQu+zPnzgZdeAh54wNzMVTn86U8ivHPOAc44Axg6VFoKXuyzj5QZNAhYty74dj7yCHDKKcB5\n5wFHHw3cfbdsz44dwZancrjjDqB+faBVK+Dcc0UWQQdLqxwGDwb231+Oz127gH/8I9jyguLyFIMA\nPAVgETMP8yj2MoAbU08zTgewJdV38TGA44jomNTkKR1SZSMhrCSCykEJI4l0OehVpkqV3EgiUw5K\nGEkElYOSK0mky0EJI4lMOShhJJEpB6U8JOHSgmgFoDOAC4jos9TPZUTUlYh0vsZZAFYAWAbJQ9gd\nAFJzG9wB4FUAiwBMYebiKCsQVBJh5aAEkUQ2OShxS8JLDkoQSYSVgxK3JLLJQQkiCS85KEEk4SUH\nJdeSsAqCmd9lZmLmJszcLPUzi5nHMvPYVBlm5h7MXJ+ZT2LmeWmfn8XMx6feeyiOSviVRFRyUPxI\nwiQHJS5J2OSg+JFEVHJQ4pKESQ6KH0nY5KD4kYRNDkouJVFpIildJRG1HBQXSbjIQYlaEq5yUFwk\nEbUclKgl4SIHxUUSrnJQXCThKgclV5KoNIIA7JKISw6KSRJ+5KBEJQm/clBMkohLDkpUkvAjB8Uk\nCb9yUEyS8CsHJReSqFSCALwlEbcclGySCCIHJawkgspBySaJuOWghJVEEDko2SQRVA5KNkkElYMS\ntyTyMqNUWFQSBQUihZIS4Igj4peD0rWrBFN98QWwbRvw+efB5KCoJPr1A0pLgbPO8i57wQXy+uij\n4eWg3HabxAxMmybxFwsXxi8H5Zxz5NVVEvfeK+XCyEFJl0T16uHkoLRqJa+DB4eXg9K5M/C3v4kk\navmaSdVOXma1btmyJdtSzhHJCWMKty4tLYts06g5L775Bnj6aQmCsslh0CD78rQcIAeBTQ4uy9yz\nR074o44Cfv97c9nnngO++kqCoExyYJaTvXt34BBLhkKtz4ABdjkMHgzcdJNdDrpMW93ffht46y17\nuR9+EDk2aGAOKgKA114D3n/fvszly4G//x249lrghBPMZQcNkjLXXmsu9/LL0hLr2dMuB7/Hm8sp\nTUSfMHNLW7kK24Lo2FHi3h9+2CyJunXtkZRr1sgVulkzaWn07Glf//33m9+fN09aDlu3Ai1aAL/9\nbfhlTpkCnH++XM3r1fMut3atRAe2aydXUi+Y5fs75RT5+447zFegDh3skZQ7d4pkb7oJmDTJWB1M\nnSr1OeggieBs2tRc/v777d/R0KHyfZ92mkjPFHHZuTPw3nvmlqXeVpx9NnDYYRIdaYq4bNRIInsb\nNfIuU1Ii+6d6dYleLSiwX0Bs9S4uBpYsKWtBRkWF7YM4/njg9tuBPn3kNiIoKoehQ+WEatwYGD48\n3Lal9zl07y59Ei+HDA+bMkXCjG++2Vxu7VppEQwdam85PPywNHkvu0x+LyyUEOGgqBx69rS3HKZO\nlfDqLl2kVTRligg1DEOHiuyuuCKaiMv0PocLLgB69ZL6BY24BH5+W9GmjfRJ/N//BY+4BEQOzzwj\ny6laNfhyslFhBQGEl0S6HPS24uKLw0kiW4dkt27hJBGnHPQef7/9wkkiiBxuu03+rlIlvCRUDuef\nL3+HjbjM1iF59NHhJJGtz+Gss8JJIk45ABVcEEBwSWSTgxJUEqanFUElkQs5KEElEUYOShhJZMpB\nCSoJ09OKoJIwdUgGlUTccgAqgSAA/5IwyUHxKwmXR5l+JZFLOSh+JRGFHJQgkvCSg+JXEi6PMv1K\nwuVphV9J5EIOQCURBOAuCRc5KK6S8BPn4CqJ8pCD4iqJKOWg+JGETQ6KqyT8xDm4SsLPo0xXSeRK\nDkAlEgRgl4QfOSg2SQQJgrJJojzloNgkEYccFBdJuMpBsUkiSBCUTRJB4hxsksilHIBKJgjAWxJB\n5KB4SSJMhKSXJPJBDoqXJOKUg2KShF85KF6SCBMh6SWJMEFQXpLItRyASigI4JeSCCMHJVMSYeSg\nZEoin+SgZEoiF3JQskkiqByUTEmEDZ8GfimJKCIkMyVRHnIAKnAk5Tnn2PsGli4Frr9eAlLefNMu\nh7Zt7UO2X39dgmUuu8wuh7lzgY8+kgPHxJgxwFNPyUFqk0OXLjLG5LXX7HI45RRg2DC7HFzqvW1b\n2Un54ot2ObRsKSHnNjn8+c8SyHZa1kylwp49Er356qsS+GWTg0t9fvhBTkAAeP55uxxOOQX4+GNz\nmVWrytIL/vOfdjm4bOe770q4+oUX2uWwbBnwySciQBuVPpJyxQrgP/8xl9m2TV537ZKrny0L9po1\n9mX+9JO8/vijlDexbJmczLZlatN0xw572fnz5YDZskV+ssEsVy/dXtsyXeq9c2fZ71u32ssDbvVZ\nuBA49lhzpCsg+SiB6Oqza1fZ79u22cszS/LYa67xLlO9etnvJSVl2xxmO3/8UV43bwa+/958QVq9\nWtYbJVZBENHTAK4AsI6Zf5FOlYjuAXBD2vJOAHAIM28kopUAtgEoBbDbxViu3HabOWGsDrwaMUJC\nX2fPtodlFxaal/nJJ3JrMXCgNAGLi+Vk9aJuXTG/aZlTpkiZ+vWBww+X/5nK//Wvsm6vUGtmGZw1\nZIgkMN2yBTjmGHMiXFu9deDV7NkyZmXePGk9mcKyCwuBI4+UWzxTItzt22UEaoMG3mUeeUTG1Pzz\nn3LCHHWUOSzbVp9du4D77pMWyYoVwKefymdMYdmFhbLeHTvMkpgyRVp2r70m+8AUlm3bzuJiaUH0\n7w+0bi0tUlMi3KZNyydp7SQArb3eZOZHNNMUgL4A5jLzxrQi56fej0wONjJHZUYRlv3JJ9JXoOnM\nL75Y4u0ffzz4dmb2OYRNhKtyOPNMua2oUkW+h7//3T5U3IvMUZl+4iSiSISbPioziohLlUOPHiJl\nP3ESt9wiracXXzSXO/po6aO5777gYdmZfQ5RJMINgkvKubcBbLSVS9ERwORQWxQSryHbYSSRKQcl\njCS8OiSDSiJTDkoYSXgN2c6VJLIN2Q4jiUw5KPkmCa8OyfKQRGRPMYjo15CWxktp/2YAbxDRJ6mJ\ncUyfv52I5hHRvPXr1wfaBls+hyCS8JKDEkQStqcVfiXhJQcliCRs+RziloQpn0MQSXjJQckXSdie\nVuRaElE+5mwL4N8ZtxdnpW492kDm9PTsT2fm8czckplbHmJLTJAF12QvfiRhk4PiRxKujzJdJWGT\ng+JHEq7JXuKShEuyFz+SsMlBKW9JuD7KzKUkohREB2TcXjDzd6nXdQCmATg1wvX9D7+ZoFwk4SoH\nxUUSrnJQbJJwlYPiIgm/maCiloSfTFAuknCVg1JekvAb55ArSUQiCCKqCeBcAEVp/9uHiPbT3wFc\nAiBgV5k3QdPEmSThVw6KSRJ+5aB4ScKvHBSTJIKmiYtKEkHSxJkk4VcOSq4lETQIKheScJlZazKA\n9wE0IKJviejWjElzAOBKAK8x8w9p/zsMwLtE9DmAjwDMZOY5UW582ByS2SQRVA5KNkkElYOSKYmg\nclCySSJsDsmwkgiTQzKbJILKQcmVJMJGSMYtCWscBDNb5hQGmHkS5HFo+v9WALAkEQtOVAlm0yVR\nXFw2i3TQ8GlAJAGIJLZvl8xBQeWgaCLcN9+UZ+OdOgUPnwbKJFFQIMFcUSSYVUn06WP//jQR7muv\nSZ06dQqXYFYlMWCASCKMHJTMRLgmbrlFcpq++qp5Qt50SXz9tcR1hA2fTk+Eq+kDo6LCRlKOHy8h\nzKawUk3IesYZEs5rYto0CZpp3tw+JfugQXJSmWYVB+SAASSJqS2noMt09cxlV6nf/EYCh8IuE5CW\nSfPmwAEHmMtNny4n3+LF5nLbtkmWZZcMyy+8IK9HHgn861/e5XRf7twJ1KjhXa5qVWDmTIkqtIU6\na9La5cvN5X78EfjLX8yBVMqUKbKPBgwwl3v/fQnDb9hQxGbCdT+OGyc/t97qVt6FCiuIwkLzSccs\nZr7lFilnCuWdN0/i/b/6Spp+gwbZs2Xv2CGZi02SuPpqYORIoGbN8Ilw9bbittvkRG7ePHwiXL2t\n6NxZhHPNNeaIS5ektevXyzp79bJL8ZFHpPW2334SBWiKuATkivvf/8pyTZIAZB/ZEuG6JK3VgVc9\ne9rrM2yY1OfOO81Ja4uLJWS8YUORzwMPmCMuAfu6Z86UVsnxx5vL+aVSjuZUOZx/vtynm0gfldmg\ngdsj0KpVpTk3aZJ9LtBLLpGDJUwiXJVDq1aSXTlsjkvg530O9eqFj7gEyuQwZAiw777msul9Dl26\nSL+F7RFovXpy8t1778/HhmQSRcQl4G9U5rBhQJMmwEUXmcsVF0vr6qGHJAQ+ikS4M2fKQLHu3YMv\nw4tKJ4igctB7Ztc4iVxJIlMOShhJZOuQDBuWHVQOSr5JIm45aJ9D2ES4ccoBqGSCCCsHJV8k4SUH\nJYgkTE8rgkoirByUfJFEruSgBJVE3HIAKpEgopKDUt6SsMlB8SMJl0eZfiURlRyU8pZEruWg+JVE\nLuQAVBJBRC0Hpbwk4SoHxUUSfuIcXCURtRyU8pJEeclBcZVEruQAVAJBxCUHJdeS8CsHxSSJIEFQ\nNknEJQcl15IobzkoNknkUg5ABRdE3HJQciWJoHJQskkiTISklyTiloOSK0nkixwUL0nkWg5ABRZE\nruSgxC2JsHJQ0iURNnwa+KUkciUHJW5J5JsclExJlIccgAqctJYI+P3v7XLo0kVex483y+HHH+UA\nmzDBvLylS+XAHzPGHkzVNTVaxbbM116TiMIbb7TLoUsXmZy2XTtzudGjJX9lv352OXTpIuXTcypm\nsmcP8Ic/yO/Dh9vl0KULcMMNdjl06SJBX926mctNmCCRh7bvcsUKOdnHjRMZeJFen8ces8uhSxfg\nuuvsctDjbexYuxy6dLHXZ9UqSbfXoAHQu7e57JdfSjCZy5D6Sp+0tkYN4JBDJFefiZtukp/jjvMu\ns2NHWebpPXskY7UXrVtL6+Bf/7LnuHz+ecmEbNpGZjmRAQlNttXnlluAt9+Wg+bgg83b2bZtmaS8\neCmV3uett6RVYJLEiy9KtGWDBsBJJ5mX26qVfD8mFiyQ1/XrJQu2KTS7dWuJ4LR9P0Wp8cSmfQhI\ny+Hoo+UEPP54oFkzc/lWrYBHHzWXKS6W10MPlZaTLTT70kvt9ZkxQ14bNZLWsini8vjjJTI0Siqs\nIPr2lVTp06YBf/yjd7mzz5asyUcemf19lcPjj0vzcedOOXBNB1j37pJOf+RIsyROOEGuIl7r1nkr\nrrtOZNeqlbQkTGHZ554r6fHHjZOmskkS7dp5rxsQOeyzj4Std+oEPPmkWRLbtwPPPivJU4880iyJ\niy4yr3vBAuCdd+TWpUEDuZLaEuHalllUJN+3Jsz1oqREWpTTpsl2vPeelDdJwrbu4mIZGzN5stzu\nPPWUPRHu6aeblzljhrRECwultTxsmIw98pLEnj2SST1KKmwfBCCjJI89Vg7sIKgceveWjMWAnHwr\nVgCzZpk/GzYRbrZJbVznAq1ZU5rRAwYET3P+0kvApk1l81Yce6yI7777fp4SPhPtk3j2WeCLL4Kt\ne8ECOZEeekiWF3RW8XSKimSEr63FlNnnUKWKtMamTJFxJkHIHLIddFbxdGbMkNZNjx7y99FHA3ff\nLX0sYcKy/VKhBQEEl0Q2OShxS8I041UuJJEpByUXksiUgxJGEkHloISRhFc+hzCSyJSDUh6SqPCC\nAPxLwiQHJS5JuEyHF6ckvOSgxCkJLzkoQSQRVg5KEEnYkr0EkYSXHJRcS8Ilo9TTRLSOiLLG1BHR\neUS0hYg+S/0MTHuvNREtIaJlRNQnyg3PxFUSLnJQopaEn7ky45CETQ5KHJKwyUHxI4mo5KD4kYRr\nJig/krDJQcmlJEJPnJPiHZ08h5kfAAAiqgpgFCSj9YkAOhKRJcVKOGyS8CMHJSpJBJlIN0pJuMpB\niVISrnJQXCQRtRwUF0n4TRPnIglXOSi5kkTUE+ekcyqAZcy8gpl3AvgHAMvT+/B4SSKIHJSwkggz\ny3YUkvArByUKSfiVg2KSRFxyUEySCJpD0iQJv3JQciGJqPogziSiBUQ0m4g0l05tAN+klfk29b+s\nRDFxjpIpiTByUIJKIowclDCSCCoHJYwkgspBySaJuOWgZJNE2ASz2SQRVA5K7JJgZusPgLoAFnq8\ntz+AfVO/Xwbgq9Tv1wCYmFauM4CRLutr0aIF2ygstBbhWbOYjz2WuX175lWr7OVdlvnkk8wOm8dL\nljBfcglzo0bMc+dGs+7XXmNu0oR5+XJzuc2bmbt1Yz70UOYJE6JZ91dfMbdpwzxpkrlcaSnzffcx\nH3MMc58+8reJ555jXrzYXGbrVubu3ZkPO4x5zBj7tnbvLj9btpjLLVvG/Mwz5jKlpcx9+zLXrct8\n773Mu3eby0+ezPzFF+Yy27fL9tWuzTxypLkss9v+WblSjnOX5TEzA5jHDudi6EApZt6a9vssIhpN\nRLUAfAcgPRPkkan/RYJLIs8dOyRwZNkyc25CP8tcu1ZS4xcWmkO3mSWEGpDXKBLMahDX+edL4IyJ\nV16Rq9SCBdEkzAVkdu/ZsyWc18TixVJmyRLJt2hb98kn2/NrvvqqfPcLF9rrM3q0pHOzTdDmmrR2\n6VJg5Upg0SL3BLMDB5qPj733ljEzxcXR7Z/p0+UnaGskG6EFQUSHA1jLzExEp0JuW0oAbAZwHBEd\nAxFDBwDXh12fYktau2OHNLsGDgRq15a/TRGXimmZa9fKAdK7t3mn6cCrUaOALVukuWsLy7atWwde\nDRokUY/16nmXnTpVvp/iYvncnXcCBx0UfN1A2fiTAw4wJ61dsEC+68aNJQrQlgi3QQMZi9GggXeZ\noiJJqrtkiSR7veMOe8bsI4+UKFNTIlyXpLXFxTJOp18/2a+2RLiNGsnt0Natsr9MkiCSJLzdu9vD\nsm3758knJbrWdoz5JYqJc64BsDA1Qc4IAB1SrZjdAO4A8CqARQCmMHNxtJufHZVDr17S5xA24hIo\nk8OQIXIv6UXmqMywEZeAv1GZU6cCGzdKn0PNmnK/PGCA/C8orqMy0/scqlWT12efDZcIt6hIvnvN\nfp2LWcWV9D6HatXcM1Odc45MtTB0qHk+jSgiLgE5ro87zj6uIwguTzE6MvMRzFydmY9k5qeYeSwz\nj029P5KZGzFzU2Y+nZnfS/vsLGY+npnrM/ND0W/+L8mUgxJGEkHloISRRFA5KGElEUQO2iFZpUo4\nSaTLQcmVJLJ1SPpJX5crScQpB6CSRFIqXnJQgkgirByUIJIIKwclqCTCyEEJKolsclDiloTpaUU+\nSSJuOQCVSBA2OSh+JBGVHBQ/kohKDopfSUQhB8WvJExyUOKShMujzHyQRC7kAFQSQbjKQXGRRNRy\nUFwkEbUcFFdJRCkHxVUSLnJQopaEnziH8pREruQAVAJB+JWDYpJEXHJQTJKISw6KTRJxyEGxScKP\nHJSoJBEkCKo8JJFLOQAVXBBB5aBkk0TcclCySSJuOShekohTDoqXJILIQQkriTARkrmURK7lAFTg\njFJh5aC0aSPBP08+mTs5KOmSqFrVnxwOOCB4+DRQJomCAlm3Hzl8/nnw8GmgTBL9+snv33wj9Qki\nB0Ul0aePPTHxbbcBEyeKJDZsEFEFDZ8GyiQxYICs2zShs4bbu0ri3nulXHnIAajgSWtvusl+Mg0a\nJGMZbMltJ06UyLbevc1y0GVecIFdDhpMVVhoLrdwoYyZ6NQJqF8/mmWuXy8RhbZyW7ZIur3GjWU2\nchOPPgr88IOcCDY5DBpkX/eePWWRibayGzZI4Jmt3LZtkpqte3d7JKV+l/372+UQdX3eflvygNrK\n/fCDfO8tWwKXX24u+8ILEjXsMibDNWmt01iMXP+4jMW49lrme+6xx8ZPnMi8erW5zPffS2x8o0bM\nI0ZYV80DB9rLfPop829/yzxlirncjh3Md97JfMYZzL16Me/aZS4/aZJ9LMa6dTIW4+677ds5ZAhz\n06ZS/5ISc9lFi5jPOYe5qMi+XJfxA9Ony/JsYzE2bJDta9qU+ZFH7Mu9+24pv369udyyZcxnncX8\n0kv2ZbrU55VXmM8+m7ldO+bPPrOXdzmOHn+c+bTT7OUWLZJzYvx4e1lm97EYFbYP4oQTZDxC374S\n0huUtWtlvMDQoRIW7PII1NaEnT9fWgQPPGC+0u7cKaNMe/aUOTJcHoHa1r1+vVyVhgyRZreJoUPl\ntqJ9e/dHoGFmFU+nqAj4/nt/ozLbt5dxG7bs0n76JG6/XdYxdaq/7c9ER2X+8Y/ufRK2ffnEExKC\nbrutWLxYkuQ+/HDw2yQvKqwggPCSSJeD3laEDcsOIge9TQoblh1EDuefL3/7iZMIKwmVg85L4UW2\nIdsXXhi9JLp0CSeJzCHbQWcVTycf5ABUcEEAwSWRTQ5KUEmEkYMSVBJh5KDkQhJh5KDkkyS88jmE\nkUS+yAGoBIIA/EvCJAfFrySikIPiVxJRyEGJUxJRyEHJB0nYkr0EkUQ+yQGoJIIA3CXhIgfFVRJR\nykFxlUSUclDikESUclDKUxKumaD8SCLf5ABUIkEAdkn4kYNik0QcclBskohDDkqUkohDDkp5SMJv\nmjgXSeSjHIBKJgjAWxJB5KB4SSJOOShekohTDkoUkohTDkouJRE0h6RJEvkqByCaeTFuSCWs/YKI\n3iOipmmcAzkUAAAgAElEQVTvrUz9/zMiMkc+RUimJMLIQcmURC7koGRKIhdyUMJIIhdyUHIhibAJ\nZrNJIp/lADhEUhLROQC2A3iGmX+RPIyIzgSwiJk3EVEbAPcz82mp91YCaMnMviZTc4mkvOgiiZgz\nsWiRpEfbbz/JaWiTw1VX2e89Z8+WE7VtW7sctGxRkV0OLuteuhT43e/kAJk71y6Hpk0lStImh7PP\nlol0TWzZIhMHn3qqTHxrYswYYOxYiWa0yaF3b4kqfOMNuxxctvPNNyXLs+1+f9s2OSmbN7f3MU2Y\nIBM13367XQ6PPirjLVq18i6zZ48Id8YMkbxNDk2bSpyMTQ7LlgHz5onQbbhGUlrHYjDz20RU1/B+\n+oTjH0CS08bO55/bZzLWK962bcC339rNu3y5fZnbt8vr1q2SBt+EJqpdv97+NMJl3ZrKvrRUrtSm\n/IN79shrnTrm5S5YIHX56itz4M6uXTIm4MMP7dv53/+WvdrKvv66hMyvW2ce4jxsmGyn6/75wx+k\n5Wbihhskj6NtmT/8IK8u9Xn3XUlIe9hh5nJbtpQt22VG7pISOd5M+2jxYjkuoiTqwVq3Apid9jcD\neIOISgGMY2bPaw8R3Q7gdgA4ymH0VY8ectX1QgdePfEE0LChXKEGDzZLorDQvMz584FPPxX7n3GG\n7FhTItzjjgM++kiyYD/0EHDggcHXvX69JC4dNUoS1r7xhj0Rrm2ZCxYAH38s9TGNw9i1SwYNzZgh\nsjUts6hIvptjjgEOPxyoXl1aW14895w9ae2wYVL3d981r7u4WATWv7+0mubPl8F8XjRtKld60zJn\nzABatJABgbVrywlqSoRbWiqDtRo18i4zYoQMknvvPblw1K9vToRbWChjf95/35wId/VqOS6iJLJO\nSiI6HyKI+9L+fRYzN4NMv9cjdbuSFWYez8wtmbnlIbZRNhYyR2VGEZad2efg+gj0gAOkydevn0xg\nEwSVw5AhwL77RpMI13XItsrhjjvsA8kyJ7Xp2lWuaK+8EmwbAZFDkyZyS2kic8j2BRdIn8RjjwVf\nd2afQ5cu0m8RJhHuiBGy/1q3jmeoeNREIggiagJgIoB2zPy/id+Y+bvU6zoA0yDT8cWK15DtMJLw\n6pDMhSQy5aCEkUQu5KCEkURQOShhJOHVIRlGEulyUPJdEqEFQURHAZgKoDMzL037/z5EtJ/+DuAS\nACESoNux5XMIIgnb04o4JeElByWIJHIpByWIJMLKQQkiCdvTiiCSyCYHJZ8lEcW8GAMBHAxgdMbj\nzMMAvJuaL+MjADOZeU4MdQDgnuzFjyRcH2XGIQmbHBQ/kigPOSh+JBGVHBQ/knB9lOlHEiY5KPkq\niSjmxbiNmQ9k5mapn5ap/69gmSujKcu8GbHNi+EqB8VFEq5yUKKUhKscFBdJlKccFBdJRC0HxUUS\nfuMcXCThIgclHyVR4SMp/cpBMUnCrxyUKCThVw6KSRL5IAfFJIm45KCYJBE0CMokCT9yUPJNEhVa\nEEHloGSTRFA5KGEkEVQOSjZJ5JMclGySiFsOSjZJhI2QzCaJIHJQ8kkSFTZpbVg5KOmSWLRITqyg\nclDSE+Feeql3OZVEQYEEzoSRg5Iuia++kvyEUcvh3HODy0Hp2lWiLd9+G/jXvyRC1EUOP/4YLsEs\nIJIARBIlJXKShZ0Ru0sXibh86y2JRbjyynAJZtMT4X75pblseiLc5s2DrzMbFVYQY8fKifDII+Zy\nrklev/5ariQnn2yesn7nTgm4YranDJswQX5cEsc+9xxwyy32sQSDBklwj+1kmjVLhNeokbk+gJxw\nLVrIuAsTTz8tGagbNDDPNl1SIqHJLkyeLK9HHCGBUF7s2QO8+KKITxPDemGaeT2dMWPk9YQTzPXZ\ntElaBC6k1+eDD8xlXbdz+nQJ/rLN3D16tLzeeqvbcl2osIL44x/ltqBnT/OBXaeOnChnnGG+mgNy\nkJgOFJ234qab7Dt36lTgrrtk/IApSlFvK269FdhrL/mMKeKyXj3grLPk1YsFC6Tl0KSJbPMDD5gP\nrtJSKXftteaU7R06AJ99Jq9e6MCrnj3tU9YPHy5S3H9/EZQp4hKQUOMdO6Q+1auby9rWPXOmjNlo\n1gzo3Nm73KZNcnLedZd9mSNGlNWnaVNzxKXLdpaWihCvuw548EHzcubMkXEtJ5xgX6cfKmwfxEEH\nyZVEm+cm7rpL4tRffTX4+vxOarNpk93k6X0ORx4ptwL9+wePuAR+3ufQsKHbI9CqVaUV8de/2puz\nJvyMyhw+HDjpJGkJuT4CPfZYGQB2331yWxSUmTOlz6F7d3M5lcNDD8ntoIn0Pocws4orKodbb5X9\naGLOHLmdNIX9B6XCCgLInSTiloP2ORxwQDhJZOuQdI2TCCuJoHJQciWJuOWghJFEvsgBqOCCAOKX\nRK7koASVhOlpRdySCCsHJW5J5EoOShBJ5JMcgEogCCA+SeRaDopfSbg8yoxLElHJQYlLErmWg+JH\nEvkmB6CSCAKIXhLlJQfFVRJ+JtKNWhJRy0GJWhLlJQfFRRL5KAegEgkCiE4S5S0HxSaJILNsRyWJ\nuOSgRCWJ8paDYpJEvsoBqGSCAMJLIl/koHhJIogclLCSiFsOSlhJ5IsclGySyGc5AJVQEEBwSeSb\nHJRMSYSRgxJUErmSgxJUEvkmByVdEvkuB8AhaW154JK0tnp1YOJE83I2bpRgmOuuk/BnE4MHA0uW\nyMlsk8PNNwNnniknmIlJkySMeORIuxxuvlnKm9i8WeQFSFSjTQ4uy1y6VATwl7+YI0NLS8tkOGKE\nXQ433wx06mSXw803S4CUKZgMkO9w3jx7fZYtk6CiZs3Kvisv/vEPOfEef9wuh5tvlgAxmxxuvlnq\n3KmTudzYsRJp2aePXQ433yyBZDY5vP++pLFbsMBcDogwaW2+0r59WQy6F/Pny+vNN5tzHu7cKVeG\nJUsk2ahtuTVqAL/5jYQ8m8YEVKkioc6XXWZe3tCh8vr55/aDYNo0uZI0bWqOuAQkitNWF71tOfdc\nc7mSEolY3bJFojhNEZeAbJ8txLu4WF47dwZaGg7V3btlDMiXX9rrs3WrvNatK9GzpojL/fcXQRxz\njLR0TDRtas8WvXixvOoFxIvS0rIWUa1a9jo1bQq88IK5zObNcpG7/HJzOb9UWEE0aiQ71gtNMDt+\nvJStUyd7OZ23YsQIGeiyY4dcVU1h2X37Sqtk/HhzItzt24Ff/9q8nUOHSojzYYcBp50m40FMkrjl\nFqBxY2k92RLhXnCBed1FRcCvfiVjRUzlSkpkjMhrr4mc5s6V8iZJtG9vXmZxsaSw/9vfJCTcq+zu\n3XLbMGSIuRwg312VKlKfTp1kbMKQId6S2LNHWiTvvWeXhK0+ixdLyv3Jk2U5XmX1tuLBB2U/7r+/\nHKumxLm2dW/eLPvn+eclQXKURDFxDhHRCCJalppAp3nae62JaEnqvT5RbriJMJPauMZJRJEIN3NS\nm1wmwnUdsp3Z51C1qvw+aVLwsGzXIdsqh65dJUO4icwh266PQKtUke149lngiy/81wWQ42XiRJGR\nqT7Z+hyCziqubN4sx8GDD9pblEFw6eKaBMB059UGwHGpn9sBjAEAIqoKYFTq/RMBdCQiS8M0PFHM\neJULSXjNeJULSQSVgxJGErmQg5ILSYSRgxJUEnHLAXBLOfc2ANOka+0gs24xM38A4AAiOgKSwXpZ\nKvXcTgD/SJWNjSinw4tTErbp8OKURFg5KEEkkUs5KHFKIgo5KH4lkQs5ANE85qwN4Ju0v79N/c/r\n/1khotuJaB4RzVu/fr3vjYhjrsw4JOE6V2YckohKDoofSZSHHJQ4JBGlHBRXSeRKDkAexUGEmTgn\nzol0o5SE34l0o5RE1HJQXCRRnnJQopREHHJQbJLIpRyAaATxHYD0ZwRHpv7n9f9IycUs21FIIugs\n21FIIi45KCZJ5IMclCgkEaccFC9J5FoOQDSCeBnAjamnGacD2MLMawB8DOA4IjqGiGoA6JAqGxm5\nkIMSRhJB5aCEkUTcclCySSKf5KCEkUQu5KBkSqI85AAAYGbjD4DJANYA2AXpR7gVQFcAXVPvE+Rp\nxXIAXwBomfbZywAsTb3Xz7Yu/WnRogXbaNuWuV8/5tJSc7lRo5ivvJL566+ti+TCQnuZxx9nPvVU\ne7kvv2Ru3Zq5ZUvmf/7TXr5/f3uZWbNkecuXm8tt2sTcrRtznTrMY8bYl9u9u/xs2WIut2gR8+TJ\n5jK7dzPfcw9zw4bM994rf5t45hnmyy9nXrrUvp3NmzOPHGkv5/JdfvUV829/yzxxorlcaSnzffcx\nn3ACc69e9vpMnszcpo18V1Fs5/jxsh+7d2feuNFcdtUq5qeesi+TmRnAPHY4F62BUszc0fI+A8jq\nc2aeBWCWXVP+eeUVYL/97NF6mjuySRP7Ml2SiK5fLzN2n3kmcMkl5rJzUvOIzZ0rP14sWCARksuW\nmSM+AQk3btfOHpo8c6YkmF20yJ5LcfRoudING2YuN326RHtqxKAXy5dLmUaN3BLM7ruvXKlNzJgh\ngW9t27olmJ05E/jtb83LnDMHePllyXVpirisUUO+xwYN3OpDBJxqmYV2zhyZiXzTJommNPHNN7Kd\ntq65L76QbO+33GIu54cKG0nZr5+kQO/USZqNXtSpI9GRREC3bvblmg6+9evlAOjVS5Z7wglmSVx9\ntURlmk5mTTBbUCC3Qh07miXhkrS2qEiyKi9aJMu0JcIF5OA75hhzWLhL0lpNTV9QILcOtkS4DRpI\nBKupzjNnAoceClxxhV0OmmC2YUOgdm1zItzOneUWaNs2Wa4tEe7OncA110gEpBeNGkl9GzXyLjNn\nDnDwwRKq/9NPMt2ATRJHHimfMSXCXb1aUu5HSd48xfBLtWpyLzhmjFx5TbRtK1+wpjkPgsrh4Yfl\ninfXXXICvvZa8GWmj8qsXl3u4SdOlDEhQSkqAr7/HvjDH/w9Au3RQ1L/zwrR3kvvc6hePZpEuEFH\nZUaddEb7JP7+d2BhiCmo00dl7refHE+FhTLxjokoEuEGocIKAsidJDLloISRRLYh29WqhZNEuhyU\nXEkiW4dk2ES4YYds55sksg3ZzndJVGhBAPFLwksOShBJmPI5BJVENjkocUvC9LQiqCSiyueQL5Iw\n5XPIZ0lUeEEA8UnCJgfFjyRckr34lYRJDkpcknB5lOlXElEneylvSbgke8lXSVQKQQDRS8JVDoqL\nJPxkgnKVhIsclKgl4WciXVdJxJEJCig/SfjJBJWPkqg0ggCik4RfOSgmSQRJE2eThB85KFFJIsgs\n2zZJxCUHJdeSCJImLt8kUakEAYSXRFA5KNkkESaHpJckgshBCSuJIHJQvCQRtxyUXEkiTA7JfJJE\nhY2DMKGSuO8+CXIxxSq0bSsHy5gx4eWg3HUX8MQTEriydasEN4VJMKuS6NtX6lNSIpmqgshBUUkU\nFNif//foAYwaJcFPmzZJiHsQOSgqiYICqc+aNZLZKm45KF27Sk7I5csltsCLdEnsvbd3OZVEQYFM\nwLx6tdQrTIJZlUSfPvZZ5G+7TS4g06bJ56Kkwgpi0CBJGWY66fbeWw6muXPtUY8TJgD/+Y8EQT36\nqH3drtsIAAMGuEd8mthrL6kPIFcXU9DQnj32qD9ADuQnnpATznYg6jb27+8WUbh7t3lW8Ro1ymat\nttWnpEQS1951lySZddlOG1pu+XJzuXXryiI9bRGXup9t9dm8Wb53G7/6lUS4HnSQ/QKj9UkiKSHR\niTt3irlNX1zdupJf0isnJSAth/XrJQqtfv3wEZeA3FbMmyeZtVu1Mue4dF1mUZGE8NaubY+4BGRg\n2cEH269kNWvaT6riYomiXLtWksHaEvEuXSpRgg8/bJbECSdIxmZTXbTlcMIJwNFHy6A7G7bvcvFi\naeF9+61EctrCsgGJELVFXDZuLFGUpkhKHXjVuLG0VHv3Nq93333t++err+SW7cILzeX8UmH7IBo3\nljDrggK5Wnphuyqmz1vRsWP4iEvg530OvXpJn0SQWcXTKSqSk7NHDznpJkywPwJ1HQVq+46KiyW5\n7EMPyRV8+XL7I1DXeTds606/rbj2Wjnxhg83f8bG4sUybcDgwRLm/N13MtbBhGufhK0+6aMyr74a\nOPlke4vVtsyvvgLGjZPjwna76JcKKwjAXRJeZJvUJmxYdrYOyT/9KZwkVA46D0e1atFLwot0OWif\nwx//GK0kvMjW53DJJeEkkS4HrU+3btFKwotsQ7YvvNBNEl6ky8HUUgtKhRYEEFwSphmvgkrC9LQi\nqCQy5aDkQhLZ5KDELQlTh2RQSWSTgxK3JEz5HIJKIm45AJVAEIB/SbhMh+dXEi6PMv1KwksOSpyS\nMMlBiUsSLk8r/ErCJAclLkm4JHvxK4lcyAGoJIIA3CXhZ65MV0n4iXNwlYRNDkocknCRgxK1JPw8\nynSVhIsclKgl4ScTlKskciUHwFEQtglwiOgeIvos9bOQiEqJ6KDUeyuJ6IvUe+YJN0Nik0SQiXRt\nkggSBGWThKsclCgl4UcOSlSSCBLnYJOEHzkoUUkiSJo4myRyKQfAbWYt6wQ4zPwIMzdj5mYA+gKY\ny8zpc2mcn3rfOlloWLwkEWaWbS9JhImQ9JKEXzkoUUgiiByUsJIIEwTlJYkgclDCSiJMDkkvSeRa\nDoBbC8LvBDgdIXksy41MSYSRg5IpiTByUDIlEVQOShhJhJGDElQSUURIZkoijByUoJKIIsFspiTK\nQw6AW6BUtglwTstWkIh+DZmm7460fzOAN4ioFMA4Zh7v8dnbIVP34aijjrJu1Jw5MoGuF1WrSgDO\ntddK+rWiIsntZ2LCBPMyjz1WTuhBgyRo6KGHzCfiBx9IHkVTyrVLL5XIwPHjRRi33y7r8GLOHImu\n27HDu8xNNwF33y0nqqk+deuW1Wf+fKnP0qXe5d95B3jxRZlt2ouLLhLpjBljXjcgs4+3by/76p13\nJOR6zRrv8lOnmpdZp05ZfT75RORgqs+nnwJTpphnFj/vPKnLk0/a63P++WUBV//+t4yV+f577/K2\n+vzmN1Kfe+4BPv5Y5PDVV97lly8HVq40b6NfonZRWwD/zri9OIuZvyOiQwG8TkSLWabz+xkpcYwH\ngJYtW7JtRR9+KFcLl8ix1q0lMamJoUMl1HrBAnO5LVvkdeNGey6At98G/vUvOchMbNtWtmzb+ouK\nJH8kW74hPdFeeMEcpVg7NdfZdddJK8KL3bvLxn7MnWu+Mp57rpxUtrroQLHSUtk/tivjF1/Yl7l5\ns7yuWWOuDyCzcc+YIa0p05iMVq1EEK7rBuyJcAG3+uh39M039lwaixbZw8b94iIIPxPgdEDG7QUz\nf5d6XUdE0yC3LL8QhF8KC2VgSq1aZkls3y7NT1Oo9dCh0pSbO9ds9AULpOUwYIDMc/Htt+aw7BNP\nlCukKWltURFwzjnS7D7wQGlim8Kyd+ywJ60dNkxup+bOlcFizZqZJbFoEXD99d7v67wVL78sTe4F\nC6S1YZLEokXm71JvK0aPllu3uXPtzWfbMhcvlitt//7ynS5caA7LbtlSwqwXLpSxJSZJ2Na9ebPc\n0o4bBxx2mNRnyBCzJGzLXLpULoT9+gEXXCAtHlNYdnklrXWaAIeIagI4F0BR2v/2IaL99HcAlwAI\nkfLz5/TqJSfsm28GX4brpDaZfQ5RJMLNnNTGdXIeE8OGASedBFx8sb8+CS8yJ7UJM6u4ktnnEDbi\nEvjlpDYXXyx9UbZHoDVryq3IgAEyICwIKgeVZtiIS0DkMHasHJ/VqokgmjcPHnEZFJfZvXdD+hRe\nBbAIwBRmLiairkSUPl/TlQBeY+Yf0v53GIB3iehzAB8BmMnMc6Lb/HCSCCoHJYwkvGa8CiOJdDko\nYSThNeNVGEl4dUiGkYTXjFe5kESmHJQwksiUg1IeknDqf2fmWcx8PDPXZ+aHUv8by8xj08pMYuYO\nGZ9bwcxNUz+N9LNRE0QSYeWgBJGEbTq8IJLIJgcliCRs0+EFkYTtaUUQSdimw4tTEl5yUIJIwksO\nSq4lUWkiKf1IIio5KH4k4TpXph9JmOSg+JGE61yZfiTh+ijTjyRc58qMQxI2OSh+JGGTg5JLSVQa\nQQBukohaDoqLJFzloLhIwkUOiosk/EykC7hJwm+cg4skXOWgRCkJVzkoLpJwlYOSK0lUKkEAZknE\nJQfFJAm/clBMkvAjB8UkCb9yUEySCBoEZZKEXzkoUUjCrxwUkyT8ykHJhSQqnSCA7JKIWw5KNkkE\nlYOSTRJB5KBkk0RQOSjZJBE2QjKbJILKQQkjiaByULJJIqgclLglUWFTztno1Qt47DGJqPvsM+Cq\nq9zksGNHuPBp4OeJcLdvl5nAg8pB0US4H30kP+3aBZODopLo00eySYeRg5KeCHf79mgSzKZLYuXK\ncHJQ9HsbPhw4LWtMsKCS6NtX6hNGDkq6JL79NpwclAsukNdHH5WYlyipsIIYNEgiEG1ZfCenwrYO\nPliCV7z47DO50tsSzO7aJSfBsmXmWcWBslZEvXpuCVldeO45ea1VS8J5Tbgkb61WTcKnW7Swf5fT\npwOffy5XcRM//CDjOm66yT3BbGGhudx33wEvvWRPmFtaKmMgiovNeSEBydQ9apQ9+vCnn6Q+nTrZ\nE8261mfdOoly1TEbXmh9bAmagbJxNUnSWkhuxp9+khPPdGDXqWOPpNSWQ5Mmdjncd5/krjzsMDlg\nTJK4+mppQpoiKUtKgIEDJeLPlmh1+HCZsv7AAyVxbNhEuHpbcd118v3YEuF26CAi7dDBu4zeVtx8\ns6Tmv+su8xW3QQPg/fclGa8tEW6DBub6lJZKS+N3vwOOOUbG4ZjGwXTuDLz3nrx6oQOvbroJ2Gcf\n4M47JcO0F40aSYh9ixbSkjRx7LH2+hQUANdcY8+K/sYbMpbnpJPM5fxSYfsgatUSW/fpUzaWIQiu\nfQ4qhx49pNnrOjmPCZXD4MHA/vubyw4fLjv/oouiibhM73No2DDcrOJKep/D0UfLa//+8c4qrqgc\nbr1VpBBmVnElfVRm3bqyzAEDZByOCdfJeUyoHG680d4SeuMNGdfhku3bLxVWEEB4SQSRQ/368j8/\nM3hlI6gclDCSyNYhGXRWcSVbh+QBB+RGEulyaNhQ/hd0VnEl25DtmjVzI4l8kQNQwQUBBJdEGDko\nQSURVg5KEEmYnlYElYTpaUXcksgmByWoJEz5HOKWRD7JAagEggD8SyIKOSh+JRGVHBQ/knB5lOlX\nEi6PMuOShEkOil9JuCR7iUsS+SYHoJIIAnCXRJRyUFwlEbUcFBdJ+IlzcJWEnziHqCXhIgfFVRJ+\nMkFFLYl8lANQiQQB2CURhxwUmyTikoNikkSQICibJIIEQUUlCT9yUGySCJImLipJ5KscgEomCMBb\nEnHKQfGSRNxyULJJIkyEpJckwkRIhpVEEDkoXpIIk0MyrCTyWQ5AJRQE8EtJ5EIOSqYkciUHJV0S\nYcOngV9KIooEs0ElEUYOSqYkokgwG1QS+S4HACC2JTeEzIsB4AkAVQFMZOaHM94/D5JJ6uvUv6Yy\n8wMun81Gy5Yted488xQa++wjYwlMbNggJwwgkXAmOezeLcEwDzxgl8MNN5RNB29bHiCysMnhhhuA\n3//eLocbbpAArTZtzOWGD5fZxe+/3y4Hv/UZNcouB5dlbt4sArj0UjlJTDz+uKSTGzDALocbbhAB\nmMKXS0vL1jlypF0OLvXZskXCqK+4QoLOTIwaJUFaBQV2Odxwg6THs8nhww+Bd9+V4QU2iOgTl2ko\nrJGUafNiXAzJaP0xEb3MzJl3cu8w8xUBP+ubCy80ZyMGyhKCHnGEHFS2kxSQ1odtuS1a2MukjwKs\nU8ftCl6zpn25bdsCr78uEX01a3qXe/ZZiTw86KBo6rNpE1CjBrBzp9THFHHpukwN2e7Y0Vy2tFT2\nIeBWn/32k/BxW47L554D/vIX+5iWUaPk9Z137CfpjTdK0l5bfQ49VH4/8EB7ferXt8tp61ZpkV1x\nhbmcX1xCrf83LwYAEJHOi+Fykof5rJHmzSWi0YsFC8Sk48fLDpgwQQ4YkyQKC4H//lcyCJsS4V5x\nhXndJSUSs//uu7KsN9+U8iZJFBbKSfXmm+ZEuNdcI7csOueDSRKFhXIwfv21OSzbVp9Nm6Q+c+fK\nbdPcuVLeJAnbMhcvlhPub3+TUbZeZfW2YsgQmSeialXZBlNY9t13SxLe8ePNkqhaVVpGpu0cNUoE\ncuCBMuhu5kyzJK64QiI5TfUpKJBbmqZNpSW8eHFZuvxsdOpk3satW2X/vPSSTF8QJS59ENnmxaid\npdyZRLSAiGYTkTaaXD8LIrqdiOYR0bz169c7bJY3mX0OtWrJrUOfPvJlmgibCDezz6F69bI+CdOc\nBoB7ZqqDDpIDrG/fsjT8XgSdVVzJ7HOIIhHu4sXAU0/JckyjMrP1OcQ9q3g6o0ZJyLhelYPOKq5k\n63NwnZzHi61b5Th44AFzVu6gRNVJ+SmAo5i5CYAnAUz3uwBmHs/MLZm55SGHHBJ4Q7w6JHMhCa8O\nyYoqCa8OyTCSCCMHJReSyJSDElQSpg7JoJKIWw6AmyCs82Iw81Zm3p76fRaA6kRUy+WzUWJ7WhGn\nJGxPKyqaJGxPK4JIIgo5KHFKwksOil9JuDyt8CuJXMgBiGheDCI6nIgo9fupqeWWuHw2KlwfZcYh\nCddHmRVFEq6PMv1IIko5KHFIwiYHxVUSfh5lukoiV3IAopsX4xoAC1PzX4wA0IGFrJ+NuhJ+08RF\nKQk/cQ5A/kvCb5yDiyTikIMSpSRc5aDYJOFHDopNErmUAxDRvBjMPDI170VTZj6dmd8zfTZKguaQ\njEISfuWg5KskggZBmSQRpxyUKCThVw6KlySCyEHxkkSu5QBU8EjKsAlmw0giqByUfJNE2AjJbJLI\nhcUt5pUAABNjSURBVByUMJIIKgclUxJh5KBkSqI85AAAYOa8+2nRogXbuPJK5j59mEtLzeUmTGBe\nvdpcZv165m7dmHv2tK6WH32UuWFD5u7dmbdsMZddsID5xRfNZXbulPVef7193S+/zHzqqczLl5vL\nlZRIfe66y77M4cOZGzSQ+mzaZC67aBHzc8+Zy+zaxdyrF3PTpsy9ezPv3m0u/8wzzG3byrJtDBxo\nLzNiBPMpp9jLLVnC3L49c/PmzK+8Es26X32V+bTTmC+/nHnhwmiWOXo0c716sn82bDCXXbWK+amn\n7MtkZgYwjx3OxQqbk3LaNPcEs2ecIZY38corkmXY1hooKZErY40aknrexHPPSevAlsD0668lIawt\n2nL3bsloXb++fZlVq0rwDJE5mAqQK37r1u4JZm0dkitXSnLb4483J5hNX2bz5uZyzLKvN2+2h0V/\n/LF8P9Jt7s301MP4Dz6QsHQTDzwgLS5TjktAwp0BCc574QXvcq+/LqHWtm3cuhVYsUJaEran/198\nAaxdmyStBSAHgCnh565dwL33SrOsWzdz0tqiIhmw8/339uzTmmC2dm1JMW6KuHRJWltcDPz4ozTv\nXRLM9u8vt0Y1apgjLgE5+HbulO21ScKWMHfECGlCH364OWmtzluhTXhbItwGDeSgrlpVbhNMbNki\n9bElwgXs9Rk1SgS6fr18/7aw7NJSSWxsS4TbqJG8b7q1eP112R8XX2zezq1b5Xu84w5ZXu3a5kS4\nq1fLoK4oqdB9EF6oHO64wz7wynVSm8w+h7ARl4DI4ZlnpJVjukfPHJXpp0/iz3+W+2Fbn4SJESOk\nNdC6tblc+qQ2e+3lnpmqTRvJ8Kxp270IM6t4OqNGSRLayy93fwRatarUZ9KkcIlwX38dWLjQPqZD\n5fDnP0ufQxSJcINQ6QSRCzkoYSQRVA5KriQRRA5aHz/p63IliXQ5KLmSRFA5KOUhiUoliFzKQQki\nibByUOKWRBg5KPkkiWxyUOKWRFg5KLmWRKURRHnIQfEjiajkoMQliSjkoOSDJExyUOKSRFRyUHIp\niUohiPKUg+IiiajloEQtiSjloJSnJFzkoEQtiajloORKEhVeEPkgB8UkibjkoEQliTjkoJSHJPzI\nQYlKEnHJQcmFJCq0IPJJDko2ScQtByWsJOKUg5JLSQSRgxJWEnHLQYlbEhU2DsKvHKpVi18OSq9e\nwGOPycGzdasEN8UtB6VtWzlY/EiiRg1/cpg3L5gcFJVE376SUckUJNWmDTB7trskCgokjD2MHJR0\nSfzqV97lVBJan+++k8CuuOWgdO0KjB0rcRB77x18OdlwSlqba1yS1hJJMk/T7NqA+3Tsu3fLOISe\nPe1yGDTIvrz0dffvbz+ZBg2SNGk2OQwaBDRubA6+AiSV24oV9u3cuFFOvoYNZZZvE08/LSn0bPXR\nKetdv3PAXnb2bBGtrdzmzRIAVa+eedZuAHjtNZlZ3LbMxYuB55+XFootmOrBB+V32zI1MvSOO+xy\n8Hu8uZzSrklry33cRbYfl7EY7dszFxTYx2KMHs189dX28QvMMn5g3Dh7ucJCe5mFC5kvu0zGBbzx\nhr28S1z+jBkS62+ry8aNErt/5532ZT7xBPOxx0r5zZvNZV3GYuzeLWMwrrvOvu7Zs+X7aduWefFi\ne3mX72jkSOb69aU+Gzeayy5bJmNBTOzezXzvvcxNmsgYk127zOWffVbGYhQX27f1jDOYhw2zl3M5\n3pYtk+N8zBh7WWb3sRgVtg+iaVO54g0YAOzZ412uRg1g6FC5Sq5YYV5m+/bAYYdJstMwFBfLFfyh\nh+S2weURqC0mf+ZMYNUqSatuIn1Upi0cWW8rbrjB/RGoaTv9jMqcM0fGqfTs6d4nYfuO9LaiUyf3\neTdMpI/KvPJKtz6J6tXd5wK95BJpDQbNcaksXy51HzJEjvcoqbCCAIAmTdwkUb26uyTatQsniXQ5\naDM8bFh2EDnYhmxn9jmEjbgMIgcdexF0VvF0Mvsc/EzOk41sQ7b9dFy6SuLii8NJIl0O1asHW4YJ\nJ0EQUWsiWkJEy4ioT5b3b0hltP6CiN4joqZp761M/f8zIrKMmfNPPkkimxyUoJLIhRyUoJIIIwcl\njCS8OiSDSsKUzyGfJBG3HAAHQaRNftMGwIkAOhJR5ni2rwGcy8wnAfgzgMxT63xmbsYunSIByAdJ\nmOSg+JVELuWg+JVEFHJQgkjC9rTCryRckr3kgyRyIQfArQXxv8lvmHknAJ385n8w83vMrF//B5Ds\n1TmlPCXhIgfFVRLlIQfFVRJRykHxIwnXR5mukvCTCao8JZErOQDRTpyj3ApgdtrfDOANIvqEiG73\n+lAUE+eUhyT8yEGxSaI85aDYJBGHHBQXSfiNc7BJIkiauPKQRC7lAETcSUlE50MEcV/av89i5maQ\nW5QeRHROts9yRBPn5FISQeSgeEkiH+SgeEkiTjkoJkkEDYLykkSYHJK5lESu5QBENHEOABBREwAT\nAbRj5v9NXcvM36Ve1wGYBrlliZVcSCKMHJRMSeSTHJRMSeRCDko2SYSNkMyURBQJZnMhifKQA+AQ\nSUlE1QAsBXAhRAwfA7ie0+a3IKKjAPwTwI2clvKeiPYBUIWZt6V+fx3AA8w8x7ROl0jKK68EHnnE\nWAQLFkjEYVGROU2Yhm1/+KHkCTRRVAT07i3LtcnhvfdkMtUePczLfOwxCZXt2dMuh27dJMb/lVfs\ncjjuOJGfTQ6dOgF//7u5zMaNwPnny+TGr7xil8Nxx8kM5DY5PPaYLNcUbr17t4QxT58ueTNtcnCp\nz+bN8r2UlMhybXK44gpgxgxzmaVLZdsmTQJatfIup1KaOtWe0fz11+WYaNvWLofly+UYHjjQvEzA\nPZLSOhaDmXcTkU5+UxXA05yaOCf1/lgAAwEcDGB0aoKt3amVHwZgWup/1QA8Z5ODK3Pm2MONS1Lt\nmHbtxNqm1PhXXSUHwAcfmJe5YYO8fv+9JEc1MWMG8NZbMmu3CW3ubtpkX7/mHCwutoc7A5Jz0bbM\nDz+0lyktFTkA0srZvNlcHhCp2Jb7zDOSE3LnTnO5b7+V1w0boquPHh+rVwPbtpnLL1kiJ55pVnHl\n+eftrco1a+TVtp3r1snrihUyW72JRYtk5vMocRqsxTLf5qyM/41N+/02ALdl+dwKAE0z/x8F990n\nVwovdODVsGFyNfv3v+3zZ6xaZV5mcbFkDu7fX076deukaemFRnuaRDZjBnD22XJF3m8/4IgjzIlw\n9+yRlsO//20fALZ0qRz4p51mHuOxbJm53npbMW0a8J//SAtm8GBzItxlyySS85BDgEsv9S5Xtaq0\nHkyJbUeNkuSu778vMjnlFHN5l/oUFEhrcNUqObFskafLlsmYiYMOMktir73sSWvHjZPjd/5883Yu\nXy5SKCgAzjvPPjo0SVrrSOaozKZNJRNzv37mPgkTmUO227cHDj00XFj2jBlygOotiOsj0MaNZSBS\nQUFZSyEbfibn8SKzz8HP5Dx33SUDnfzOKp5Oep9DmFnFlcw+Bz+ZqVwn5zExbpxkBm/Xzlxu+XJg\n5EjpI6tePXzEZVAqnSC8hmyHkYRXPocwksiUg5JPkvDqkMyVJLJ1SIaRhFeHZK4kEVQOSnlIolIJ\nwpbPIYgkbMlegkjCSw5KPkjC9rQibkmYnlYEkYTtaUXckggrByXXkqg0gnBN9uJHEq6ZoPxIwiYH\npTwl4fooMy5JuDzK9CMJ10eZcUkiKjkouZREpRCE30xQLpJwlYPiIglXOSjlIQm/E+lGLQk/cQ4u\nkvAb5xC1JKKWg5IrSVR4QQRNE2eShF85KCZJ+JWDkktJBJ1lOypJBAmCMkkiaBBUVJKISw5KLiRR\noQURNodkNkkElYOSTRJB5aDkQhJB5aCElUSYCMlskggbIRlWEnHLQYlbEhVWEFElmE2XxMKF4eSg\npEsirByUOCWxeHE4OShBJRFFgtl0SXz5ZfjwaSC4JHIlByVOSVTYrNYjR0qQ0rBh5nKDBskXeOaZ\n5nJTpki0Wr9+blPWr1ghiVFt2wgARx9tnsV50yYZM7FtmwRLmXjmGfmxJTFdvhx46SW76GrUkMi/\nFi0kI7OJ6dOBzz+XE9vEDz+IeGwzigMS4QqIXGyRqZqU1USNGsALL8gJs/fe8rsXU6ZIkNTy5fZl\nPvGEPYUfIGHWgNRn/nzvcswSDl5QUJa4Nxt79sjxWFxsl91TT0m06y232LfTlQoriG7dJCqyd29g\n3329y9WpI7n/6tTxLlNcLBF6S5bIjhs40BxxuWuX7IgbbzRL4uqrJZrRFEmpA686dQJ++klEYpJE\nvXoSSnz22eaIS0CWZRKT3lZcfbVIrFMnc5bwDh2kFdOhg3eZzZtFsnfeaV43AIweLfkw999fBH7J\nJebygL0+BQUyTqd+faBjR3PEZefOMqDv4IPtY0YOOMBen3Hjyupz6qn2sGwi8zL37JHvsl072e/X\nXmseU3TLLUkk5f849FA5kfv0AbZvD76c9D6HZs2iHwVqIn1UZv36coXs08c+LiBsjkvg530OjRuX\n3W6EieVXOTz4oP1qO3q0SOnyy+V2Y9EiSUMflPQ+hyZNop+cx0b6bUWPHsDXX4eLuFQ5dOwox6Xr\nKNCoqbCCAMJLIluHZBxDxbORbch2rVq5kUS2Dslq1cJJIqgclDCSyNYhGccMXl5k63MII4l0OTRp\nIv/zM1Q8Siq0IIDgkjA9rYhbEqZ8DnFLwvS0IqgkwspBCSIJ09OKXEjC1CEZRBLZ5KCUhyQqvCAA\n/5JweZQZlyRckr3EJQmXR5l+JRGVHBQ/knB5lBmnJFyeVviRhEkOSq4lUSkEAbhLwk+cQ9SS8JMJ\nKmpJ+IlzcJVE1HJQXCThJ84hDkm4PsoE3CThIgcll5KoNIIA7JIIEgQVlST8yEGJShJBgqBskohL\nDopJEkGCoKKUhB85KCZJ+JGDkitJRDVxDhHRiNT7C4iouetno8ZLEmEiJMNKIogclLCSCBMh6SWJ\nuOWgZJNEmAjJKCQRRA5KNkkEkYOSC0lENXFOGwDHpX5uBzDGx2cjJ1MSYcOngeCSCCMHJagkwoZP\nA7+URK7koKRLIooEs2EkEUYOSrokwshBiVsSLklrzwBwPzNfmvq7LwAw8+C0MuMA/IuZJ6f+XgLg\nPAB1bZ/NhkvS2rPOkoPHxLp1Mr06AEyebJfDtddKdJ2JBQvk5Hj+efNksrt2SdAMICHXNjm4rHvD\nBklg2rMncMYZ5rKPPgp89JFMMW+Tg8u6d+8Grr9efh83zi6Ha68FbrrJLochQyTgq6Ulferjj0sS\n4MJCuxz81GfvvcuiOb2YPRv4y19kf9rkMGwYcPrp9sjdkSOBt9+Wi4dNDi71KS2VQLerrrKXBdyT\n1roI4hoArVN5J0FEnQGcxsx3pJWZAeBhZn439febkLkx6to+m7aM2yGtDxx11FEtVq1aZdyuDz6w\nhyUDciVt2dIt7NclnBWQHXtO1tk9fs6aNTIO4/TTo1v3Z58BxxzjVp9p0ySqMKp1r1sng7tMGZsV\n1+9o9WqRp8t4mqjrs2ULsHKljMeJat3ffCN1KY/989NPctEyZQhXKpwg0nFpQSQkJAQnsrT3cJs4\nx6tMdYfPJiQk5CkuTzE+BnAcER1DRDUAdADwckaZlwHcmHqacTqALcy8xvGzCQkJeUpUE+fMAnAZ\ngGUA/gvg96bPxlKThISEyLH2QZQHSR9EQkK8uPZBVKpIyoSEhGhJBJGQkOBJIoiEhARPEkEkJCR4\nkpedlES0HoA5lBKoBWBDDjYnbipDPSpDHYDKUQ/XOhzNzIfYCuWlIFwgonkuvbD5TmWoR2WoA1A5\n6hF1HZJbjISEBE8SQSQkJHhSkQXhMJd2haAy1KMy1AGoHPWItA4Vtg8iISEhfipyCyIhISFmEkEk\nJCR4kveCCJMwN59wqMd5RLSFiD5L/Qwsj+30goieJqJ1RLTQ4/2Ksh9s9cjr/QAARFSHiN4ioi+J\nqJiIfpF8MbL9wcx5+wMZIr4cQD0ANQB8DuDEjDKXAZgNgACcDuDD8t7ugPU4D8CM8t5WQx3OAdAc\nwEKP9/N+PzjWI6/3Q2objwDQPPX7fgCWxnVe5HsL4lQAy5h5BTPvBPAPAJlpQ9sBeIaFDwAcQERH\n5HpDLbjUI69h5rcBbDQUqQj7waUeeQ8zr2HmT1O/bwOwCEDtjGKR7I98F0RtAN+k/f0tfvlFuJQp\nb1y38cxUc3A2EQVM7F5uVIT94EqF2Q9EVBfAyQA+zHgrkv3hkpMyITd8CuAoZt5ORJcBmA6ZZyQh\nt1SY/UBE+wJ4CcCfmHlrHOvI9xZEmIS5+YR1G5l5KzNvT/0+C0B1IqqVu00MTUXYD1Yqyn4gouoQ\nOTzLzFOzFIlkf+S7IMIkzM0nrPUgosOJZCoeIjoVsm9Kcr6lwakI+8FKRdgPqe17CsAiZh7mUSyS\n/ZHXtxgcImFuPuFYj2sAdCOi3QB+BNCBU93R+QARTYb08Nciom8BFEKmNagw+wFwqkde74cUrQB0\nBvAFEX2W+l8BgKOAaPdHEmqdkJDgSb7fYiQkJJQjiSASEhI8SQSRkJDgSSKIhIQETxJBJCQkeJII\nIiEhwZNEEAkJCZ78PxPEltvx/oKtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xb6d444cc>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import matplotlib.tri as tri\n",
    "import math\n",
    "import numpy.random as rd\n",
    "\n",
    "\n",
    "Nx=18\n",
    "Ny=18\n",
    "\n",
    "x_m=0\n",
    "x_M=2\n",
    "y_m=0\n",
    "y_M=2\n",
    "\n",
    "x=np.linspace(x_m,x_M,Nx+2)\n",
    "y=np.linspace(y_m,y_M,Ny+2)\n",
    "\n",
    "X,Y=np.meshgrid(x,y)\n",
    "\n",
    "X=X.flatten()\n",
    "Y=Y.flatten()\n",
    "\n",
    "triang = tri.Triangulation(X, Y)\n",
    "\n",
    "NTri=np.shape(triang.triangles)[0]\n",
    "NSom=np.shape(triang.x)[0]\n",
    "\n",
    "#Tableau avec coordonnes des noeuds\n",
    "TabSom=np.zeros([NSom,2])\n",
    "TabSom[:,0]=triang.x\n",
    "TabSom[:,1]=triang.y\n",
    "\n",
    "# Tableau avec noeuds des triangles\n",
    "TabTri=triang.triangles\n",
    "\n",
    "plt.figure(1)\n",
    "plt.gca().set_aspect('equal')\n",
    "plt.triplot(X,Y,triang.triangles, 'b-', lw=0.5)\n",
    "plt.title('maillage')\n",
    "\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Les matrices élémentaires **\n",
    "\n",
    "Pour construire les matrices de masse M et de rigidité K, on commence par définir des matrices élémentaires qui permettent, localement sur un triangle $T_N$ de noeuds $S_{I_1},\\ S_{I_2}$ et $S_{I_3},$ de calculer \n",
    "$$\n",
    "\\int_{T_N}\\nabla\\Phi_{I_i}\\dot\\nabla\\Phi_{I_j}\\ \\ \\ \\textrm{et}\\ \\ \\ \\ \\int_{T_N}\\Phi_{I_i}\\dot\\Phi_{I_j} \n",
    "$$\n",
    "pour $i,\\ j=1,2,3.$ Ces intégrales vont contribuer respectivement à l'élément $K_{I_i,I_j}$ et $M_{I_i,I_j}$ des matrices globales.\n",
    "\n",
    "** Q3)** Construire deux fonctions $M_{elem}(\\mathcal{S1},\\mathcal{S2},\\mathcal{S3})$ et $K_{elem}(\\mathcal{S1},\\mathcal{S2},\\mathcal{S3})$ calculant les matrices de masse et de rigidité élémentaires sur un triangle $T$ de sommets $\\mathcal{S1}=(x_1,y_1),\\ \\mathcal{S2}=(x_2,y_2)$ et $\\mathcal{S3}=(x_3,y_3).$  \n",
    "Pour calculer la matrice de masse élémentaire sur le triangle $T,$ on peut utiliser les formules ci-dessous, que l'on vérifiera. \n",
    "Les fonctions de base associées à chaque sommet de $T$ sont localement données par\n",
    "$$\n",
    "\\lambda_1(x,y)=\\frac{1}{D}(y_{23}(x-x_3)-x_{23}(y-y_3)),\\ \\ \\ \\lambda_2(x,y)=\\frac{1}{D}(y_{31}(x-x_1)-x_{31}(y-y_1)),\\ \\ \\ \\lambda_3(x,y)=\\frac{1}{D}(y_{12}(x-x_2)-x_{12}(y-y_2)),\n",
    "$$\n",
    "où $x_{ij}=x_i-x_j,\\ y_{ij}=y_i-y_j$ et $D=x_{23}y_{31}-x_{31}y_{23}.$ $|D|$ est égal à deux fois la surface du triangle.\n",
    "\n",
    "Pour $k_1,\\ k_2,\\ k_3\\in\\{0,1,2\\},$ on a \n",
    "$$\n",
    "\\int_{T}\\lambda_1^{k_1}\\lambda_2^{k_2}\\lambda_3^{k_3}=2\\frac{k_1!k_2!k_3!}{(k_1+k_2+k_3+2)!}aire(T).\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# a compléter\n",
    "#def M_elem(S1,S2,S3):\n",
    "    #...\n",
    "    \n",
    "#def K_elem(S1,S2,S3):\n",
    "    #...\n",
    "\n",
    "# A LA FIN : vérifier le calcul des matrices élémentaires sur le triangle de référence de noeuds (0,0), (0,1) et (1,0). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Assemblage des matrices $M$ et $K.$ **\n",
    "\n",
    "Dans cette partie on met en oeuvre un algorithme pour construire les matrices $M$ et $K.$ Pour ce faire on remarque que, par exemple, \n",
    "$$\n",
    "M_{I,J}=\\sum_{N=1}^{NTri}\\int_{T_N}\\Phi_I\\Phi_J=\\sum_{N : S_I,S_J\\in T_N}\\int_{T_N}\\Phi_I\\Phi_J.\n",
    "$$\n",
    "\n",
    "L'algoritme de construction des matrices consiste alors à faire la boucle suivante :\n",
    "\n",
    "    Pour N= 1...NTri\n",
    "         Détérmination des sommets S_I1, S_I2 et S_I3 du triangle T_N\n",
    "         Calcul des matrices élémentaires associées au triangle T_N\n",
    "             Pour i=1..3\n",
    "                Pour j=1...3 \n",
    "                    M(Ii,Ij)=M(Ii,Ij)+M_elem(i,j)\n",
    "\n",
    "**Q4)** Compléter dans le programme l'assemblage des matrices M et K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "K = np.zeros([NSom,NSom])\n",
    "M = np.zeros([NSom,NSom])\n",
    "# A compléter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Calcul du second membre $F.$**\n",
    "\n",
    "On peut calculer exactement les composantes $F_I$ du vecteur $F$ ou approcher ses valeurs en utilisant par exemple des formules de quadrature.\n",
    "\n",
    "Ici on fait le choix de remplacer $f$ par son interpolé $P^1$ aux points du maillage, *i.e.* on approche $f$ par la fonction $\\displaystyle{\\sum_{I=1}^{NSom}f(S_I)\\Phi_I}.$ \n",
    "\n",
    "**Q5)** En approchant $f$ par son interpolé, donner une approximation du second membre $F$ faisant intervenir la matrice de masse.  \n",
    "\n",
    "On admet que cette approximation ne change pas l'approximation par éléments finis du problème. \n",
    "\n",
    "**Validation : calcul d'une solution connue**\n",
    "\n",
    "On considère $f$ tel que la fonction\n",
    "$$\n",
    "u(x,y)=\\cos(\\pi x)\\cos(2\\pi y)\n",
    "$$\n",
    "est solution du problème (PN) dans $\\Omega=[0,2]\\times[0,2].$ \n",
    "\n",
    "**Q6)** Construire dans le programme une fonction $f(x,y)$ définissant le second membre $f$ et calculer l'approximation du vecteur $F$ obtenue comme expliqué dessus. Calculer le vecteur $U$ des coefficients de la solution approchée donnée par la méthode des élements finis et utiliser le code suivant pour la visualiser."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#plt.figure(2)\n",
    "#plt.gca().set_aspect('equal')\n",
    "#  #plt.tripcolor(triang, U, shading='flat')\n",
    "#plt.tripcolor(triang.x,triang.y,triang.triangles, U, shading='flat')\n",
    "#plt.colorbar()\n",
    "#plt.title('solution approchée par EF P1')\n",
    "\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q7)** Soit $\\Pi_h u=\\displaystyle{\\sum_I=1^{NSom} u(S_I)\\Psi_i}$ l'interpolé $P1$ de la solution exacte $u$ de $(PN)$ aux points $S_I.$ Calculer, pour différentes valeurs de $h,$ la norme $L^2$ de l'erreur $u_h-\\Pi_h u.$ Remarquer que cette intégrale peut se calculer en utilisant la matrice de masse $M.$ Evaluer l'ordre de précision de la méthode en norme $L^2.$ \n",
    "\n",
    "**Q8) Une autre structure de maillage** \n",
    "Remplacez le debut de votre programme par le code ci-dessus, pour générer un maillage non structuré du carré."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Nx=18\n",
    "Ny=18\n",
    "\n",
    "x_m=0.\n",
    "x_M=2.\n",
    "y_m=0.\n",
    "y_M=2.\n",
    "\n",
    "Rx=rd.random([Nx+2,Nx+2])\n",
    "Rx[:,0]=0\n",
    "Rx[:,-1]=0\n",
    "Ry=rd.random([Ny+2,Ny+2])\n",
    "Ry[0,:]=0\n",
    "Ry[-1,:]=0\n",
    "x=np.linspace(x_m,x_M,Nx+2)#+4*((x_M-x_m)/(5*(Nx+2)))*Rx\n",
    "y=np.linspace(y_m,y_M,Ny+2)#+4*((y_M-y_m)/(5*(Ny+2)))*Ry\n",
    "\n",
    "X,Y=np.meshgrid(x,y)\n",
    "X=X+0.9*((x_M-x_m)/((Nx+2)))*Rx\n",
    "Y=Y+0.9*((y_M-y_m)/((Ny+2)))*Ry\n",
    "print X\n",
    "print Y\n",
    "\n",
    "X=X.flatten()\n",
    "Y=Y.flatten()\n",
    "\n",
    "triang = tri.Triangulation(X, Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Q9) Domaine coronne** On considère le problème (PN) dans une coronne. Pour cela remplacer le debut de votre code par le code ci-dessous, qui crée le domaine de calcul et le maillage. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def triangulation_couronne(nrad=8,nang=36):\n",
    "    radii = np.linspace(0.25, 1, nrad)\n",
    "    angles = np.linspace(0, 2*np.pi, nang, endpoint=False)\n",
    "    angles = np.repeat(angles[..., np.newaxis], nrad, axis=1)\n",
    "    angles[:, 1::2] += np.pi/nang\n",
    "    x = (radii*np.cos(angles)).flatten()\n",
    "    y = (radii*np.sin(angles)).flatten()\n",
    "\n",
    "    # On construit une triangulation de Delaunay et on filtre les triangles proches de l'origine\n",
    "    triang = tri.Triangulation(x, y)\n",
    "    xmid = x[triang.triangles].mean(axis=1)\n",
    "    ymid = y[triang.triangles].mean(axis=1)\n",
    "    mask = xmid*xmid + ymid*ymid >= .25*.25\n",
    "    T = triang.triangles[mask,:]\n",
    "    return x,y,T\n",
    "\n",
    "[x,y,T] = triangulation_couronne(8,36)\n",
    "\n",
    "\n",
    "NTri=np.shape(T)[0]\n",
    "NSom=np.shape(x)[0]\n",
    "\n",
    "TabSom=np.zeros([NSom,2])\n",
    "TabSom[:,0]=x\n",
    "TabSom[:,1]=y\n",
    "\n",
    "# noeuds des triangles\n",
    "TabTri=T\n",
    "\n",
    "# .... Le programme\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}