added velocity vectors

growing_domain/full_model/viz_growing_domain.py

@@ -145,8 +145,14 @@ def visualize(params, DIR='', offscreen=False):
         # heat map
         n_cmap_vals = 16
         scalar_cmap = 'viridis'
+        vector_color = 'w'
+        vector_scale= 0.1
         geometry = np.dstack((geometry, np.zeros(geometry.shape[0:2])))
 
+        v = np.dstack((v, np.zeros(v.shape[0:2])))
+        vmag = np.linalg.norm(v, axis=2)
+        vmax = np.max(vmag)
+        v = v/vmax
         rhomin, rhomax = np.min(rho), np.max(rho)
         if params['morphogen']:
             cmin, cmax = np.min(c), np.max(c)
@@ -163,10 +169,19 @@ def visualize(params, DIR='', offscreen=False):
                     )
         plotter += scalar_actor
 
+        vector_actor = vd.shapes.Arrows(geometry[0],  geometry[0]+ vector_scale * v[0])
+        vector_actor.color(vector_color)
+        plotter+=vector_actor
         def update(idx):
+            nonlocal plotter, vector_actor
             scalar_actor.points(pts=geometry[idx], transformed=False)
             scalar_actor.pointdata['density'] = rho[idx]
-        
+            plotter.remove(vector_actor)
+            vector_actor = vd.shapes.Arrows(geometry[idx],  geometry[idx]+ vector_scale * v[idx])
+            vector_actor.color(vector_color)
+            plotter += vector_actor
+
+       
         def slider_update(widget, event):
             value = widget.GetRepresentation().GetValue()
             idx = (abs(times-value)).argmin()
@@ -176,7 +191,7 @@ def visualize(params, DIR='', offscreen=False):
                             xmin=times[0], xmax=times.max(), 
                             value=times[0], title=r"$t/\tau$")
     
-        vd.show(interactive=(not offscreen), zoom=0.2)
+        vd.show([scalar_actor, vector_actor],interactive=(not offscreen), zoom=0.8)
         if offscreen:
             FPS = 5
             movFile = '%s_cell_density.mov' % params['timestamp']
@@ -203,7 +218,7 @@ def visualize(params, DIR='', offscreen=False):
             scalar_actor1 = vd.mesh.Mesh(poly1)
             #scalar_actor1.computeNormals(points=True, cells=True)
             scalar_actor1.pointdata['concentration'] = c[0]
-            scalar_actor1.cmap(scalar_cmap, c[0], vmin=cmin, vmax=cmax, n=n_cmap_vals)
+            scalar_actor1.cmap('plasma', c[0], vmin=cmin, vmax=cmax, n=n_cmap_vals)
             scalar_actor1.add_scalarbar(title = r'$c/c_0$', 
                         pos=(0.8, 0.04), nlabels=2,
                         # titleYOffset=15, titleFontSize=28, size=(100, 600)
@@ -245,19 +260,30 @@ def visualize(params, DIR='', offscreen=False):
                 tmp_dir1.cleanup()
 
         # R(t) vs t
+        plt.rcParams.update({'font.size': 22})
+
         radius = np.array(np.zeros(len(times)))
         for j in range(len(times)):
             radius[j] = np.max(radial_coordinate[j])
+
+        figradius3, axradius3 = plt.subplots(1,1)
+        axradius3.set_xlabel(r'$\log t$')
+        # axradius3.set_xlim(np.min(np.log(times)), np.max(np.log(times)))
+        # axradius3.set_ylim(np.min(np.log(radius)), np.max(np.log(radius)))
+        axradius3.set_ylabel(r'$\log R(t)$')
+        axradius3.loglog(times, radius)
+
         figradius1, axradius1 = plt.subplots(1,1)
         axradius1.set_xlabel(r'$t$')
-        axradius1.set_xlim(np.min(times), np.max(times))
-        axradius1.set_ylim(np.min(radius), np.max(radius)+1)
+        # axradius1.set_xlim(np.min(times), np.max(times))
+        # axradius1.set_ylim(np.min(np.log(radius)), np.max(np.log(radius)))
         axradius1.set_ylabel(r'$\log R(t)$')
         axradius1.semilogy(times, radius)
+
         figradius2, axradius2 = plt.subplots(1,1)
         axradius2.set_xlabel(r'$t$')
         axradius2.set_xlim(np.min(times), np.max(times))
-        # axradius2.set_ylim(np.min(radius), np.max(radius))
+        axradius2.set_ylim(np.min(radius), np.max(radius))
         axradius2.set_ylabel(r'$R(t)$')
         axradius2.plot(times, radius)
         plt.show()