/[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 1809 - (hide annotations)
Thu Sep 25 06:43:44 2008 UTC (11 years ago) by ksteube
File MIME type: text/x-python
File size: 48777 byte(s)
Copyright updated in all python files

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

  ViewVC Help
Powered by ViewVC 1.1.26