/[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 263 - (hide annotations)
Tue Nov 29 09:55:16 2005 UTC (13 years, 10 months ago) by cochrane
File MIME type: text/x-python
File size: 47697 byte(s)
Added visualisation of bathymmetry data, the coastline and the wave height
data as per earlier tsunami movies.  Also autogenerates the .params file for
ppmtompeg and generates the movie automatically.

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

  ViewVC Help
Powered by ViewVC 1.1.26