Skip to content

Jigyasa Watwani / growth-pattern-control

Go to a project
Commit 8b2e5a8c by Vijay Kumar Krishnamurthy
Options

Visualization script. Need to get profiles in 2D/3D

parent 8ac363b5
Inline Side-by-side
Showing with 295 additions and 0 deletions
growing_domain/viz_growing_domains.py 0 → 100644
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']
# Read mesh geometry from h5 file
var = 'concentration'
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)
# 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):
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):
n_cmap_vals = 16
scalar_cmap = 'viridis'
(times, geometry, topology, concentration) = get_data(params, DIR)
if params['dimension']==2:
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]
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)],
xmin=times[0], xmax=times.max(),
value=times[0], title=r"$t/\tau$")
vd.show(interactive=(not offscreen), zoom=0.8)
# make movie
if offscreen:
FPS = 10
movFile = '%s.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()
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
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)
(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
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or sign in to comment