turing patterns on a growing domain

growing_domain/growing_domain_patterns/viz_growing_domain_patterns.py

+import dolfin as df
+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=''):
+    savesteps = int(params['maxtime']/params['savetime'])
+    times = np.arange(savesteps+1) * params['savetime']    
+
+    # Read mesh geometry from h5 file
+    var = 'c1'
+    h5 = h5py.File(os.path.join(DIR, '%s_%s.h5' % (params['timestamp'], var)), "r")
+    # should be in the loop if remeshing
+    topology = np.array(h5['%s/%s_0/mesh/topology'%(var,var)])
+    geometry = []
+    for i in range(len(times)):
+        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'])
+    else:
+        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[0]))
+        editor.init_cells(len(topology))
+        for i in range(len(geometry[0])):
+            editor.add_vertex(i, geometry[0][i])
+        for i in range(len(topology)):
+            editor.add_cell(i, topology[i])
+        editor.close()
+    
+    # Read concentration data
+    concentration = np.zeros((len(times), len(geometry[0])))
+    cFile = os.path.join(DIR, '%s_%s.xdmf' % (params['timestamp'],var))
+    with df.XDMFFile(cFile) as cFile:
+        for steps in range(savesteps+1):
+            if params['dimension']==1:
+                mesh.coordinates()[:] = geometry[steps]
+            VFS = df.FunctionSpace(mesh, 'P', 1)
+            c = df.Function(VFS)
+            cFile.read_checkpoint(c, var, steps)
+            concentration[steps] = c.compute_vertex_values(mesh)
+
+    return (times, geometry, topology, concentration)
+
+def visualize(params, DIR='', offscreen=False):
+    (times, geometry, topology, concentration) = get_data(params, DIR)
+
+    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.add_scalarbar(title = r'$c_{numerical}$', 
+                    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)
+        
+        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)
+        if offscreen:
+            FPS = 10
+            movFile = '%s_numerical.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.io.screenshot(fr)
+            os.system(command + " " + str(fps)
+                        + " -i " + tmp_dir.name + os.sep
+                        + "%03d.png " + options + " " + movFile)
+            tmp_dir.cleanup()
+               
+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))