Removed mshr dependence, superfluous viz code, added standard disk/sphere meshes

growing_domain/diffusion_on_growing_domain.py

 import dolfin as df
-import mshr as ms
 import numpy as np
 import progressbar
 import os
@@ -28,11 +27,9 @@ class GrowthDiffusion(object):
         if self.dimension==1:
             self.mesh = df.IntervalMesh(self.resolution, 0, self.system_size)
         elif self.dimension==2: 
-            geometry = ms.Circle(df.Point(0, 0), self.system_size)
-            self.mesh = ms.generate_mesh(geometry, self.resolution)     
+            self.mesh = df.Mesh('disk.xml.gz')
         elif self.dimension==3:
-            geometry = ms.Sphere(df.Point(0, 0, 0), self.system_size)
-            self.mesh = ms.generate_mesh(geometry, self.resolution)
+            self.mesh = df.Mesh('sphere.xml.gz')
 
         # create mesh, function space, define function, test function
         self.SFS = df.FunctionSpace(self.mesh, 'P', 1)
@@ -91,6 +88,8 @@ if __name__ == '__main__':
     # parse parameters
     assert params['dimension'] in (1,2,3)
     assert params['growth'] in ('none','linear','exponential'), 'Unknown growth model'
+    if params['dimension'] in (2,3):
+        del params['resolution']
     
     timestamp = datetime.datetime.now().strftime("%d%m%y-%H%M%S")
     params['timestamp'] = timestamp

growing_domain/viz_growing_domains.py

 import dolfin as df
 import numpy as np
 import vedo as vd
-import matplotlib.pyplot as plt
-from mpl_toolkits.axes_grid1.inset_locator import inset_axes
-import progressbar
 import os
-import json
 import h5py
 from tempfile import TemporaryDirectory
 
-
 def get_data(params, DIR=''):
     savesteps = int(params['maxtime']/params['savetime'])
     times = np.arange(savesteps+1) * params['savetime']    
@@ -77,7 +72,6 @@ def visualize(params, DIR='', offscreen=False):
                 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]
@@ -90,7 +84,6 @@ def visualize(params, DIR='', offscreen=False):
     slider = plotter.addSlider2D(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)
     
@@ -113,171 +106,7 @@ def visualize(params, DIR='', offscreen=False):
                     + " -i " + tmp_dir.name + os.sep
                     + "%03d.png " + options + " " + movFile)
         tmp_dir.cleanup()
-    
-def visualize1(params, DIR='', offscreen=False):
-    n_cmap_vals = 10
-    scalar_cmap = 'viridis'
-    vector_color = '#ffffff'
-    vector_scale = 0.1
-    savesteps = int(params['maxtime']/params['savetime'])
-    times = params['alpha'] * np.arange(savesteps+1) * params['savetime']
-    halftime = int(len(times)/2)
-
-    mesh = get_mesh_from_xdmf(DIR, params)
-
-    print('Reading polarity')
-    polarity = np.zeros((len(times), mesh.num_vertices(), 3))
-    polarity_functions = []
-    VFS = df.VectorFunctionSpace(mesh, 'P', 1)
-    p = df.Function(VFS)
-    with df.XDMFFile(os.path.join(DIR,'%s_polarity.xdmf' % params['timestamp'])) as pFile:
-        for steps in progressbar.progressbar(range(savesteps+1)):
-            pFile.read_checkpoint(p, 'polarity', steps)
-            polarity_functions.append(p)
-            vector_values = p.compute_vertex_values(mesh)
-            polarity[steps] = vector_values.reshape(3, int(vector_values.shape[0]/3)).T
-    pmag = np.linalg.norm(polarity, axis=2)
-    pmin, pmax = 0, np.max(pmag)
-    polarity = polarity / pmax
-
-
-    # load previously computed convergence data
-    data = np.loadtxt(os.path.join(DIR, '%s.dat' % params['timestamp']))
-
-    fig1, ax = plt.subplots(1)
-    fig1.subplots_adjust(left=0.15, bottom=0.15, right=0.9, top=0.85)
-    ax.plot(times, data[:,1])
-    ax.set_xlabel(r'$\alpha t$')
-    ax.set_ylabel(r'$\int_{\Gamma} (\mathbf{p}_{t+\Delta t}-\mathbf{p}_{t})^2$')
-    axins = inset_axes(ax, width="50%", height="40%")
-    axins.plot(times[halftime:], data[halftime:,1])
-    fig1.savefig(os.path.join(DIR, '%s_pnorm.pdf' % params['timestamp']), transparent=True)
-
-
-    VSH_labels = [r'$\mathbf{\Psi}_{10}$',
-                r'$\mathbf{\Phi}_{10}$',
-                r'$\mathbf{\Psi}_{11}$',
-                r'$\mathbf{\Phi}_{11}$',
-                r'$\mathbf{\Psi}_{20}$',
-                r'$\mathbf{\Phi}_{20}$',
-                r'$\mathbf{\Psi}_{21}$',
-                r'$\mathbf{\Phi}_{21}$',
-                r'$\mathbf{\Psi}_{22}$',
-                r'$\mathbf{\Phi}_{22}$'
-    ]
-
-    fig2, ax = plt.subplots(1)
-    fig2.subplots_adjust(left=0.15, bottom=0.15, right=0.9, top=0.85)
-    ax.set_xlim(times[halftime:].min(), times.max())
-    for i in range(len(VSH_labels)):
-        ax.plot(times[halftime:], data[halftime:, i+2], label=VSH_labels[i], lw=1)
-    ax.set_xlabel(r'$\alpha t$')
-    ax.set_ylabel(r'$\sqrt{p^{\Psi}_{lm} (p^{\Psi}_{lm})^{\ast}}, \ \sqrt{p^{\Phi}_{lm} (p^{\Phi}_{lm})^{\ast}}$')
-    ax.legend(loc=0, ncol=5, fontsize='small', bbox_to_anchor=(0.1, 1))
-    fig2.savefig(os.path.join(DIR, '%s_modes.pdf' % params['timestamp']), transparent=True)
-    if not offscreen:
-        fig1.show()
-        fig2.show()
-
-
-    (_, _, _, _, normal, _, _) = compute_geometric_quantities(params['meshfile'])
-    projector = df.Identity(3) - df.outer(normal, normal)
-    TFS = df.TensorFunctionSpace(mesh, 'P', 1)
-    p = polarity_functions[-1]
-    grad_p = df.project(projector * df.grad(p) * projector, TFS)
-    div_p = df.tr(grad_p)
-    omega = df.skew(grad_p)
-    Theta = df.assemble(df.inner(div_p, div_p) * df.dx(mesh))
-    Phi   = df.assemble(df.inner(omega, omega) * df.dx(mesh))
-
-    defect_locations, defect_eigvals = find_defects(mesh, p, grad_p, pmag[-1])
-    assert len(defect_locations)==2
-
-    params['Theta_norm'] = Theta
-    params['Phi_norm'] = Phi
-
-    # clean up
-    for key in ["converged", "dtpnorm_ave"]:
-        if key in params.keys():
-            del params[key]
-            
-    with open(os.path.join(DIR, params['timestamp'] + '_polarity.json'), "w") as fp:
-        json.dump(params, fp, indent=4)
-    
-
-    plotter = vd.plotter.Plotter(axes=4)
-    poly = vd.utils.buildPolyData(mesh.coordinates(), mesh.cells())
-    scalar_actor = vd.mesh.Mesh(poly)
-    scalar_actor.computeNormals(points=True, cells=True)
-    
-    scalar_actor.cmap(scalar_cmap, pmag[0], vmin=pmin, vmax=pmax, n=n_cmap_vals)
-    scalar_actor.addScalarBar(
-        title = r'$\sqrt{\left(\frac{\beta}{\alpha}\right)} \; \left\vert\mathbf{p}\right\vert$', 
-        pos=(0.8, 0.1), nlabels=2,
-        titleYOffset=15, titleFontSize=28, size=(100, 600))
-    scalar_actor.scalarbar.SetLabelFormat("%0.1f")
-    plotter.add(scalar_actor)
-    
-    defects = vd.shapes.Spheres(defect_locations, c='red', r=0.02)
-    plotter.add(defects)
-    vec = defect_locations[1] - defect_locations[0]
-    uvec = vec / np.linalg.norm(vec)
-    line = vd.Line(-1.1*uvec, 1.1*uvec, lw=6, c='black')
-    plotter.add(line)
-
-    endPoints = mesh.coordinates() + vector_scale * polarity[0]
-    vector_actor = vd.shapes.Arrows(mesh.coordinates(), endPoints)
-    vector_actor.color(vector_color)
-    plotter.add(vector_actor)
-    def update(idx):
-        nonlocal vector_actor
-        scalar_actor.pointdata['polarity'] = pmag[idx]
-        scalar_actor.cmap(scalar_cmap, pmag[idx], vmin=pmin, vmax=pmax, n=n_cmap_vals)
-        plotter.remove(vector_actor)
-        endPoints = mesh.coordinates() + vector_scale * polarity[idx]
-        vector_actor = vd.shapes.Arrows(mesh.coordinates(), endPoints)
-        vector_actor.color(vector_color)
-        plotter.add(vector_actor)
-    def slider_update(widget, event):
-        value = widget.GetRepresentation().GetValue()
-        idx = (abs(times-value)).argmin()
-        update(idx)
-        
-    slider = plotter.addSlider2D(slider_update, pos=[(0.25,0.05), (0.75,0.05)],
-                            xmin=times[0], xmax=times.max(), 
-                            value=times[0], title=r"$\alpha \, t$")
-    slider.GetRepresentation().SetLabelFormat('%0.1f')
-    
-    # optimized camera position
-    cam = dict(pos=(4.952, 3.107, 1.804),
-           focalPoint=(0.07537, 0.1585, -0.07462),
-           viewup=(-0.2810, -0.1402, 0.9494),
-           distance=6.000,
-           clippingRange=(2.836, 10.26))
-
-    vd.show(interactive=(not offscreen), camera=cam, zoom=1.4)
 
-
-    # make movie
-    if offscreen:
-        print('Saving movie...')
-        FPS = 10
-        movFile = os.path.join(DIR, '%s_polarity.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))):
-            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
@@ -288,7 +117,6 @@ if __name__ == '__main__':
 
     with open(args.jsonfile) as jsonFile:
         params = json.load(jsonFile)
-        
-    (times, geometry, topology, concentration) = get_data(params, DIR='')
+
     visualize(params, DIR=os.path.dirname(args.jsonfile), offscreen=(not args.onscreen))
     
\ No newline at end of file