Skip to content
Snippets Groups Projects
Half_Frap_Phi_u.ipynb 3.2 KiB
Newer Older
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#### READ ME ###\n",
    "\n",
    "################ When load this in Paraview !!!!!!!!!!!! => Work the first time then bug !!!!!!!!!!!\n",
    "\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",
    "dt = 0.01\n",
    "Radius_Droplet = 5\n",
    "\n",
    "\n",
    "mesh = generate_mesh(Circle(Point(0,0),20.0),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_notbleach = 0.35\n",
    "\n",
    "Phi_u0 = interpolate(Expression('sqrt(x[1]*x[1]+x[0]*x[0]) <= Radius_Droplet && atan2(x[1],x[0]) < 0  ? Phi_u0_int : sqrt(x[1]*x[1]+x[0]*x[0]) <= Radius_Droplet && atan2(x[1],x[0]) >= 0 ?  Phi_u0_notbleach : Phi_u0_ext' , degree=2,tol=tol,Phi_u0_int = Phi_u0_int,Phi_u0_ext=Phi_u0_ext, Phi_u0_notbleach = Phi_u0_notbleach, Radius_Droplet=Radius_Droplet),V)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\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_Droplet ? Phi_tot_int : Phi_tot_ext', degree=2,tol=tol,Phi_tot_int=Phi_tot_int,Phi_tot_ext=Phi_tot_ext,Radius_Droplet=Radius_Droplet),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",
    "\n",
    "vtkfile = File('Disk/Phi0.pvd')\n",
    "vtkfile << Phi_u0\n",
    "\n",
    "\n",
    "\n",
    "t = 0\n",
    "cFile = XDMFFile('test.xmdf')\n",
    "cFile.write(Phi_u0, t)\n",
    "\n",
    "ti = time.time()\n",
    "for i in range(200):\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"
   ]
  },
  {
   "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
}