/[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 1989 - (hide annotations)
Fri Nov 7 04:11:33 2008 UTC (11 years ago) by caltinay
File MIME type: text/x-python
File size: 50580 byte(s)
modellib: completely revised tsunami script.
-updated finley usage (due to changed interface)
-moved VTK stuff to one place making it easy to disable (and use dumps instead)
-removed dependency on non-existing data server - script now depends on saved coastline/bathymetry data files OR calls scripts to create the data on the fly.
-PDE and triangulation code is unchanged except indentation updates

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

  ViewVC Help
Powered by ViewVC 1.1.26