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

  ViewVC Help
Powered by ViewVC 1.1.26