u, v, rho, 1D and 2D, initial condition dependence for 2D problem

u_v_rho/1D/compare_growth_rates.py

+import numpy as np
+import matplotlib.pyplot as plt
+
+l1 = np.load()
+l10 =
+l100 = 
+
+times = np.arange(0, 500.1, 0.1)
+
+plt.plot(times, l1, lw =3, label='$l_0=1$')
+plt.plot(times, l10, lw =3, label='$l_0=10$')
+plt.plot(times, l100, lw =3, label='$l_0=100$')
+
+plt.show()
\ No newline at end of file

u_v_rho/1D/kymographs.py

+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.colors import SymLogNorm
+
+u1 = np.load('060625-155017_u.npy')
+geometry1 = np.load('060625-155017_geometry.npy')
+times1 = np.load('060625-155017_times.npy')
+v1 = np.load('060625-155017_v.npy')
+rho1 = np.load('060625-155017_rho.npy')
+print(times1)
+
+u2 = np.load('060625-162043_u.npy')
+geometry2 = np.load('060625-162043_geometry.npy')
+times2 = np.load('060625-162043_times.npy')
+v2 = np.load('060625-162043_v.npy')
+rho2 = np.load('060625-162043_rho.npy')
+
+estress1 = np.load('060625-155017_elastic_stress.npy')
+vstress1 = np.load('060625-155017_viscous_stress.npy')
+astress1 = np.load('060625-155017_active_stress.npy')
+
+estress2 = np.load('060625-162043_elastic_stress.npy')
+vstress2 = np.load('060625-162043_viscous_stress.npy')
+astress2 = np.load('060625-162043_active_stress.npy')
+min_u1 = min(min(sublist) for sublist in u1)
+max_u1 = max(max(sublist) for sublist in u1)
+min_u2 = min(min(sublist) for sublist in u2)
+max_u2 = max(max(sublist) for sublist in u2)
+min_v1 = min(min(sublist) for sublist in v1)
+max_v1 = max(max(sublist) for sublist in v1)
+min_v2 = min(min(sublist) for sublist in v2)
+max_v2 = max(max(sublist) for sublist in v2)
+min_rho1 = min(min(sublist) for sublist in rho1)
+max_rho1 = max(max(sublist) for sublist in rho1)
+min_rho2 = min(min(sublist) for sublist in rho2)
+max_rho2 = max(max(sublist) for sublist in rho2)
+min_u = min(min_u1, min_u2)
+max_u = max(max_u1, max_u2)
+min_v = min(min_v1, min_v2)
+max_v = max(max_v1, max_v2)
+min_rho = min(min_rho1, min_rho2)
+max_rho = max(max_rho1, max_rho2)
+
+min_estress1 = min(min(sublist) for sublist in estress1)
+max_estress1 = max(max(sublist) for sublist in estress1)
+min_vstress1 = min(min(sublist) for sublist in vstress1)
+max_vstress1 = max(max(sublist) for sublist in vstress1)
+min_astress1 = min(min(sublist) for sublist in astress1)
+max_astress1 = max(max(sublist) for sublist in astress1)
+min_estress2 = min(min(sublist) for sublist in estress2)
+max_estress2 = max(max(sublist) for sublist in estress2)
+min_vstress2 = min(min(sublist) for sublist in vstress2)
+max_vstress2 = max(max(sublist) for sublist in vstress2)
+min_astress2 = min(min(sublist) for sublist in astress2)
+max_astress2 = max(max(sublist) for sublist in astress2)
+min_stress = min(min_estress1, min_estress2, min_vstress1, min_vstress2, min_astress1, min_astress2)
+max_stress = max(max_estress1, max_estress2, max_vstress1, max_vstress2, max_astress1, max_astress2)
+
+fig1, ax1 = plt.subplots(1, 2, figsize=(8, 4))
+norm = SymLogNorm(linthresh=1, vmin=min_u, vmax=max_u, base=10)
+
+im0 = ax1[0].imshow(u1, aspect='auto', extent=[min(geometry1), max(geometry1), 
+                                                min(times1), max(times1)], 
+                                                origin='lower', cmap='PiYG')
+im1 = ax1[1].imshow(u2, aspect='auto', extent=[min(geometry2), max(geometry2),
+                                                min(times2), max(times2)], 
+                                                origin='lower', cmap='PiYG')
+cbar = fig1.colorbar(im0, ax=ax1, orientation='vertical', fraction=0.02, pad=0.04)
+cbar.set_label(r'$u$', labelpad=10)
+
+plt.savefig('u.png', dpi=1000, bbox_inches='tight')
+
+fig2, ax2 = plt.subplots(1, 2, figsize=(8, 4))
+norm = SymLogNorm(linthresh=1, vmin=min_v, vmax=max_v, base=10)
+im0 = ax2[0].imshow(v1, aspect='auto', extent=[min(geometry1), max(geometry1), 
+                                                min(times1), max(times1)], 
+                                                origin='lower', cmap='seismic')
+im1 = ax2[1].imshow(v2, aspect='auto', extent=[min(geometry2), max(geometry2),
+                                                min(times2), max(times2)], 
+                                                origin='lower', cmap='seismic')
+cbar = fig2.colorbar(im0, ax=ax2, orientation='vertical', fraction=0.02, pad=0.04)
+cbar.set_label(r'$v$', labelpad=10)
+plt.savefig('v.png', dpi=1000, bbox_inches='tight')
+
+fig3, ax3 = plt.subplots(1, 2, figsize=(8, 4))
+norm = SymLogNorm(linthresh=1, vmin=min_rho, vmax=max_rho, base=10)
+im0 = ax3[0].imshow(rho1, aspect='auto', extent=[min(geometry1), max(geometry1), 
+                                                min(times1), max(times1)], 
+                                                origin='lower', cmap='Oranges')
+
+im1 = ax3[1].imshow(rho2, aspect='auto', extent=[min(geometry2), max(geometry2),
+                                                min(times2), max(times2)], 
+                                                origin='lower', cmap='Oranges')
+
+cbar = fig3.colorbar(im0, ax=ax3, orientation='vertical', fraction=0.02, pad=0.04)
+cbar.set_label(r'$\rho$', labelpad=10)
+plt.savefig('rho.png', dpi=1000, bbox_inches='tight')
+
+fig4, ax4 = plt.subplots(2, 3, figsize=(16, 6))
+im0 = ax4[0, 0].imshow(estress1, aspect='auto', extent=[min(geometry1), max(geometry1),
+                                                min(times1), max(times1)], 
+                                                origin='lower', vmin=min_stress, vmax=max_stress, cmap='coolwarm_r')
+im1 = ax4[0, 1].imshow(vstress1, aspect='auto', extent=[min(geometry1), max(geometry1),
+                                                min(times1), max(times1)], 
+                                                origin='lower', vmin=min_stress, vmax=max_stress, cmap='coolwarm_r')
+im2 = ax4[0, 2].imshow(astress1, aspect='auto', extent=[min(geometry1), max(geometry1),
+                                                min(times1), max(times1)], 
+                                                origin='lower', vmin=min_stress, vmax=max_stress, cmap='coolwarm_r')
+im3 = ax4[1, 0].imshow(estress2, aspect='auto', extent=[min(geometry2), max(geometry2),
+                                                min(times2), max(times2)], 
+                                                origin='lower', vmin=min_stress, vmax=max_stress, cmap='coolwarm_r')
+im4 = ax4[1, 1].imshow(vstress2, aspect='auto', extent=[min(geometry2), max(geometry2),
+                                                min(times2), max(times2)], 
+                                                origin='lower', vmin=min_stress, vmax=max_stress, cmap='coolwarm_r')
+im5 = ax4[1, 2].imshow(astress2, aspect='auto', extent=[min(geometry2), max(geometry2),
+                                                min(times2), max(times2)], 
+                                                origin='lower', vmin=min_stress, vmax=max_stress, cmap='coolwarm_r')
+cbar = fig4.colorbar(im0, ax=ax4, orientation='vertical', fraction=0.02, pad=0.04)
+cbar.set_label(r'$\sigma_{elastic}$', labelpad=10)
+plt.savefig('stress.png', dpi=1000, bbox_inches='tight')
+plt.show()
+
+

u_v_rho/1D/parameters.json

+{
+    "active_death": true,
+    "active_stress_setpoint": 4,
+    "average_rho": 1.0,
+    "diffusion_rho": 0.1,
+    "elasticity": 0.5,
+    "friction": 1.0,
+    "lamda": 1.0,
+    "maxtime": 50.0,
+    "noise_level": 0.0,
+    "resolution": 128,
+    "saturation_rho": 1.0,
+    "savetime": 0.1,
+    "system_size": 1.0,
+    "timestep": 0.01,
+    "turnover_rho": 1,
+    "viscosity": 1,
+    "zero_rho_boundary": false
+}
\ No newline at end of file

u_v_rho/1D/run_tissue.py

+import json, datetime
+import os
+from tissue import OneDimTissue
+from viz_tissue 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)])[:, 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 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 = OneDimTissue(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='')
\ No newline at end of file

u_v_rho/1D/scan_tissue.py

+from __future__ import print_function
+import json
+import itertools
+import os
+import datetime
+import time
+import numpy as np
+
+from tissue import OneDimTissue
+
+assert os.path.isfile("parameters.json"), 'parameters.json file not found'
+with open("parameters.json") as jsonFile:
+    parameters = json.load(jsonFile)
+
+#=========================================
+# parameters to scan
+
+rho_hat = np.arange(0.5, 5.5, 0.5)
+k = np.arange(0.5, 5.5, 0.5 )
+
+plist = list(itertools.product(rho_hat, k))
+
+#=========================================
+
+# scan
+for p in plist:
+    rho_hat, k = p
+    print('\n------------------\n')
+    print('active_stress_setpoint = ',rho_hat, ', elasticity = ', k)
+    
+    timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
+    params = parameters.copy()
+    params['timestamp'] = timestamp
+    params['active_stress_setpoint'] = float(rho_hat)
+    params['elasticity'] = float(k)
+
+    t = OneDimTissue(params)
+    t.solve()
+
+    parametersFile = timestamp + '__parameters.json'
+
+    with open(parametersFile, "w") as fp:
+        json.dump(params, fp, indent=4, sort_keys=True, default=str)
+
+    from viz_tissue import visualize
+    with open(parametersFile) as jFile:
+        parameters = json.load(jFile)
+    visualize(parameters)
+    time.sleep(2)
\ No newline at end of file

u_v_rho/1D/tissue.py

+import numpy as np
+import dolfin as df
+import progressbar
+
+df.set_log_level(df.LogLevel.ERROR)
+df.parameters['form_compiler']['optimize'] = True
+
+class OneDimTissue(object):
+    def __init__(self, parameters, mesh=None):
+        # read in parameters
+        for key in parameters:
+            setattr(self, key, parameters[key])
+
+        # 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)
+        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
+        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, self.mixed_element)
+        
+        # boundary condition on rho
+        if self.zero_rho_boundary:
+            self.rho_boundary = 0.0
+        else:
+            self.rho_boundary = self.average_rho
+
+        self.bc = df.DirichletBC(self.function_space.sub(2), 
+                                 df.Constant(self.rho_boundary), 
+                                 'on_boundary')
+        
+        # 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)
+
+    def active_stress(self, rho):
+        if self.active_death:
+            return (self.lamda * rho * ((rho - self.active_stress_setpoint)/(rho**2 + self.saturation_rho**2)) * df.Identity(1))
+        else:
+            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.05:  # if the width of the cell > threshold, 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 
+        viscous_stress = self.viscosity * eps_v
+        return (elastic_stress + viscous_stress)
+    
+    def stress(self, u, v, rho):
+        return (self.passive_stress(u, v) + self.active_stress(rho))
+
+    def diffusion_reaction_rho(self, rho, trho):
+        return (self.diffusion_rho * df.inner(df.nabla_grad(rho), 
+                                              df.nabla_grad(trho)) 
+                - self.turnover_rho * df.inner(rho*(1 - rho/self.average_rho), 
+                                               trho)
+                )
+
+    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:
+            base_rho = self.average_rho
+        rho0 = df.interpolate(df.Expression(
+                            'base_rho + 0.1 * cos(PI*x[0]/L)*cos(PI*x[0]/L)', 
+                            L=self.system_size,
+                            base_rho=base_rho,
+                            degree=1, PI=np.pi), 
+                            self.function_space.sub(2).collapse())
+
+        
+        # velocity from u0 and rho0
+        VFS = self.function_space.sub(1).collapse()
+        v0  = df.Function(VFS)
+        tv0 = df.TestFunction(VFS)
+        v0form = (self.friction * df.inner(v0, tv0) 
+                    + df.inner(self.stress(u0, v0, rho0), 
+                               self.epsilon(tv0))
+                ) * df.dx
+        
+        df.solve(v0form == 0, v0)
+        df.assign(self.function0, [u0, v0, rho0])
+
+    def setup_weak_forms(self):
+        u0, v0, rho0  = df.split(self.function0)
+        u,  v, rho    = df.split(self.function)
+        tu, tv, trho  = df.TestFunctions(self.function_space)
+                  
+        uform = (df.inner((u - u0)/self.timestep, tu)
+                - df.inner(v0, tu))
+        
+        vform = (self.friction * df.inner(v, tv) 
+                + df.inner(self.stress(u, v, rho), 
+                           self.epsilon(tv)))
+
+        rhoform = (df.inner((rho - rho0)/self.timestep, trho)
+                + self.advection(rho0, v0, trho)
+                + 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)
+        
+        # time-variables
+        self.time = initTime
+        savesteps = round(self.savetime/self.timestep)
+        maxsteps = round(self.maxtime/self.timestep)
+
+        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)
+
+        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)
+
+        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)
+
+            # 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)
+
+            # 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 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()
+        
\ No newline at end of file

u_v_rho/1D/viz_tissue.py

+import dolfin as df
+import numpy as np
+import os
+import h5py
+from matplotlib.widgets import Slider
+import progressbar
+import matplotlib.pyplot as plt
+from tempfile import TemporaryDirectory
+import scipy
+
+plt.rcParams['font.size'] = 15
+
+def visualize(params, DIR=''):
+    savesteps = round(params['maxtime'] / params['savetime'])
+    times = np.arange(savesteps + 1) * params['savetime']
+    
+    # Read mesh geometry and topology from h5 file
+    var = 'density'
+    h5 = h5py.File(os.path.join(DIR, '%s_%s.h5' % (params['timestamp'], var)), "r")
+    
+    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)])[:, 0])
+        topology.append(np.array(h5['%s/%s_%d/mesh/topology' % (var, var, i)]))
+
+    h5.close()
+
+    # Initialize field storage
+    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']))
+
+    for k in progressbar.ProgressBar()(range(len(times))):
+        # Create and process the mesh
+        mesh = df.Mesh()
+        editor = df.MeshEditor()
+        editor.open(mesh, 'interval', 1, 1)
+
+        # Initialize 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()
+
+        # Read field values for the current mesh
+        SFS = df.FunctionSpace(mesh, 'P', 1)
+        VFS = df.VectorFunctionSpace(mesh, 'P', 1)
+        ui, vi, rhoi = df.Function(VFS), df.Function(VFS), df.Function(SFS)
+
+        uFile.read_checkpoint(ui, 'displacement', k)
+        vFile.read_checkpoint(vi, 'velocity', k)
+        rhoFile.read_checkpoint(rhoi, 'density', k)
+
+        u.append(ui.compute_vertex_values(mesh))
+        v.append(vi.compute_vertex_values(mesh))
+        rho.append(rhoi.compute_vertex_values(mesh))
+
+        e_stress = params['elasticity'] * ui.dx(0)
+        e_stress = df.project(e_stress, VFS)
+        elastic_stress.append(e_stress.compute_vertex_values(mesh))
+
+        v_stress = params['viscosity'] * vi.dx(0)
+        v_stress = df.project(v_stress, VFS)
+        viscous_stress.append(v_stress.compute_vertex_values(mesh))
+
+        if params['active_death'] == False:
+            a_stress = params['lamda'] * rhoi / (rhoi + params['saturation_rho'])
+        else:
+            a_stress = (params['lamda'] * rhoi *(rhoi - params['active_stress_setpoint']) 
+                                   / (rhoi**2 + params['saturation_rho']**2))
+
+        active_stress.append(df.project(a_stress, SFS).compute_vertex_values(mesh))
+
+        del mesh
+
+    uFile.close()
+    vFile.close()
+    rhoFile.close()
+
+    # Sort the geometry and fields
+    for k in range(len(times)):
+        sorted_indices = np.argsort(geometry[k])
+        geometry[k] = geometry[k][sorted_indices]
+        u[k] = u[k][sorted_indices]
+        v[k] = v[k][sorted_indices]
+        rho[k] = rho[k][sorted_indices]
+        elastic_stress[k] = elastic_stress[k][sorted_indices]
+        viscous_stress[k] = viscous_stress[k][sorted_indices]
+        active_stress[k] = active_stress[k][sorted_indices]
+
+
+    # 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, min_geometry = (map(lambda x: min(min(sublist) for sublist in x), 
+                                                [elastic_stress, viscous_stress, active_stress, geometry])
+                                            )
+    max_estress, max_vstress, max_astress, max_geometry = (map(lambda x: max(max(sublist) for sublist in x), 
+                                                [elastic_stress, viscous_stress, active_stress, geometry])
+                                            )    
+    # interactive plot
+    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/\ell \kappa$')
+    axes[2].set_ylabel(r'$u/\ell$')
+    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], ms=3, color='g', lw = 2)
+    velplot, = axes[1].plot(geometry[0], v[0], ms=3, color='r', lw = 2)
+    uplot, = axes[2].plot(geometry[0], u[0], ms=3, color='b', lw = 2)
+    
+    def update(value):
+        ti = np.abs(times-value).argmin()
+        rhoplot.set_ydata(rho[ti])
+        rhoplot.set_xdata(geometry[ti])
+        velplot.set_ydata(v[ti])
+        velplot.set_xdata(geometry[ti])
+        uplot.set_ydata(u[ti])
+        uplot.set_xdata(geometry[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()
+
+    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()
+
+    figss, axss = plt.subplots(3, 1, sharex=True, figsize=(8,8))
+    axss[0].set_xlabel(r'$x$')
+    axss[0].set_ylabel(r'$u$')
+    axss[1].set_ylabel(r'$v$')
+    axss[2].set_ylabel(r'$\rho$')
+    axss[0].set_xlim(min_geometry-5, max_geometry+5)
+    axss[0].set_ylim(min_u-0.1, max_u+0.1)
+    axss[1].set_ylim(min_v-0.1, max_v+0.1)
+    axss[2].set_ylim(min_rho-0.1, max_rho+0.1)
+    axss[0].plot(geometry[-1], u[-1], ms=1, color='b', lw = 2)
+    axss[1].plot(geometry[-1], v[-1], ms=1, color='r', lw = 2)
+    axss[2].plot(geometry[-1], rho[-1], ms=1, color='g', lw = 2)
+    plt.savefig("%sss_fields.png" %params['timestamp'])
+
+    slope = 1
+    figss_u, axss_u = plt.subplots(1, 1, sharex=True, figsize=(8,8))
+    axss_u.set_xlabel(r'$x$')
+    axss_u.set_ylabel(r'$u/\ell$')
+    axss_u.set_xlim(min_geometry-5, max_geometry+5)
+    axss_u.set_ylim(min_u-0.1, max_u+0.1)
+    axss_u.plot(geometry[-1], u[-1], ms=1, color='b', lw = 1, label=r'$u_{numerical}$')
+    axss_u.plot(geometry[-1], (slope) * geometry[-1], color='r', linestyle=':', label=r'$u_{guess}$')
+    # slope 1 for UG, slope for CG in terms of parameters
+    plt.grid(True)
+    plt.legend()
+    plt.savefig("%sss_u_with_u_guess.png" %params['timestamp'])
+
+    fig_u_near_bdry, ax_u_near_bdry = plt.subplots(1, 1, sharex=True, figsize=(8,8))
+    ax_u_near_bdry.set_xlabel(r'$x$')
+    ax_u_near_bdry.set_ylabel(r'$u$')
+    # axss[2].set_ylim(min_u-0.1, max_u+0.1)
+    ax_u_near_bdry.plot(geometry[-1][-10:], u[-1][-10:], ms=3, color='b', lw=2)    
+    plt.savefig("%sss_u_near_bdry.png" %params['timestamp'])
+
+    fig_rho_near_bdry, ax_rho_near_bdry = plt.subplots(1, 1, sharex=True, figsize=(8,8))
+    ax_rho_near_bdry.set_xlabel(r'$x$')
+    ax_rho_near_bdry.set_ylabel(r'$\rho$')
+    # axss[2].set_ylim(min_u-0.1, max_u+0.1)
+    ax_rho_near_bdry.plot(geometry[-1][-20:], rho[-1][-20:], ms=3, color='b', lw=2)    
+    plt.savefig("%sss_rho_near_bdry_2.png" %params['timestamp'])
+
+    fig_rho_near_bdry2, ax_rho_near_bdry2 = plt.subplots(1, 1, sharex=True, figsize=(8,8))
+    ax_rho_near_bdry2.set_xlabel(r'$x$')
+    ax_rho_near_bdry2.set_ylabel(r'$\rho$')
+    # axss[2].set_ylim(min_u-0.1, max_u+0.1)
+    ax_rho_near_bdry2.plot(geometry[-1][-5:], rho[-1][-5:], ms=3, color='b', lw=2)    
+    plt.savefig("%sss_rho_near_bdry4.png" %params['timestamp'])
+
+    fig_v_near_bdry, ax_v_near_bdry = plt.subplots(1, 1, sharex=True, figsize=(8,8))
+    ax_v_near_bdry.set_xlabel(r'$x$')
+    ax_v_near_bdry.set_ylabel(r'$v$')
+    # axss[2].set_ylim(min_u-0.1, max_u+0.1)
+    ax_v_near_bdry.plot(geometry[-1][-10:], v[-1][-10:], ms=3, color='b', lw=2)    
+    plt.savefig("%sss_v_near_bdry.png" %params['timestamp'])
+
+    # plotting the stresses
+    # figure, axes1 = plt.subplots(3, 1, sharex=True, figsize=(8,8))
+    # axes1[-1].set_xlabel(r'$x$')
+    # axes1[2].set_ylabel('Elastic stress')
+    # axes1[1].set_ylabel('Viscous stress')
+    # axes1[0].set_ylabel('Active stress')
+    # 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], lw = 2, color='b')
+    # v_stress_plot, = axes1[1].plot(geometry[0], viscous_stress[0], lw = 2, color='r')
+    # a_stress_plot, = axes1[0].plot(geometry[0], active_stress[0], lw = 2, color='g')
+    
+    # def update2(value):
+    #     ti = np.abs(times-value).argmin()
+
+    #     e_stress_plot.set_ydata(elastic_stress[ti])
+    #     e_stress_plot.set_xdata(geometry[ti])
+
+    #     v_stress_plot.set_ydata(viscous_stress[ti])
+    #     v_stress_plot.set_xdata(geometry[ti])
+
+    #     a_stress_plot.set_ydata(active_stress[ti])
+    #     a_stress_plot.set_xdata(geometry[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(update2)
+
+    # figss_stress, axss_stress = plt.subplots(3, 1, sharex=True, figsize=(8,8))
+    # axss_stress[-1].set_xlabel(r'$x$')
+    # axss_stress[0].set_ylabel('Active stress')
+    # axss_stress[1].set_ylabel('Viscous stress')
+    # axss_stress[2].set_ylabel('Elastic stress')
+    # axss_stress[0].set_xlim(min_geometry-0.5, max_geometry+0.5)
+    # axss_stress[0].set_ylim(min_astress-0.1, max_astress+0.1)
+    # axss_stress[1].set_ylim(min_vstress-0.1, max_vstress+0.1)
+    # axss_stress[2].set_ylim(min_estress-0.1, max_estress+0.1)
+    # axss_stress[0].plot(geometry[-1], active_stress[-1], ms=3, color='g', lw = 2)
+    # axss_stress[1].plot(geometry[-1], viscous_stress[-1], ms=3, color='r', lw = 2)
+    # axss_stress[2].plot(geometry[-1], elastic_stress[-1], ms=3, color='b', lw = 2)
+    # plt.savefig("%sss_stress.png" %params['timestamp'])
+
+    # # make movie
+    # 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 = 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])
+    print('length_numerics=', length[-1])
+
+
+    # #NOTE: Alternative length calculation
+    integrated_strain = []
+
+    for k in range(len(times)):
+        strain = elastic_stress[k] / params['elasticity']
+
+    # Integrate strain over the geometry using the trapezoidal rule
+    integration = np.trapezoid(strain, geometry[k])
+    integrated_strain.append(integration)
+
+    print('length_calc=', integrated_strain[-1] + params['system_size'])
+
+
+    fig1, ax1 = plt.subplots(1, 1, figsize=(8, 8))
+    ax1.set_xlabel('$t$')
+    ax1.set_ylabel('$L(t)$')
+    ax1.set_xlim(np.min(times)-1, np.max(times))
+    ax1.set_ylim(np.min(length), np.max(length)+0.1)
+    # ax1.tick_params(axis='both', which='both')
+    ax1.plot(times, length, lw = 3, color='darkblue')
+    # ax1.grid(True)
+    plt.savefig("%s_length.png" %params['timestamp'])
+    np.save('%s_length.npy' %params['timestamp'], length)
+   
+    # # dldt vs t
+    dldt = np.gradient(length, times)
+    print('dldt[-1]:', dldt[-1])
+    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')
+    ax31.semilogy(times, np.abs(dldt), lw = 3, color='darkblue')
+    # 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)
+    args = parser.parse_args()
+
+    with open(args.jsonfile) as jsonFile:
+        params = json.load(jsonFile)
+
+    visualize(params, DIR=os.path.dirname(args.jsonfile))
\ No newline at end of file

u_v_rho/2D/circle_to_ellipse.py

+import numpy as np
+import dolfin as df
+import matplotlib.pyplot as plt
+from matplotlib.animation import FuncAnimation
+from meshpy.triangle import MeshInfo, build
+
+def triangulate_with_meshpy(boundary_points, max_volume=0.005):
+    mesh_info = MeshInfo()
+    mesh_info.set_points(boundary_points)
+    mesh_info.set_facets([(i, (i+1) % len(boundary_points)) for i in range(len(boundary_points))])
+    mesh = build(mesh_info, max_volume=max_volume)
+    points = np.array(mesh.points)
+    cells = np.array(mesh.elements)
+    return points, cells
+
+def generateEllipseMesh(r0=1.0, scale=0.1, nsegments=36, max_volume=0.005, time=0.0):
+    # Generate the vertices of the ellipse
+    theta = np.linspace(0, 2 * np.pi, nsegments, endpoint=False)
+    vertices = []
+    for th in theta:
+        # Deform the circle into an ellipse over time
+        x = r0 * (1 + scale * np.cos(2 * np.pi * time)) * np.cos(th)
+        y = r0 * (1 - scale * np.sin(2 * np.pi * time)) * np.sin(th)
+        vertices.append((x, y))
+    vertices = np.array(vertices)
+
+    points, cells = triangulate_with_meshpy(vertices, max_volume=max_volume)
+    # Convert the mesh to a FEniCS mesh
+    mesh = df.Mesh()
+    editor = df.MeshEditor()
+    editor.open(mesh, "triangle", 2, 2)
+    editor.init_vertices(len(points))
+    editor.init_cells(len(cells))
+    for i, point in enumerate(points):
+        editor.add_vertex(i, point)
+    for i, cell in enumerate(cells):
+        editor.add_cell(i, cell)
+    editor.close()
+
+    return mesh
+
+# Parameters
+radius = 1.0
+scale = 0.3
+nsegments = 36
+max_volume = 0.005
+frames = 100  # Number of frames in the animation
+
+# Create a figure for the animation
+fig, ax = plt.subplots()
+
+# Initialize the plot
+mesh = generateEllipseMesh(r0=radius, scale=scale, nsegments=nsegments, max_volume=max_volume, time=0.0)
+plot = df.plot(mesh)  # Plot directly without passing ax
+ax.set_title("Deforming Circle to Ellipse")
+ax.set_xlim(-1.5, 1.5)
+ax.set_ylim(-1.5, 1.5)
+
+# Update function for the animation
+def update(frame):
+    ax.clear()
+    ax.axis('off')
+    mesh = generateEllipseMesh(r0=radius, scale=scale, nsegments=nsegments, max_volume=max_volume, time=frame / frames)
+    df.plot(mesh)  # Plot directly without passing ax
+    # ax.set_title(f"Time: {frame / frames:.2f}")
+    ax.set_xlim(-1.5, 1.5)
+    ax.set_ylim(-1.5, 1.5)
+
+# Create the animation
+ani = FuncAnimation(fig, update, frames=frames, interval=50)
+
+# Save the animation as a movie
+ani.save("circle_to_ellipse.mp4", writer="ffmpeg")
+
+plt.show()
\ No newline at end of file

u_v_rho/2D/generate_ellipse_mesh.py

+# import dolfin as df
+# import mshr as ms
+# from matplotlib import pyplot as plt
+
+# center = df.Point(0.0, 0.0)
+# semi_major_axis = 1.25
+# semi_minor_axis = 1.0
+# nsegments = 60
+
+# domain = ms.Ellipse(center, semi_major_axis, semi_minor_axis, nsegments)
+
+# mresolution = nsegments/3
+# mesh = ms.generate_mesh(domain, mresolution)
+# df.plot(mesh)
+# plt.show()
+# mname = 'ellipse_majoraxis_%3.2f_minoraxis_%3.2f' % (semi_major_axis, semi_minor_axis)
+
+# mFile = df.File(mname + '.xml.gz')
+# mFile << mesh
+
+import dolfin as df
+import mshr as ms
+from matplotlib import pyplot as plt
+
+# # Define the corner points of the square
+# p1 = df.Point(0.0, 0.0)
+# p2 = df.Point(1, 1.0)  # Adjust these points to change the size of the square
+
+# # Create a rectangular domain
+# domain = ms.Rectangle(p1, p2)
+# Define the center and radii of the ellipse
+center = df.Point(0.0, 0.0)  # Center of the ellipse
+radius_x = 2  # Radius along the x-axis
+radius_y = 1.0  # Radius along the y-axis
+
+# Create an elliptical domain
+domain = ms.Ellipse(center, radius_x, radius_y)
+
+# Define the mesh resolution
+mresolution = 15  # Adjust this value to change the mesh resolution
+
+# Generate the mesh
+mesh = ms.generate_mesh(domain, mresolution)
+
+# Plot the mesh
+df.plot(mesh)
+plt.show()
+
+# Save the mesh to a file
+mname = 'ellipse_a_%s_b_' \
+'%s_resolution_15' %(radius_x, radius_y)
+mFile = df.File(mname + '.xml.gz')
+mFile << mesh
\ No newline at end of file

u_v_rho/2D/generate_meshes.py

+import numpy as np
+import dolfin as df
+import argparse
+from meshpy.triangle import MeshInfo, build
+import matplotlib.pyplot as plt
+def triangulate_with_meshpy(boundary_points, max_volume=0.005):
+    mesh_info = MeshInfo()
+    mesh_info.set_points(boundary_points)
+    mesh_info.set_facets([(i, (i+1) % len(boundary_points)) for i in range(len(boundary_points))])
+    mesh = build(mesh_info, max_volume=max_volume)
+    points = np.array(mesh.points)
+    cells = np.array(mesh.elements)
+    return points, cells
+
+def generateMesh(r0=1.0, scale=0.1, nmode=3, nsegments=36, max_volume=0.005):
+    # Generate the vertices of the polygon
+    theta = np.linspace(0, 2 * np.pi, nsegments, endpoint=False)
+    vertices = []
+    for th in theta:
+        r = r0 * (1 + scale * np.cos(nmode * th))
+        x, y = r * np.cos(th), r * np.sin(th)
+        vertices.append((x, y))
+    vertices = np.array(vertices)
+
+    points, cells = triangulate_with_meshpy(vertices, max_volume=max_volume)
+    # Convert the mesh to a FEniCS mesh
+    mesh = df.Mesh()
+    editor = df.MeshEditor()
+    editor.open(mesh, "triangle", 2, 2)
+    editor.init_vertices(len(points))
+    editor.init_cells(len(cells))
+    for i, point in enumerate(points):
+        editor.add_vertex(i, point)
+    for i, cell in enumerate(cells):
+        editor.add_cell(i, cell)
+    editor.close()
+
+    return mesh
+
+
+parser = argparse.ArgumentParser()
+parser.add_argument('-radius', type=float, help='radius of circle')
+parser.add_argument('-mode', type=int, help='frequency of mode')
+parser.add_argument('-scale', type=float, 
+        help='scale for spharm deformation [0,1)')
+args = parser.parse_args()
+assert(args.radius is not None), 'Need radius'
+assert(args.mode is not None), 'Need mode'
+assert(args.scale is not None), 'Need scale'
+assert(args.radius>=0.0), 'Radius should greater 1'
+assert(0.0<=args.scale<1.0), 'Scale between 0 and 1'
+
+scale = args.scale
+mode = args.mode
+radius = args.radius
+
+mesh2 = generateMesh(r0=radius, scale=scale, nmode=mode, nsegments=20, max_volume=0.01)
+mesh3 = generateMesh(r0=radius, scale=scale, nmode=mode, nsegments=3*30, max_volume=0.005)
+mesh4 = generateMesh(r0=radius, scale=scale, nmode=mode, nsegments=4*30, max_volume=0.002)
+
+mname = 'mesh_radius%3.2f_mode%d_scale%3.2f__' % (radius, mode, scale)
+
+mFile = df.File(mname+'2.xml.gz')
+mFile<<mesh2
+# mFile = df.File(mname+'3.xml.gz')
+# mFile<<mesh3
+# mFile = df.File(mname+'4.xml.gz')
+# mFile<<mesh4
+# mFile = df.File(mname+'5.xml.gz')
+# mFile<<mesh5
+# mFile = df.File(mname+'6.xml.gz')
+# mFile<<mesh6
+# mFile = df.File(mname+'7.xml.gz')
+# mFile<<mesh7
+# mFile = df.File(mname+'8.xml.gz')
+# mFile<<mesh8
+# mFile = df.File(mname+'9.xml.gz')
+# mFile<<mesh9
+# mFile = df.File(mname+'10.xml.gz')
+# mFile<<mesh10
+
+plt.figure()
+df.plot(mesh3)
+plt.title('Mesh 2')
+plt.savefig(mname + '2_plot.png')
+
+# Plot and save mesh3
+# plt.figure()
+# df.plot(mesh3)
+# plt.title('Mesh 3')
+# plt.savefig(mname + '3_plot.png')
+
+# # Plot and save mesh4
+# plt.figure()
+# df.plot(mesh4)
+# plt.title('Mesh 4')
+# plt.savefig(mname + '4_plot.png')

u_v_rho/2D/growing_domain.py

+import dolfin as df
+import numpy as np
+import progressbar
+import meshzoo
+
+df.set_log_level(df.LogLevel.ERROR)
+df.parameters['form_compiler']['optimize'] = True
+
+class Growing_Domain(object):
+    def __init__(self, parameters, mesh = None):
+        # read in parameters
+        for key in parameters:
+            setattr(self, key, parameters[key])
+
+
+        # if no mesh is given as initial condition, create one
+        if mesh == 'None':
+            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()
+
+        else:
+            self.mesh = df.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])
+
+        # define function space with this mixed element
+        self.function_space = df.FunctionSpace(self.mesh, mixed_element)
+
+        # dirichlet boundaries for rho
+        if self.zero_rho_boundary:
+            bc_rho = df.Constant(0.0)
+        else:
+            bc_rho = df.Constant(self.average_rho)
+
+        self.bc = df.DirichletBC(self.function_space.sub(2), bc_rho, 'on_boundary')
+
+        # define the reqd functions on this function space
+        self.f  = df.Function(self.function_space)                # f at current time
+        self.f0 = df.Function(self.function_space)                # f at t =0
+
+    def advection(self, conc, vel, tconc):
+        return (-df.inner(vel * conc, df.nabla_grad(tconc)))
+    
+    def active_stress(self, rho):               
+        return (self.lamda * rho * (rho - self.stress_setpoint_rho)/(rho**2 + self.saturation_rho**2) 
+                    * df.Identity(self.dimension)
+                    )
+        
+    def epsilon(self, v):
+        return df.sym(df.nabla_grad(v))
+    
+    def passive_stress(self, u, v):
+        eps_u, eps_v = self.epsilon(u), self.epsilon(v)
+
+        elastic_stress = (2 * self.shear_elasticity * eps_u + 
+                    (self.bulk_elasticity - 2 * self.shear_elasticity/self.dimension) * 
+                        df.tr(eps_u) * df.Identity(self.dimension))
+        
+        viscous_stress = (2 * self.shear_viscosity * eps_v + 
+                    (self.bulk_viscosity - 2 * self.shear_viscosity/self.dimension) * 
+                        df.tr(eps_v) * df.Identity(self.dimension))
+        
+        return (elastic_stress + viscous_stress)
+    
+    def stress(self, u, v, rho):
+        return (self.passive_stress(u, v) + self.active_stress(rho))
+        
+    def reaction_rho(self, rho, trho):
+        return -self.turnover_rho * df.inner(rho * (1 - rho/self.average_rho), trho)
+    
+    def diffusion_reaction_rho(self, rho, trho):
+        return (self.diffusion_rho * df.inner(df.nabla_grad(rho), df.nabla_grad(trho)) 
+                + 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.05:
+                cell_markers[cell] = True
+        mesh = df.refine(mesh, cell_markers)
+
+        return mesh
+
+    def setup_initial_conditions(self):
+        # ic for u
+        self.zero_vector = df.Constant((0.0, 0.0))
+        u0   = df.interpolate(self.zero_vector, self.function_space.sub(0).collapse())
+        noise_u = self.noise_level * (2*np.random.random(u0.vector().size()) - 1)
+        u0.vector()[:] = u0.vector()[:] + noise_u 
+
+        # ic for rho
+        rho0 = '(a/2)*(1-tanh((sqrt(x[0]*x[0]+x[1]*x[1]) - r0)/w))'
+
+        # NOTE: For ellipse IC
+        # rho0 = '(a/2)*(1 - tanh((sqrt(x[0]*x[0]+x[1]*x[1]) - 1/(sqrt(3*sin(atan2(x[1], x[0]))*sin(atan2(x[1], x[0]))+1)) )/w))' 
+                # for ellipse IC
+
+        # NOTE: For rho=rho_0 BC
+        # rho0 = '(a/2)*(1+tanh((-x[0] - r0)/w))'
+                
+        rho0 = df.interpolate(df.Expression(rho0, a=self.average_rho, r0=0.5, w=0.1, degree=1), 
+                                            self.function_space.sub(2).collapse())
+        noise_rho = self.noise_level * (2 * np.random.random(rho0.vector().size()) - 1)
+        rho0.vector()[:] = rho0.vector()[:] + noise_rho
+        
+        VFS = self.function_space.sub(1).collapse()
+        v0  = df.Function(VFS)
+        tv = df.TestFunction(VFS)
+
+        vform = (self.friction * df.inner(v0, tv) 
+                 + df.inner(self.stress(u0, v0, rho0), self.epsilon(tv))
+                ) * df.dx
+
+        df.solve(vform == 0, v0)
+
+        df.assign(self.f0, [u0, v0, rho0])
+
+    def setup_weak_forms(self):
+        u0, v0, rho0  = df.split(self.f0)
+        u,  v, rho    = df.split(self.f)
+        tu, tv, trho = df.TestFunctions(self.function_space)
+                  
+        uform = (df.inner((u - u0)/self.timestep, tu)
+                - df.inner(v0, tu)
+                )
+    
+        vform = (self.friction * df.inner(v, tv) 
+                + df.inner(self.stress(u, v, rho), self.epsilon(tv))
+                )
+
+        rhoform = (df.inner((rho - rho0)/self.timestep, trho)
+                + self.advection(rho0, v0, trho)
+                + 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 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
+        self.time = initTime
+        savesteps = round(self.savetime/self.timestep)
+        maxsteps = round(self.maxtime/self.timestep)
+
+        if not extend:
+            # find the initial velocity
+            self.setup_initial_conditions()
+
+            # 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:
+            # assign fields to f0
+            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.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)
+        
+        self.setup_weak_forms()
+
+        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)
+
+            # save the fields
+            u, v, rho  = self.f.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)
+
+            # assign the calculated func to a newly-defined function to be used later
+            old_function = df.Function(self.function_space)
+            old_function.assign(self.f)
+
+            # move 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
+            old_function.set_allow_extrapolation(True)
+
+            self.f0 = df.project(old_function, self.function_space) 
+
+            # new weak form
+            self.setup_weak_forms()
+
+        self.uFile.close()
+        self.vFile.close()
+        self.rhoFile.close() 

u_v_rho/2D/parameters.json

+{
+    "average_rho": 1.0,
+    "bulk_elasticity": 1,
+    "bulk_viscosity": 1,
+    "diffusion_rho": 1,
+    "dimension": 2,
+    "friction": 1.0,
+    "lamda": 1,
+    "maxtime": 50.0,
+    "mesh": "mesh_radius1.00_mode0_scale0.00__5.xml.gz",
+    "noise_level": 0.0,
+    "saturation_rho": 1.0,
+    "savetime": 0.1,
+    "shear_elasticity": 1,
+    "shear_viscosity": 1,
+    "stress_setpoint_rho": 5,
+    "timestep": 0.01,
+    "turnover_rho": 1.0,
+    "zero_rho_boundary": true
+}
\ No newline at end of file

u_v_rho/2D/plot_meshes.py

+import dolfin as df
+import matplotlib.pyplot as plt
+
+# Load the mesh
+mesh = df.Mesh("mesh_radius1.00_mode2_scale0.50_asymm0.80__10.xml.gz")
+
+# Plot the mesh
+df.plot(mesh)
+# plt.title("FEniCS Mesh")
+# plt.gca().set_aspect("equal")
+plt.show()
\ No newline at end of file

u_v_rho/2D/quantify_bdry_shapes.py

+import os
+import numpy as np
+import h5py
+import dolfin as df
+import progressbar
+import matplotlib.pyplot as plt
+from scipy.spatial.distance import euclidean
+from scipy.spatial import ConvexHull
+
+def get_mesh(params, DIR=''):
+    savesteps = round(params['maxtime'] / params['savetime'])
+    times = np.arange(savesteps + 1) * params['savetime']
+
+    # Read mesh geometry from h5 file
+    var = 'density'
+    h5 = h5py.File(os.path.join(DIR, '%s_%s.h5' % (params['timestamp'], var)), "r")
+
+    # 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 = []
+    boundary_points  = []
+
+    for i in range(len(times)):
+        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()
+
+    for steps in progressbar.ProgressBar()(range(savesteps + 1)):
+        # Create a mesh for the current timestep
+        mesh = df.Mesh()
+        editor = df.MeshEditor()
+        editor.open(mesh,
+                    type='triangle' if params['dimension'] == 2 else 'tetrahedron',
+                    tdim=params['dimension'], gdim=params['dimension'])
+        editor.init_vertices(len(geometry[steps]))
+        editor.init_cells(len(topology[steps]))
+        for j in range(len(geometry[steps])):
+            editor.add_vertex(j, geometry[steps][j])
+        for j in range(len(topology[steps])):
+            editor.add_cell(j, topology[steps][j])
+        editor.close()
+
+        boundary = df.BoundaryMesh(mesh, "exterior", False)
+        boundary_coords = boundary.coordinates()
+        boundary_coords = boundary_coords[np.argsort(np.arctan2(boundary_coords[:, 1] - np.mean(boundary_coords[:, 1]),
+                                        boundary_coords[:, 0] - np.mean(boundary_coords[:, 0])))]
+        boundary_points.append(boundary_coords)
+
+    return (times, boundary_points)
+
+def quantifify_bdry(params, DIR='', offscreen=False):
+    times, bdry = get_mesh(params, DIR)
+
+
+    # Plot the initial, intermediate, final boundary 
+    plt.figure(figsize=(8, 8))
+    plt.plot(bdry[0][:, 0], bdry[0][:, 1], c='blue', label='Initial Boundary', alpha=0.7)
+    plt.plot(bdry[-1][:, 0], bdry[-1][:, 1], c='red', label='Final Boundary', alpha=0.7)
+    plt.xlabel('X')
+    plt.ylabel('Y')
+    plt.axis('equal')
+    plt.legend()
+    plt.grid(True)
+    plt.show()
+    plt.savefig(os.path.join(DIR, '%s_boundary_shapes.png' %params['timestamp']), dpi=1000)
+
+    area = []
+    aspect_ratio = []
+    perimeter = []
+    roundness = []
+    for i in range(len(times)):
+        x = bdry[i][:, 0]
+        y = bdry[i][:, 1]
+
+        # Area by shoelace formula
+        ar = 0.5 * np.abs(np.dot(x, np.roll(y, 1)) - np.dot(y, np.roll(x, 1)))
+        area.append(ar)
+        
+        # Aspect ratio by max distance to min distance from center of shape
+
+        # NOTE: The following may not work is remeshing is not suficient
+        # asp_rat = (np.max(np.sqrt((x - np.mean(x))**2 + (y - np.mean(y))**2)) 
+        #         / np.min(np.sqrt((x - np.mean(x))**2 + (y - np.mean(y))**2)))
+        # aspect_ratio.append(asp_rat)
+
+        hull = ConvexHull(bdry[i])
+        hull_points = bdry[i][hull.vertices]
+
+        x_min, x_max = np.min(hull_points[:, 0]), np.max(hull_points[:, 0])
+        y_min, y_max = np.min(hull_points[:, 1]), np.max(hull_points[:, 1])
+
+        width = x_max - x_min
+        height = y_max - y_min
+
+        asp_rat = max(width, height) / min(width, height)
+        aspect_ratio.append(asp_rat)
+
+        # Perimeter calculation
+        per = sum(euclidean(bdry[i][j], bdry[i][(j + 1) % len(bdry[i])]) for j in range(len(bdry[i])))
+        perimeter.append(per)
+        
+        # Roundness calculation
+        roundness.append((4 * np.pi * ar) / (per**2))
+
+    fig, ax = plt.subplots(1, 2, figsize=(6, 3))
+    ax[0].grid(True)
+    ax[1].grid(True)
+    ax[0].plot(times, area)
+    ax[0].set_xlabel('Time')
+    ax[0].set_ylabel('Area')
+    ax[1].plot(times, aspect_ratio)
+    ax[1].set_xlabel('Time')
+    ax[1].set_ylabel('Aspect Ratio')
+    plt.tight_layout()
+    # # ax[1, 0].plot(times, perimeter)
+    # # ax[1, 0].set_xlabel('Time')
+    # # ax[1, 0].set_ylabel('Perimeter')
+    # # ax[1, 1].plot(times, roundness)
+    # # ax[1, 1].set_xlabel('Time')
+    # # ax[1, 1].set_ylabel('Roundness')
+    plt.show()
+    plt.savefig(os.path.join(DIR, '%s_boundary_quantities.png' %params['timestamp']), dpi=1000)
+
+    fig1, ax1 = plt.subplots(1, 1, figsize=(6, 3))
+    ax1.grid(True)
+    ax1.semilogy(times, np.gradient(area), label='$dA/dt$')
+    ax1.set_xlabel('Time')
+    ax1.set_ylabel('Rate of Change of Area')
+    plt.tight_layout()
+    plt.show()
+    plt.savefig(os.path.join(DIR, '%s_boundary_area_rate.png' %params['timestamp']), dpi=1000)
+    
+if __name__ == '__main__':
+    import argparse, json
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-j', '--jsonfile', help='parameters file', required=True)
+    args = parser.parse_args()
+
+    with open(args.jsonfile) as jsonFile:
+        params = json.load(jsonFile)
+
+    quantifify_bdry(params, DIR=os.path.dirname(args.jsonfile))
\ No newline at end of file

u_v_rho/2D/run_growing_domain.py

+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 = parameters['mesh']
+
+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='')

u_v_rho/2D/viz_growing_domain.py

+import dolfin as df
+import numpy as np
+import vedo as vd
+import os
+import h5py
+import progressbar
+import matplotlib.pyplot as plt
+from tempfile import TemporaryDirectory
+
+def get_data(params, DIR=''):
+    savesteps = round(params['maxtime'] / params['savetime'])
+    times = np.arange(savesteps + 1) * params['savetime']
+
+    # Read mesh geometry from h5 file
+    var = 'density'
+    h5 = h5py.File(os.path.join(DIR, '%s_%s.h5' % (params['timestamp'], var)), "r")
+
+    # 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 = []
+
+    for i in range(len(times)):
+        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()
+
+    # Read field values
+    u = []
+    v = []
+    rho = []
+
+    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']))
+    bar = progressbar.ProgressBar(maxval=savesteps)
+
+    print('Processing meshes and field values...')
+    for steps in bar(range(savesteps + 1)):
+        # Create a mesh for the current timestep
+        mesh = df.Mesh()
+        editor = df.MeshEditor()
+        editor.open(mesh,
+                    type='triangle' if params['dimension'] == 2 else 'tetrahedron',
+                    tdim=params['dimension'], gdim=params['dimension'])
+        editor.init_vertices(len(geometry[steps]))
+        editor.init_cells(len(topology[steps]))
+        for j in range(len(geometry[steps])):
+            editor.add_vertex(j, geometry[steps][j])
+        for j in range(len(topology[steps])):
+            editor.add_cell(j, topology[steps][j])
+        editor.close()
+    
+        # Process the mesh directly without saving it
+        SFS = df.FunctionSpace(mesh, 'P', 1)
+        VFS = df.VectorFunctionSpace(mesh, '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.append(u_vec.reshape(params['dimension'], int(u_vec.shape[0] / params['dimension'])).T)
+
+        v_vec = vi.compute_vertex_values(mesh)
+        v.append(v_vec.reshape(params['dimension'], int(v_vec.shape[0] / params['dimension'])).T)
+
+        rho.append(rhoi.compute_vertex_values(mesh))
+
+    uFile.close()
+    vFile.close()
+    rhoFile.close()
+    return (times, topology, geometry, u, v, rho)
+
+
+def visualize(params, DIR='', offscreen=False):
+    (times, topology, geometry, u, v, rho) = get_data(params, DIR)
+    print(v[-1].max())
+
+    n_cmap_vals = 16
+    scalar_cmap = 'viridis'
+    cmin, cmax = min(min(sublist) for sublist in rho), max(max(sublist) for sublist in rho)
+    plotter = vd.plotter.Plotter(axes=0)
+
+    poly = vd.utils.buildPolyData(geometry[0], topology[0])
+    scalar_actor = vd.mesh.Mesh(poly)
+    scalar_actor.pointdata['density'] = rho[0]
+    scalar_actor.cmap(scalar_cmap, rho[0], vmin=cmin, vmax=cmax)
+    scalar_actor.add_scalarbar(title=r'$\rho/\rho_0$',
+                               pos=(0.8, 0.04), nlabels=2)
+
+    plotter += scalar_actor
+
+    # plot mesh
+    mesh_actor = vd.mesh.Mesh(poly, c='gray', alpha=0.3)
+    mesh_actor.wireframe()
+    plotter += mesh_actor
+
+    # plot velocity
+    # arrow_scale = 1
+    # velocity_vectors = vd.Arrows(geometry[0] - arrow_scale * v[0]/2, 
+    #                              geometry[0] + arrow_scale * v[0]/2, c='k')
+    # plotter += velocity_vectors
+
+    def update(idx):
+        nonlocal plotter, scalar_actor, mesh_actor
+        # nonlocal plotter, scalar_actor, velocity_vectors
+        poly = vd.utils.buildPolyData(geometry[idx], topology[idx])
+
+        new_scalar_actor = vd.mesh.Mesh(poly)
+        new_scalar_actor.pointdata['density'] = rho[idx]
+
+        new_scalar_actor.cmap(scalar_cmap, rho[idx], vmin=cmin, vmax=cmax)
+        new_scalar_actor.add_scalarbar(title=r'$c/c_0$',
+                                       pos=(0.8, 0.04), nlabels=2)
+        plotter.remove(scalar_actor)
+        plotter += new_scalar_actor
+        # new_velocity_vectors = vd.Arrows(geometry[idx] - arrow_scale * v[idx]/2, 
+        #                                  geometry[idx] + arrow_scale * v[idx]/2, c='k')
+        # plotter.remove(velocity_vectors)
+        # plotter += new_velocity_vectors
+
+        scalar_actor = new_scalar_actor
+        # velocity_vectors = new_velocity_vectors
+
+        new_mesh_actor = vd.mesh.Mesh(poly, c='gray', alpha=0.3)
+        new_mesh_actor.wireframe()
+        plotter.remove(mesh_actor)
+        plotter += new_mesh_actor
+        mesh_actor = new_mesh_actor
+
+
+    def slider_update(widget, event):
+        value = widget.GetRepresentation().GetValue()
+        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, mesh_actor], interactive=False, zoom=0.5)
+
+    # if offscreen:
+    FPS = 5
+    movFile = '%s_density.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 range(len(times)):
+        idx = (abs(times - times[tt])).argmin()
+        update(idx)
+        slider.GetRepresentation().SetValue(times[tt])
+        fr = get_filename("%03d.png" % tt)
+        vd.screenshot(fr)
+    os.system(command + " " + str(fps)
+              + " -i " + tmp_dir.name + os.sep
+              + "%03d.png " + options + " " + movFile)
+    tmp_dir.cleanup()
+    
+    # # radial coordinate
+    r = np.zeros(len(times))
+    for i in range(len(times)):
+        r[i] = np.max(np.sqrt(np.sum(np.square(geometry[i]), axis=1)))
+
+    print(r[0], r[-1])
+    fig, ax = plt.subplots()
+    ax.loglog(times, r)
+    ax.set_xlabel('time')
+    ax.set_ylabel('radius')
+    np.save('%s_radius.npy' % params['timestamp'], r)
+    plt.savefig('%s_radius.png' % params['timestamp'])
+
+    plt.show()
+
+    # drdt = np.gradient(r, times)
+    # fig1, ax1 = plt.subplots()
+    # ax1.loglog(times, drdt)
+    # ax1.set_xlabel('time')
+    # ax1.set_ylabel('dr/dt')
+    # plt.savefig('%s_dr_dt.png' % params['timestamp'])
+
+    # plt.show()
+
+
+if __name__ == '__main__':
+    import argparse, json
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-j', '--jsonfile', help='parameters file', required=True)
+    args = parser.parse_args()
+
+    with open(args.jsonfile) as jsonFile:
+        params = json.load(jsonFile)
+
+    visualize(params, DIR=os.path.dirname(args.jsonfile))
\ No newline at end of file