imposed boundary stress, no activity, what are the response curves

u_v_rho/1D/check_response_curves_without_activity/parameters.json

+{
+    "average_rho": 1.0,
+    "diffusion_rho": 0.1,
+    "elasticity": 1,
+    "friction": 1.0,
+    "lamda": 2.0,
+    "maxtime": 100.0,
+    "noise_level": 0.0,
+    "resolution": 128,
+    "savetime": 0.1,
+    "system_size": 1.0,
+    "timestep": 0.01,
+    "turnover_rho": 1,
+    "viscosity": 1,
+    "applied_boundary_stress": 2.0,
+    "zero_rho_boundary": true
+}
\ No newline at end of file

u_v_rho/1D/check_response_curves_without_activity/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/check_response_curves_without_activity/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 Left(df.SubDomain):
+    def __init__(self, x_left, tol):
+        super().__init__()  # call init of parent class
+        self.x_left = x_left
+        self.tol = tol
+
+    def inside(self, x, on_boundary):
+        return on_boundary and df.near(x[0], self.x_left, self.tol)
+
+class Right(df.SubDomain):
+    def __init__(self, x_right, tol):
+        super().__init__()
+        self.x_right =  x_right
+        self.tol = tol
+
+    def inside(self, x, on_boundary):
+        return on_boundary and df.near(x[0], self.x_right, self.tol)
+
+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)
+
+        # Define boundary markers
+        boundaries = df.MeshFunction("size_t", self.mesh, self.mesh.topology().dim() - 1)          # size_t: non-negative integer values will be stored
+                                                                                                   # mesh.topology().dim() - 1: the facets (edges in 2D, faces in 3D)
+        boundaries.set_all(0)
+        tol = 1E-14
+        xmin = self.mesh.coordinates().min()
+        xmax = self.mesh.coordinates().max()
+
+        left = Left(xmin, tol)
+        right = Right(xmax, tol)
+
+        left.mark(boundaries, 1)
+        right.mark(boundaries, 2)
+        
+        self.ds = df.Measure("ds", domain = self.mesh, subdomain_data = boundaries)
+
+    def imposed_boundary_stress_normal_component(self):
+        return df.Constant((self.applied_boundary_stress, ))
+
+    def advection(self, conc, vel, tconc):
+        return df.inner(df.div(vel * conc), tconc)
+
+    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 * 2 * eps_u 
+        viscous_stress = self.viscosity * 2 * eps_v
+        return (elastic_stress + viscous_stress)
+
+    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.passive_stress(u0, v0), 
+                               self.epsilon(tv0))
+                ) * df.dx
+                + df.inner(tv0, self.imposed_boundary_stress_normal_component()) * self.ds(1)
+                - df.inner(tv0, self.imposed_boundary_stress_normal_component()) * self.ds(2)
+        )
+        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))* df.dx
+        
+        vform = ((self.friction * df.inner(v, tv) 
+                + df.inner(self.passive_stress(u, v), 
+                           self.epsilon(tv))) * df.dx
+                + df.inner(tv, self.imposed_boundary_stress_normal_component()) * self.ds(1)
+                - df.inner(tv, self.imposed_boundary_stress_normal_component()) * self.ds(2))
+
+        rhoform = (df.inner((rho - rho0)/self.timestep, trho)
+                + self.advection(rho0, v0, trho)
+                + self.diffusion_reaction_rho(rho, trho))* df.dx
+        
+        self.form = (uform + vform + rhoform)
+
+    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()  # Builds space on moved, unrefined mesh
+            self.function0 = df.project(self.function0, self.function_space)
+
+
+        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/check_response_curves_without_activity/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
+from scipy.interpolate import interp1d
+
+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 = []
+
+    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))
+
+        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]
+
+    # 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_geometry = min(min(sublist) for sublist in geometry)
+    max_geometry = max(max(sublist) for sublist in 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, marker = 'o')
+    velplot, = axes[1].plot(geometry[0], v[0], ms=3, color='r', lw = 2, marker = 'o')
+    uplot, = axes[2].plot(geometry[0], u[0], ms=3, color='b', lw = 2, marker = 'o')
+    
+    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()
+
+    
+
+    # # # 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])
+
+    # #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(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(params['elasticity'], params['lamda'], dldt[-1])
+
+    # print(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'])
+    np.save('%s_dldt.npy' %params['timestamp'], dldt)
+
+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))