{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "from fenics import *\n",
    "from dolfin import *\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import time\n",
    "from math import *\n",
    "\n",
    "set_log_level(40)\n",
    "\n",
    "def point_expr(c, p_list, Phi_tot_int, Phi_tot_ext, Phi_tot_init,\n",
    "               rad_drop, tol):\n",
    "    s = ('sqrt(pow((x[0]-'+str(c[0])+'),2)+pow((x[1]-'+str(c[1])+'),2)+'+\n",
    "            'pow((x[2]-'+str(c[2])+'),2))')\n",
    "    s = s + '<= '+str(rad_drop)+' ? '+str(Phi_tot_init)+' : '\n",
    "    for p in p_list:\n",
    "        y = ('sqrt(pow((x[0]-'+str(p[0])+'),2)+pow((x[1]-'+str(p[1])+'),2)+'+\n",
    "            'pow((x[2]-'+str(p[2])+'),2))')\n",
    "        s = s + y + '<= '+str(rad_drop)+' ? '+str(Phi_tot_int)+' : '\n",
    "    s = s + str(Phi_tot_ext)\n",
    "    return Expression(s, degree=2, tol=tol)\n",
    "\n",
    "tol = 1E-14\n",
    "dt = 0.1\n",
    "rad_drop = 0.25\n",
    "\n",
    "Phi_tot_int = 0.99\n",
    "Phi_tot_ext = 0.001\n",
    "Phi_tot_initial = 0\n",
    "\n",
    "m = ['Meshes/multi_drop_gauss.xml', 'Meshes/multi_drop_gauss_med.xml',\n",
    "     'Meshes/multi_drop_gauss_far.xml']\n",
    "point_lists = [[[4, 4.5, 0.5], [4, 3.5, 0.5], [3.5, 4, 0.5], [4.5, 4, 0.5]],\n",
    "               [[4, 5, 0.5], [4, 3, 0.5], [3, 4, 0.5], [5, 4, 0.5]],\n",
    "               [[4, 5.5, 0.5], [4, 2.5, 0.5], [2.5, 4, 0.5], [5.5, 4, 0.5]]]\n",
    "cs = []\n",
    "for j in range(len(m)):\n",
    "    mesh = Mesh(m[j])\n",
    "    V = FunctionSpace(mesh, 'CG', 1)\n",
    "    Phi_u = Function(V)\n",
    "    v = TestFunction(V)\n",
    "\n",
    "    Phi_u0 = Function(V)\n",
    "    Phi_tot = Function(V)\n",
    "\n",
    "    Phi_u0 = point_expr([4,4,0.5], point_lists[j],\n",
    "                        Phi_tot_int, Phi_tot_ext, Phi_tot_initial, rad_drop,\n",
    "                        tol)\n",
    "    Phi_u0 = interpolate(Phi_u0, V)\n",
    "\n",
    "    Phi_tot = point_expr([4,4,0.5], point_lists[j],\n",
    "                        Phi_tot_int, Phi_tot_ext, Phi_tot_int, rad_drop, tol)\n",
    "    Phi_tot = interpolate(Phi_tot,V)\n",
    "\n",
    "    Form = (inner((Phi_u-Phi_u0)/dt, v) -\n",
    "            inner(((1-Phi_tot)/Phi_tot)*Phi_u*grad(Phi_tot)- \n",
    "            (1-Phi_tot)*grad(Phi_u), grad(v)))*dx\n",
    "\n",
    "    t = 0\n",
    "    cFile = XDMFFile('4_neighbors_'+str(j)+'.xdmf')\n",
    "    cFile.write(Phi_u0, t)\n",
    "    \n",
    "    cs.append([])\n",
    "    ti = time.time()\n",
    "    for i in range(100):\n",
    "        cs[j].append([Phi_u([x, 4, 0.5]) for x in np.linspace(4, 4.49, 100)])\n",
    "        solve(Form == 0, Phi_u)\n",
    "        assign(Phi_u0, Phi_u)\n",
    "        t += dt\n",
    "        cFile.write(Phi_u0, t)\n",
    "        print(time.time() - ti)\n",
    "\n",
    "    cFile.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "font = {'family' : 'normal',\n",
    "        'weight' : 'normal',\n",
    "        'size'   : 12}\n",
    "import matplotlib\n",
    "matplotlib.rc('font', **font)\n",
    "\n",
    "def make_graph_pretty(xlab, ylab, loc=0, markerfirst=True, handlelength=None):\n",
    "    color1 = 'black'#'darkorange'\n",
    "    color2 = (0, 0, 0)\n",
    "    color2 = (1, 1, 1)\n",
    "    ax = plt.gca()\n",
    "    ax.tick_params(axis='x', colors=color1, which='both')\n",
    "    ax.tick_params(axis='y', colors=color1, which='both')\n",
    "    plt.xlabel(xlab, color=color1)\n",
    "    plt.ylabel(ylab, color=color1)\n",
    "    ax.spines['left'].set_color(color1)\n",
    "    ax.spines['bottom'].set_color(color1)\n",
    "    # Hide the right and top spines\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.set_facecolor(color2)\n",
    "    plt.gcf().subplots_adjust(bottom=0.17, left=0.15)\n",
    "    plt.legend(frameon=False, loc=loc, labelspacing=0.1,\n",
    "               markerfirst=markerfirst, handlelength=handlelength)\n",
    "\n",
    "for i in range(len(cs[0])):\n",
    "    plt.plot([np.mean(x[1:49])/0.8 for x in cs[0][0:i]], label='dist=0.5')\n",
    "    plt.plot([np.mean(x[1:49])/0.8 for x in cs[1][0:i]], label='dist=1')\n",
    "    plt.plot([np.mean(x[1:49])/0.8 for x in cs[2][0:i]], label='dist=1.5')\n",
    "    plt.plot(range(0, 100), np.ones(100), linestyle='--', color='k')\n",
    "    plt.ylim((0, 1.1))\n",
    "    plt.xlim((0, 100))\n",
    "    make_graph_pretty('t (a.u.)', 'intensity (a.u)', loc=4)\n",
    "    plt.savefig('img_'+str(i)+'.png', dpi=300)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}