with morphogen and nematic, most minimal model

u_v_c_Q/compare_with_expt.py

+import numpy as np
+import matplotlib.pyplot as plt
+plt.rcParams['font.size'] = 18
+savesteps = round(30/0.1)
+times = np.arange(savesteps+1) * 0.1
+l1=np.load('060725-161644_radius.npy')
+l2 =np.load('060725-163030_radius.npy')
+plt.plot(times, l1/l1[0], lw = 2, label='More alignment', color='black')
+plt.plot(times, l2/l2[0], lw = 2, label='Less alignment', color='orange')
+plt.xlabel(r'$\kappa t$')
+plt.ylabel(r'$L(t)/L(0)$')
+plt.legend()
+plt.savefig('compare_with_expt_lol.png', dpi=300, bbox_inches='tight')
+plt.show()
\ No newline at end of file

u_v_c_Q/growing_domain.py

+import dolfin as df
+import numpy as np
+import progressbar
+import meshzoo
+from ufl.constantvalue import PermutationSymbol as epsilon
+
+df.set_log_level(df.LogLevel.ERROR)
+df.parameters['form_compiler']['optimize'] = True
+
+class Growing_Domain(object):
+    def __init__(self, parameters, mesh = None):
+        print('running the eqns with q tensor..')
+
+        # 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:
+            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 = mesh
+
+            
+    def setup(self):
+        scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
+        vector_element = df.VectorElement('P', self.mesh.ufl_cell(), 1)
+        tensor_element = df.TensorElement('P', self.mesh.ufl_cell(), 1)
+
+        # u, v, c, Q, l, L
+        mixed_element = df.MixedElement([vector_element, vector_element, 
+                                         scalar_element, tensor_element, 
+                                         scalar_element, scalar_element, 
+                                        ])
+
+        # define function space with this mixed element
+        self.function_space = df.FunctionSpace(self.mesh, mixed_element)
+
+        # 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 active_stress(self, c, q):               
+        return self.lamda * (c/(c + self.saturation_c)) * q   
+        # return self.lamda * q       
+    
+    def passive_stress(self, u, v, q, c):
+        eps_u, eps_v = df.sym(df.nabla_grad(u)), df.sym(df.nabla_grad(v))
+        
+        elastic_stress = (self.shear_elasticity * 2 * eps_u + 
+                    (self.bulk_elasticity - self.shear_elasticity) * 
+                        df.tr(eps_u) * df.Identity(2))
+        
+        viscous_stress = (self.shear_viscosity * 2 * eps_v + 
+                    (self.bulk_viscosity - self.shear_viscosity) * 
+                        df.tr(eps_v) * df.Identity(2))
+        
+
+        nematic_stress = (df.dot(q, self.H_polynomial_terms(q, c) - self.Kq * df.div(df.nabla_grad(q))) 
+                          - df.dot(self.H_polynomial_terms(q, c) - self.Kq * df.div(df.nabla_grad(q)), q) 
+                          - self.stretch  * (self.H_polynomial_terms(q, c) - self.Kq * df.div(df.nabla_grad(q)))
+                        )        
+        
+        if self.neglect_orientational_stress:
+            # print('neglecting orientational stress')
+            return (elastic_stress + viscous_stress)
+        else:
+            # print('including orientational stress')
+            return (elastic_stress + viscous_stress + nematic_stress)
+    
+    def stress(self, u, v, c, q):
+        return (self.passive_stress(u, v, q, c) + self.active_stress(c, q))
+
+    def H_polynomial_terms(self, q, c):
+        return (self.A * q * c**2
+                + self.B * df.dot(q, q) 
+                + self.C * q * df.inner(q, q) 
+                + self.D * df.dot(df.dot(q, q), q))
+        
+   
+    def refine_mesh(self, mesh):
+        cell_markers = df.MeshFunction("bool", mesh, mesh.topology().dim())
+        cell_markers.set_all(False)
+        for cell in df.cells(mesh):
+            if cell.volume() > 2 * 0.005:
+                cell_markers[cell] = True
+        mesh = df.refine(mesh, cell_markers)
+
+        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_u * (2*np.random.random(u0.vector().size())-1)
+        u0.vector()[:] = u0.vector()[:] + noise_u 
+
+        # ic for c
+        c0 = 'b * (1 + tanh((x[0] - a)/w))'
+        if self.c_scaling_profile:
+            c0 = df.interpolate(df.Expression(c0, b = self.c_amplitude, 
+                                              a =  self.c_origin * max([np.sqrt(vertex.x(0)**2 + vertex.x(1)**2) 
+                                                    for vertex in df.vertices(self.mesh)]), 
+                                            w = self.c_width, degree=1), 
+                                            self.function_space.sub(2).collapse())
+        else:
+            c0 = df.interpolate(df.Expression(c0, b = self.c_amplitude, 
+                                              a = self.c_origin, w =self.c_width, degree=1), 
+                                         self.function_space.sub(2).collapse()) 
+
+        # ic for q
+        a = np.random.randn(2, 2)
+        a = (a + a.T)/2
+        a -= np.trace(a) * np.eye(2)/2
+        a = self.noise_level_q * a
+
+        #NOTE: If random initial conditions for q
+        q0 = df.Constant(((0.01 , 0.1), (0.1, -0.01)))
+
+        #NOTE: Polarised initial conditions for q
+        # q0 = 0.01 * (df.outer(df.nabla_grad(c0), df.nabla_grad(c0)) - df.Identity(2) * df.tr(df.outer(df.nabla_grad(c0), df.nabla_grad(c0))) / 2
+        # )
+        # q0_noise = df.Constant(((1 , 0), (0, -1)))
+        # q0 = df.project(q0+ q0_noise, self.function_space.sub(3).collapse())
+        q0 = df.project(q0, self.function_space.sub(3).collapse())
+
+        l0 = df.interpolate(df.Constant(0), self.function_space.sub(4).collapse())
+        L0 = df.interpolate(df.Constant(0), self.function_space.sub(5).collapse())
+
+        VFS = df.VectorFunctionSpace(self.mesh, 'P', 1)
+        vi = df.Function(VFS)
+        tvi = df.TestFunction(VFS) 
+
+        viform = (self.friction * df.inner(vi, tvi) 
+                    + df.inner(self.stress(u0, vi, c0, q0), df.nabla_grad(tvi))
+                ) * df.dx
+  
+        df.solve(viform == 0, vi)
+
+        df.assign(self.f0, [u0, vi, c0, q0, l0, L0])
+
+    def antisymm_gradv(self, v):
+        return df.skew(df.nabla_grad(v))
+
+    def symm_gradv(self, v):
+        return df.sym(df.nabla_grad(v))
+        
+    def corotation_stretch_term(self, v, q):
+        return (  df.dot(self.antisymm_gradv(v), q)
+                - df.dot(q, self.antisymm_gradv(v))
+                #NOTE: sign of following imp
+                - self.stretch * df.sqrt(df.inner(q, q)) * (self.symm_gradv(v) - df.div(v) * df.Identity(2)/2)
+                
+             )                
+
+    def setup_weak_forms(self):
+        u0, v0, c0, q0, l0, L0 = df.split(self.f0)
+        u, v, c, q, l, L = df.split(self.f)
+        tu, tv, tc, tq, tl, tL = df.TestFunctions(self.function_space)
+                  
+        uform = (df.inner((u - u0)/self.timestep, tu)
+                - df.inner(v, tu)
+                ) * df.dx
+        
+        vform = (self.friction * df.inner(v, tv) 
+                + df.inner(self.stress(u, v, c, q), df.nabla_grad(tv))
+                ) * df.dx
+
+        cform = (df.inner((c - c0)/self.timestep, tc)
+                #  + df.inner(df.dot(v0, df.nabla_grad(c0)), tc)
+                ) * df.dx
+
+        qform = (df.inner((q - q0)/self.timestep, tq)
+                 + df.inner(df.dot(v0, df.nabla_grad(q0)), tq) 
+                 + df.inner(self.corotation_stretch_term(v0, q0), tq)
+                 + self.mobility * df.inner(self.H_polynomial_terms(q, c), tq)
+                 + self.mobility * self.Kq * df.inner(df.nabla_grad(q), df.nabla_grad(tq))
+                 + l * df.inner(df.Identity(2), tq)
+                 + L * df.inner(epsilon(2), tq)) * df.dx           
+        
+        lform = df.inner(tl, df.tr(q)) * df.dx
+
+        Lform = df.inner(tL, df.inner(epsilon(2), q)) * df.dx
+
+                
+        self.form = (uform + vform + qform + cform + lform + Lform )
+                    
+        
+    def solve(self, initu=None, initq = 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.cFile = df.XDMFFile('%s_c.xdmf' % self.timestamp)
+        self.qFile = df.XDMFFile('%s_q.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, c0, q0, l0, L0  = self.f0.split(deepcopy=True)
+            self.uFile.write_checkpoint(u0, 'displacement', self.time)
+            self.vFile.write_checkpoint(v0, 'velocity', self.time)
+            self.cFile.write_checkpoint(c0, 'c', self.time)
+            self.qFile.write_checkpoint(q0, 'q', 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())
+            q0 = df.project(initq, self.function_space.sub(3).collapse())
+            l0 = df.Constant(0.0)
+            L0 = df.Constant(0.0)
+
+            df.assign(self.f0, [u0, v0, c0, q0, l0, L0])
+
+            # 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)
+
+            # save the fields
+            u, v, c, q, l, L  = 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.cFile.write_checkpoint(c, 'c', self.time, append=True)
+                self.qFile.write_checkpoint(q, 'q', 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.cFile.close()
+        self.qFile.close()

u_v_c_Q/parameters.json

+{
+    "A": -0.7,
+    "B": 0,
+    "C": 1,
+    "D": 0,
+    "Kq": 1,
+    "bulk_elasticity": 1,
+    "bulk_viscosity": 3,
+    "c_amplitude": 1,
+    "c_origin": 0.5,
+    "c_scaling_profile": false,
+    "c_width": 0.5,
+    "dim": 2,
+    "friction": 1,
+    "lamda": -2,
+    "maxtime": 30.0,
+    "mobility": 1.0,
+    "neglect_orientational_stress": false,
+    "noise_level_q": 0.01,
+    "noise_level_u": 0.0,
+    "saturation_c": 1.0,
+    "savetime": 0.1,
+    "shear_elasticity": 0.5,
+    "shear_viscosity": 1,
+    "stretch": 0,
+    "timestep": 0.01
+}
\ No newline at end of file

u_v_c_Q/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_c_Q/run_growing_domain.py

+import json, datetime
+import os
+from growing_domain import Growing_Domain
+from viz_growing_domain_2 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)
+parser.add_argument('--onscreen', action=argparse.BooleanOptionalAction)
+
+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
+    initTime = 0.0
+    initq = None
+    mesh = None
+
+else:
+    extend = True
+    print('Extending run %s...' % parameters['timestamp'])
+    parametersFile = args.jsonfile
+
+    savesteps = round(parameters['maxtime']/parameters['savetime'])
+
+    # read geometry and topology at the last timepoint from h5 file
+    var = 'density'
+    h5 = h5py.File( '%s_%s.h5' % (parameters['timestamp'], var,), 'r+')
+    geometry = np.array(h5['%s/%s_%d/mesh/geometry'%(var,var,savesteps)])
+    topology = np.array(h5['%s/%s_%d/mesh/topology'%(var,var,savesteps)])
+
+    # create mesh with this geometry and topology
+    mesh = df.Mesh()
+    editor = df.MeshEditor()
+    editor.open(mesh, 
+                type='triangle' if parameters['dimension']==2 else 'tetrahedron', 
+                tdim=parameters['dimension'], gdim=parameters['dimension'])
+    editor.init_vertices(len(geometry))
+    editor.init_cells(len(topology))
+    for j in range(len(geometry)):
+        editor.add_vertex(j, geometry[j])
+    for j in range(len(topology)):
+        editor.add_cell(j, topology[j])
+
+    editor.close()
+    # Read field values at the last time point
+    SFS = df.FunctionSpace(mesh, 'P', 1)
+    VFS = df.VectorFunctionSpace(mesh, 'P', 1)
+
+    initu, initv, initrho = df.Function(VFS), df.Function(VFS), df.Function(SFS)
+    uFile   = df.XDMFFile('%s_displacement.xdmf' % parameters['timestamp'])
+    vFile   = df.XDMFFile('%s_velocity.xdmf' % parameters['timestamp'])
+    rhoFile = df.XDMFFile('%s_density.xdmf' % parameters['timestamp'])
+    uFile.read_checkpoint(initu, 'displacement', savesteps)
+    vFile.read_checkpoint(initv, 'velocity', savesteps)
+    rhoFile.read_checkpoint(initrho, 'density', savesteps)
+
+    uFile.close()
+    vFile.close()
+    rhoFile.close()
+    initTime = parameters['maxtime']
+    parameters['maxtime'] = args.time
+
+tissue = Growing_Domain(parameters, mesh)
+tissue.solve(initu, initq, 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='', offscreen=(not args.onscreen))

u_v_c_Q/viz_growing_domain _old.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 = 'c'
+    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_c.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, 'c', 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

u_v_c_Q/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 = 'c'
+    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
+    v = []
+    conc = []
+
+    vFile = df.XDMFFile(os.path.join(DIR, '%s_velocity.xdmf' % params['timestamp']))
+    cFile = df.XDMFFile(os.path.join(DIR, '%s_c.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['dim'] == 2 else 'tetrahedron',
+                    tdim=params['dim'], gdim=params['dim'])
+        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)
+        vi, ci = df.Function(VFS), df.Function(SFS)
+
+        vFile.read_checkpoint(vi, 'velocity', steps)
+        cFile.read_checkpoint(ci, 'c', steps)
+
+        v_vec = vi.compute_vertex_values(mesh)
+        v.append(v_vec.reshape(params['dim'], int(v_vec.shape[0] / params['dim'])).T)
+
+        conc.append(ci.compute_vertex_values(mesh))
+
+    vFile.close()
+    return (times, topology, geometry, v, conc)
+
+
+def visualize(params, DIR='', offscreen=False):
+    (times, topology, geometry, v, conc) = get_data(params, DIR)
+
+    n_cmap_vals = 16
+    scalar_cmap = 'Greens'
+    cmin, cmax = min(min(sublist) for sublist in conc), max(max(sublist) for sublist in conc)
+    plotter = vd.plotter.Plotter(axes=0)
+
+    poly = vd.utils.buildPolyData(geometry[0], topology[0])
+    scalar_actor = vd.mesh.Mesh(poly)
+    scalar_actor.pointdata['c'] = conc[0]
+    scalar_actor.cmap(scalar_cmap, conc[0], vmin=cmin, vmax=cmax)
+    # scalar_actor.add_scalarbar(title=r'$c/c_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, velocity_vectors
+        poly = vd.utils.buildPolyData(geometry[idx], topology[idx])
+
+        new_scalar_actor = vd.mesh.Mesh(poly)
+        new_scalar_actor.pointdata['density'] = conc[idx]
+
+        new_scalar_actor.cmap(scalar_cmap, conc[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
+
+
+    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, velocity_vectors], interactive=False, zoom=1)
+
+    # if offscreen:
+    FPS = 20
+    movFile = '%s_v.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()
+
+    def plot_velocity_at_time(params, DIR='', target_time=0.9):
+        (times, topology, geometry, v, conc) = get_data(params, DIR)
+
+        # Find the index corresponding to the target time
+        idx = (abs(times - target_time)).argmin()
+
+        n_cmap_vals = 16
+        scalar_cmap = 'Greens'
+        cmin, cmax = min(min(sublist) for sublist in conc), max(max(sublist) for sublist in conc)
+        plotter = vd.plotter.Plotter(axes=0)
+
+        # Build the polydata and scalar actor for the target time
+        poly = vd.utils.buildPolyData(geometry[idx], topology[idx])
+        scalar_actor = vd.mesh.Mesh(poly)
+        scalar_actor.pointdata['density'] = conc[idx]
+        scalar_actor.cmap(scalar_cmap, conc[idx], vmin=cmin, vmax=cmax)
+
+        # Add scalar actor to the plotter
+        plotter += scalar_actor
+
+        # Plot velocity vectors for the target time
+        arrow_scale = 1
+        velocity_vectors = vd.Arrows(geometry[idx] - arrow_scale * v[idx]/2, 
+                                    geometry[idx] + arrow_scale * v[idx]/2, c='k')
+        plotter += velocity_vectors
+
+        # Show the plot
+        vd.show([scalar_actor, velocity_vectors], interactive=True, zoom=1)
+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_c_Q/viz_growing_domain_2.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 = 'c'
+    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 = []
+    q = []
+
+    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_c.xdmf' % params['timestamp']))
+    qFile  = df.XDMFFile(os.path.join(DIR, '%s_q.xdmf' % params['timestamp']))
+    bar = progressbar.ProgressBar(maxval=savesteps)
+
+    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['dim'] == 2 else 'tetrahedron',
+                    tdim=params['dim'], gdim=params['dim'])
+        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)
+        TFS = df.TensorFunctionSpace(mesh, 'P', 1)
+        ui, vi, rhoi, qi = df.Function(VFS), df.Function(VFS), df.Function(SFS), df.Function(TFS)
+
+        uFile.read_checkpoint(ui, 'displacement', steps)
+        vFile.read_checkpoint(vi, 'velocity', steps)
+        rhoFile.read_checkpoint(rhoi, 'c', steps)
+        qFile.read_checkpoint(qi, 'q', steps)
+        
+        u_vec = ui.compute_vertex_values(mesh)
+        u.append(u_vec.reshape(2, int(u_vec.shape[0]/2)).T)
+
+        v_vec = vi.compute_vertex_values(mesh)
+        v.append(v_vec.reshape(2, int(v_vec.shape[0]/2)).T)
+
+        rho.append(rhoi.compute_vertex_values(mesh))
+
+        qvals = qi.compute_vertex_values(mesh)
+        qvals = qvals.reshape(2, 2, int(qvals.shape[0]/2**2)).T
+        q.append(qvals.transpose(0, 2, 1))
+
+    uFile.close()
+    vFile.close()
+    rhoFile.close()
+    qFile.close()
+    return (times, topology, geometry, u, v, rho, q)
+
+
+def visualize(params, DIR='', offscreen=False):
+    (times, topology, geometry, u, v, rho, q) = get_data(params, DIR)
+    # for i in range(len(times)):
+        # print(f"AS part of q at time {times[i]}:", q[i][:, 0, 1] - q[i][:, 1, 0])
+        # print(f"Trace of q at time {times[i]}:", q[i][:, 0, 0] + q[i][:, 1, 1])
+    qmag = []
+    director = []
+    for i in range(len(times)):
+        qmag.append(np.sqrt(np.einsum('ijk, ijk -> i', q[i], q[i])))
+        dir = np.zeros((len(qmag[i]), 2))
+        for j in range(len(qmag[i])):
+            l, v = np.linalg.eig(q[i][j])
+            idx = np.argmax(l)
+            dir[j] = np.real(v[:, idx]) * qmag[i][j]
+        director.append(dir)
+
+    qmin, qmax = min(min(sublist) for sublist in qmag), max(max(sublist) for sublist in qmag)
+    rhomin, rhomax = min(min(sublist) for sublist in rho), max(max(sublist) for sublist in rho)
+
+    def create_movie(scalar_field, scalar_title, scalar_min, scalar_max, filename):
+        plotter = vd.plotter.Plotter(axes=0)
+
+        poly = vd.utils.buildPolyData(geometry[0], topology[0])
+        scalar_actor = vd.mesh.Mesh(poly)
+        if scalar_field is rho:
+            scalar_actor.pointdata['c'] = rho[0]
+            scalar_actor.cmap('Greens', rho[0], vmin=rhomin, vmax=rhomax)
+            scalar_actor.add_scalarbar( 
+                        pos=[[0.8, 0.04], [0.9, 0.4]], nlabels=2,
+
+                        # titleYOffset=15, titleFontSize=28, size=(100, 600)
+                        )
+        else:
+            scalar_actor.pointdata['q'] = qmag[0]
+            scalar_actor.cmap('Greens', qmag[0], vmin=qmin, vmax=qmax)
+            scalar_actor.add_scalarbar(title = r'$|Q|$', 
+                        pos=[[0.8, 0.04], [0.9, 0.4]], nlabels=2,
+                        # titleYOffset=15, titleFontSize=28, size=(100, 600)
+                        )
+        plotter += scalar_actor
+
+        scale =  0.1
+        points = geometry[0]
+        startpoints = points - scale * director[0]/2
+        endpoints = points + scale * director[0]/2
+        lines = vd.shapes.Lines(startpoints, endpoints, lw=3, c='red')
+        plotter += lines
+
+        def update(idx):
+            nonlocal plotter, scalar_actor, lines
+            poly = vd.utils.buildPolyData(geometry[idx], topology[idx])
+
+            new_scalar_actor = vd.mesh.Mesh(poly)
+            if scalar_field is q:
+                new_scalar_actor.pointdata['q'] = qmag[idx]
+                new_scalar_actor.cmap('Greens', qmag[idx], vmin=qmin, vmax=qmax)
+                new_scalar_actor.add_scalarbar(title = r'$|Q|$', 
+                        pos=[[0.8, 0.04], [0.9, 0.4]], nlabels=2,
+                        # titleYOffset=15, titleFontSize=28, size=(100, 600)
+                        )
+            else:
+                new_scalar_actor.pointdata['c'] = rho[idx] 
+                new_scalar_actor.cmap('Greens', rho[idx], vmin=rhomin, vmax=rhomax)  
+                new_scalar_actor.add_scalarbar( 
+                        pos=[[0.8, 0.04], [0.9, 0.4]], nlabels=2,
+                        # titleYOffset=15, titleFontSize=28, size=(100, 600)
+                        ) 
+
+            
+            plotter.remove(scalar_actor)
+            plotter += new_scalar_actor
+
+            scalar_actor = new_scalar_actor
+
+            points = geometry[idx]
+            startpoints = points - scale * director[idx]/2
+            endpoints = points + scale * director[idx]/2
+            new_lines = vd.shapes.Lines(startpoints, endpoints, lw=3, c='red')
+            plotter.remove(lines)
+            plotter += new_lines
+
+            lines = new_lines
+
+        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$")
+        
+        slider.GetRepresentation().SetLabelFormat('%0.1f')
+
+        vd.show([scalar_actor], interactive=False, zoom=0.6)
+        if offscreen:
+            if scalar_field is q:
+                FPS = 20
+                movFile = '%s_q_on_qmag.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()
+            else:
+                FPS = 20
+                movFile = '%s_q_on_c.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()
+    create_movie(q, r'$|Q|$', qmin, qmax, '%s_q_on_qmag.mov' % params['timestamp'])
+    create_movie(rho, r'$c/c_0$', rhomin, rhomax, '%s_q_on_c.mov' % params['timestamp'])
+    
+    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.plot(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()
+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:
+        params = json.load(jsonFile)
+
+    visualize(params, DIR=os.path.dirname(args.jsonfile), offscreen=(not args.onscreen))

u_v_c_Q/viz_growing_domain_old.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