Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit a8b9a7bb by Jigyasa Watwani
Options

fully working 1-D problem with remeshing and extend functionality

parent 3a76deec
Inline Side-by-side
Showing with 206 additions and 141 deletions
basic_model_only_density/parameters.json
......@@ -15,6 +15,6 @@
"noise_level": 0.0,
"timestep": 0.01,
"savetime": 0.1,
"maxtime": 100
"maxtime": 10
}
\ No newline at end of file
basic_model_only_density/run_tissue.py
......@@ -35,18 +35,29 @@ else:
parametersFile = args.jsonfile
savesteps = round(parameters['maxtime']/parameters['savetime'])
mesh = df.IntervalMesh(parameters['resolution'],
- parameters['system_size']/2,
parameters['system_size']/2)
# read geometry from h5 file
# 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)])
mesh.coordinates()[:,0] = geometry[:,0]
geometry = np.array(h5['%s/%s_%d/mesh/geometry'%(var,var,savesteps)])[:, 0]
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, 'interval', 1, 1)
# initialise vertices and cells
editor.init_vertices(len(geometry))
editor.init_cells(len(topology))
# add vertices
for j in range(len(geometry)):
editor.add_vertex(j, [geometry[j]])
# add cells
for j in range(len(topology)):
editor.add_cell(j, topology[j])
editor.close()
# Read data
# 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)
......@@ -71,5 +82,4 @@ if extend:
with open(parametersFile, "w") as jsonFile:
json.dump(parameters, jsonFile, indent=4, sort_keys=True)
from viz_tissue import visualize
visualize(parameters, DIR='')
\ No newline at end of file
basic_model_only_density/tissue.py
......@@ -13,22 +13,20 @@ class Tissue(object):
# if no mesh is given as initial condition, create one
if mesh is None:
self.mesh = df.IntervalMesh(self.resolution,
-self.system_size/2,
self.system_size/2)
self.mesh = df.IntervalMesh(self.resolution, -self.system_size/2, self.system_size/2)
else:
self.mesh = mesh
def setup(self):
scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
vector_element = df.VectorElement('P', self.mesh.ufl_cell(), 1)
# u, v, rho
mixed_element = df.MixedElement([vector_element,
vector_element,
scalar_element])
self.mixed_element = df.MixedElement([vector_element, vector_element, scalar_element])
# define function space with this mixed element
self.function_space = df.FunctionSpace(self.mesh, mixed_element)
self.function_space = df.FunctionSpace(self.mesh, self.mixed_element)
# boundary condition on rho
if self.zero_rho_boundary:
......@@ -43,7 +41,6 @@ class Tissue(object):
# define functions on function space
self.function0 = df.Function(self.function_space)
self.function = df.Function(self.function_space)
def advection(self, conc, vel, tconc):
return df.inner(df.div(vel * conc), tconc)
......@@ -56,10 +53,18 @@ class Tissue(object):
return (-self.lamda * rho /(rho + self.saturation_rho)
* df.Identity(1))
def epsilon(self, v):
return df.sym(df.nabla_grad(v))
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.h() > 0.1: # if the width of the cell> 0.1, refine
cell_markers[cell] = True
mesh = df.refine(mesh, cell_markers)
return mesh
def passive_stress(self, u, v):
eps_u, eps_v = self.epsilon(u), self.epsilon(v)
elastic_stress = self.elasticity * eps_u
......@@ -77,10 +82,11 @@ class Tissue(object):
)
def setup_initial_conditions(self):
# ic for u
zero_vector = df.Constant((0.0,))
u0 = df.interpolate(zero_vector,
self.function_space.sub(0).collapse())
# ic for rho
if self.zero_rho_boundary:
base_rho = 0.0
else:
......@@ -122,63 +128,83 @@ class Tissue(object):
+ self.diffusion_reaction_rho(rho, trho))
self.form = (uform + vform + rhoform) * df.dx
def solve(self, initu=None, initv=None, initrho=None, initTime=0, extend=False):
# create function space, functions, bc
self.setup()
# create files if they don't exist, open them if they do 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)
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
self.time = initTime
savesteps = round(self.savetime/self.timestep)
maxsteps = round(self.maxtime/self.timestep)
# if not extend, find the initial velocity, save the ic, move the mesh with this initial velocity
if not extend:
# find the initial velocity
self.setup_initial_conditions()
# save the ic
u0, v0, rho0 = self.function0.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)
# if extend, move the mesh with this initial velocity
else:
# assign fields to func0
u0 = df.project(initu, self.function_space.sub(0).collapse())
v0 = df.project(initv, self.function_space.sub(1).collapse())
rho0 = df.project(initrho, self.function_space.sub(2).collapse())
df.assign(self.function0, [u0, v0, rho0])
# move the mesh with initial v read from file
dr0 = df.project(v0 * self.timestep, self.function_space.sub(0).collapse())
df.ALE.move(self.mesh, dr0)
df.assign(self.function0, [u0, v0, rho0])
self.setup_weak_forms()
for steps in progressbar.progressbar(range(1, maxsteps + 1)):
self.time += self.timestep
# solve to get fields at next timestep
df.solve(self.form == 0, self.function, self.bc)
self.function0.assign(self.function)
u, v, rho = self.function0.split(deepcopy=True)
# save the fields
u, v, rho = self.function.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.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)
# assign the calculated function to a newly-defined function, to be used later
old_function = df.Function(self.function_space)
old_function.assign(self.function)
# move mesh
# move the mesh
dr = df.project(v * self.timestep, self.function_space.sub(0).collapse())
df.ALE.move(self.mesh, dr)
# refine the mesh
self.mesh = self.refine_mesh(self.mesh)
# new function space, bc, functions
self.setup()
# new function0 should be the old function in the old function space projected on the new function space
self.function0 = df.project(old_function, self.function_space)
# new weak form
self.setup_weak_forms()
self.uFile.close()
self.vFile.close()
self.rhoFile.close()
basic_model_only_density/viz_tissue.py
......@@ -11,89 +11,116 @@ def visualize(params, DIR=''):
savesteps = round(params['maxtime']/params['savetime'])
times = np.arange(savesteps+1) * params['savetime']
# Read mesh geometry from h5 file
# Read mesh geometry and topology from h5 file
var = 'density'
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][:, 0] is a numpy array of shape (Num of vertices, ) for timestep k
# NOTE: topology[k] is a numpy array of shape (Num of cells, 2) for timestep k
geometry = []
topology = []
print('Reading geometry and topology of the mesh..')
for i in progressbar.progressbar(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)])[:, 0])
topology.append(np.array(h5['%s/%s_%d/mesh/topology'%(var,var,i)]))
h5.close()
geometry = np.array(geometry)
geometry, zeros = np.dsplit(geometry, 2)
mesh = df.IntervalMesh(params['resolution'],
- params['system_size']/2,
params['system_size']/2)
# Read data
u = np.zeros((len(times), mesh.num_vertices(), 1))
v = np.zeros_like(u)
rho = np.zeros((len(times), mesh.num_vertices()))
elastic_stress = np.zeros((len(times), mesh.num_vertices()))
viscous_stress = np.zeros_like(elastic_stress)
active_stress = np.zeros_like(elastic_stress)
uFile = df.XDMFFile(os.path.join(DIR,
'%s_displacement.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']))
# create mesh at every timestep with this geometry and topology
meshes = []
print('Making the mesh..')
for k in progressbar.progressbar(range(len(times))):
# initialise mesh and mesh editor
mesh = df.Mesh()
editor = df.MeshEditor()
editor.open(mesh, 'interval', 1, 1)
# initialise vertices and cells
editor.init_vertices(len(geometry[k]))
editor.init_cells(len(topology[k]))
# add vertices
for j in range(len(geometry[k])):
editor.add_vertex(j, [geometry[k][j]])
# add cells
for j in range(len(topology[k])):
editor.add_cell(j, topology[k][j])
editor.close()
meshes.append(mesh)
# Read field values
u = []
v = []
rho = []
elastic_stress = []
viscous_stress = []
active_stress = []
uFile = df.XDMFFile(os.path.join(DIR, '%s_displacement.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']))
# Reading data
print('Reading data...')
print('Reading the fields..')
for steps in progressbar.progressbar(range(savesteps+1)):
mesh.coordinates()[:] = geometry[steps]
SFS = df.FunctionSpace(mesh, 'P', 1)
VFS = df.VectorFunctionSpace(mesh, 'P', 1)
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)
uFile.read_checkpoint(ui, 'displacement', steps)
vFile.read_checkpoint(vi, 'velocity', steps)
rhoFile.read_checkpoint(rhoi, 'density', steps)
u_vec = ui.compute_vertex_values(mesh)
u[steps] = u_vec.reshape(1, int(u_vec.shape[0])).T
v_vec = vi.compute_vertex_values(mesh)
v[steps] = v_vec.reshape(1, int(v_vec.shape[0])).T
rho[steps] = rhoi.compute_vertex_values(mesh)
u.append(ui.compute_vertex_values(meshes[steps]))
v.append(vi.compute_vertex_values(meshes[steps]))
rho.append(rhoi.compute_vertex_values(meshes[steps]))
e_stress = params['elasticity'] * ui.dx(0)
e_stress = df.project(e_stress, VFS)
elastic_stress[steps] = e_stress.compute_vertex_values(mesh)
elastic_stress.append(e_stress.compute_vertex_values(meshes[steps]))
v_stress = params['viscosity'] * vi.dx(0)
v_stress = df.project(v_stress, VFS)
viscous_stress[steps] = v_stress.compute_vertex_values(mesh)
viscous_stress.append(v_stress.compute_vertex_values(meshes[steps]))
a_stress = -params['lamda'] * rhoi/(rhoi + params['saturation_rho'])
a_stress = df.project(a_stress, SFS)
active_stress[steps] = a_stress.compute_vertex_values(mesh)
active_stress.append(a_stress.compute_vertex_values(meshes[steps]))
uFile.close()
vFile.close()
rhoFile.close()
# minima and maxima of the fields
min_rho, min_u, min_v = map(lambda x: min(min(sublist) for sublist in x), [rho, u, v])
max_rho, max_u, max_v = map(lambda x: max(max(sublist) for sublist in x), [rho, u, v])
min_estress, min_vstress, min_astress = (map(lambda x: min(min(sublist) for sublist in x),
[elastic_stress, viscous_stress, active_stress])
)
max_estress, max_vstress, max_astress = (map(lambda x: max(max(sublist) for sublist in x),
[elastic_stress, viscous_stress, active_stress])
)
min_geometry, max_geometry = min(min(sublist) for sublist in geometry), max(max(sublist) for sublist in geometry)
# interactive plot
plt.rcParams.update({'font.size': 15})
fig, axes = plt.subplots(3, 1, sharex=True, figsize=(8,8))
axes[-1].set_xlabel(r'$x$')
axes[0].set_ylabel(r'$\rho/\rho_0$')
axes[1].set_ylabel(r'$v/v_0$')
axes[2].set_ylabel(r'$u/u_0$')
axes[0].set_xlim(np.min(geometry), np.max(geometry))
axes[0].set_ylim(np.min(rho)-0.1, np.max(rho)+0.1)
axes[1].set_ylim(np.min(v)-0.1, np.max(v)+0.1)
axes[2].set_ylim(np.min(u)-0.1, np.max(u)+0.1)
axes[0].set_xlim(min_geometry, max_geometry)
axes[0].set_ylim(min_rho-0.1, max_rho+0.1)
axes[1].set_ylim(min_v-0.1, max_v+0.1)
axes[2].set_ylim(min_u-0.1, max_u+0.1)
rhoplot, = axes[0].plot(geometry[0], rho[0], 'g-', ms=3)
velplot, = axes[1].plot(geometry[0], v[0], 'r-', ms=3)
uplot, = axes[2].plot(geometry[0], u[0], 'b-', ms=3)
rhoplot, = axes[0].plot(geometry[0], rho[0], ms=3, color='g', marker='.', linestyle='None')
velplot, = axes[1].plot(geometry[0], v[0], ms=3, color='r', marker='.', linestyle='None')
uplot, = axes[2].plot(geometry[0], u[0], ms=3, color='b', marker='.', linestyle='None')
def update(value):
ti = np.abs(times-value).argmin()
......@@ -113,22 +140,22 @@ def visualize(params, DIR=''):
slider.on_changed(update)
plt.show()
# print('Saving movie-...')
# FPS = 50
# movFile = os.path.join(DIR, '%s_fields.mov' % params['timestamp'])
# fps = float(FPS)
# command = "ffmpeg -y -r"
# options = "-b:v 3600k -qscale:v 4 -vcodec mpeg4"
# tmp_dir = TemporaryDirectory()
# get_filename = lambda x: os.path.join(tmp_dir.name, x)
# for tt in progressbar.progressbar(range(len(times))):
# slider.set_val(times[tt])
# fr = get_filename("%03d.png" % tt)
# fig.savefig(fr, facecolor=fig.get_facecolor(), dpi=100)
# os.system(command + " " + str(fps)
# + " -i " + tmp_dir.name + os.sep
# + "%03d.png " + options + " " + movFile)
# tmp_dir.cleanup()
print('Saving movie for fields..')
FPS = 50
movFile = os.path.join(DIR, '%s_fields.mov' % params['timestamp'])
fps = float(FPS)
command = "ffmpeg -y -r"
options = "-b:v 3600k -qscale:v 4 -vcodec mpeg4"
tmp_dir = TemporaryDirectory()
get_filename = lambda x: os.path.join(tmp_dir.name, x)
for tt in progressbar.progressbar(range(len(times))):
slider.set_val(times[tt])
fr = get_filename("%03d.png" % tt)
fig.savefig(fr, facecolor=fig.get_facecolor(), dpi=100)
os.system(command + " " + str(fps)
+ " -i " + tmp_dir.name + os.sep
+ "%03d.png " + options + " " + movFile)
tmp_dir.cleanup()
# plotting the stresses
figure, axes1 = plt.subplots(3, 1, sharex=True, figsize=(8,8))
......@@ -136,14 +163,14 @@ def visualize(params, DIR=''):
axes1[2].set_ylabel('Elastic stress')
axes1[1].set_ylabel('Viscous stress')
axes1[0].set_ylabel('Active stress')
axes1[0].set_xlim(np.min(geometry), np.max(geometry))
axes1[2].set_ylim(np.min(elastic_stress)-0.1, np.max(elastic_stress)+0.1)
axes1[1].set_ylim(np.min(viscous_stress)-0.1, np.max(viscous_stress)+0.1)
axes1[0].set_ylim(np.min(active_stress)-0.1, np.max(active_stress)+0.1)
axes1[0].set_xlim(min_geometry, max_geometry)
axes1[2].set_ylim(min_estress-0.1, max_estress+0.1)
axes1[1].set_ylim(min_vstress-0.1, max_vstress+0.1)
axes1[0].set_ylim(min_astress-0.1, max_astress+0.1)
e_stress_plot, = axes1[2].plot(geometry[0], elastic_stress[0], 'g-', ms=3)
v_stress_plot, = axes1[1].plot(geometry[0], viscous_stress[0], 'r-', ms=3)
a_stress_plot, = axes1[0].plot(geometry[0], -active_stress[0], 'b-', ms=3)
e_stress_plot, = axes1[2].plot(geometry[0], elastic_stress[0], ms=3, color='b', marker='.', linestyle='None')
v_stress_plot, = axes1[1].plot(geometry[0], viscous_stress[0], ms=3, color='r', marker='.', linestyle='None')
a_stress_plot, = axes1[0].plot(geometry[0], active_stress[0], ms=3, color='g', marker='.', linestyle='None')
def update2(value):
ti = np.abs(times-value).argmin()
......@@ -167,34 +194,42 @@ def visualize(params, DIR=''):
plt.show()
# make movie
# print('Saving movie-...')
# FPS = 50
# movFile = os.path.join(DIR, '%s_stresses.mov' % params['timestamp'])
# fps = float(FPS)
# command = "ffmpeg -y -r"
# options = "-b:v 3600k -qscale:v 4 -vcodec mpeg4"
# tmp_dir = TemporaryDirectory()
# get_filename = lambda x: os.path.join(tmp_dir.name, x)
# for tt in progressbar.progressbar(range(len(times))):
# slider.set_val(times[tt])
# fr = get_filename("%03d.png" % tt)
# figure.savefig(fr, facecolor=figure.get_facecolor(), dpi=100)
# os.system(command + " " + str(fps)
# + " -i " + tmp_dir.name + os.sep
# + "%03d.png " + options + " " + movFile)
# tmp_dir.cleanup()
print('Saving movie for stresses..')
FPS = 50
movFile = os.path.join(DIR, '%s_stresses.mov' % params['timestamp'])
fps = float(FPS)
command = "ffmpeg -y -r"
options = "-b:v 3600k -qscale:v 4 -vcodec mpeg4"
tmp_dir = TemporaryDirectory()
get_filename = lambda x: os.path.join(tmp_dir.name, x)
for tt in progressbar.progressbar(range(len(times))):
slider.set_val(times[tt])
fr = get_filename("%03d.png" % tt)
figure.savefig(fr, facecolor=figure.get_facecolor(), dpi=100)
os.system(command + " " + str(fps)
+ " -i " + tmp_dir.name + os.sep
+ "%03d.png " + options + " " + movFile)
tmp_dir.cleanup()
# L(t) vs t
length = geometry[:,-1,0]-geometry[:,0,0]
# print(length[-1])
length = np.zeros(len(times))
# NOTE: Since after remeshing, the new points are added at the end of the array, it doesn't make
# NOTE: sense to find the length by subtracting the coordinate of the last and first space point.
# NOTE: Rather, subtract the minima and maxima of the corrdinates of the re-meshed mesh
for i in range(len(times)):
length[i] = np.max(geometry[i])-np.min(geometry[i][0])
fig1, ax1 = plt.subplots(1, 1, figsize=(8, 8))
ax1.set_xlabel('$t$')
ax1.set_ylabel('$L(t)$')
ax1.set_xlim(np.min(times), np.max(times))
ax1.set_ylim(np.min(length), np.max(length)+0.01)
ax1.tick_params(axis='both', which='both', labelsize=10)
ax1.grid(True)
ax1.plot(times, length, marker='o', ms=3, linestyle='None')
ax1.grid(True)
plt.savefig("%s.png" %params['timestamp'])
np.save('%s_length.npy' %params['timestamp'], length)
......@@ -202,18 +237,16 @@ def visualize(params, DIR=''):
ax2.set_xlabel('$t$')
ax2.set_ylabel('$L(t)$')
ax2.tick_params(axis='both', which='both', labelsize=10)
ax2.grid(True)
ax2.loglog(times, length, marker='o', ms=3, linestyle='None')
ax2.grid(True)
plt.savefig("%s_loglog.png" %params['timestamp'])
fig3, ax3 = plt.subplots(1, 1, figsize=(8, 8))
ax3.set_xlabel('$t$')
ax3.set_ylabel('$L(t)$')
ax3.tick_params(axis='both', which='both', labelsize=10)
ax3.grid(True)
ax3.semilogy(times, length, marker='o', ms=3, linestyle='None')
ax3.grid(True)
plt.savefig("%s_semilog.png" %params['timestamp'])
# dldt vs t
......@@ -224,34 +257,30 @@ def visualize(params, DIR=''):
ax11.set_xlim(np.min(times), np.max(times))
ax11.set_ylim(np.min(dldt), np.max(dldt)+0.01)
ax11.tick_params(axis='both', which='both', labelsize=10)
ax11.grid(True)
ax11.plot(times, dldt, marker='o', ms=3, linestyle='None')
ax11.grid(True)
plt.savefig("%s_dldt.png" %params['timestamp'])
fig21, ax21 = plt.subplots(1, 1, figsize=(8, 8))
ax21.set_xlabel('$t$')
ax21.set_ylabel('$dL/dt$')
ax21.tick_params(axis='both', which='both', labelsize=10)
ax21.grid(True)
ax21.loglog(times, dldt, marker='o', ms=3, linestyle='None')
ax21.grid(True)
plt.savefig("%s_loglog_dldt.png" %params['timestamp'])
fig31, ax31 = plt.subplots(1, 1, figsize=(8, 8))
ax31.set_xlabel('$t$')
ax31.set_ylabel('$dL/dt$')
ax31.tick_params(axis='both', which='both', labelsize=10)
ax31.grid(True)
ax31.semilogy(times, dldt, marker='o', ms=3, linestyle='None')
ax31.grid(True)
plt.savefig("%s_semilog_dldt.png" %params['timestamp'])
if __name__ == '__main__':
import argparse, json
parser = argparse.ArgumentParser()
parser.add_argument('-j','--jsonfile', help='parameters file', required=True)
parser.add_argument('--onscreen', action=argparse.BooleanOptionalAction)
args = parser.parse_args()
with open(args.jsonfile) as jsonFile:
......
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