/[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 2881 - (hide annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 9 months ago) by jfenwick
File MIME type: text/x-python
File size: 50526 byte(s)
Don't panic.
Updating copyright stamps

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 gross 247 from esys.escript import *
25     from esys.escript.linearPDEs import LinearPDE
26     from esys.escript.modelframe import Model
27 jfenwick 2455 import numpy
28 gross 247 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 jfenwick 2455 x_array = numpy.array(x)
66 caltinay 1989 x_shape0 = x_array.shape[0]
67     return_data_object = False
68 gross 247
69 jfenwick 2455 data=numpy.zeros(x_shape0, numpy.float_)
70 caltinay 1989 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 jfenwick 2455 x_array = x.getTupleForDataPoint(i)
78 caltinay 1989 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 jfenwick 2455 data_grd=numpy.array(data_grd_list)
763     x_grd=numpy.array(x_grd_list)
764     y_grd=numpy.array(y_grd_list)
765 caltinay 1989 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 2067 # create output directory if necessary
806     if not os.path.isdir(WORKDIR): os.mkdir(WORKDIR)
807    
808 caltinay 1989 print "Initializing ocean region..."
809     if hasattr(self.source,"readline"):
810     f=self.source
811     data_name=f.geturl()
812     else:
813     f=open(self.source,"r")
814     data_name=self.source
815 gross 247
816 caltinay 1989 segs=[]
817     name=""
818     line=f.readline()
819     while not line=="":
820     line=line.strip()
821     if line[:7]=="#region":
822     data=line[line.index("[")+1:line.index("]")].split(",")
823     reg=RegionOnEarthSurface(PointOnEarthSurface(lat=self.south,long=self.west),PointOnEarthSurface(lat=self.north,long=self.east),self.resolution)
824     self.coastline=Coastline(region=reg,name=data_name)
825     elif line.find("Shore Bin")>-1:
826     self.coastline.append(Polyline(segs,name))
827     name=line[2:]
828     segs=[]
829     if not (line=="" or line[0]=="#" or line[0]==">") :
830     x=line.split()
831     segs.append(PointOnEarthSurface(long=float(x[0]),lat=float(x[1])))
832     line=f.readline()
833     self.coastline.append(Polyline(segs,name))
834     d=self.bathymetry_data.interpolate([[self.east,self.south],
835     [self.west,self.south],
836     [self.east,self.north],
837     [self.west,self.north]])
838     self.domain = self.coastline.makeTriangulation(
839     east_south_is_water=d[0]<=0,
840     west_south_is_water=d[1]<=0,
841     east_north_is_water=d[2]<=0,
842     west_north_is_water=d[3]<=0
843     ).getFinleyDomain()
844     self.domain.dump(os.path.join(WORKDIR, "coastline.nc"))
845     self.bathymetry = maximum(-self.bathymetry_data.interpolate(Function(self.domain).getX()),0.)
846 gross 247
847    
848 caltinay 1989 #=============================================================================
849 gross 247 class TsunamiSource(Model):
850 caltinay 1989 """
851     Defines a wave in Gaussian form between start and end.
852     """
853     GAMMA=0.05
854     def __init__(self,**kwargs):
855     Model.__init__(self,**kwargs)
856     self.declareParameter(domain = None,
857     start_lat = -10.,
858     start_long = 110.,
859     end_lat = -12.,
860     end_long = 120.,
861     width = 5.,
862     decay_zone = 0.1,
863     amplitude = 1.,
864     wave_height = 1.)
865 gross 247
866 caltinay 1989 def doInitialization(self):
867     """
868     Sets the initial wave form
869     """
870     beta=math.sqrt(-math.log(self.GAMMA))/self.decay_zone
871     x=self.domain.getX()
872     x_long = x[0]
873     x_lat = x[1]
874     mid_long = (self.start_long+self.end_long)/2
875     mid_lat = (self.start_lat+self.end_lat)/2
876     dist = math.sqrt((mid_long-self.end_long)**2+(mid_lat-self.end_lat)**2)
877     a = (self.start_long-mid_long)/dist
878     b = (self.start_lat-mid_lat)/dist
879     self.trace("source length = %s" % (dist*2.))
880     x_long_hat = a*(x_long-mid_long)+b*(x_lat-mid_lat)
881     x_lat_hat = -b*(x_long-mid_long)+a*(x_lat-mid_lat)
882     # x_lat-direction
883     s = abs(x_lat_hat)-self.width
884     m = whereNegative(s)
885     f1 = (1.-m)*exp(-(s*beta)**2)+m
886 gross 247
887 caltinay 1989 # x_long-direction
888     s = abs(x_long_hat)-dist
889     m = whereNegative(s)
890     f0 = (1.-m)*exp(-(s*beta)**2)+m
891     self.wave_height = f1*f0*self.amplitude
892    
893    
894     #=============================================================================
895 gross 247 class TsunamiInDeepWater(Model):
896 caltinay 1989 """
897     Runs the deep water tsunami model based on a simplified form of the
898     shallow water equation.
899 gross 247
900 jfenwick 2625 *d^2 h/dt^2 =div(c grad(h))*
901 gross 247
902 caltinay 1989 where h is the wave height above sea level, and c=sqrt(g*H),
903     with H - depth of the water level, g - gravity constant
904 gross 247
905 caltinay 1989 The simulation uses the Verlet scheme.
906     """
907     def __init__(self,**kwargs):
908     Model.__init__(self,**kwargs)
909     self.declareParameter(domain = None,
910     wave_height = 1.,
911     wave_velocity = 0.,
912     initial_time_step = None,
913     bathymetry = 1.,
914     safety_factor = 1.)
915 gross 247
916 caltinay 1989 def doInitialization(self):
917     """
918     Initializes the time integration scheme
919     """
920 gross 247
921 caltinay 1989 self.__pde = LinearPDE(self.domain)
922     self.__pde.setSolverMethod(self.__pde.LUMPING)
923     self.__pde.setValue(D=1.)
924 jfenwick 2455 self.__c2 = RegionOnEarthSurface.GRAVITY*self.bathymetry/(RegionOnEarthSurface.RADIUS*2*numpy.pi/360.)**2
925 caltinay 1989 c_max = math.sqrt(Lsup(self.__c2))
926     self.__dt = self.safety_factor*inf(self.domain.getSize()/(sqrt(self.__c2)+EPS*c_max))
927     if self.initial_time_step==None:
928     self.initial_time_step=self.__dt
929     self.trace("Maximum wave velocity %s m/sec" % c_max)
930     self.trace("Time step size is %s sec" % self.__dt)
931 gross 247
932 caltinay 1989 def getSafeTimeStepSize(self, dt):
933     """
934     Returns new step size
935 gross 247
936 jfenwick 2625 :param dt: last time step size used
937     :type dt: ``float``
938     :return: time step size that can safely be used
939     :rtype: ``float``
940 caltinay 1989 """
941     return self.__dt
942 gross 247
943 caltinay 1989 def doStepPostprocessing(self,dt):
944     """
945     Performs the time step using the Verlet scheme
946 gross 247
947 jfenwick 2625 :param dt: time step size to be used
948     :type dt: ``float``
949 caltinay 1989 """
950     self.__pde.setValue(X=-self.__c2*grad(self.wave_height))
951 gross 247
952 caltinay 1989 new_height = self.wave_height+dt*self.wave_velocity+dt*(self.initial_time_step+dt)/2*self.__pde.getSolution()
953 gross 247
954 caltinay 1989 self.wave_velocity = (new_height-self.wave_height)/dt
955     self.wave_height = new_height
956     self.initial_time_step = dt
957     self.trace("Wave height range is %e %e" % (inf(self.wave_height), sup(self.wave_height)))
958 gross 259
959    
960 caltinay 1989 #=============================================================================
961 gross 247 class SurfMovie(Model):
962 caltinay 1989 """
963     movie from a wave propagation on the sea
964 gross 247
965 jfenwick 2625 :ivar time: current time
966     :ivar bathymetry: scalar data set
967     :ivar wave_height: vector data set
968     :ivar filename: name of the movie file
969 caltinay 1989 """
970     def __init__(self,**kwargs):
971     Model.__init__(self,**kwargs)
972     self.declareParameter(bathymetry = 1.,
973     wave_height = 1.,
974     coastline = None,
975     t = 0.,
976     dt = 1.,
977     south = 2.,
978     north = 5.,
979     east = 3.,
980     west = 15.,
981     max_height = 1.,
982     filename = "tsunami.mpg")
983 gross 247
984 caltinay 1989 def doInitialization(self):
985     """
986     Initializes the time integration scheme
987     """
988     self.__frame_name="tsunami"
989     self.__next_t=self.dt
990     self.__fileIndex=0
991 gross 247
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 jfenwick 2455 bathXData = numpy.zeros(numXPoints, numpy.float_)
1012     bathYData = numpy.zeros(numYPoints, numpy.float_)
1013 caltinay 1989 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 jfenwick 2455 height = numpy.zeros(numColours, numpy.float_)
1079     red = numpy.zeros(numColours, numpy.float_)
1080     green = numpy.zeros(numColours, numpy.float_)
1081     blue = numpy.zeros(numColours, numpy.float_)
1082 caltinay 1989 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 jfenwick 2625 :param dt:
1233 gross 247 """
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 jfenwick 2625 :param dt: last time step size used
1254     :type dt: ``float``
1255     :return: time step size that can savely be used
1256     :rtype: ``float``
1257 caltinay 1989 """
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 caltinay 2067 sq.t_end = 30000.
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