/[escript]/trunk/modellib/py_src/tsunami.py
ViewVC logotype

Annotation of /trunk/modellib/py_src/tsunami.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3432 - (hide annotations)
Fri Jan 7 01:32:07 2011 UTC (8 years, 8 months ago) by jfenwick
File MIME type: text/x-python
File size: 50558 byte(s)
Made import statements a bit more specific to clean up the epydoc
1 cochrane 263
2 ksteube 1809 ########################################################
3 ksteube 1312 #
4 jfenwick 2881 # Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1809 # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 ksteube 1312 #
8 ksteube 1809 # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11 ksteube 1312 #
12 ksteube 1809 ########################################################
13 gross 247
14 jfenwick 2881 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 ksteube 1809 Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
17     Primary Business: Queensland, Australia"""
18 elspeth 628 __license__="""Licensed under the Open Software License version 3.0
19 ksteube 1809 http://www.opensource.org/licenses/osl-3.0.php"""
20 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 elspeth 628
22 cochrane 263 import os, sys
23     import vtk
24 jfenwick 3432 #from esys.escript import *
25 gross 247 from esys.escript.linearPDEs import LinearPDE
26     from esys.escript.modelframe import Model
27 caltinay 3346 from esys.weipa import saveVTK
28 jfenwick 2455 import numpy
29 gross 247 import math
30 caltinay 1989 from tempfile import mkstemp
31 gross 247
32 caltinay 1989 WORKDIR="./work/"
33 gross 323 EPS=1.e-8
34 gross 247
35 caltinay 1989 #=============================================================================
36    
37 gross 247 class GridData:
38 caltinay 1989 """
39     This object is used to store data on a grid.
40     It will be replaced by Bruce at a later stage.
41 gross 247
42 caltinay 1989 data[i,j] are x=j*s[0]+o[0] and y=i*s[1]+o[1]
43 gross 247
44 caltinay 1989 for 0<=j<n[0] and 0<=i<n[1]
45     """
46     def __init__(self, s, o, n, data):
47     self.s=tuple(s)
48     self.o=tuple(o)
49     self.n=tuple(n)
50     self.data=data
51 gross 247
52 caltinay 1989 def getOrigin(self):
53     return self.o
54 gross 247
55 caltinay 1989 def getSpacing(self):
56     return self.s
57 gross 247
58 caltinay 1989 def getData(self):
59     return self.data
60 gross 247
61 caltinay 1989 def interpolate(self,x):
62     if hasattr(x, "getNumberOfDataPoints"):
63     x_shape0 = x.getNumberOfDataPoints()
64     return_data_object = True
65     else:
66 jfenwick 2455 x_array = numpy.array(x)
67 caltinay 1989 x_shape0 = x_array.shape[0]
68     return_data_object = False
69 gross 247
70 jfenwick 2455 data=numpy.zeros(x_shape0, numpy.float_)
71 caltinay 1989 ox,oy = self.getOrigin()
72     dx,dy = self.getSpacing()
73     data_array = self.getData()
74     i_dx = 1
75     i_dy = 1
76     for i in range(x_shape0):
77     if return_data_object:
78 jfenwick 2455 x_array = x.getTupleForDataPoint(i)
79 caltinay 1989 x_long = x_array[0]-ox
80     x_lat = x_array[1]-oy
81     else:
82     x_long = x_array[i,0]-ox
83     x_lat = x_array[i,1]-oy
84     j0 = min(max(int(x_long/dx),0), data_array.shape[1]-1-i_dy)
85     j1 = min(max(int(x_lat/dy),0), data_array.shape[0]-1-i_dx)
86     f01 = (x_long-j0*dx)/dx
87     f00 = 1.-f01
88     f11 = (x_lat-j1*dy)/dy
89     f10 = 1.-f11
90     H00 = data_array[j1,j0]
91     H01 = data_array[j1,j0+i_dx]
92     H11 = data_array[j1+i_dy,j0+i_dx]
93     H10 = data_array[j1+i_dy,j0]
94     data[i] = (H00*f00+H01*f01)*f10 + (H10*f00+H11*f01)*f11
95    
96     if return_data_object:
97     dataObj = Scalar(0, x.getFunctionSpace(), True)
98     for i in range(x_shape0):
99     dataObj.setValueOfDataPoint(i, data[i])
100     return dataObj
101     else:
102     return data
103    
104    
105 gross 247 class PointOnEarthSurface:
106 caltinay 1989 """
107     Coordinates of a point on the surface of the Earth
108     """
109     def __init__(self, long=0, lat=0):
110     self.long=long
111     self.lat=lat
112 gross 247
113 caltinay 1989 def __str__(self):
114     return "(%s,%s)"%(self.long,self.lat)
115 gross 247
116 caltinay 1989 def __sub__(self, other):
117     return self.dist(other)
118 gross 247
119 caltinay 1989 def split(self,p,t):
120     return PointOnEarthSurface(long=self.long+t*(p.long-self.long),
121     lat=self.lat+t*(p.lat-self.lat) )
122 gross 247
123 caltinay 1989 def midPoint(self,other):
124     return PointOnEarthSurface(long=(self.long+other.long)/2,
125     lat=(self.lat+other.lat)/2 )
126 gross 247
127 caltinay 1989 def dist(self,other):
128     return math.sqrt((self.long-other.long)**2+(self.lat-other.lat)**2)
129    
130    
131 gross 247 class RegionOnEarthSurface:
132 caltinay 1989 """
133     Defines a region by a south-east and north-west corner
134     """
135     RADIUS=6378.8e3
136     GRAVITY=9.81
137     WEST=0
138     NORTH=1
139     EAST=2
140     SOUTH=3
141     def __init__(self, west_south=PointOnEarthSurface(), east_north=PointOnEarthSurface(), resolution=1.):
142     if resolution <= 0:
143     raise ValueError, "resolution must be positive"
144     if west_south.long >= east_north.long:
145     raise ValueError, "south-west corner must be west of north-east corner"
146     if west_south.lat >= east_north.lat:
147     raise ValueError, "south-east corner must be south of north-west corner"
148     if east_north.lat-west_south.lat < resolution/2:
149     raise ValueError, "latitude length of region must be at least 2*larger than the resolution"
150     if east_north.long-west_south.long < resolution/2:
151     raise ValueError, "longitude length of region must be at least 2*larger than the resolution"
152 gross 247
153 caltinay 1989 self.west_south = west_south
154     self.east_north = east_north
155     self.resolution = resolution
156 gross 247
157 caltinay 1989 def __str__(self):
158     return "RegionOnEarthSurface between %s and %s" % (str(self.west_south), str(self.east_north))
159 gross 247
160 caltinay 1989 def isOnFace(self, p):
161     return self.isOnThisFace(p,self.SOUTH) or \
162     self.isOnThisFace(p,self.NORTH) or \
163     self.isOnThisFace(p,self.WEST) or \
164     self.isOnThisFace(p,self.EAST)
165    
166     def isOnThisFace(self, p, face):
167 gross 247 if face==self.WEST:
168 caltinay 1989 return self.west_south.long==p.long
169 gross 247 if face==self.EAST:
170 caltinay 1989 return p.long==self.east_north.long
171 gross 247 if face==self.SOUTH:
172 caltinay 1989 return self.west_south.lat==p.lat
173 gross 247 if face==self.NORTH:
174 caltinay 1989 return p.lat==self.east_north.lat
175 gross 247
176 caltinay 1989 def isBoxVertex(self, p):
177     return ( self.isOnThisFace(p,self.WEST) and self.isOnThisFace(p,self.SOUTH) ) or \
178     ( self.isOnThisFace(p,self.WEST) and self.isOnThisFace(p,self.NORTH) ) or \
179     ( self.isOnThisFace(p,self.EAST) and self.isOnThisFace(p,self.NORTH) ) or \
180     ( self.isOnThisFace(p,self.EAST) and self.isOnThisFace(p,self.SOUTH) )
181 gross 247
182 caltinay 1989 def getFace(self, p):
183 gross 247 # order is critical
184 caltinay 1989 if self.west_south.long == p.long: return self.WEST
185     if p.lat == self.east_north.lat: return self.NORTH
186     if p.long == self.east_north.long: return self.EAST
187     if self.west_south.lat == p.lat: return self.SOUTH
188 gross 247
189 caltinay 1989 def comparePointsOnFace(self, p0, p1):
190     f0 = self.getFace(p0)
191     f1 = self.getFace(p1)
192 gross 247
193 caltinay 1989 if f0 < f1:
194 gross 247 return -1
195 caltinay 1989 elif f0 > f1:
196 gross 247 return 1
197 caltinay 1989 else:
198 gross 247 if f0 == self.WEST:
199 caltinay 1989 if p0.lat < p1.lat:
200     return -1
201     elif p0.lat > p1.lat:
202     return 1
203 gross 247 else:
204 caltinay 1989 return 0
205 gross 247 elif f0 == self.EAST:
206 caltinay 1989 if p0.lat < p1.lat:
207     return 1
208     elif p0.lat > p1.lat:
209     return -1
210 gross 247 else:
211 caltinay 1989 return 0
212 gross 247 elif f0 == self.NORTH:
213 caltinay 1989 if p0.long < p1.long:
214     return -1
215     elif p0.long > p1.long:
216     return 1
217 gross 247 else:
218 caltinay 1989 return 0
219 gross 247 else:
220 caltinay 1989 if p0.long < p1.long:
221     return 1
222     elif p0.long > p1.long:
223     return -1
224 gross 247 else:
225 caltinay 1989 return 0
226 gross 247
227 caltinay 1989 def isInRegion(self, p):
228     return self.west_south.long <= p.long \
229     and p.long <= self.east_north.long \
230     and self.west_south.lat <= p.lat \
231     and p.lat <= self.east_north.lat
232 gross 247
233 caltinay 1989 def cutSegment(self, p0, p1):
234     t = None
235     p = None
236     tmp = self.interceptOfSegment(p0, p1, d=0, v=self.west_south.long)
237     if not tmp == None:
238     p_tmp = p0.split(p1, tmp)
239     if self.isInRegion(p_tmp) and t < tmp:
240     t = tmp
241     p = p_tmp
242 gross 247
243 caltinay 1989 tmp = self.interceptOfSegment(p0, p1, d=0, v=self.east_north.long)
244     if not tmp == None:
245     p_tmp = p0.split(p1, tmp)
246     if self.isInRegion(p_tmp) and t < tmp:
247     t = tmp
248     p = p_tmp
249 gross 247
250 caltinay 1989 tmp = self.interceptOfSegment(p0, p1, d=1, v=self.west_south.lat)
251     if not tmp == None:
252     p_tmp = p0.split(p1, tmp)
253     if self.isInRegion(p_tmp) and t < tmp:
254     t = tmp
255     p = p_tmp
256 gross 247
257 caltinay 1989 tmp = self.interceptOfSegment(p0, p1, d=1, v=self.east_north.lat)
258     if not tmp == None:
259     p_tmp = p0.split(p1, tmp)
260     if self.isInRegion(p_tmp) and t < tmp:
261     t = tmp
262     p = p_tmp
263     return p
264 gross 247
265 caltinay 1989 def interceptOfSegment(self, p0, p1, d=0, v=0.):
266     """
267     Finds t in [0,1] such that p0+t*(p0-p1) cuts x[d]=v. If it does not
268     exist None is returned.
269     """
270     if d == 0:
271     a = p0.long
272     b = p1.long
273     else:
274     a = p0.lat
275     b = p1.lat
276     if b == a:
277     if a == v:
278     t = 0
279     else:
280     t = None
281     else:
282     t = (v-a)/(b-a)
283     if not (0<=t and t<=1):
284     t = None
285     return t
286 gross 247
287    
288     class Polyline:
289 caltinay 1989 """
290     Defines a set of segments through a list of coordinates
291     """
292     def __init__(self, list_of_coordinates=[], name="none"):
293     c=[]
294     if len(list_of_coordinates)>0:
295     for i in range(len(list_of_coordinates)-1):
296     if list_of_coordinates[i]-list_of_coordinates[i+1] > 0:
297     c.append(list_of_coordinates[i])
298     c.append(list_of_coordinates[-1])
299     self.list_of_coordinates = c
300     self.name = name
301    
302     def getDiameter(self):
303     out = 0.
304     for p in self.list_of_coordinates:
305 gross 247 for q in self.list_of_coordinates:
306 caltinay 1989 out = max(out, p-q)
307     return out
308 gross 247
309 caltinay 1989 def isLoop(self):
310     if len(self) > 0:
311     return not self[0]-self[-1]>EPS
312     else:
313     return True
314 gross 247
315 caltinay 1989 def insert(self, index, coordinate):
316     """
317     Inserts point before index
318     """
319     if self.list_of_coordinates[index]-coordinate < EPS:
320     return index
321     elif self.list_of_coordinates[index-1]-coordinate < EPS:
322     return index-1
323     else:
324     self.list_of_coordinates.insert(index, coordinate)
325     return index
326 gross 247
327 caltinay 1989 def split(self, index):
328     """
329     Splits the polyline at point index
330     """
331     return Polyline(self.list_of_coordinates[:index], self.name), \
332     Polyline(self.list_of_coordinates[index:], self.name)
333 gross 247
334 caltinay 1989 def __getitem__(self, s):
335     return self.list_of_coordinates.__getitem__(s)
336    
337     def __iter__(self):
338     return iter(self.list_of_coordinates)
339    
340     def __str__(self):
341     if self.isLoop():
342     out = "loop["
343     else:
344     out = "["
345     for i in self:
346     if out[-1] == "[":
347     out += "%s" % str(i)
348 gross 247 else:
349 caltinay 1989 out += ",%s" % str(i)
350     return out+"]"
351 gross 247
352 caltinay 1989 def __len__(self):
353     return len(self.list_of_coordinates)
354 gross 247
355 caltinay 1989 def orientation(self):
356     """
357     Returns the orientation of the polyline
358     """
359     if not self.isLoop():
360     raise TypeError, "polyline is not a loop"
361 gross 247
362 caltinay 1989 integ = 0.
363     for i in range(len(self.list_of_coordinates)-1):
364     p0 = self.list_of_coordinates[i]
365     p1 = self.list_of_coordinates[i+1]
366     integ += (p1.lat-p0.lat)*(p1.long+p0.long)-(p1.long-p0.long)*(p1.lat+p0.lat)
367     if integ >= 0.:
368 gross 247 return 1.
369 caltinay 1989 else:
370 gross 247 return -1.
371    
372 caltinay 1989 def givePositiveOrientation(self):
373     if self.orientation() < 0:
374     self.list_of_coordinates.reverse()
375 gross 247
376 caltinay 1989
377 gross 247 class Coastline:
378 caltinay 1989 """
379     Defines a coast line by a Polyline within a RegionOnEarthSurface
380     """
381     def __init__(self, region, name="none"):
382     self.region = region
383     self.name = name
384     self.polylines = []
385    
386     def __str__(self):
387     out = self.name + " in " + str(self.region)
388     for pl in self.polylines:
389     out += "\n" + str(pl)
390     return out
391    
392     def makeTriangulation(self, west_south_is_water=True, east_south_is_water=True, west_north_is_water=True, east_north_is_water=True):
393     self.clean()
394     vertices=[]
395     segments=[]
396     holes=[]
397     vertices_on_face=[]
398     for pl in self.polylines:
399     if pl.getDiameter() > self.region.resolution:
400     short_pl = [pl[0]]
401     for i in range(1, len(pl)):
402     if pl[i]-short_pl[-1] > 0*EPS+self.region.resolution/10:
403     short_pl.append(pl[i])
404     elif i == len(pl)-1:
405     short_pl[-1] = pl[i]
406 gross 247 if pl.isLoop():
407 caltinay 1989 if len(short_pl) > 3:
408     offset = len(vertices)
409     v_tmp = [short_pl[0]]
410     s_tmp = []
411     for i in range(1, len(short_pl)):
412     if i == len(short_pl)-1:
413     s_tmp.append((len(v_tmp)-1+offset, offset))
414     else:
415     s_tmp.append((len(v_tmp)-1+offset, len(v_tmp)+offset))
416     v_tmp.append(short_pl[i])
417     vertices += v_tmp
418     segments += s_tmp
419     # find a point in the loop:
420     d_long = v_tmp[1].long-v_tmp[0].long
421     d_lat = v_tmp[1].lat-v_tmp[0].lat
422     l = math.sqrt(d_long**2+d_lat**2)
423     mid = v_tmp[0].midPoint(v_tmp[1])
424     n_long = -d_lat/l
425     n_lat = d_long/l
426     s = self.region.resolution
427     for seg in s_tmp:
428     p0 = vertices[seg[0]]
429     p1 = vertices[seg[1]]
430     a_long = p1.long-p0.long
431     a_lat = p1.lat-p0.lat
432     d = a_lat*n_long-a_long*n_lat
433     if abs(d) > 0.:
434     t = ((mid.lat-p0.lat)*n_long-(mid.long-p0.long)*n_lat)/d
435     if 0<=t and t<=1:
436     s_tmp = ((mid.lat-p0.lat)*a_long-(mid.long-p0.long)*a_lat)/d
437     if s_tmp > EPS:
438     s = min(s,s_tmp)
439     h = PointOnEarthSurface(long=mid.long+s/2*n_long,lat=mid.lat+s/2*n_lat)
440     holes.append(h)
441 gross 247 else:
442 caltinay 1989 if len(short_pl) > 1:
443     if self.region.isOnFace(short_pl[0]):
444     vertices_on_face.append(short_pl[0])
445     if self.region.isOnFace(short_pl[-1]):
446     vertices_on_face.append(short_pl[-1])
447     vertices.append(short_pl[0])
448     for i in range(1, len(short_pl)):
449     segments.append((len(vertices)-1, len(vertices)))
450     vertices.append(short_pl[i])
451     # put into the bounding box:
452     new_vertices = []
453     if west_south_is_water:
454     new_vertices.append(PointOnEarthSurface(long=self.region.west_south.long, lat=self.region.west_south.lat))
455     if east_south_is_water:
456     new_vertices.append(PointOnEarthSurface(long=self.region.east_north.long, lat=self.region.west_south.lat))
457     if west_north_is_water:
458     new_vertices.append(PointOnEarthSurface(long=self.region.west_south.long, lat=self.region.east_north.lat))
459     if east_north_is_water:
460     new_vertices.append(PointOnEarthSurface(long=self.region.east_north.long, lat=self.region.east_north.lat))
461 gross 247
462 caltinay 1989 # add new vertices if they don't exist yet
463     for q in new_vertices:
464     for p2 in vertices_on_face:
465     if p2-q < EPS:
466     q = None
467     raise ValueError, "coast line crosses boundary box vertex. This case is currently not supported."
468     if not q == None:
469 gross 247 vertices.append(q)
470     vertices_on_face.append(q)
471 caltinay 1989 vertices_on_face.sort(self.region.comparePointsOnFace)
472     index = 0
473     walking_on_water = True
474     l = len(vertices_on_face)
475     while index < l:
476     p1 = vertices_on_face[(index+1)%l]
477     p0 = vertices_on_face[index]
478     if walking_on_water:
479     segments.append((vertices.index(p0), vertices.index(p1)))
480     walking_on_water = False
481     else:
482     if self.region.isBoxVertex(p0):
483     segments.append((vertices.index(p0), vertices.index(p1)))
484     else:
485     walking_on_water = True
486     index += 1
487     return EarthTriangulation(vertices, segments, holes, self.region.resolution)
488 gross 247
489 caltinay 1989 def clean(self):
490 gross 247 """
491 caltinay 1989 Cleans up the coastline by joining polylines to loops or connecting
492     faces of the region
493 gross 247 """
494 caltinay 1989 # find poylines that are linked
495 gross 247 while True:
496 caltinay 1989 k0 = None
497 gross 247 for pl in self.polylines:
498 caltinay 1989 if not pl.isLoop():
499     for k in [0,-1]:
500     for ql in self.polylines:
501     if not (ql==pl or ql.isLoop()):
502     for k2 in [0,-1]:
503     if ql[k2]-pl[k] < EPS:
504     pl0 = pl
505     pl1 = ql
506     k0 = k
507     k1 = k2
508     break
509     if not k0==None: break # ql
510     if not k0==None: break # k
511     if not k0==None: break # pl
512 gross 247
513 caltinay 1989 if k0 == None:
514 gross 247 break
515     else:
516     self.polylines.remove(pl0)
517     self.polylines.remove(pl1)
518 caltinay 1989 pl0c = pl0.list_of_coordinates
519     pl1c = pl1.list_of_coordinates
520     if k0 == 0: pl0c.reverse()
521     if k1 == -1: pl1c.reverse()
522     pl = Polyline(pl0c+pl1c[1:], pl0.name+" + "+pl1.name)
523 gross 247 self.append(pl)
524    
525 caltinay 1989 # find a polyline that is not a loop and has an end or start point
526     # not on the face of the region:
527 gross 247 while True:
528 caltinay 1989 pl = None
529     k = None
530 gross 247 for pl2 in self.polylines:
531 caltinay 1989 if not pl2.isLoop():
532     pl = pl2
533     if not self.region.isOnFace(pl[0]): k=0
534     if not self.region.isOnFace(pl[-1]): k=-1
535     if not k == None: break
536     if k == None: break
537 gross 247 self.polylines.remove(pl)
538 caltinay 1989 d_min = 50000.
539     k_min = None
540 gross 247 for pl2 in self.polylines:
541 caltinay 1989 if not pl2.isLoop():
542     for k2 in [0,-1]:
543     if not self.region.isOnFace(pl2[k2]):
544     d2 = pl2[k2]-pl[k]
545     if d2 < d_min:
546     d_min = d2
547     pl_min = pl2
548     k_min = k2
549     if k_min == None:
550     raise ValueError, "cannot link coastline %s to any other coastline." % pl.name
551     plc = pl.list_of_coordinates
552     plc_min = pl_min.list_of_coordinates
553     if k == 0: plc.reverse()
554     if k_min == -1: plc_min.reverse()
555     if d_min < EPS:
556     new_pl = Polyline(plc+plc_min[1:], pl.name+" + "+pl_min.name)
557 gross 247 else:
558 caltinay 1989 new_pl = Polyline(plc+plc_min, pl.name+" + "+pl_min.name)
559 gross 247 self.polylines.remove(pl_min)
560     self.append(new_pl)
561     # give positive orientation to loops:
562     for pl in self.polylines:
563 caltinay 1989 if pl.isLoop(): pl.givePositiveOrientation()
564 gross 247
565 caltinay 1989 def append(self, polyline=Polyline()):
566     """
567     Appends a polyline
568     """
569     if len(polyline) > 1:
570     pl = []
571     outside_region = None
572     for i in range(len(polyline)):
573     if not self.region.isInRegion(polyline[i]):
574     outside_region = i
575     break
576     #pl.append(self.region.nudgeToFace(polyline[i]))
577     pl.append(polyline[i])
578     if not outside_region == None:
579     if outside_region == 0:
580     for i in range(outside_region+1, len(polyline)):
581     if self.region.isInRegion(polyline[i]):
582     polyline.insert(i,self.region.cutSegment( \
583     polyline[i-1], polyline[i]))
584     pl1 = polyline.split(i)[1]
585     self.append(pl1)
586     break
587     else:
588     # split polyline in two parts: first is fully within the
589     # region, the other starts with point outside the region
590     c = self.region.cutSegment(polyline[outside_region-1], polyline[outside_region])
591     i = polyline.insert(outside_region, c)
592     pl0,pl1 = polyline.split(i+1)
593     self.append(pl0)
594     self.append(pl1)
595     else:
596     if len(pl) > 1:
597     pply = Polyline(pl,polyline.name)
598     self.polylines.append(pply)
599 gross 247
600 caltinay 1989
601 gross 247 class EarthTriangulation:
602 caltinay 1989 """
603     Generates earth mesh and triangulates it
604     """
605     GENERATOR = "triangle -pqa%g %s"
606 gross 247
607 caltinay 1989 def __init__(self,vertices=[], segments=[], holes=[], resolution=1.):
608     (fd, self.fn) = mkstemp(suffix=".poly")
609     self.fn = self.fn.replace(".poly", "")
610     f=os.fdopen(fd, "w")
611     # write triangle input file
612     f.writelines("%d %d %d %d\n" % (len(vertices),2,0,0))
613     for i in range(len(vertices)):
614     f.writelines("%d %e %e\n" % (i,vertices[i].long,vertices[i].lat))
615     f.writelines("%d %d\n"%(len(segments),0))
616     for i in range(len(segments)):
617     f.writelines("%d %d %d\n" % (i,segments[i][0],segments[i][1]))
618     f.writelines("%d\n"%(len(holes)))
619     for i in range(len(holes)):
620     f.writelines("%d %e %e\n"%(i,holes[i].long,holes[i].lat))
621     f.close()
622 gross 247
623 caltinay 1989 # create mesh using generator tool
624     os.system(self.GENERATOR % (resolution**2, self.fn))
625 gross 247
626 caltinay 1989 # read the resulting mesh file
627     self.node_coordinates=[]
628     self.node_tags=[]
629     self.node_ids=[]
630     self.triangles_nodes=[]
631     self.triangles_id=[]
632     node_file=open("%s.1.node" % self.fn, "r")
633     nodes=node_file.readline().strip().split()
634     nn=int(nodes[0])
635     for i in range(nn):
636     nodes=node_file.readline().strip().split()
637     self.node_coordinates.append((float(nodes[1]),float(nodes[2])))
638     self.node_tags.append(int(nodes[3]))
639     self.node_ids.append(int(nodes[0]))
640     node_file.close()
641     ele_file=open("%s.1.ele" % self.fn, "r")
642     elem=ele_file.readline().strip().split()
643     ne=int(elem[0])
644     for i in range(ne):
645     elem=ele_file.readline().strip().split()
646     self.triangles_id.append(int(elem[0]))
647     self.triangles_nodes.append((int(elem[1]),int(elem[2]),int(elem[3])))
648     ele_file.close()
649 gross 247
650 caltinay 1989 # clean up
651     os.remove("%s.poly" % self.fn)
652     os.remove("%s.1.poly" % self.fn)
653     os.remove("%s.1.node" % self.fn)
654     os.remove("%s.1.ele" % self.fn)
655 gross 247
656 caltinay 1989 def getFinleyDomain(self):
657     from esys.finley import ReadMesh
658     finley_file = open("%s.msh" % self.fn, "w")
659     finley_file.writelines("%s\n2D-Nodes %d\n" % (self.fn,
660     len(self.node_coordinates)))
661     for i in range(len(self.node_coordinates)):
662     finley_file.writelines("%s %s %s %e %e\n" % (self.node_ids[i],
663     self.node_ids[i],self.node_tags[i],
664     self.node_coordinates[i][0],
665     self.node_coordinates[i][1]))
666 gross 247
667 caltinay 1989 finley_file.writelines("Tri3 %d\n" % len(self.triangles_nodes))
668     for i in range(len(self.triangles_nodes)):
669     finley_file.writelines("%s 0 %s %s %s\n" % (
670     self.triangles_id[i],
671     self.triangles_nodes[i][0],
672     self.triangles_nodes[i][1],
673     self.triangles_nodes[i][2]))
674     finley_file.writelines("Line2 0\n")
675     finley_file.writelines("Line2_Contact 0\n")
676     finley_file.writelines("Point1 0\n")
677     finley_file.close()
678     # read the mesh with finley
679     out=ReadMesh("%s.msh" % self.fn)
680     os.remove("%s.msh" % self.fn)
681     return out
682    
683    
684 gross 247 #=================================
685 caltinay 1989 # Model interfaces #
686     #=================================
687 gross 247 class OceanRegionCollector(Model):
688 caltinay 1989 """
689     Opens input streams for coastline and bathymetry data (generated by GMT)
690     """
691 gross 247
692 caltinay 1989 def __init__(self,**kwargs):
693     Model.__init__(self,**kwargs)
694     self.declareParameter(coastline_source = "coastline_default.dat",
695     bathymetry_source = "bathymetry_default.dat",
696     resolution = 1.,
697     south = 0.,
698     north = 10.,
699     east = 0.,
700     west = 20.,
701     range360 = True,
702     coastline_stream = None,
703     bathymetry_stream = None)
704 gross 247
705 caltinay 1989 def doInitialization(self):
706     """
707     Initializes the ocean region by reading coordinate files
708     """
709     # We presume that %%'s in the source name requires calling a command
710     # otherwise it is a filename containing the required data
711     if self.coastline_source.find("%%") != -1:
712     c=self.__mergeParameters(self.coastline_source)
713     self.coastline_stream=os.popen(c).read()
714     else:
715     self.coastline_stream=self.coastline_source
716 gross 247
717 caltinay 1989 if self.bathymetry_source.find("%%") != -1:
718     b=self.__mergeParameters(self.bathymetry_source)
719     self.bathymetry_stream=os.popen(b).read()
720     else:
721     self.bathymetry_stream=self.bathymetry_source
722 gross 247
723 caltinay 1989 def __mergeParameters(self,txt):
724     return txt.replace("%%west%%", str(self.west)) \
725     .replace("%%east%%", str(self.east)) \
726     .replace("%%south%%", str(self.south)) \
727     .replace("%%north%%", str(self.north)) \
728     .replace("%%resolution%%", str(self.resolution)) \
729     .replace("%%range360%%", str(self.range360).lower())
730 gross 247
731    
732 caltinay 1989 #=============================================================================
733 gross 247 class Bathymetry(Model):
734 caltinay 1989 """
735     Generates the bathymetry data within a region on the earth
736     """
737     def __init__(self,**kwargs):
738     Model.__init__(self,**kwargs)
739     self.declareParameter(source = "none",
740     bathymetry = 1.)
741 gross 247
742 caltinay 1989 def doInitialization(self):
743     """
744     Initializes the bathymetry grid
745     """
746 gross 247
747 caltinay 1989 print "Intializing bathymetry..."
748     if hasattr(self.source,"readline"):
749     f=self.source
750     else:
751     f=open(self.source,"r")
752     x_grd_list=[]
753     y_grd_list=[]
754     data_grd_list=[]
755     line=f.readline().strip()
756     while line != "":
757     v=line.split()
758     x_grd_list.append(float(v[0]))
759     y_grd_list.append(float(v[1]))
760     data_grd_list.append(float(v[2]))
761     line=f.readline().strip()
762     self.trace("%s lines have been read from %s." % (len(data_grd_list), self.source))
763 jfenwick 2455 data_grd=numpy.array(data_grd_list)
764     x_grd=numpy.array(x_grd_list)
765     y_grd=numpy.array(y_grd_list)
766 caltinay 1989 if len(x_grd)<2:
767     raise ValueError,"%s: data base is too small"%str(self)
768     ox=x_grd[0]
769     oy=y_grd[0]
770     diam=max(abs(x_grd[len(x_grd)-1]-ox),abs(y_grd[len(y_grd)-1]-oy))
771     dx=x_grd[1]-ox
772     nx=1
773     nx=1
774     while abs(x_grd[nx]-ox)>EPS*diam:
775     nx+=1
776     dy=y_grd[nx]-oy
777     ny=len(x_grd)/nx
778     data_grd.resize((ny,nx))
779     self.bathymetry=GridData(s=[dx,dy],o=[ox,oy],n=[nx,ny],data=data_grd)
780     self.trace("%sx%s grid with %sx%s spacing." % (nx,ny,dx,dy))
781 gross 247
782 caltinay 1989
783     #=============================================================================
784 gross 247 class OceanRegion(Model):
785 caltinay 1989 """
786     Generates the ocean region with a coast line and a bathymetry
787     """
788     def __init__(self,**kwargs):
789     Model.__init__(self,**kwargs)
790     self.declareParameter(domain = None,
791     resolution = 1.,
792     south = 0.,
793     north = 10.,
794     east = 0.,
795     west = 20.,
796     bathymetry = None,
797     bathymetry_data = None,
798     coastline = None,
799     source = "none")
800 gross 247
801 caltinay 1989 def doInitialization(self):
802     """
803     Initializes the ocean region
804     """
805 gross 247
806 caltinay 2067 # create output directory if necessary
807     if not os.path.isdir(WORKDIR): os.mkdir(WORKDIR)
808    
809 caltinay 1989 print "Initializing ocean region..."
810     if hasattr(self.source,"readline"):
811     f=self.source
812     data_name=f.geturl()
813     else:
814     f=open(self.source,"r")
815     data_name=self.source
816 gross 247
817 caltinay 1989 segs=[]
818     name=""
819     line=f.readline()
820     while not line=="":
821     line=line.strip()
822     if line[:7]=="#region":
823     data=line[line.index("[")+1:line.index("]")].split(",")
824     reg=RegionOnEarthSurface(PointOnEarthSurface(lat=self.south,long=self.west),PointOnEarthSurface(lat=self.north,long=self.east),self.resolution)
825     self.coastline=Coastline(region=reg,name=data_name)
826     elif line.find("Shore Bin")>-1:
827     self.coastline.append(Polyline(segs,name))
828     name=line[2:]
829     segs=[]
830     if not (line=="" or line[0]=="#" or line[0]==">") :
831     x=line.split()
832     segs.append(PointOnEarthSurface(long=float(x[0]),lat=float(x[1])))
833     line=f.readline()
834     self.coastline.append(Polyline(segs,name))
835     d=self.bathymetry_data.interpolate([[self.east,self.south],
836     [self.west,self.south],
837     [self.east,self.north],
838     [self.west,self.north]])
839     self.domain = self.coastline.makeTriangulation(
840     east_south_is_water=d[0]<=0,
841     west_south_is_water=d[1]<=0,
842     east_north_is_water=d[2]<=0,
843     west_north_is_water=d[3]<=0
844     ).getFinleyDomain()
845     self.domain.dump(os.path.join(WORKDIR, "coastline.nc"))
846     self.bathymetry = maximum(-self.bathymetry_data.interpolate(Function(self.domain).getX()),0.)
847 gross 247
848    
849 caltinay 1989 #=============================================================================
850 gross 247 class TsunamiSource(Model):
851 caltinay 1989 """
852     Defines a wave in Gaussian form between start and end.
853     """
854     GAMMA=0.05
855     def __init__(self,**kwargs):
856     Model.__init__(self,**kwargs)
857     self.declareParameter(domain = None,
858     start_lat = -10.,
859     start_long = 110.,
860     end_lat = -12.,
861     end_long = 120.,
862     width = 5.,
863     decay_zone = 0.1,
864     amplitude = 1.,
865     wave_height = 1.)
866 gross 247
867 caltinay 1989 def doInitialization(self):
868     """
869     Sets the initial wave form
870     """
871     beta=math.sqrt(-math.log(self.GAMMA))/self.decay_zone
872     x=self.domain.getX()
873     x_long = x[0]
874     x_lat = x[1]
875     mid_long = (self.start_long+self.end_long)/2
876     mid_lat = (self.start_lat+self.end_lat)/2
877     dist = math.sqrt((mid_long-self.end_long)**2+(mid_lat-self.end_lat)**2)
878     a = (self.start_long-mid_long)/dist
879     b = (self.start_lat-mid_lat)/dist
880     self.trace("source length = %s" % (dist*2.))
881     x_long_hat = a*(x_long-mid_long)+b*(x_lat-mid_lat)
882     x_lat_hat = -b*(x_long-mid_long)+a*(x_lat-mid_lat)
883     # x_lat-direction
884     s = abs(x_lat_hat)-self.width
885     m = whereNegative(s)
886     f1 = (1.-m)*exp(-(s*beta)**2)+m
887 gross 247
888 caltinay 1989 # x_long-direction
889     s = abs(x_long_hat)-dist
890     m = whereNegative(s)
891     f0 = (1.-m)*exp(-(s*beta)**2)+m
892     self.wave_height = f1*f0*self.amplitude
893    
894    
895     #=============================================================================
896 gross 247 class TsunamiInDeepWater(Model):
897 caltinay 1989 """
898     Runs the deep water tsunami model based on a simplified form of the
899     shallow water equation.
900 gross 247
901 jfenwick 2625 *d^2 h/dt^2 =div(c grad(h))*
902 gross 247
903 caltinay 1989 where h is the wave height above sea level, and c=sqrt(g*H),
904     with H - depth of the water level, g - gravity constant
905 gross 247
906 caltinay 1989 The simulation uses the Verlet scheme.
907     """
908     def __init__(self,**kwargs):
909     Model.__init__(self,**kwargs)
910     self.declareParameter(domain = None,
911     wave_height = 1.,
912     wave_velocity = 0.,
913     initial_time_step = None,
914     bathymetry = 1.,
915     safety_factor = 1.)
916 gross 247
917 caltinay 1989 def doInitialization(self):
918     """
919     Initializes the time integration scheme
920     """
921 gross 247
922 caltinay 1989 self.__pde = LinearPDE(self.domain)
923     self.__pde.setSolverMethod(self.__pde.LUMPING)
924     self.__pde.setValue(D=1.)
925 jfenwick 2455 self.__c2 = RegionOnEarthSurface.GRAVITY*self.bathymetry/(RegionOnEarthSurface.RADIUS*2*numpy.pi/360.)**2
926 caltinay 1989 c_max = math.sqrt(Lsup(self.__c2))
927     self.__dt = self.safety_factor*inf(self.domain.getSize()/(sqrt(self.__c2)+EPS*c_max))
928     if self.initial_time_step==None:
929     self.initial_time_step=self.__dt
930     self.trace("Maximum wave velocity %s m/sec" % c_max)
931     self.trace("Time step size is %s sec" % self.__dt)
932 gross 247
933 caltinay 1989 def getSafeTimeStepSize(self, dt):
934     """
935     Returns new step size
936 gross 247
937 jfenwick 2625 :param dt: last time step size used
938     :type dt: ``float``
939     :return: time step size that can safely be used
940     :rtype: ``float``
941 caltinay 1989 """
942     return self.__dt
943 gross 247
944 caltinay 1989 def doStepPostprocessing(self,dt):
945     """
946     Performs the time step using the Verlet scheme
947 gross 247
948 jfenwick 2625 :param dt: time step size to be used
949     :type dt: ``float``
950 caltinay 1989 """
951     self.__pde.setValue(X=-self.__c2*grad(self.wave_height))
952 gross 247
953 caltinay 1989 new_height = self.wave_height+dt*self.wave_velocity+dt*(self.initial_time_step+dt)/2*self.__pde.getSolution()
954 gross 247
955 caltinay 1989 self.wave_velocity = (new_height-self.wave_height)/dt
956     self.wave_height = new_height
957     self.initial_time_step = dt
958     self.trace("Wave height range is %e %e" % (inf(self.wave_height), sup(self.wave_height)))
959 gross 259
960    
961 caltinay 1989 #=============================================================================
962 gross 247 class SurfMovie(Model):
963 caltinay 1989 """
964     movie from a wave propagation on the sea
965 gross 247
966 jfenwick 2625 :ivar time: current time
967     :ivar bathymetry: scalar data set
968     :ivar wave_height: vector data set
969     :ivar filename: name of the movie file
970 caltinay 1989 """
971     def __init__(self,**kwargs):
972     Model.__init__(self,**kwargs)
973     self.declareParameter(bathymetry = 1.,
974     wave_height = 1.,
975     coastline = None,
976     t = 0.,
977     dt = 1.,
978     south = 2.,
979     north = 5.,
980     east = 3.,
981     west = 15.,
982     max_height = 1.,
983     filename = "tsunami.mpg")
984 gross 247
985 caltinay 1989 def doInitialization(self):
986     """
987     Initializes the time integration scheme
988     """
989     self.__frame_name="tsunami"
990     self.__next_t=self.dt
991     self.__fileIndex=0
992 gross 247
993 caltinay 1989 self.firstFrame = True
994     self.imageFiles = []
995 cochrane 326
996 caltinay 1989 factGraphics = vtk.vtkGraphicsFactory()
997     factGraphics.SetUseMesaClasses(1)
998     factImage = vtk.vtkImagingFactory()
999     factImage.SetUseMesaClasses(1)
1000 cochrane 326
1001 caltinay 1989 print "Preparing bathymetry data..."
1002     # get depth (z), origin, dx and dy of the bathymetry
1003     bathZData = self.bathymetry.getData()
1004     bathOrigin = self.bathymetry.getOrigin()
1005     bathSpacing = self.bathymetry.getSpacing()
1006 cochrane 263
1007 caltinay 1989 # now construct the x and y data from the spacing and the origin
1008     numXPoints = bathZData.shape[1]
1009     numYPoints = bathZData.shape[0]
1010     numPoints = numXPoints*numYPoints
1011 cochrane 263
1012 jfenwick 2455 bathXData = numpy.zeros(numXPoints, numpy.float_)
1013     bathYData = numpy.zeros(numYPoints, numpy.float_)
1014 caltinay 1989 for i in range(numXPoints):
1015     bathXData[i] = bathOrigin[0] + i*bathSpacing[0]
1016 cochrane 263
1017 caltinay 1989 for i in range(numYPoints):
1018     bathYData[i] = bathOrigin[1] + i*bathSpacing[1]
1019 cochrane 263
1020 caltinay 1989 # prepare the bathymetry data set
1021     bathPoints = vtk.vtkPoints()
1022     bathPoints.SetNumberOfPoints(numPoints)
1023     bathData = vtk.vtkFloatArray()
1024     bathData.SetNumberOfComponents(1)
1025     bathData.SetNumberOfTuples(numPoints)
1026 cochrane 263
1027 caltinay 1989 # add the points and data values
1028     count = 0
1029     for i in range(numXPoints):
1030     for j in range(numYPoints):
1031     bathPoints.InsertPoint(count, bathXData[i], bathYData[j], 0.0)
1032     bathData.InsertTuple1(count, bathZData[j,i])
1033     count += 1
1034    
1035     # create an unstructured grid for the bathymetry
1036     bathGrid = vtk.vtkUnstructuredGrid()
1037     bathGrid.SetPoints(bathPoints)
1038     bathGrid.GetPointData().SetScalars(bathData)
1039 cochrane 263
1040 caltinay 1989 # do a delaunay on the grid
1041     bathDelaunay = vtk.vtkDelaunay2D()
1042     bathDelaunay.SetInput(bathGrid)
1043     bathDelaunay.SetTolerance(0.001)
1044     bathDelaunay.SetAlpha(2.5)
1045 cochrane 263
1046 caltinay 1989 # save static bathymetry file
1047     writer = vtk.vtkPolyDataWriter()
1048     writer.SetInput(bathDelaunay.GetOutput())
1049     writer.SetScalarsName("bathymetry")
1050     writer.SetFileName(os.path.join(WORKDIR, "bathymetry.vtk"))
1051     writer.Write()
1052 cochrane 263
1053 caltinay 1989 # create the bathymetry colourmap
1054     data = []
1055     data.append([-8000, 0, 0, 0])
1056     data.append([-7000, 0, 5, 25])
1057     data.append([-6000, 0, 10, 50])
1058     data.append([-5000, 0, 80, 125])
1059     data.append([-4000, 0, 150, 200])
1060     data.append([-3000, 86, 197, 184])
1061     data.append([-2000, 172, 245, 168])
1062     data.append([-1000, 211, 250, 211])
1063     data.append([0, 16, 48, 16])
1064     data.append([300, 70, 150, 50])
1065     data.append([500, 146, 126, 60])
1066     data.append([1000, 198, 178, 80])
1067     data.append([1250, 250, 230, 100])
1068     data.append([1500, 250, 234, 126])
1069     data.append([1750, 252, 238, 152])
1070     data.append([2000, 252, 243, 177])
1071     data.append([2250, 253, 249, 216])
1072     data.append([2500, 255, 255, 255])
1073 cochrane 263
1074 caltinay 1989 # the amount to scale the data by
1075     scale = 255.0
1076     numColours = len(data)
1077 cochrane 263
1078 caltinay 1989 # convert the colourmap into something vtk is more happy with
1079 jfenwick 2455 height = numpy.zeros(numColours, numpy.float_)
1080     red = numpy.zeros(numColours, numpy.float_)
1081     green = numpy.zeros(numColours, numpy.float_)
1082     blue = numpy.zeros(numColours, numpy.float_)
1083 caltinay 1989 for i in range(numColours):
1084     (h, r, g, b) = data[i]
1085     height[i] = float(h)
1086     red[i] = float(r)/scale
1087     green[i] = float(g)/scale
1088     blue[i] = float(b)/scale
1089 cochrane 263
1090 caltinay 1989 transFunc = vtk.vtkColorTransferFunction()
1091     transFunc.SetColorSpaceToRGB()
1092     for i in range(numColours):
1093     transFunc.AddRGBPoint(height[i], red[i], green[i], blue[i])
1094     h = height[i]
1095     while i < numColours-1 and h < height[i+1]:
1096 cochrane 263 h += 1
1097     transFunc.AddRGBPoint(h, red[i], green[i], blue[i])
1098    
1099 caltinay 1989 # set up the mapper
1100     bathMapper = vtk.vtkDataSetMapper()
1101     bathMapper.SetInput(bathDelaunay.GetOutput())
1102     bathMapper.SetLookupTable(transFunc)
1103 cochrane 263
1104 caltinay 1989 # set up the actor
1105     self.bathActor = vtk.vtkActor()
1106     self.bathActor.SetMapper(bathMapper)
1107 cochrane 263
1108 caltinay 1989 print "Done preparing bathymetry data"
1109 cochrane 263
1110 caltinay 1989 ### prepare and add the coastline ###
1111 cochrane 263
1112 caltinay 1989 print "Preparing coastline data..."
1113     # create the coastline grid
1114     coastGrid = vtk.vtkUnstructuredGrid()
1115     coastPoints = vtk.vtkPoints()
1116     totalCoastPoints = 0
1117     for polyline in self.coastline.polylines:
1118     numPoints = len(polyline)
1119     coastLine = vtk.vtkPolyLine()
1120     coastLine.GetPointIds().SetNumberOfIds(numPoints)
1121     j = 0
1122     for point in polyline:
1123     coastLine.GetPointIds().SetId(j, j+totalCoastPoints)
1124     coastPoints.InsertNextPoint(point.long, point.lat, 0.0)
1125     j += 1
1126     coastGrid.InsertNextCell(coastLine.GetCellType(),
1127     coastLine.GetPointIds())
1128     totalCoastPoints += numPoints
1129 cochrane 263
1130 caltinay 1989 coastGrid.SetPoints(coastPoints)
1131     coastMapper = vtk.vtkDataSetMapper()
1132     coastMapper.SetInput(coastGrid)
1133 cochrane 263
1134 caltinay 1989 # create an actor for the coastline
1135     self.coastActor = vtk.vtkActor()
1136     self.coastActor.SetMapper(coastMapper)
1137     self.coastActor.GetProperty().SetColor(0,0,0)
1138 cochrane 263
1139 caltinay 1989 print "Done preparing coastline data"
1140 cochrane 263
1141 caltinay 1989 # set up the lookup table for the wave data
1142     refLut = vtk.vtkLookupTable()
1143     self.lutTrans = vtk.vtkLookupTable()
1144     refLut.Build()
1145     alpha = 0.7 # alpha channel value
1146     for i in range(256):
1147     (r,g,b,a) = refLut.GetTableValue(255-i)
1148     if g == 1.0 and b < 0.5 and r < 0.5:
1149     self.lutTrans.SetTableValue(i, r, g, b, 0.0)
1150     else:
1151     self.lutTrans.SetTableValue(i, r, g-0.2, b, alpha)
1152 cochrane 263
1153 caltinay 1989 # create an actor for the wave
1154     self.waveActor = vtk.vtkActor()
1155 cochrane 263
1156 caltinay 1989 # calculate window size
1157     dataWidth = max(bathXData) - min(bathXData)
1158     dataHeight = max(bathYData) - min(bathYData)
1159     self._winHeight = 600
1160     self._winWidth = int(dataWidth*float(self._winHeight)/dataHeight)
1161 cochrane 263
1162 caltinay 1989 def saveImage(self, prefix):
1163     """
1164     Saves current state as PNM image by saving a VTK file with the
1165     wave height first and reading it back afterwards.
1166     """
1167 cochrane 263
1168 caltinay 1989 # do stuff here that should only be done once
1169     if self.firstFrame:
1170     ### set up the renderer and the render window ###
1171     self.ren = vtk.vtkRenderer()
1172     self.renWin = vtk.vtkRenderWindow()
1173     self.renWin.AddRenderer(self.ren)
1174     self.renWin.SetSize(self._winWidth, self._winHeight)
1175     # add actors to the renderer
1176     self.ren.AddActor(self.bathActor)
1177     self.ren.AddActor(self.coastActor)
1178     self.ren.AddActor(self.waveActor)
1179     self.ren.GetActiveCamera().Zoom(2.0)
1180 cochrane 263
1181 caltinay 1989 # use interactor...
1182     self.iren = vtk.vtkRenderWindowInteractor()
1183     self.iren.SetRenderWindow(self.renWin)
1184     self.iren.Initialize()
1185     self.iren.Start()
1186 cochrane 263
1187 caltinay 1989 # ...or make sure rendering is offscreen
1188     #self.renWin.OffScreenRenderingOn()
1189 cochrane 263
1190 caltinay 1989 self.firstFrame = False
1191 cochrane 263
1192 caltinay 1989 vtuFile = prefix+".vtu"
1193     imgFile = prefix+".pnm"
1194 cochrane 263
1195 caltinay 1989 # save the VTK file first
1196     print "Writing", vtuFile
1197     saveVTK(vtuFile, h=self.wave_height)
1198 cochrane 263
1199 caltinay 1989 # make a reader for the data
1200     waveReader = vtk.vtkXMLUnstructuredGridReader()
1201     waveReader.SetFileName(vtuFile)
1202     waveReader.Update()
1203 cochrane 263
1204 caltinay 1989 waveGrid = waveReader.GetOutput()
1205     waveGrid.Update()
1206 cochrane 263
1207 caltinay 1989 waveMapper = vtk.vtkDataSetMapper()
1208     waveMapper.SetInput(waveGrid)
1209     waveMapper.SetLookupTable(self.lutTrans)
1210     waveMapper.SetScalarRange(-self.max_height*0.3, self.max_height*0.3)
1211 cochrane 263
1212 caltinay 1989 self.waveActor.SetMapper(waveMapper)
1213    
1214     # now render the window and save the image file
1215     print "Rendering to", imgFile
1216     self.renWin.Render()
1217    
1218     # convert the render window to an image
1219     win2imgFilter = vtk.vtkWindowToImageFilter()
1220     win2imgFilter.SetInput(self.renWin)
1221    
1222     # save the image to file
1223     outWriter = vtk.vtkPNMWriter()
1224     outWriter.SetInput(win2imgFilter.GetOutput())
1225     outWriter.SetFileName(imgFile)
1226     outWriter.Write()
1227     self.imageFiles.append(imgFile)
1228    
1229     def doStepPostprocessing(self, dt):
1230 gross 247 """
1231     Does any necessary postprocessing after each step
1232    
1233 jfenwick 2625 :param dt:
1234 gross 247 """
1235 caltinay 1989 if self.t >= self.__next_t:
1236     prefix = os.path.join(WORKDIR, \
1237     "%s%04d" % (self.__frame_name, self.__fileIndex))
1238 cochrane 263
1239 caltinay 1989 zMin = inf(self.wave_height)
1240     zMax = sup(self.wave_height)
1241     print "T = %7.1f: Wave height range=%f...%f" % (self.t, zMin, zMax)
1242 cochrane 408
1243 caltinay 1989 #self.saveImage(prefix)
1244     self.wave_height.dump(os.path.join(WORKDIR, \
1245     "waveheight%04d.nc" % self.__fileIndex))
1246 cochrane 263
1247 caltinay 1989 self.__next_t += self.dt
1248     self.__fileIndex += 1
1249 cochrane 263
1250 caltinay 1989 def getSafeTimeStepSize(self,dt):
1251     """
1252     returns new step size
1253 cochrane 263
1254 jfenwick 2625 :param dt: last time step size used
1255     :type dt: ``float``
1256     :return: time step size that can savely be used
1257     :rtype: ``float``
1258 caltinay 1989 """
1259     return self.__next_t-self.t
1260 cochrane 263
1261 caltinay 1989 def doFinalization(self):
1262     """
1263     Finalises the visualisation. For instance, makes a movie of the
1264     image files.
1265     """
1266 cochrane 263
1267 caltinay 1989 if len(self.imageFiles) == 0:
1268     return
1269 cochrane 263
1270 caltinay 1989 # set up the movie parameters
1271     paramsFileString = "REFERENCE_FRAME DECODED\n"
1272     paramsFileString += "FRAME_RATE 24\n"
1273     paramsFileString += "OUTPUT %s\n" % self.filename
1274     paramsFileString += "PATTERN IBBPBBPBB\n"
1275     paramsFileString += "FORCE_ENCODE_LAST_FRAME\n"
1276     paramsFileString += "GOP_SIZE 20\n"
1277     paramsFileString += "BSEARCH_ALG CROSS2\n"
1278     paramsFileString += "PSEARCH_ALG TWOLEVEL\n"
1279     paramsFileString += "IQSCALE 10\n"
1280     paramsFileString += "PQSCALE 11\n"
1281     paramsFileString += "BQSCALE 16\n"
1282     paramsFileString += "RANGE 8\n"
1283     paramsFileString += "SLICES_PER_FRAME 1\n"
1284     paramsFileString += "BASE_FILE_FORMAT PNM\n"
1285     paramsFileString += "INPUT_DIR .\n"
1286     paramsFileString += "INPUT_CONVERT *\n"
1287     paramsFileString += "INPUT\n"
1288     for fname in self.imageFiles:
1289     paramsFileString += "%s\n" % fname
1290     paramsFileString += "END_INPUT\n"
1291     paramsFileString += "PIXEL HALF\n"
1292     paramsFileString += "ASPECT_RATIO 1\n"
1293 cochrane 263
1294 caltinay 1989 # write the parameter file
1295     fp = open("%s.params" % self.filename, "w")
1296     fp.write(paramsFileString)
1297     fp.close()
1298 cochrane 263
1299 caltinay 1989 print "Performing conversion to mpeg"
1300     convertString = "ppmtompeg %s.params" % self.filename
1301     # write the movie into self.filename
1302     result = os.system(convertString)
1303     if result != 0:
1304     print "An error occurred in mpeg conversion"
1305 cochrane 276
1306 caltinay 1989 # now clean up the image files
1307     if result == 0:
1308     print "Removing temporary image files"
1309     os.unlink("%s.params" % self.filename)
1310     for fname in self.imageFiles:
1311     os.unlink(fname)
1312 cochrane 276
1313 cochrane 263
1314 caltinay 1989 def main():
1315     from esys.escript.modelframe import Link,Simulation
1316     from esys.modellib.input import Sequencer
1317 gross 247
1318 caltinay 1989 oc = OceanRegionCollector()
1319     oc.range360 = True
1320     oc.west = 102.9
1321     oc.east = 232.6
1322     oc.south = -71.3
1323     oc.north = 26.7
1324     oc.resolution = 0.25
1325 gross 247
1326 caltinay 1989 # Uncomment and use the following if GMT is available and the scripts work
1327     #oc.coastline_source="tsunami_create_coast.py west=%%west%% east=%%east%% south=%%south%% north=%%north%% resolution=%%resolution%% range360=%%range360%% river=false city=false citytyp=wcity.dat volcano=false hotspot=false shoreline=true state=false WSMdata=false plate=false earthquake=false"
1328     #oc.bathymetry_source = "tsunami_create_bath.py west=%%west%% east=%%east%% south=%%south%% north=%%north%% resolution=%%resolution%%"
1329 gross 247
1330 caltinay 1989 # ...otherwise use pre-created data files (range parameters are ignored)
1331     oc.coastline_source = "coastline_default.dat"
1332     oc.bathymetry_source = "bathymetry_default.dat"
1333 gross 247
1334 caltinay 1989 b = Bathymetry()
1335     b.source = Link(oc,"bathymetry_stream")
1336 gross 247
1337 caltinay 1989 oreg = OceanRegion()
1338     oreg.source = Link(oc,"coastline_stream")
1339     oreg.resolution = Link(oc, "resolution")
1340     oreg.south = Link(oc, "south")
1341     oreg.north = Link(oc, "north")
1342     oreg.east = Link(oc, "east")
1343     oreg.west = Link(oc, "west")
1344     oreg.bathymetry_data = Link(b, "bathymetry")
1345 cochrane 263
1346 caltinay 1989 src = TsunamiSource()
1347     src.domain = Link(oreg, "domain")
1348     src.decay_zone = 0.01
1349     src.end_long = 185.
1350     src.end_lat = -37.
1351     src.width = 0.5
1352     src.start_long = 174.
1353     src.start_lat = -15.
1354     src.amplitude = 3
1355 cochrane 264
1356 caltinay 1989 ts = TsunamiInDeepWater()
1357     ts.domain = Link(oreg, "domain")
1358     ts.wave_height = Link(src, "wave_height")
1359     ts.wave_velocity = 0.
1360     ts.bathymetry = Link(oreg, "bathymetry")
1361 gross 247
1362 caltinay 1989 sq = Sequencer()
1363 caltinay 2067 sq.t_end = 30000.
1364 gross 247
1365 caltinay 1989 sm = SurfMovie()
1366     sm.bathymetry = Link(b, "bathymetry")
1367     sm.wave_height = Link(ts, "wave_height")
1368     sm.coastline = Link(oreg, "coastline")
1369     sm.t = Link(sq,"t")
1370     sm.filename = "tsunami.mpg"
1371     sm.north = 8.7
1372     sm.west = 138.9
1373     sm.dt = 50.
1374     sm.east = 196.6
1375     sm.south = -53.3
1376     sm.max_height = Link(src, "amplitude")
1377 gross 247
1378 caltinay 1989 s = Simulation([sq,oc,b,oreg,src,ts,sm])
1379     #s.writeXML()
1380     s.run()
1381 cochrane 263
1382 elspeth 281 if __name__=="__main__":
1383 caltinay 1989 #from esys.modellib import tsunami
1384     main()
1385 elspeth 281
1386 cochrane 263 # vim: expandtab shiftwidth=4:

  ViewVC Help
Powered by ViewVC 1.1.26