Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit edf05acf by Jigyasa Watwani
Options

added flags for morphogen-mediated active stress, turnover, average density

parent 164b33e9
Inline Side-by-side
Showing with 103 additions and 27 deletions
growing_domain/full_model/growing_domain.py
......@@ -2,6 +2,8 @@ import dolfin as df
import numpy as np
import os
import progressbar
import meshzoo
import matplotlib.pyplot as plt
df.set_log_level(df.LogLevel.ERROR)
df.parameters['form_compiler']['optimize'] = True
......@@ -19,14 +21,18 @@ class Growing_Domain(object):
assert self.bulk_elasticity == self.shear_elasticity
assert self.bulk_viscosity == self.shear_viscosity
elif self.dimension==2:
if hasattr(self, 'meshfile'):
print('Using %s' % str(self.meshfile))
self.mesh = df.Mesh(str(self.meshfile))
else:
self.mesh = df.RectangleMesh(df.Point(0,0), df.Point(L,L),
self.resolution, self.resolution)
# self.mesh = df.EllipseMesh(df.Point(0.0, 0.0), [3.0, 1.0], 0.2)
points, cells = meshzoo.disk(10,8)
self.mesh = df.Mesh()
e = df.MeshEditor()
e.open(self.mesh, type='triangle', tdim=2, gdim=2)
e.init_vertices(len(points))
e.init_cells(len(cells))
for i in range(len(points)):
e.add_vertex(i, points[i])
for i in range(len(cells)):
e.add_cell(i, cells[i])
e.close()
scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
vector_element = df.VectorElement('P', self.mesh.ufl_cell(), 1)
if self.morphogen:
......@@ -48,11 +54,14 @@ class Growing_Domain(object):
self.f0 = df.Function(self.function_space) # f at t =0
def advection(self, conc, vel, tconc):
return (df.inner(df.div(vel * conc), tconc))
return df.inner(df.div(vel * conc), tconc)
def active_stress(self, rho, c):
return self.lamda * rho/(rho + self.saturation_rho) * c * df.Identity(self.dimension)
if self.morphogen_mediated_active_stress:
return self.lamda * rho/(rho + self.saturation_rho) * (c-self.average_c) * df.Identity(self.dimension)
else:
return self.lamda * rho/(rho + self.saturation_rho)* df.Identity(self.dimension)
def epsilon(self, v):
return df.sym(df.nabla_grad(v))
......@@ -73,7 +82,12 @@ class Growing_Domain(object):
return (self.passive_stress(u, v) + self.active_stress(rho, c))
def reaction_rho(self, rho, trho):
return self.turnover_rho * df.inner(rho - self.average_rho, trho)
if self.morphogen_mediated_turnover_rho:
return self.average_c * df.inner(rho - self.average_rho, trho)
elif self.morphogen_mediated_average_rho:
return self.turnover_rho * df.inner(rho - self.average_c, trho)
else:
return self.turnover_rho * df.inner(rho - self.average_rho, trho)
def reaction_c(self, c, tc):
return self.turnover_c * df.inner(c - self.average_c, tc)
......@@ -99,18 +113,34 @@ class Growing_Domain(object):
# ic for rho
r2 = "+".join(['x['+str(i)+']*x['+str(i)+']' for i in range(self.dimension)])
# rho0 = df.interpolate(df.Constant(0.0), self.function_space.sub(2).collapse())
rho0 = 'a*x[0]*x[0]*(1-tanh((sqrt(%s) - r0)/w))' %r2
rho0 = df.interpolate(df.Expression(rho0, a=self.average_rho, r0=0.5, w=0.2, degree=1), self.function_space.sub(2).collapse())
if self.dimension==1:
rho0='1-cos(2*pi*x[0]/L)'
else:
# circular domain --> circular domain
# rho0 = '1 + cos(pi*%s/R)' %r2
rho0 = 'a*(1-tanh((sqrt(%s) - r0)/w))' %r2
# circular domain --> elliptical domain (w=0.3, lambda = -30, k=1)
# rho0 = 'a*x[0]*x[0]*(1-tanh((sqrt(%s) - r0)/w))' %r2
# circular domain --> dumbbell shaped domain (w = 0.8, k = 0.01, lambda = -80)
# rho0 = '0.25 + a*x[0]*x[1]*(1-tanh((sqrt(%s) - r0)/w))' %r2
# circular --> kidney (w=0.8, r0=0.5, k=0.1, lambda=-30)
# rho0 = '0.6 + x[0]*(1-tanh((sqrt(%s) - r0)/w))' %r2
rho0 = df.interpolate(df.Expression(rho0, pi=np.pi, R=self.system_size, a=self.average_rho, r0=0.5, w=0.3, degree=1), self.function_space.sub(2).collapse())
# add noise
noise_rho = self.noise_level * (2*np.random.random(rho0.vector().size())-1)
noise_rho = self.noise_level * (2 * np.random.random(rho0.vector().size()) - 1)
rho0.vector()[:] = rho0.vector()[:] + noise_rho
# ic for c
if self.morphogen:
c0 = 'a*(1 - tanh((sqrt(%s) - r0)/w))/2.0' % r2
# c0 = df.interpolate(df.Expression(c0, L=self.system_size, pi=np.pi, degree=1), self.function_space.sub(3).collapse())
c0 = df.interpolate(df.Expression(c0, a=self.average_rho, r0=0.5, w=0.1, degree=1), self.function_space.sub(3).collapse())
c0 = '1 + 0.1*a*cos(2.0*pi*%s/R)' %r2
# c0 = df.interpolate(df.Constant(0.0), self.function_space.sub(3).collapse())
c0 = df.interpolate(df.Expression(c0,a = self.average_c, R=self.system_size, pi=np.pi, degree=1), self.function_space.sub(3).collapse())
# add noise
noise_c = self.noise_level * (2*np.random.random(c0.vector().size())-1)
c0.vector()[:] = c0.vector()[:] + noise_c
......@@ -188,6 +218,10 @@ class Growing_Domain(object):
self.uFile.write_checkpoint(u, 'displacement', self.time, append=True)
self.vFile.write_checkpoint(v, 'velocity', self.time, append=True)
self.rhoFile.write_checkpoint(rho, 'density', self.time, append=True)
radius_array = np.zeros(maxsteps+1)
radius_array[0] = np.max(np.sqrt(self.mesh.coordinates()[:][0]**2 + self.mesh.coordinates()[:][1]**2))
for steps in progressbar.progressbar(range(1, maxsteps + 1)):
# solve
df.solve(self.form == 0, self.f, self.bc)
......@@ -209,15 +243,16 @@ class Growing_Domain(object):
# move mesh
dr = df.project(u, self.function_space.sub(0).collapse())
df.ALE.move(self.mesh, dr)
radius_array[steps] = np.max(np.sqrt(self.mesh.coordinates()[:][0]**2 + self.mesh.coordinates()[:][1]**2))
# update time
self.time += self.timestep
self.uFile.close()
self.vFile.close()
self.rhoFile.close() # def save_data_uvrho(self):
if self.morphogen:
self.cFile.close()
return (radius_array)
if __name__ == '__main__':
import json, datetime
......@@ -235,7 +270,18 @@ if __name__ == '__main__':
params['timestamp'] = timestamp
gd = Growing_Domain(params)
gd.solve()
radius_array = gd.solve()
times = np.linspace(0, params['maxtime'], int(params['maxtime']/params['timestep']) + 1 )
fig,ax = plt.subplots(1,1)
ax.set_yscale("log")
ax.set_ylabel("$log R(t)$")
ax.set_xlabel("$t$")
ax.scatter(times, radius_array, facecolors='none', edgecolors='b')
plt.show()
with open(params['timestamp'] + '_parameters.json', "w") as fp:
json.dump(params, fp, indent=4)
\ No newline at end of file
json.dump(params, fp, indent=4)
from viz_growing_domain import visualize
visualize(params, DIR="", offscreen=False)
\ No newline at end of file
growing_domain/full_model/viz_growing_domain.py
......@@ -131,7 +131,8 @@ def visualize(params, DIR='', offscreen=False):
geometry = np.dstack((geometry, np.zeros(geometry.shape[0:2])))
rhomin, rhomax = np.min(rho), np.max(rho)
if params['morphogen']:
cmin, cmax = np.min(c), np.max(c)
plotter = vd.plotter.Plotter(axes=0)
poly = vd.utils.buildPolyData(geometry[0], topology)
......@@ -139,7 +140,7 @@ def visualize(params, DIR='', offscreen=False):
#scalar_actor.computeNormals(points=True, cells=True)
scalar_actor.pointdata['density'] = rho[0]
scalar_actor.cmap(scalar_cmap, rho[0], vmin=rhomin, vmax=rhomax, n=n_cmap_vals)
scalar_actor.add_scalarbar(title = r'$\rho$',
scalar_actor.add_scalarbar(title = r'$\rho/\rho_0$',
pos=(0.8, 0.04), nlabels=2,
# titleYOffset=15, titleFontSize=28, size=(100, 600)
)
......@@ -160,6 +161,35 @@ def visualize(params, DIR='', offscreen=False):
vd.show(interactive=(not offscreen), zoom=0.8)
if params['morphogen']:
plotter1 = vd.plotter.Plotter(axes=0)
poly1 = vd.utils.buildPolyData(geometry[0], topology)
scalar_actor1 = vd.mesh.Mesh(poly1)
#scalar_actor1.computeNormals(points=True, cells=True)
scalar_actor1.pointdata['concentration'] = c[0]
scalar_actor1.cmap(scalar_cmap, c[0], vmin=cmin, vmax=cmax, n=n_cmap_vals)
scalar_actor1.add_scalarbar(title = r'$c/c_0$',
pos=(0.8, 0.04), nlabels=2,
# titleYOffset=15, titleFontSize=28, size=(100, 600)
)
plotter1 += scalar_actor1
def update1(idx):
scalar_actor1.points(pts=geometry[idx], transformed=False)
scalar_actor1.pointdata['concentration'] = c[idx]
def slider_update1(widget, event):
value = widget.GetRepresentation().GetValue()
idx = (abs(times-value)).argmin()
update1(idx)
slider1 = plotter1.add_slider(slider_update1, pos=[(0.1,0.94), (0.5,0.94)],
xmin=times[0], xmax=times.max(),
value=times[0], title=r"$t/\tau$")
vd.show(interactive=(not offscreen), zoom=0.8)
if offscreen:
FPS = 10
movFile = '%s_numerical.mov' % params['timestamp']
......@@ -199,7 +229,7 @@ def visualize(params, DIR='', offscreen=False):
fc='#999999')
slider1.drawon = False
slider1.on_changed(update)
plt.show()
# plt.show()
if __name__ == '__main__':
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or sign in to comment