{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.384448766708374\n",
      "2.8050124645233154\n",
      "4.201659917831421\n",
      "5.629311561584473\n",
      "7.02629828453064\n",
      "8.605925798416138\n",
      "10.40172266960144\n",
      "12.074671268463135\n",
      "13.842711925506592\n",
      "15.555825471878052\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "from __future__ import print_function\n",
    "from fenics import *\n",
    "from mshr import *\n",
    "from dolfin import *\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "from math import *\n",
    "\n",
    "set_log_level(40)\n",
    "\n",
    "\n",
    "def solver_disk(Radius_Drop, Radius_Disk):\n",
    "    \n",
    "    dt = 0.5\n",
    "    mesh = generate_mesh(Circle(Point(0,0),Radius_Disk),120)\n",
    "    V = FunctionSpace(mesh, 'CG', 1)\n",
    "\n",
    "    \n",
    "    Phi_u = Function(V)\n",
    "    v = TestFunction(V)\n",
    "\n",
    "\n",
    "    Phi_u0 = Function(V)\n",
    "    Phi_tot = Function(V)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "    Phi_tot = interpolate(Expression('0.5-0.4*tanh(sqrt(x[0]*x[0]+x[1]*x[1]) - 5) ', degree=2), V)\n",
    "\n",
    "    Phi_u0 = interpolate(Expression('0.0505 - (0.0495)*tanh(-sqrt(x[0]*x[0]+x[1]*x[1]) + 5) ', degree=2), V)\n",
    "\n",
    "    \n",
    "    def boundary(x, on_boundary):\n",
    "        tol = 1E-14\n",
    "        return on_boundary\n",
    "\n",
    "    u_D = Constant(0.1)\n",
    "\n",
    "    bcC = DirichletBC(V, u_D, boundary)\n",
    "\n",
    "\n",
    "    Form = (inner((Phi_u-Phi_u0)/dt, v) -  inner(((1-Phi_tot)/Phi_tot)*Phi_u*grad(Phi_tot) - (1-Phi_tot)*grad(Phi_u), grad(v))) *dx\n",
    "\n",
    "    \n",
    "    \n",
    "\n",
    "\n",
    "    t = 0\n",
    "    cFile = XDMFFile('Phi_u_Disk.xmdf')\n",
    "    cFile.write(Phi_u0, t)\n",
    "\n",
    "    ti = time.time()\n",
    "    for i in range(100):\n",
    "        solve(Form == 0, Phi_u, bcC)\n",
    "        assign(Phi_u0, Phi_u)\n",
    "        t += dt\n",
    "        cFile.write(Phi_u0, t)\n",
    "        print(time.time() - ti)\n",
    "\n",
    "    cFile.close()\n",
    "    \n",
    "    return 0\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def solver_sphere(Radius_Drop, Radius_Sphere):\n",
    "\n",
    "\n",
    "    dt = 0.15\n",
    "\n",
    "\n",
    "   # mesh = generate_mesh(Sphere(Point(0,0,0),Radius_Sphere),30)\n",
    "    mesh = Mesh('Meshes_Sphere.xml')    \n",
    "    P1 = FiniteElement(\"Lagrange\", mesh.ufl_cell(), 1)\n",
    "    V = FunctionSpace(mesh, P1)\n",
    "\n",
    "    # V = FunctionSpace(mesh, 'CG', 1)\n",
    "    Phi_u = Function(V)\n",
    "    v = TestFunction(V)\n",
    "\n",
    "\n",
    "    Phi_u0 = Function(V)\n",
    "    Phi_tot = Function(V)\n",
    "\n",
    "    tol = 1E-14\n",
    "\n",
    "\n",
    "    Phi_tot = interpolate(Expression('0.5-0.4*tanh(sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]) - Radius_Drop) ', Radius_Drop = Radius_Drop, degree=2), V)\n",
    "\n",
    "    Phi_u0 = interpolate(Expression('0.0505 - (0.0495)*tanh(-sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]) + Radius_Drop) ',Radius_Drop = Radius_Drop , degree=2), V)\n",
    "\n",
    "    u_D = Constant(0.01)\n",
    "\n",
    "\n",
    "    def boundary(x, on_boundary):\n",
    "        tol = 1E-14\n",
    "        return on_boundary\n",
    "\n",
    "\n",
    "\n",
    "    bc = DirichletBC(V, u_D, boundary)\n",
    "\n",
    "\n",
    "    Form = (inner((Phi_u-Phi_u0)/dt, v) -  inner(((1-Phi_tot)/Phi_tot)*Phi_u*grad(Phi_tot) - (1-Phi_tot)*grad(Phi_u), grad(v))) *dx\n",
    "\n",
    "\n",
    "\n",
    "    t = 0\n",
    "    cFile = XDMFFile('Phi_u_Sphere.xmdf')\n",
    "    cFile.write(Phi_u0, t)\n",
    "\n",
    "    ti = time.time()\n",
    "    for i in range(10):\n",
    "        #solve(Form == 0, Phi_u, bcC)\n",
    "        solve(Form == 0, Phi_u, bc,solver_parameters={'newton_solver': {'linear_solver': 'gmres'}})\n",
    "        assign(Phi_u0, Phi_u)\n",
    "        t += dt\n",
    "        cFile.write(Phi_u0, t)\n",
    "        print(time.time() - ti)\n",
    "\n",
    "    cFile.close()\n",
    "\n",
    "    return 0\n",
    "\n",
    "\n",
    "\n",
    "def solver_half_frap_disk(Radius_Drop, Radius_Disk):\n",
    "\n",
    "    dt = 0.01\n",
    "    mesh = generate_mesh(Circle(Point(0,0),Radius_Disk),100)\n",
    "\n",
    "    V = FunctionSpace(mesh, 'CG', 1)\n",
    "    Phi_u = Function(V)\n",
    "    v = TestFunction(V)\n",
    "\n",
    "\n",
    "    Phi_u0 = Function(V)\n",
    "    Phi_tot = Function(V)\n",
    "\n",
    "    tol = 1E-14\n",
    "\n",
    "\n",
    "    Phi_u0_int = 0.01\n",
    "    Phi_u0_ext = 0.05\n",
    "    Phi_u0_int_notbleach = 0.35\n",
    "\n",
    "    Phi_u0 = interpolate(Expression('sqrt(x[1]*x[1]+x[0]*x[0]) <= Radius_Drop && atan2(x[1],x[0]) < 0  ? Phi_u0_int : sqrt(x[1]*x[1]+x[0]*x[0]) <= Radius_Drop && atan2(x[1],x[0]) >= 0 ?  Phi_u0_int_notbleach : Phi_u0_ext' , degree=2,tol=tol,Phi_u0_int = Phi_u0_int,Phi_u0_ext=Phi_u0_ext, Phi_u0_int_notbleach = Phi_u0_int_notbleach, Radius_Drop=Radius_Drop),V)\n",
    "\n",
    "    \n",
    "    Phi_tot_int = 0.35\n",
    "    Phi_tot_ext = 0.05\n",
    "\n",
    "    Phi_tot = interpolate(Expression('sqrt(x[1]*x[1]+x[0]*x[0]) <= Radius_Drop ? Phi_tot_int : Phi_tot_ext', degree=2,tol=tol,Phi_tot_int=Phi_tot_int,Phi_tot_ext=Phi_tot_ext,Radius_Drop=Radius_Drop),V)\n",
    "\n",
    "\n",
    "    u_D = Constant(0.05)\n",
    "\n",
    "    def boundary(x, on_boundary):\n",
    "        tol = 1E-14\n",
    "        return on_boundary\n",
    "\n",
    "    bc = DirichletBC(V, u_D, boundary)\n",
    "\n",
    "    Form = (inner((Phi_u-Phi_u0)/dt, v) - inner(((1-Phi_tot)/Phi_tot)*Phi_u*grad(Phi_tot) - (1-Phi_tot)*grad(Phi_u), grad(v))) *dx\n",
    "\n",
    "\n",
    "    t = 0\n",
    "    cFile = XDMFFile('Phi_u_half_frap_disk.xmdf')\n",
    "    cFile.write(Phi_u0, t)\n",
    "\n",
    "    ti = time.time()\n",
    "    for i in range(50):\n",
    "        solve(Form == 0, Phi_u, bc)\n",
    "        assign(Phi_u0, Phi_u)\n",
    "        t += dt\n",
    "        cFile.write(Phi_u0, t)\n",
    "        print(time.time() - ti)\n",
    "\n",
    "    cFile.close()\n",
    "\n",
    "    \n",
    "    \n",
    "    \n",
    "def solver_half_frap_sphere(Radius_Drop, Radius_Sphere):\n",
    "\n",
    "    dt = 0.01\n",
    "    #mesh = generate_mesh(Sphere(Point(0,0,0),Radius_Sphere),50)\n",
    "    mesh = Mesh('Meshes_Sphere.xml')    \n",
    "    P1 = FiniteElement(\"Lagrange\", mesh.ufl_cell(), 1)\n",
    "    V = FunctionSpace(mesh, P1)\n",
    "    Phi_u = Function(V)\n",
    "    v = TestFunction(V)\n",
    "\n",
    "\n",
    "    Phi_u0 = Function(V)\n",
    "    Phi_tot = Function(V)\n",
    "\n",
    "    tol = 1E-14\n",
    "\n",
    "\n",
    "    Phi_u0_int = 0.01\n",
    "    Phi_u0_ext = 0.05\n",
    "    Phi_u0_int_notbleach = 0.35\n",
    "\n",
    "    Phi_u0 = interpolate(Expression('sqrt(x[1]*x[1]+x[0]*x[0]+x[2]*x[2]) <= Radius_Drop && atan2(x[1],x[0]) < 0  ? Phi_u0_int : sqrt(x[1]*x[1]+x[0]*x[0]+x[2]*x[2]) <= Radius_Drop && atan2(x[1],x[0]) >= 0 ?  Phi_u0_int_notbleach : Phi_u0_ext' , degree=2,tol=tol,Phi_u0_int = Phi_u0_int,Phi_u0_ext=Phi_u0_ext, Phi_u0_int_notbleach = Phi_u0_int_notbleach, Radius_Drop=Radius_Drop),V)\n",
    "\n",
    "    \n",
    "    Phi_tot_int = 0.35\n",
    "    Phi_tot_ext = 0.05\n",
    "\n",
    "    Phi_tot = interpolate(Expression('sqrt(x[1]*x[1]+x[0]*x[0]+x[2]*x[2]) <= Radius_Drop ? Phi_tot_int : Phi_tot_ext', degree=2,tol=tol,Phi_tot_int=Phi_tot_int,Phi_tot_ext=Phi_tot_ext,Radius_Drop=Radius_Drop),V)\n",
    "\n",
    "\n",
    "    u_D = Constant(0.05)\n",
    "\n",
    "    def boundary(x, on_boundary):\n",
    "        tol = 1E-14\n",
    "        return on_boundary\n",
    "\n",
    "    bc = DirichletBC(V, u_D, boundary)\n",
    "\n",
    "    Form = (inner((Phi_u-Phi_u0)/dt, v) - inner(((1-Phi_tot)/Phi_tot)*Phi_u*grad(Phi_tot) - (1-Phi_tot)*grad(Phi_u), grad(v))) *dx\n",
    "\n",
    "\n",
    "    t = 0\n",
    "    cFile = XDMFFile('Phi_u_half_frap_sphere.xmdf')\n",
    "    cFile.write(Phi_u0, t)\n",
    "\n",
    "    ti = time.time()\n",
    "    for i in range(50):\n",
    "        #solve(Form == 0, Phi_u, bc)\n",
    "        solve(Form == 0, Phi_u, bc,solver_parameters={'newton_solver': {'linear_solver': 'gmres'}})\n",
    "\n",
    "        assign(Phi_u0, Phi_u)\n",
    "        t += dt\n",
    "        cFile.write(Phi_u0, t)\n",
    "        print(time.time() - ti)\n",
    "\n",
    "    cFile.close()\n",
    "    \n",
    "    return 0\n",
    "    \n",
    "\n",
    "#solver_disk(3,10)\n",
    "\n",
    "#solver_sphere(0.2,1)\n",
    "\n",
    "#solver_half_frap_disk(2,10)\n",
    "\n",
    "solver_half_frap_sphere(0.2,1)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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
}