{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import dolfin as df\n",
    "import matplotlib.pyplot as plt\n",
    "import mshr as ms\n",
    "import numpy as np\n",
    "import time\n",
    "df.set_log_level(40)\n",
    "\n",
    "domain = ms.Sphere(df.Point(0, 0, 0), 1.0)\n",
    "# mesh = ms.generate_mesh(domain, 35)\n",
    "mesh = df.UnitSquareMesh(50,50)\n",
    "dt = 0.01\n",
    "\n",
    "F = df.FunctionSpace(mesh, 'CG', 1)\n",
    "\n",
    "c0 = df.Function(F)\n",
    "c_tot = df.Function(F)\n",
    "c = df.Function(F)\n",
    "tc = df.TestFunction(F)\n",
    "\n",
    "# Interpolate c_tot and initial conditions\n",
    "# c_tot.interpolate(df.Expression('-0.4*tanh((sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))/0.1-0.5)) ', degree=1))\n",
    "c_tot.interpolate(df.Expression('0.4*tanh(-35*(sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))-0.2)) +0.5', degree=1))\n",
    "# tanh from left to right initial conditions\n",
    "# c_tot.interpolate(df.Expression('0.4*tanh(5*(x[0]-0.5)) - 1', degree=1))\n",
    "# Step function left to right initial condition\n",
    "# c0.interpolate(df.Expression('(x[0]<0.5) ? 0 : 1.0', degree=1))\n",
    "# Somewhat realistic half FRAP\n",
    "c0.interpolate(df.Expression(('(x[0]<0.5) && sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))<0.2 ? 0 :'\n",
    "                              '0.4*tanh(-35*(sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))-0.2)) +0.5'),\n",
    "                             degree=1))\n",
    "\n",
    "# Weak form\n",
    "form = ((df.inner((c-c0)/dt, tc) + (1-c_tot)* df.inner(df.grad(c), df.grad(tc))) -\n",
    "        (1-c_tot)/c_tot* c*df.inner(df.grad(c_tot), df.grad(tc))) * df.dx\n",
    "\n",
    "# Open file for writing\n",
    "cFile = df.XDMFFile('conc.xdmf')\n",
    "t = 0\n",
    "cFile.write(c0, t)\n",
    "\n",
    "# Solve in time\n",
    "ti = time.time()\n",
    "for i in range(80):\n",
    "    print(t)\n",
    "    df.solve(form == 0, c)\n",
    "    df.assign(c0, c)\n",
    "    t += dt\n",
    "    cFile.write(c0, t)\n",
    "cFile.close()\n",
    "print(time.time() - ti)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# c_tot.interpolate(df.Expression('-0.4*tanh((sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))/0.1-0.1)) ', degree=1))\n",
    "c_tot.interpolate(df.Expression('0.4*tanh(-35*(sqrt((x[0]-0.5)*(x[0]-0.5)+(x[1]-0.5)*(x[1]-0.5))-0.2)) +0.5', degree=1))\n",
    "plt.plot(np.linspace(0, 1, 100), [c_tot([x, 0.5]) for x in np.linspace(0, 1, 100)])"
   ]
  }
 ],
 "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}