Half_Frap_Phi_u.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "dt = 0.1\n",
    "Radius_Droplet = 0.25\n",
    "mesh = Mesh('Meshes/multi_drop_gauss.xml')\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",
    "tol = 1E-14\n",
    "\n",
    "Phi_tot_int = 0.99\n",
    "Phi_tot_ext = 0.001\n",
    "\n",
    "Phi_u0 = Expression('sqrt((x[1]-2)*(x[1]-2)+(x[0]-2)*(x[0]-2)+(x[2]-0.5)*(x[2]-0.5))'+\n",
    "                    '<= Radius_Droplet ? 0 : '+\n",
    "                    'sqrt((x[1]-2.5)*(x[1]-2.5)+(x[0]-2.5)*(x[0]-2.5)+(x[2]-0.5)*(x[2]-0.5))'+\n",
    "                    '<= Radius_Droplet ? Phi_tot_ext : Phi_tot_ext',\n",
    "                    Phi_tot_int=Phi_tot_int, degree=2, tol=tol, Phi_tot_ext=Phi_tot_ext,\n",
    "                    Radius_Droplet=Radius_Droplet)\n",
    "Phi_u0 = interpolate(Phi_u0, V)\n",
    "\n",
    "Phi_tot = Expression('sqrt((x[1]-2)*(x[1]-2)+(x[0]-2)*(x[0]-2)+(x[2]-0.5)*(x[2]-0.5))'+\n",
    "                     ' <= Radius_Droplet ? Phi_tot_int : '+\n",
    "                     'sqrt((x[1]-2.5)*(x[1]-2.5)+(x[0]-2.5)*(x[0]-2.5)+(x[2]-0.5)*(x[2]-0.5))'+\n",
    "                     '<= Radius_Droplet ? Phi_tot_ext: Phi_tot_ext', \n",
    "                     degree=2, tol=tol, Phi_tot_int=Phi_tot_int, Phi_tot_ext=Phi_tot_ext,\n",
    "                     Radius_Droplet=Radius_Droplet)\n",
    "Phi_tot = interpolate(Phi_tot,V)\n",
    "\n",
    "# Gamma_int = 1/2.5/10;\n",
    "# Gamma_ext = 1/10;\n",
    "Gamma0 = Expression('sqrt((x[1]-2)*(x[1]-2)+(x[0]-2)*(x[0]-2)+(x[2]-0.5)*(x[2]-0.5))'+\n",
    "                     ' <= Radius_Droplet ? Gamma_int : '+\n",
    "                     'sqrt((x[1]-2.5)*(x[1]-2.5)+(x[0]-2.5)*(x[0]-2.5)+(x[2]-0.5)*(x[2]-0.5))'+\n",
    "                     '<= Radius_Droplet ? Gamma_int: Gamma_ext', \n",
    "                     degree=2, tol=tol, Gamma_int=Gamma_int, Gamma_ext=Gamma_ext,\n",
    "                     Radius_Droplet=Radius_Droplet)\n",
    "Gamma0 = interpolate(Gamma0,V)\n",
    "\n",
    "Form = (inner((Phi_u-Phi_u0)/dt, v) - 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('multi_drop_no_neighbor.xdmf')\n",
    "cFile.write(Phi_u0, t)\n",
    "\n",
    "# c_in_b = []\n",
    "c_in_b1 = []\n",
    "ti = time.time()\n",
    "for i in range(100):\n",
    "    c_in_b1.append(Phi_u([2, 2, 0.5]))\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": [
    "plt.plot(c_in_b)\n",
    "plt.plot(c_in_b1)"
   ]
  },
  {
   "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
}