5-Diffusion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "caf83a53-b6e5-4d2c-9d2f-a9e92f688daf",
   "metadata": {},
   "source": [
    "# Your first transient solver- diffusion equation\n",
    "\n",
    "- [x] time stepping\n",
    "- [x] time integration\n",
    "- [x] initial condition\n",
    "- [x] flux/surface integration\n",
    "- [x] 1d/2d diffusion\n",
    "\n",
    "Author: Yang Bai @ M3 Group\n",
    "\n",
    "Date  : 2022.04.15\n",
    "\n",
    "QQ group: 628204857"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb41194c-1b90-42d9-80a0-d0a78c47aa22",
   "metadata": {},
   "source": [
    "# Governing equation\n",
    "\n",
    "$\\frac{\\partial c}{\\partial t}=-\\nabla\\cdot\\mathbf{j}$, with $\\mathbf{j}=-D\\nabla c$\n",
    "\n",
    "$\\dot{c}=\\nabla\\cdot(D\\nabla c)+\\dot{r}(c,\\nabla c,...)=D\\nabla\\cdot(\\nabla c)=D\\Delta c=f(c)$\n",
    "\n",
    "$c$: concentration of species or normalized concentration(molar fraction)\n",
    "\n",
    "$D$: diffusion coefficient"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af56e62b-fc7a-4af2-8c62-138acb573d98",
   "metadata": {},
   "source": [
    "# Time integration\n",
    "\n",
    "$c^{n}$: solution of current time step (unknown)\n",
    "\n",
    "$c^{n-1}$: solution of previous time step (known)\n",
    "\n",
    "$\\dot{c}=f(c)$ with $\\dot{c}=\\frac{c^{n}-c^{n-1}}{\\Delta t}$\n",
    "\n",
    "===> $\\dot{c}^{n}=f(c^{n})(\\mathrm{implicit, backward-euler~method~(BE)})$\n",
    "\n",
    "===> $\\dot{c}^{n}=f(c^{n-1})(\\mathrm{explicit, forward-euler~method, CFL~condition})$\n",
    "\n",
    "## weighted integration\n",
    "$\\int_{\\Omega}\\dot{c}\\delta c dV=\\int_{\\partial\\Omega}D\\nabla c\\cdot\\vec{n}\\delta cdS-\\int_{\\Omega}D\\nabla c\\nabla\\delta c dV$\n",
    "\n",
    "## time discretization\n",
    "### use BE\n",
    "$\\int_{\\Omega}\\frac{c^{n}-c^{n-1}}{\\Delta t}\\delta c dV=\\int_{\\partial\\Omega}D\\nabla c^{n}\\cdot\\vec{n}\\delta cdS-\\int_{\\Omega}D\\nabla c^{n}\\nabla\\delta c dV$, $-D\\nabla c^{n}\\cdot\\vec{n}=j_{0}$\n",
    "\n",
    "### use FE\n",
    "$\\int_{\\Omega}\\frac{c^{n}-c^{n-1}}{\\Delta t}\\delta c dV=\\int_{\\partial\\Omega}D\\nabla c^{n-1}\\cdot\\vec{n}\\delta cdS-\\int_{\\Omega}D\\nabla c^{n-1}\\nabla\\delta c dV$, $-D\\nabla c^{n-1}\\cdot\\vec{n}=j_{0}$\n",
    "\n",
    "## K matrix and right hand side F (BE)\n",
    "$K=\\int_{\\Omega}\\frac{1}{\\Delta t}N^{J}N^{I}dV+\\int_{\\Omega}D\\nabla N^{J}\\nabla N^{I}dV$\n",
    "\n",
    "$F=\\int_{\\Omega}\\frac{c^{n-1}}{\\Delta t}N^{I}dV-\\int_{\\partial\\Omega}j_{0}N^{I}dS$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1414a01-2b70-44fb-8562-ee02451bd31e",
   "metadata": {},
   "source": [
    "# Diffusion in 1d case"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6e409032-9d4b-4f57-a6fc-5aeb617fbac7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEKCAYAAAA1qaOTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfBklEQVR4nO3dfXhU9Z338fc3BBCKPAgkRh4SdBEhgIABY3chWpooaqu9a1vxiUq9EbzFWguI1VXr9lqt1VXbWl1rrbZ3r7q2PtGKN2tRwV4VDWwFRIogooaAoICwQcjDfO8/ZsgOkzOTQX55Kp/Xdc1Fzjnf+Z3fb35zzmfmzCSYuyMiInK4ctq6AyIi8vdBgSIiIkEoUEREJAgFioiIBKFAERGRIHLbugNtqV+/fl5UVNTW3RAR6VBWrFjxkbv3T11/RAdKUVERy5cvb+tuiIh0KGb2XtR6XfISEZEgFCgiIhKEAkVERIJQoIiISBAKFBERCUKBIiIiQShQREQkCAWKiIgEoUAREZEgFCgiIhKEAkVERIJQoIiISBAKFBERCUKBIiIiQShQREQkCAWKiIgEoUAREZEgFCgiIhKEAkVERIJQoIiISBAKFBERCUKBIiIiQShQREQkCAWKiIgEoUAREZEg2lWgmNlZZrbOzDaY2fyI7WZmP05sX2Vm41K2dzKzv5rZH1uv1yIiAu0oUMysE3A/MAUYAUw1sxEpZVOAoYnbDOCBlO3fBta2cFdFRCRCuwkUYAKwwd03unst8DhwXkrNecCvPG4Z0NvMCgDMbCBwDvBwa3ZaRETi2lOgDAA+SFquSqzLtuZeYB4Qy7QTM5thZsvNbPn27dsPq8MiIvI/2lOgWMQ6z6bGzM4Ftrn7iuZ24u4PuXuJu5f079//s/RTREQitKdAqQIGJS0PBKqzrPlH4Mtmton4pbIvmNn/bbmuiohIqvYUKJXAUDMbYmZdgAuBBSk1C4DLEt/2KgU+cfct7n6Duw9096LE/V5090tatfciIke43LbuwAHuXm9mVwOLgE7AI+6+xsxmJrY/CCwEzgY2AHuBy9uqvyIicjBzT/2Y4shRUlLiy5cvb+tuiIh0KGa2wt1LUte3p0teIiLSgSlQREQkCAWKiIgEoUAREZEgFCgiIhKEAkVERIJQoIiISBAKFBERCUKBIiIiQShQREQkCAWKiIgEoUAREZEgFCgiIhKEAkVERIJQoIiISBAKFBERCUKBIiIiQShQREQkCAWKiIgEoUAREZEgFCgiIhKEAkVERIJQoIiISBAKFBERCUKBIiIiQShQREQkCAWKiIgEoUAREZEgFCgiIhKEAkVERIJoV4FiZmeZ2Toz22Bm8yO2m5n9OLF9lZmNS6wfZGYvmdlaM1tjZt9u/d6LiBzZ2k2gmFkn4H5gCjACmGpmI1LKpgBDE7cZwAOJ9fXAd919OFAK/J+I+4qISAtqN4ECTAA2uPtGd68FHgfOS6k5D/iVxy0DeptZgbtvcff/AnD3PcBaYEBrdl5E5EjXngJlAPBB0nIVTUOh2RozKwLGAq+F76KIiKTTngLFItb5odSYWQ/gSeBad98duROzGWa23MyWb9++/TN3VkREDtaeAqUKGJS0PBCozrbGzDoTD5PfuPtT6Xbi7g+5e4m7l/Tv3z9Ix0VEpH0FSiUw1MyGmFkX4EJgQUrNAuCyxLe9SoFP3H2LmRnwC2Ctu/9b63ZbREQActu6Awe4e72ZXQ0sAjoBj7j7GjObmdj+ILAQOBvYAOwFLk/c/R+BS4HVZvZGYt333H1hKw5BROSIZu6pH1McOUpKSnz58uVt3Q0RkQ7FzFa4e0nq+vZ0yUtERDowBYqIiAShQBERkSAUKCIiEoQCRUREglCgiIhIEAoUEREJQoEiIiJBKFBERCQIBYqIiAShQBERkSAUKCIiEoQCRUREglCgiIhIEAoUEREJQoEiIiJBKFBERCQIBYqIiAShQBERkSAUKCIiEoQCRUREglCgiIhIEAoUEREJQoEiIiJBKFBERCQIBYqIiAShQBERkSAUKIdh+vTp5OXlMXLkyIx1H3zwAWeccQbDhw+nuLiY++67L7Ju3759TJgwgZNPPpni4mJuueWWjO02NDQwduxYzj333Ix1RUVFjBo1ijFjxlBSUpK2bteuXVxwwQWcdNJJDB8+nFdffbVJzbp16xgzZkzjrWfPntx7772R7d1zzz0UFxczcuRIpk6dyr59+yLr7rvvPkaOHElxcXHatiD68d6xYwfl5eUMHTqU8vJydu7cGVn3u9/9juLiYnJycli+fHna9ubOnctJJ53E6NGj+cpXvsKuXbsi6/75n/+Z0aNHM2bMGCoqKqiurs74fLjrrrswMz766KPIultvvZUBAwY0Pq4LFy5M295PfvIThg0bRnFxMfPmzUs7lm984xuN7RUVFTFmzJjIujfeeIPS0tLG58frr78eWbdy5UpOO+00Ro0axZe+9CV2796d9rmdOi+rV6+OrEudl3Ttpc7Lm2++GVmXOi+VlZUZj70D87Jy5crIuqh5yXQ8J8/NzJkzI+tS52XEiBGRdanz8oc//CGyLnVetm3bFnkeSZ2TLVu2RNZFHStZc/cj9nbKKaf44ViyZImvWLHCi4uLM9ZVV1f7ihUr3N199+7dPnToUF+zZk2Tulgs5nv27HF399raWp8wYYK/+uqradu9++67ferUqX7OOedk3H9hYaFv3769ueH4ZZdd5j//+c/d3X3//v2+c+fOjPX19fWen5/vmzZtarKtqqrKi4qKfO/eve7u/rWvfc1/+ctfNqlbvXq1FxcXe01NjdfV1fnkyZP97bffjtxf1OM9d+5cv/32293d/fbbb/d58+ZF1r311lv+t7/9zcvKyryysjJte4sWLfK6ujp3d583b17a9j755JPGn++77z6/8sor0z4f3n//fa+oqPDBgwf79u3bI+tuueUW/9GPftTseF988UWfPHmy79u3z93dP/zww7S1ya677jr//ve/H1lXXl7uCxcudHf35557zsvKyiLrSkpK/OWXX3Z391/84hd+0003pX1up87LVVddFVmXOi/p2kudl3Ttpc7LJZdckvbYS56X1atXR9ZFzUu6PqbOzapVq5o97q+77jqfM2dOZF3qvJSWlkbWpc7LjTfeGHkeSZ2TuXPnRtZFHSupgOUecU7VO5TDMGnSJI455phm6woKChg3bhwARx99NMOHD2fz5s1N6syMHj16AFBXV0ddXR1mFtlmVVUVzz33HFdcccVhjOB/7N69m6VLl/Ktb30LgC5dutC7d++M91m8eDEnnHAChYWFkdvr6+v59NNPqa+vZ+/evRx33HFNatauXUtpaSndu3cnNzeXsrIynn766cj2oh7vZ599lmnTpgEwbdo0nnnmmci64cOHM2zYsGbbq6ioIDc3F4DS0lKqqqoi63r27Nn4c01NDWaW9vnwne98hzvvvLNxLrN93kTVPfDAA8yfP5+uXbsCkJeX12yb7s4TTzzB1KlTI+vMjN27dwPwySefcNxxx0XWrVu3jkmTJgFQXl7Ok08+mfa5nTovf/rTnyLrUuclXXup87Jr167IutR56dGjR9pjL3lejj322KyO0Ux9TJ2bUaNGZWzzwLzMmDEjsi51XoYMGRJZlzovTz31VOR5JHVOnn322ci6qGMlWwqUVrZp0yb++te/cuqpp0Zub2hoYMyYMeTl5VFeXp627tprr+XOO+8kJ6f5KTQzKioqOOWUU3jooYciazZu3Ej//v25/PLLGTt2LFdccQU1NTUZ23388ceZOnVq5LYBAwYwZ84cBg8eTEFBAb169aKioqJJ3ciRI1m6dCkff/wxe/fubbykkK0PP/yQgoICIH6gb9u2Lev7NueRRx5hypQpabffeOONDBo0iN/85jfcdtttkTULFixgwIABnHzyyc3u76c//SmjR49m+vTp7Ny5M7Lm7bff5pVXXuHUU0+lrKyMysrKZtt95ZVXyM/PZ+jQoZHb7733XubOncugQYOYM2cOt99+e2TdyJEjWbBgARC/LJI6T8nP7Uzz0twx0Fxd6ryk1qWbl+S6TPOS2l6meUmuzTQ3UWOJmpfkukzzklwXNS9R55GoOcn2fJO1qLctB25AaabtoW/AWcA6YAMwP2K7AT9ObF8FjMv2vlG3z3zJq7rafdIk9y1b/N13301/ySupzt19z549Pm7cOH/yyScz1rm779y5008//XRfvXp1k9o/FBf7rGnT3N39pZdeir7kldTm5s2b3T1+eWT06NG+ZMmSJnWVzz/vnTp18mXLlrm7+zXXXOM33XRT2j7u37/f+/bt61u3bo3c7461a/2MM87wbdu2eW1trZ933nn+61//OrK9hx9+2MeOHesTJ070K6+80q+99tq0+019vHv16nVQbe/c3IzzUlZW5pXPP9/s/P3gBz/w8886y2MTJzY7z/96ww1+8+DBTepqamp8woQJvmvXLnd3Lxw40Lefdlpke1u3bvX6+npvaGjw711zjV+enx9ZV1xc7LNnz/ZYLOavPfecF3Xt6rHq6sjH5oCZM2f6XTffnHbMs2fP9t///vfu7v4fDz7ok3v3jqxbu3atl5eX+7hx4/zW737Xj0k81u5Nn9uN85KYv96J5XTHQNlpp3nl2LHNHis/uP56P79v38Yxpz2mUuYluS7TvKS2l25eovbdODebN/trY8Z40aBBjZexo/o487LL/K7jj0/bXuO8VFf7fwwf7pMnToysa5yXUaP81sJCP6ZPn8Z9JJ9H0h0rqXWNc/IZLnk1d4KvBf4FyM1UF+IGdALeAY4HugArgREpNWcDzyeCpRR4Ldv7Rt0+c6DMmuWek+M+a1bmQEmqq62t9YqKCr/77rsz1iW79dZbm1y/9VmzfD74gM99zgsLCz0/P9+7devmF198cVZtNrkmnKjbMm2aFxYWNq5eunSpn3322Wnbe+aZZ7y8vDztWJ6oqPDp06c3rn7sscd8VnJf0vTvhhtu8Pvvvz9tXerjfeKJJ3p14uRSfdllfiJknJeysjKv/OpXM87fo48+6qWlpV5zxRVZzfOmSy7x4oj9rlq1yvv37++FhYVeWFjoncx8EPiWadMytvfuxRdHtufufuaZZ/pLL73U+NgcD77tm9+MfGzc3evq6jwvL88/uPTStGPp2bOnx2Ixd3ePzZzpRzfzGLq7r5s61ccn6qKe243zMmuWV5v5ib16ZTwGyo47zivNMh4rjz76qJfm53tNM3UHHJiX2hkzDqpLNy/vX3ppxvaS5yVq341zk3jOHt+zp2/evDmyzbq6Os/r1s0/yDCWxnmZNctjZn50587NnkfWmfn4vLyDVh84j6Q7VlLrGufkMwSKxbdFM7PJwC+AHcAl7v7W4bwbysTMTgNudfczE8s3ALj77Uk1/w687O6/TSyvA04Hipq7b5SSkhI/pG8xdOsGKd9U2gScA/QvK2tct+iVV+gaizUuOzAN6AWszlC3HegMdMvJ4Yuf/zyrVq1i8ODB9O3bt0ntAf9pxtQ+fRg1alRkmzVADOiS1GZhYSGVa9Y0aW8i8DMzZk+axKZNm2hoaGBDdXXkfr8OrBo2jGOPPTZyv68B04E/m3H+xImsW7eOo48+mjc3bmzS3jagV04OZePHs2rVKsaOHcuLy5ZF7nedGeO6dWP8+PEAvPPOO1y5eTM3uXMH8SfqnaSflzNjMe4CDnzXLbVux44d1K1ezVKgf9J+U+v27t3LKytW0DUW4yfAEuD3GfbbNRajCFgO9Iuo279/Py+//jpdYzHuSTx+j0fUVVdXc8H69fwr8DYwGXif+Cus1McG4LtLl3KXO0syjKWyspKnP/2UCncWA/OAFRF1tbW1vPTaa3SOxfgm8QPvcqKf2+OWLCEPmA/cAXwMfBhRd+CxOR24CzglTXvfXbqU6xPj6E/6Yyp1Xl4GPgccA/wwJ4czJ06MnJdK4LsRdVHz8ts0+75g6VK2uXNbYm6+kLhFjfmlWIzbiT9v0o1l65IlPJh4nBcDc4GREXW/XrqUQe7EgG8mPYap55Fdu3YddKxUAT8AegOfdu1KxfjxXH/99Y3fGj399NO56667Ir8ZamYr3L3JhoyBkrjj0cC9wFTgJnf/t4x3+IzM7ALgLHe/IrF8KXCqu1+dVPNH4A53/3NieTFwPfFAyXjfpDZmADMABg8efMp7772XfSe3bIE5c6h94gm61NfzdWBxTg473encuTNFRUUUFBRwzP79zNq4kYkffcRRsRiLzfiiOz27d6ch8cHskCFDGNqjx0F1lWZM7dSJ7V26UA/079+foqIigCZt7svJ4ZV+/bitXz/e2rq1MVBS69aacW6nTnyUaDMvL4/CwsLI9h7t1YsbamupA4466iiGDRtGfizWpO6Fvn356s6dTDj11MYPSqPa+99HHcUzsRixnBx69OjBsGHD6FdX16RufKdOvJebSywnhxNOOIE+ffpEtjelc2dejcWora9vfLz79evH+jffpHNNDSc2NPAkMDNiXjp37sw769fTUFtLH2AM8QMzte7999/HYjEKYjF61NdzGvGQSq3bsWMH+2tq6Fdby/CGBh4CvtPM8+GKbdtYAVwVUbdr1y727dlD39paRjQ08HPg2xF1+fn5bFq7lpydO+nT0MDdwOfTPDYFBQW8u2YN5+/fzw9rajgqFot8znbv3p1N69dzzP79HFNfz4PADyPqGhoa2Lp5M73r6riwoYEfAS+meW4P7tqV6lWr2FNXRyEwG/gKNKn7XG0t72zYwO5YjN7AYOKhm1r37vr1dKur47hYDAMGAX+MqNu6detB83IF8L+AIZ06saNrV+rNGDJkCH379j1oXu4HvhRRt23btibz8laaMZ/QvTt7Vq5k6/79dAUuIh4CUcf9f7/xBl/av5/Z7mnPD31iMba8/Tbd6uvpRjwsropoz/bsYUdVFd0aGvhqYrwX5eY2OY/U1dUddKzcBlxjRkOvXsTy8/n6RRdx88038/TTTzN79my2b99O7969GTNmDIsWLTroVJguUA7lktRXgXrgv4Hdybds22im/a8BDyctXwr8JKXmOeCfkpYXEw/kZu8bdftMl7xmzoxfOjjqqMhLNi1W15b7bu91HaGPemw05vZQd6i1aXA4Xxs2sxLi747WE3/hNDvlFkIV8RceBwwEqrOsyea+YXz4IcycCcuWxf/durV16tpy3+29riP0UY+Nxtwe6g619hA19xlKLnAL8ctKPyP+7al9wfbedF8HLgtvJn5Z8yJ3X5NUcw5wNfEP508FfuzuE7K5b5RD/gxFRETSXvLKbeZ+lcQ/p5ri7otbpGcJ7l5vZlcDi4h/a+sRd19jZjMT2x8EFhIPkw3AXuKfCaa9b0v2V0REDtbcO5RfA1e7+yet16XWo3coIiKH7jO9Q3H3S1uuSyIi8vdEf3pFRESCUKCIiEgQChQREQlCgSIiIkEoUEREJAgFioiIBKFAERGRIBQoIiIShAJFRESCUKCIiEgQChQREQlCgSIiIkEoUEREJAgFioiIBKFAERGRIBQoIiIShAJFRESCUKCIiEgQChQREQlCgSIiIkEoUEREJAgFioiIBKFAERGRIBQoIiIShAJFRESCUKCIiEgQChQREQlCgSIiIkEoUEREJIh2EShmdoyZvWBm6xP/9klTd5aZrTOzDWY2P2n9j8zsb2a2ysyeNrPerdZ5EREB2kmgAPOBxe4+FFicWD6ImXUC7gemACOAqWY2IrH5BWCku48G3gZuaJVei4hIo/YSKOcBjyV+fgw4P6JmArDB3Te6ey3weOJ+uPt/unt9om4ZMLBluysiIqnaS6Dku/sWgMS/eRE1A4APkparEutSTQeeD95DERHJKLe1dmRmfwKOjdh0Y7ZNRKzzlH3cCNQDv8nQjxnADIDBgwdnuWsREWlOqwWKu38x3TYz+9DMCtx9i5kVANsiyqqAQUnLA4HqpDamAecCk93dScPdHwIeAigpKUlbJyIih6a9XPJaAExL/DwNeDaiphIYamZDzKwLcGHifpjZWcD1wJfdfW8r9FdERFK0l0C5Ayg3s/VAeWIZMzvOzBYCJD50vxpYBKwFnnD3NYn7/xQ4GnjBzN4wswdbewAiIke6VrvklYm7fwxMjlhfDZydtLwQWBhR9w8t2kEREWlWe3mHIiIiHZwCRUREglCgiIhIEAoUEREJQoEiIiJBKFBERCQIBYqIiAShQBERkSAUKCIiEoQCRUREglCgiIhIEAoUEREJQoEiIiJBKFBERCQIBYqIiAShQBERkSAUKCIiEoQCRUREglCgiIhIEAoUEREJQoEiIiJBKFBERCQIBYqIiAShQBERkSAUKCIiEoQCRUREglCgiIhIEAoUEREJQoEiIiJBKFBERCQIBYqIiATRLgLFzI4xsxfMbH3i3z5p6s4ys3VmtsHM5kdsn2Nmbmb9Wr7XIiKSrF0ECjAfWOzuQ4HFieWDmFkn4H5gCjACmGpmI5K2DwLKgfdbpcciInKQ9hIo5wGPJX5+DDg/omYCsMHdN7p7LfB44n4H3APMA7wF+ykiImm0l0DJd/ctAIl/8yJqBgAfJC1XJdZhZl8GNrv7yuZ2ZGYzzGy5mS3fvn374fdcREQAyG2tHZnZn4BjIzbdmG0TEevczLon2qjIphF3fwh4CKCkpETvZkREAmm1QHH3L6bbZmYfmlmBu28xswJgW0RZFTAoaXkgUA2cAAwBVprZgfX/ZWYT3H1rsAGIiEhG7eWS1wJgWuLnacCzETWVwFAzG2JmXYALgQXuvtrd89y9yN2LiAfPOIWJiEjrai+BcgdQbmbriX9T6w4AMzvOzBYCuHs9cDWwCFgLPOHua9qovyIikqLVLnll4u4fA5Mj1lcDZyctLwQWNtNWUej+iYhI89rLOxQREengFCgiIhKEAkVERIJQoIiISBAKFBERCUKBIiIiQShQREQkCAWKiIgEoUAREZEgFCgiIhKEAkVERIJQoIiISBAKFBERCUKBIiIiQShQREQkCAWKiIgEoUAREZEgFCgiIhKEAkVERIJQoIiISBAKFBERCUKBIiIiQShQREQkCAWKiIgEYe7e1n1oM2a2HXjvM969H/BRwO50BBrzkUFjPjIczpgL3b1/6sojOlAOh5ktd/eStu5Ha9KYjwwa85GhJcasS14iIhKEAkVERIJQoHx2D7V1B9qAxnxk0JiPDMHHrM9QREQkCL1DERGRIBQoIiIShAKlGWZ2lpmtM7MNZjY/YruZ2Y8T21eZ2bi26GdIWYz54sRYV5nZX8zs5LboZ0jNjTmpbryZNZjZBa3Zv9CyGa+ZnW5mb5jZGjNb0tp9DC2L53UvM/uDma1MjPnytuhnSGb2iJltM7M302wPe/5yd93S3IBOwDvA8UAXYCUwIqXmbOB5wIBS4LW27ncrjPnzQJ/Ez1OOhDEn1b0ILAQuaOt+t/Ac9wbeAgYnlvPaut+tMObvAT9M/Nwf2AF0aeu+H+a4JwHjgDfTbA96/tI7lMwmABvcfaO71wKPA+el1JwH/MrjlgG9zaygtTsaULNjdve/uPvOxOIyYGAr9zG0bOYZYDbwJLCtNTvXArIZ70XAU+7+PoC7HwljduBoMzOgB/FAqW/dbobl7kuJjyOdoOcvBUpmA4APkparEusOtaYjOdTxfIv4K5yOrNkxm9kA4CvAg63Yr5aSzRyfCPQxs5fNbIWZXdZqvWsZ2Yz5p8BwoBpYDXzb3WOt0702E/T8lXvY3fn7ZhHrUr9nnU1NR5L1eMzsDOKB8k8t2qOWl82Y7wWud/eG+AvYDi2b8eYCpwCTgW7Aq2a2zN3fbunOtZBsxnwm8AbwBeAE4AUze8Xdd7dw39pS0POXAiWzKmBQ0vJA4q9eDrWmI8lqPGY2GngYmOLuH7dS31pKNmMuAR5PhEk/4Gwzq3f3Z1qlh2Fl+7z+yN1rgBozWwqcDHTUQMlmzJcDd3j8w4UNZvYucBLweut0sU0EPX/pkldmlcBQMxtiZl2AC4EFKTULgMsS35YoBT5x9y2t3dGAmh2zmQ0GngIu7cCvWJM1O2Z3H+LuRe5eBPweuKqDhglk97x+FphoZrlm1h04FVjbyv0MKZsxv0/8HRlmlg8MAza2ai9bX9Dzl96hZODu9WZ2NbCI+LdEHnH3NWY2M7H9QeLf+Dkb2ADsJf4qp8PKcsw3A32BnyVesdd7B/5LrVmO+e9GNuN197Vm9v+AVUAMeNjdI7962hFkOcf/AjxqZquJXwq63t079J+0N7PfAqcD/cysCrgF6Awtc/7Sn14REZEgdMlLRESCUKCIiEgQChQREQlCgSIiIkEoUEREJAgFioiIBKFAEWkHzCzHzJaaWeovkXZP/Mn1B9qqbyLZUqCItAOJP0L4TeALZjY9adMPif8C8py26JfIodAvNoq0I4nf3L4TGAX8A/Hf7D7d3f/cph0TyYICRaSdMbNFxP/CbxHwuLvPa9seiWRHgSLSzpjZEOL/u+A7wEh339/GXRLJij5DEWl/pgOfEv9T4se3cV9EsqZ3KCLtiJmNB/4CfBmYBeQDn3f3hjbtmEgW9A5FpJ0ws6OAXwGPuvvzwAziH8zrMxTpEPQORaSdMLN7gPOB0e6+J7HuQuAx4JSO/P+RyJFBgSLSDpjZJOBF4Ivu/nLKtieIf5ZS6u71bdA9kawoUEREJAh9hiIiIkEoUEREJAgFioiIBKFAERGRIBQoIiIShAJFRESCUKCIiEgQChQREQni/wMlRPyp3f5BiwAAAABJRU5ErkJggg==\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",
    "from FEToy.postprocess.Result import ResultIO # to csv file\n",
    "import numpy as np\n",
    "\n",
    "# define your mesh\n",
    "my1dmesh=mesh1d(nx=10,xmax=1.0,meshtype='edge4')\n",
    "my1dmesh.createmesh()\n",
    "my1dmesh.plotmesh(withnode=True,withnodeid=True)\n",
    "\n",
    "#define your shape fun and gauss points\n",
    "shp=shape1d(meshtype='edge4')\n",
    "shp.update()\n",
    "\n",
    "qpoints=gausspoint1d(ngp=3)\n",
    "qpoints.creategausspoint()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "653a4d02-cdff-4011-9e84-9705d1b208a3",
   "metadata": {},
   "source": [
    "$$K=\\int_{\\Omega}\\frac{1}{\\Delta t}N^{J}N^{I}dV+\\int_{\\Omega}D\\nabla N^{J}\\nabla N^{I}dV$$\n",
    "$$F=\\int_{\\Omega}\\frac{c^{n-1}}{\\Delta t}N^{I}dV-\\int_{\\partial\\Omega}j_{0}N^{I}dS$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a75b79a0-af5d-4890-a098-7bc39847eba2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "write result to diff1d-000000.csv\n",
      "step=    1 finished!\n",
      "step=    2 finished!\n",
      "step=    3 finished!\n",
      "step=    4 finished!\n",
      "step=    5 finished!\n",
      "step=    6 finished!\n",
      "step=    7 finished!\n",
      "step=    8 finished!\n",
      "step=    9 finished!\n",
      "step=   10 finished!\n",
      "step=   11 finished!\n",
      "step=   12 finished!\n",
      "step=   13 finished!\n",
      "step=   14 finished!\n",
      "step=   15 finished!\n",
      "step=   16 finished!\n",
      "step=   17 finished!\n",
      "step=   18 finished!\n",
      "step=   19 finished!\n",
      "step=   20 finished!\n"
     ]
    }
   ],
   "source": [
    "# define your parameters\n",
    "D=0.5        # diffusion coefficient\n",
    "dt=1.0e-3    # time step size\n",
    "totalstep=20 # the total time step\n",
    "j0=0.005      # flux\n",
    "\n",
    "nDofs=my1dmesh.nodes*1\n",
    "K=np.zeros((nDofs,nDofs))\n",
    "F=np.zeros(nDofs)\n",
    "cold=np.zeros(nDofs) # solution in previous step\n",
    "c=np.zeros(nDofs)    # solution in current step\n",
    "\n",
    "# apply initial condition\n",
    "cold[:]=0.001\n",
    "\n",
    "# define a csv output\n",
    "result=ResultIO('diff1d')\n",
    "result.save2csv(mesh=my1dmesh,solution=cold,varnamelist=['c'],step=0) # to csv file\n",
    "\n",
    "for step in range(totalstep):\n",
    "    # init your k and rhs\n",
    "    K[:,:]=0.0;F[:]=0.0\n",
    "    for e in range(my1dmesh.elements):\n",
    "        elconn=my1dmesh.elementconn[e,:]  # local element connectivity\n",
    "        nodes=my1dmesh.nodecoords[elconn] # node coordinates of local element\n",
    "        # now we can do the integration\n",
    "        for gp in range(qpoints.ngp):\n",
    "            xi=qpoints.gpcoords[gp,1] # 1->xi, 2->eta\n",
    "            w =qpoints.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",
    "            # gradc=0.0\n",
    "            # for i in range(my1dmesh.nodesperelement):\n",
    "            #     iInd=elconn[i]\n",
    "            #     gradc+=shp_grad[i]*c[iInd]\n",
    "            # now we can calculate the F term\n",
    "            for i in range(my1dmesh.nodesperelement):\n",
    "                iInd=elconn[i]\n",
    "                # assemble to global F\n",
    "                F[iInd]+=(cold[iInd]/dt)*shp_val[i]*JxW\n",
    "                # now we can calculate the K matrix\n",
    "                for j in range(my1dmesh.nodesperelement):\n",
    "                    jInd=elconn[j]\n",
    "                    K[iInd,jInd]+=(1.0/dt)*shp_val[j]*shp_val[i]*JxW\\\n",
    "                                 +D*shp_grad[j]*shp_grad[i]*JxW\n",
    "    # apply the flux, in 1d case(bc elemt is 0-d), you don't need surface integration\n",
    "    iInd=my1dmesh.bcnodeids['right']\n",
    "    F[iInd]+=j0\n",
    "    \n",
    "    # solve the equation\n",
    "    c=np.linalg.solve(K,F) # current solution\n",
    "    cold[:]=c[:] # update solution\n",
    "    # result.save2csv(mesh=my1dmesh,solution=cold,varnamelist=['c'],step=step+1)\n",
    "    # result.save2vtu(mesh=my1dmesh,solution=cold,varnamelist=['c'],step=step+1)\n",
    "    print('step=%5d finished!'%(step+1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4cc5e11-7706-418f-aa49-15e0d9a6a4d9",
   "metadata": {},
   "source": [
    "# Diffusion in 2d case"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c3d44691-691e-42e7-aa62-c5bc4b2992ed",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "write result to diff2d-000000.vtu\n",
      "write result to diff2d-000001.vtu\n",
      "step=    1 finished!\n",
      "write result to diff2d-000002.vtu\n",
      "step=    2 finished!\n",
      "write result to diff2d-000003.vtu\n",
      "step=    3 finished!\n",
      "write result to diff2d-000004.vtu\n",
      "step=    4 finished!\n",
      "write result to diff2d-000005.vtu\n",
      "step=    5 finished!\n",
      "write result to diff2d-000006.vtu\n",
      "step=    6 finished!\n",
      "write result to diff2d-000007.vtu\n",
      "step=    7 finished!\n",
      "write result to diff2d-000008.vtu\n",
      "step=    8 finished!\n",
      "write result to diff2d-000009.vtu\n",
      "step=    9 finished!\n",
      "write result to diff2d-000010.vtu\n",
      "step=   10 finished!\n",
      "write result to diff2d-000011.vtu\n",
      "step=   11 finished!\n",
      "write result to diff2d-000012.vtu\n",
      "step=   12 finished!\n",
      "write result to diff2d-000013.vtu\n",
      "step=   13 finished!\n",
      "write result to diff2d-000014.vtu\n",
      "step=   14 finished!\n",
      "write result to diff2d-000015.vtu\n",
      "step=   15 finished!\n",
      "write result to diff2d-000016.vtu\n",
      "step=   16 finished!\n",
      "write result to diff2d-000017.vtu\n",
      "step=   17 finished!\n",
      "write result to diff2d-000018.vtu\n",
      "step=   18 finished!\n",
      "write result to diff2d-000019.vtu\n",
      "step=   19 finished!\n",
      "write result to diff2d-000020.vtu\n",
      "step=   20 finished!\n",
      "write result to diff2d-000021.vtu\n",
      "step=   21 finished!\n",
      "write result to diff2d-000022.vtu\n",
      "step=   22 finished!\n",
      "write result to diff2d-000023.vtu\n",
      "step=   23 finished!\n",
      "write result to diff2d-000024.vtu\n",
      "step=   24 finished!\n",
      "write result to diff2d-000025.vtu\n",
      "step=   25 finished!\n",
      "write result to diff2d-000026.vtu\n",
      "step=   26 finished!\n",
      "write result to diff2d-000027.vtu\n",
      "step=   27 finished!\n",
      "write result to diff2d-000028.vtu\n",
      "step=   28 finished!\n",
      "write result to diff2d-000029.vtu\n",
      "step=   29 finished!\n",
      "write result to diff2d-000030.vtu\n",
      "step=   30 finished!\n",
      "write result to diff2d-000031.vtu\n",
      "step=   31 finished!\n",
      "write result to diff2d-000032.vtu\n",
      "step=   32 finished!\n",
      "write result to diff2d-000033.vtu\n",
      "step=   33 finished!\n",
      "write result to diff2d-000034.vtu\n",
      "step=   34 finished!\n",
      "write result to diff2d-000035.vtu\n",
      "step=   35 finished!\n",
      "write result to diff2d-000036.vtu\n",
      "step=   36 finished!\n",
      "write result to diff2d-000037.vtu\n",
      "step=   37 finished!\n",
      "write result to diff2d-000038.vtu\n",
      "step=   38 finished!\n",
      "write result to diff2d-000039.vtu\n",
      "step=   39 finished!\n",
      "write result to diff2d-000040.vtu\n",
      "step=   40 finished!\n",
      "write result to diff2d-000041.vtu\n",
      "step=   41 finished!\n",
      "write result to diff2d-000042.vtu\n",
      "step=   42 finished!\n",
      "write result to diff2d-000043.vtu\n",
      "step=   43 finished!\n",
      "write result to diff2d-000044.vtu\n",
      "step=   44 finished!\n",
      "write result to diff2d-000045.vtu\n",
      "step=   45 finished!\n",
      "write result to diff2d-000046.vtu\n",
      "step=   46 finished!\n",
      "write result to diff2d-000047.vtu\n",
      "step=   47 finished!\n",
      "write result to diff2d-000048.vtu\n",
      "step=   48 finished!\n",
      "write result to diff2d-000049.vtu\n",
      "step=   49 finished!\n",
      "write result to diff2d-000050.vtu\n",
      "step=   50 finished!\n",
      "write result to diff2d-000051.vtu\n",
      "step=   51 finished!\n",
      "write result to diff2d-000052.vtu\n",
      "step=   52 finished!\n",
      "write result to diff2d-000053.vtu\n",
      "step=   53 finished!\n",
      "write result to diff2d-000054.vtu\n",
      "step=   54 finished!\n",
      "write result to diff2d-000055.vtu\n",
      "step=   55 finished!\n",
      "write result to diff2d-000056.vtu\n",
      "step=   56 finished!\n",
      "write result to diff2d-000057.vtu\n",
      "step=   57 finished!\n",
      "write result to diff2d-000058.vtu\n",
      "step=   58 finished!\n",
      "write result to diff2d-000059.vtu\n",
      "step=   59 finished!\n",
      "write result to diff2d-000060.vtu\n",
      "step=   60 finished!\n",
      "write result to diff2d-000061.vtu\n",
      "step=   61 finished!\n",
      "write result to diff2d-000062.vtu\n",
      "step=   62 finished!\n",
      "write result to diff2d-000063.vtu\n",
      "step=   63 finished!\n",
      "write result to diff2d-000064.vtu\n",
      "step=   64 finished!\n",
      "write result to diff2d-000065.vtu\n",
      "step=   65 finished!\n",
      "write result to diff2d-000066.vtu\n",
      "step=   66 finished!\n",
      "write result to diff2d-000067.vtu\n",
      "step=   67 finished!\n",
      "write result to diff2d-000068.vtu\n",
      "step=   68 finished!\n",
      "write result to diff2d-000069.vtu\n",
      "step=   69 finished!\n",
      "write result to diff2d-000070.vtu\n",
      "step=   70 finished!\n",
      "write result to diff2d-000071.vtu\n",
      "step=   71 finished!\n",
      "write result to diff2d-000072.vtu\n",
      "step=   72 finished!\n",
      "write result to diff2d-000073.vtu\n",
      "step=   73 finished!\n",
      "write result to diff2d-000074.vtu\n",
      "step=   74 finished!\n",
      "write result to diff2d-000075.vtu\n",
      "step=   75 finished!\n",
      "write result to diff2d-000076.vtu\n",
      "step=   76 finished!\n",
      "write result to diff2d-000077.vtu\n",
      "step=   77 finished!\n",
      "write result to diff2d-000078.vtu\n",
      "step=   78 finished!\n",
      "write result to diff2d-000079.vtu\n",
      "step=   79 finished!\n",
      "write result to diff2d-000080.vtu\n",
      "step=   80 finished!\n",
      "write result to diff2d-000081.vtu\n",
      "step=   81 finished!\n",
      "write result to diff2d-000082.vtu\n",
      "step=   82 finished!\n",
      "write result to diff2d-000083.vtu\n",
      "step=   83 finished!\n",
      "write result to diff2d-000084.vtu\n",
      "step=   84 finished!\n",
      "write result to diff2d-000085.vtu\n",
      "step=   85 finished!\n",
      "write result to diff2d-000086.vtu\n",
      "step=   86 finished!\n",
      "write result to diff2d-000087.vtu\n",
      "step=   87 finished!\n",
      "write result to diff2d-000088.vtu\n",
      "step=   88 finished!\n",
      "write result to diff2d-000089.vtu\n",
      "step=   89 finished!\n",
      "write result to diff2d-000090.vtu\n",
      "step=   90 finished!\n",
      "write result to diff2d-000091.vtu\n",
      "step=   91 finished!\n",
      "write result to diff2d-000092.vtu\n",
      "step=   92 finished!\n",
      "write result to diff2d-000093.vtu\n",
      "step=   93 finished!\n",
      "write result to diff2d-000094.vtu\n",
      "step=   94 finished!\n",
      "write result to diff2d-000095.vtu\n",
      "step=   95 finished!\n",
      "write result to diff2d-000096.vtu\n",
      "step=   96 finished!\n",
      "write result to diff2d-000097.vtu\n",
      "step=   97 finished!\n",
      "write result to diff2d-000098.vtu\n",
      "step=   98 finished!\n",
      "write result to diff2d-000099.vtu\n",
      "step=   99 finished!\n",
      "write result to diff2d-000100.vtu\n",
      "step=  100 finished!\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAZ7UlEQVR4nO3df4zc9X3n8eerm1h3+cGtCYFYtjm7e3vhrBZcOjWo6XU3pWltI2WDVCRDz7hA5LgXN0FXdLj5o43O9wdFpJzQEa82iVWQSnzoAseetIlrud3lIiDndbo1tqnLxiV4sWNvoLtkyx2nxe/7Y76bfDV8dz2zn5nvsOzrIY1mvp/35/OZ92dtzXu+3+/MfBURmJmZLdbPtTsBMzNb2lxIzMwsiQuJmZklcSExM7MkLiRmZpbkfe1OoB2uuOKKWLduXbvTMDNbUo4ePfrjiPhobfuyLCTr1q1jdHS03WmYmS0pkn5Y1O5DW2ZmlsSFxMzMkriQmJlZEhcSMzNL4kJiZmZJSi0kkjZLOiVpXNKegvjvSjqW3Z6VdN2lxkq6XNIhSS9l9yvLWo+ZmZVYSCR1AI8AW4ANwG2SNtR0+wegJyKuBfYCA3WM3QMcjohu4HC2bWZmJSnzeySbgPGIOA0g6QDQB5yc6xARz+b6Pw+sqWNsH9Cb9XsUGAbua8UCBgYGePzxx+eNnz17FkmsWrWqoVhqvF1jl2NerZzbeTmvVs8NcPvtt7Nz585544tRZiFZDZzJbU8ANyzQ/27g23WMvSoizgFExDlJVxZNJmknsBPg6quvbjh5gMcff5yxsTE2btxYGD99+jRA4T/iQrHUeLvGLse8Wjm383JerZ57bGwMoOmFRGVd2ErSrcBvR8Rns+3twKaI+IOCvp8Evgr8WkS8ttBYSVMR0Zkb+48RseB5kkqlEov5Zntvby8Aw8PDhfHOzmoaU1NTDcVS4+0auxzzauXczst5tXruS72GXYqkoxFRqW0v82T7BLA2t70GOFvbSdK1wNeBvoh4rY6x5yWtysauAi40OW8zM1tAmYXkCNAtab2kFcA2YDDfQdLVwJPA9oj4+zrHDgI7ssc7gKdbuAYzM6tR2jmSiJiVtBs4CHQA+yPihKRdWbwf+GPgI8BXJQHMRkRlvrHZ1PcDT0i6G3gFuLWsNZmZWcm//hsRQ8BQTVt/7vFngc/WOzZrfw24qbmZmplZvfzNdjMzS+JCYmZmSVxIzMwsiQuJmZklcSExM7MkLiRmZpbEhcTMzJK4kJiZWRIXEjMzS+JCYmZmSVxIzMwsiQuJmZklcSExM7MkLiRmZpbEhcTMzJK4kJiZWZJSC4mkzZJOSRqXtKcgfo2k5yS9JeneXPvHJY3lbm9IuieLfVnSq7nY1hKXZGa27JV2hURJHcAjwKeACeCIpMGIOJnr9jrwBeAz+bERcQrYmJvnVeCpXJeHIuLBliVvZmbzKvNSu5uA8Yg4DSDpANAH/LSQRMQF4IKkmxeY5ybgBxHxw1YmW2RkZASAzs7Owvj09PS88YViqfF2jV2OebVybuflvFo998zMDF1dXYWxFGUe2loNnMltT2RtjdoGfLOmbbekY5L2S1pZNEjSTkmjkkYnJycX8bRmZktfRDR9zjL3SFTQ1tCKJK0APg38Ua55H7A3m2sv8BXgrnc8UcQAMABQqVQW9Zfs6ekBYHh4uDA+9y5gamqqoVhqvF1jl2NerZzbeTmvVs/d29tb2J6qzD2SCWBtbnsNcLbBObYA34+I83MNEXE+It6OiIvA16geQjMzs5KUWUiOAN2S1md7FtuAwQbnuI2aw1qSVuU2bwGOJ2VpZmYNKe3QVkTMStoNHAQ6gP0RcULSrizeL+ljwChwGXAx+4jvhoh4Q9IHqH7i63M1Uz8gaSPVQ1svF8TNzKyFyjxHQkQMAUM1bf25xz+iesiraOybwEcK2rc3OU0zM2uAv9luZmZJXEjMzCyJC4mZmSVxITEzsyQuJGZmlsSFxMzMkriQmJlZEhcSMzNL4kJiZmZJXEjMzCyJC4mZmSVxITEzsyQuJGZmlsSFxMzMkriQmJlZEhcSMzNLUmohkbRZ0ilJ45L2FMSvkfScpLck3VsTe1nSC5LGJI3m2i+XdEjSS9n9yjLWYmZmVaUVEkkdwCPAFmADcJukDTXdXge+ADw4zzSfjIiNEVHJte0BDkdEN3A42zYzs5KUeandTcB4RJwGkHQA6ANOznWIiAvABUk3NzBvH9CbPX4UGAbua0K+7zAyMgJAZ2dnYXx6enre+EKx1Hi7xi7HvFo5t/NyXq2ee2Zmhq6ursJYijIPba0GzuS2J7K2egXwl5KOStqZa78qIs4BZPdXFg2WtFPSqKTRycnJBlM3M3tviIimz1nmHokK2hpZ0Sci4qykK4FDkv4uIp6pd3BEDAADAJVKZVF/yZ6eHgCGh4cL43PvAqamphqKpcbbNXY55tXKuZ2X82r13L29vYXtqcrcI5kA1ua21wBn6x0cEWez+wvAU1QPlQGcl7QKILu/0JRszcysLmUWkiNAt6T1klYA24DBegZK+qCkD889Bn4LOJ6FB4Ed2eMdwNNNzdrMzBZU2qGtiJiVtBs4CHQA+yPihKRdWbxf0seAUeAy4KKke6h+wusK4ClJczk/HhHfyaa+H3hC0t3AK8CtZa3JzMzKPUdCRAwBQzVt/bnHP6J6yKvWG8B188z5GnBTE9M0M7MG+JvtZmaWxIXEzMySuJCYmVkSFxIzM0viQmJmZklcSMzMLIkLiZmZJXEhMTOzJC4kZmaWxIXEzMySuJCYmVkSFxIzM0viQmJmZklcSMzMLIkLiZmZJSm1kEjaLOmUpHFJewri10h6TtJbku7Nta+V9NeSXpR0QtIXc7EvS3pV0lh221rWeszMrMQLW0nqAB4BPkX1+u1HJA1GxMlct9eBLwCfqRk+C/xhRHw/u+TuUUmHcmMfiogHW7sCMzMrUuYVEjcB4xFxGkDSAaAP+GkhiYgLwAVJN+cHRsQ54Fz2+CeSXgRW58eWYWRkBIDOzs7C+PT09LzxhWKp8XaNXY55tXJu5+W8Wj33zMwMXV1dhbEUZR7aWg2cyW1PZG0NkbQO+CXge7nm3ZKOSdovaeU843ZKGpU0Ojk52ejTmpm9J0RE0+csc49EBW0NrUjSh4BvAfdExBtZ8z5gbzbXXuArwF3veKKIAWAAoFKpLOov2dPTA8Dw8HBhfO5dwNTUVEOx1Hi7xi7HvFo5t/NyXq2eu7e3t7A9VZl7JBPA2tz2GuBsvYMlvZ9qEfmLiHhyrj0izkfE2xFxEfga1UNoZmZWkjILyRGgW9J6SSuAbcBgPQMlCfgG8GJE/FlNbFVu8xbgeJPyNTOzOpR2aCsiZiXtBg4CHcD+iDghaVcW75f0MWAUuAy4KOkeYANwLbAdeEHSWDbllyJiCHhA0kaqh7ZeBj5X1prMzKzccyRkL/xDNW39ucc/onrIq9Z3KT7HQkRsb2aOZmbWGH+z3czMkriQmJlZEhcSMzNL4kJiZmZJXEjMzCyJC4mZmSVxITEzsyQuJGZmlsSFxMzMkriQmJlZEhcSMzNL4kJiZmZJXEjMzCyJC4mZmSVxITEzsyQuJGZmlmTBQiLpxmY+maTNkk5JGpe0pyB+jaTnJL0l6d56xkq6XNIhSS9l9yubmbOZmS3sUnskz0jaKyn5SoqSOoBHgC1UL597m6QNNd1eB74APNjA2D3A4YjoBg5n22ZmVpJLFYgtwDeAmyX9u4g4mfBcm4DxiDgNIOkA0Af8dM6IuABckHRzA2P7gN6s36PAMHBfQp7zGhkZAaCzs7MwPj09PW98oVhqvF1jl2NerZzbeTmvVs89MzNDV1dXYSzFgnskEXEY+EXgb4BRSf8h4blWA2dy2xNZW+rYqyLiXJbvOeDKogkk7ZQ0Kml0cnKyocTNzN4rIqLpc17ykFVE/AS4W9IQ8N8k/SfgYk2fy+p4LhVNX1eWaWOrnSMGgAGASqWyqL9kT08PAMPDw4XxuXcBU1NTDcVS4+0auxzzauXczst5tXru3t7ewvZUdZ37kFQB/jPwEtXzF7OLeK4JYG1uew1wtgljz0taFRHnJK0CLiwiNzMzW6QFC0l2kv1PqJ5z+CqwJyL+7yKf6wjQLWk98CqwDbi9CWMHgR3A/dn904vMz8zMFuFSeyRHgMuBLdn5kkWLiFlJu4GDQAewPyJOSNqVxfslfQwYBS4DLkq6B9gQEW8Ujc2mvh94QtLdwCvArSl5mplZYy5VSI4DuyNiuhlPFhFDwFBNW3/u8Y+oHraqa2zW/hpwUzPyMzOzxi1YSCJie1mJmJnZ0uSfSDEzsyQuJGZmlsSFxMzMkriQmJlZEhcSMzNL4kJiZmZJXEjMzCyJC4mZmSVxITEzsyQuJGZmlsSFxMzMkriQmJlZEhcSMzNL4kJiZmZJXEjMzCxJqYVE0mZJpySNS9pTEJekh7P4MUnXZ+0flzSWu72RXT0RSV+W9GoutrXMNZmZLXeXukJi00jqAB4BPgVMAEckDUbEyVy3LUB3drsB2AfcEBGngI25eV4FnsqNeygiHmz5IszM7B1KKyTAJmA8Ik4DSDoA9AH5QtIHPBYRATwvqVPSqog4l+tzE/CDiPhhWYnPGRkZAaCzs7MwPj09PW98oVhqvF1jl2NerZzbeTmvVs89MzNDV1dXYSxFmYe2VgNnctsTWVujfbYB36xp250dCtsvaWXRk0vaKWlU0ujk5GTj2ZuZvQdU36c3V5l7JCpoq13Rgn0krQA+DfxRLr4P2Jv12wt8BbjrHZNEDAADAJVKZVF/yZ6eHgCGh4cL43PvAqamphqKpcbbNXY55tXKuZ2X82r13L29vYXtqcrcI5kA1ua21wBnG+yzBfh+RJyfa4iI8xHxdkRcBL5G9RCamZmVpMxCcgTolrQ+27PYBgzW9BkE7sg+vXUjMF1zfuQ2ag5rSVqV27wFON781M3MbD6lHdqKiFlJu4GDQAewPyJOSNqVxfuBIWArMA68Cdw5N17SB6h+4utzNVM/IGkj1UNbLxfEzcyshco8R0JEDFEtFvm2/tzjAD4/z9g3gY8UtG9vcppmZtYAf7PdzMySuJCYmVkSFxIzM0viQmJmZklcSMzMLIkLiZmZJXEhMTOzJC4kZmaWxIXEzMySuJCYmVkSFxIzM0viQmJmZklcSMzMLIkLiZmZJXEhMTOzJKUWEkmbJZ2SNC5pT0Fckh7O4sckXZ+LvSzpBUljkkZz7ZdLOiTppex+ZVnrMTOzEguJpA7gEarXXd8A3CZpQ023LUB3dtsJ7KuJfzIiNkZEJde2BzgcEd3A4WzbzMxKUuYVEjcB4xFxGkDSAaAPOJnr0wc8ll0p8XlJnZJW1Vy3vVYf0Js9fhQYBu5rcu4AjIyMANDZ2VkYn56enje+UCw13q6xyzGvVs7tvJxXq+eemZmhq6urMJaizENbq4Ezue2JrK3ePgH8paSjknbm+lw1V2iy+yuLnlzSTkmjkkYnJycTlmFmtnRV36c3V5l7JCpoq13RQn0+ERFnJV0JHJL0dxHxTL1PHhEDwABApVJZ1F+yp6cHgOHh4cL43LuAqamphmKp8XaNXY55tXJu5+W8Wj13b29vYXuqMvdIJoC1ue01wNl6+0TE3P0F4Cmqh8oAzktaBZDdX2h65mZmNq8yC8kRoFvSekkrgG3AYE2fQeCO7NNbNwLTEXFO0gclfRhA0geB3wKO58bsyB7vAJ5u9ULMzOxnSju0FRGzknYDB4EOYH9EnJC0K4v3A0PAVmAceBO4Mxt+FfCUpLmcH4+I72Sx+4EnJN0NvALcWtKSzMyMcs+REBFDVItFvq0/9ziAzxeMOw1cN8+crwE3NTdTMzOrl7/ZbmZmSVxIzMwsiQuJmZklcSExM7MkLiRmZpbEhcTMzJK4kJiZWRIXEjMzS+JCYmZmSVxIzMwsiQuJmZklcSExM7MkLiRmZpbEhcTMzJK4kJiZWRIXEjMzS1JqIZG0WdIpSeOS9hTEJenhLH5M0vVZ+1pJfy3pRUknJH0xN+bLkl6VNJbdtpa5JjOz5a60KyRK6gAeAT4FTABHJA1GxMlcty1Ad3a7AdiX3c8CfxgR38+u3X5U0qHc2Ici4sGy1mJmZj9T5qV2NwHj2WVzkXQA6APyhaQPeCy75O7zkjolrYqIc8A5gIj4iaQXgdU1Y1tuZGQEgM7OzsL49PT0vPGFYqnxdo1djnm1cm7n5bxaPffMzAxdXV2FsRRlHtpaDZzJbU9kbQ31kbQO+CXge7nm3dmhsP2SVhY9uaSdkkYljU5OTi5yCWZmS1v1fXpzlblHooK22hUt2EfSh4BvAfdExBtZ8z5gb9ZvL/AV4K53TBIxAAwAVCqVRf0le3p6ABgeHi6Mz70LmJqaaiiWGm/X2OWYVyvndl7Oq9Vz9/b2FranKnOPZAJYm9teA5ytt4+k91MtIn8REU/OdYiI8xHxdkRcBL5G9RCamZmVpMxCcgTolrRe0gpgGzBY02cQuCP79NaNwHREnJMk4BvAixHxZ/kBklblNm8BjrduCWZmVqu0Q1sRMStpN3AQ6AD2R8QJSbuyeD8wBGwFxoE3gTuz4Z8AtgMvSBrL2r4UEUPAA5I2Uj209TLwuVIWZGZmQLnnSMhe+Idq2vpzjwP4fMG471J8/oSI2N7kNM3MrAH+ZruZmSVxITEzsyQuJGZmlsSFxMzMkriQmJlZEhcSMzNL4kJiZmZJXEjMzCyJC4mZmSVxITEzsyQuJGZmlsSFxMzMkriQmJlZEhcSMzNL4kJiZmZJXEjMzCxJqYVE0mZJpySNS9pTEJekh7P4MUnXX2qspMslHZL0Una/sqz1mJlZiYVEUgfwCLAF2ADcJmlDTbctQHd22wnsq2PsHuBwRHQDh7NtMzMrSZmX2t0EjEfEaQBJB4A+4GSuTx/wWHbJ3ecldUpaBaxbYGwf0JuNfxQYBu5rxQJGRkYA6O3tLYxPT0/PG18olhpv19jlmFcr53ZezqvVc4+NjbFx48bCWIoyD22tBs7ktieytnr6LDT2qog4B5DdX1n05JJ2ShqVNDo5ObnoRZiZLVUbN27k9ttvb/q8Ze6RqKAt6uxTz9gFRcQAMABQqVQaGpubYzHDzMze08rcI5kA1ua21wBn6+yz0Njz2eEvsvsLTczZzMwuocxCcgTolrRe0gpgGzBY02cQuCP79NaNwHR2uGqhsYPAjuzxDuDpVi/EzMx+prRDWxExK2k3cBDoAPZHxAlJu7J4PzAEbAXGgTeBOxcam019P/CEpLuBV4Bby1qTmZmBluNx/0qlEqOjo+1Ow8xsSZF0NCIqte3+ZruZmSVxITEzsyQuJGZmlsSFxMzMkizLk+2SJoEfLnL4FcCPm5jOUuA1Lw9e8/KQsuZ/GREfrW1cloUkhaTRok8tvJd5zcuD17w8tGLNPrRlZmZJXEjMzCyJC0njBtqdQBt4zcuD17w8NH3NPkdiZmZJvEdiZmZJXEjMzCyJC8k8JG2WdErSuKR3XAc++6n7h7P4MUnXtyPPZqpjzb+brfWYpGclXdeOPJvpUmvO9fsVSW9L+p0y82u2etYrqVfSmKQTkkbKzrHZ6vh//S8k/U9Jf5ut+c525NlMkvZLuiDp+Dzx5r5+RYRvNTeqP1X/A+DngRXA3wIbavpsBb5N9eqNNwLfa3feJaz5V4GV2eMty2HNuX5/RfUyB7/T7rxb/G/cCZwErs62r2x33iWs+UvAn2aPPwq8Dqxod+6J6/514Hrg+Dzxpr5+eY+k2CZgPCJOR8T/Aw4AfTV9+oDHoup5oHPuSo1L1CXXHBHPRsQ/ZpvPU71S5VJWz78zwB8A32LpX32znvXeDjwZEa8ARMRyWHMAH5Yk4ENUC8lsuWk2V0Q8Q3Ud82nq65cLSbHVwJnc9kTW1mifpaTR9dxN9R3NUnbJNUtaDdwC9JeYV6vU82/8r4GVkoYlHZV0R2nZtUY9a/6vwL+hevnuF4AvRsTFctJrm6a+fpV2hcQlRgVttZ+TrqfPUlL3eiR9kmoh+bWWZtR69az5vwD3RcTb1TesS1o9630f8MvATcA/B56T9HxE/H2rk2uRetb828AY8BtAF3BI0v+KiDdanFs7NfX1y4Wk2ASwNre9huq7lUb7LCV1rUfStcDXgS0R8VpJubVKPWuuAAeyInIFsFXSbET8j1IybK56/1//OCL+CfgnSc8A1wFLtZDUs+Y7gfujevJgXNI/ANcA/7ucFNuiqa9fPrRV7AjQLWm9pBXANmCwps8gcEf26YcbgemIOFd2ok10yTVLuhp4Eti+hN+h5l1yzRGxPiLWRcQ64L8D/36JFhGo7//108C/lfQ+SR8AbgBeLDnPZqpnza9Q3QND0lXAx4HTpWZZvqa+fnmPpEBEzEraDRyk+qmP/RFxQtKuLN5P9RM8W4Fx4E2q72qWrDrX/MfAR4CvZu/QZ2MJ/3JqnWt+z6hnvRHxoqTvAMeAi8DXI6LwI6RLQZ3/xnuBP5f0AtVDPvdFxJL+aXlJ3wR6gSskTQB/ArwfWvP65Z9IMTOzJD60ZWZmSVxIzMwsiQuJmZklcSExM7MkLiRmZpbEhcTMzJK4kJi1kaSfk/SMpNovf34g++nzfe3KzaxeLiRmbZT9OODvAb8h6a5c6E+pfmH43nbkZdYIfyHR7F0g+6b1A8AvAv+K6jexeyPiu21NzKwOLiRm7xKSDlL9xd11wIGI+I/tzcisPi4kZu8SktZTvZrfD4BfiIi32pySWV18jsTs3eMu4P9Q/Unvn29zLmZ18x6J2buApF8BngU+Dfw+cBXwqxHxdlsTM6uD90jM2kzSPwMeA/48Ir4N7KR6wt3nSGxJ8B6JWZtJegj4DHBtRPwka9sGPAr88lK+HogtDy4kZm0k6deBvwJ+MyKGa2JPUD1XcmNEzLYhPbO6uJCYmVkSnyMxM7MkLiRmZpbEhcTMzJK4kJiZWRIXEjMzS+JCYmZmSVxIzMwsiQuJmZkl+f8MQDwpn4afMwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from FEToy.mesh.lagrange2dmesh import mesh2d\n",
    "from FEToy.fe.shapefun import shape1d\n",
    "from FEToy.fe.shapefun import shape2d\n",
    "from FEToy.fe.gaussrule import gausspoint1d\n",
    "from FEToy.fe.gaussrule import gausspoint2d\n",
    "from FEToy.postprocess.Result import ResultIO\n",
    "import numpy as np\n",
    "\n",
    "# define your mesh\n",
    "my2dmesh=mesh2d(nx=40,ny=8,xmax=1.0,ymax=0.2,meshtype='quad9')\n",
    "my2dmesh.createmesh()\n",
    "my2dmesh.plotmesh()\n",
    "\n",
    "#define your shape fun and gauss points\n",
    "shp1d=shape1d(meshtype='edge2')\n",
    "shp1d.update()\n",
    "\n",
    "# 4--7--3\n",
    "# |  |  |\n",
    "# 8--9--6 \n",
    "# |  |  |\n",
    "# 1--5--2\n",
    "shp2d=shape2d(meshtype='quad9')\n",
    "shp2d.update()\n",
    "\n",
    "qpoints1d=gausspoint1d(ngp=2)\n",
    "qpoints1d.creategausspoint()\n",
    "\n",
    "qpoints2d=gausspoint2d(ngp=2) # for quad9, you need ngp>=3!\n",
    "qpoints2d.creategausspoint()\n",
    "\n",
    "############################################\n",
    "# define your parameters\n",
    "D=0.5        # diffusion coefficient\n",
    "dt=1.0e-2    # time step size\n",
    "totalstep=100 # the total time step\n",
    "j0=0.005      # flux\n",
    "\n",
    "nDofs=my2dmesh.nodes*1\n",
    "K=np.zeros((nDofs,nDofs))\n",
    "F=np.zeros(nDofs)\n",
    "cold=np.zeros(nDofs) # solution in previous step\n",
    "c=np.zeros(nDofs)    # solution in current step\n",
    "\n",
    "# apply initial condition\n",
    "cold[:]=0.001\n",
    "# cold= np.random.normal(0.5, 0.02, nDofs)\n",
    "\n",
    "# define a csv output\n",
    "result=ResultIO('diff2d')\n",
    "# result.save2csv(mesh=my2dmesh,solution=cold,varnamelist=['c'],step=0)\n",
    "result.save2vtu(mesh=my2dmesh,solution=cold,varnamelist=['c'],step=0)\n",
    "for step in range(totalstep):\n",
    "    # init your k and F\n",
    "    K[:,:]=0.0;F[:]=0.0\n",
    "    for e in range(my2dmesh.elements):\n",
    "        elconn=my2dmesh.elementconn[e,:]\n",
    "        nodes=my2dmesh.nodecoords[elconn,:]\n",
    "        # now we can do the integration\n",
    "        for gp in range(qpoints2d.ngp2):\n",
    "            xi =qpoints2d.gpcoords[gp,1] # 1->xi, 2->eta\n",
    "            eta=qpoints2d.gpcoords[gp,2] # 1->xi, 2->eta\n",
    "            w  =qpoints2d.gpcoords[gp,0] # weight\n",
    "        \n",
    "            shp_val,shp_grad,j=shp2d.calc(xi,eta,nodes[:,0],nodes[:,1]) # xi,eta,x,y\n",
    "            JxW=j*w # the weight times the determinate of jacobian transformation\n",
    "            # now we can calculate the F term\n",
    "            for i in range(my2dmesh.nodesperelement):\n",
    "                iInd=elconn[i]\n",
    "                # assemble to global F\n",
    "                F[iInd]+=(cold[iInd]/dt)*shp_val[i]*JxW\n",
    "                # now we can calculate the K matrix\n",
    "                for j in range(my2dmesh.nodesperelement):\n",
    "                    jInd=elconn[j]\n",
    "                    K[iInd,jInd]+=(1.0/dt)*shp_val[j]*shp_val[i]*JxW\\\n",
    "                                 +D*(shp_grad[j,0]*shp_grad[i,0]\\\n",
    "                                    +shp_grad[j,1]*shp_grad[i,1])*JxW\n",
    "#     # apply the flux, now you need the surface integration\n",
    "#     bcelmt=my2dmesh.bcconn['right'] # we want to apply the flux on the right side\n",
    "#     for e in range(np.shape(bcelmt)[0]):\n",
    "#         elconn=bcelmt[e,:]\n",
    "#         nodes=my2dmesh.nodecoords[elconn,:]\n",
    "#         for gp in range(qpoints1d.ngp):\n",
    "#             xi =qpoints1d.gpcoords[gp,1] # 1->xi, 2->eta\n",
    "#             w  =qpoints1d.gpcoords[gp,0] # weight\n",
    "        \n",
    "#             shp_val,shp_grad,j=shp1d.calc(xi,nodes[:,1]) # for left side, you need 'y'\n",
    "#             JxW=j*w # the weight times the determinate of jacobian transformation\n",
    "#             # now we can calculate the F term\n",
    "#             for i in range(len(elconn)):\n",
    "#                 iInd=elconn[i]\n",
    "#                 # assemble to global F\n",
    "#                 F[iInd]+=j0*shp_val[i]*JxW\n",
    "    ####\n",
    "    # Appyl boundary conditions\n",
    "    y2=5.0\n",
    "    iInd=my2dmesh.bcnodeids['left']\n",
    "    K[iInd,iInd]+=1.0e16\n",
    "    F[iInd]+=y2*1.0e16\n",
    "\n",
    "    iInd=my2dmesh.bcnodeids['right']\n",
    "    K[iInd,iInd]+=1.0e16\n",
    "    F[iInd]+=y2*1.0e16\n",
    "\n",
    "    iInd=my2dmesh.bcnodeids['bottom']\n",
    "    K[iInd,iInd]+=1.0e16\n",
    "    F[iInd]+=y2*1.0e16\n",
    "\n",
    "    iInd=my2dmesh.bcnodeids['top']\n",
    "    K[iInd,iInd]+=1.0e16\n",
    "    F[iInd]+=y2*1.0e16\n",
    "        \n",
    "    \n",
    "    # solve the equation\n",
    "    c=np.linalg.solve(K,F)\n",
    "    cold[:]=c[:] # update solution\n",
    "    # result.save2csv(mesh=my2dmesh,solution=cold,varnamelist=['c'],step=step+1)\n",
    "    result.save2vtu(mesh=my2dmesh,solution=cold,varnamelist=['c'],step=step+1)\n",
    "    print('step=%5d finished!'%(step+1))"
   ]
  }
 ],
 "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
}