Skip to content
Snippets Groups Projects
Commit 02a3d2fb authored by Lars Hubatsch's avatar Lars Hubatsch
Browse files

Adding Gamma_0 / mobility to FP equation.

parent 88afc875
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Half_Frap_Phi_u.ipynb
+ 37
− 9
View file @ 02a3d2fb
%% Cell type:markdown id: tags:
### FRAP geometries
%% Cell type:code id: tags:
``` python
from fem_sol import frap_solver
import matplotlib.pyplot as plt
import numpy as np
```
%% Cell type:code id: tags:
``` python
font = {'family' : 'normal',
'weight' : 'normal',
'size' : 12}
import matplotlib
matplotlib.rc('font', **font)
def make_graph_pretty(xlab, ylab, loc=0, markerfirst=True, handlelength=None):
color1 = 'black'#'darkorange'
color2 = (0, 0, 0)
color2 = (1, 1, 1)
ax = plt.gca()
ax.tick_params(axis='x', colors=color1, which='both')
ax.tick_params(axis='y', colors=color1, which='both')
plt.xlabel(xlab, color=color1)
plt.ylabel(ylab, color=color1)
ax.spines['left'].set_color(color1)
ax.spines['bottom'].set_color(color1)
# Hide the right and top spines
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.set_facecolor(color2)
plt.gcf().subplots_adjust(bottom=0.17, left=0.15)
plt.legend(frameon=False, loc=loc, labelspacing=0.1,
markerfirst=markerfirst, handlelength=handlelength)
```
%% Cell type:code id: tags:
``` python
m = ['Meshes/multi_drop_gauss.xml', 'Meshes/multi_drop_gauss_med.xml',
'Meshes/multi_drop_gauss_far.xml']
point_lists = [[[4, 4.5, 0.5], [4, 3.5, 0.5], [3.5, 4, 0.5], [4.5, 4, 0.5]],
[[4, 5, 0.5], [4, 3, 0.5], [3, 4, 0.5], [5, 4, 0.5]],
[[4, 5.5, 0.5], [4, 2.5, 0.5], [2.5, 4, 0.5], [5.5, 4, 0.5]]]
f_i = []
for j in range(len(m)):
f = frap_solver([4, 4, 0.5], m[j], name=str(j), point_list=point_lists[j],
T=3)
f.solve_frap()
f_i.append(f)
```
%% Cell type:code id: tags:
``` python
for i in range(len(cs[0])):
plt.plot([np.mean(x[1:49])/0.8 for x in cs[0][0:i]], label='dist=0.5')
plt.plot([np.mean(x[1:49])/0.8 for x in cs[1][0:i]], label='dist=1')
plt.plot([np.mean(x[1:49])/0.8 for x in cs[2][0:i]], label='dist=1.5')
plt.plot(range(0, 100), np.ones(100), linestyle='--', color='k')
plt.ylim((0, 1.1))
plt.xlim((0, 100))
make_graph_pretty('t (a.u.)', 'intensity (a.u)', loc=4)
# plt.savefig('img_'+str(i)+'.png', dpi=300)
plt.show()
```
%% Cell type:code id: tags:
``` python
f_coverslip = frap_solver([4, 4, 0.5], 'Meshes/single_drop_coverslip.xml',
name='FRAP_coverslip', T=300)
f_coverslip.solve_frap()
f_symmetric = frap_solver([4, 4, 4], 'Meshes/single_drop_symmetric.xml',
name='FRAP_symmetric', T=300)
f_symmetric.solve_frap()
f_0 = frap_solver([4, 4, 0.5], 'Meshes/single_drop_coverslip.xml',
name='FRAP_coverslip', T=300, phi_tot_ext=0.001)
f_0.solve_frap()
f_1 = frap_solver([4, 4, 1.5], 'Meshes/single_drop_1_5.xml',
name='FRAP_symmetric', T=300, phi_tot_ext=0.001)
f_1.solve_frap()
f_2 = frap_solver([4, 4, 4], 'Meshes/single_drop_symmetric.xml',
name='FRAP_symmetric', T=300, phi_tot_ext=0.001)
f_2.solve_frap()
```
%% Cell type:code id: tags:
``` python
f_3 = frap_solver([4, 4, 0.5], 'Meshes/single_drop_coverslip.xml',
name='FRAP_symmetric', T=100, phi_tot_ext=0.01)
f_3.solve_frap()
f_4 = frap_solver([4, 4, 1.5], 'Meshes/single_drop_1_5.xml',
name='FRAP_symmetric', T=100, phi_tot_ext=0.01)
f_4.solve_frap()
f_5 = frap_solver([4, 4, 4], 'Meshes/single_drop_symmetric.xml',
name='FRAP_symmetric', T=100, phi_tot_ext=0.01)
f_5.solve_frap()
```
%% Cell type:code id: tags:
``` python
plt.plot([np.mean(x[1:49])/0.8 for x in f_coverslip.cs], label='coverslip')
plt.plot([np.mean(x[1:49])/0.8 for x in f_symmetric.cs], label='symmetric')
plt.plot([np.mean(x[1:49])/0.99 for x in f_0.cs], label='dist=0.5')
plt.plot([np.mean(x[1:49])/0.99 for x in f_1.cs], label='dist=1.5')
plt.plot([np.mean(x[1:49])/0.99 for x in f_2.cs], label='dist=4')
plt.plot([np.mean(x[1:49])/0.99 for x in f_3.cs], label='dist=0.5, G_out=10')
plt.plot([np.mean(x[1:49])/0.99 for x in f_4.cs], label='dist=1.5, G_out=10')
plt.plot([np.mean(x[1:49])/0.99 for x in f_5.cs], label='dist=4, G_out=10')
plt.plot(range(0, 300), np.ones(300), linestyle='--', color='k')
plt.ylim((0, 1.1))
plt.xlim((0, 300))
make_graph_pretty('t (a.u.)', 'intensity (a.u)', loc=4)
```
%% Cell type:code id: tags:
``` python
f_coverslip.cs
plt.plot(np.transpose(f_3.cs))
```
%% Cell type:code id: tags:
``` python
```
......
Meshes/single_drop_1_5.geo 0 → 100644
+ 14
− 0
View file @ 02a3d2fb
SetFactory("OpenCASCADE");
Box(1) = {0, 0, 0, 8, 8, 8};
Physical Surface("surfacedomain") = {1};
Physical Volume("volumedomain") = {1};
Field[3] = MathEval;
Field[3].F = ".5 - (.5-0.01)*exp(-(1./0.025) *((sqrt((x-4)*(x-4)+(y-4)*(y-4)+(z-1.5)*(z-1.5))-0.25))*((sqrt((x-4)*(x-4)+(y-4)*(y-4)+(z-1.5)*(z-1.5))-0.25)))";
Field[8] = Min;
Field[8].FieldsList = {3};
Mesh.Algorithm = 1;
Background Field = 8;
\ No newline at end of file
Meshes/single_drop_coverslip.geo
+ 1
− 1
View file @ 02a3d2fb
......@@ -10,5 +10,5 @@ Field[3].F = ".5 - (.5-0.01)*exp(-(1./0.025) *((sqrt((x-4)*(x-4)+(y-4)*(y-4)+(z-
Field[8] = Min;
Field[8].FieldsList = {3};
Mesh.Algorithm = 5;
Mesh.Algorithm = 1;
Background Field = 8;
\ No newline at end of file
fem_sol.py
+ 16
− 7
View file @ 02a3d2fb
......@@ -6,7 +6,7 @@ import time
set_log_level(40)
def point_expr(p_c, p_list, phi_tot_int, phi_tot_ext, Phi_tot_init,
def point_expr(p_c, p_list, phi_tot_int, phi_tot_ext, phi_init,
rad_drop, tol):
'''
Create dolfin expression for initial condition of Phi or
......@@ -16,7 +16,7 @@ def point_expr(p_c, p_list, phi_tot_int, phi_tot_ext, Phi_tot_init,
'''
s = ('sqrt(pow((x[0]-'+str(p_c[0])+'),2)+pow((x[1]-'+str(p_c[1])+'),2)+'+
'pow((x[2]-'+str(p_c[2])+'),2))')
s = s + '<= '+str(rad_drop)+' ? '+str(Phi_tot_init)+' : '
s = s + '<= '+str(rad_drop)+' ? '+str(phi_init)+' : '
for p in p_list:
y = ('sqrt(pow((x[0]-'+str(p[0])+'),2)+pow((x[1]-'+str(p[1])+'),2)+'+
'pow((x[2]-'+str(p[2])+'),2))')
......@@ -27,9 +27,9 @@ def point_expr(p_c, p_list, phi_tot_int, phi_tot_ext, Phi_tot_init,
class frap_solver:
def __init__(self, cent_poin, mesh, tol=1E-14, dt=0.1, name='FRAP',
phi_tot_int=0.99, phi_tot_ext=0.001, phi_tot_initial=0,
point_list=[], rad_drop=0.25, T=100):
def __init__(self, cent_poin, mesh, tol=1E-14, dt=0.1, G_in=1,
G_out=1, name='FRAP', phi_tot_int=0.99, phi_tot_ext=0.001,
phi_tot_initial=0, point_list=[], rad_drop=0.25, T=100):
'''
Initialise solver
'''
......@@ -38,6 +38,8 @@ class frap_solver:
self.name = name # name of output file
self.tol = tol # Tolerance for expressions
self.dt = dt # Time step
self.G_in = G_in # Gamma inside, mobility
self.G_out = G_out # Gamma outside, mobility
self.phi_tot_int = phi_tot_int # Value inside droplet
self.phi_tot_ext = phi_tot_ext # Value outside droplet
self.phi_tot_initial = phi_tot_initial # Initial value inside droplet
......@@ -67,10 +69,17 @@ class frap_solver:
self.phi_tot_ext, self.phi_tot_int, self.rad_drop,
self.tol)
Phi_tot = interpolate(Phi_tot, V)
Gamma_0 = interpolate(Expression('(x[0]-xc)*(x[0]-xc)+'
'(x[1]-yc)*(x[1]-yc)+'+
'(x[2]-zc)*(x[2]-zc) <= rad_drop*rad_drop ? '+
'G_in : G_out' ,
degree=2, rad_drop=self.rad_drop, G_in=self.G_in,
G_out=self.G_out, xc=self.cent_poin[0],
yc=self.cent_poin[1], zc=self.cent_poin[2]), V)
Form = (inner((Phi_u-Phi_u0)/self.dt, v) -
inner(((1-Phi_tot)/Phi_tot)*Phi_u*grad(Phi_tot)-
(1-Phi_tot)*grad(Phi_u), grad(v)))*dx
inner(Gamma_0*((1-Phi_tot)/Phi_tot)*Phi_u*grad(Phi_tot)-
Gamma_0*(1-Phi_tot)*grad(Phi_u), grad(v)))*dx
t = 0
cFile = XDMFFile(self.name+'.xdmf')
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment