{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For $\\phi_{tot}$ we obtain \n",
    "\n",
    "$\\partial_\\mathrm{t}\\phi_\\mathrm{tot}=k_\\mathrm{B}T\\Gamma_\\mathrm{0} \\,\\nabla\\cdot\\left[\\left[1-2\\chi\\phi_\\mathrm{tot}(1-\\phi_\\mathrm{tot})\\right]\\nabla \\phi_\\mathrm{tot}-\\phi_\\mathrm{tot}(1-\\phi_\\mathrm{tot})\\kappa\\nabla\\nabla^2\\phi_\\mathrm{tot}\\right]$\n",
    "\n",
    "We use the [Fenics Cahn-Hilliard example](https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/cahn-hilliard/demo_cahn-hilliard.py.rst#rst-header-id1) (the maths is explained [here](https://fenicsproject.org/docs/dolfinx/dev/python/demos/cahn-hilliard/demo_cahn-hilliard.py.html)) and replace $\\mu$:\n",
    "\n",
    "\n",
    "$\\mu = \\left[1-2\\chi\\phi_\\mathrm{tot}(1-\\phi_\\mathrm{tot})\\right]\\nabla \\phi_\\mathrm{tot}-\\phi_\\mathrm{tot}(1-\\phi_\\mathrm{tot})\\kappa\\nabla\\nabla^2\\phi_\\mathrm{tot}$\n",
    "\n",
    "\n",
    "Cast into the weak form we obtain:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extended Standard Cahn-Hilliard Example to Binary Flory Huggins \n",
    "# as discussed on 28/11/2019\n",
    "# Example can be found at \n",
    "# https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/\n",
    "# documented/cahn-hilliard/demo_cahn-hilliard.py.rst#rst-header-id1\n",
    "# Runs with fenics 2019.01\n",
    "# The resulting .pvd file can be opened using default settings in ParaView\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "import random\n",
    "from dolfin import *\n",
    "# Class representing the intial conditions\n",
    "class InitialConditions(UserExpression):\n",
    "    def __init__(self, **kwargs):\n",
    "        random.seed(2 + MPI.rank(MPI.comm_world))\n",
    "        super().__init__(**kwargs)\n",
    "    def eval(self, values, x):\n",
    "        if (x[0] > 0.45):# & (x[0] < 0.55):# & (x[1] > 0.45) & (x[1] < 0.55):\n",
    "#             values[0] = 0.63 + 0.02*(0.5 - random.random())\n",
    "            values[0] = 0.6\n",
    "        else:\n",
    "            values[0] = 0.35                \n",
    "        values[1] = 0.0\n",
    "#         values[0] = 0.63 + 0.02*(0.5 - random.random())\n",
    "        values[1] = 0.0\n",
    "    def value_shape(self):\n",
    "        return (2,)\n",
    "# Class for interfacing with the Newton solver\n",
    "class CahnHilliardEquation(NonlinearProblem):\n",
    "    def __init__(self, a, L):\n",
    "        NonlinearProblem.__init__(self)\n",
    "        self.L = L\n",
    "        self.a = a\n",
    "    def F(self, b, x):\n",
    "        assemble(self.L, tensor=b)\n",
    "    def J(self, A, x):\n",
    "        assemble(self.a, tensor=A)\n",
    "# Model parameters\n",
    "lmbda  = 10e-04  # surface parameter\n",
    "dt     = 10.0e-04  # time step\n",
    "X = 2.2 # Chi Flory Huggins parameter\n",
    "# time stepping family, e.g.: \n",
    "# theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson\n",
    "theta  = 0.5      \n",
    "# Form compiler optionsmai\n",
    "parameters[\"form_compiler\"][\"optimize\"]     = True\n",
    "parameters[\"form_compiler\"][\"cpp_optimize\"] = True\n",
    "# Create mesh and build function space\n",
    "x_mesh = 96\n",
    "y_mesh = 1\n",
    "mesh = UnitSquareMesh.create(x_mesh, y_mesh, CellType.Type.quadrilateral)\n",
    "P1 = FiniteElement(\"Lagrange\", mesh.ufl_cell(), 1)\n",
    "ME = FunctionSpace(mesh, P1*P1)\n",
    "# Define trial and test functions\n",
    "du    = TrialFunction(ME)\n",
    "q, v  = TestFunctions(ME)\n",
    "# Define functions\n",
    "u   = Function(ME)  # current solution\n",
    "u0  = Function(ME)  # solution from previous converged step\n",
    "\n",
    "# Split mixed functions\n",
    "dc, dmu = split(du)\n",
    "c,  mu  = split(u)\n",
    "c0, mu0 = split(u0)\n",
    "# Create intial conditions and interpolate\n",
    "u_init = InitialConditions(degree=1)\n",
    "u.interpolate(u_init)\n",
    "u0.interpolate(u_init)\n",
    "# Compute the chemical potential df/dc\n",
    "c = variable(c)\n",
    "dfdc = ln(c/(1-c))+(1-2*X*c)\n",
    "# mu_(n+theta)\n",
    "mu_mid = (1.0-theta)*mu0 + theta*mu\n",
    "# Weak statement of the equations\n",
    "L0 = c*q*dx - c0*q*dx + dt*c*(1-c)*dot(grad(mu_mid), grad(q))*dx\n",
    "L1 = mu*v*dx - (ln(c/(1-c))+X*(1-2*c))*v*dx - lmbda*dot(grad(c), grad(v))*dx\n",
    "L = L0 + L1\n",
    "# Compute directional derivative about u in the direction of du (Jacobian)\n",
    "a = derivative(L, u, du)\n",
    "# Create nonlinear problem and Newton solver\n",
    "problem = CahnHilliardEquation(a, L)\n",
    "solver = NewtonSolver()\n",
    "solver.parameters[\"linear_solver\"] = \"lu\"\n",
    "solver.parameters[\"convergence_criterion\"] = \"incremental\"\n",
    "solver.parameters[\"relative_tolerance\"] = 1e-6\n",
    "# Output file\n",
    "file = File(\"output.pvd\", \"compressed\")\n",
    "\n",
    "# Step in time\n",
    "t = 0.0\n",
    "T = 500*dt\n",
    "while (t < T):\n",
    "    file << (u.split()[0], t)\n",
    "    t += dt\n",
    "    u0.vector()[:] = u.vector()\n",
    "    solver.solve(problem, u.vector())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1279f4cf8>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD4CAYAAAANbUbJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxV9Z3/8dfn3ux7QhZCFhJCIIRdAggKikvFWpcutmLVrlo62qrTaaudX2c603amnS5jF1trXVutVq2tVq24FVDZIwgkIZCdkD0h+3aX7++PG51ICSQh4eTe83k+HveRnHPPvffzJeHNl+/5nu8RYwxKKaUCm8PqApRSSk0+DXullLIBDXullLIBDXullLIBDXullLKBIKs+ODEx0WRlZVn18Uop5ZcKCwtbjDFJY32dZWGflZXFnj17rPp4pZTySyJSPZ7X6TCOUkrZgIa9UkrZgIa9UkrZgIa9UkrZgIa9UkrZgIa9UkrZgIa9UkrZgGXz7JVSY+PxGho6+zna1kvt8T66+l0MuL30uzw4RAgPdhIe4iQmPJjk6FBSYsKYHhNGeIjT6tLVFKBhr9QU5fUa3qk5zvbyVnZUtlJYfZx+l3fM75MSE0p2YiTZiVHMS40mPzWGvNQYokL1r7+d6E9bqSmmtXuAZwpr+cOuGqpbewGYlxrDdcszmZMSTUZCOOnxEcRHBBMa5CQkyIHXGPpcHvoHPbT3uWjqHKCxs5/6jj4qW3qpbOnmpQP1PLGr5v3PmZ0cxZKMOJZmxlEwM4Hc5CgcDrGq2WqSadgrNUV0D7j531cP8/vt1Qx6vKzITuDOS+ZwwZwk4iNDTvlaJ0Kw00FMWDDJMWHMSYn+h2OM8Q0DFdd1UlzXybu17bxxqIlnCmsBiI8IZkV2AqtzElmTm0h2YiQiGv6BQsNeKYsZY3jxQD3ffaGYpq4Brl2Wzs1rZpF7ksA+EyJCamw4qbHhXDwv5f3PrmnrZVdlGzsq2thZ2cqmokYA0uPDWTsniUvnpbAqZxphwTr278/EqnvQFhQUGF0ITdld36CHrz29j5cONDB/Rgzfu2YBSzPjLa2purWHrUda2Hq4mbfLWugd9BAR4mRtbhIfXpTKxXnJROp4v2VEpNAYUzDm140m7EVkPfAzwAk8YIz5wQnPfx349NBmEDAPSDLGtI30nhr2yu6auvq5+dE97D/WwTfX53Hzmlk4p9iYeb/Lw46KVl4vaeKV4gYaOwcIDXJwUV4y1yxNY93cZEKCdAb32TRpYS8iTuAwcClQC+wGNhhjikc4/krgTmPMRad6Xw17ZWelDV18/pHdtPUM8vMNS7k0P8Xqkk7L6zUU1hznxf31vLC/npbuARIiQ7hq8Qw2rMhk7vSJHXZSJzfesB/N/8VWAGXGmIqhD3oSuBo4adgDG4AnxlqIUnZR09rL9b/dgdMhPPWlVSxMj7W6pFFxOITlWQksz0rg/10xjzePtPhmDe2s4ZFtVazISuCGVTNZP3+69vanoNGEfRpwdNh2LbDyZAeKSASwHrhthOdvAW4ByMzMHFOhSgWCjl4Xn3tkF26v4emNq5iVFGV1SeMS5HSwLi+ZdXnJHO8Z5OnCozy2o4avPrGX6TFhfHFNNtetyNS5/FPIaP75Pdkg4khjP1cCb480Vm+Mud8YU2CMKUhKGvNdtZTya4NuLxsfK6SmrZf7b1zmt0F/ovjIEG5Zm8Pmf7mQhz+7nKzECL73Ygnn/eAN7nntMF39LqtLVIyuZ18LZAzbTgfqRjj2OnQIR6mT+vZfDrK9opWffnIxK2dNs7qcCedwyPu9/b01x/n15nLuee0Ij26r4tZ1s7nh3Jk6fdNCo+nZ7wZyRSRbRELwBfrzJx4kIrHABcBzE1uiUv7v5YP1/HHPUW5bN5uPnZNudTmTbmlmPPffVMDzt53HgrRYvvdiCRf/ZAsvH2zAqunednfasDfGuPGNwW8CSoCnjDFFIrJRRDYOO/SjwCvGmJ7JKVUp/9TR6+LbzxWRnxrD7ZfkWl3OWbUoPY7ff2Elf/jiSqLDgtj4WCGfe2Q3VS0aE2ebXlSl1CT7+tPv8uzeYzx3q6+Xa1duj5dHtlVxz2tHGPR4+cZlc/n8edm6Hs8YjXfqpc6PUmoSbT3czNOFtXxp7SxbBz34ZvB8cc0sXv/aBazNTeR7L5bw6Qd2cqy9z+rSbEHDXqlJ0jfo4e5nDzArKZKvXmyv4ZtTSYkJ47c3FfCDjy1kf2076/93Ky8fbLC6rICnYa/UJHlsRzXH2vv4/jULdRbKCUSE61Zk8rfb1zIrKZKNjxXyk1dK8Xj15O1k0bBXahL0Drq5b0s5a3ITWZUTeNMsJ0rmtAj++KVVXLssnV+8UcYXH91NR5/Oy58MGvZKTYLfba+mtWeQOy6ZY3UpU15YsJP/+cQivnvNAt4qa+FTv9lOU2e/1WUFHA17pSZY94Cb32wp54I5SSybae1yxf5CRLjx3Jk8/NkV1LT18vH7tlHdqtMzJ5KGvVIT7NFtVRzvdXHnpdqrH6vzcxN54uZz6e538/Ffb6e4rtPqkgKGhr1SE6ir38X9Wyu4KC+ZJRlxVpfjlxZnxPH0xtWEOIVPP7CDI41dVpcUEDTslZpAT+2ppaPPxe061fKMzE6O4g83n0uQ08END+6kZujG62r8NOyVmiDGGB7fWc3SzDgWa6/+jGUlRvLYF1Yy4PZy/QM7aOjQk7ZnQsNeqQmyo6KNiuYePr1yptWlBIy506N59HMraO91cdNDO3W55DOgYa/UBHl8ZzUxYUF8ZFGq1aUElMUZcfzmxmWUN/dwx5P79MKrcdKwV2oCtHQPsKmogU8sy9CrZSfBebMT+c6V+bx+qIn/2XTI6nL8kt4zTKkJ8PSeWlwew/Ur9Xabk+XGVVkcbuzmN1sqmJMczceXBf59ASaS9uyVOkNer+EPu6o5d1YCs5MD41aDU9W/XZnP6pxp3P3nAzoHf4y0Z6/UGXqzrIWjbX1847K8M3sjY2CgC/raYKAbBoce7gFw9/u+egbB4wKvB7xuMB7f98bjuzO08f7fA+N7z/duGf2Be1eYD36unwgGHkzz8FT9Ud595Clyl6UT7PDDPmvW+ZB76Vn9SA17pc7Qn9+pJS4imMvmTz/9wV4vtFXAsUJoKYXjVXC8GjrroLfFF+YTRkBk2Nehfe8/LR881k+EAzdicA0Y2CHgjzc/cQRp2CvlT/pdHl4raeKKhamEBI3Qw+xrh0MvQvFzcHQH9Hf49osT4jIgPgtyLoLIRIiYBhEJEBoNIVG+R3AYBIWBM2ToEewLC4fT9x7vfRXH0ENOCPLA4wB+vukQ9/69nF9sWMqVi2dYXdKUp2Gv1BnYeriZ7gE3V5xsumXdXtj6Yzjyiq/HHpsJ8z8Kact8j8S54NS/guN1xyVz2F7eyreePcDSzDjS4yOsLmlK0980pc7ASwfqiYsI/uCa9W2V8Mb34OAzEB4Py2+GBR/zBXyA97jPpmCng59dt5TL7tnKt/58kEc/txzRP98RadgrNU79Lg+vFjdy5eIZBDuHhnD2PAwvfd03zLLma3De7RBm73vPTqaMhAi+cdlcvvPXYp5955hOxzwFPzyNrdTUsOVwMz2DHt8QjtcDm/4VXrgDZl0IX30HLv43Dfqz4KZVWRTMjOc/XyimqUvXzxmJhr1S4/TSgXriI4JZlR4GT90E238JK26BDU9CjJ4wPFscDuGHn1hEn8vDvz9XZHU5U5aGvVLj0O/y8FpxI+vnpxD0/Jeh9CVY/0P48I/0pKsFcpKiuOOSXP52sIFNRQ1WlzMladgrNQ6bS31DOF8MeQVK/gqX/Aecu9HqsmztljWzmJsSzXdfKKbf5bG6nClHw16pcXj5YD1rwquYtfcHMPcKWP0Vq0uyvSCng3+/Kp/a4338ZkuF1eVMORr2So2Rx2vYW1rBz5w/Q2JmwDX36pTKKWJ1TiJXLEzlV5vLqD2ud7caTsNeqTHaX9vO7e6HiPW2wbWP+ubSqynjW1fMQwS+/2KJ1aVMKRr2So3R/n27uMbxNoMFGyHtHKvLUSdIiwvnny6czd8ONvB2WYvV5UwZGvZKjVF20b0MSCjhF9xpdSlqBLesnUVaXDj/9VIJXr2zFaBhr9SYtFcf4Pz+rRxM/xRETjv9C5QlwoKd/Mtlcyiq6+Sv++usLmdK0LBXagy6X/kvegkl7II7rC5FncbVi9OYlxrDjzaVMuDWqZga9kqNVlMJM479jT/K5eTnZFtdjToNh0O46/I8ao/38fiOGqvLsZyGvVKjZLb+mD5CKZ/9GZz+eMMMG1qbm8h5s6fxizeO0NnvsrocS2nYKzUafe2Y4ud5yn0By+fnWl2NGiUR4a718zje6+K3W+19oZWGvVKjcegFHN5BnvOex9rcJKurUWOwMD2WKxal8tBblbT1TORtH/3LqMJeRNaLSKmIlInIXSMcc6GI7BORIhHZMrFlKmWxA8/Q4JyON/UcpkWFWl2NGqM7L8ml1+Xhfhv37k8b9iLiBO4FLgfygQ0ikn/CMXHAr4CrjDHzgWsnoValrNHdjKncwrOD53K+9ur90uzkaK5aPINHt1XR0j1gdTmWGE3PfgVQZoypMMYMAk8CV59wzPXAs8aYGgBjTNPElqmUhYr/ghgvf3Gv4txZOrfeX3314lwG3B5+s6Xc6lIsMZqwTwOODtuuHdo33BwgXkQ2i0ihiNw0UQUqZbmDf6I5PIcKyWTZTF0Hx1/lJEVxzdI0fre9mqZO+93RajRhf7I5ZidefxwELAOuAC4Dvi0ic/7hjURuEZE9IrKnubl5zMUqdda1H4Wa7bziOI9F6bFEhuqNSfzZVy/Kxe01/Gqz/Xr3own7WiBj2HY6cOL1x7XAy8aYHmNMC7AVWHziGxlj7jfGFBhjCpKSdOxT+YGiZwF48Pg5OoQTALISI/nY0jSe2FVju/vVjibsdwO5IpItIiHAdcDzJxzzHLBGRIJEJAJYCej6osr/HfwTXdMWUeFNZqWGfUC4dd1sXB4vD7xZaXUpZ9Vpw94Y4wZuAzbhC/CnjDFFIrJRRDYOHVMCvAzsB3YBDxhjDk5e2UqdBV0NUP8u70SuxekQCnS8PiBkJUZy1eIZPLaj2lbz7kc1z94Y85IxZo4xJscY8/2hffcZY+4bdsyPjDH5xpgFxph7Jqtgpc6aqrcAeLE7V8frA8yt62bTO+jhobfs07vXK2iVGknlVkxoDH9pSNTx+gCTmxLN5Qum8+i2Kjr67LFmjoa9UiOpeou2xAIGvaJhH4Buu2g2XQNufretyupSzgoNe6VOprMO2srZH7RQx+sD1PwZsVycl8yDb1fSPeC2upxJp2Gv1MlUvgnAi12zdbw+gN120Wzae108sTPw17vXsFfqZKrexITF8VxDAiuzdQgnUC3NjGd1zjR++2YF/a7AvpuVhr1SJ1P1JseTluPyCsuzdAgnkN22bjZNXQM8U1hrdSmTSsNeqRO1H4XjVRSHLgJ8vT8VuFblTGNJRhz3bSnH7fFaXc6k0bBX6kRD8+tf65vLrKRIEiJDLC5ITSYR4dZ1s6k93sfz7564Ekzg0LBX6kRVb2LCE3i+LpZl2qu3hYvzksmbHs2vNpfj9Z64zmNg0LBX6kSVb9KTei5tfR5d0tgmHA7hyxfmUNbUzaaiBqvLmRQa9koN114DHTUcDvct2lqgJ2dt4yOLZpCdGMkv/16GMYHXu9ewV2q4Y+8AsK0/m9jwYGYlRllckDpbnEO9+6K6TjaXBt79NjTslRqufh84gnipKYFzMuNwOE527x4VqD66NI20uHB+8caRgOvda9grNVzdPjyJ8yhuHtTxehsKdjrYeGEO79S0s72i1epyJpSGvVLvMQbq99EYlQfAORr2tnTtsnSSo0P55RtlVpcyoTTslXpPew30HafIZON0CEsy4qyuSFkgLNjJLWtnsa28lcLqNqvLmTAa9kq9p34fAFu60shPjSEiRBc/s6vrV2aSGBXCPa8dsbqUCaNhr9R76vZhHEH8tTFBx+ttLiIkiC+tzeHNIy3srgqM3r2GvVLvqd9Hf9wcOlxOlmbqEI7d3XDuTBKjQvnfVw9bXcqE0LBXCnwnZ+v2cSx8LgDn6DIJthce4uTLF+awrbyVHQEwM0fDXimAjqPQ18a7niymRYaQHh9udUVqCvj0ykySowOjd69hrxRAne/k7BudqSzJiENEL6ZSvpk5X74wh52VbWwra7G6nDOiYa8UQP0+jDh5rS1Zp1yqD9iwIpPU2DB++PIhv14RU8NeKYC6ffTE5jJACEv05KwaJizYyT9fOod3azt44UC91eWMm4a9UkNXztaE5gKwKF3DXn3Qx85JJ296ND/adIgBt3/eq1bDXqmOWuht5R1XFjlJkcSGB1tdkZpinA7hWx+ex9G2Pn6/vdrqcsZFw16p+ncBeK19OksydMqlOrm1c5JYk5vIL94oo6PXZXU5Y6Zhr1RTMQC7elN1vF6d0t2Xz6Oz38Uv3vC/ZRQ07JVqKqYnMoNewliqM3HUKeTPiOFTBRk8vK2K4rpOq8sZEw17pZpKqA3KIjTIwdzp0VZXo6a4uy7PIy48mLv/fACPH03F1LBX9uYegNYyDrhmsDAtlmCn/pVQpxYXEcK3P5LPu0fbeWyH/5ys1d9sZW8tR8Dr5u0uvZhKjd7VS2awJjeRH20qpaGj3+pyRkXDXtlbUwkARe50PTmrRk1E+N41C3B5vHz7uYN+cb9aDXtlb03FeCWISpPKYr2YSo3BzGmRfP2yubxa3Mgj26qsLue0NOyVvTUV0xiSQWxUhK50qcbsC+dnc2l+Ct9/sYTC6uNWl3NKGvbK3pqKOeRJZ3G6rnSpxk5E+PG1i0mNC+O2P7xDW8+g1SWNSMNe2ddAF7TXUNifymI9OavGKTY8mF9/ehmtPYPc/uReBt1eq0s6qVGFvYisF5FSESkTkbtO8vyFItIhIvuGHv828aUqNcGaSwEo9aZr2KszsiAtlu9dvYA3j7TwT48XTsnF0k4b9iLiBO4FLgfygQ0ikn+SQ980xiwZevznBNep1MRrLAKg1GSwOD3W4mKUv/vk8gy+e/V8Xitp4ubfFdLvmlqBP5qe/QqgzBhTYYwZBJ4Erp7cspQ6C5pKGJAwnPEziYsIsboaFQBuXJXFDz++kDePNPO5h3fT3DVgdUnvG03YpwFHh23XDu070SoReVdE/iYi80/2RiJyi4jsEZE9zc3N4yhXqQnUVEwZ6SzKTLC6EhVAPrU8k59+cjF7qtu46CebefjtStwe68fxg0ZxzMmmKJx4BcE7wExjTLeIfBj4C5D7Dy8y5n7gfoCCgoKpfxWCCmiexiKKXPP1ylk14T66NJ1F6XF85/ki/uOvxfxx91GuXOxbkmNhWizxkWf/f5KjCftaIGPYdjpQN/wAY0znsO9fEpFfiUiiMca/79CrAld3M87eFkpNOldo2KtJkJMUxe8+v4JNRQ389NXD/GhT6fvPbbwgh7suzzur9Ywm7HcDuSKSDRwDrgOuH36AiEwHGo0xRkRW4Bseap3oYpWaMM2+ZRLKyCQ/NcbiYlSgEhHWL0hl/YJUOvpcFB3rYP+xDhamnf0JAacNe2OMW0RuAzYBTuAhY0yRiGwcev4+4BPAl0XEDfQB1xl/WCxC2dfQmjgmOZ+wYKfFxSg7iA0PZvXsRFbPTrTk80fTs8cY8xLw0gn77hv2/S+BX05saUpNHtNYTKeJJDMzy+pSlDor9ApaZUv9dUWUmnQW6z1nlU1o2Cv7MQZnaylHvOk6E0fZhoa9sp+uBkJcnVQ7M5mVFGV1NUqdFRr2yn6GZuJ4k/JwOnSlS2UPGvbKdlwNxQDEZS60uBKlzp5RzcZRKpC0V+3HaaKYM2uW1aUoddZo2Cvb8TaVUGHSWTJTZ+Io+9BhHGUvxhDTVc6x4CySo8Osrkaps0bDXtlLZx3h3h5c0+ZYXYlSZ5WGvbKV9pr9AESl68lZZS8a9spWGsreBSBtzlKLK1Hq7NKwV7YyWFdEq4khL0dn4ih70bBXthLefoRjwTN1pUtlOxr2yjY8Hi/TB6vojfuHm6gpFfA07JVtVFUeIVr6CE096S2SlQpoGvbKNmpLCwFInr3E4kqUOvs07JVtdNceBCA1R2fiKPvRsFe24Ww5RLsjHkfUNKtLUeqs07BXttDR62LGQAWdMXpyVtmThr2yhb3VzcyVWpypeuWssicNe2ULFYf2EyouEnOWWV2KUpbQsFe20F2zD4DQ9EUWV6KUNTTsVcBzebyEtRbjkSBInGt1OUpZQsNeBbxD9V3kmiq6o3MgKMTqcpSyhIa9Cnh7qtvIcxwlOE1Pzir70rBXAe9QRRWp0kZE+mKrS1HKMhr2KuD11PjWsGf6AmsLUcpCGvYqoB1r7yOl94hvI0WHcZR9adirgFZYfZw8qcEVngRRSVaXo5RlNOxVQCusamO+s0avnFW2p2GvAtreqmZy5RgOHa9XNqdhrwJWR5+LgcZSgnHBdO3ZK3vTsFcBa1dlG3Op8W2kaM9e2ZuGvQpYOypaWRBUg3GGQKIubazsTcNeBazt5a2sDK9DkuaCM9jqcpSylIa9CkjtvYOUNHSQ6ymH6XrlrFKjCnsRWS8ipSJSJiJ3neK45SLiEZFPTFyJSo3dzso2Mmgiwt0O6QVWl6OU5U4b9iLiBO4FLgfygQ0ikj/CcT8ENk10kUqN1fbyVlYGl/k20pdbW4xSU8BoevYrgDJjTIUxZhB4Erj6JMd9BfgT0DSB9Sk1LjsqWrk0ugZCoiB5ntXlKGW50YR9GnB02Hbt0L73iUga8FHgvlO9kYjcIiJ7RGRPc3PzWGtValTaegY51NDFEjkCaeeAw2l1SUpZbjRhLyfZZ07Yvgf4pjHGc6o3Msbcb4wpMMYUJCXpOiVqcuyqbCWMAZJ6jugQjlJDgkZxTC2QMWw7Hag74ZgC4EkRAUgEPiwibmPMXyakSqXGYHt5KwXB1YjxQPoKq8tRakoYTdjvBnJFJBs4BlwHXD/8AGNM9nvfi8gjwAsa9MoqOyra+FxCLXSgM3GUGnLaYRxjjBu4Dd8smxLgKWNMkYhsFJGNk12gUmPR1NlPaWMXK4LLIWEWRCZaXZJSU8JoevYYY14CXjph30lPxhpjPnvmZSk1PpsPNwOGjJ4iyF1ndTlKTRl6Ba0KKFtKm1kc3UlwX5OenFVqGA17FTBcHi9bjzTzqZQG3w4Ne6Xep2GvAsY71cfp6nezOqwSgsIhZb7VJSk1ZWjYq4Cx+XAzQQ4ho+cgzFiqK10qNYyGvQoYfz/UxHmZ4TgbD0CGDuEoNZyGvQoIDR39HGro4rqkSvC6YJbOxFFqOA17FRA2l/rW31vlKfQtfjZztcUVKTW1aNirgLC5tJkZMaHE1m6GWRdCUKjFFSk1tWjYK7836PbyVlkLn5zZjXTWQu6HrC5JqSlHw175vR0VrXQPuLki7IBvR+6l1hak1BSkYa/83gv764gKDSKn/W2YvhBiZlhdklJTjoa98muDbi8vH2zgqrkROGp36RCOUiPQsFd+7a2yZjr73WyYVgbGA7mXWV2SUlOShr3yay+8W09MWBD53TsgPF7Xr1dqBBr2ym/1uzy8UtzI+vxknOWvwexL9H6zSo1Aw175rS2Hm+kecHP9jDrobdEhHKVOQcNe+a0X9tcTHxHMovpnICwW8q6wuiSlpiwNe+WX+gY9vF7SyCfmhuIo+Sssvh5CIqwuS6kpS8Ne+aXXShrpHfRwQ+hW38JnBZ+3uiSlpjQNe+WXHt9ZTWZcCJmVf4SsNZA0x+qSlJrSNOyV3znS2MWOija+mXsM6TiqvXqlRkHDXvmdx3fWEOJ0cGnPCxCZDHkfsbokpaY8DXvlV3oG3PypsJYb5hpCKl6Dc26CoBCry1JqyguyugClxuK5fXV0Dbi5lT/67jG77LNWl6SUX9CevfIbxhge21HNJ6dVMq3iOTjvDojLsLospfyC3/Xsy5u7ea24EadDCHIITqeDiGAnkaFBRIY6iQsPIT4ymITIECJC/K556hTeqWnnSH0bT0x7EOJmwpp/trokpfyG36VhSX0n//23Q6M6Nio0iOmxYaTGhjEjNpyZiRHMTIgkKzGCnKQowoJ1HRV/cv/Wcr4cuonYngrY8EcIDre6JKX8ht+F/eULUin+z2TcXoPHY3B5vPS5PHQPuOkZ8NDeO8jx3kFaewZp7hqgoaOf+o5+Xq9vpKV78P33EYHMhAhyk6PJT40mf0YsC9JiSIsLR0QsbKE6mQO1HewvKuaXEX+C3Ctg7nqrS1LKr/hd2DsdMu7hme4BNzWtvVS0dFPW1M2Rpm5KG7p441AjXuM7ZlpkCIsz4liSEceymfEsyYgjMtTv/pgCzr2b9vLbsHsIcgis/2+ry1HK79gqxaJCg8ifEUP+jJgP7O8b9HCooZODdZ28e7SdfUfb+XtpE8b4/nGZPyOG5VkJrJo1jRWzEogJC7aoBfZUWFbPjVXfIj+oCrn2cYifaXVJSvkdMcZY8sEFBQVmz549lnz2aHT2u9hb086eqjZ2Vbax92g7g24vDoEFabGcPzuR83MTWTYzntAgHfufNB43u390Jcv7tzFw1a8JPed6qytSylIiUmiMGfNdejTsR6nf5WFvTTvbK1rZVtbC3qPteLyG8GAnq3KmccGcJC6Yk0RWYqTVpQaO3jZan7iFaUdfZcfcb3Luhm9ZXZFSltOwP8u6+l3srGhj65Fmth5upqq1F4DsxEjWzU3morxkVmQnEBKklzKMy+FXMM/fhru7lXudN7Lxrp/o7CmlGH/Y22rMfiJFhwVzSX4Kl+SnAFDV0sOWw828caiJx3ZW89DblUSFBrF2TiIX56WwLi+ZhEi9rP+0Govh7Z/B/idpjZjNjQN3cOcNH9OgV+oMadhPkKzESLISI/nM6ix6B91sK2vl9UONvF7SxEsHGnAIFGQlcOm8FD40P4WZ03S4530D3VDxd9j1W6jcAkHhtCy5lbW7VnLp4kw+NH+61RUq5fd0GGeSeb2GA8c6eK2kkVeLGznU0AXA3JRoLs33Bf/CtFj7zO33eqG9GpoPQcMBqNgMR3f5brSsik4AAAk6SURBVEASkwYrbsa1+EauebiExs5+Xr3zAuL1f0RKvW9Sx+xFZD3wM8AJPGCM+cEJz18NfBfwAm7gDmPMW6d6z3GHvccFrt6xv26KqD3ex99Lm3jjUBPv1LTj9RpSYkJZl5fMRXNTKMiKI9g5znH+EX+W5tTHvL/PDH0/9NV4h33v8W173/vq9j08g76fiXsA3P2+n42rDwa6oK8d+jugtxW66qCzHjqPffDnl7oYZq2DnHUw8zxwBvPz14/w01cPc98Ny1i/QHv1Sg03aWEvIk7gMHApUAvsBjYYY4qHHRMF9BhjjIgsAp4yxuSd6n3HHfZFf4anPzv21ylrBEdCRAJEp0JMKsSk++4qlZwPSXN9Nwof5tXiRr70+z18ZNEMfr5hqUVFKzV1TeYJ2hVAmTGmYuiDngSuBt4Pe2NM97DjI/lAV3KCpSyEy/5r0t7eKoMeLxXNPZQ2dHGkqYu+QQ9Oh5CdGMnc6dHMSYkmalRX8o4wHPSBYSI5yX754LYIiMO3XwTECQ6nb58jaOh7JzhDfOvJO0N8a9UER/i+hkT5gtw5+gvQ3j3azleeeIcFabH84OMLR/06pdTpjSY90oCjw7ZrgZUnHiQiHwX+G0gGrpiQ6k4mcbbvEWBCgLyhh9vjZXfVcV4tbuTe4gZq9/YhAksz4rgkP4VL56UwOzkqoMb5a1p7+cKju0mKDuXBzyzXFUuVmmCjGca5FrjMGPPFoe0bgRXGmK+McPxa4N+MMZec5LlbgFsAMjMzl1VXV59h+YHPGMOhhi5eLfad4D1wrAOAmdMiuDgvhUvmJVOQ5d/z+Wtae/nMw7s43jvIn768mpykKKtLUmrKmswx+1XAd4wxlw1t3w1gjBlxNSoRqQSWG2NaRjrGLrNxJlp9Rx+vlzTxekkjb5e3Muj2Eh0axJo5iaybm8wFc5NIjg6zusxR21XZxsbHCvF4DQ99toBlMxOsLkmpKW0ywz4I3wnai4Fj+E7QXm+MKRp2zGygfOgE7TnAX4F0c4o317A/cz0Dbt4ua3l/dk9j5wAA82fEcMGcJNbOSeKczPgp2+t/es9RvvXnA2TER/DgZ5eTrUtNKHVak3aC1hjjFpHbgE34pl4+ZIwpEpGNQ8/fB3wcuElEXEAf8KlTBb2aGJGhQXxo/nQ+NH86xhhK6rvYfLiJzaXN/GZrBb/aXE5EiJNzZ01jdc40Vuckkjc9GofD2rH+mtZevvtiMa8WN3L+7ETuvf4cYiN0JVGlJpNeVBWgOvtd7Chv5c0jLbx55P/W7kmIDGFFVgLLsxNYmZ1A3vRogsY7r3+MOvpcPPhWJfdtKSfIIXzloly+uCZ7/NcVKGVDujaO+oCYsOD3e/0Ax9r72F7eyrbyFnZXtfFyUQMA4cFOFqbHsjQjjkXpccyfEUNmQsSE9f6NMRRWH+eJXUd58UAd/S4vVy+Zwd2Xz2N6rP+cW1DK32nP3qbqO/p86/TXtLP3aDvFdR24PL7fhcgQJ3OnRzMrKYpZSZFkT4skNS6c6TFhJEWH4hzhHwK3x0tj1wB17X0UHetgV1UbuyqP09I9QFRoEFctmcH1KzJZkBZ70tcrpU5PlzhWZ6Tf5eFIYzfF9R0U13VyqKGLipYemrsGPnCcCESGBBEZ6iQiJAivMbjcXgY9Xo73uvB4/+/3KS0unJXZCayencjlC6br7R2VmgA6jKPOSNjQcM7C9A/2urv6XVS39tLQ0U9DZz9Nnf10D3joGXDTM+jG6RCCnQ6CnQ6mRYaQFh/OjLhwZidHkRYXblFrlFIn0rBXpxQdFsyCtFgdelHKz+k0CKWUsgENe6WUsgENe6WUsgENe6WUsgENe6WUsgENe6WUsgENe6WUsgENe6WUsgHLlksQkWZgvLeqSgRGvDGKDdi5/XZuO9i7/dp2n5nGmKSxvoFlYX8mRGTPeNaGCBR2br+d2w72br+2/czarsM4SillAxr2SillA/4a9vdbXYDF7Nx+O7cd7N1+bfsZ8Msxe6WUUmPjrz17pZRSY6Bhr5RSNuB3YS8i60WkVETKROQuq+uZTCKSISJ/F5ESESkSkduH9ieIyKsicmToa7zVtU4WEXGKyF4ReWFo205tjxORZ0Tk0NDvwCq7tF9E7hz6nT8oIk+ISFggt11EHhKRJhE5OGzfiO0VkbuHMrBURC4bzWf4VdiLiBO4F7gcyAc2iEi+tVVNKjfwNWPMPOBc4Nah9t4FvG6MyQVeH9oOVLcDJcO27dT2nwEvG2PygMX4/hwCvv0ikgZ8FSgwxiwAnMB1BHbbHwHWn7DvpO0dyoDrgPlDr/nVUDaekl+FPbACKDPGVBhjBoEngastrmnSGGPqjTHvDH3fhe8vexq+Nj86dNijwDXWVDi5RCQduAJ4YNhuu7Q9BlgLPAhgjBk0xrRjk/bju2VquIgEARFAHQHcdmPMVqDthN0jtfdq4EljzIAxphIow5eNp+RvYZ8GHB22XTu0L+CJSBawFNgJpBhj6sH3DwKQbF1lk+oe4BuAd9g+u7R9FtAMPDw0jPWAiERig/YbY44BPwZqgHqgwxjzCjZo+wlGau+4ctDfwl5Osi/g546KSBTwJ+AOY0yn1fWcDSLyEaDJGFNodS0WCQLOAX5tjFkK9BBYwxYjGhqbvhrIBmYAkSJyg7VVTSnjykF/C/taIGPYdjq+/94FLBEJxhf0jxtjnh3a3SgiqUPPpwJNVtU3ic4DrhKRKnzDdReJyGPYo+3g+12vNcbsHNp+Bl/426H9lwCVxphmY4wLeBZYjT3aPtxI7R1XDvpb2O8GckUkW0RC8J2keN7imiaNiAi+MdsSY8xPhz31PPCZoe8/Azx3tmubbMaYu40x6caYLHw/5zeMMTdgg7YDGGMagKMiMndo18VAMfZofw1wrohEDP0duBjf+So7tH24kdr7PHCdiISKSDaQC+w67bsZY/zqAXwYOAyUA/9qdT2T3Nbz8f33bD+wb+jxYWAavrPzR4a+Jlhd6yT/OVwIvDD0vW3aDiwB9gz9/P8CxNul/cB/AIeAg8DvgdBAbjvwBL7zEy58PfcvnKq9wL8OZWApcPloPkOXS1BKKRvwt2EcpZRS46Bhr5RSNqBhr5RSNqBhr5RSNqBhr5RSNqBhr5RSNqBhr5RSNvD/AX+Q9X6PKbJhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def tanh_fit(x, a, b, c, d):\n",
    "     return np.tanh((x-a)/b)*c+d\n",
    "xdata = np.linspace(0, x_mesh, x_mesh+1)\n",
    "popt, pcov = curve_fit(tanh_fit, xdata, u.compute_vertex_values()[0:x_mesh+1])\n",
    "plt.plot(xdata, u.compute_vertex_values()[0:x_mesh+1])\n",
    "a, b, c, d = popt\n",
    "plt.plot(xdata, tanh_fit(xdata, a, b, c, d))# np.tanh((xdata-52)/12.4)*0.36+0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x127c054a8>]"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAANR0lEQVR4nO3dQYic533H8e+vUgQNSWMTbUIqyZVa5CQ62MWZOKY0rdPQWnIPIuCD7RATExCmdsjRptDk4EtzKIRgO0IYYXKJDo1JlKLEFErigutWK7Bly8ZmKxNro4DXcUjBORjZ/x5mUqbr2Z135Xd3Nc9+P7Cw7/s+2vk/rPj69WhnJ1WFJGn2/d5mDyBJ6odBl6RGGHRJaoRBl6RGGHRJasT2zXrgnTt31t69ezfr4SVpJp05c+b1qpqbdG3Tgr53717m5+c36+ElaSYl+flK13zKRZIaYdAlqREGXZIaYdAlqREGXZIaMTXoSY4neS3J8ytcT5JvJ1lIcjbJDf2PKUmapssd+mPAwVWuHwL2jz6OAN9572NJktZqatCr6kngjVWWHAa+W0NPA1cl+VhfA0qSuunjOfRdwIWx48XRuXdJciTJfJL5paWlHh5akvQ7fQQ9E85NfNeMqjpWVYOqGszNTXzlqiTpMvUR9EVgz9jxbuBiD19XkrQGfQT9JHDX6KddbgJ+U1W/7OHrSpLWYOov50ryPeBmYGeSReAbwPsAquoocAq4FVgAfgvcvV7DSpJWNjXoVXXHlOsF3NvbRJKky+IrRSWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEZ2CnuRgkpeSLCR5YML1DyX5UZJnk5xLcnf/o0qSVjM16Em2AQ8Dh4ADwB1JDixbdi/wQlVdD9wM/FOSHT3PKklaRZc79BuBhao6X1VvASeAw8vWFPDBJAE+ALwBXOp1UknSqroEfRdwYex4cXRu3EPAJ4GLwHPA16rqneVfKMmRJPNJ5peWli5zZEnSJF2CngnnatnxLcAzwB8Cfwo8lOQP3vWHqo5V1aCqBnNzc2seVpK0si5BXwT2jB3vZngnPu5u4PEaWgBeAT7Rz4iSpC66BP00sD/JvtE/dN4OnFy25lXg8wBJPgp8HDjf56CSpNVtn7agqi4luQ94AtgGHK+qc0nuGV0/CjwIPJbkOYZP0dxfVa+v49ySpGWmBh2gqk4Bp5adOzr2+UXgb/odTZK0Fr5SVJIaYdAlqREGXZIaYdAlqREGXZIaYdAlqREGXZIaYdAlqREGXZIaYdAlqREGXZIaYdAlqREGXZIaYdAlqREGXZIaYdAlqREGXZIaYdAlqREGXZIaYdAlqREGXZIaYdAlqREGXZIaYdAlqREGXZIaYdAlqREGXZIaYdAlqREGXZIaYdAlqRGdgp7kYJKXkiwkeWCFNTcneSbJuSQ/63dMSdI026ctSLINeBj4a2AROJ3kZFW9MLbmKuAR4GBVvZrkI+s1sCRpsi536DcCC1V1vqreAk4Ah5etuRN4vKpeBaiq1/odU5I0TZeg7wIujB0vjs6Nuxa4OslPk5xJctekL5TkSJL5JPNLS0uXN7EkaaIuQc+Ec7XseDvwKeBvgVuAf0hy7bv+UNWxqhpU1WBubm7Nw0qSVjb1OXSGd+R7xo53AxcnrHm9qt4E3kzyJHA98HIvU0qSpupyh34a2J9kX5IdwO3AyWVrfgh8Nsn2JO8HPgO82O+okqTVTL1Dr6pLSe4DngC2Acer6lySe0bXj1bVi0l+ApwF3gEerarn13NwSdL/l6rlT4dvjMFgUPPz85vy2JI0q5KcqarBpGu+UlSSGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGtEp6EkOJnkpyUKSB1ZZ9+kkbye5rb8RJUldTA16km3Aw8Ah4ABwR5IDK6z7JvBE30NKkqbrcod+I7BQVeer6i3gBHB4wrqvAt8HXutxPklSR12Cvgu4MHa8ODr3f5LsAr4AHF3tCyU5kmQ+yfzS0tJaZ5UkraJL0DPhXC07/hZwf1W9vdoXqqpjVTWoqsHc3FzXGSVJHWzvsGYR2DN2vBu4uGzNADiRBGAncGuSS1X1g16mlCRN1SXop4H9SfYBvwBuB+4cX1BV+373eZLHgH8x5pK0saYGvaouJbmP4U+vbAOOV9W5JPeMrq/6vLkkaWN0uUOnqk4Bp5admxjyqvryex9LkrRWvlJUkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhph0CWpEZ2CnuRgkpeSLCR5YML1LyY5O/p4Ksn1/Y8qSVrN1KAn2QY8DBwCDgB3JDmwbNkrwF9W1XXAg8CxvgeVJK2uyx36jcBCVZ2vqreAE8Dh8QVV9VRV/Xp0+DSwu98xJUnTdAn6LuDC2PHi6NxKvgL8eNKFJEeSzCeZX1pa6j6lJGmqLkHPhHM1cWHyOYZBv3/S9ao6VlWDqhrMzc11n1KSNNX2DmsWgT1jx7uBi8sXJbkOeBQ4VFW/6mc8SVJXXe7QTwP7k+xLsgO4HTg5viDJNcDjwJeq6uX+x5QkTTP1Dr2qLiW5D3gC2AYcr6pzSe4ZXT8KfB34MPBIEoBLVTVYv7ElSculauLT4etuMBjU/Pz8pjy2JM2qJGdWumH2laKS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1IhOQU9yMMlLSRaSPDDhepJ8e3T9bJIb+h9VkrSaqUFPsg14GDgEHADuSHJg2bJDwP7RxxHgOz3PKUmaossd+o3AQlWdr6q3gBPA4WVrDgPfraGngauSfKznWSVJq+gS9F3AhbHjxdG5ta4hyZEk80nml5aW1jqrJGkVXYKeCefqMtZQVceqalBVg7m5uS7zSZI66hL0RWDP2PFu4OJlrJEkraMuQT8N7E+yL8kO4Hbg5LI1J4G7Rj/tchPwm6r6Zc+zSpJWsX3agqq6lOQ+4AlgG3C8qs4luWd0/ShwCrgVWAB+C9y9fiNLkiaZGnSAqjrFMNrj546OfV7Avf2OJklaC18pKkmNMOiS1AiDLkmNMOiS1IgM/z1zEx44WQJ+fpl/fCfweo/jzAL3vDW4563hvez5j6pq4iszNy3o70WS+aoabPYcG8k9bw3ueWtYrz37lIskNcKgS1IjZjXoxzZ7gE3gnrcG97w1rMueZ/I5dEnSu83qHbokaRmDLkmNuKKDvhXfnLrDnr842uvZJE8luX4z5uzTtD2Prft0kreT3LaR862HLntOcnOSZ5KcS/KzjZ6xbx3+bn8oyY+SPDva80z/1tYkx5O8luT5Fa7336+quiI/GP6q3v8G/hjYATwLHFi25lbgxwzfMekm4D83e+4N2POfAVePPj+0FfY8tu7fGP7Wz9s2e+4N+D5fBbwAXDM6/shmz70Be/574Jujz+eAN4Admz37e9jzXwA3AM+vcL33fl3Jd+hb8c2pp+65qp6qql+PDp9m+O5Qs6zL9xngq8D3gdc2crh10mXPdwKPV9WrAFU16/vusucCPpgkwAcYBv3Sxo7Zn6p6kuEeVtJ7v67koPf25tQzZK37+QrD/8LPsql7TrIL+AJwlDZ0+T5fC1yd5KdJziS5a8OmWx9d9vwQ8EmGb1/5HPC1qnpnY8bbFL33q9MbXGyS3t6ceoZ03k+SzzEM+p+v60Trr8uevwXcX1VvD2/eZl6XPW8HPgV8Hvh94D+SPF1VL6/3cOuky55vAZ4B/gr4E+Bfk/x7Vf3Peg+3SXrv15Uc9K345tSd9pPkOuBR4FBV/WqDZlsvXfY8AE6MYr4TuDXJpar6wcaM2Luuf7dfr6o3gTeTPAlcD8xq0Lvs+W7gH2v4BPNCkleATwD/tTEjbrje+3UlP+WyFd+ceuqek1wDPA58aYbv1sZN3XNV7auqvVW1F/hn4O9mOObQ7e/2D4HPJtme5P3AZ4AXN3jOPnXZ86sM/4+EJB8FPg6c39ApN1bv/bpi79BrC745dcc9fx34MPDI6I71Us3wb6rruOemdNlzVb2Y5CfAWeAd4NGqmvjjb7Og4/f5QeCxJM8xfDri/qqa2V+rm+R7wM3AziSLwDeA98H69cuX/ktSI67kp1wkSWtg0CWpEQZdkhph0CWpEQZdkhph0CWpEQZdkhrxv0JmifRqw5HQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 50\n",
    "plt.plot(u.compute_vertex_values()[n*(x_mesh+1):(n+1)*(x_mesh+1)])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}