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

  ViewVC Help
Powered by ViewVC 1.1.26