3-PoissonEquation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "101f55ab-c566-41aa-a4ce-4d3888910e3a",
   "metadata": {},
   "source": [
    "# Your first FEM solver- poisson equation solver\n",
    "\n",
    "- [x] isoparametric element(mesh)\n",
    "- [x] dirichlet bc\n",
    "- [x] elemen loop\n",
    "- [x] matrix/vector assemble\n",
    "- [x] 1d/2d poisson equation\n",
    "\n",
    "Author: Yang Bai @ M3 Group\n",
    "\n",
    "Date  : 2022.01.23\n",
    "\n",
    "QQ group: 628204857"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0dd37a11-fb46-41e4-9fd5-b57ed1e1a937",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## Poisson equation\n",
    "\n",
    "$-\\Delta\\phi=f$ or in 1d case: $-\\frac{\\partial^{2}\\phi}{\\partial x^{2}}=f$\n",
    "\n",
    "Boundary conditions:\n",
    "\n",
    "$\\phi(x=x_{1})=y_{1}$ and $\\phi(x=x_{2})=y_{2}$\n",
    "\n",
    "The analytical solution in 1D case can be read as follows:\n",
    "\n",
    "$\\phi(x)=-\\frac{f}{2}x^{2}+ax+b$, with $a=\\frac{y_{1}-y_{2}}{x_{1}-x_{2}}+\\frac{f}{2}(x_{1}+x_{2})$ and $b=y_{1}-ax_{1}+\\frac{f}{2}x_{1}^{2}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8051026-fd0c-47b6-b4ba-d61009d20cd4",
   "metadata": {},
   "source": [
    "# isoparametric mesh\n",
    "\n",
    "$\\int_{\\Omega}f(x,y,z)dxdydz=\\int_{\\Omega}f(x,y,z)d\\xi d\\eta d\\zeta J$, with $J=\\mathrm{det}(\\frac{\\partial\\mathbf{x}}{\\partial\\mathbf{\\xi}})$\n",
    "\n",
    "In 1D case, the jacobian matrix looks like: $J=\\frac{\\partial x}{\\partial\\xi}$.\n",
    "\n",
    "In 2D case, it should be:\n",
    "$\\mathbf{J}=\\begin{bmatrix}\n",
    "\\frac{\\partial x}{\\partial\\xi} & \\frac{\\partial y}{\\partial\\xi} \\\\\n",
    "\\frac{\\partial x}{\\partial\\eta} & \\frac{\\partial y}{\\partial\\eta}\n",
    "\\end{bmatrix}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "31c260af-f4f6-4e16-a2c9-1320c0c52c48",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f1c307c43a0>]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAPbUlEQVR4nO3db4hl9X3H8ffHmSgdkiYhOwnt6s7YYv5IiSVOjC1N3TS0We0DCeSBZolUEgbXGPJQW2hSCEsbQiEE3W4GESkO+qCRxIRNpFB2U7DbOoJRN2LYGnfdGuqYhBQMKKvfPpibZJy9M/eMe2Zm5+f7BQNzzvntPb+fs7z3eO69c1NVSJK2v/O2egKSpH4YdElqhEGXpEYYdElqhEGXpEaMb9WJd+zYUdPT01t1eknalh555JEXqmpy2LEtC/r09DQLCwtbdXpJ2paSnFjtmLdcJKkRBl2SGmHQJakRBl2SGmHQJakRI4Oe5K4kzyd5YpXjSfK1JMeTPJbkA/1Pc8n8/DzT09Ocd955TE9PMz8/v1GnkqTebXTDulyh3w3sWeP41cAlg69Z4J/Oflpnmp+fZ3Z2lhMnTlBVnDhxgtnZWaMuaVvYjIaly6/PTTINfKeq/mDIsa8Dh6vq3sH2U8DuqvrJWo85MzNT63kd+vT0NCdOnPnyywsuuIArr7yy8+NI0lY4evQoL7300hn7p6ameOaZZzo/TpJHqmpm2LE+7qHvBJ5dtn1qsG/YRGaTLCRZWFxcXNdJTp48OXT/sP9AknSuWa1Vq7Xt9ejjnaIZsm/oZX9VzQFzsHSFvp6T7Nq1a+gV+tTUFIcPH17PQ0nSplvtLsOuXbt6O0cfV+ingIuWbV8IPNfD477G/v37mZiYeM2+iYkJ9u/f3/epJKl3m9GwPoL+AHDD4NUuVwK/GHX//PXYu3cvc3NzXHDBBcDSlfnc3Bx79+7t+1SS1LvNaNjIJ0WT3AvsBnYA/wt8EXgTQFUdTBLgdpZeCfNL4MaqGvls53qfFP2V3bt3A3ibRdK2dLYNW+tJ0ZH30Kvq+hHHC/js65qZJKk3vlNUkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhrRKehJ9iR5KsnxJLcNOf7WJN9O8oMkx5Lc2P9UJUlrGRn0JGPAHcDVwKXA9UkuXTHss8APq+oyYDfwj0nO73mukqQ1dLlCvwI4XlVPV9XLwH3AtSvGFPCWJAHeDPwMON3rTCVJa+oS9J3As8u2Tw32LXc78D7gOeBx4PNV9erKB0oym2QhycLi4uLrnLIkaZguQc+QfbVi+2PAo8DvAn8I3J7kt8/4Q1VzVTVTVTOTk5PrnKokaS1dgn4KuGjZ9oUsXYkvdyNwfy05DvwYeG8/U5QkddEl6A8DlyS5ePBE53XAAyvGnAQ+CpDkXcB7gKf7nKgkaW3jowZU1ekktwAPAmPAXVV1LMlNg+MHgS8Bdyd5nKVbNLdW1QsbOG9J0gojgw5QVYeAQyv2HVz2/XPAX/Q7NUnSevhOUUlqhEGXpEYYdElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEYYdElqhEGXpEZ0CnqSPUmeSnI8yW2rjNmd5NEkx5Ic6XeakqRRxkcNSDIG3AH8OXAKeDjJA1X1w2Vj3gYcAPZU1ckk79yg+UqSVtHlCv0K4HhVPV1VLwP3AdeuGPNJ4P6qOglQVc/3O01J0ihdgr4TeHbZ9qnBvuXeDbw9yeEkjyS5YdgDJZlNspBkYXFx8fXNWJI0VJegZ8i+WrE9DlwO/CXwMeBvk7z7jD9UNVdVM1U1Mzk5ue7JSpJWN/IeOktX5Bct274QeG7ImBeq6kXgxSTfBy4DftTLLCVJI3W5Qn8YuCTJxUnOB64DHlgx5lvAh5OMJ5kAPgQ82e9UJUlrGXmFXlWnk9wCPAiMAXdV1bEkNw2OH6yqJ5N8D3gMeBW4s6qe2MiJS5Jeq8stF6rqEHBoxb6DK7a/Anylv6lJktbDd4pKUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiM6BT3JniRPJTme5LY1xn0wyStJPtHfFCVJXYwMepIx4A7gauBS4Pokl64y7svAg31PUpI0Wpcr9CuA41X1dFW9DNwHXDtk3OeAbwDP9zg/SVJHXYK+E3h22fapwb5fS7IT+DhwcK0HSjKbZCHJwuLi4nrnKklaQ5egZ8i+WrH9VeDWqnplrQeqqrmqmqmqmcnJyY5TlCR1Md5hzCngomXbFwLPrRgzA9yXBGAHcE2S01X1zT4mKUkarUvQHwYuSXIx8D/AdcAnlw+oqot/9X2Su4HvGHNJ2lwjg15Vp5PcwtKrV8aAu6rqWJKbBsfXvG8uSdocXa7QqapDwKEV+4aGvKr+6uynJUlaL98pKkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1IhOQU+yJ8lTSY4nuW3I8b1JHht8PZTksv6nKklay8igJxkD7gCuBi4Frk9y6YphPwauqqr3A18C5vqeqCRpbV2u0K8AjlfV01X1MnAfcO3yAVX1UFX9fLB5FLiw32lKkkbpEvSdwLPLtk8N9q3m08B3hx1IMptkIcnC4uJi91lKkkbqEvQM2VdDByYfYSnotw47XlVzVTVTVTOTk5PdZylJGmm8w5hTwEXLti8Enls5KMn7gTuBq6vqp/1MT5LUVZcr9IeBS5JcnOR84DrggeUDkuwC7gc+VVU/6n+akqRRRl6hV9XpJLcADwJjwF1VdSzJTYPjB4EvAO8ADiQBOF1VMxs3bUnSSl1uuVBVh4BDK/YdXPb9Z4DP9Ds1SdJ6+E5RSWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWpEp6An2ZPkqSTHk9w25HiSfG1w/LEkH+h/qnDzzTdz5MgRjhw5wvj4ODfffPNGnEaSNsT8/DxHjx7lyJEjTE9PMz8/3+8JqmrNL2AM+G/g94DzgR8Al64Ycw3wXSDAlcB/jnrcyy+/vNZj3759BZzxtW/fvnU9jiRthXvuuacmJiZe06+JiYm655571vU4wEKt0tUsHV9dkj8C/q6qPjbY/uvBPwR/v2zM14HDVXXvYPspYHdV/WS1x52ZmamFhYXO//CMj4/zyiuvDD121VVXdX4cSdoKR48e5aWXXjpj/9TUFM8880znx0nySFXNDDvW5ZbLTuDZZdunBvvWO4Yks0kWkiwsLi52OPVvrBZzSdoOhsUc4OTJk72dY7zDmAzZt/KyvssYqmoOmIOlK/QO5/61sbGxoVEfGxvj8OHD63koSdp009PTnDhx4oz9u3bt6u0cXa7QTwEXLdu+EHjudYw5K7Ozs+vaL0nnkv379zMxMfGafRMTE+zfv7+3c3QJ+sPAJUkuTnI+cB3wwIoxDwA3DF7tciXwi7Xun78eBw4cYN++fYyNjQFLV+b79u3jwIEDfZ5GkjbE3r17mZubY2pqiiRMTU0xNzfH3r17ezvHyCdFAZJcA3yVpVe83FVV+5PcBFBVB5MEuB3YA/wSuLGq1nzGc71PikqS1n5StMs9dKrqEHBoxb6Dy74v4LNnM0lJ0tnxnaKS1AiDLkmNMOiS1AiDLkmN6PQqlw05cbIInPkq+252AC/0OJ3twDW/MbjmN4azWfNUVU0OO7BlQT8bSRZWe9lOq1zzG4NrfmPYqDV7y0WSGmHQJakR2zXoc1s9gS3gmt8YXPMbw4aseVveQ5cknWm7XqFLklYw6JLUiHM66OfKh1Nvpg5r3jtY62NJHkpy2VbMs0+j1rxs3AeTvJLkE5s5v43QZc1Jdid5NMmxJEc2e4596/B3+61Jvp3kB4M137gV8+xLkruSPJ/kiVWO99+v1T5sdKu/2KAPpz6Xvzqu+Y+Btw++v/qNsOZl4/6Npd/6+Ymtnvcm/JzfBvwQ2DXYfudWz3sT1vw3wJcH308CPwPO3+q5n8Wa/xT4APDEKsd779e5fIV+BXC8qp6uqpeB+4BrV4y5FvjnWnIUeFuS39nsifZo5Jqr6qGq+vlg8yhLnw61nXX5OQN8DvgG8PxmTm6DdFnzJ4H7q+okQFVt93V3WXMBbxl8vsKbWQr66c2dZn+q6vssrWE1vffrXA56bx9OvY2sdz2fZulf+O1s5JqT7AQ+DhykDV1+zu8G3p7kcJJHktywabPbGF3WfDvwPpY+vvJx4PNV9ermTG9L9N6vTh9wsUV6+3DqbaTzepJ8hKWg/8mGzmjjdVnzV4Fbq+qVpYu3ba/LmseBy4GPAr8F/EeSo1X1o42e3AbpsuaPAY8Cfwb8PvCvSf69qv5vg+e2VXrv17kc9HPiw6k3Waf1JHk/cCdwdVX9dJPmtlG6rHkGuG8Q8x3ANUlOV9U3N2WG/ev6d/uFqnoReDHJ94HLgO0a9C5rvhH4h1q6wXw8yY+B9wL/tTlT3HS99+tcvuVyTnw49SYbueYku4D7gU9t46u15UauuaourqrpqpoG/gW4eRvHHLr93f4W8OEk40kmgA8BT27yPPvUZc0nWfo/EpK8C3gP8PSmznJz9d6vc/YKvapOJ7kFeJDffDj1seUfTs3SKx6uAY4z+HDqrZpvHzqu+QvAO4ADgyvW07WNf1NdxzU3pcuaq+rJJN8DHgNeBe6sqqEvf9sOOv6cvwTcneRxlm5H3FpV2/bX6ia5F9gN7EhyCvgi8CbYuH751n9JasS5fMtFkrQOBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakR/w9KwKxYOnauRwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "x=np.asarray([0.0,1.0,1.0,0.0,0.0])\n",
    "y=np.asarray([0.0,0.0,1.0,1.0,0.0])\n",
    "\n",
    "plt.plot(x,y,'k-o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0b4b6a96-6d8e-4151-b32f-7c1615cbc3f0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f1c1dc05ee0>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAvx0lEQVR4nO3deXxU9b3/8dcnCQECYUlYJAkkoLiAFMEIogiZY+tF1Fqv9oKiCAITK1av1yuC1EJB+ih4qz7qAhlWBdxwRaXU9k5YKgoEkSUg14CEJYgkYZUl2/f3Ryb8piGBCZnMmeXzfDzyyMxZ5nyOMnnP9yyfEWMMSimlIleU3QUopZSylwaBUkpFOA0CpZSKcBoESikV4TQIlFIqwsXYXcDFaNOmjUlLS7O7DKWUCikbNmwoNMa0rT49JIMgLS2NnJwcu8tQSqmQIiL5NU3XQ0NKKRXhNAiUUirCaRAopVSE0yBQSqkIp0GglFIRzi9BICLzRORHEdlay/wMETkqIt94fn7vNW+QiOwQkTwRGe+PepRSCmDx4sWkpaURFRVFWloaixcvtrukoOSvy0cXAK8Ab5xnmdXGmNu9J4hINPAq8AtgH7BeRJYaY7b5qS6lVIRavHgxTqeTkydPApCfn4/T6QRg2LBhdpYWdPwyIjDGrAKKL2LVPkCeMWaXMaYEeBu40x81KaUi28SJE8+GQJWTJ08yceJEmyoKXoE8R9BPRDaJyF9FpLtnWjKw12uZfZ5p5xARp4jkiEjOoUOHGrpWpVSIMsawcuVK8vNrvHeKPXv2BLii4BeoIPgaSDXG9AReBj7yTJcalq3xm3KMMS5jTLoxJr1t23PukFZKRbiioiJefPFFunXrRkZGBiI1/XmBTp06Bbiy4BeQIDDGHDPGnPA8XgY0EpE2VI4AOnotmgIUBKImpVToM8awevVq7r//fpKTk/mv//ovWrVqxfz585k7dy5xcXH/snxcXBzTpk2zqdrgFZBeQyJyCXDQGGNEpA+VAVQEHAG6ikhnYD8wFLgvEDUppUJXcXExb7zxBi6Xi+3bt9OiRQtGjx6N0+nkZz/72dnlYmNjGTduHAUFBbRu3ZqXX35ZTxTXwC9BICJvARlAGxHZB0wCGgEYY2YB9wC/EZEy4BQw1FR+WXKZiDwK/A2IBuYZY3L9UZNSKrwYY/jiiy9wuVwsWbKE06dP07dvX+bOncuQIUNo1qzZOesMGzaM++67j9TUVPr27ashUAu/BIEx5t4LzH+FystLa5q3DFjmjzqUUuHn8OHDLFy4kKysLLZt20Z8fDwjR44kMzOTnj17XnB9EcGyLD799FMqKiqIitL7aKvT/yJKqaBjjGHNmjU8+OCDJCUl8fjjj9OsWTPmzJnDgQMHeO2113wKgSqWZVFUVMSWLVsasOrQFZLfR6CUCk9Hjhxh4cKFuFwutm7dSnx8PCNGjMDpdNKrV6+Lfl2HwwGA2+2uU4BECh0RKKVsZYzhyy+/ZOTIkSQlJfHYY4/RpEkTXC4XBQUFzJw5s14hANCxY0cuu+wysrOz/VR1eNERgVLKFkePHmXRokVkZWWxZcsWmjdvzvDhw3E6nfTu3dvv27Msi7fffpuysjJiYvRPnzcdESilAsYYw9q1a3nooYfo0KEDjz76KI0aNSIrK4uCggJmzZrVICEAlUFw7Ngxvv766wZ5/VCmsaiUanBHjx5l8eLFZGVlsXnzZpo1a8b9999PZmYm1157bUBqyMjIACA7O5s+ffoEZJuhQkcESqkGYYxh3bp1jBo1iqSkJMaOHUt0dDSzZs2ioKAAl8sVsBAAaN++Pd27d8ftdgdsm6FCRwRKKb86duwYb775JllZWXzzzTfExcVx7733kpmZSXp6eq09gALBsizmzp1LSUkJsbGxttURbHREoJTyi5ycHMaMGUNSUhK/+c1vMMbw2muvceDAAebMmcN1111nawhA5WWkJ0+eZN26dbbWEWx0RKCUumjHjx/nzTffxOVy8fXXXxMXF8fQoUPJzMwMij/81Q0cOBARwe12079/f7vLCRo6IlBK1dmGDRvIzMwkKSmJhx9+mLKyMl599VUKCgqYO3cuffr0CboQAEhISKBXr156nqAaHREopXxy/Phx3nrrLVwuFxs2bKBp06YMHToUp9NJ3759g/IPf00cDgcvv/wyp06domnTpnaXExR0RKCUOq+NGzfy8MMPk5SURGZmJmfOnOHll1+moKCAefPmcf3114dMCEDlCeOSkhLWrFljdylBQ0cESqlznDhxgrfffpusrCxycnJo0qQJQ4YMITMzM+T+8Fd30003ER0djdvt5uabb7a7nKCgQaCUOuubb74hKyuLxYsXc/z4cbp3785f/vIX7r//flq3bm13eX4RHx9Pnz599DyBFw0CpSLcTz/9xNtvv43L5WLdunU0adKE//iP/8DpdHLDDTeE9Kf/2jgcDqZPn87x48eJj4+3uxzb6TkCpSLUpk2bGDt2LElJSYwePZrjx4/z0ksvsX//fl5//XVuvPHGsAwBqDxPUF5ezurVq+0uJSjoiECpCHLy5EneeecdsrKyWLt2LY0bN+bXv/41mZmZYf2Hv7obbriB2NhY3G43gwcPtrsc22kQKBUBtmzZQlZWFosWLeLo0aNceeWVvPjiiwwfPpyEhAS7ywu4pk2b0q9fP/1+Ag8NAqXC1MmTJ3n33XdxuVx8+eWXNG7cmHvuuQen08lNN90UMZ/+a2NZFpMnT6a4uDgiw9CbX84RiMg8EflRRLbWMn+YiGz2/KwRkZ5e83aLyBYR+UZEcvxRj1KRLDc3l8cee4zk5GRGjhxJcXExf/7zn9m/fz+LFi1iwIABER8CUBkExhhWrlxpdym289eIYAHwCvBGLfO/BwYaYw6LyK2AC+jrNd9hjCn0Uy1KRZxTp06xZMkSsrKyWLNmDbGxsdx9991kZmbqH/5a9OnTh7i4OLKzs7nrrrvsLsdWfgkCY8wqEUk7z3zvW/i+AlL8sV2lIt22bdvIysrijTfe4MiRI1x++eX8z//8Dw8++CBt2rSxu7ygFhsbS//+/fV+Auy5fHQU8Fev5wb4XEQ2iIiztpVExCkiOSKSc+jQoQYvUqlgderUKRYuXMhNN91E9+7dmTlzJoMGDSI7O5tvv/2WJ598UkPAR5ZlkZuby8GDB+0uxVYBDQIRcVAZBE97Tb7RGNMbuBUYKyIDalrXGOMyxqQbY9Lbtm0bgGqVCi7bt2/niSeeIDk5meHDh/PDDz8wY8YM9u/fz1tvvUVGRoYeAqojy7IAWLFihb2F2CxgVw2JyM+AOcCtxpiiqunGmALP7x9F5EOgD7AqUHUpFcxOnz7N+++/T1ZWFqtXr6ZRo0bcddddZGZmkpGRQVSU3hNaH7169aJFixa43W6GDBlidzm2CUgQiEgn4APgAWPM/3lNbwZEGWOOex7fAkwJRE1KBbNvv/0Wl8vF66+/TnFxMZdeeinTp09nxIgRtGvXzu7ywkZMTAwDBw6M+PMEfgkCEXkLyADaiMg+YBLQCMAYMwv4PZAIvOYZupYZY9KB9sCHnmkxwJvGmOX+qEmpUHPmzBnef/99XC4XK1euJCYm5uynf4fDoZ/+G4hlWXzyySfs3buXjh072l2OLfx11dC9F5g/Ghhdw/RdQM9z11AqcuzYsYPZs2ezYMECioqK6NKlC3/6058YMWIE7du3t7u8sOdwOADIzs5m+PDhNldjD72zWCkbnDlzhg8//JCsrCxWrFhBTEwMd955J5mZmdx888366T+AevToQWJiIm63W4NAKdXwvvvuO1wuFwsWLKCwsJDOnTvzxz/+kZEjR3LJJZfYXV5EioqKwuFw4Ha7McZE5JVXGgRKNbCSkhI+/PBDXC4Xbreb6Ojos5/+f/7zn+un/yDgcDh477332LVrF5deeqnd5QScBoFSDSQvL4/Zs2czf/58Dh06RFpaGtOmTWPkyJF06NDB7vKUl6r7CdxutwaBUqp+SkpK+Oijj3C5XPzv//4v0dHR/PKXv8TpdHLLLbfop/8gdcUVV9ChQwfcbjdjxoyxu5yA0yBQyg927tx59tP/jz/+SKdOnZg6dSoPPfQQSUlJdpenLkBEsCyLf/zjHxF5nkCDQKmLVFpayscff0xWVhb/+Mc/iI6O5vbbbyczM5NbbrmF6Ohou0tUdeBwOFi8eDHbt2+nW7dudpcTUBoEStXRrl27mDNnDvPmzePgwYN07NiRKVOm8NBDD5GcnGx3eeoieZ8n0CBQSp2jtLSUTz75hKysLD7//HOioqK47bbbyMzMZNCgQfrpPwx07tyZtLQ03G43jz76qN3lBJQGgVLnsXv3bmbPns28efP44YcfSElJYfLkyYwaNYqUFP1ajXDjcDj46KOPqKioiKgT+5Gzp0r5qLS0lA8//JBBgwadbfeQnp7OJ598wu7du5k0aZKGQJiyLIvDhw+zadMmu0sJKB0RKOWRn59/9tP/gQMHSE5O5ve//z2jRo2K2GZkkaaq75Db7aZXr142VxM4GgQqopWVlfHpp5/icrlYvryy8e3gwYNxOp0MHjyYmBh9i0SS5ORkLr/8crKzs3nyySftLidg9F+5ikh79uxhzpw5zJ07l4KCApKSkvjd737H6NGj6dSpk93lKRtZlsWiRYsoLS2lUaNGdpcTEHqOQEWMsrIyli5dym233UZaWhrPPfccPXv25KOPPiI/P58pU6ZoCCgsy+LEiRNs2LDB7lICRkcEKuzt3bv37Kf//fv306FDByZOnMioUaNIS0uzuzwVZDIyMoDK7ye4/vrr7S0mQHREoMJSeXk5n3zyCXfccQdpaWlMnTqVHj168MEHH5Cfn8/UqVM1BFSN2rZtS48ePSLq6yt1RKDCyr59+5g7dy5z5sxh3759XHLJJYwfP57Ro0fTuXNnu8tTIcKyLLKysjhz5gyNGze2u5wGpyMCFfLKy8v57LPP+OUvf0lqaiqTJ0+mW7duvP/+++zZs4dp06ZpCKg6sSyL06dPs3btWrtLCQgdEaiQtX///rOf/vfu3Uv79u15+umnGT16NF26dLG7PBXCBgwYQFRUFG63mwEDBthdToPzy4hAROaJyI8isrWW+SIifxGRPBHZLCK9veYNEpEdnnnj/VGPCl/l5eUsW7aMX/3qV6SmpjJp0iSuvPJKlixZwp49e/jjH/+oIaDqrVWrVvTu3TtizhP469DQAmDQeebfCnT1/DiBmQAiEg286pnfDbhXRCKr7Z/ySUFBAc899xxdunThtttu48svv+S///u/ycvL4/PPP+eee+4hNjbW7jJVGLEsi6+++oqTJ0/aXUqD80sQGGNWAcXnWeRO4A1T6SuglYh0APoAecaYXcaYEuBtz7JKUV5ezvLly7nrrrvo1KkTzz77LF27duXdd99l7969/OlPf4rIrxVUgeFwOCgtLeWLL76wu5QGF6hzBMnAXq/n+zzTapret6YXEBEnlaMJveknzB04cIB58+Yxe/Zs8vPzadu2LU8++SRjxozhsssus7s8FSH69+9PTEwMbrebX/ziF3aX06ACFQQ1fe+bOc/0cyca4wJcAOnp6TUuo0JXRUUFf//733G5XCxdupSysjIsy2L69On86le/iohL+FRwad68OX379o2I8wSBCoJ9gHf7xhSgAIitZbqKED/88APz589n9uzZfP/997Rp04b//M//ZMyYMVx++eV2l6cinGVZTJs2jaNHj9KyZUu7y2kwgbqPYCkw3HP10PXAUWPMAWA90FVEOotILDDUs6wKY1Wf/u+55x46duzIM888Q2pqKm+99Rb79u3j+eef1xBQQcHhcFBRUcHq1avtLqVB+WVEICJvARlAGxHZB0wCGgEYY2YBy4DBQB5wEhjpmVcmIo8CfwOigXnGmFx/1KSCz8GDB89++t+1axeJiYk8/vjjjBkzhiuuuMLu8pQ6R79+/WjcuDFut5vbb7/d7nIajF+CwBhz7wXmG2BsLfOWURkUKgxVVFTgdrvJysrio48+oqysjIEDB/Lcc89x11130aRJE7tLVKpWTZo04cYbbyQ7O9vuUhqU3lmsGsSPP/7IggULcLlc7Ny5k4SEBH7729/idDq58sor7S5PKZ85HA6effZZioqKSExMtLucBqG9hpTfVH36HzJkCCkpKTz99NMkJSWxaNEi9u/fzwsvvKAhoEKOZVkArFixwt5CGpCOCFS9HTp06Oyn/7y8PFq3bs3YsWNxOp1cddVVdpenVL1cd911NGvWDLfbzd133213OQ1Cg0BdFGMMK1aswOVy8cEHH1BSUkL//v2ZNGkSd999N02bNrW7RKX8olGjRtx0001hfZ5Ag0DVSWFh4dlP/9999x2tWrXiN7/5DWPGjKF79+52l6dUg7Asi3HjxnHgwAE6dOhgdzl+p+cI1AUZY1i5ciX33XcfycnJPPXUU7Rr147XX3+dgoICXnrpJQ0BFdaqzhOE66hARwSqVkVFRbz++uu4XC527NhBy5YtyczMxOl0cvXVV9tdnlIBc80119CqVSuys7O577777C7H7zQI1L8wxrB69WqysrJ47733KCkpoV+/fixYsIBf//rXxMXF2V2iUgEXHR3NwIEDw7bvkAaBAqC4uJg33ngDl8vF9u3badGiBWPGjCEzM5MePXrYXZ5StrMsi48//pj8/HxSU1PtLsev9BxBBDPG8M9//pMHHniApKQknnjiCVq2bMm8efMoKCjglVde0RBQyiOczxPoiCACHT58+Oyn/23bttGiRQtGjRqF0+mkZ8+edpenVFDq3r07bdu2xe12M2LECLvL8SsNgghhjGHNmjVkZWWxZMkSTp8+TZ8+fZg7dy5DhgyhWbNmdpeoVFATERwOB263G2MMIjV9nUpo0iAIc4cPH2bhwoW4XC5yc3OJj49n5MiROJ1OrrnmGrvLUyqkWJbFu+++S15eHl27drW7HL/RIAhDxhi++uorsrKyeOeddzh9+jTXXXcds2fPZujQoTRv3tzuEpUKSQ6HAwC3261BoILTkSNHWLRoES6Xiy1bttC8eXMefPBBnE4nvXv3trs8pUJe165dSU5Oxu12k5mZaXc5fqNBEIIWL17MxIkT2bNnDx07duShhx5i9+7dvPPOO5w6dYprr70Wl8vF0KFDiY+Pt7tcpcKGiGBZFsuXLw+r8wQaBCFm8eLFOJ1OTp48CcCePXuYPHkyjRs3Pvvp/9prr7W5SqXCl8PhYOHCheTm5obNHfYaBCFm4sSJZ0PAW7t27cjKyrKhIqUiS9X9BG63O2yCQG8oCzF79uypcfq+ffsCXIlSkSk1NZUuXbqEVbsJDYIQk5KSUuP0Tp06BbgSpSKXZVmsXLmS8vJyu0vxC78EgYgMEpEdIpInIuNrmP+UiHzj+dkqIuUikuCZt1tEtnjm5fijnnDWq1evc6bFxcUxbdo0G6pRKjI5HA6OHDnCN998Y3cpflHvIBCRaOBV4FagG3CviHTzXsYY87wx5hpjzDXABGClMabYaxGHZ356fesJZ/v27ePzzz+nX79+pKamIiKkpqbicrkYNmyY3eUpFTG87ycIB/4YEfQB8owxu4wxJcDbwJ3nWf5e4C0/bDfiTJo0iYqKChYvXszu3bupqKhg9+7dGgJKBViHDh246qqrwqYBnT+CIBnY6/V8n2faOUQkDhgEvO812QCfi8gGEXHWthERcYpIjojkHDp0yA9lh5atW7eyYMECxo4dS+fOne0uR6mI53A4WLVqFaWlpXaXUm/+CIKa7qgwtSx7B/BFtcNCNxpjelN5aGmsiAyoaUVjjMsYk26MSW/btm39Kg5B48ePJz4+nokTJ9pdilKKyhPGP/30E+vXr7e7lHrzRxDsAzp6PU8BCmpZdijVDgsZYwo8v38EPqTyUJPysmLFCj777DOeeeYZEhMT7S5HKQVkZGQA4fH9BP4IgvVAVxHpLCKxVP6xX1p9IRFpCQwEPvaa1kxE4qseA7cAW/1QU9ioqKjgqaeeIiUlhd/+9rd2l6OU8khMTKRnz55hccK43ncWG2PKRORR4G9ANDDPGJMrIg975s/yLHoX8Lkx5iev1dsDH3r6dcQAbxpjlte3pnCyZMkScnJyWLBgAU2bNrW7HKWUF8uyeO211zh9+jRNmjSxu5yLJsbUdjg/eKWnp5ucnPC/5aCkpISrrrqKZs2asXHjRqKjo+0uSSnl5dNPP+WOO+7A7XafvaQ0mInIhpou09c7i4PYrFmz2LVrFzNmzNAQUCoIDRgwgOjo6JA/T6BBEKSOHj3KlClTuPnmm/m3f/s3u8tRStWgRYsWXHvttSF/nkCDIEjNmDGDoqIipk+fHjY9z5UKR5ZlsXbtWk6cOGF3KRdNgyAI7d+/nxdffJH77rtPv1tAqSBnWRZlZWV88cUXdpdy0TQIgtCkSZMoLy/nueees7sUpdQF3HjjjTRq1CikDw9pEASZrVu3Mn/+fG0loVSIiIuL4/rrr9cgUP6jrSSUCj2WZfH1119z5MgRu0u5KBoEQaSqlcSECRO0lYRSIcThcFBRUcGqVavsLuWiaBAECWMM48aNIyUlhccee8zucpRSdXD99dfTpEmTkD08pF9eHySWLFnC+vXrmT9/vraSUCrENG7cmP79+4dsEOiIIAiUlJQwYcIEevTowQMPPGB3OUqpi2BZFlu2bCEUvy9FgyAIaCsJpUJfVa+hFStW2FvIRdAgsNnRo0eZOnUqlmVpKwmlQlh6ejrx8fEheXhIg8BmM2bMoLCwkBkzZmgrCaVCWExMDAMGDAjJBnQaBDaqaiVx7733aisJpcKAw+Fgx44d7N+/3+5S6kSDwEaTJk2irKyMadOm2V2KUsoPLMsCQu/rKzUIbJKbm8v8+fN59NFHtZWEUmGiZ8+etG7dWoNA+UZbSSgVfqKiosjIyAi5E8YaBDZYuXIln376qbaSUCoMWZbF7t27+f777+0uxWcaBAFmjOGpp57SVhJKhalQPE/glyAQkUEiskNE8kRkfA3zM0TkqIh84/n5va/rhpuqVhJTp07VVhJKhaGrrrqK9u3bh9ThoXr3GhKRaOBV4BfAPmC9iCw1xmyrtuhqY8ztF7luWNBWEkqFPxHB4XDgdrsxxoTE/UH+GBH0AfKMMbuMMSXA28CdAVg35GRlZbFr1y6mT5+urSSUCmOWZXHgwAF27Nhhdyk+8UcQJAN7vZ7v80yrrp+IbBKRv4pI9zquG/KOHTvGlClTsCyLQYMG2V2OUqoBhdp5An8EQU3jHlPt+ddAqjGmJ/Ay8FEd1q1cUMQpIjkikhOK3f20lYRSkaNLly507NgxZM4T+CMI9gEdvZ6nAAXeCxhjjhljTngeLwMaiUgbX9b1eg2XMSbdGJPetm1bP5QdOPv37+eFF17QVhJKRQgRwbIssrOzqaiosLucC/JHEKwHuopIZxGJBYYCS70XEJFLxPMxWET6eLZb5Mu64WDy5MnaSkKpCGNZFkVFRWzdutXuUi6o3kFgjCkDHgX+BmwH3jXG5IrIwyLysGexe4CtIrIJ+Asw1FSqcd361hRMcnNzmTdvHmPHjtVWEkpFkKrvJwiFw0NiTI2H5INaenq6ycnJsbsMn9xxxx2sXr2anTt36l3ESkWYrl27ctVVV7F0aXAc6BCRDcaY9OrT9c7iBlTVSmL8+PEaAkpFIMuyWLlyJWVlZXaXcl4aBA3EGMO4ceNISUnh8ccft7scpZQNLMvi2LFjbNy40e5SzkuDoIG89957rFu3TltJKBXBMjIygOA/T6DnCBpASUkJ3bp1Iy4ujo0bN+pdxEpFsKuvvpqUlBSWL19udym1niOod68hda6srCx27tzJsmXLNASUinCWZTF37lxKSkqIjY21u5wa6aEhP6tqJeFwOLSVhFIKh8PByZMnWbdund2l1EqDwM+0lYRSytvAgQMRkaA+T6BB4EcFBQVnW0mkp59zGE4pFYESEhLo1atXUDeg0yDwo0mTJmkrCaXUORwOB2vWrOHUqVN2l1IjDQI/2bZtm7aSUErVyLIsSkpKWLNmjd2l1EiDwE/Gjx9P8+bNmThxot2lKKWCzE033UR0dHTQHh7SIPCDVatW8cknnzBhwgTatGljdzlKqSATHx9Pnz59gvaEsQZBPRljeOqpp7SVhFLqvBwOB+vWreP48eN2l3IODYJ6qmolMWXKFG0loZSqlWVZlJeX889//tPuUs6hQVAPJSUlTJgwgauvvprhw4fbXY5SKojdcMMNxMbGBuXhIW0xUQ8ul4udO3fy2WefaSsJpdR5NW3alH79+gVlEOiI4CIdO3aMP/zhDzgcDm699Va7y1FKhQDLsti4cSPFxcV2l/IvNAgu0vPPP6+tJJRSdWJZFsYYVq1aZXcp/0KD4CIUFBTw5z//maFDh2orCaWUz/r06UNcXFzQHR7SILgI2kpCKXUxYmNj6d+/f3gGgYgMEpEdIpInIuNrmD9MRDZ7ftaISE+vebtFZIuIfCMiwfttMx5VrSQeeeQRunTpYnc5SqkQY1kWubm5HDx40O5Szqp3EIhINPAqcCvQDbhXRLpVW+x7YKAx5mfAVMBVbb7DGHNNTd+cE2yqWkn87ne/s7sUpVQIsiwLgBUrVthbiBd/jAj6AHnGmF3GmBLgbeBO7wWMMWuMMYc9T78CUvyw3YDTVhJKqfrq1asXLVq0CKrDQ/4IgmRgr9fzfZ5ptRkF/NXruQE+F5ENIuL0Qz0NoqqVRHJysraSUEpdtJiYGAYOHBhUDej8EQQ1XTtpalxQxEFlEDztNflGY0xvKg8tjRWRAbWs6xSRHBHJOXToUH1rrrOqVhJTp07VVhJKqXqxLIvvvvuOvXv3XnjhAPBHEOwDOno9TwEKqi8kIj8D5gB3GmOKqqYbYwo8v38EPqTyUNM5jDEuY0y6MSa9bdu2fijbdyUlJTzzzDPaSkIp5RcOhwMgaEYF/giC9UBXEeksIrHAUGCp9wIi0gn4AHjAGPN/XtObiUh81WPgFmCrH2ryK5fLRV5eHtOnT9dWEkqpeuvRoweJiYlBEwT17jVkjCkTkUeBvwHRwDxjTK6IPOyZPwv4PZAIvOa5C7fMc4VQe+BDz7QY4E1jzPL61uRPx44dY8qUKdpKQinlN1FRUTgcDtxuN8YY27sT+KXpnDFmGbCs2rRZXo9HA6NrWG8X0LP69GDy/PPPc+jQIW0loZTyK4fDwXvvvceuXbu49NJLba1F7yw+D20loZRqKFX3EwTDZaQaBOcxefJkbSWhlGoQV1xxBR06dAiK8wQaBLXYtm0bc+fO1VYSSqkGISJYlnX2PIGdNAhqMWHCBG0loZRqUA6Hg4MHD7J9+3Zb69AgqMGqVatYunQp48eP11YSSqkGU3WewO7DQxoE1RhjGDdunLaSUEo1uM6dO5OWlmb7CWP9zuJq3n//fdauXcvcuXOJi4uzuxylVJhzOBx8/PHHVFRUEBVlz2dzHRF4KS0tZcKECXTv3p0HH3zQ7nKUUhHAsiyKi4vZvHmzbTVoEHipaiUxY8YMbSWhlAqIqr5Ddh4e0iDwOHbsGH/4wx/IyMjQVhJKqYBJTk7m8ssv1yAIBtpKQillF8uyWLVqFWVlZbZsX4OAylYSL7zwAkOGDOG6666zuxylVISxLIvjx4+zYcMGW7avQUBlK4nS0lJtJaGUskVGRgZg33mCiA+C7du3n20lYXcHQKVUZGrbti09evTQILDL+PHjtZWEUsp2lmXxxRdfcObMmYBvO6KDYPXq1dpKQikVFCzL4tSpU6xduzbg247YIDDG8NRTT2krCaVUUBgwYABRUVG2HB6K2CCoaiUxZcoUbSWhlLJdq1at6N27ty0N6CIyCLSVhFIqGFmWxZdffsnJkycDut2IDIKqVhLTp0/XVhJKqaDhcDgoLS3liy++COh2Iy4IvFtJDB482O5ylFLqrP79+xMTExPww0N+CQIRGSQiO0QkT0TG1zBfROQvnvmbRaS3r+v6yyOPPEJMTAwtW7bk0KFDJCYmaisJpVRQad68OX379g34CeN6B4GIRAOvArcC3YB7RaRbtcVuBbp6fpzAzDqsW2+PPPIIM2fOpLy8/Oy0999/n0ceecTfm1JKqXqxLIv169dz9OjRgG3TH19M0wfIM8bsAhCRt4E7gW1ey9wJvGEqv6H5KxFpJSIdgDQf1q03l8tV4/SZM2eyZcsWEhMTadOmzTm/vR+3atVKzycopRqcw+Fg6tSprF69mttvvz0g2/RHECQDe72e7wP6+rBMso/rAiAiTipHE3Tq1KlOBXqPBKpr1KgRu3btYv369RQWFlJSUlLjciJCQkLCecOi+u+EhAQND6VUnfTr14/GjRuTnZ0dUkFQ04F24+MyvqxbOdEYF+ACSE9Pr3GZ2kRHR9cYBtHR0f9yLM4Yw08//URhYSGFhYUUFRXV+nvPnj18/fXXFBYW1npLuIjQqlWr84ZF9WkJCQk0atSoLrunlAojTZo04cYbbwzoeQJ/BME+oKPX8xSgwMdlYn1Yt96cTiczZ86scbo3EaF58+Y0b96ctLQ0n17bGMPJkyfPCYuaAmT//v1s2rSJoqKi814n3LJlS59HHVUjj9jY2Dr9N1FKBS+Hw8Gzzz5LUVERiYmJDb49fwTBeqCriHQG9gNDgfuqLbMUeNRzDqAvcNQYc0BEDvmwbr299tprQOW5gvLycqKjo3E6nWen14eI0KxZM5o1a1anQ1ZV4XG+UUdhYSE//PADubm5FBYW8tNPP9X6ei1atPB51FH1u3HjxvXef6WU/1mWxbPPPsvKlSv593//9wbfnlSev63ni4gMBl4CooF5xphpIvIwgDFmllRep/kKMAg4CYw0xuTUtu6Ftpeenm5ycnLqXXeoOX369L+Ex4UOXxUWFnL8+PFaX6958+bnHW3UFCRNmjQJ4B4rFZlKS0tp3bo1I0aM4JVXXvHb64rIBmNM+jnT/REEgRapQXAxzpw5Q3FxcY0hUVuAnO+ytbi4OJ8PW1X91l5OStXdrbfeSn5+Ptu2+e8iytqCwB+HhlQQa9y4MR06dKBDhw4+r1NaWno2PC406vj+++8pLCzkyJEjtb5e06ZN6zTqSExMpFmzZnrDn4polmUxbtw4fvjhBy655JIG3ZYGgTpHo0aNaN++Pe3bt/d5nbKyMoqLi30adWzcuJHCwkIOHz5MbSPSxo0b12nU0aZNG5o3b67hocKGZVkAZGdnc++99zbotjQIlF/ExMTQrl072rVr5/M65eXlHD582KdzHZs3b6aoqIji4mIqKipqfL1GjRrVadTRpk0bWrRooeGhgtI111xDq1atcLvdGgQqfEVHR5/9A+2riooKjhw54tO5jm3btp09uV7bTYUxMTEXHG1UD5KWLVsSFRVx/RpVgEVHRzNw4MCANKDTIFAhJSoqioSEBBISEnxep6KigqNHj15w1FFYWMiOHTvOTisrK6vx9aKjo0lISKjTSfPWrVtreKg6syyLjz/+mPz8fFJTUxtsOxoEKuxFRUXRunVrWrduzWWXXebTOsYYjh075tOlujt37mTt2rUUFhZSWlpaaw01tSg53+Gr1q1ba4uSCOd9nmDEiBENth29fFQpPzHGcOLEiQuOOqpPO1+LktatW9fppHlCQgIxMfr5LlwYY2jfvj2DBg3ijTfeqPfr6eWjSjUwESE+Pp74+Hg6d+7s0zpVLUp8uVR37969Z6+4On36dK2vWVt/q9oCJDExUftbBSkRweFwkJ2djTGmwS5s0CBQykbeLUrqcgy4pv5WNY06Dhw4wJYtWygqKrpgixJf+lp5T9P+VoFhWRbvvvsueXl5dO3atUG2oUGgVAiKi4sjLi6Ojh07Xnhhj1OnTvnU3+rHH39k+/btFBYWcuLEiVpfLz4+vk5t2RMTE7VFyUVwOBwAuN1uDQKlVP00bdqUlJQUUlJSfF7nzJkzNYZHTQHy3XffUVhYyLFjx2p9vWbNmtX5Xo+mTZv6Y/dDVteuXUlOTiY7O5vMzMwG2YYGgVKqVo0bNyYpKYmkpCSf1ykpKam1v1X1INm5cydFRUUXbFFSl7bsiYmJxMXFhc2NgiKCZVksX768wc4TaBAopfwqNjaWSy65pE79capalPhyo2B+fv7ZFiW1adKkSZ3bsgdzixKHw8HChQvJzc3l6quv9vvraxAopWx3MS1KysrKOHz4sE8nzTdt2kRhYSHFxcW19reKjY2t06ijTZs2xMfHByQ8vO8naIgg0PsIlFIRo7y8nCNHjtTpOz2KiorO29+q6hJcXw9ftWzZ8qLCo127dpw4cYLTp0/TqVMnpk2bxrBhw+r0GnofgVIq4kVHR5/9w3355Zf7tE5VixJfbhD89ttvzz4/X3+r6i1KLhQgy5Yto7i4+Oxr5ufnn/2q3bqGQU10RKCUUn5W1aLE11FH1ePaWpTUJjU1ld27d/u8vI4IlFIqQESEli1b0rJlSy699FKf1jHGcPz48RpD4oknnqhxnT179vilXg0CpZQKAiJCixYtaNGixTktSl566SXy8/PPWadTp05+2bb2xVVKqSA3bdq0c777Oy4ujmnTpvnl9esVBCKSICJ/F5HvPL9b17BMRxHJFpHtIpIrIo97zZssIvtF5BvPz+D61KOUUuFo2LBhuFwuUlNTERFSU1NxuVx+OVEM9TxZLCIzgGJjzJ9EZDzQ2hjzdLVlOgAdjDFfi0g8sAH4lTFmm4hMBk4YY/6nLtvVk8VKKVV3tZ0sru+hoTuB1z2PXwd+VX0BY8wBY8zXnsfHge1Acj23q5RSyk/qGwTtjTEHoPIPPnDe2wJFJA3oBaz1mvyoiGwWkXk1HVryWtcpIjkiknPo0KF6lq2UUqrKBYNARP4hIltr+LmzLhsSkebA+8B/GmOq2hPOBC4FrgEOAH+ubX1jjMsYk26MSW/btm1dNq2UUuo8Lnj5qDHm57XNE5GDItLBGHPAcy7gx1qWa0RlCCw2xnzg9doHvZaZDXxal+KVUkrVX30PDS0FHvQ8fhD4uPoCUtlUYy6w3RjzQrV5Hbye3gVsrWc9Siml6qi+Vw0lAu8CnYA9wK+NMcUikgTMMcYMFpH+wGpgC1DVuekZY8wyEVlI5WEhA+wGMqvOOVxgu4eAc++u8E0boPAi1w02ui/BJ1z2A3RfglV99iXVGHPOsfWQ7DVUHyKSU9PlU6FI9yX4hMt+gO5LsGqIfdE7i5VSKsJpECilVISLxCBw2V2AH+m+BJ9w2Q/QfQlWft+XiDtHoJRS6l9F4ohAKaWUFw0CpZSKcGEbBCIySER2iEiepzNq9fkiIn/xzN8sIr3tqNMXPuzLMM8+bBaRNSLS0446L+RC++G13HUiUi4i9wSyvrrwZV9EJMPTXj1XRFYGukZf+fDvq6WIfCIimzz7MtKOOi/E06/sRxGp8cbUEHvPX2hf/PueN8aE3Q8QDewEugCxwCagW7VlBgN/BQS4Hlhrd9312JcbqGwBDnBrMO6LL/vhtZwbWAbcY3fd9fh/0grYBnTyPG9nd9312JdngOmex22BYiDW7tpr2JcBQG9gay3zQ+I97+O++PU9H64jgj5AnjFmlzGmBHibypbZ3u4E3jCVvgJaVWt5ESwuuC/GmDXGmMOep18BKQGu0Re+/D8B+C2Vfalq7FsVJHzZl/uAD4wxewCMMcG6P77siwHiPe1imlMZBGWBLfPCjDGrqKytNqHynr/gvvj7PR+uQZAM7PV6vo9zvwPBl2WCQV3rHEXlp55gc8H9EJFkKntOzQpgXRfDl/8nlwOtRWSFiGwQkeEBq65ufNmXV4CrgAIqW8U8boypIPSEynu+rur9ng/XL6+XGqZVv07Wl2WCgc91ioiDyn8U/Ru0oovjy368BDxtjCmv/PAZtHzZlxjgWuBmoCnwpYh8ZYz5v4Yuro582Zd/A74BLCrbxv9dRFab/99OPlSEynveZ/56z4drEOwDOno9T6Hy00xdlwkGPtUpIj8D5gC3GmOKAlRbXfiyH+nA254QaAMMFpEyY8xHAanQd77++yo0xvwE/CQiq4CeQLAFgS/7MhL4k6k8IJ0nIt8DVwLrAlOi34TKe94n/nzPh+uhofVAVxHpLCKxwFAqW2Z7WwoM91xJcD1w1PjQ+dQGF9wXEekEfAA8EISfOKtccD+MMZ2NMWnGmDTgPeCRIAwB8O3f18fATSISIyJxQF8qv6Y12PiyL3uoHNkgIu2BK4BdAa3SP0LlPX9B/n7Ph+WIwBhTJiKPAn+j8qqIecaYXBF52DN/FpVXpQwG8oCTVH7qCTo+7svvgUTgNc+n6TITZJ0WfdyPkODLvhhjtovIcmAzle3X5xhjgu77Nnz8/zIVWCAiW6g8vPK0MSboWjqLyFtABtBGRPYBk4BGEFrvefBpX/z6ntcWE0opFeHC9dCQUkopH2kQKKVUhNMgUEqpCKdBoJRSEU6DQCmlIpwGgVJKRTgNAqWUinD/D7bemXp6nrrWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x=np.asarray([0.0,1.2,1.0,0.2,0.0])\n",
    "y=np.asarray([0.0,-0.2,1.5,0.8,0.0])\n",
    "\n",
    "plt.plot(x,y,'k-o')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2a75bfa-2c0a-402f-abf7-d8cbfa8eddb9",
   "metadata": {},
   "source": [
    "For isoparametric mesh, one can define $\\xi,\\eta,\\zeta\\in\\{-1,1\\}$, then we can do the maping between $(x,y,z)$ and $(\\xi,\\eta,\\zeta)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d0c3c112-0cc6-446c-9fe2-7279048142e0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.         -0.57735027]\n",
      " [ 1.          0.57735027]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEKCAYAAAA1qaOTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAXHElEQVR4nO3df5BddZnn8ffTncSAgiGSjjGh0xmmB5MOAklPYJzhh2KAACXKjyniD7IgRYVdLKxdRRyrli3dKtDVBdxBWVQc1KlJ6YxIZgyyCKtRCUInYCCGQIwZ0iZAxCiSiEmHZ//oC9u0t9M35Htv37bfr6pb3eec77nP86Wb87nnntM3kZlIknSgWka6AUnSnwYDRZJUhIEiSSrCQJEkFWGgSJKKGDfSDYykww8/PDs6Oka6DUkaVVavXv2rzJwyeP2YDpSOjg56enpGug1JGlUi4t+rrfctL0lSEQaKJKkIA0WSVISBIkkqwkCRJBVhoEiSijBQJElFGCiSpCIMFElSEQaKJKkIA0WSVISBIkkqwkCRJBVhoEiSijBQJElFGCiSpCIMFElSEQaKJKkIA0WSVISBIkkqwkCRJBVhoEiSijBQJElFGCiSpCIMFElSEU0VKBFxRkRsiIiNEXF1le0REZ+rbF8bEfMGbW+NiIci4t8a17UkCZooUCKiFbgJWATMARZHxJxBwxYBnZXHZcAXBm2/Elhf51YlSVU0TaAAC4CNmbkpM3cDy4BzBo05B/hq9rsfmBQR0wAiYgZwFvClRjYtSerXTIEyHdgyYLm3sq7WMTcAVwEv7qtIRFwWET0R0bN9+/YDaliS9P81U6BElXVZy5iIOBt4JjNXD1ckM2/JzO7M7J4yZcqr6VOSVEUzBUovcMSA5RnA1hrH/DXwzojYTP9bZW+PiK/Xr1VJ0mDNFCgPAp0RMSsiJgAXAssHjVkOXFS52+sE4LeZuS0zP5aZMzKzo7LfvZn5voZ2L0lj3LiRbuAlmdkXEVcAdwGtwK2ZuS4illa23wysAM4ENgK7gItHql9J0itF5uDLFGNHd3d39vT0jHQbkjSqRMTqzOwevL6Z3vKSJI1iBookqQgDRZJUhIEiSSrCQJEkFWGgSJKKMFAkSUUYKJKkIgwUSVIRBookqQgDRZJUhIEiSSrCQJEkFWGgSJKKMFAkSUUYKJKkIgwUSVIRBookqQgDRZJUhIEiSSrCQJEkFWGgSJKKMFAkSUUYKJKkIgwUSVIRBookqQgDRZJUhIEiSSrCQJEkFWGgSJKKaKpAiYgzImJDRGyMiKurbI+I+Fxl+9qImFdZf0RE/N+IWB8R6yLiysZ3L0ljW9MESkS0AjcBi4A5wOKImDNo2CKgs/K4DPhCZX0f8F8yczZwAvCfquwrSaqjpgkUYAGwMTM3ZeZuYBlwzqAx5wBfzX73A5MiYlpmbsvMNQCZ+TtgPTC9kc1L0ljXTIEyHdgyYLmXPw6FYcdERAdwHPCT8i1KkobSTIESVdbl/oyJiNcB/wJ8KDOfq1ok4rKI6ImInu3bt7/qZiVJr9RMgdILHDFgeQawtdYxETGe/jD5x8z81lBFMvOWzOzOzO4pU6YUaVyS1FyB8iDQGRGzImICcCGwfNCY5cBFlbu9TgB+m5nbIiKALwPrM/N/NrZtSRLAuJFu4CWZ2RcRVwB3Aa3ArZm5LiKWVrbfDKwAzgQ2AruAiyu7/zXwfuCRiHi4su7vMnNFA6cgSWNaZA6+TDF2dHd3Z09Pz0i3IUmjSkSszszuweub6S0vSdIoZqBIkoowUCRJRRgokqQiDBRJUhEGiiSpCANFklSEgSJJKsJAkSQVYaBIkoowUCRJRRgokqQiDBRJUhEGiiSpCANFklSEgSJJKsJAkSQVYaBIkoowUCRJRRgokqQiDBRJUhEGiiSpCANFklSEgSJJKsJAkSQVYaBIkoowUCRJRRgoB+CSSy6hra2NuXPnNrTuli1beNvb3sbs2bPp6urixhtvbFjtF154gQULFnDMMcfQ1dXFNddc07DaAHv37uW4447j7LPPbmjdjo4Ojj76aI499li6u7sbVvc3v/kN559/Pm9+85uZPXs2q1atakjdDRs2cOyxx778OPTQQ7nhhhsaUvv666+nq6uLuXPnsnjxYl544YWG1L3xxhuZO3cuXV1dDZsrVD+O/PrXv2bhwoV0dnaycOFCduzY0bDa3/zmN+nq6qKlpYWenp79e8LMHLOP+fPn54H4wQ9+kKtXr86urq4Dep79tXXr1ly9enVmZj733HPZ2dmZ69ata0jtF198MX/3u99lZubu3btzwYIFuWrVqobUzsz87Gc/m4sXL86zzjqrYTUzM2fOnJnbt29vaM3MzIsuuii/+MUvZmbmH/7wh9yxY0fDe+jr68upU6fm5s2b616rt7c3Ozo6cteuXZmZecEFF+RXvvKVutd95JFHsqurK3fu3Jl79uzJU089NR9//PG6182sfhz5yEc+ktdee21mZl577bV51VVXNaz2z372s3zsscfy5JNPzgcffLDqfkBPVjmmeoZyAE466SQmT57c8LrTpk1j3rx5ABxyyCHMnj2bX/7ylw2pHRG87nWvA2DPnj3s2bOHiGhI7d7eXr7zne9w6aWXNqTeSHvuuedYuXIlH/jABwCYMGECkyZNangf99xzD0ceeSQzZ85sSL2+vj5+//vf09fXx65du3jTm95U95rr16/nhBNO4OCDD2bcuHGcfPLJ3H777XWvC9WPI3fccQdLliwBYMmSJXz7299uWO3Zs2dz1FFHvarnM1BGuc2bN/PQQw9x/PHHN6zm3r17OfbYY2lra2PhwoUNq/2hD32IT3/607S0NP7XNiI47bTTmD9/PrfccktDam7atIkpU6Zw8cUXc9xxx3HppZeyc+fOhtQeaNmyZSxevLghtaZPn86HP/xh2tvbmTZtGq9//es57bTT6l537ty5rFy5kmeffZZdu3axYsUKtmzZUve6Q3n66aeZNm0a0P8C8plnnhmxXvbHPv/PjIgTGtVIpd4ZEbEhIjZGxNVVtkdEfK6yfW1EzKt136K2bYOTT4annqprmeFqP//885x33nnccMMNHHrooQ2r29raysMPP0xvby8PPPAAjz76aN3r/tvXvkZbWxvz58+vX60havPUU/z4xz9mzZo13Hnnndx0002sXLmy7nX7nnqKNWvWcPnll/PQQw/x2te+luuuu65+dQfUful3e/fu3SxfvpwLLrigIXV3PPYYd9xxB7/4xS/YunUrO3fu5Otf/3rda89eupSPLl3KwoULOeOMMzjmmGMYN25c3es2w3GkpOFe6q2MiE9GRJ3/y0JEtAI3AYuAOcDiiJgzaNgioLPyuAz4wn7sW84nPwk/+hF84hN1KzFc7T3XXMN5553He9/7Xs4999yG1R0450mTJnHKKafw3e9+t+51f3z99SxfvpyOjg4uvPBC7r33Xt73vvfVr+6A2nziEy+/7dLW1sa73/1uHnjggbrXnbFsGTNmzHj5DPD8889nzZo19as7oPZLP+c777yTefPmMXXq1IbU/d6VVzJr1iymTJnC+PHjOffcc7nvvvsaUvsDGzeyZs0aVq5cyeTJk+ns7GxI3WrHkalTp7Jt2zYAtm3bRltbW8NqH4jov74yxMaIU4EvA78G3peZPyta/ZW1/gr4b5l5emX5YwCZee2AMf8b+H5m/lNleQNwCtAx3L7VdHd3537dxXDQQTDojpPNwFnAlJNPrv15XoW7fvhDXvPiiwAksASYDHyqpYXTTzyxIXUBtgPjgYNaWnjHW9/K2rVraW9v5w1veENd677kDy0tHH/00WzZsoWjjz66aM2hau8EXgQmDJjzzJkzi18/qzbnE4HPR/DBk05i8+bN7N27lyOPPLJo3aFqA/wtsPaoo3jjG99YvGa1uj8BLgF+FMG7TjyRDRs2cMghhzB9+vS6134GaAOeiOCYiRM57rjjGD9+fN3rwh8fR37+858zfvx42tvbefLJJ9mzZ0+Rn/tQtc+O4NFB60855RQ+85nPVL2rMSJWZ+YfbdjnGUpm3gMcDTwE9ETEf97fCeyH6cDANy17K+tqGVPLvgBExGUR0RMRPdu3b9+/Djdtgve8h92VU+G/Bea3tLA+glWrVr38iqIeFi9YwN1tbbzQ0sKPga8By1tbmTpxIj09PTz77LN1rwuwOYLuceNomziRNWvWcNhhhxUPk2p1X2hp4e62NhY34HrN4NpPRnDsgDlPnjy5LjdjVJvzew87jJMOPpienh6ef/552tvbi9cdqva/TpnCt8eN4/DDD69LzWp1j2lpYd7BB9P+mtfQ09NDZr58LaHetd8NzGxtpfugg+js7KxLmFSrW+040t7ezo4dO3jggQfYsWNHsZ/74NqLW1v5q4kT2dDayowZM/jyl7/M7bffzowZM1i1ahVnnXUWp59+eu0Fqt36Ve0BnAf0Ac8Dzw181Pocwzz/BcCXBiy/H/hfg8Z8B/ibAcv3APNr2bfa41XdNrx0aWZLS+bEif1fL798/5/j1Rqp2mOt7kjWds7OeRTU5kBuG46IbuC/A08AVwIfHPQooRc4YsDyDGBrjWNq2beMp5+GpUvh/vv7vzbygtpI1R5rdUeytnN2zqO49nDXUMYB1wAfBT4PXJ2Zdfmz1Uqtx4FTgV8CDwLvycx1A8acBVwBnAkcD3wuMxfUsm81+30NRZI05DWU4e7eepD+a7+Lsv96St1kZl9EXAHcBbQCt2bmuohYWtl+M7CC/jDZCOwCLt7XvvXsV5L0SsOdoXwNuCIzf9u4lhrHMxRJ2n+v6gwlM99fv5YkSX9K/OgVSVIRBookqQgDRZJUhIEiSSrCQJEkFWGgSJKKMFAkSUUYKJKkIgwUSVIRBookqQgDRZJUhIEiSSrCQJEkFWGgSJKKMFAkSUUYKJKkIgwUSVIRBookqQgDRZJUhIEiSSrCQJEkFWGgSJKKMFAkSUUYKJKkIgwUSVIRBookqQgDRZJUhIEiSSrCQJEkFdEUgRIRkyPi7oh4ovL1sCHGnRERGyJiY0RcPWD9/4iIxyJibUTcHhGTGta8JAlokkABrgbuycxO4J7K8itERCtwE7AImAMsjog5lc13A3Mz8y3A48DHGtK1JOllzRIo5wC3Vb6/DXhXlTELgI2ZuSkzdwPLKvuRmf8nM/sq4+4HZtS3XUnSYM0SKFMzcxtA5WtblTHTgS0Dlnsr6wa7BLizeIeSpH0a16hCEfE94I1VNn281qeosi4H1fg40Af84z76uAy4DKC9vb3G0pKk4TQsUDLzHUNti4inI2JaZm6LiGnAM1WG9QJHDFieAWwd8BxLgLOBUzMzGUJm3gLcAtDd3T3kOEnS/mmWt7yWA0sq3y8B7qgy5kGgMyJmRcQE4MLKfkTEGcBHgXdm5q4G9CtJGqRZAuU6YGFEPAEsrCwTEW+KiBUAlYvuVwB3AeuBb2Tmusr+fw8cAtwdEQ9HxM2NnoAkjXUNe8trXzLzWeDUKuu3AmcOWF4BrKgy7s/r2qAkaVjNcoYiSRrlDBRJUhEGiiSpCANFklSEgSJJKsJAkSQVYaBIkoowUCRJRRgokqQiDBRJUhEGiiSpCANFklSEgSJJKsJAkSQVYaBIkoowUCRJRRgokqQiDBRJUhEGiiSpCANFklSEgSJJKsJAkSQVYaBIkoowUCRJRRgokqQiDBRJUhEGiiSpCANFklSEgSJJKsJAkSQVYaBIkopoikCJiMkRcXdEPFH5etgQ486IiA0RsTEirq6y/cMRkRFxeP27liQN1BSBAlwN3JOZncA9leVXiIhW4CZgETAHWBwRcwZsPwJYCDzZkI4lSa/QLIFyDnBb5fvbgHdVGbMA2JiZmzJzN7Csst9LrgeuArKOfUqShtAsgTI1M7cBVL62VRkzHdgyYLm3so6IeCfwy8z86XCFIuKyiOiJiJ7t27cfeOeSJADGNapQRHwPeGOVTR+v9SmqrMuIOLjyHKfV8iSZeQtwC0B3d7dnM5JUSMMCJTPfMdS2iHg6IqZl5raImAY8U2VYL3DEgOUZwFbgSGAW8NOIeGn9mohYkJlPFZuAJGmfmuUtr+XAksr3S4A7qox5EOiMiFkRMQG4EFiemY9kZltmdmRmB/3BM88wkaTGapZAuQ5YGBFP0H+n1nUAEfGmiFgBkJl9wBXAXcB64BuZuW6E+pUkDdKwt7z2JTOfBU6tsn4rcOaA5RXAimGeq6N0f5Kk4TXLGYokaZQzUCRJRRgokqQiDBRJUhEGiiSpCANFklSEgSJJKsJAkSQVYaBIkoowUCRJRRgokqQiDBRJUhEGiiSpCANFklSEgSJJKsJAkSQVYaBIkoowUCRJRRgokqQiDBRJUhEGiiSpCANFklSEgSJJKsJAkSQVEZk50j2MmIjYDvz7q9z9cOBXBdsZDZzz2OCcx4YDmfPMzJwyeOWYDpQDERE9mdk90n00knMeG5zz2FCPOfuWlySpCANFklSEgfLq3TLSDYwA5zw2OOexoficvYYiSSrCMxRJUhEGiiSpCANlGBFxRkRsiIiNEXF1le0REZ+rbF8bEfNGos+SapjzeytzXRsR90XEMSPRZ0nDzXnAuL+MiL0RcX4j+yutlvlGxCkR8XBErIuIHzS6x9Jq+L1+fUT8a0T8tDLni0eiz5Ii4taIeCYiHh1ie9njV2b6GOIBtAI/B/4MmAD8FJgzaMyZwJ1AACcAPxnpvhsw57cCh1W+XzQW5jxg3L3ACuD8ke67zj/jScDPgPbKcttI992AOf8d8KnK91OAXwMTRrr3A5z3ScA84NEhthc9fnmGsm8LgI2ZuSkzdwPLgHMGjTkH+Gr2ux+YFBHTGt1oQcPOOTPvy8wdlcX7gRkN7rG0Wn7OAB8E/gV4ppHN1UEt830P8K3MfBIgM8fCnBM4JCICeB39gdLX2DbLysyV9M9jKEWPXwbKvk0HtgxY7q2s298xo8n+zucD9L/CGc2GnXNETAfeDdzcwL7qpZaf8V8Ah0XE9yNidURc1LDu6qOWOf89MBvYCjwCXJmZLzamvRFT9Pg17oDb+dMWVdYNvs+6ljGjSc3ziYi30R8of1PXjuqvljnfAHw0M/f2v4Ad1WqZ7zhgPnAqcBCwKiLuz8zH691cndQy59OBh4G3A0cCd0fEDzPzuTr3NpKKHr8MlH3rBY4YsDyD/lcv+ztmNKlpPhHxFuBLwKLMfLZBvdVLLXPuBpZVwuRw4MyI6MvMbzekw7Jq/b3+VWbuBHZGxErgGGC0Bkotc74YuC77Ly5sjIhfAG8GHmhMiyOi6PHLt7z27UGgMyJmRcQE4EJg+aAxy4GLKndLnAD8NjO3NbrRgoadc0S0A98C3j+KX7EONOycM3NWZnZkZgfwz8B/HKVhArX9Xt8BnBgR4yLiYOB4YH2D+yypljk/Sf8ZGRExFTgK2NTQLhuv6PHLM5R9yMy+iLgCuIv+u0Ruzcx1EbG0sv1m+u/4ORPYCOyi/1XOqFXjnP8r8Abg85VX7H05ij+ptcY5/8moZb6ZuT4ivgusBV4EvpSZVW89HQ1q/Bl/EviHiHiE/reCPpqZo/oj7SPin4BTgMMjohe4BhgP9Tl++dErkqQifMtLklSEgSJJKsJAkSQVYaBIkoowUCRJRRgokqQiDBSpCURES0SsjIjBf0R6cOUj178wUr1JtTJQpCZQ+RDC/wC8PSIuGbDpU/T/AfKHR6IvaX/4h41SE6n85fangaOBP6f/L7tPycwfjWhjUg0MFKnJRMRd9H/CbwewLDOvGtmOpNoYKFKTiYhZ9P/rgj8H5mbmH0a4JakmXkORms8lwO/p/yjxPxvhXqSaeYYiNZGI+EvgPuCdwOXAVOCtmbl3RBuTauAZitQkImIi8FXgHzLzTuAy+i/Mew1Fo4JnKFKTiIjrgXcBb8nM31XWXQjcBswfzf8eicYGA0VqAhFxEnAv8I7M/P6gbd+g/1rKCZnZNwLtSTUxUCRJRXgNRZJUhIEiSSrCQJEkFWGgSJKKMFAkSUUYKJKkIgwUSVIRBookqYj/BzEBY8awzGuUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from FEToy.mesh.lagrange1dmesh import mesh1d\n",
    "from FEToy.fe.shapefun import shape1d\n",
    "from FEToy.fe.gaussrule import gausspoint1d\n",
    "\n",
    "mymesh=mesh1d(nx=10,meshtype='edge2')\n",
    "mymesh.createmesh()\n",
    "mymesh.plotmesh(withnode=True,withnodeid=True)\n",
    "\n",
    "gpoints=gausspoint1d(ngp=2)\n",
    "gpoints.creategausspoint()\n",
    "print(gpoints.gpcoords)\n",
    "\n",
    "shp=shape1d(meshtype='edge2')\n",
    "shp.update()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2ab1e26-333d-49d2-86d0-796b87049e09",
   "metadata": {},
   "source": [
    "$\\phi(x)=-\\frac{f}{2}x^{2}+ax+b$, with $a=\\frac{y_{1}-y_{2}}{x_{1}-x_{2}}+\\frac{f}{2}(x_{1}+x_{2})$ and $b=y_{1}-ax_{1}+\\frac{f}{2}x_{1}^{2}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f3a6b330-b988-47a7-960c-3592b29e55cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here are your parameters\n",
    "f=200.0\n",
    "x1=0.0;y1=0.5 # left\n",
    "x2=1.0;y2=5.0 # right\n",
    "\n",
    "def Analytical(x):\n",
    "    a=(y1-y2)/(x1-x2)+0.5*f*(x1+x2)\n",
    "    b=y1-a*x1+0.5*f*x1*x1\n",
    "    return -0.5*f*x*x+a*x+b"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad6e3848-4efc-4af8-988c-80a47832c2d0",
   "metadata": {},
   "source": [
    "The K matrix and the right hand side vector can be defined as follows:\n",
    "\n",
    "$K=\\int_{\\Omega}\\nabla N^{J}\\nabla N^{I}dV$, $F=\\int_{\\Omega}fN^{I}dV$\n",
    "\n",
    "where $I$ and $J$ are the I-th and J-th node of the current element, respectively. $N$ and $\\nabla N$ represent the shape function value and its derivatives, respectively. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "525073ab-eb32-4bf0-9d46-3fa386efa201",
   "metadata": {},
   "outputs": [],
   "source": [
    "nDofs=mymesh.nodes*1 # total degrees of freedom\n",
    "K=np.zeros((nDofs,nDofs)) # Your K matrix\n",
    "F=np.zeros(nDofs)         # Your RHS\n",
    "\n",
    "# do element loop\n",
    "for e in range(mymesh.elements):\n",
    "    elconn=mymesh.elementconn[e,:]  # local element connectivity\n",
    "    nodes=mymesh.nodecoords[elconn] # node coordinates of local element\n",
    "    # now we can do the integration\n",
    "    for gp in range(gpoints.ngp):\n",
    "        xi=gpoints.gpcoords[gp,1] # 1->xi, 2->eta\n",
    "        w =gpoints.gpcoords[gp,0] # weight\n",
    "        \n",
    "        shp_val,shp_grad,j=shp.calc(xi,nodes) # j->mapping from x->xi\n",
    "        JxW=j*w # the weight times the determinate of jacobian transformation\n",
    "        \n",
    "        # now we can calculate the F term\n",
    "        for i in range(mymesh.nodesperelement*1):\n",
    "            iInd=elconn[i]\n",
    "            # assemble to global F\n",
    "            F[iInd]+=f*shp_val[i]*JxW\n",
    "            # now we can calculate the K matrix\n",
    "            for j in range(mymesh.nodesperelement*1):\n",
    "                jInd=elconn[j]\n",
    "                K[iInd,jInd]+=shp_grad[j]*shp_grad[i]*JxW\n",
    "#\n",
    "#print(K)\n",
    "#print(F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8d7f3bf8-d57e-4ef9-a101-e315e22bad09",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we should apply our boundary conditions\n",
    "iInd=mymesh.bcnodeids['left']\n",
    "K[iInd,iInd]+=1.0e16 # here we use the penalty method, 1.0e16 is the 'large' number for penalty\n",
    "F[iInd]+=y1*1.0e16\n",
    "\n",
    "iInd=mymesh.bcnodeids['right']\n",
    "K[iInd,iInd]+=1.0e16\n",
    "F[iInd]+=y2*1.0e16\n",
    "\n",
    "#print(K)\n",
    "\n",
    "# now we solve K*u=F by u=K\\F\n",
    "phi=np.linalg.solve(K,F)\n",
    "#print(phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "51b0917f-976a-4b78-a985-fc5d71e8cb92",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f1c1db0dca0>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAoqElEQVR4nO3dd3hUVf7H8fc3M8kkpAEpEAIkQOggLXRUEBuIDcuiqIgiupbFdYvlt7oWXHV31V0rFuzYZS3YFulSZBOEhA4hJJSQBALpPef3R+IuKkhIZnLnznxfz5Mnk8lk5nMFPzmcufccMcaglFLKfgKsDqCUUqpptMCVUsqmtMCVUsqmtMCVUsqmtMCVUsqmnC35YtHR0SYxMbElX1IppWwvNTX1oDEm5qf3t2iBJyYmkpKS0pIvqZRSticiWce6X6dQlFLKprTAlVLKprTAlVLKprTAlVLKprTAlVLKprTAlVLKprTAlSo+AK9OgOJcq5ModVJa9DxwpbxR9eJHcWat5sCnD7Cu/70UV1RTXFFDRXUtdQZqjcEhQnBgACFBDsJcTqLCXESFBhEb4SImzIWIWH0Yyg9pgSu/UVxRzcZ9RWzcV8iOvGJmbxxPENUENnw/bsc8ztsxjwoTSK/K1xv9vMGBAXRu24qEqFB6tQ+nd1wEfeIiSIhqpcWuPEoLXPmsvOIK1uwqYHXGIb7LPMSu/NL/fi86zMVtsa9xU9Vr9CtaQWBdBbWOYEq6nkv52AdJj+pAcKCDABECBOoMVFTXUlFdS1FFDQWllRwsqeJAYQV7CsrILihj18FSFm/No7aufpOUqNAgBie0YWhiG0YnRdMnLkILXbmVFrjyGcYYNu4r4pstuSzamsvGfUUAhLmcDE1sw8UD4+nfMZL+8ZFEhbnqf+izRbCuCpzBOGqriIyMIjI+4WfP7RAIdTkJbZg+6RIdeswMFdW17MwrIX1fISm7D7Mu+zALN9fPrUeHuTitRzRn92nP2J4xBAc66uffP5wOl74G4e088t9F+S4tcGV723OL+XT9fj5L20/WoTJEYHDnNvzhnJ6MSYqmb4cInI7jvF9fmgdDpkPydEh5FUqa90ZmcKCDfvGR9IuP5IphnQHILapg+fZ8lu84yOKtecxft4+QQAdn9IrlrroX6Zi9Bln2GEx6olmvrfyPtOSemMnJyUYXs1LuUFJZw2cb9vPu2mw27C0kQGB0UjSTTonjzN7t/jfC9jLVtXV8t6uA4W/3JtBU/fwBThf8Ka/lgymvJiKpxpjkn96vI3BlK7sPlvLqykw+TN1LaVUtPdqFcd+kPpw/oAMx4d5Z2kcLdAQwpns03JFO3df/h9myAEdtBeUmiK9qh7Iw7lYu3HSAM3u3wxGg8+Xql2mBK1tYl32Y55dm8M2WXJwBwvmndGDqiAQGd25tzzcGw9sT4IqAuvr59+CaKpI6xfG3AhdfvJlKYlQrZp7WjUuGxONyOqxOq7yUFrjyamszC3h68Q5W7DhI61aB3DI2iWtGJhAbEWx1tOY7av5dUl6lf0kuy2eM49+bc5mzLIN7/pXOP77Zzo2nd2Pq8M71b3oqdRSdA1fe46gzMjYWBfPYV1tZseMg0WFBzDytK1OHJxDq8o8xhzGGVRmHeHbJTlZlHKJ9RDC3npHE5cmdCHLqBdT+5nhz4FrgynssuAOT8grfRp7P1blTaNMqkFvGJTF1eAIhQf47+lyVcZDH/72d1KzDdIkO5e4JvTirTzt7Th2pJtECV95rdizUVP7sbuNwIffqGRlQPyJfvDWPv3yxhYz8UkZ2jeK+8/vQOy7C6miqBRyvwPXfYspSxhgWn72QrwNOpdwE1d/nDIH+lyG3p1ucznuICON7t+Or20/jgQv6suVAEZOe/paHFmympLLG6njKIlrgyjI5heXc8EYq183fS7UzjGCpBmcwUlsJrgi9MvEYAh0BTBuVyNLfj+Xy5E7M/TaTMx9fxlcbD1gdTVlAC1y1OGMMb67J4qwnlvPtznzumdiLiV0cSPJ1MOOb+jMzmnlFpK9r3SqIRyb3Z/7No2gTGsRNb6Vy69vrOFTy86ko5bt0Dly1qJzCcv74YRordhzk1O7RPHxRfzpHtbI6lq1V19YxZ2kGTy3eQURwIA9d1I+J/eOsjqXcSOfAleU+/n4fZz+5nNSswzx8cT/euG6YlrcbBDoCuG18dxbcdirxbUK4ed46fv/BBp0b9wP+cVKtslRJZQ33fbyR+d/vIzmhDY9fPoCEqGOv5qearmf7cD769SieXrSDZ5bsZG1mAf+YMpDBndtYHU15yAlH4CLSSUSWiMgWEdkkIrMa7r9fRPaJyPqGj4mej6vsZuO+Qs5/+ls+Xr+P28/szns3jtTy9qBARwB3nN2T924cSZ0xXD5nNS+v2EVLTpWqlnPCOXARiQPijDHrRCQcSAUuAi4HSowxf2/si+kcuH95d202932yibahQfxzykCGd42yOpJfKaqo5o8fpPHVpgOc1acdf790AJGtAk/8g8rrNHkO3BiTY4xZ13C7GNgCxLs/ovIVFdW1/PHDDdw1P53hXdvyxaxTtbwtEBEcyPNXDea+SX1YsjWPSc+sYEtOkdWxlBud1JuYIpIIDAK+a7jrVhFJE5FXREQn2hT7j5Rz6ZxVvJ+yl9vOSOK16cNoGxpkdSy/JSJcN6YL7980kuoaw+TnVvF5Wk79ujOvToBiPV3Tzhpd4CISBnwE3G6MKQKeB7oBA4Ec4PHj/NxMEUkRkZT8/PzmJ1Zea132YS54ZiW7D5bx0jXJ/O7snrqmtZcY3LkNn942mj4dIrjl7XV8/+bdmOw1sOwxq6OpZmjUeeAiEggsAL42xvxs36eGkfkCY0y/X3oenQP3XfPX7eWuj9JpHxnMy9OS6dEu3OpI6hjM7FjkGOvO6E5A3q3Jc+BSv+TZXGDL0eXd8ObmDy4GNrojqLIXYwxPLtzOHe9vYHBCaz65ZbSWtxeTWWmYfpdSE1C/nnoFLsp7XQKzdN0ZO2rMeeCjgauBdBFZ33DfPcAVIjIQMMBu4EYP5FNerLq2jrvnp/Nh6l4uGdyRRyb317WqvV14e8QVgdNUURvgIqi2io+3l9B/XCu66+9d2zlhgRtjvgWONZH5hfvjKLsoqazh12+lsmLHQWaN787tZ3bX9antomEnIEfydA4te4HWW7dz6ZzVvHJtMkMS2lqdTp0EXQtFnbSC0iqufXUtm/YX8cjF/bl8aCerI6lm2FNQxjWvrGX/kXKevXIwZ/bRVSC9ja6Fotxi/5FyLpuzim0HinnhqiFa3j6gU9tWfHjTSHq1D+fGt1KZv26v1ZFUI2mBq0bLPFjKZXNWk1dUyRvXDdORmg+JCnPx9g0jGNG1LXe8v4E312RZHUk1gha4apSdecVc/sJqyqtreWfmCL2y0geFupzMnTaUM3vHcu/HG3lhWYbVkdQJaIGrE9p6oIhfvbAGY+DdmSPoFx9pdSTlIcGBDp6/agiTTonjkS+38tSiHVZHUr9Al5NVv2jz/iKmvryGIGcAb98wgm4xYVZHUh4W6Ajgn1MGEeQM4ImF2wH4zfjuFqdSx6IFro5r24Fipr68huBAB+/OHKHLwPoRR4Dwt0sHAGiJezEtcPVjxQfgw+lkjnuGqW/tIsgZwDs3aHn7o5+WuCNAuGVcksWp1NG0wNWPLfsrJms1Ka/fjThn8vYNI0iM1vL2Vz+UeF2d4W9fb6NVkIPpo7tYHUs10AJX9WbHQsMiRwJcZr7msuqv4QVd5MjfOQKEv182gPLqWh74bDOtghz8amhnq2Mp9CwU9YNZaVT0mkwFLgDqHMHQ/zJd5EgB4HQE8NQVgzi9Rwx3zU/n0w37rY6k0AJXDYoCo1i0q5wgU0Wdw0VAXRW4IiBcL9ZR9VxOBy9cPYShiW353fvrWb5d1/e3mha4oqK6lhteT8FZfpCc7lcScMMiGDIdSnS3FvVjwYEOXp6WTFJsODe9lcr6PUesjuTXtMD9XF2d4Xfvb+C7zAIqLnmd+Kueg/b9YdITMGWe1fGUF4oIDuT164YSHeZi+qtr2ZlXYnUkv6UF7seMMTy4YDOfp+fwfxN7c+FA3ataNU5seDBvXj8MR4Bw7atryS8+xi4/yuO0wP3YSyt28dqq3Vw/pgs3nNbV6jjKZhKiQnnl2qEcKqni+tf/Q1lVjdWR/I4WuJ/6Mj2Hv3yxlfP6x/F/E3tbHUfZ1CkdW/P0FYPYuK+Q297+npraOqsj+RUtcD/0ffZhbn9vPYM7t+bxywcQoDvHq2Y4s087HrigL4u25vHggs1Wx/EreiGPn9lTUMYNb6TQLiKYl65JJjjQYXUk5QOuHplIdkEZL63IJCk2jGtGJlodyS9ogfuRksoaZryeQlVNHe/OHEpUmMvqSMqH3DWhN7vyS3ngs80kRoVyWo8YqyP5PJ1C8RN1dYbfvreenfklPDd1CEmxuiysci9HgPDPKwbRPTaMW+atY2desdWRfJ4WuJ94fOE2Fm7O5d7zejOme7TVcZSPCnM5eXlaMq5ABzNeT6GwrNrqSD5NC9wPfLJ+H88uyeCKYZ2ZNirR6jjKx3Vs04o5Vw1m35Fybnv3e2rrjNWRfJYWuI/btL+QOz9KY1hiWx64oC8iesaJ8rzkxLY8cEE/lm/P569fb7U6js/SNzF92JGyKm56K5XWIUE8O3UwQU79fa1azpXDO7NpfyEvLNtFn7gIvdLXA/T/aB9VW2e47Z3vyS2sZM7VQ4gJ1zNOVMv78/l9GZrYhrs+SmfbAX1T0920wH3UEwu3sWLHQR68sC8DO7W2Oo7yU0HOAJ69cjBhwU5ueiuVogp9U9OdtMB90KItuTy7JIMpQzsxZZjunKKsFRsRzLNXDia7oIw/fLABU5QDr06AYl2uuLlOWOAi0klElojIFhHZJCKzGu5vKyILRWRHw+c2no+rTmRPQRm/fW89fTtEcP8Ffa2OoxQAw7q05e4Jvfh6Uy6b3r0XstfAssesjmV7Yswvn+IjInFAnDFmnYiEA6nARcC1QIEx5lERuQtoY4y585eeKzk52aSkpLgluPq5yppaLn1+NbsPlbLgtjG6k7zyKmZ2LFJzjGVnnbrv6omISKoxJvmn959wBG6MyTHGrGu4XQxsAeKBC4HXGx72OvWlriw0e8EW0vcV8vfLBmh5K68js9Ko7nMJFQQBYJwhuu9qM53UHLiIJAKDgO+AdsaYHKgveSD2OD8zU0RSRCQlP1/30POUL9JzeHNNFjPGdOGcvu2tjqPUz4W3JzAkEhfVVJhATE0FJihc911thkYXuIiEAR8Btxtjihr7c8aYF40xycaY5JgYXdzGE/YUlHHnR2kM6NSaP57by+o4Sh1faR6SfB3fjH6bt2rGk5mVaXUiW2vUhTwiEkh9ec8zxsxvuDtXROKMMTkN8+Q6iWWB6to6bnvnezDw9JRBerGO8m4N+6yeZwy35kfzwKYDfJh9mEGd9RyIpmjMWSgCzAW2GGOeOOpbnwLTGm5PAz5xfzx1Io//ezvr9xzh0UtOoXNUK6vjKNUoIsJfJvenfUQwv3n3e4r1/PAmacxwbTRwNXCGiKxv+JgIPAqcJSI7gLMavlYtaOXOg8xZVr9I1XmnxFkdR6mTEhkSyFNXDGT/kQru/Xij1XFs6YRTKMaYb4HjrYA03r1xVGMdLq3id+9voFtMKPdN6mN1HKWaZEhCW2aN784TC7dzWo8YJg/uaHUkW9EJUxsyxnD3/HQOlVbyzymDCAnSbdGUfd0yLolhXdpy78cbyT5UZnUcW9ECt6H3U/bw1aYD/OGcnvSLj7Q6jlLN4ggQnvzVQAJEuOP99bp++EnQAreZrEP1ew6O6hbFjDFdrY6jlFvEtw7hwYv6kpJ1mDnLMqyOYxta4DZSW2e44/0NOAKExy8fQECAbs6gfMdFA+M575Q4nly4nfS9hVbHsQUtcBt5cfkuUrMO89CF/YiLDLE6jlJuJSI8fFE/osNc3P7e91RU11odyetpgdvE5v1FPLFwGxP7t+fCgR2sjqOUR7RuFcTfLjuFjPxS/v71NqvjeD0tcBuorKnljvfXExkSxOyL+uu+lsqnndo9hqnDOzN3ZSZrMwusjuPVtMBt4JnFO9l6oJjHLulP29Agq+Mo5XH3TOxNxzYh/OHDDZRV1Vgdx2tpgXu59L2FPLc0g8mD4xnfW1dtU/4h1OXk75cOILugjEe/1F3tj0cL3ItV1tTy+w82EB0WxJ8n6e46yr8M7xrF9FFdeGN1FqszDlkdxytpgXur4gPkPTWegtw9PDK5P5GtAq1OpFSL+8M5PUmMasWdH6XpVMoxaIF7qUNfPER84XqejPuaM3rp1InyTyFBDh695BSyC8p4/N/brY7jdRq1HrhqQbNjoaaSKACBMYc/gfsjdd9A5bdGdI3iqhGdeWVlJuedEsdgXTv8v3QE7m1mpbE99lzKTcPZJrpvoFLcNaE3HSJD+MMHG/QCn6NogXuZzMpwUg7U4JJqcAZDbSW4InTfQOXXwlxO/jK5Pxn5pTy3VNdK+YEWuBepqzPc9VEasQFFVAyYBjO+gSHToSTX6mhKWe70HjFcNLADzy/dyfbcYqvjeAUtcC/yfsoevsss4ODEubS6+J/Qvj9MeuK/+wgq5e/undSHMJeTu+enU6fLzmqBe4v84kr+8sUWhnVpy6+GdrI6jlJeKSrMxZ/O60Nq1mHeXpttdRzLaYF7iYc/30x5dS1/uVjXOlHql0weHM+YpGge+3IruUUVVsexlBa4F1ixI5+P1+/n12OTSIoNszqOUl5NRHj44n5U1dbx4ILNVsexlBa4xSqqa/m/f22ka3QoN4/tZnUcpWwhISqUW8cl8XlaDku3+e/1EVrgFnt2yU6yC8qYfVE/ggN1c2KlGmvm6V3pGhPKfZ9s8ttzw7XALbQrv4QXlu3iooEdGJUUbXUcpWzF5XQw+8J+ZBeU8eySnVbHsYQWuEWMMfz50024nAHcc15vq+MoZUujkqK5eFA8c5ZlsDOvxOo4LU4L3CJfpB9gxY6D/P6cnsSGB1sdRynbumdib4IDHdz/6SaM8a9zw7XALVBSWcODCzbRt0MEV41IsDqOUrYWE+7id2f14NudB/ly4wGr47QoLXALPL1oB7lFlcy+qB+OAD3nW6nmumpEAr3ahzN7wWa/Wjf8hAUuIq+ISJ6IbDzqvvtFZJ+IrG/4mOjZmL4jI7+EV1ZmcnlyRwbpsphKuYXTEcBDF/Vjf2EFzyz2nzc0GzMCfw049xj3P2mMGdjw8YV7Y/kmYwwPfLaZ4EAHfzy3l9VxlPIpQxPbMnlQPC+t2MWufP94Q/OEBW6MWQ4UtEAWn7dwcy7Lt+fz2zN7EB3msjqOUj7nrom9CHY6eMhPrtBszhz4rSKS1jDFcty5ABGZKSIpIpKSn5/fjJezt4rqWh5csJke7cK4eqS+camUJ8SGB/Ob8d1Zsi2fJVt9/wrNphb480A3YCCQAzx+vAcaY140xiQbY5JjYmKa+HL299LyXew9XM795/cl0KHvHSvlKdNGJdI1JpSHFmymqqbO6jge1aQmMcbkGmNqjTF1wEvAMPfG8i0HCit4bmkGE/q11ysulfKwIGcA907qw66Dpby2KtPqOB7VpAIXkbijvrwY2Hi8xyr461dbqTWGeybqFZdKtYRxPWM5o1csTy3aSV6x7y4525jTCN8BVgM9RWSviFwP/FVE0kUkDRgH/NbDOW3r++zDzP9+HzPGdKFT21ZWx1HKb/zpvN5U1tTy+NfbrY7iMc4TPcAYc8Ux7p7rgSw+p66u/rTBmHAXN49LsjqOUn6la0wY00YmMndlJteMSqBvh0irI7mdvpvmQZ9u2M/6PUf44zk9CXOd8HelUsrNbhvfndYhgTz8+RafXCdFC9xDKqpr+etXW+kXH8ElgztaHUcpvxQZEsjtZ/ZgVcYhlq/bCK9OgOJcq2O5jRa4h8z9NpP9hRX86bw+BOh6J0pZ5srhnekWE8rhL2djstfAssesjuQ2+u96D8gvruS5JTs5q087RnSNsjqOUn4t8JH2LKqp/N8dKXPrP5wu+JO9L/bREbgHPPnNdipr6rh7gq53opTlZqVh+l1KpdQvX2GcIdD/MpiVbnGw5tMCd7PtucW8uzabq0Yk0DVGd5hXynLh7RFXBEGmigoTiKmpAFcEhLezOlmzaYG72aNfbiXU5eQ347tbHUUp9YPSPCT5Op7p9gJv151J2eH9VidyC50Dd6PVGYdYvDWPuyb0om1okNVxlFI/mDIPgKmF5Yz9WytSAtvzD2sTuYWOwN2krs7wyJdbiIsM5tpRiVbHUUodQ1xkCNeP6cLH6/eTvrfQ6jjNpgXuJp+n55C2t5A7zupBcKDD6jhKqeP49dhutA0N4pEv7X9xjxa4G1TV1PG3r7fRq304k/WiHaW8WnhwILedkcSqjEOs2HHQ6jjNogXuBvO+yyK7oIw7J/TSTYqVsoErh3emY5sQHvtqK3V19h2Fa4E3U0llDc8s3smIrm0Z28N/N6xQyk5cTge/O7sHm/YX8Vmafc9I0QJvprkrMjlUWsWd5/ZCREffStnFhQPi6R0XweP/3m7bnXu0wJuhoLSKl1bs4py+7RjU+bjbgiqlvFBAgPDHc3uSXVDGO2uzrY7TJFrgzfDckp2UVdXw+7N7Wh1FKdUEY3vEMLxLW55eXP//st1ogTfRviPlvLEmi8mDO9K9XbjVcZRSTSBSPwo/WFLJa6t2Wx3npGmBN9E/v9kOBm4/Uy+ZV8rOhiS0ZXyvWOYszaCwvNrqOCdFC7wJduWX8NG6fQ2nIuk+l0rZ3e/O7klRRQ0vLd9ldZSTogXeBP/4ZgdBjgBu0X0ulfIJfTpEcP6ADryyMpP84soT/4CX0AI/SVsP1J83eu3oRGLCXVbHUUq5yW/P7E5lTR3PLd1pdZRG0wI/SU8u3E5YkJMbT+tqdRSllBt1jQnj0sEdmfddNjmF5VbHaRQt8JOQvreQrzflMuPUrrRupcvFKuVrbj0jCWMMzy6xxyhcC/wkPL5wG21aBXLdmESroyilPKBT21ZcntyJ9/6zh72Hy6yOc0Ja4I2UmnWYpdvymXlaN8KDA62Oo5TykFvPSEIQnlns/aNwLfBG+sc324kKDWLaqASroyilPCguMoQrh3fmg9S9ZB0qtTrOL9ICb4TUrAJW7DjIjad3pVWQ7kKnlK+7eWw3nAHCPxftsDrKLzphgYvIKyKSJyIbj7qvrYgsFJEdDZ99dyWn4gO0mncBPUNLuWqEjr6V8gexEcFcPSKBj7/fR+ZB7x2FN2YE/hpw7k/uuwtYZIzpDixq+Non5S54iJ6VG3my/b919K2UH5l5eleCnAFePRd+wgI3xiwHCn5y94XA6w23Xwcucm8sLzA7Fu6PpN22twgQQ599H8D9kfX3K6V8Xmx4MFcNT+Dj9fvY7aWj8KbOgbczxuQANHw+bquJyEwRSRGRlPz8/Ca+nAVmpXGwywWUm4bzvZ0h0P8ymJVubS6lVIuZeXpXnAHCM156XrjH38Q0xrxojEk2xiTHxNhoy7Hw9mzIq8Ul1RinC2orwRUB4e2sTqaUaiGx4cFcNSKBf32/zyvPSGlqgeeKSBxAw+c890XyDqlZh6kuymNL/GXIjEUwZDqU5FodSynVwm5sGIU/7YVz4U0t8E+BaQ23pwGfuCeO93h68Q7uCbqTLtOeh/b9YdITMGWe1bGUUi0sNjyYqcPrR+F7Crzr6szGnEb4DrAa6Ckie0XkeuBR4CwR2QGc1fC1z9iw5whLt+Vzw6l63rdSqn4U7hDhuaUZVkf5kRO2kzHmiuN8a7ybs3iNpxbtoHWrQK4eqed9K6WgXUQwlw/tyHv/2cNtZyTRoXWI1ZEAvRLzZzbtL2TR1jyuH92FMJeOvpVS9W46vRvGwItetGuPFvhPPLckg3CXk2tGJVodRSnlRTq2acXkwfG8szabvOIKq+MAWuA/sjOvhC825nDNqAQiQ3TFQaXUj908Nonq2jpeXpFpdRRAC/xHnlu6k2Cng+tGd7E6ilLKCyVGh3LBgA68uTqLgtIqq+Nogf9gT0EZn6zfzxXDOhMVpntdKqWO7eZxSZRX1/Laqt1WR9EC/8GcZRk4RJipe10qpX5Bj3bhnN2nHa+tzKSkssbSLFrgQF5RBR+k7OWSIR1pHxlsdRyllJe7eVwSRRU1vP1dlqU5tMCBud9mUlNXx02n6+hbKXViAzu1ZnRSFC+tyKSiutayHH5f4IVl1by1JovzTulAQlSo1XGUUjZxy9gk8osr+TB1r2UZ/L7A31yzm9KqWn59ejeroyilbGRktygGdGrNC8szqKmtsySDXxd4eVUtr67czdieMfTpEGF1HKWUjYgIN4/txp6Ccj5Pz7Ekg18X+PspezhUWsXNY5OsjqKUsqGzerejW0woc5btwhjT4q/vtwVeXVvHi8t3kZzQhmFd2lodRyllQwEBwo2ndWNLThHLtrf8jmN+W+Cfp+Ww70g5N+nct1KqGS4c1IH2EcHMWdbyS836ZYEbY3hh+S66x4ZxRi/dpFgp1XQup4Prx3Rhza4C1u850qKv7ZcFvmLHQbbkFHHDaV0JCBCr4yilbO6K4Z2JCHYyp4U3fPDLAn9heQbtIlxcOLCD1VGUUj4gzOXk6pEJfL35ALvyS1rsdf2uwDfuK2TlzkNcN7oLLqfD6jhKKR9x7aguBAYEMPfblltq1u8K/IXluwh3OblieGeroyilfEhMuIvJg+P5MHUvh0oqW+Q1/arA9xSU8Xnafq4c3pmIYN2wQSnlXjNO7UJlTR1vrG6ZRa78qsBfWZlJgAjTdcMGpZQHJMWGc2bvWN5ck0V5lecXufKbAi8sr+b9/+zhggEddMlYpZTH3HBqVwpKq/honecXufKbAn9nbTalVbXMOFWXjFVKec6wLm0Z0DGSud9mUlvn2cvr/aLAq2rqeHVlJqOTonTRKqWUR4kIN5zWlcyDpXyzJdejr+UXBb4gbT+5RZU6+lZKtYhz+7YnvnUIcz28e73PF7gxhpdWZNI9NoyxPWKsjqOU8gNORwDTRyeydncBaXuPeOx1fL7AV2UcYktOETNO7YKIXjavlGoZvxraiTCX06MX9jSrwEVkt4iki8h6EUlxVyh3mvttJtFhQVw4MN7qKEopPxIeHMiUoZ34PC2HA/t2w6sToNi9c+LuGIGPM8YMNMYku+G53Cojv4TFW/OYOjyB4EC9bF4p1bKuHZ1InTHs++QByF4Dyx5z6/M73fpsXua1lbsJcgRw1YgEq6MopfxQx2cT2eWqhLyGO1Lm1n84XfCnvF/82cZo7gjcAP8WkVQRmXmsB4jITBFJEZGU/PyW27HiSFkVH6bu5cKBHYgJd7XY6yql1H/NSqOg6wWUm6D6r50h0P8ymJXulqdvboGPNsYMBiYAt4jIaT99gDHmRWNMsjEmOSam5c4CeWftHsqra7n+VL1sXillkfD2tG0TTbBUgzMYaivBFQHh7dzy9M0qcGPM/obPecC/gGHuCNVc1bV1vL5qN6OToujVXi/cUUpZqDQPSb4OZnwDQ6ZDifveyGzyHLiIhAIBxpjihttnAw+6LVkzfLnxAAeKKnj44n5WR1FK+bsp8/53e9ITbn3q5ryJ2Q74V8O51U7gbWPMV25J1UyvrcwkMaoV43rqfpdKKd/V5AI3xuwCBrgxi1ts2HOEddlH+PP5fXS/S6WUT/O5KzFfW7WbMJeTS4d0tDqKUkp5lE8VeF5xBQvS9nPpkI6E6447Sikf51MFPm9NNjV1hmmjEq2OopRSHuczBV5ZU8u877IZ1zOWLtGhVsdRSimP85kC/yI9h4MllVyro2+llJ/wmQJ/fVUWXWNCGZMUbXUUpZRqET5R4Bv2HGH9niNMG5mopw4qpfyGTxT466t3ExrkYPJgXfNbKeU/bF/gh0oqWbAhh0v01EGllJ+xfYG/+589VNXWcc1IXfNbKeVfbF3gNbV1zFuTxZikaJJiw62Oo5RSLcrWBf7Nllz2F1bo6Fsp5ZdsXeBvrskivnUI43u7Z3F0pZSyE9sW+M68ElbuPMSVwzvj0FMHlVJ+yLYFPu+7LAIdwq+GdrI6ilJKWcKWBV5WVcOHqXuZ2D+O6DDdsFgp5Z9sWeCfrN9PcUUNV4/QNy+VUv7LdgVujOHN1Vn0ah/OkIQ2VsdRSinL2K7A12UfYXNOEVeNSKBhP06llPJLtivweWuyCHM5uWiQrnuilPJvtirww6VVLEjP4eJB8YS5mrwfs1JK+QRbFfhH6/ZSVVPHlcM7Wx1FKaUsZ5sCN8Yw77tshiS0oXdchNVxlFLKcrYp8NUZh8g8WMpUHX0rpRRglwIvPkDMhxfTLaSEif3jrE6jlFJewRYFXvbNI3QrT+ex6C8JDnRYHUcppbyCd5/KMTsWaippBSCQnP8vuP9f4HTBn/KsTqeUUpZq1ghcRM4VkW0islNE7nJXqP+alQb9LqMmILj+a2cI9L8MZqW7/aWUUspumlzgIuIAngUmAH2AK0Skj7uCARDeHlzhOE0VOIOhthJcERCu638rpVRzRuDDgJ3GmF3GmCrgXeBC98Q6SmkeDJkOM76p/1yS6/aXUEopO2rOHHg8sOeor/cCw5sX5ximzPvf7UlPuP3plVLKrpozAj/WSlLmZw8SmSkiKSKSkp+f34yXU0opdbTmFPhe4OjtcDoC+3/6IGPMi8aYZGNMckxMTDNeTiml1NGaU+D/AbqLSBcRCQKmAJ+6J5ZSSqkTafIcuDGmRkRuBb4GHMArxphNbkumlFLqFzXrQh5jzBfAF27KopRS6iTY4lJ6pZRSPyfG/OzEEc+9mEg+kNXEH48GDroxjh3oMfsHPWb/0JxjTjDG/OwskBYt8OYQkRRjTLLVOVqSHrN/0GP2D544Zp1CUUopm9ICV0opm7JTgb9odQAL6DH7Bz1m/+D2Y7bNHLhSSqkfs9MIXCml1FG0wJVSyqa8rsBPtMuP1Huq4ftpIjLYipzu1IhjntpwrGkiskpEBliR050au5uTiAwVkVoRubQl87lbY45XRMaKyHoR2SQiy1o6o7s14u91pIh8JiIbGo55uhU53UlEXhGRPBHZeJzvu7e/jDFe80H9mioZQFcgCNgA9PnJYyYCX1K/nO0I4Durc7fAMY8C2jTcnuAPx3zU4xZTv1zDpVbn9vCfcWtgM9C54etYq3O3wDHfAzzWcDsGKACCrM7ezOM+DRgMbDzO993aX942Am/MLj8XAm+YemuA1iIS19JB3eiEx2yMWWWMOdzw5Rrql+61s8bu5nQb8BFg9x2sG3O8VwLzjTHZAMYYfzhmA4SLiABh1Bd4TcvGdC9jzHLqj+N43Npf3lbgx9rlJ74Jj7GTkz2e66n/DW5nJzxmEYkHLgbmtGAuT2nMn3EPoI2ILBWRVBG5psXSeUZjjvkZoDf1+wikA7OMMXUtE88ybu2vZq1G6AGN2eWnUTsB2Uijj0dExlFf4GM8msjzGnPM/wDuNMbU1g/QbK0xx+sEhgDjgRBgtYisMcZs93Q4D2nMMZ8DrAfOALoBC0VkhTGmyMPZrOTW/vK2Am/MLj+N2gnIRhp1PCJyCvAyMMEYc6iFsnlKY445GXi3obyjgYkiUmOM+bhFErpXY/9eHzTGlAKlIrIcGADYtcAbc8zTgUdN/eTwThHJBHoBa1smoiXc2l/eNoXSmF1+PgWuaXg3dwRQaIzJaemgbnTCYxaRzsB84Gobj8iOdsJjNsZ0McYkGmMSgQ+Bm21a3tC4v9efAKeKiFNEWlG/QfiWFs7pTo055mzq/8WBiLQDegK7WjRly3Nrf3nVCNwcZ5cfEbmp4ftzqD8jYSKwEyij/re4bTXymO8DooDnGkakNcbGK7k18ph9RmOO1xizRUS+AtKAOuBlY8wxT0Wzg0b+GT8EvCYi6dRPLdxpjLH1ErMi8g4wFogWkb3An4FA8Ex/6aX0SillU942haKUUqqRtMCVUsqmtMCVUsqmtMCVUsqmtMCVUsqmtMCVUsqmtMCVUsqm/h9LSCxfFofRhgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x=np.linspace(x1,x2,100)\n",
    "sol=Analytical(x)\n",
    "\n",
    "plt.plot(x,sol,label='analytical')\n",
    "plt.plot(mymesh.nodecoords,phi,'*',label='numerical')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "dddc3e36-01f1-4565-8101-bd0eb1c4735a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L2 error is: 7.53331e-14\n"
     ]
    }
   ],
   "source": [
    "sum=0.0\n",
    "for i in range(nDofs):\n",
    "    sum+=(phi[i]-Analytical(mymesh.nodecoords[i]))**2\n",
    "\n",
    "print('L2 error is: %g'%(np.sqrt(sum)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2c809f0b-6821-4ff8-a515-6f100f7ce407",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "save result to result.jpg\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEKCAYAAAAW8vJGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAASnElEQVR4nO3df6zd9V3H8efLItFFI3NsCylgKzK1UWi2OyRGXd38QfnDarI/oMuIuKVBhzExZhD/cH/MROePaBbZmmYhuD+wMUq0m53EaAAN4iixY5SFee0iXLtkRR375Y+Uvf3jHOR8rrftbfl+v/d77/f5SJr0nO/3vnh/W3JfPefc8z6pKiRJesk3bPQAkqRxsRgkSQ2LQZLUsBgkSQ2LQZLUuGSjB3ilLr/88tqxY8dGjyFJm8oTTzzxfFW9dq1jm74YduzYwbFjxzZ6DEnaVJL8y9mO+VSSJKlhMUiSGhaDJKlhMUiSGhaDJKkxWDEkuTfJF5I8dZbjSfLBJMtJnkzyxqFmkyS9bMhHDPcBN53j+F7g2vmvA8CHB5hJkrTKYO9jqKpHkuw4xyn7gI/WbA/4Y0kuS3JFVX2+j3kOHTrE/fff32nmqVOnSMIVV1wxibw+MqeW10fm2PP6yJxa3kv279/PgQMHOs2Ecb3BbTvw3MLtlfl9/68Ykhxg9qiCq6+++qL+Y/fffz/Hjx9n9+7dF/X1azl58iRAZ3/5Y8/rI3NqeX1kjj2vj8yp5QEcP34coJdiyJAf1DN/xPDxqvq+NY79BfAbVfV389t/Dby3qp44V+bS0lJdzDuf9+zZA8BDDz10wV97NpdddhkAX/ziFyeR10fm1PL6yBx7Xh+ZU8uDV/49LMkTVbW01rEx/VTSCnDVwu0rgVMbNIskTdaYiuEIcNv8p5NuBF7o6/UFSdLZDfYaQ5I/AvYAlydZAd4HfCNAVR0EjgI3A8vA14Dbh5pNkvSyIX8q6dbzHC/gPQONI0k6izE9lSRJGgGLQZLUsBgkSQ2LQZLUsBgkSQ2LQZLUGNOupEE9/PDDwMtvVe/CCy+80Gnm2PP6yJxaXh+ZY8/rI3NqeQBf+cpXuOaaazrLW+QjBknapPradTfZRwxvectbAJfojS1zanl9ZI49r4/MqeXBy0v0+uAjBklSw2KQJDUsBklSw2KQJDUsBklSw2KQJDUsBklSw2KQJDUsBklSw2KQJDUsBklSY7K7ktyu2o2xzzj2vD4yx57XR+bU8sDtqpKkNbhdtWNuV+3G2Gcce14fmWPP6yNzanngdlVJ0oAsBklSw2KQJDUsBklSw2KQJDUsBklSw2KQJDUsBklSY9BiSHJTkmeSLCe5e43j35bkY0k+leREktuHnE+SNGAxJNkG3APsBXYBtybZteq09wBPV9X1wB7gd5NcOtSMkqRhV2LcACxX1UmAJIeBfcDTC+cU8K1JAnwL8O/AmT6GcYleN8Y+49jz+sgce14fmVPLg62zRG878NzC7ZX5fYv+APhe4BTwaeCXqurrq4OSHEhyLMmx06dP9zWvJI3aVliilzXuW31VPwkcB94KXAP8VZK/raovNV9UdQg4BLC0tHRRfzIu0evG2Gcce14fmWPP6yNzanmwdZborQBXLdy+ktkjg0W3Aw/UzDLwOeB7BppPksSwxfA4cG2SnfMXlG8Bjqw651ngbQBJXg98N3BywBklafIGeyqpqs4kuRN4ENgG3FtVJ5LcMT9+EHg/cF+STzN76umuqnp+qBklSQN/UE9VHQWOrrrv4MLvTwE/MeRMkqSW73yWJDUsBklSw2KQJDUsBklSw2KQJDUsBklSY9AfVx0Tl+h1Y+wzjj2vj8yx5/WRObU82DpL9CRJHdoKS/RGxSV63Rj7jGPP6yNz7Hl9ZE4tD7bOEj1J0iZgMUiSGhaDJKlhMUiSGhaDJKlhMUiSGhaDJKlhMUiSGhaDJKlhMUiSGpNdieESvW6Mfcax5/WROfa8PjKnlgcu0ZMkrcEleh1ziV43xj7j2PP6yBx7Xh+ZU8sDl+hJkgZkMUiSGhaDJKlhMUiSGhaDJKlhMUiSGhaDJKlhMUiSGoMWQ5KbkjyTZDnJ3Wc5Z0+S40lOJHl4yPkkSQO+8znJNuAe4MeBFeDxJEeq6umFcy4DPgTcVFXPJnndUPNJkmaGXIlxA7BcVScBkhwG9gFPL5yzH3igqp4FqKov9DWMS/S6MfYZx57XR+bY8/rInFoebJ0letuB5xZur8zvW/QG4NVJHkryRJLb1gpKciDJsSTHTp8+3dO4kjRuW2GJXta4b/VVXQK8CXgb8M3A3yd5rKo+23xR1SHgEMDS0tJF/cm4RK8bY59x7Hl9ZI49r4/MqeVBv0v0hiyGFeCqhdtXAqfWOOf5qvoq8NUkjwDXA59FkjSIIZ9Kehy4NsnOJJcCtwBHVp3z58APJ7kkyauAHwA+M+CMkjR5gz1iqKozSe4EHgS2AfdW1Ykkd8yPH6yqzyT5S+BJ4OvAR6rqqaFmlCQN/EE9VXUUOLrqvoOrbv828NtDziVJepnvfJYkNSwGSVLDYpAkNSwGSVLDYpAkNSwGSVJj0B9XHROX6HVj7DOOPa+PzLHn9ZE5tTzYOkv0JEkd2gpL9EbFJXrdGPuMY8/rI3PseX1kTi0P+l2i5yMGSVLDYpAkNc5ZDEluHGoQSdI4nO8RwyNJ3p9ksq9FSNLUnK8Y9gLvBD6ZZNcA80iSNtg5i6Gq/hr4fuAfgWNJfnmQqSRJG+a8Lz5X1Zer6l3MHjn8VpKvJPnS4q/+x5QkDWVdrx0kWQJ+Hfgn4HeAM30OJUnaOOcshvmLzu8D7gI+BNxdVf81xGCSpI1xvkcMjwPfDuydv94gSdrizvcaw1PAdZaCJE3HOR8xVNU7hxpkaG5X7cbYZxx7Xh+ZY8/rI3NqeeB2VUnSGtyu2jG3q3Zj7DOOPa+PzLHn9ZE5tTxwu6okaUAWgySpYTFIkhoWgySpYTFIkhoWgySpYTFIkhoWgySpMWgxJLkpyTNJlpPcfY7z3pzkxSRvH3I+SdKAxZBkG3APs48L3QXcutbHhc7P+wDw4FCzSZJeNuRKjBuA5ao6CZDkMLAPeHrVeb8I/Cnw5j6HcYleN8Y+49jz+sgce14fmVPLg62zRG878NzC7ZX5ff8nyXbgZ4CD5wpKciDJsSTHTp8+3fmgkrQZbIUlelnjvtVX9fvAXVX1YrLW6fMvqjoEHAJYWlq6qD8Zl+h1Y+wzjj2vj8yx5/WRObU86HeJ3pDFsAJctXD7SuDUqnOWgMPzUrgcuDnJmar6s0EmlCQNWgyPA9cm2Qn8K3ALsH/xhKra+dLvk9wHfNxSkKRhDVYMVXUmyZ3MftpoG3BvVZ1Icsf8+DlfV5AkDWPQD+qpqqPA0VX3rVkIVfWzQ8wkSWr5zmdJUsNikCQ1LAZJUsNikCQ1LAZJUsNikCQ1Bv1x1TFxiV43xj7j2PP6yBx7Xh+ZU8uDrbNET5LUoa2wRG9UXKLXjbHPOPa8PjLHntdH5tTyoN8lej5ikCQ1LAZJUsNikCQ1LAZJUsNikCQ1LAZJUsNikCQ1LAZJUsNikCQ1LAZJUmOyKzFcoteNsc849rw+Msee10fm1PLAJXqSpDW4RK9jLtHrxthnHHteH5ljz+sjc2p54BI9SdKALAZJUsNikCQ1LAZJUsNikCQ1LAZJUsNikCQ1LAZJUmPQYkhyU5JnkiwnuXuN4+9I8uT816NJrh9yPknSgMWQZBtwD7AX2AXcmmTXqtM+B7ylqq4D3g8cGmo+SdLMkCsxbgCWq+okQJLDwD7g6ZdOqKpHF85/DLiyr2FcoteNsc849rw+Msee10fm1PJg6yzR2w48t3B7ZX7f2bwL+MRaB5IcSHIsybHTp093OKIkbR5bYYle1rhvzatK8qPMiuGH1jpeVYeYP820tLR0UX8yLtHrxthnHHteH5ljz+sjc2p50O8SvSGLYQW4auH2lcCp1ScluQ74CLC3qv5toNkkSXNDPpX0OHBtkp1JLgVuAY4snpDkauAB4J1V9dkBZ5MkzQ32iKGqziS5E3gQ2AbcW1UnktwxP34Q+DXgNcCHkgCcqaqloWaUJA38QT1VdRQ4uuq+gwu/fzfw7iFnkiS1fOezJKlhMUiSGhaDJKlhMUiSGhaDJKlhMUiSGhaDJKkx6PsYxsTtqt0Y+4xjz+sjc+x5fWROLQ+2znZVSVKHtsJ21VFxu2o3xj7j2PP6yBx7Xh+ZU8uDfrer+ohBktSwGCRJDYtBktSwGCRJDYtBktSwGCRJDYtBktSwGCRJDYtBktSwGCRJjcmuxHCJXjfGPuPY8/rIHHteH5lTywOX6EmS1uASvY65RK8bY59x7Hl9ZI49r4/MqeWBS/QkSQOyGCRJDYtBktSwGCRJDYtBktSwGCRJDYtBktSwGCRJjUGLIclNSZ5Jspzk7jWOJ8kH58efTPLGIeeTJA1YDEm2AfcAe4FdwK1Jdq06bS9w7fzXAeDDQ80nSZoZciXGDcByVZ0ESHIY2Ac8vXDOPuCjNVsA8liSy5JcUVWf73qYl5bodfm28pcWZXWVOfa8PjKnltdH5tjz+sicWh7A8ePH2b17d2d5i4Z8Kmk78NzC7ZX5fRd6DkkOJDmW5Njp06c7H1SSxm737t3s37+/l+whHzFkjftWrwZczzlU1SHgEMDS0tJFrRfsayuhJG12Qz5iWAGuWrh9JXDqIs6RJPVoyGJ4HLg2yc4klwK3AEdWnXMEuG3+00k3Ai/08fqCJOnsBnsqqarOJLkTeBDYBtxbVSeS3DE/fhA4CtwMLANfA24faj5J0sygH9RTVUeZffNfvO/gwu8LeM+QM0mSWr7zWZLUsBgkSQ2LQZLUsBgkSY1s9jd6JTkN/MtFfvnlwPMdjrMZeM3T4DVPwyu55u+oqteudWDTF8MrkeRYVS1t9BxD8pqnwWuehr6u2aeSJEkNi0GS1Jh6MRza6AE2gNc8DV7zNPRyzZN+jUGS9P9N/RGDJGkVi0GS1JhEMSS5KckzSZaT3L3G8ST54Pz4k0neuBFzdmkd1/yO+bU+meTRJNdvxJxdOt81L5z35iQvJnn7kPP1YT3XnGRPkuNJTiR5eOgZu7aO/7e/LcnHknxqfs2bektzknuTfCHJU2c53v33r6ra0r+Yrfj+Z+A7gUuBTwG7Vp1zM/AJZp8gdyPwDxs99wDX/IPAq+e/3zuFa14472+Ybfl9+0bPPcDf82XMPlf96vnt12303ANc868CH5j//rXAvwOXbvTsr+CafwR4I/DUWY53/v1rCo8YbgCWq+pkVf0PcBjYt+qcfcBHa+Yx4LIkVww9aIfOe81V9WhV/cf85mPMPi1vM1vP3zPALwJ/CnxhyOF6sp5r3g88UFXPAlTVZr/u9VxzAd+aJMC3MCuGM8OO2Z2qeoTZNZxN59+/plAM24HnFm6vzO+70HM2kwu9nncx+xfHZnbea06yHfgZ4CBbw3r+nt8AvDrJQ0meSHLbYNP1Yz3X/AfA9zL7WOBPA79UVV8fZrwN0fn3r0E/qGeDZI37Vv+M7nrO2UzWfT1JfpRZMfxQrxP1bz3X/PvAXVX14uwfk5veeq75EuBNwNuAbwb+PsljVfXZvofryXqu+SeB48BbgWuAv0ryt1X1pZ5n2yidf/+aQjGsAFct3L6S2b8kLvSczWRd15PkOuAjwN6q+reBZuvLeq55CTg8L4XLgZuTnKmqPxtkwu6t9//t56vqq8BXkzwCXA9s1mJYzzXfDvxmzZ6AX07yOeB7gE8OM+LgOv/+NYWnkh4Hrk2yM8mlwC3AkVXnHAFum7+6fyPwQlV9fuhBO3Tea05yNfAA8M5N/K/HRee95qraWVU7qmoH8CfAL2ziUoD1/b/958APJ7kkyauAHwA+M/CcXVrPNT/L7BESSV4PfDdwctAph9X5968t/4ihqs4kuRN4kNlPNNxbVSeS3DE/fpDZT6jcDCwDX2P2L45Na53X/GvAa4APzf8FfaY28WbKdV7zlrKea66qzyT5S+BJ4OvAR6pqzR973AzW+ff8fuC+JJ9m9jTLXVW1addxJ/kjYA9weZIV4H3AN0J/379ciSFJakzhqSRJ0gWwGCRJDYtBktSwGCRJDYtBktSwGCRJDYtB6kiSb0jySJLVbyZ81XxN9Ic3ajbpQlgMUkfmi9p+Fnhrkp9bOPQBZm8m/ZWNmEu6UL7BTerY/F24vwV8P/BdzN6lu6eq/m5DB5PWyWKQepDkQWbbTHcAh6vqvRs7kbR+FoPUgyQ7mX3S2D8D31dV/73BI0nr5msMUj9+DvhPZiuQv3ODZ5EuiI8YpI4leTPwKPBTwM8Drwd+sKpe3NDBpHXyEYPUoSTfBHwUuK+qPgEcYPYCtK8xaNPwEYPUoSS/B/w0cF1VfXl+3y3AHwJv2syfhaDpsBikjiT5EeBvgB+rqodWHftjZq813FhVZzZgPGndLAZJUsPXGCRJDYtBktSwGCRJDYtBktSwGCRJDYtBktSwGCRJDYtBktT4XwqWQF3kgyKlAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbcAAAEoCAYAAADbp799AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAz40lEQVR4nO2de9Bk5V3nP7+5RcgMMMwMKJfJgAKGZYMmbwJrjEvuMGLYWHETchMSl8IAoqtrom7kjy23ko3rBiuJUyOZUFQoMBJMWB0hqVJES0CGSLhKMksE3kzMzHDLABPm9ts/uht6+u3u85xznufc+vupmpq3u59zznO6T59P/37PzdwdIYQQokssqrsCQgghRGwkNyGEEJ1DchNCCNE5JDchhBCdQ3ITQgjROSQ3IYQQnWNJXQdevXq1r1u3rq7DCyFEZ7j77rt3uvuaWPt7yxt/xJ948kCube65d+8t7n52rDqUpTa5rVu3ji1bttR1eCGE6Axm9mjM/T3x5AH+9q+PzrXNymPnV8esQ1mUlhRCCNE5JDchhBCdI1NuZrbJzLab2f0TXjcz+2Mz22pm95rZq+NXUwghhAgnJHK7GpjWSHgOcFL/30XAn5SvlhBCCFGcTLm5+23Ak1OKnAdc4z3uAI4wsx+LVUEhhBAiLzF6Sx4LPD70eL7/3PdGC5rZRfSiO9auXRvh0NXx9pd/sO4qCCEq4pbnrqm7CqIkMTqU2Jjnxq6j4+4b3X3O3efWrIk2JEMIIYQ4iBiR2zxw/NDj44BtEfYrOsqiVUfWXYVOcuCJaa0HQswWMeR2E3CpmV0PnAE84+4LUpIiLhKEGKXN14TE3H3MbBNwLrDd3U/rP3c6sAFYDvwr8D53/8HIdscD1wA/ChwANrr7lVnHCxkKcB1wO3CKmc2b2YfN7GIzu7hfZDPwCLAV+FPgIyEnKvKzaNWRL/4Tokvo2p4JrmZhz/urgI+5+78H/gL4b2O22wf8pru/EjgTuMTMTs06WGbk5u7nZ7zuwCVZ+xHF0Jc9HF+5ou4qBGFP7aq7Co1m+JpXRNcd3P02M1s38vQpwG39v78O3AJ8fGS779HvoOjuu8zsIXqdFh+cdrza5pYU0+mC1Noim6qp+n1ps0wH3wNJrrPcD7wD+CrwSxzcf2MBfTn+NHBn1o4lt4pom6wkpu5Q1WeZUqKxvz+S5XR+6Ma39y3Ou9lqMxueDX+ju2/M2OZDwB+b2e/T67+xZ1JBM1sOfBn49dF2uXFIbhWQUmySkGgKRa7FNkeVYgE73X0uzwbu/i/A2wDM7GTg58eVM7Ol9MR2rbvfGLJvyS0xKcTWBaHtXXlI3VXoFEuf2l13FQoxfC1LdLOHmR3l7tvNbBHw3+n1nBwtY8DngYfc/Y9C9y25JSSm2JoiNEmpmaT8XKoSp0TXbfo978+il76cB64AlpvZoEPijcAX+mWPAa5y9/XA64EPAPeZ2T39sr/r7punHU9yS0QMsdUlNAlMDDPtekglPomue0zpeb9gzJq7bwPW9//+B8bPhDUVyS0BZcUWQ2pdFNTew5fWXYVKWPrM3rqrEEzIdVZWgIPvgyQn8iC5Raao2LoqtFkRUkxivmdNEOXodVlUdormRB4kt4gUEVsZqVUtM4mqfRT9zFJKMYbsFM2JLCS3SOQRW5OF1hWBvXB47jE6SXnZM/vrrkIupl0HscU3fE3nFZ2iOTEJyS0CoWIrKrXYQmuKwJomoJSkPNeqxTnp+okhvRiik+QESG6VUURsMaVWtdBmSVx1E/pep5bg8DUWS3RFU5ZZglu06kjNUtJxJLeSpJp9JIbYUgpN8mof0z6z2OIbXHtlJVdUcCFIcN1GcquAPFFb06TWFontWZF7GEytLNs1drH62hj+nGOKLkY0V0RwIdGb6DaSWwlCorZQsZWVWlmhNU1ibZNVXvKeX5UyHHctxBBeGdGlEpyit+4iudVMnVKrQ2hdl1Yq6pZhbOEVSVumTFGKg9njS3hkz5qcWz2apC5FkdwSkmr6rDYITRKrl6z3P4b8YggvbzSXV3BKT84uklsiQsSWN2prktSaLK89zZhjegHLGnSPHff5xRJe0YguNJqT4EQIkltN5BFbne1pTZBYU2WVlzznUYcIhz/rMqIbXG9lJBdbcGL2kNwSECsdWYfUqpJZV4SVirzvT2wZxhBdGcnFFpyit9lDcotMjHRk1VJLLTSJLD1Z73EZ+ZUVXVHJSXCiDJJbh2iC1Jomsr3L6z3+0mfrPf6A0c+lqOxGr5k8sisiOaUoRVEkt4pJEbXlkVoModUlsLpFVYQida5CiOM+wyLCKyK7vJKT4EQRJLcWU4XUqhRZG+WVgqz3IZX8hj/rspFdbMlJcCIvkluFxIra2io1ySsOIe9jWQGWFV0KyUlwIg+SW4tI3aYWU2hNFNm+ww5UcpwlP1hUyXGmMfr+l5FdGdHtWWHB7XIhkgsRnGgmZrYJOBfY7u6n9Z/7M+CUfpEjgKfd/afGbHsEcBVwGuDAh9z99mnHk9wqomzUljJaiyG1qmRWlaDKULSOKaU47vMpIrwiossTxUH2QPAswSl6ayxXA58Brhk84e7vHvxtZv8beGbCtlcCN7v7u8xsGXBo1sEktwYQS2x5pFZGaKlF1gaBpSDrvGPLr2x0N7iGUkiurOBE83D328xs3bjXzMyA/wy8acxrhwE/B1zQ388eYE/W8SS3jlCF2GJLbVYlVpRJ71cs6e1dXi6ayyO5pi35I2rnDcD33f3bY147EdgBfMHMTgfuBi539+em7VByq5myUVvbojUJLT6j72kZ2Q0+59SSC4niFL3Vxw8PLGXrC0fn3Wy1mW0ZerzR3TcGbns+cN2E15YArwYuc/c7zexK4GPAx6ftUHKrkTaILYbU6hTa/uVxV5jOw+Jn61kjL4bsykouVhQnwbWKne4+l3cjM1sC/CLwmglF5oF5d7+z//gGenKbiuRWAUXWbIsltjqkVoXM6pRWKHnqmFKEw59HXtEVlVzeKK6M4MahTiWt4i3Av7j7/LgX3f3fzOxxMzvF3R8G3gw8mLVTya0mys4fGUIRsRWVWkqhtUFkZZl2jjHFV1R0ZSTXpKV+RH2Y2XXAWfTSl/PAFe7+eeA9jKQkzewY4Cp3X99/6jLg2n5PyUeAC7OOJ7k1kBhRW16xFZFaCqHNgsjyMu49iSG8IqIrIrkQwZWJ3pSabAfufv6E5y8Y89w2YP3Q43uAXClPya1hlBVbaqnFEpokVo7R96+s7AafayrJhaQpU6QnxewiuTWIKsVWR6RWp9CWLs8cFhOFvc8uq+Q4o8SS3b7DDiRNV2ZFcRomIGIhuSVmXGeSFO1toWLrmtSqklYoeeqTUoRlZJc3ioNyvSvzMCl6G5eaVKeS2UZyawkxlqqpMgWZQmhNE1lZss4npvwGn0cVkssSnKI3UQWSW0RCVuEuQtl0ZFXRWkyhdU1kRRh9D2LIrqjkuiI4rcY9O0huFTMpJTmpva1qseWVWiyhSWbZxJTd8OcWIroUHU6KDhNQxxIRguRWkEWrjswsU2TwdkzyiC2P1MoKrU6RHb4iXRvMM7uq/bxjyS5PNFdEckXb4fJGbxoSIIYJkpuZnU1vyYHF9AbWfWLk9cOBLwJr+/v8Q3f/QuS6ihFirb8WKrYyUqtCaCnFFev4KQU4/B4XEV1eycWYsDn1IO+QTiWLVh3JgSeeTFcJUQuZcjOzxcBngbfSm+PrLjO7yd2Hpz+5BHjQ3X/BzNYAD5vZtf2lCUQGRVKSMdKRVURrsaVWt8DKklX/WPIrI7r9y/dHjeKyordpglPnElGUkMjtdcBWd38EwMyuB87j4Lm9HFjRX5NnOfAksC9yXRtNSGeSWEMAyoottdRiCq3tMsvL6PnGkN3g88gjudhRXBnBjSPPkIBR1Kkkmz2+hEd3r6q7GqUIkduxwONDj+eBM0bKfAa4CdgGrADe7e4L7qBmdhFwEcDatWuL1FeUJKXYYkht1mSWRUzZFZVcirY4IVITciWOy42N5gneDtwDHAP8FPCZ/uqpB2/kvtHd59x9bs2aNTmr2i5SdSZJ0e1/HHnEtnT5nlJiO3zF7hf/ienEeK/yfl6xesRmXZuTru0YYzzF7BESuc0Dxw89Po5ehDbMhcAn3N2BrWb2HeAngX+KUsuGEdJTsiwpvtAhUVteqRWlLpEd9fK0U2hsfy7ycuVTGH4Pi0R0eSK50DRlrI4mIWhIgJhGiNzuAk4ysxOA79JbnuC9I2Ueo7fGzt+b2dHAKfSWJRAZZM0nOUzqqK0KsVUhtdQCK3rslOIbvK9FJReaqgxNU04idtubEJPIlJu77zOzS4Fb6A0F2OTuD5jZxf3XNwD/A7jazO6jl8b8qLvvTFjv1lHF+m3TyIraQsXWFKnVKbCiTKpzTOkVlVxMwaWI3ib1mlT0JiYRNM7N3TcDm0ee2zD09zbgbXGr1h5STbuVh2lRW51iiym1NgothNHziiG7IpLLm6YsKrgyA7tD0GBuAZqhpDXEGrA9SiqxxZJaV4U2jeFzLiu6opKLIbiiTEpNFh3zptUBZhPJrYHk7UxSNGpLIbYYUptFoU0ilujySi6G4OqM3oSQ3BIQOgwgT2eSicfqiNgktGxiiO7wFbujC06IJiK5tYBUKclpdFFqxx/6dCXHefz5I5IfY/C+FZFcXsHB9Ha4FNFb7NSkqB8z2wScC2x399OGnr8MuJTerFZ/5e6/PWbbqfMbj0Nyaxh5UpKporbUYksttKoklvf4KaRXNJorM3RgHEUFFwP1mGwNV9ObzeqawRNm9kZ60zm+yt1fMLOjRjcKnN94AZJbBdQ9DGCYmAuK5hVbbKnVLbG8pJZekWguNIoLSVEW6WBSJHoT7cTdbzOzdSNP/yq9CUBe6JfZPmbTkPmNFyC51URoe9uklGSZrv/TCIna6pJa22QWyuh5lZVdXsmFRnFl2uDUuURM4GTgDWb2B8APgd9y97tGyoTMb7wAya1BxJhyK3U6sg6xdVVqkxicbwzJxY7isgSXanhAXjTWrRx7Diwucv2tNrMtQ483uvvGjG2WACuBM4HXAl8ysxP7UzkOCJnfeOyORcsoMs1W28Q2a0Ibx/B7UFR0TRJcrQO7NdatCna6+1zObeaBG/sy+yczOwCsBnaMlMma33gBWp+iweTtJVk0HRlbbEe9/NnCYjv+0KcltjGUeV/yfh4hn3XRadjyXqPjvgMxMhxNmFVIAPAV4E0AZnYysAwYnbrxxfmNzWwZvfmNb8rasSK3GhjX3hb6hU0RtWWRV2x5qUJmrzjkieTHGCbVQo9lork8UVye4QLjiN25JIuyPSYXrTqSA088WXh7kY2ZXQecRS99OQ9cAWwCNpnZ/cAe4Jfd3c3sGHpd/tdPmt8463iSW0dIFbWlFFtKqVUts6zjp5Bdkba5vGnKaRTtYFLlsjiiObj7+RNeev+YstuA9UOPF8xvnIXk1nGmRW0xVs4eULfY6pZZFilld/yhTycRXNn2t6Z0LhGzieTWUMa1NcRaZTuU0KitDrE1XWZZxJZdEcFB9nCBsunJSSh6E6nR1dUAyjaQF0lJxkpHSmxxiHFOqTrjZF0L066lvO29Vf+AE91FchOFqVpsrzjkiU6KbUCs88vzPmvCatFVJLeKibESQCiTfjXHiNqqFFvXpTZKjPONLbgy0dskQjMOdUwcLtqP5NYSqkrXNFFss0rbBDeJmPOZDhiXyq/yh6NoPpJbYopMmpznl+qkX79Fo7YsqhLbrEVrkyj7PlSZoozV+zbWD7kmTVguqke9JcWLZP0yr0JsqYX2Ey/7ftL9b33h6CT7HbwvRXpV5u1JOY2ivSc1LEBUjeRWMzGmEgol5ri2LPKKLZXUUsss63ixZVdUcqGCKzvIO+/Abg0JaCb79i+KNti/LiS3FpM3JTmNmFFbHrHFllrVMssilexeccgTtQku1di3aWhtN5EXya0F1D32p+lia5rQpjFc17KiKxLFVRHBTYrelJoUVSK5RWbvymp/0YYyLSVZtBfcKKFiiyG1NgltErFEl1dyMdrgYkZv41KTWsBUlEXJ7oYR2lMyZkpyGqFRW1Vi+4mXfb8TYhslxnnleW9DPq+mD/DO217d1B+eIg2SW4XUNQ4nddRWpdi6TpWCK8u0a6fKzktCjENy6xB1RG1ViK2r0dokyp5v6HtdR/RW5hoNyWpoILcYILm1kKJrt41SNmpLLbZZk9ooVZx72fk+Y7XXxrqmhRgguTWcunpKZv1iT7169ixLbZiigo/d/iZE21BvyZL4ymbP6lp320eRqC212E5ctiPJfh/ZsybJfqH3nuTtUVlkLNw4ig4LKLpS94CQHpN7VhjLdnnhY4juosgtIVlz28WcnSRvW0aZdFLKdGRKsZ24bEcysVWx/5TvTZnoLe+1lGIiZSFGkdzEAurqAp7q5p1aOlUeL+97pMmnJ9P0rEvXMLNNZrbdzO4f89pvmZmb2eop2y82s382s78MOZ7k1jLqbnhPFbWlEFvVUpt0/Nh1SPUjYNpn2/Qxb5PQygCN4mrg7NEnzex44K3AYxnbXw48FHowya3DTGpvi9XDLRaxb9Z1S20cseuU5z2rInqbdE3lafMN+eGmhUvbi7vfBjw55qX/A/w2MLHx1MyOA34euCr0eOpQ0gGa1oaR52YaS2xNk9kkBvWM0fmkSCcTIULYf2BRkenVVpvZlqHHG91947QNzOwdwHfd/ZtmU/sgfJqeAIN/3khu4iCmpZ9CUpJVi60tUhslluRCBRfac3LavJNll8MZRpMod5Kd7j4XWtjMDgV+D3hbRrlzge3ufreZnRW6f6UlG0QVKZempCTLiq2JqcciVHkeqdOTKa6tulfEEEn5ceAE4Jtm9q/AccA3zOxHR8q9HnhHv8z1wJvM7ItZO5fcOkod49uq6pnXBamNUuacmj7gve6xlqKZuPt97n6Uu69z93XAPPBqd/+3kXK/4+7H9cu8B/gbd39/1v4ltwZT9a/WsinJUMrcjLsotgFVCC7kB0gXe02K+jGz64DbgVPMbN7MPjyl7DFmtrnM8dTm1iLqHgYwjdCoTWKbzonLdiSd6aTpjFvbTXQDdz8/4/V1Q39vA9aPKXMrcGvI8XQVVUQTZiuf1iZSVdRWlFkQ24Ci5xozepvGtGulKW26QgTJzczONrOHzWyrmX1sQpmzzOweM3vAzP4ubjXFJJowDCB11DZLYhuQWnBZVPGDpgnXrugumXIzs8XAZ4FzgFOB883s1JEyRwCfA97h7v8O+KX4VRWzyCyKbUDKc+/ytFxNyJKI+glpc3sdsNXdHwEws+uB84AHh8q8F7jR3R8DcPftsSvaNWJOmjxK3t5pZVKSKaO21GI7aUmcyOHb+9LdTIu0wcUa+xZzzFvZFQKmoZUBxDhC0pLHAo8PPZ7vPzfMycBKM7vVzO42sw+O25GZXWRmW8xsy44ds/uLfByaVqg6TlqyP5rYBvsTQjSLELmNCzFGfyYtAV5Db+6vtwMfN7OTF2zkvtHd59x9bs2a2e0RVgd1NvQ3KWpLJaKUgpvl1KwQRQmR2zxw/NDj44BtY8rc7O7PuftO4Dbg9DhVFE2lbe02qSOsJkVwVfWcFKKphMjtLuAkMzvBzJbRGyF+00iZrwJvMLMl/fnCziDH0gRiIbM87VCKSKUq8aQ6Tpuit9RZAqXwRQiZHUrcfZ+ZXQrcAiwGNrn7A2Z2cf/1De7+kJndDNwLHACucvcFC9KJ5jELM05UHVGdtGR/0k4mVVLVRMqiWfh+S9YBqCqCZihx983A5pHnNow8/hTwqXhVE1mkHicUY6xT3va22BFKXanCFIKb9dlLhMiDZihpCU2beqsNbTV1t4HVfXy1u4lZRnLrGG2egT1m1Fa3WAbErkeb2t5CCMk+zHL7syiO5CY6R1PENqBp9Wk6KbMUe1fmXl1atBTJTSSjjnXGmiqSmPXKE701fa03IVIhuYmxTOtMkqKNJka6raliG9Dk+k37TJuwKoQQeZHcZoBJ4466NAygyeIYpi31DCXPNdTm9mDRPiQ3kYQ86bCyUVvbhBGjvk3vWKJ13UTdSG4RUWN19bRNbAOqrLfa3cQsIrmJWikTgbRVbAPK1r/p0ZsQdSK5iVw0ZcBv28U2oEnn0ZTPVnQTM9tkZtvN7P6h5z5lZv9iZvea2V/0F74et+1vmNkDZna/mV1nZj+SdTzJTSygbO+40DSYIo8eZQQX+h4qNSkawNXA2SPPfR04zd1fBXwL+J3RjczsWODXgDl3P43eHMfvyTqY5CaEECI57n4b8OTIc19z9339h3fQW1JtHEuAQ8xsCXAoC5ddG7uBEEII8RIHjMXP5p74e7WZbRl6vNHdN+bY/kPAn40+6e7fNbM/BB4DdgNfc/evZe1MchNCCBGDne4+V2RDM/s9YB9w7ZjXVgLnAScATwN/bmbvd/cvTtun0pJCCCFqw8x+GTgXeJ+7+5gibwG+4+473H0vcCPwM1n7ldyEEELUgpmdDXwUeIe7Pz+h2GPAmWZ2qJkZ8Gbgoax9S25CCCGSY2bXAbcDp5jZvJl9GPgMsAL4upndY2Yb+mWPMbPNAO5+J3AD8A3gPnreymzLU5ubEEKI5Lj7+WOe/vyEstuA9UOPrwCuyHM8RW5CCCE6h+QmhBCic0huQgghOofkJoQQonNIbkIIITqH5CaEEKJzSG5CCCE6h+QmhBCic2gQtxBCiIOw/bDkB+2OfdpdeyGEEGIMkptYwOPPH1HJcR7Zs6aS43SZ0Pdw6wtHJ66JEM1CchPRqeJG+u19uRdSbCxNOZdHd6+quwrJWfrU7rqrICpCchO5aNINsClSKEMXzkGIJiK5iVopm5pssxzK1l1pXSEmI7lFRCmPl6iyjaeNgquyzmpvE7OI5DbDbH9ued1VAOJEIG0SXIy6NiVqm3QNPbPrkIprIsTBSG4zQJEbTVU9JmPSBsE1tY7T2lLbeC0IIbmJ3KToVBIrEmmqPCBe3fK8V0pJillFchPJqOvG2kTBNbFOVbP32WVB5do+M4ZoBrqKOkboDaTrNEkmMevSlLa2Kln6bN01EDEws01mtt3M7h967pfM7AEzO2BmcxO2O97M/tbMHuqXvTzkeJJbS5iFX7Oxb9xNEFwT6pBFk8Yuik5zNXD2yHP3A78I3DZlu33Ab7r7K4EzgUvM7NSsg3X/jimSEHpDrLvNp0651C22GO996s4ki59tvvxFHNz9NuDJkececveHM7b7nrt/o//3LuAh4Nis40luLSbGjWHacIA6esmlSLvVIZkUx2xaSrIpQ0lEfOxALx2c5x+w2sy2DP27KHq9zNYBPw3cmVU2SG5mdraZPWxmW83sY1PKvdbM9pvZu4JrK8aidoa4VCm4uiO2rrNsV901EBPY6e5zQ/82xty5mS0Hvgz8urv/IKt8ptzMbDHwWeAc4FTg/HH5zn65TwK35K20EMOkilCqkE6qY6R6T1K0t2kAt4iNmS2lJ7Zr3f3GkG1CIrfXAVvd/RF33wNcD5w3ptxl/YNvD6yvqJAUN5y2tLsNk1JwTYrYmvSeC1EGMzPg88BD7v5HoduFyO1Y4PGhx/OMNOaZ2bHAO4ENGZW8aJCP3bFjR2gdZwKlWg4mZftSCgmlFFvT2tqEKIKZXQfcDpxiZvNm9mEze6eZzQP/AfgrM7ulX/YYM9vc3/T1wAeAN5nZPf1/67OOtySkTmOe85HHnwY+6u77e5IdTz8HuxFgbm5udB8iEnufXcbS5XvqrkZpHtmzhhOXpfkR1KQoaxopxZYVeU/rUJS3M4nGXwp3P3/CS38xpuw2YH3/739gvIemEhK5zQPHDz0+Dtg2UmYOuN7M/hV4F/A5M/tPeSszSyzb1Ry3l+kxmTo1OctRS9Fzn7WUZJO+S6I5hMjtLuAkMzvBzJYB7wFuGi7g7ie4+zp3XwfcAHzE3b8Su7JiIbMwTmgWBZdabE0YuJ3q2n3ZM/uT7Fe0i0y5ufs+4FJ6vSAfAr7k7g+Y2cVmdnHqCnaFJnzhpnUqaXL0BrMpuDopmpJUT0nRFELa3HD3zcDmkefGdh5x9wvKV0u0jUd3r+IVhzyRWW7rC0fzEy/7fqFjpGyDaxKzELVNYhammRPVoCupReT54ne1Ab/LEdwje9YkP786xNbVa1E0G8mtwVQ9S0nZqbiqGvfWRcGVPaeYnUhi9pKcxCy0FYt6kdw6QBtvFBLcS1QlttRRW4r2Nk1DJ4oiuTWIKgZyp27wz3MDjSG4NksuRv1jd/uvY7JsIVIQ1KFEtJMig7m3P7eco14+/ufy488fwfGHPh2hZi9RpoPJgIEg2tLZJJaQ84gtRtRWJCWp9rZ2YvvbP2uSIjcRnbw30ljRR9MjuZj1m7WB2kLkRXJrGTG6SpdJTaZKW8W8WTdNcrHrk/e9Cv2xUeazjZHuDrm22x5NiOpQWlIsYFpqMpTQcW/DxEhRDjMslKpTlqnkmkpsWcRcmLTKDlBLn9lb2bFEs1DklpAqv1hN7DFZ5MaaKt1WVTSX8jhKRcbFnlIY2GUktxqpYsLXSQ36TUxNDkh5E08ln9TyLPKe5PlxkSIlWUVnEk2aLCahtGRJ7Kld+MoVyfa/9FnYGy8jFEyM1CQUS09C/BTlKE1qk8sidcSWJbaYKck8aIybKIMiN1GIPL/0i7b7KA1X/D1o8vyR49CckiI2uqJaSKwbQdkeblUJbhYlV+a8q0pHQv5rqIltw6KbSG4dIvaNIyQdVYXgYLaiuDLnGltsdaUkxxEyDKAJS0uJZiC5VUhdX7xpDfsxxidVKbguS67s+VUttmnXjmYmEXUjuTWM0EGqVbVRhP5yr0pw0D3JxTifKlORsRl3LYd0JsnbU3LpU7tzlRdxMbPLzex+M3vAzH59QpmzzOyefpm/K3M89ZaMzNKndrN3ZfNWI542z+Qzuw7h8BWTv/ihPSfzzD1ZtBflMKNCSNm7MiaxxZxCbCmiNrW3zS5mdhrwX4DXAXuAm83sr9z920NljgA+B5zt7o+Z2VFljim5tYA8wwEWP7uY/cvra3fIKzigtOQGNFV2qaLMtvWIFDPNK4E73P15gH5U9k7gfw2VeS9wo7s/BuDu28scUHITQLzoDfKvHhAjihvHsFSqFl3KtGlRqdUZtZVFc0pWS29VgNwD5Feb2ZahxxvdfWP/7/uBPzCzVcBuYD2wZWT7k4GlZnYrsAK40t2vyV35PpJbzSzb5exZYYW2XfKDRew77EBw+SJL4AyTWnAQL4obJXVUV0UbYJlILZbYYqPxbZ1ip7vPjXvB3R8ys08CXweeBb4J7BsptgR4DfBm4BDgdjO7w92/VaQyklsHKZqazIre8lJk/bdUUdwobeqQUjb9GFNsRXvX5mlvizUziSZNbhbu/nng8wBm9j+B+ZEi8/QE+RzwnJndBpwOFJKbfjYlpsgXbFwKJtYXPittlHXzyvvLvkjPvEd3r1J7Up+qxBZC1rVRd/d/jXFrNoMOIma2FvhF4LqRIl8F3mBmS8zsUOAM4KGix1Pk1lFSdizJO+9k0RW8U6cqm0zVUkuZjizbS3Lcjz1NmNxKvtxvc9sLXOLuT5nZxQDuvqGfurwZuBc4AFzl7vcXPZjkVjEve2Y/Lxwer0t03nY3yG57C0lPViU4mC3JxYhYU4gtRdSm9rbZwt3fMOa5DSOPPwV8KsbxdHV1mDK/mEPaVqpIUQ7T9XRlW8U2DY1tE3UhuTWAcSmWPF2fi/wCjtU+UrXg4CXJdUF0Mc+lrlRkzKhNy9yIWCgt2SKqXtsttPdklSnKUUal0PTUZSohpxJb0zuRCDEJya3jTOtYEjLuLaXggGiSG9A02aWOLotEwrHEloVSkqJOJLeOUKRjCdQrOEgnuQFVy66qVGnR9G5MsRWN2tSRRFSB5FYD43pMjpupZNku2LPi4G2LpCazhgWUnblkmMHNs6jkIJ3ooN3zMZZpr4zd1T9LbEWitnHtbaHDAMqOcTvwxJOlthfNQ3JrIZMENy16Kyu4vLOXFIniBlQlujZQtgNOEamlbGdrStRmT2myyq4juXWMounJEIoIDvJHccOkTls2lRi9SvOKLcbCtRAvaouJ1nKbPSS3miiTmoT60pODG2BVUdyAWYjmYk2VlSJaG1AmHZk3aiuTktS8kuVYtN9bP52Z5NZBykRvoe1vdURxA0Yl0FbZpVgRO2W0lqrbv8a2iRRIbhWw9Jm97D18afz9ToneJgkuZM7JPIKD/FEcxJHcgHGSaJrwUohsmNQpyBCxqeu/aBKSW8PIk5osSkzBQbGlcoZvxjFFN2CaTFKJL7XARinaA7IOsaVKSQoxCcmt5RSJ3qAZghuQWnSjVC2hWJTtzl+kw0hqsSklKVIhuSVg6VO72bsy+0YSe4WAvKQQHORLU45SteiaTozxaUV7QdYRseWl7Z0eRDoktwaSNzVZNHoLJe8g7xiSg9kWXaxB16mitRhMitpipyQ1DGA2kdxmgLLpSXjphleH5GA2RBdzFpGU0dqAuqO2SWgYgADJLQr21C585fQeH6l6TL64/xIrBuRZtbvIVF0xJQcLJdBW2aVY/brMQOwqxaa2NpGaILmZ2dnAlcBiekt/f2Lk9fcBH+0/fBb4VXf/ZsyKinK9JrPSk4ObVWgUV2QuytiSGzBJEk2RXgqJjaMJ0VoI08SWZx1DIaaRKTczWwx8FngrMA/cZWY3ufuDQ8W+A/xHd3/KzM4BNgJnpKhw14jZqSTGem8p05QDUklulGlSiS2+qgQ2Stkps1KILXY6clJ7mzqTtAsz+w3gVwAH7gMudPcfDr1u9IKo9cDzwAXu/o2ixwuJ3F4HbHX3R/oVuB44D3hRbu7+j0Pl7wCOK1qhNnDgiSdZtOrIqWVCe0xOYlynkt7z6aK3AXnTlFBOcgNSy26YumQUgxhzQObtNBJLbEpHziZmdizwa8Cp7r7bzL4EvAe4eqjYOcBJ/X9nAH9CiSApRG7HAo8PPZ7POOCHgb8e94KZXQRcBLB27drAKoo8ZEVvKQQH5SQ3YPimXaXo2kCsSY2L9ISsSmxKSXaeJcAhZrYXOBTYNvL6ecA17u7AHWZ2hJn9mLt/r+jBslgYPvTCyoUFzd5IT24/O+51d99IL2XJ3NzczE03EKtTSVb0FlNwENYONyCG5ECiiyWzASmlBhKbmI67f9fM/hB4DNgNfM3dvzZSbFwgdSyQTG7zwPFDj49joXExs1cBVwHnuHva5Y47xqR2t0mpyRgMbkYpojg4+GYaU3QDuia82DIbUHTMWsz2tTJiy9vepmEAcbD9XuS9XG1mW4Yeb+wHNJjZSnqR2QnA08Cfm9n73f2Lw4cds8/CQVCI3O4CTjKzE4Dv0suTvne4gJmtBW4EPuDu3ypamVkmb8eSstHbgJRR3IBY0dwwoTKoW4KppDWNMoOwY3caUcQ2U+x097kJr70F+I677wAwsxuBnwGG5RYUSIWSKTd332dmlwK30BsKsMndHzCzi/uvbwB+H1gFfK7X4YV9U05ypsmbmpwWvcUUHIRHcVBOchBXdNOoQy51UHZWkZgpyAESmxjiMeBMMzuUXlryzcCWkTI3AZf2Oy2eATxTtL0NAse5uftmYPPIcxuG/v4Vel08RQmKpCdjCQ6KSW5AmbQlVCe7rhBriqwUUoM4YlNKsju4+51mdgPwDWAf8M/AxpEgaTO9YQBb6Q0FuLDMMTVDyQyQd/xbkfkoy0R0MP5mLeH1iD3XY5GB2E0Rm2gv7n4FcMXI08NBkgOXxDqe5FYD01KTKaI3KCY4CIvihhm+cRYV3YBZFF7KSYtTSw2UihTNQXJLSNmB3HkIFRxUIzmIK7oBoTf/pkmwqpn2Ryk6XVaRWUZiiS1m1KYVAWYXya0mYkdvuY5dUHJlls4pm7bMS10yqZsYcz/GjtYgntjU3iZCkdw6xOAGEjo9V17JlYniBqSI5madGEKDNFKD9GLLiz2l3OgsILlFImTZmzyUid7yzj9ZVHIQT3QDJLzpxBLZMKmkBvWkIgcoJTnbSG6JmdbuVnQ6rlDBQVrJQTzRDZh085416aWQ2DAp2tRGiSm2aVGbUpJiHJJbg5k2a0lo+1uRVQSKSA4W3jBjyG5A1s2+bfJLLa9xFF2KpshM/lWJTYhJSG4RKZKaTL1CNxRfJqfs+nCxo7pp5JVFbBnWIassyq6rVnR5mipTkXmjNrW3zQ6SWwWUGRKQFb0BwRFcr2y+4xeN4kapUnQhNFFGZYm5SGhTxFY0alN7m5DcOkCeIQJlojgoLzlIm76cJWKveA3lFhOtWmyx2toOPPFklP10Cdvvrf+BILkVJGQ17lCyUpMhKwZUEcXBwTe/GKKD6TfpWRdfCoGNUnZ17DyzjtQpNqUkZwvJLTKT2t2yUpMxBAf5o7gBZUUH8WQ3TMjNva0CrEJco5QV2YC802hVOVdk2yMOEQfJrQQxozeoR3AvbfPS30VEB9XIbhx1SKItxJLZgCJzQ+YVW1XpSFBKsstIbgkoGr1BXMFBWJpy4bYv/V1UdDD+xlqV8GaR2CIbUHSy4zxSS93dXynJ2UNyq5gqBQflJNfb/qW/y4huQOgNWBI8mFTimkSZ2ftTSi0ralNKUgyQ3EpSJDUZS3BAbslBHNH19lNoN0Fk3cy7Jr+q5TUg5hI0dUdqecWmlGS3kdwSkTWgO4bgIF8UN2DcTahs+rK3j9y7KExdMmgrqdZRi92eNo2ibW1KSc4mkluN1Cm4UVJEdr19Fa2RKEoVC4JWKTWIn45U1NZ9JLeEhEzHFVNwEJ6mnEYM0b20r/CyEmE2da5kXaQ7f2qpZaGobXaR3BpALMHBwTeTpoku+1hh5boqwTrFNUqMcWlVpSDViaQ9mNkRwFXAaYADH3L328eUey1wB/Bud7+hyLEkt4YQU3ADUooO0stucj1qOWznSDWwuo52NdEargRudvd3mdky4NDRAma2GPgkcEuZA0luicmzUkAKwQ0YveGkkB3UJzyxkCpnBYF6pDYtalNKslmY2WHAzwEXALj7HmDPmKKXAV8GXlvmeJJbw0gpuGHG3YhSCW8UCTA+VYtsQJ1RmtKRreNEYAfwBTM7HbgbuNzdnxsUMLNjgXcCb0Jyaz5513kLFRwQdS24VMIbpeiNeFakWJeosog1Nq0KqSlqK8m+/UXew9VmtmXo8UZ33zj0eAnwauAyd7/TzK4EPgZ8fKjMp4GPuvt+s3Lfd8mtoYSuATd6o4i98GmKdGZRYt/0Y8myqTIqQ+xB1rHa0hStNZqd7j435fV5YN7d7+w/voGe3IaZA67vi201sN7M9rn7V/JWRnJrMEUWOR2+iaRY4btJsitLF6WUh9TzOcbsHCKptR93/zcze9zMTnH3h4E3Aw+OlDlh8LeZXQ38ZRGxgeRWGXlTkwPKrOKdOqqDsBtkmwXYVlKLaxp1S00pyUZzGXBtv6fkI8CFZnYxgLtviHkgyS0CsZe+GaWM4A7aTwWyG0eRG62EeDB1ymoaqbruK1LrJu5+D73U4zBjpebuF5Q5luQWiRDBDX5RFo3ghkkhu4P2X5H4JtHUm/ksUtXYs7JCC43YNPXWbCC5RSQ0giuaohxm3I0ghvBe3H/ADa1uAYpy1D1gOmZ0JrGJUSS3yOQRHBSL4iaRWngLjpfj5igRpqduWYUQO92Yp31NYpstJLcE5GmDixHFTSNFOrNQPWq68VYt1TYIpkpStZ3l7TQisc0eklsi8goO4kZxk5h2s6lLfCmRbMrRxI4dEpsIQXJLSN5elFVKbhxFbmRdFGLXaKKgilCki7/ENrtIbokpMkwgxjidqgRZ9Y1zlmXaFUmNI8XYNIlttpHcKiD1OLhxjN4s6ooGY9PlG/wskXqgtcQmJLeKqENww0y6mXRFeiI+bZ3pQ2ITILlVSt2CG0esG5gk2QzaKqRYSGyR2L+v9e+l5FYxgwumaZIry6zfVEW9tP1GLOKzqO4KzCr6MgoRB32XxDiCIjczOxu4ElgMXOXunxh53fqvrweeBy5w929ErmvnqOtL2bWoUbQbyUmkIFNuZrYY+CzwVnqLzd1lZje5+/A6POcAJ/X/nQH8Sf9/0UB0MxFCdJ2QtOTrgK3u/oi77wGuB84bKXMecI33uAM4wsx+LHJdhRBCiCBC5HYs8PjQ4/n+c3nLYGYXmdkWM9uyY8eOvHUVQgghgghpc7Mxz3mBMrj7RmAjwNzc3ILXm8wtz11TdxWEEEIEEhK5zQPHDz0+DthWoIwQQghRCSFyuws4ycxOMLNlwHuAm0bK3AR80HqcCTzj7t+LXFchhBAiiEy5ufs+4FLgFuAh4Evu/oCZXWxmF/eLbQYeAbYCfwp8JFF9hRBCtBQzW2xm/2xmfznmtcPN7P+a2TfN7AEzu7DMsYLGubn7ZnoCG35uw9DfDlxSpiJCCCE6z+X0gqTDxrx2CfCgu/+Cma0BHjaza/u99HOjGUqEEEIkx8yOA34euGpCEQdW9CcFWQ48CewrejzNLSmEEKIKPg38NjBplvXP0Ou/sa1f5t3ufqDowayXUaweM9sBPFrLwYuzGthZdyUioPNoFjqPZtHG83iFu6+JtTMzu5ne+5CHHwF+OPR4Y3/4F2Z2LrDe3T9iZmcBv+Xu544c813A64H/Cvw48HXgdHf/QZFzqC1yi/lBVIWZbXH3ubrrURadR7PQeTSLrpxHGdz97Mi7fD3wDjNbT0+Ch5nZF939/UNlLgQ+0e/DsdXMvgP8JPBPRQ6oNjchhBBJcfffcffj3H0dveFkfzMiNoDHgDcDmNnRwCn0euEXQm1uQgghamEwnKzf+/5/AFeb2X30Zr36qLsXTg9LbvnYWHcFIqHzaBY6j2bRlfNoJO5+K3Br/+/hIWXbgLfFOk5tHUqEEEKIVKjNTQghROeQ3MZgZmeb2cNmttXMPjbmdTOzP+6/fq+ZvbqOemYRcB7v69f/XjP7RzM7vY56ZpF1HkPlXmtm+/tdihtHyHmY2Vlmdk9/+qG/q7qOIQRcV1GnUUqBmW0ys+1mdv+E11vxHRdTcHf9G/oHLAb+H3AisAz4JnDqSJn1wF/Ta/Q8E7iz7noXPI+fAVb2/z6nrecxVO5v6E0T9666613w8zgCeBBY2398VN31Lngevwt8sv/3GnozTSyru+4jdfw54NXA/RNeb/x3XP+m/1PktpCurDyeeR7u/o/u/lT/4R30lipqGiGfB8BlwJeB7VVWLgch5/Fe4EZ3fwzA3Zt4LiHnEXUapRS4+2306jWJNnzHxRQkt4VEW3m8ZvLW8cP0fqk2jczzMLNjgXcCG2guIZ/HycBKM7vVzO42sw9WVrtwQs7jM8Ar6U2jdB9wuZeYRqkm2vAdF1PQUICFRFt5vGaC62hmb6Qnt59NWqNihJzHp+mNidnfCxYaSch5LAFeQ28g6yHA7WZ2h7t/K3XlchByHm8H7gHeRH8aJTP7ey84jVJNtOE7LqYguS2kKyuPB9XRzF5Fb5buc9z9iYrqloeQ85gDru+LbTWw3sz2uftXKqlhGKHX1U53fw54zsxuA04HmiS3kPOIOo1STbThOy6moLTkQrqy8njmeZjZWuBG4AMNiw6GyTwPdz/B3dd5b2qfG4CPNExsEHZdfRV4g5ktMbNDgTPorX3VJELOI+o0SjXRhu+4mIIitxHcfZ+ZDVYeXwxs8v7K4/3XN9Drkbee3srjz9P7pdooAs/j94FVwOf6Uc8+b9iEsYHn0XhCzsPdH7LebOz3AgeAq9x9bFf1ugj8PKJOo5QCM7sOOAtYbWbzwBXAUmjPd1xMRzOUCCGE6BxKSwohhOgckpsQQojOIbkJIYToHJKbEEKIziG5CSGE6BySmxBCiM4huQkhhOgckpsQQojO8f8BeCIMYOJQe7IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from FEToy.mesh.lagrange2dmesh import mesh2d\n",
    "from FEToy.fe.gaussrule import gausspoint2d\n",
    "from FEToy.fe.shapefun import shape2d\n",
    "from FEToy.postprocess.PlotResult import Plot2D\n",
    "\n",
    "mymesh=mesh2d(nx=15,ny=15,meshtype='quad4')\n",
    "mymesh.createmesh()\n",
    "mymesh.plotmesh()\n",
    "\n",
    "gpoints=gausspoint2d(ngp=2)\n",
    "gpoints.creategausspoint()\n",
    "\n",
    "shp=shape2d(meshtype='quad4')\n",
    "shp.update()\n",
    "\n",
    "nDofs=mymesh.nodes*1\n",
    "K=np.zeros((nDofs,nDofs))\n",
    "F=np.zeros(nDofs)\n",
    "\n",
    "# element loop\n",
    "for e in range(mymesh.elements):\n",
    "    elconn=mymesh.elementconn[e,:]\n",
    "    nodes=mymesh.nodecoords[elconn,:]\n",
    "    # do integration\n",
    "    for gp in range(gpoints.ngp2):\n",
    "        xi =gpoints.gpcoords[gp,1]\n",
    "        eta=gpoints.gpcoords[gp,2]\n",
    "        w  =gpoints.gpcoords[gp,0]\n",
    "        \n",
    "        shp_val,shp_grad,j=shp.calc(xi,eta,nodes[:,0],nodes[:,1]) # xi,eta,x,y\n",
    "        JxW=j*w\n",
    "        for i in range(mymesh.nodesperelement):\n",
    "            iInd=elconn[i]\n",
    "            F[iInd]+=f*shp_val[i]*JxW\n",
    "            for j in range(mymesh.nodesperelement):\n",
    "                jInd=elconn[j]\n",
    "                K[iInd,jInd]+=(shp_grad[j,0]*shp_grad[i,0]+shp_grad[j,1]*shp_grad[i,1])*JxW\n",
    "####\n",
    "# Appyl boundary conditions\n",
    "iInd=mymesh.bcnodeids['left']\n",
    "K[iInd,iInd]+=1.0e16\n",
    "F[iInd]+=y2*1.0e16\n",
    "\n",
    "iInd=mymesh.bcnodeids['right']\n",
    "K[iInd,iInd]+=1.0e16\n",
    "F[iInd]+=y2*1.0e16\n",
    "\n",
    "iInd=mymesh.bcnodeids['bottom']\n",
    "K[iInd,iInd]+=1.0e16\n",
    "F[iInd]+=y2*1.0e16\n",
    "\n",
    "iInd=mymesh.bcnodeids['top']\n",
    "K[iInd,iInd]+=1.0e16\n",
    "F[iInd]+=y2*1.0e16\n",
    "\n",
    "### solve the equaiton\n",
    "phi=np.linalg.solve(K,F)\n",
    "\n",
    "### plot result\n",
    "Plot2D.Contour2D(mymesh,phi,savefig=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "265a0327-68de-4fc8-9e1e-63c7911982de",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "write result to poisson-000000.vtu\n"
     ]
    }
   ],
   "source": [
    "from FEToy.postprocess.Result import ResultIO\n",
    "\n",
    "# result=ResultIO('result')\n",
    "# result.save2csv(mesh=mymesh,solution=phi,varnamelist=['phi'],step=0)\n",
    "\n",
    "result=ResultIO('poisson')\n",
    "result.save2vtu(mesh=mymesh,solution=phi,varnamelist=['phi'],step=0)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}