MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
plot.py
1 #!/usr/bin/env python
2 # -*- coding: utf-8 -*-
3 
8 from __future__ import division,unicode_literals
9 import numpy
10 from math import radians,cos,sin,tan,pi,sqrt
11 import matplotlib.pyplot as plt
12 from mpl_toolkits.mplot3d import Axes3D
13 from matplotlib import cm
14 import matplotlib.animation as animation
15 path='./'
16 plots=[1,1,0,1]
17 saveplots=[0,0,0,0]
18 fontsize=25
19 linewidth=4
20 data=numpy.loadtxt(path+'filmthickness.csv')
21 time,dr,np,vt,fs,hmin,hloc,hcenter,havg=\
22  numpy.loadtxt(path+'results_1d.csv',skiprows=1,unpack=True)
23 points=int(np[0]) #discretization points
24 radius=[]
25 for i in range(points):
26  radius.append(dr[0]*(0.5+i))
27 xran=len(radius)
28 radius=numpy.array(radius)
29 if plots[0]:
30  nlines=5
31  xpart=1
32  xran=int(len(radius)*xpart)
33  fig = plt.figure(figsize=(10.0,8.0))
34  plt.rc('xtick', labelsize=fontsize)
35  plt.rc('ytick', labelsize=fontsize)
36  for i in range(nlines):
37  j=int((len(time)-1)/(nlines-1)*i)
38  for k in range(xran):
39  radius[k]=dr[j]*(0.5+k)
40  line1,=plt.plot(radius[0:xran]*1e6,data[j][0:xran]*1e6,
41  label='t={0:g} s'.format(time[j]),lw=linewidth)
42 # plt.ylim(0,1)
43 # plt.xlim(0,40)
44  plt.legend(loc=2, fontsize=fontsize)
45  plt.xlabel('Distance from center (μm)', fontsize=fontsize)
46  plt.ylabel('Film half thickness (μm)', fontsize=fontsize)
47  if saveplots[0]:
48  plt.savefig('profiles.pdf')
49  plt.show()
50 
51 if plots[1]:
52  colors=['b','g','m']
53  linewidth=2
54  times=[0,int(len(data)/10),len(data)-1]
55  fig = plt.figure(figsize=(10.0,10.0))
56  plt.rc('xtick', labelsize=fontsize)
57  plt.rc('ytick', labelsize=fontsize)
58  j=numpy.argmax(dr)
59  for k in range(xran):
60  radius[k]=dr[j]*(0.5+k)
61  Rc=radius[-1]+dr[j]/2
62  plt.xlim(-Rc/sqrt(3.0)*1e6,Rc*1e6)
63  plt.ylim(-Rc*1e6,Rc*1e6)
64  for i,j in enumerate(times):
65  xran=len(radius)
66  for k in range(xran):
67  radius[k]=dr[j]*(0.5+k)
68  Rc=radius[-1]
69  x0=radius[0:xran]
70  y0=data[times[i]][0:xran]
71  a=tan(pi/2)
72  c=0
73  d=(x0+a*(y0-c))/(1+a**2)
74  x1=2*d-x0+Rc+data[times[i]][-1]*tan(pi/6)
75  y1=2*d*a-y0+2*c
76  a=tan(pi/3)
77  c=0
78  d=(x1+a*(y1-c))/(1+a**2)
79  x2=2*d-x1
80  y2=2*d*a-y1+2*c
81  x2=x2[::-1]
82  y2=y2[::-1]
83  x1=numpy.append(x1,x2)
84  y1=numpy.append(y1,y2)
85  if i==0:
86  line1,=plt.plot(x1*1e6,y1*1e6,color=colors[i],
87  label='t={0:g} s'.format(time[times[i]]),lw=linewidth)
88  else:
89  line1,=plt.plot(x1*1e6,y1*1e6,color=colors[i],lw=linewidth)
90  a=tan(-pi/3)
91  c=0
92  d=(x1+a*(y1-c))/(1+a**2)
93  x2=2*d-x1
94  y2=2*d*a-y1+2*c
95  if i==1:
96  line2,=plt.plot(x2*1e6,y2*1e6,color=colors[i],
97  label='t={0:g} s'.format(time[times[i]]),lw=linewidth)
98  else:
99  line2,=plt.plot(x2*1e6,y2*1e6,color=colors[i],lw=linewidth)
100  a=tan(0)
101  c=0
102  d=(x1+a*(y1-c))/(1+a**2)
103  x3=2*d-x1
104  y3=2*d*a-y1+2*c
105  if i==2:
106  line3,=plt.plot(x3*1e6,y3*1e6,color=colors[i],
107  label='t={0:g} s'.format(time[times[i]]),lw=linewidth)
108  else:
109  line3,=plt.plot(x3*1e6,y3*1e6,color=colors[i],lw=linewidth)
110  plt.legend(loc=1, fontsize=fontsize)
111  plt.xlabel('Dimensions in (μm)', fontsize=fontsize)
112  plt.ylabel('Dimensions in (μm)', fontsize=fontsize)
113  if saveplots[1]:
114  plt.savefig('profiles2d.pdf')
115  plt.show()
116 
117 if plots[2]: #not to scale when bubble grows
118  fig = plt.figure(figsize=(5.0,4.0))
119  ax = fig.gca(projection='3d')
120  radius3d, time3d = numpy.meshgrid(radius, time)
121  surf = ax.plot_surface(time3d, radius3d*1e6, data*1e6, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
122  #surf = ax.plot_surface(time3d, radius3d, data, rstride=2, cstride=2, cmap=cm.coolwarm, linewidth=0, antialiased=False)
123  # ax.set_xlabel('time')
124  # ax.set_ylabel('radius (um)')
125  # ax.set_zlabel('thickness (um)')
126  # plt.rc('xtick', labelsize=10)
127  # plt.rc('ytick', labelsize=10)
128  plt.show()
129  if saveplots[2]:
130  plt.savefig('profiles3d.png')
131 
132 if plots[3]:
133  j=numpy.argmax(dr)
134  for k in range(xran):
135  radius[k]=dr[j]*(0.5+k)
136  fig,ax = plt.subplots(figsize=(10.0,8.0))
137  ax.set_xlim(radius[0]*1e6, (radius[-1]+dr[0]/2)*1e6)
138  ax.set_ylim(0, max(map(max,data))*1e6)
139  line,=ax.plot(radius[0:xran]*1e6,data[0][0:xran]*1e6,
140  label='t={0:.1e} s'.format(time[0]),lw=linewidth)
141  plt.rc('xtick', labelsize=fontsize)
142  plt.rc('ytick', labelsize=fontsize)
143  plt.legend(loc=2, fontsize=fontsize)
144  plt.xlabel('position', fontsize=fontsize)
145  plt.ylabel('half thickness (um)', fontsize=fontsize)
146  iii=[]
147  iii.append(-1)
148 
149 def update(x):
150  iii[0]=iii[0]+1
151  if iii[0]>len(time)-1:
152  iii[0]=0
153  plt.legend(['t={0:.1e} s'.format(time[iii[0]])], loc=2, fontsize=fontsize)
154  for k in range(xran):
155  radius[k]=dr[iii[0]]*(0.5+k)
156  line.set_xdata(radius[0:xran]*1e6)
157  line.set_ydata(data[iii[0]][0:xran]*1e6)
158  return line,
159 
160 def data_gen():
161  while True: yield data[0][0:xran]*1e6
162 
163 if plots[3]:
164  try:
165  ani=animation.FuncAnimation(fig,update,data_gen,interval=100)
166  plt.show()
167  except:
168  pass
169  if saveplots[3]:
170  print 'saving of animation not implemented'