Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit 162e68e0 by Jigyasa Watwani
Options

visualisation for isotropic stress / areal expansion

parent 2982497b
Inline Side-by-side
Showing with 58 additions and 6 deletions
growing_domain/full_model/viz_growing_domain.py
...@@ -48,6 +48,8 @@ def get_data(params, DIR=''): ...@@ -48,6 +48,8 @@ def get_data(params, DIR=''):
u = np.zeros((len(times), mesh.num_vertices(), params['dimension'])) u = np.zeros((len(times), mesh.num_vertices(), params['dimension']))
v = np.zeros_like(u) v = np.zeros_like(u)
rho = np.zeros((len(times), mesh.num_vertices())) rho = np.zeros((len(times), mesh.num_vertices()))
iso_stress = np.zeros((len(times), mesh.num_vertices()))
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']))
...@@ -55,6 +57,7 @@ def get_data(params, DIR=''): ...@@ -55,6 +57,7 @@ def get_data(params, DIR=''):
ci = df.Function(SFS) ci = df.Function(SFS)
c = np.zeros_like(rho) c = np.zeros_like(rho)
cFile = df.XDMFFile(os.path.join(DIR, '%s_concentration.xdmf' % params['timestamp'])) cFile = df.XDMFFile(os.path.join(DIR, '%s_concentration.xdmf' % params['timestamp']))
for steps in progressbar.progressbar(range(savesteps+1)): 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)
...@@ -67,28 +70,43 @@ def get_data(params, DIR=''): ...@@ -67,28 +70,43 @@ def get_data(params, DIR=''):
v[steps] = v_vec.reshape(params['dimension'], int(v_vec.shape[0]/params['dimension'])).T v[steps] = v_vec.reshape(params['dimension'], int(v_vec.shape[0]/params['dimension'])).T
rho[steps] = rhoi.compute_vertex_values(mesh) rho[steps] = rhoi.compute_vertex_values(mesh)
i_stress = (params['bulk_elasticity'] * df.div(ui)
+ params['bulk_viscosity']*df.div(vi)
+ 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']: if params['morphogen']:
cFile.read_checkpoint(ci, 'concentration', steps) cFile.read_checkpoint(ci, 'concentration', steps)
c[steps] = ci.compute_vertex_values(mesh) c[steps] = ci.compute_vertex_values(mesh)
uFile.close() uFile.close()
vFile.close() vFile.close()
rhoFile.close() rhoFile.close()
if params['morphogen']: if params['morphogen']:
cFile.close() cFile.close()
return (times, topology, geometry, u, v, rho, c) return (times, topology, geometry, u, v, rho, c, iso_stress)
else: else:
return (times, topology, geometry, u, v, rho) return (times, topology, geometry, u, v, rho, iso_stress)
def visualize(params, DIR='', offscreen=False): def visualize(params, DIR='', offscreen=False):
if params['morphogen']: if params['morphogen']:
(times, topology, geometry, u, v, rho, c) = get_data(params, DIR) (times, topology, geometry, u, v, rho, c, iso_stress) = get_data(params, DIR)
else: else:
(times, topology, geometry, u, v, rho) = get_data(params, DIR) (times, topology, geometry, u, v, rho, iso_stress) = get_data(params, DIR)
radial_coordinate = np.linalg.norm(geometry, axis=2) 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: if params['dimension']==1:
# c(x,t) and rho(x,t) vs x and t # c(x,t) and rho(x,t) vs x and t
if params['morphogen']: if params['morphogen']:
...@@ -141,8 +159,10 @@ def visualize(params, DIR='', offscreen=False): ...@@ -141,8 +159,10 @@ def visualize(params, DIR='', offscreen=False):
axlength.plot(times, length) axlength.plot(times, length)
plt.show() plt.show()
###################################### 2D ##############################################################################################################
else: else:
# heat map # rho and v
n_cmap_vals = 16 n_cmap_vals = 16
scalar_cmap = 'viridis' scalar_cmap = 'viridis'
vector_color = 'w' vector_color = 'w'
...@@ -172,6 +192,7 @@ def visualize(params, DIR='', offscreen=False): ...@@ -172,6 +192,7 @@ def visualize(params, DIR='', offscreen=False):
vector_actor = vd.shapes.Arrows(geometry[0], geometry[0]+ vector_scale * v[0]) vector_actor = vd.shapes.Arrows(geometry[0], geometry[0]+ vector_scale * v[0])
vector_actor.color(vector_color) vector_actor.color(vector_color)
plotter+=vector_actor plotter+=vector_actor
def update(idx): def update(idx):
nonlocal plotter, vector_actor nonlocal plotter, vector_actor
scalar_actor.points(pts=geometry[idx], transformed=False) scalar_actor.points(pts=geometry[idx], transformed=False)
...@@ -181,7 +202,6 @@ def visualize(params, DIR='', offscreen=False): ...@@ -181,7 +202,6 @@ def visualize(params, DIR='', offscreen=False):
vector_actor.color(vector_color) vector_actor.color(vector_color)
plotter += vector_actor plotter += vector_actor
def slider_update(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()
...@@ -192,6 +212,38 @@ def visualize(params, DIR='', offscreen=False): ...@@ -192,6 +212,38 @@ def visualize(params, DIR='', offscreen=False):
value=times[0], title=r"$t/\tau$") value=times[0], title=r"$t/\tau$")
vd.show([scalar_actor, vector_actor],interactive=(not offscreen), zoom=0.8) 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,
# titleYOffset=15, titleFontSize=28, size=(100, 600)
)
plotter_i_stress += scalar_actor_i_stress
def update_i_stress(idx):
scalar_actor_i_stress.points(pts=geometry[idx], transformed=False)
scalar_actor_i_stress.pointdata['iso_stress'] = iso_stress[idx]
def slider_update_i_stress(widget, event):
value = widget.GetRepresentation().GetValue()
idx = (abs(times-value)).argmin()
update_i_stress(idx)
slider_i_stress = plotter_i_stress.add_slider(slider_update_i_stress, 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_i_stress,interactive=(not offscreen), zoom=0.8)
if offscreen: if offscreen:
FPS = 5 FPS = 5
movFile = '%s_cell_density.mov' % params['timestamp'] movFile = '%s_cell_density.mov' % params['timestamp']
......
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