/[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 918 - (hide annotations)
Wed Jan 3 06:30:00 2007 UTC (12 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 48332 byte(s)
fixes on ESySXML to get runmodel going.

    * object ids are now local for read and write of XML
    * ParameterSets are now written with class name
    * ParameterSets linked by other ParameterSets (previously only Models) are written now 
    * list are now lists of str (rather than bools), lists with bool, int or float are mapped to numarray
    * numarray are writen with generic types Bool, Int, Float (portability!)


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

  ViewVC Help
Powered by ViewVC 1.1.26