/[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 264 - (hide annotations)
Tue Nov 29 10:06:39 2005 UTC (13 years, 10 months ago) by cochrane
File MIME type: text/x-python
File size: 47890 byte(s)
Removing image files after generation of movie.

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 cochrane 264 self.imageFiles = []
907    
908 gross 247 def doInitialization(self):
909     """
910 cochrane 263 Initializes the time integration scheme
911 gross 247 """
912     self.__fn=os.tempnam()+".xml"
913     self.__frame_name=os.tempnam()
914     self.__next_t=self.dt
915     # self.coastline.getVTK()
916     # self.bathymetry.getVTK()
917     # wndow(south,west,north,east)
918    
919 cochrane 263 # the bathymmetry colourmap
920     data = []
921     data.append([-8000, 0, 0, 0])
922     data.append([-7000, 0, 5, 25])
923     data.append([-6000, 0, 10, 50])
924     data.append([-5000, 0, 80, 125])
925     data.append([-4000, 0, 150, 200])
926     data.append([-3000, 86, 197, 184])
927     data.append([-2000, 172, 245, 168])
928     data.append([-1000, 211, 250, 211])
929     data.append([0, 16, 48, 16])
930     data.append([300, 70, 150, 50])
931     data.append([500, 146, 126, 60])
932     data.append([1000, 198, 178, 80])
933     data.append([1250, 250, 230, 100])
934     data.append([1500, 250, 234, 126])
935     data.append([1750, 252, 238, 152])
936     data.append([2000, 252, 243, 177])
937     data.append([2250, 253, 249, 216])
938     data.append([2500, 255, 255, 255])
939    
940     # the amount to scale the data by
941     scale = 255.0
942     numColours = len(data)
943    
944     # convert the colourmap into something vtk is more happy with
945     height = numarray.zeros(numColours, numarray.Float)
946     red = numarray.zeros(numColours, numarray.Float)
947     green = numarray.zeros(numColours, numarray.Float)
948     blue = numarray.zeros(numColours, numarray.Float)
949     for i in range(numColours):
950     (h, r, g, b) = data[i]
951     height[i] = float(h)
952     red[i] = float(r)/scale
953     green[i] = float(g)/scale
954     blue[i] = float(b)/scale
955    
956     print "Loading bathymmetry data..."
957     # grab the z data
958     bathZData = self.bathymetry.getData()
959    
960     # get the origin
961     bathOrigin = self.bathymetry.getOrigin()
962     # get the delta_x and delta_y
963     bathSpacing = self.bathymetry.getSpacing()
964    
965     # now construct the x and y data from the spacing and the origin
966     numXPoints = bathZData.shape[1]
967     numYPoints = bathZData.shape[0]
968     numPoints = numXPoints*numYPoints
969    
970     bathXData = numarray.zeros(numXPoints, numarray.Float)
971     bathYData = numarray.zeros(numYPoints, numarray.Float)
972     for i in range(numXPoints):
973     bathXData[i] = bathOrigin[0] + i*bathSpacing[0]
974    
975     for i in range(numYPoints):
976     bathYData[i] = bathOrigin[1] + i*bathSpacing[1]
977    
978     # calculate the appropriate window size
979     dataWidth = max(bathXData) - min(bathXData)
980     dataHeight = max(bathYData) - min(bathYData)
981     winHeight = 600
982     winWidth = int(dataWidth*float(winHeight)/dataHeight)
983    
984     print "Done loading bathymmetry data"
985    
986     ### now do the vtk stuff
987    
988     # make sure rendering is offscreen
989     factGraphics = vtk.vtkGraphicsFactory()
990     factGraphics.SetUseMesaClasses(1)
991    
992     factImage = vtk.vtkImagingFactory()
993     factImage.SetUseMesaClasses(1)
994    
995     # make the bathymmetry colourmap
996     transFunc = vtk.vtkColorTransferFunction()
997     transFunc.SetColorSpaceToRGB()
998     for i in range(numColours):
999     transFunc.AddRGBPoint(height[i], red[i], green[i], blue[i])
1000     h = height[i]
1001     while i < numColours-1 and h < height[i+1]:
1002     h += 1
1003     transFunc.AddRGBPoint(h, red[i], green[i], blue[i])
1004    
1005     # set up the lookup table for the wave data
1006     refLut = vtk.vtkLookupTable()
1007     self.lutTrans = vtk.vtkLookupTable()
1008     refLut.Build()
1009     alpha = 0.4 # alpha channel value
1010     for i in range(256):
1011     (r,g,b,a) = refLut.GetTableValue(255-i)
1012     if g == 1.0 and b < 0.5 and r < 0.5:
1013     self.lutTrans.SetTableValue(i, r, g, b, 0.0)
1014     else:
1015     self.lutTrans.SetTableValue(i, r, g-0.2, b, alpha)
1016    
1017     print "Generating the bathymmetry vtk object..."
1018    
1019     # set up the renderer and the render window
1020     self.ren = vtk.vtkRenderer()
1021     self.renWin = vtk.vtkRenderWindow()
1022     self.renWin.AddRenderer(self.ren)
1023     self.renWin.SetSize(winWidth, winHeight)
1024     self.renWin.OffScreenRenderingOn()
1025    
1026     # make an unstructured grid
1027     bathGrid = vtk.vtkUnstructuredGrid()
1028    
1029     # make the points for the vtk data
1030     bathPoints = vtk.vtkPoints()
1031     bathPoints.SetNumberOfPoints(numPoints)
1032    
1033     # make the vtk bathymmetry data set
1034     bathData = vtk.vtkFloatArray()
1035     bathData.SetNumberOfComponents(1)
1036     bathData.SetNumberOfTuples(numPoints)
1037    
1038     # add the points and data
1039     count = 0
1040     for i in range(numXPoints):
1041     for j in range(numYPoints):
1042     bathPoints.InsertPoint(count, bathXData[i], bathYData[j], 0.0)
1043     bathData.InsertTuple1(count, bathZData[j,i])
1044     count += 1
1045    
1046     # set the data to the grid
1047     bathGrid.SetPoints(bathPoints)
1048     bathGrid.GetPointData().SetScalars(bathData)
1049    
1050     # do a delaunay on the grid
1051     bathDelaunay = vtk.vtkDelaunay2D()
1052     bathDelaunay.SetInput(bathGrid)
1053     bathDelaunay.SetTolerance(0.001)
1054     bathDelaunay.SetAlpha(2.5)
1055    
1056     # set up the mapper
1057     bathMapper = vtk.vtkDataSetMapper()
1058     bathMapper.SetInput(bathDelaunay.GetOutput())
1059     bathMapper.SetLookupTable(transFunc)
1060    
1061     # set up the actor
1062     bathActor = vtk.vtkActor()
1063     bathActor.SetMapper(bathMapper)
1064    
1065     self.ren.AddActor(bathActor)
1066    
1067     ### now add the coastline
1068     print "Loading the coastline data..."
1069    
1070     # make the coastline points
1071     coastPoints = vtk.vtkPoints()
1072    
1073     # make the coastline grid
1074     coastGrid = vtk.vtkUnstructuredGrid()
1075    
1076     # now point the points and lines into the grid
1077     totalCoastPoints = 0
1078     for polyline in self.coastline.polylines:
1079     numPoints = len(polyline)
1080     coastLine = vtk.vtkPolyLine()
1081     coastLine.GetPointIds().SetNumberOfIds(numPoints)
1082     j = 0
1083     for point in polyline:
1084     coastLine.GetPointIds().SetId(j, j+totalCoastPoints)
1085     coastPoints.InsertNextPoint(point.long, point.lat, 0.0)
1086     j += 1
1087     coastGrid.InsertNextCell(coastLine.GetCellType(),
1088     coastLine.GetPointIds())
1089     totalCoastPoints += numPoints
1090    
1091     coastGrid.SetPoints(coastPoints)
1092    
1093     # make the coast's mapper
1094     coastMapper = vtk.vtkDataSetMapper()
1095     coastMapper.SetInput(coastGrid)
1096    
1097     # make its actor
1098     coastActor = vtk.vtkActor()
1099     coastActor.SetMapper(coastMapper)
1100     coastActor.GetProperty().SetColor(0,0,0)
1101    
1102     # add the actor to the renderer
1103     self.ren.AddActor(coastActor)
1104    
1105     # set up the actor for the wave
1106     self.waveActor = vtk.vtkActor()
1107    
1108     # add the actor to the renderer
1109     self.ren.AddActor(self.waveActor)
1110    
1111     print "Done loading the coastline data"
1112    
1113 gross 247 def doStepPostprocessing(self, dt):
1114     """
1115     Does any necessary postprocessing after each step
1116    
1117     @param dt:
1118     """
1119     if self.t>=self.__next_t:
1120     print self.t,"write ",Lsup(self.wave_height)
1121     saveVTK(self.__fn,h=self.wave_height)
1122     # vtkobj=...
1123     # save(self.__frame_name)
1124 cochrane 263
1125     # get the wave data
1126     waveDomain = self.wave_height.getDomain().getX()
1127     waveX = waveDomain[0].convertToNumArray()
1128     waveY = waveDomain[1].convertToNumArray()
1129     waveZ = self.wave_height.convertToNumArray()
1130    
1131     numPoints = len(waveZ)
1132    
1133     # make the points
1134     wavePoints = vtk.vtkPoints()
1135     wavePoints.SetNumberOfPoints(numPoints)
1136    
1137     # make the vtk data array
1138     waveData = vtk.vtkFloatArray()
1139     waveData.SetNumberOfComponents(1)
1140     waveData.SetNumberOfTuples(numPoints)
1141     waveData.SetName("data")
1142    
1143     # put the data into the points and array
1144     for i in range(numPoints):
1145     wavePoints.InsertPoint(i, waveX[i], waveY[i], 0.0)
1146     waveData.InsertTuple1(i, waveZ[i])
1147    
1148     # make the grid
1149     waveGrid = vtk.vtkUnstructuredGrid()
1150     waveGrid.SetPoints(wavePoints)
1151     waveGrid.GetPointData().AddArray(waveData)
1152     waveGrid.GetPointData().SetActiveScalars("data")
1153    
1154     (zMin, zMax) = waveGrid.GetPointData().GetScalars().GetRange()
1155     print "Wave height range %f - %f" % (zMin, zMax)
1156    
1157     # do a delaunay on the data
1158     waveDelaunay = vtk.vtkDelaunay2D()
1159     waveDelaunay.SetInput(waveGrid)
1160     waveDelaunay.SetTolerance(0.001)
1161     waveDelaunay.SetAlpha(2.5)
1162    
1163     # make a mapper for the grid
1164     waveMapper = vtk.vtkDataSetMapper()
1165     waveMapper.SetInput(waveDelaunay.GetOutput())
1166     waveMapper.SetLookupTable(self.lutTrans)
1167     waveMapper.SetScalarRange(zMin, zMax)
1168    
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     self.paramsFileString += "%s\n" % imgFname
1194 cochrane 264 self.imageFiles.append(imgFname)
1195 cochrane 263
1196 gross 247 self.__next_t+=self.dt
1197    
1198     def getSafeTimeStepSize(self,dt):
1199     """
1200     returns new step size
1201    
1202     @param dt: last time step size used
1203     @type dt: C{float}
1204     @return: time step size that can savely be used
1205     @rtype: C{float}
1206     """
1207     return self.__next_t-self.t
1208    
1209     def doFinalization(self):
1210     """
1211     Finalises the visualisation. For instance, makes a movie of the
1212     image files.
1213     """
1214     # make the movie into self.filename
1215 cochrane 263 self.paramsFileString += "END_INPUT\n"
1216     self.paramsFileString += "PIXEL HALF\n"
1217     self.paramsFileString += "ASPECT_RATIO 1\n"
1218 gross 247
1219 cochrane 263 # write the params string to file
1220     fp = open("%s.params" % self.filename, "w")
1221     fp.write(self.paramsFileString + '\n')
1222     fp.close()
1223 gross 247
1224 cochrane 263 print "Performing conversion to mpeg"
1225     convertString = "ppmtompeg %s.params" % self.filename
1226     result = os.system(convertString)
1227     if result != 0:
1228     print "An error occurred in mpeg conversion"
1229    
1230 cochrane 264 # now clean up the image files
1231     print "Removing temporary image files"
1232     os.unlink("%s.params" % self.filename)
1233     for fname in self.imageFiles:
1234     os.unlink(fname)
1235    
1236 gross 247 if __name__=="__main__":
1237     from esys.escript.modelframe import Link,Simulation
1238     from esys.modellib.input import Sequencer
1239    
1240     oc=OceanRegionCollector()
1241 gross 251 oc.resolution=2
1242 gross 247 oc.south=-45.5
1243     oc.north=-5.4
1244     oc.east=161.1
1245     oc.west=108.2
1246     oc.range360=True
1247    
1248    
1249     b=Bathymetry()
1250     b.source=Link(oc,"bathymetry_stream")
1251    
1252     oreg=OceanRegion()
1253     oreg.source=Link(oc,"coastline_stream")
1254     oreg.resolution=Link(oc,"resolution")
1255     oreg.south=Link(oc,"south")
1256     oreg.north=Link(oc,"north")
1257     oreg.east=Link(oc,"east")
1258     oreg.west=Link(oc,"west")
1259     oreg.bathymetry_data=Link(b,"bathymetry")
1260    
1261     src=TsunamiSource()
1262     src.domain=Link(oreg,"domain")
1263     src.start_lat=-10.
1264 gross 251 src.end_lat=-12.
1265 gross 247 src.start_long=110.
1266     src.end_long=120.
1267     src.width=0.1
1268     src.decay_zone=0.01
1269     src.amplitude=1.
1270    
1271     ts=TsunamiInDeepWater()
1272     ts.domain=Link(oreg,"domain")
1273     ts.wave_height=Link(src,"wave_height")
1274     ts.wave_velocity=0.
1275     ts.bathymetry=Link(oreg,"bathymetry")
1276    
1277     sq=Sequencer()
1278 gross 251 sq.t_end=100000.
1279 gross 247
1280     sm=SurfMovie()
1281     sm.bathymetry=Link(b,"bathymetry")
1282     sm.wave_height=Link(ts,"wave_height")
1283     sm.coastline=Link(oreg,"coastline")
1284     sm.t=Link(sq,"t")
1285 gross 251 sm.dt=5000.
1286 gross 247 sm.filename="movie.mpg"
1287    
1288     s=Simulation([sq,oc,b,oreg,src,ts,sm])
1289 gross 259 # s.writeXML()
1290 gross 247 s.run()
1291 cochrane 263
1292     # vim: expandtab shiftwidth=4:

  ViewVC Help
Powered by ViewVC 1.1.26