1-D plotting with matplotlib

growing_domain/viz_growing_domains.py

@@ -3,6 +3,9 @@ import numpy as np
 import vedo as vd
 import os
 import h5py
+from matplotlib.widgets import Slider
+
+import matplotlib.pyplot as plt
 from tempfile import TemporaryDirectory
 
 def get_data(params, DIR=''):
@@ -19,7 +22,9 @@ def get_data(params, DIR=''):
         geometry.append(np.array(h5['%s/%s_%d/mesh/geometry'%(var,var,i)]))
     h5.close()
     geometry = np.array(geometry)
-    
+    if params['dimension']==1:
+        geometry, zeros= np.dsplit(geometry, 2)
+ 
     # create a mesh
     if params['dimension']==1:
         mesh = df.IntervalMesh(params['resolution'], 0, params['system_size'])
@@ -42,7 +47,8 @@ def get_data(params, DIR=''):
     cFile = os.path.join(DIR, '%s_%s.xdmf' % (params['timestamp'],var))
     with df.XDMFFile(cFile) as cFile:
         for steps in range(savesteps+1):
-            mesh.coordinates()[:] = geometry[steps]
+            if params['dimension']==1:
+                mesh.coordinates()[:] = geometry[steps]
             VFS = df.FunctionSpace(mesh, 'P', 1)
             c = df.Function(VFS)
             cFile.read_checkpoint(c, var, steps)
@@ -51,41 +57,90 @@ def get_data(params, DIR=''):
     return (times, geometry, topology, concentration)
 
 def visualize(params, DIR='', offscreen=False):
-    n_cmap_vals = 16
-    scalar_cmap = 'viridis'
-    
     (times, geometry, topology, concentration) = get_data(params, DIR)
-    if params['dimension']==2:
+
+    radius = np.linalg.norm(geometry, axis=2)
+
+    if params['dimension']==1:
+        fig, axc = plt.subplots(1,1, figsize=(8,8))
+        axc.set_xlabel(r'$x$')
+        axc.set_xlim(np.min(radius), np.max(radius))
+        axc.set_ylim(np.min(concentration), np.max(concentration))
+        axc.set_ylabel(r'$c(x,t)$')
+        cplot, = axc.plot(radius[0], concentration[0])
+
+        def update(value):
+            ti = np.abs(times-value).argmin()
+            cplot.set_ydata(concentration[ti])
+            cplot.set_xdata(radius[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()
+
+    else:
+        # heat map
+        n_cmap_vals = 16
+        scalar_cmap = 'viridis'
         geometry = np.dstack((geometry, np.zeros(geometry.shape[0:2])))
 
-    cmin, cmax = np.min(concentration), np.max(concentration)
-    
-    plotter = vd.plotter.Plotter(axes=0)
-    
-    poly = vd.utils.buildPolyData(geometry[0], topology)
-    scalar_actor = vd.mesh.Mesh(poly)
-    #scalar_actor.computeNormals(points=True, cells=True)
-    scalar_actor.pointdata['concentration'] = concentration[0]
-    scalar_actor.cmap(scalar_cmap, concentration[0], vmin=cmin, vmax=cmax, n=n_cmap_vals)
-    scalar_actor.addScalarBar(title = r'$c$', 
-                pos=(0.8, 0.04), nlabels=2,
-                titleYOffset=15, titleFontSize=28, size=(100, 600))
-    plotter += scalar_actor
-    
-    def update(idx):
-        scalar_actor.points(pts=geometry[idx], transformed=False)
-        scalar_actor.pointdata['concentration'] = concentration[idx]
+        cmin, cmax = np.min(concentration), np.max(concentration)
+        
+        plotter = vd.plotter.Plotter(axes=0)
+        
+        poly = vd.utils.buildPolyData(geometry[0], topology)
+        scalar_actor = vd.mesh.Mesh(poly)
+        #scalar_actor.computeNormals(points=True, cells=True)
+        scalar_actor.pointdata['concentration'] = concentration[0]
+        scalar_actor.cmap(scalar_cmap, concentration[0], vmin=cmin, vmax=cmax, n=n_cmap_vals)
+        scalar_actor.add_scalarbar(title = r'$c$', 
+                    pos=(0.8, 0.04), nlabels=2,
+                    # titleYOffset=15, titleFontSize=28, size=(100, 600)
+                    )
+        plotter += scalar_actor
+
+        def update(idx):
+            scalar_actor.points(pts=geometry[idx], transformed=False)
+            scalar_actor.pointdata['concentration'] = concentration[idx]
         
-    def slider_update(widget, event):
-        value = widget.GetRepresentation().GetValue()
-        idx = (abs(times-value)).argmin()
-        update(idx)
+        def slider_update(widget, event):
+            value = widget.GetRepresentation().GetValue()
+            idx = (abs(times-value)).argmin()
+            update(idx)
         
-    slider = plotter.addSlider2D(slider_update, pos=[(0.1,0.94), (0.5,0.94)],
+        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(interactive=(not offscreen), zoom=0.8)
+        vd.show(interactive=(not offscreen), zoom=0.8)
+
+        
+        # c(r,t) vs r
+        fig1, axc1 = plt.subplots(1,1, figsize=(8,8))
+        axc1.set_xlabel(r'$r$')
+        axc1.set_xlim(np.min(radius), np.max(radius))
+        axc1.set_ylim(np.min(concentration), np.max(concentration))
+        axc1.set_ylabel(r'$c(r,t)$')
+        cplot, = axc1.plot(radius[0], concentration[0],'o', ms=3)
+
+        def update(value):
+            ti = np.abs(times-value).argmin()
+            cplot.set_ydata(concentration[ti])
+            cplot.set_xdata(radius[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()
     
     # make movie
     if offscreen:
@@ -119,4 +174,3 @@ if __name__ == '__main__':
         params = json.load(jsonFile)
 
     visualize(params, DIR=os.path.dirname(args.jsonfile), offscreen=(not args.onscreen))
-    
\ No newline at end of file