Newer
Older
{
"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",
Lars Hubatsch
committed
"def four_points(point, dist, rad):\n",
" e = Expression('sqrt((x[1]-4)*(x[1]-4)+(x[0]-4)*(x[0]-4)+(x[2]-0.5)*(x[2]-0.5))'+\n",
" '<= Radius_Droplet ? 0 : '+\n",
" 'sqrt((x[1]-5)*(x[1]-5)+(x[0]-4)*(x[0]-4)+(x[2]-0.5)*(x[2]-0.5))'+\n",
" '<= Radius_Droplet ? Phi_tot_int : '+\n",
" 'sqrt((x[1]-3)*(x[1]-3)+(x[0]-4)*(x[0]-4)+(x[2]-0.5)*(x[2]-0.5))'+\n",
" '<= Radius_Droplet ? Phi_tot_int : '+\n",
" 'sqrt((x[1]-4)*(x[1]-4)+(x[0]-5)*(x[0]-5)+(x[2]-0.5)*(x[2]-0.5))'+\n",
" '<= Radius_Droplet ? Phi_tot_int : '+\n",
" 'sqrt((x[1]-4)*(x[1]-4)+(x[0]-3)*(x[0]-3)+(x[2]-0.5)*(x[2]-0.5))'+\n",
" '<= Radius_Droplet ? Phi_tot_int : 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",
" return e\n",
"\n",
Lars Hubatsch
committed
"dt = 0.1\n",
"Radius_Droplet = 0.25\n",
Lars Hubatsch
committed
"mesh = Mesh('Meshes/multi_drop_gauss_med.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",
Lars Hubatsch
committed
"Phi_tot_int = 0.99\n",
"Phi_tot_ext = 0.001\n",
"\n",
Lars Hubatsch
committed
"Phi_u0 = Expression('sqrt((x[1]-4)*(x[1]-4)+(x[0]-4)*(x[0]-4)+(x[2]-0.5)*(x[2]-0.5))'+\n",
Lars Hubatsch
committed
" '<= Radius_Droplet ? 0 : '+\n",
Lars Hubatsch
committed
" 'sqrt((x[1]-5)*(x[1]-5)+(x[0]-4)*(x[0]-4)+(x[2]-0.5)*(x[2]-0.5))'+\n",
" '<= Radius_Droplet ? Phi_tot_int : '+\n",
" 'sqrt((x[1]-3)*(x[1]-3)+(x[0]-4)*(x[0]-4)+(x[2]-0.5)*(x[2]-0.5))'+\n",
" '<= Radius_Droplet ? Phi_tot_int : '+\n",
" 'sqrt((x[1]-4)*(x[1]-4)+(x[0]-5)*(x[0]-5)+(x[2]-0.5)*(x[2]-0.5))'+\n",
" '<= Radius_Droplet ? Phi_tot_int : '+\n",
" 'sqrt((x[1]-4)*(x[1]-4)+(x[0]-3)*(x[0]-3)+(x[2]-0.5)*(x[2]-0.5))'+\n",
" '<= Radius_Droplet ? Phi_tot_int : Phi_tot_ext',\n",
Lars Hubatsch
committed
" 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",
Lars Hubatsch
committed
"Phi_tot = Expression('sqrt((x[1]-4)*(x[1]-4)+(x[0]-4)*(x[0]-4)+(x[2]-0.5)*(x[2]-0.5))'+\n",
Lars Hubatsch
committed
" ' <= Radius_Droplet ? Phi_tot_int : '+\n",
Lars Hubatsch
committed
" 'sqrt((x[1]-5)*(x[1]-5)+(x[0]-4)*(x[0]-4)+(x[2]-0.5)*(x[2]-0.5))'+\n",
" '<= Radius_Droplet ? Phi_tot_int: '+\n",
" 'sqrt((x[1]-3)*(x[1]-3)+(x[0]-4)*(x[0]-4)+(x[2]-0.5)*(x[2]-0.5))'+\n",
" '<= Radius_Droplet ? Phi_tot_int: '+\n",
" 'sqrt((x[1]-4)*(x[1]-4)+(x[0]-5)*(x[0]-5)+(x[2]-0.5)*(x[2]-0.5))'+\n",
" '<= Radius_Droplet ? Phi_tot_int: '+\n",
" 'sqrt((x[1]-4)*(x[1]-4)+(x[0]-3)*(x[0]-3)+(x[2]-0.5)*(x[2]-0.5))'+\n",
" '<= Radius_Droplet ? Phi_tot_int: Phi_tot_ext', \n",
Lars Hubatsch
committed
" 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",
Lars Hubatsch
committed
"# 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",
Lars Hubatsch
committed
"\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",
Lars Hubatsch
committed
"cFile = XDMFFile('multi_drop_6_neighbors_med.xdmf')\n",
"cFile.write(Phi_u0, t)\n",
"\n",
Lars Hubatsch
committed
"# c_in_b = []\n",
Lars Hubatsch
committed
"# c_in_b1 = []\n",
"c_in_6n_med = []\n",
"ti = time.time()\n",
Lars Hubatsch
committed
"for i in range(100):\n",
Lars Hubatsch
committed
" c_in_6n_med.append([Phi_u([x, 4, 0.5]) for x in np.linspace(4, 4.49, 100)])\n",
Lars Hubatsch
committed
" 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",
Lars Hubatsch
committed
"cFile.close()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
Lars Hubatsch
committed
"source": [
Lars Hubatsch
committed
"# plt.plot(c_in_b)\n",
"# plt.plot(c_in_b1)\n",
"plt.plot(np.transpose(c_in_6n), c='r')\n",
"plt.plot(np.transpose(c_in_6n_med), c='k')\n",
"plt.plot(np.transpose(c_in_6n_far), c='g')\n",
"plt.plot(np.transpose(c_in_0n), c='b')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"plt.plot([np.mean(x[1:49]) for x in c_in_0n])\n",
"plt.plot([np.mean(x[1:49]) for x in c_in_6n])\n",
"plt.plot([np.mean(x[1:49]) for x in c_in_6n_far])\n",
"plt.plot([np.mean(x[1:49]) for x in c_in_6n_med])"
Lars Hubatsch
committed
]
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
},
{
"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
}