Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit e38db8a2 by Jigyasa Watwani
Options

fully working 2D problem with remeshing and extend functionality

parent 82d74da7
Inline Side-by-side
Showing with 313 additions and 510 deletions
growing_domain/full_model/growing_domain.py
import dolfin as df import dolfin as df
import numpy as np import numpy as np
import os
import progressbar import progressbar
import meshzoo import meshzoo
import matplotlib.pyplot as plt
df.set_log_level(df.LogLevel.ERROR) df.set_log_level(df.LogLevel.ERROR)
df.parameters['form_compiler']['optimize'] = True df.parameters['form_compiler']['optimize'] = True
class Growing_Domain(object): class Growing_Domain(object):
def __init__(self, parameters): def __init__(self, parameters, mesh = None):
# read in parameters # read in parameters
for key in parameters: for key in parameters:
setattr(self, key, parameters[key]) setattr(self, key, parameters[key])
# create mesh, mixed element and function space
L = self.system_size # if no mesh is given as initial condition, create one
if self.dimension==1: if mesh is None:
self.mesh = df.IntervalMesh(self.resolution, 0.0, L)
assert self.bulk_elasticity == self.shear_elasticity
assert self.bulk_viscosity == self.shear_viscosity
elif self.dimension==2:
points, cells = meshzoo.disk(10,8) points, cells = meshzoo.disk(10,8)
self.mesh = df.Mesh() self.mesh = df.Mesh()
e = df.MeshEditor() e = df.MeshEditor()
...@@ -33,13 +27,13 @@ class Growing_Domain(object): ...@@ -33,13 +27,13 @@ class Growing_Domain(object):
e.add_cell(i, cells[i]) e.add_cell(i, cells[i])
e.close() e.close()
else:
self.mesh = mesh
def setup(self):
scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1) scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
vector_element = df.VectorElement('P', self.mesh.ufl_cell(), 1) vector_element = df.VectorElement('P', self.mesh.ufl_cell(), 1)
if self.morphogen:
# u, v, rho, c
mixed_element = df.MixedElement([vector_element, vector_element,
scalar_element, scalar_element])
else:
# u, v, rho # u, v, rho
mixed_element = df.MixedElement([vector_element, vector_element, scalar_element]) mixed_element = df.MixedElement([vector_element, vector_element, scalar_element])
...@@ -47,12 +41,7 @@ class Growing_Domain(object): ...@@ -47,12 +41,7 @@ class Growing_Domain(object):
self.function_space = df.FunctionSpace(self.mesh, mixed_element) self.function_space = df.FunctionSpace(self.mesh, mixed_element)
# dirichlet boundaries for rho # dirichlet boundaries for rho
if self.dirichlet_boundary_rho:
if self.zero_rho_boundary:
self.bc = df.DirichletBC(self.function_space.sub(2), df.Constant(0.0), 'on_boundary') self.bc = df.DirichletBC(self.function_space.sub(2), df.Constant(0.0), 'on_boundary')
elif self.rho0_boundary:
self.bc = df.DirichletBC(self.function_space.sub(2), df.Constant(self.average_rho), 'on_boundary')
# define the reqd functions on this function space # define the reqd functions on this function space
self.f = df.Function(self.function_space) # f at current time self.f = df.Function(self.function_space) # f at current time
...@@ -61,21 +50,19 @@ class Growing_Domain(object): ...@@ -61,21 +50,19 @@ class Growing_Domain(object):
def advection(self, conc, vel, tconc): 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, v): def qtensor(self):
if self.saturating_density_active_stress: q = df.Expression((('q_xx', 'q_xy'),
return self.lamda * rho/(rho + self.saturation_rho)* df.Identity(self.dimension) ('q_yx', 'q_yy')),
q_xx = self.q_xx, q_xy = self.q_xy,
elif self.difference_in_morphogen_active_stress: q_yx = self.q_xy, q_yy = -self.q_xx, degree=1)
return self.lamda * (c-self.average_c) * df.Identity(self.dimension) return q
elif self.difference_in_density_active_stress:
return self.lamda * (c-self.average_rho) * df.Identity(self.dimension)
elif self.saturating_density_and_difference_in_density_active_stress:
return self.lamda * (rho*(rho-self.average_rho))/(rho + self.saturation_rho) * df.Identity(self.dimension)
elif self.saturating_density_and_difference_in_morphogen_active_stress: def active_stress(self, rho, u, v):
return self.lamda * (rho*(c-self.average_c))/(rho + self.saturation_rho) * df.Identity(self.dimension) return -(self.lamda * rho * (rho - self.stress_setpoint_rho)/(rho**2 + self.saturation_rho**2)
# self.qtensor()
* df.Identity(self.dimension)
)
def epsilon(self, v): def epsilon(self, v):
return df.sym(df.nabla_grad(v)) return df.sym(df.nabla_grad(v))
...@@ -93,42 +80,28 @@ class Growing_Domain(object): ...@@ -93,42 +80,28 @@ class Growing_Domain(object):
return (elastic_stress + viscous_stress) return (elastic_stress + viscous_stress)
def stress(self, u, v, rho, c): def stress(self, u, v, rho):
return (self.passive_stress(u, v) + self.active_stress(rho, c, v)) return (self.passive_stress(u, v) + self.active_stress(rho, u, v))
def reaction_rho(self, rho, trho, c):
if self.morphogen_mediated_turnover_rho: def reaction_rho(self, rho, trho):
return self.turnover_rho * df.inner((c - self.average_c/2)*(rho - self.average_rho), trho) return -self.turnover_rho * df.inner(rho * (1 - rho/self.average_rho), trho)
elif self.morphogen_mediated_average_rho: def diffusion_reaction_rho(self, rho, trho):
return self.turnover_rho * df.inner(rho - self.average_c, trho)
elif self.density_dependent_turnover_rho:
# to model that cells are born or die only where the density is non-zero
# tanh approximation: Theta = 1+ tanh = 1 + rho - rho^3/3! + . . .
return self.turnover_rho * df.inner((1 + rho) * (rho - self.average_rho), 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)
def diffusion_reaction_rho(self, rho, trho, c):
return (self.diffusion_rho * df.inner(df.nabla_grad(rho), df.nabla_grad(trho)) return (self.diffusion_rho * df.inner(df.nabla_grad(rho), df.nabla_grad(trho))
+ self.reaction_rho(rho, trho, c)) + self.reaction_rho(rho, trho))
def refine_mesh(self, mesh):
cell_markers = df.MeshFunction("bool", mesh, mesh.topology().dim())
cell_markers.set_all(False)
for cell in df.cells(mesh):
if cell.volume() > 2 * 0.005:
cell_markers[cell] = True
mesh = df.refine(mesh, cell_markers)
def diffusion_reaction_c(self, c, tc): return mesh
return (self.diffusion_c * df.inner(df.nabla_grad(c), df.nabla_grad(tc))
+ self.reaction_c(c, tc))
def setup_initial_conditions(self): def setup_initial_conditions(self):
# ic for u # ic for u
if self.dimension==1:
self.zero_vector = df.Constant((0.0,))
elif self.dimension==2:
self.zero_vector = df.Constant((0.0, 0.0)) self.zero_vector = df.Constant((0.0, 0.0))
u0 = df.interpolate(self.zero_vector, self.function_space.sub(0).collapse()) u0 = df.interpolate(self.zero_vector, self.function_space.sub(0).collapse())
# add noise # add noise
...@@ -136,181 +109,124 @@ class Growing_Domain(object): ...@@ -136,181 +109,124 @@ class Growing_Domain(object):
u0.vector()[:] = u0.vector()[:] + noise_u u0.vector()[:] = u0.vector()[:] + noise_u
# ic for rho # ic for rho
rho0 = '(a/2)*(1-tanh((sqrt(x[0]*x[0]+x[1]*x[1]) - r0)/w))'
r2 = "+".join(['x['+str(i)+']*x['+str(i)+']' for i in range(self.dimension)]) rho0 = df.interpolate(df.Expression(rho0, a=self.average_rho, r0=0.5, w=0.1, degree=1),
self.function_space.sub(2).collapse())
if self.dimension==1:
rho0='1-cos(2*pi*x[0]/R)'
else:
# for rho=rho_0 boundary condition
# rho0 = "%s" % str(self.average_rho)
# rho0 = '(a/2)*(1 + tanh((sqrt(%s) - r0)/w))' %r2
# grad rho=0 boundary condition
rho0 = '1 + 0.1 * a * cos(2.0*pi*%s/R)' %r2
# rho = 0 boundary condition
# circular domain --> circular domain
# rho0 = '(a/2)*(1 - tanh((sqrt(%s) - r0)/w))' %r2
# rho0 = 'a*(1-tanh((sqrt(%s) - r0)/w))' %r2 # grows to a different size
# circular domain --> elliptical domain (w=0.19, lambda = -10, 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 = '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.19, degree=1), self.function_space.sub(2).collapse())
# add noise # 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 rho0.vector()[:] = rho0.vector()[:] + noise_rho
# ic for c
if self.morphogen:
# c0 = "%s" % str(5.0)
# c0 = '1 + 0.1*a*cos(2.0*pi*sqrt(%s)/R)' %r2 # shrinks and stops
c0 = '1 + 0.1*a*cos(2.0*pi*%s/R)' %r2 # grows and stops
# c0 = '1 + 0.1*a*cos(2.0*pi*x[0]/R)' # grows very slightly, shrinks even more slightly, changes shape
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
else:
c0 = df.Constant(1.0)
VFS = self.function_space.sub(1).collapse() VFS = self.function_space.sub(1).collapse()
v0 = df.Function(VFS) v0 = df.Function(VFS)
tv = df.TestFunction(VFS) tv = df.TestFunction(VFS)
vform = (self.friction * df.inner(v0, tv) vform = (self.friction * df.inner(v0, tv)
+ df.inner(self.stress(u0, v0, rho0, c0), self.epsilon(tv)) + df.inner(self.stress(u0, v0, rho0), self.epsilon(tv))
) * df.dx ) * df.dx
df.solve(vform == 0, v0) df.solve(vform == 0, v0)
if self.morphogen:
df.assign(self.f0, [u0, v0, rho0, c0])
else:
df.assign(self.f0, [u0, v0, rho0]) df.assign(self.f0, [u0, v0, rho0])
def setup_weak_forms(self): def setup_weak_forms(self):
if self.morphogen:
u0, v0, rho0, c0 = df.split(self.f0)
u, v, rho, c = df.split(self.f)
tu, tv, trho, tc = df.TestFunctions(self.function_space)
else:
u0, v0, rho0 = df.split(self.f0) u0, v0, rho0 = df.split(self.f0)
u, v, rho = df.split(self.f) u, v, rho = df.split(self.f)
tu, tv, trho = df.TestFunctions(self.function_space) tu, tv, trho = df.TestFunctions(self.function_space)
c = df.Constant(1.0)
uform = (df.inner((u - u0)/self.timestep, tu) uform = (df.inner((u - u0)/self.timestep, tu)
- df.inner(v0, tu) - df.inner(v0, tu)
) )
vform = (self.friction * df.inner(v, tv) vform = (self.friction * df.inner(v, tv)
+ df.inner(self.stress(u, v, rho, c), self.epsilon(tv)) + df.inner(self.stress(u, v, rho), self.epsilon(tv))
) )
rhoform = (df.inner((rho - rho0)/self.timestep, trho) rhoform = (df.inner((rho - rho0)/self.timestep, trho)
+ self.advection(rho0, v0, trho) + self.advection(rho0, v0, trho)
+ self.diffusion_reaction_rho(rho, trho, c) + self.diffusion_reaction_rho(rho, trho)
) )
if self.morphogen:
cform = (df.inner((c - c0)/self.timestep, tc)
+ self.advection(c0, v0, tc)
+ self.diffusion_reaction_c(c, tc)
)
self.form = (uform + vform + rhoform + cform) * df.dx
else:
self.form = (uform + vform + rhoform) * df.dx self.form = (uform + vform + rhoform) * df.dx
def solve(self,DIR=''): def solve(self, initu=None, initv=None, initrho=None, initTime = 0, extend=False):
# for saving data
self.uFile = df.XDMFFile(os.path.join(DIR, '%s_displacement.xdmf' % self.timestamp))
self.vFile = df.XDMFFile(os.path.join(DIR, '%s_velocity.xdmf' % self.timestamp))
self.rhoFile = df.XDMFFile(os.path.join(DIR, '%s_density.xdmf' % self.timestamp))
if self.morphogen:
self.cFile = df.XDMFFile(os.path.join(DIR, '%s_concentration.xdmf' % self.timestamp))
self.setup_initial_conditions() # create function space, functions, bc
self.setup_weak_forms() self.setup()
# create files if they don't exist, open them if they exist
self.uFile = df.XDMFFile('%s_displacement.xdmf' % self.timestamp)
self.vFile = df.XDMFFile('%s_velocity.xdmf' % self.timestamp)
self.rhoFile = df.XDMFFile('%s_density.xdmf' % self.timestamp)
# time-variables # time-variables
self.time = 0.0 self.time = initTime
savesteps = int(self.savetime/self.timestep) savesteps = round(self.savetime/self.timestep)
maxsteps = int(self.maxtime/self.timestep) maxsteps = round(self.maxtime/self.timestep)
if self.morphogen:
u, v, rho, c = self.f0.split(deepcopy=True) if not extend:
self.uFile.write_checkpoint(u, 'displacement', self.time, append=True) # find the initial velocity
self.vFile.write_checkpoint(v, 'velocity', self.time, append=True) self.setup_initial_conditions()
self.rhoFile.write_checkpoint(rho, 'density', self.time, append=True)
self.cFile.write_checkpoint(c, 'concentration', self.time, append=True) # save the ic
u0, v0, rho0 = self.f0.split(deepcopy=True)
self.uFile.write_checkpoint(u0, 'displacement', self.time)
self.vFile.write_checkpoint(v0, 'velocity', self.time)
self.rhoFile.write_checkpoint(rho0, 'density', self.time)
# move the mesh with this initial velocity
dr0 = df.project(v0 * self.timestep, self.function_space.sub(0).collapse())
df.ALE.move(self.mesh, dr0)
else: else:
u, v, rho = self.f0.split(deepcopy=True) # assign fields to f0
self.uFile.write_checkpoint(u, 'displacement', self.time, append=True) u0 = df.project(initu, self.function_space.sub(0).collapse())
self.vFile.write_checkpoint(v, 'velocity', self.time, append=True) v0 = df.project(initv, self.function_space.sub(1).collapse())
self.rhoFile.write_checkpoint(rho, 'density', self.time, append=True) rho0 = df.project(initrho, self.function_space.sub(2).collapse())
df.assign(self.f0, [u0, v0, rho0])
# move the mesh with this initial velocity
dr0 = df.project(v0 * self.timestep, self.function_space.sub(0).collapse())
df.ALE.move(self.mesh, dr0)
for steps in progressbar.progressbar(range(1, maxsteps + 1)): self.setup_weak_forms()
# solve
if self.dirichlet_boundary_rho: bar = progressbar.ProgressBar(maxval=maxsteps)
for steps in bar(range(1, maxsteps + 1)):
self.time += self.timestep
# solve to get fields at next timestep
df.solve(self.form == 0, self.f, self.bc) df.solve(self.form == 0, self.f, self.bc)
elif self.neumann_boundary_rho:
if self.zero_grad_rho_boundary: # save the fields
df.solve(self.form==0, self.f) u, v, rho = self.f.split(deepcopy=True)
# update
self.f0.assign(self.f)
if self.morphogen:
u, v, rho, c = self.f0.split(deepcopy=True)
if steps % savesteps == 0:
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)
self.cFile.write_checkpoint(c, 'concentration', self.time, append=True)
else:
u, v, rho = self.f0.split(deepcopy=True)
if steps % savesteps == 0: if steps % savesteps == 0:
self.uFile.write_checkpoint(u, 'displacement', self.time, append=True) self.uFile.write_checkpoint(u, 'displacement', self.time, append=True)
self.vFile.write_checkpoint(v, 'velocity', self.time, append=True) self.vFile.write_checkpoint(v, 'velocity', self.time, append=True)
self.rhoFile.write_checkpoint(rho, 'density', self.time, append=True) self.rhoFile.write_checkpoint(rho, 'density', self.time, append=True)
# move mesh
dr = df.project(u, self.function_space.sub(0).collapse())
df.ALE.move(self.mesh, dr)
# update time
self.time += self.timestep
self.uFile.close() # assign the calculated func to a newly-defined function to be used later
self.vFile.close() old_function = df.Function(self.function_space)
self.rhoFile.close() # def save_data_uvrho(self): old_function.assign(self.f)
if self.morphogen: # move mesh
self.cFile.close() dr = df.project(v * self.timestep, self.function_space.sub(0).collapse())
df.ALE.move(self.mesh, dr)
if __name__ == '__main__': # refine the mesh
import json, datetime self.mesh = self.refine_mesh(self.mesh)
assert os.path.isfile('parameters.json'), 'parameters.json file not found' # assert that parameters.json is a valid file, otherwise
# give an error message parameters.json file not found
# load the parameters # new function space, bc, functions
with open('parameters.json') as jsonFile: self.setup()
params = json.load(jsonFile)
# parse parameters # new function0 should be the old function in the old function space projected on the new function space
assert params['dimension'] in (1,2) old_function.set_allow_extrapolation(True)
timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S") self.f0 = df.project(old_function, self.function_space)
params['timestamp'] = timestamp
gd = Growing_Domain(params) # new weak form
gd.solve() self.setup_weak_forms()
with open(params['timestamp'] + '_parameters.json', "w") as fp:
json.dump(params, fp, indent=4)
from viz_growing_domain import visualize self.uFile.close()
visualize(params, DIR="", offscreen=False) self.vFile.close()
\ No newline at end of file self.rhoFile.close()
growing_domain/full_model/parameters.json
{ {
"morphogen" : false, "morphogen" : false,
"dimension" : 2, "dimension" : 2,
"resolution" : 100,
"system_size" : 1,
"meshfile" : "disk.xml.gz",
"dirichlet_boundary_rho": false,
"zero_rho_boundary":false,
"rho0_boundary": false,
"neumann_boundary_rho": true,
"zero_grad_rho_boundary": true,
"difference_in_density_active_stress": false,
"difference_in_morphogen_active_stress": false,
"saturating_density_active_stress": true,
"saturating_density_and_difference_in_morphogen_active_stress": false,
"saturating_density_and_difference_in_density_active_stress":false,
"morphogen_mediated_turnover_rho": false, "resolution" : 128,
"density_dependent_turnover_rho": true, "system_size" : 1,
"morphogen_mediated_average_rho": false, "bulk_elasticity" : 1,
"shear_elasticity": 1,
"bulk_elasticity" : 1.0,
"shear_elasticity": 1.0,
"shear_viscosity" : 1.0, "shear_viscosity" : 1.0,
"bulk_viscosity" : 1.0, "bulk_viscosity" : 3.0,
"friction" : 1.0, "friction" : 1.0,
"lamda": -1.0, "lamda": -1,
"diffusion_rho" : 1.0,
"turnover_rho" : 1.0, "diffusion_rho" : 0.1,
"turnover_rho" : 1,
"average_rho" : 1.0, "average_rho" : 1.0,
"saturation_rho" : 1.0, "saturation_rho" : 1.0,
"diffusion_c" : 1.0, "stress_setpoint_rho": 2.0,
"turnover_c" : 1.0,
"average_c" : 1.0,
"noise_level": 0.0, "noise_level": 0.0,
"timestep" : 0.1, "timestep" : 0.01,
"savetime": 0.2, "savetime": 0.1,
"maxtime" : 100.0 "maxtime" : 100.0,
"q_xx" : 1.0,
"q_xy" : 0.0
} }
\ No newline at end of file
growing_domain/full_model/run_growing_domain.py 0 → 100644
import json, datetime
import os
from growing_domain import Growing_Domain
from viz_growing_domain import visualize
import argparse
import dolfin as df
import h5py
import numpy as np
parser = argparse.ArgumentParser()
parser.add_argument('-j','--jsonfile', help='data file',
default='parameters.json')
parser.add_argument('-t','--time', help='time to run', type=float)
args = parser.parse_args()
assert os.path.isfile(args.jsonfile), '%s file not found' % args.jsonfile
with open(args.jsonfile) as jsonFile:
parameters = json.load(jsonFile)
if args.jsonfile=='parameters.json':
extend = False
print('Fresh run...')
timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
parameters['timestamp'] = timestamp
parametersFile = timestamp + '_parameters.json'
initu = None
initv = None
initrho = None
initTime = 0.0
mesh = None
else:
extend = True
print('Extending run %s...' % parameters['timestamp'])
parametersFile = args.jsonfile
savesteps = round(parameters['maxtime']/parameters['savetime'])
# read geometry and topology at the last timepoint from h5 file
var = 'density'
h5 = h5py.File( '%s_%s.h5' % (parameters['timestamp'], var,), 'r+')
geometry = np.array(h5['%s/%s_%d/mesh/geometry'%(var,var,savesteps)])
topology = np.array(h5['%s/%s_%d/mesh/topology'%(var,var,savesteps)])
# create mesh with this geometry and topology
mesh = df.Mesh()
editor = df.MeshEditor()
editor.open(mesh,
type='triangle' if parameters['dimension']==2 else 'tetrahedron',
tdim=parameters['dimension'], gdim=parameters['dimension'])
editor.init_vertices(len(geometry))
editor.init_cells(len(topology))
for j in range(len(geometry)):
editor.add_vertex(j, geometry[j])
for j in range(len(topology)):
editor.add_cell(j, topology[j])
editor.close()
# Read field values at the last time point
SFS = df.FunctionSpace(mesh, 'P', 1)
VFS = df.VectorFunctionSpace(mesh, 'P', 1)
initu, initv, initrho = df.Function(VFS), df.Function(VFS), df.Function(SFS)
uFile = df.XDMFFile('%s_displacement.xdmf' % parameters['timestamp'])
vFile = df.XDMFFile('%s_velocity.xdmf' % parameters['timestamp'])
rhoFile = df.XDMFFile('%s_density.xdmf' % parameters['timestamp'])
uFile.read_checkpoint(initu, 'displacement', savesteps)
vFile.read_checkpoint(initv, 'velocity', savesteps)
rhoFile.read_checkpoint(initrho, 'density', savesteps)
uFile.close()
vFile.close()
rhoFile.close()
initTime = parameters['maxtime']
parameters['maxtime'] = args.time
tissue = Growing_Domain(parameters, mesh)
tissue.solve(initu, initv, initrho, initTime, extend)
if extend:
parameters['maxtime'] = initTime + parameters['maxtime']
with open(parametersFile, "w") as jsonFile:
json.dump(parameters, jsonFile, indent=4, sort_keys=True)
visualize(parameters, DIR='')
growing_domain/full_model/viz_growing_domain.py
...@@ -3,246 +3,156 @@ import numpy as np ...@@ -3,246 +3,156 @@ import numpy as np
import vedo as vd import vedo as vd
import os import os
import h5py import h5py
from matplotlib.widgets import Slider
import progressbar import progressbar
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from tempfile import TemporaryDirectory from tempfile import TemporaryDirectory
def get_data(params, DIR=''): def get_data(params, DIR=''):
savesteps = int(params['maxtime']/params['savetime']) savesteps = round(params['maxtime']/params['savetime'])
times = np.arange(savesteps+1) * params['savetime'] times = np.arange(savesteps+1) * params['savetime']
# Read mesh geometry from h5 file # Read mesh geometry from h5 file
var = 'density' var = 'density'
h5 = h5py.File(os.path.join(DIR, '%s_%s.h5' % (params['timestamp'], var)), "r") h5 = h5py.File(os.path.join(DIR, '%s_%s.h5' % (params['timestamp'], var)), "r")
# should be in the loop if remeshing
topology = np.array(h5['%s/%s_0/mesh/topology'%(var,var)]) # NOTE: geometry and topology are lists of length len(times)
# NOTE: geometry[k] is a numpy array of shape (Num of vertices, 2) for timestep k
# NOTE: topology[k] is a numpy array of shape (Num of cells, 3) for timestep k
topology = []
geometry = [] geometry = []
for i in range(len(times)): for i in range(len(times)):
geometry.append(np.array(h5['%s/%s_%d/mesh/geometry'%(var,var,i)])) geometry.append(np.array(h5['%s/%s_%d/mesh/geometry'%(var,var,i)]))
topology.append(np.array(h5['%s/%s_%d/mesh/topology'%(var,var,i)]))
h5.close() h5.close()
geometry = np.array(geometry)
if params['dimension']==1: # create a mesh at every timestep with this geometry and topology
geometry, zeros= np.dsplit(geometry, 2) meshes = []
print('Making the mesh..')
# create a mesh
if params['dimension']==1: for k in range(len(times)):
mesh = df.IntervalMesh(params['resolution'], 0, params['system_size'])
else:
mesh = df.Mesh() mesh = df.Mesh()
editor = df.MeshEditor() editor = df.MeshEditor()
editor.open(mesh, editor.open(mesh,
type='triangle' if params['dimension']==2 else 'tetrahedron', type='triangle' if params['dimension']==2 else 'tetrahedron',
tdim=params['dimension'], gdim=params['dimension']) tdim=params['dimension'], gdim=params['dimension'])
editor.init_vertices(len(geometry[0])) editor.init_vertices(len(geometry[k]))
editor.init_cells(len(topology)) editor.init_cells(len(topology[k]))
for i in range(len(geometry[0])): for j in range(len(geometry[k])):
editor.add_vertex(i, geometry[0][i]) editor.add_vertex(j, geometry[k][j])
for i in range(len(topology)): for j in range(len(topology[k])):
editor.add_cell(i, topology[i]) editor.add_cell(j, topology[k][j])
editor.close() editor.close()
meshes.append(mesh)
# Read concentration data
SFS = df.FunctionSpace(mesh, 'P', 1) print(meshes[0].num_vertices(), meshes[0].num_cells())
VFS = df.VectorFunctionSpace(mesh, 'P', 1) print(meshes[-1].num_vertices(), meshes[-1].num_cells())
ui, vi, rhoi = df.Function(VFS), df.Function(VFS), df.Function(SFS)
u = np.zeros((len(times), mesh.num_vertices(), params['dimension'])) # df.plot(meshes[0])
v = np.zeros_like(u) # plt.show()
rho = np.zeros((len(times), mesh.num_vertices())) # df.plot(meshes[1])
iso_stress = np.zeros((len(times), mesh.num_vertices())) # plt.show()
# df.plot(meshes[2])
# plt.show()
# df.plot(meshes[3])
# plt.show()
# df.plot(meshes[4])
# plt.show()
# df.plot(meshes[5])
# plt.show()
# df.plot(meshes[6])
# plt.show()
# df.plot(meshes[7])
# plt.show()
# df.plot(meshes[8])
# plt.show()
# Read field values
u = []
v = []
rho = []
uFile = df.XDMFFile(os.path.join(DIR, '%s_displacement.xdmf' % params['timestamp'])) uFile = df.XDMFFile(os.path.join(DIR, '%s_displacement.xdmf' % params['timestamp']))
vFile = df.XDMFFile(os.path.join(DIR, '%s_velocity.xdmf' % params['timestamp'])) vFile = df.XDMFFile(os.path.join(DIR, '%s_velocity.xdmf' % params['timestamp']))
rhoFile = df.XDMFFile(os.path.join(DIR, '%s_density.xdmf' % params['timestamp'])) rhoFile = df.XDMFFile(os.path.join(DIR, '%s_density.xdmf' % params['timestamp']))
if params['morphogen']: bar = progressbar.ProgressBar(maxval=savesteps)
ci = df.Function(SFS)
c = np.zeros_like(rho) for steps in bar(range(savesteps+1)):
cFile = df.XDMFFile(os.path.join(DIR, '%s_concentration.xdmf' % params['timestamp'])) SFS = df.FunctionSpace(meshes[steps], 'P', 1)
VFS = df.VectorFunctionSpace(meshes[steps], 'P', 1)
ui, vi, rhoi = df.Function(VFS), df.Function(VFS), df.Function(SFS)
for steps in progressbar.progressbar(range(savesteps+1)):
uFile.read_checkpoint(ui, 'displacement', steps) uFile.read_checkpoint(ui, 'displacement', steps)
vFile.read_checkpoint(vi, 'velocity', steps) vFile.read_checkpoint(vi, 'velocity', steps)
rhoFile.read_checkpoint(rhoi, 'density', steps) rhoFile.read_checkpoint(rhoi, 'density', steps)
u_vec = ui.compute_vertex_values(mesh) u_vec = ui.compute_vertex_values(meshes[steps])
u[steps] = u_vec.reshape(params['dimension'], int(u_vec.shape[0]/params['dimension'])).T u.append(u_vec.reshape(params['dimension'], int(u_vec.shape[0]/params['dimension'])).T)
v_vec = vi.compute_vertex_values(mesh)
v[steps] = v_vec.reshape(params['dimension'], int(v_vec.shape[0]/params['dimension'])).T
rho[steps] = rhoi.compute_vertex_values(mesh)
i_stress = (params['bulk_elasticity'] * df.div(ui) v_vec = vi.compute_vertex_values(meshes[steps])
+ params['bulk_viscosity']*df.div(vi) v.append(v_vec.reshape(params['dimension'], int(v_vec.shape[0]/params['dimension'])).T)
+ params['lamda'] * rhoi*(rhoi-params['average_rho'])/(rhoi + params['saturation_rho'])
)
i_stress = df.project(i_stress, SFS)
iso_stress[steps] = i_stress.compute_vertex_values(mesh)
if params['morphogen']: rho.append(rhoi.compute_vertex_values(meshes[steps]))
cFile.read_checkpoint(ci, 'concentration', steps)
c[steps] = ci.compute_vertex_values(mesh)
uFile.close() uFile.close()
vFile.close() vFile.close()
rhoFile.close() rhoFile.close()
if params['morphogen']: return (times, topology, geometry, u, v, rho)
cFile.close()
return (times, topology, geometry, u, v, rho, c, iso_stress)
else:
return (times, topology, geometry, u, v, rho, iso_stress)
def visualize(params, DIR='', offscreen=False): def visualize(params, DIR='', offscreen=False):
if params['morphogen']: (times, topology, geometry, u, v, rho) = get_data(params, DIR)
(times, topology, geometry, u, v, rho, c, iso_stress) = get_data(params, DIR)
else:
(times, topology, geometry, u, v, rho, iso_stress) = get_data(params, DIR)
radial_coordinate = np.linalg.norm(geometry, axis=2)
# print(geometry[-1,1,:] - geometry[-1,0,:])
# print(geometry[-1,1,:]+ v[-2, 1,:] * params['timestep'])
# print(geometry[-1,0,:]+ v[-2, 0,:] * params['timestep'])
##################################################### 1-D ############################################################################################################
if params['dimension']==1:
# c(x,t) and rho(x,t) vs x and t
if params['morphogen']:
fig, (axrho, axc) = plt.subplots(2,1, sharex=True, figsize=(8,8))
axc.set_xlabel(r'$x$')
axc.set_ylabel(r"$c(x,t)$")
axc.set_xlim(np.min(radial_coordinate), np.max(radial_coordinate))
axc.set_ylim(np.min(c), np.max(c))
else:
fig, axrho = plt.subplots(1,1, sharex=True, figsize=(8,8))
axrho.set_xlabel(r'$x$')
axrho.set_ylabel(r"$\rho(x,t)$")
axrho.set_xlim(np.min(radial_coordinate), np.max(radial_coordinate))
axrho.set_ylim(np.min(rho), np.max(rho))
rhoplot, = axrho.plot(radial_coordinate[0], rho[0],'o',ms=3)
if params['morphogen']:
cplot, = axc.plot(radial_coordinate[0], c[0])
def update(value):
ti = np.abs(times-value).argmin()
rhoplot.set_ydata(rho[ti])
rhoplot.set_xdata(radial_coordinate[ti])
if params['morphogen']:
cplot.set_ydata(c[ti])
cplot.set_xdata(radial_coordinate[ti])
plt.draw()
sax = plt.axes([0.1, 0.92, 0.7, 0.02])
slider = Slider(sax, r'$t/\tau$', min(times), max(times),
valinit=min(times), valfmt='%3.1f',
fc='#999999')
slider.drawon = False
slider.on_changed(update)
plt.show()
# heat map
# L(t) vs t
length = np.array(np.zeros(len(times)))
for j in range(len(times)):
length[j] = np.max(radial_coordinate[j])
figlength, axlength = plt.subplots(1, 1,)
axlength.set_xlabel('$t$')
axlength.set_ylabel('$L(t)$')
axlength.set_xlim(np.min(times), np.max(times))
axlength.set_ylim(np.min(length), np.max(length))
axlength.plot(times, length)
plt.show()
###################################### 2D ##############################################################################################################
else:
# rho and v
n_cmap_vals = 16 n_cmap_vals = 16
scalar_cmap = 'viridis' scalar_cmap = 'viridis'
vector_color = 'w' rhomin, rhomax = min(min(sublist) for sublist in rho), max(max(sublist) for sublist in rho)
vector_scale= 0.1
geometry = np.dstack((geometry, np.zeros(geometry.shape[0:2])))
v = np.dstack((v, np.zeros(v.shape[0:2])))
vmag = np.linalg.norm(v, axis=2)
vmax = np.max(vmag)
v = v/vmax
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) plotter = vd.plotter.Plotter(axes=0)
poly = vd.utils.buildPolyData(geometry[0], topology) poly = vd.utils.buildPolyData(geometry[0], topology[0])
scalar_actor = vd.mesh.Mesh(poly) scalar_actor = vd.mesh.Mesh(poly)
#scalar_actor.computeNormals(points=True, cells=True)
scalar_actor.pointdata['density'] = rho[0] scalar_actor.pointdata['density'] = rho[0]
scalar_actor.cmap(scalar_cmap, rho[0], vmin=rhomin, vmax=rhomax, n=n_cmap_vals) scalar_actor.cmap(scalar_cmap, rho[0], vmin=rhomin, vmax=rhomax)
scalar_actor.add_scalarbar(title = r'$\rho/\rho_0$', scalar_actor.add_scalarbar(title = r'$\rho/\rho_0$',
pos=(0.8, 0.04), nlabels=2, pos=(0.8, 0.04), nlabels=2,
# titleYOffset=15, titleFontSize=28, size=(100, 600) # titleYOffset=15, titleFontSize=28, size=(100, 600)
) )
plotter += scalar_actor plotter += scalar_actor
# Create a wireframe of the scalar_actor and add it to the plotter
vector_actor = vd.shapes.Arrows(geometry[0], geometry[0]+ vector_scale * v[0]) wireframe_actor = scalar_actor.clone().wireframe(True).lw(0.1).c('black')
vector_actor.color(vector_color) plotter += wireframe_actor
plotter+=vector_actor
def update(idx): def update(idx):
nonlocal plotter, vector_actor nonlocal plotter, scalar_actor, wireframe_actor
scalar_actor.points(pts=geometry[idx], transformed=False) poly = vd.utils.buildPolyData(geometry[idx], topology[idx])
scalar_actor.pointdata['density'] = rho[idx] new_scalar_actor = vd.mesh.Mesh(poly)
plotter.remove(vector_actor) new_scalar_actor.pointdata['density'] = rho[idx]
vector_actor = vd.shapes.Arrows(geometry[idx], geometry[idx]+ vector_scale * v[idx])
vector_actor.color(vector_color)
plotter += vector_actor
def slider_update(widget, event): new_scalar_actor.cmap(scalar_cmap, rho[idx], vmin=rhomin, vmax=rhomax)
value = widget.GetRepresentation().GetValue() new_scalar_actor.add_scalarbar(title = r'$\rho/\rho_0$',
idx = (abs(times-value)).argmin()
update(idx)
slider = plotter.add_slider(slider_update, pos=[(0.1,0.94), (0.5,0.94)],
xmin=times[0], xmax=times.max(),
value=times[0], title=r"$t/\tau$")
vd.show([scalar_actor, vector_actor],interactive=(not offscreen), zoom=0.8)
# isotropic stress
iso_stressmin, iso_stressmax = np.min(iso_stress), np.max(iso_stress)
plotter_i_stress = vd.plotter.Plotter(axes=0)
poly_i_stress = vd.utils.buildPolyData(geometry[0], topology)
scalar_actor_i_stress = vd.mesh.Mesh(poly_i_stress)
#scalar_actor.computeNormals(points=True, cells=True)
scalar_actor_i_stress.pointdata['iso_stress'] = iso_stress[0]
scalar_actor_i_stress.cmap(scalar_cmap, iso_stress[0], vmin=iso_stressmin, vmax=iso_stressmax, n=n_cmap_vals)
scalar_actor_i_stress.add_scalarbar(title = r'$\rho/\rho_0$',
pos=(0.8, 0.04), nlabels=2, pos=(0.8, 0.04), nlabels=2,
# titleYOffset=15, titleFontSize=28, size=(100, 600) # titleYOffset=15, titleFontSize=28, size=(100, 600)
) )
plotter_i_stress += scalar_actor_i_stress plotter.remove(scalar_actor)
plotter.remove(wireframe_actor)
plotter += new_scalar_actor
new_wireframe_actor = new_scalar_actor.clone().wireframe(True).lw(0.1).c('black')
plotter += new_wireframe_actor
def update_i_stress(idx): scalar_actor = new_scalar_actor
scalar_actor_i_stress.points(pts=geometry[idx], transformed=False) wireframe_actor = new_wireframe_actor # Update the reference to the wireframe_actor
scalar_actor_i_stress.pointdata['iso_stress'] = iso_stress[idx] # plotter.render()
def slider_update_i_stress(widget, event): def slider_update(widget, event):
value = widget.GetRepresentation().GetValue() value = widget.GetRepresentation().GetValue()
idx = (abs(times-value)).argmin() idx = (abs(times-value)).argmin()
update_i_stress(idx) update(idx)
slider_i_stress = plotter_i_stress.add_slider(slider_update_i_stress, pos=[(0.1,0.94), (0.5,0.94)], slider = plotter.add_slider(slider_update, pos=[(0.1,0.94), (0.5,0.94)],
xmin=times[0], xmax=times.max(), xmin=times[0], xmax=times.max(),
value=times[0], title=r"$t/\tau$") value=times[0], title=r"$t/\tau$")
vd.show([scalar_actor], interactive=(not offscreen), zoom=0.5)
vd.show(scalar_actor_i_stress,interactive=(not offscreen), zoom=0.8)
if offscreen: if offscreen:
FPS = 5 FPS = 5
...@@ -257,124 +167,29 @@ def visualize(params, DIR='', offscreen=False): ...@@ -257,124 +167,29 @@ def visualize(params, DIR='', offscreen=False):
update(idx) update(idx)
slider.GetRepresentation().SetValue(times[tt]) slider.GetRepresentation().SetValue(times[tt])
fr = get_filename("%03d.png" % tt) fr = get_filename("%03d.png" % tt)
vd.io.screenshot(fr) vd.screenshot(fr)
os.system(command + " " + str(fps) os.system(command + " " + str(fps)
+ " -i " + tmp_dir.name + os.sep + " -i " + tmp_dir.name + os.sep
+ "%03d.png " + options + " " + movFile) + "%03d.png " + options + " " + movFile)
tmp_dir.cleanup() tmp_dir.cleanup()
if params['morphogen']: # radial coordinate
plotter1 = vd.plotter.Plotter(axes=0) r = np.zeros(len(times))
for i in range(len(times)):
poly1 = vd.utils.buildPolyData(geometry[0], topology) r[i] = np.max(np.sqrt(np.sum(np.square(geometry[i]), axis=1)))
scalar_actor1 = vd.mesh.Mesh(poly1)
#scalar_actor1.computeNormals(points=True, cells=True)
scalar_actor1.pointdata['concentration'] = c[0]
scalar_actor1.cmap('plasma', 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: plt.plot(times, r)
FPS1 = 5 plt.xlabel('time')
movFile1 = '%s_morphogen.mov' % params['timestamp'] plt.ylabel('radius')
fps1 = float(FPS1)
command1 = "ffmpeg -y -r"
options1 = "-b:v 3600k -qscale:v 4 -vcodec mpeg4"
tmp_dir1 = TemporaryDirectory()
get_filename1 = lambda x: os.path.join(tmp_dir1.name, x)
for tt in range(len(times)):
idx1 = (abs(times-times[tt])).argmin()
update1(idx1)
slider1.GetRepresentation().SetValue(times[tt])
fr1 = get_filename1("%03d.png" % tt)
vd.io.screenshot(fr1)
os.system(command1+ " " + str(fps1)
+ " -i " + tmp_dir1.name + os.sep
+ "%03d.png " + options1 + " " + movFile1)
tmp_dir1.cleanup()
# R(t) vs t
plt.rcParams.update({'font.size': 22})
radius = np.array(np.zeros(len(times)))
for j in range(len(times)):
radius[j] = np.max(radial_coordinate[j])
figradius3, axradius3 = plt.subplots(1,1)
axradius3.set_xlabel(r'$\log t$')
# axradius3.set_xlim(np.min(np.log(times)), np.max(np.log(times)))
# axradius3.set_ylim(np.min(np.log(radius)), np.max(np.log(radius)))
axradius3.set_ylabel(r'$\log R(t)$')
axradius3.loglog(times, radius)
figradius1, axradius1 = plt.subplots(1,1)
axradius1.set_xlabel(r'$t$')
# axradius1.set_xlim(np.min(times), np.max(times))
# axradius1.set_ylim(np.min(np.log(radius)), np.max(np.log(radius)))
axradius1.set_ylabel(r'$\log R(t)$')
axradius1.semilogy(times, radius)
figradius2, axradius2 = plt.subplots(1,1)
axradius2.set_xlabel(r'$t$')
axradius2.set_xlim(np.min(times), np.max(times))
axradius2.set_ylim(np.min(radius), np.max(radius))
axradius2.set_ylabel(r'$R(t)$')
axradius2.plot(times, radius)
plt.show() plt.show()
# c(r,t) vs r drdt = np.gradient(r, times)
if params['morphogen']: print(drdt)
figradial, (axcradial, axrhoradial) = plt.subplots(2,1, figsize=(8,8)) plt.plot(times, drdt)
axcradial.set_xlabel(r'$r$') plt.xlabel('time')
axcradial.set_xlim(np.min(radial_coordinate), np.max(radial_coordinate)) plt.ylabel('dr/dt')
axcradial.set_ylim(np.min(c), np.max(c))
axcradial.set_ylabel(r'$c(r,t)$')
cplot, = axcradial.plot(radial_coordinate[0], c[0],'o', ms=3)
else:
figradial, axrhoradial = plt.subplots(1,1, figsize=(8,8))
axrhoradial.set_xlabel(r'$r$')
axrhoradial.set_xlim(np.min(radial_coordinate), np.max(radial_coordinate))
axrhoradial.set_ylim(np.min(rho), np.max(rho))
axrhoradial.set_ylabel(r'$\rho(r,t)$')
rhoplot, = axrhoradial.plot(radial_coordinate[0], rho[0], 'o', ms=3)
def update(value):
ti = np.abs(times-value).argmin()
rhoplot.set_ydata(rho[ti])
rhoplot.set_xdata(radial_coordinate[ti])
if params['morphogen']:
cplot.set_ydata(c[ti])
cplot.set_xdata(radial_coordinate[ti])
plt.draw()
sax = plt.axes([0.1, 0.92, 0.7, 0.02])
slider1 = Slider(sax, r'$t/\tau$', min(times), max(times),
valinit=min(times), valfmt='%3.1f',
fc='#999999')
slider1.drawon = False
slider1.on_changed(update)
plt.show() plt.show()
if __name__ == '__main__': if __name__ == '__main__':
import argparse, json import argparse, json
parser = argparse.ArgumentParser() parser = argparse.ArgumentParser()
......
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