/[escript]/trunk/modellib/py_src/tsunami.py
ViewVC logotype

Contents of /trunk/modellib/py_src/tsunami.py

Parent Directory Parent Directory | Revision Log Revision Log

Revision 326 - (show annotations)
Wed Dec 7 03:55:19 2005 UTC (17 years, 2 months ago) by cochrane
File MIME type: text/x-python
File size: 47914 byte(s) 
Fixed bug with not getting the movie parameter correctly.  Also no longer
cleaning up the movie files if the movie generation failed.
1 #!/usr/bin/env python
2
3 # $Id$
4
5 import os, sys
6 import vtk
7 from esys.escript import *
8 from esys.escript.linearPDEs import LinearPDE
9 from esys.escript.modelframe import Model
10 import numarray
11 import math
12 #import urllib
13 import urllib2
14
15 EPS=1.e-8
16
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 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 f00=1.-f01
60 f11=(x_lat-j1*dy)/dy
61 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 def __init__(self,long=0,lat=0):
83 self.long=long
84 self.lat=lat
85 def __str__(self):
86 return "(%s,%s)"%(self.long,self.lat)
87
88 def __sub__(self,other):
89 return self.dist(other)
90
91 def split(self,p,t):
92 return PointOnEarthSurface(long=self.long+t*(p.long-self.long),lat=self.lat+t*(p.lat-self.lat))
93
94 def midPoint(self,other):
95 return PointOnEarthSurface(long=(self.long+other.long)/2,lat=(self.lat+other.lat)/2)
96
97 def dist(self,other):
98 return math.sqrt((self.long-other.long)**2+(self.lat-other.lat)**2)
99
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 def __init__(self,west_south=PointOnEarthSurface(),east_north=PointOnEarthSurface(),resolution=1.):
111 if resolution<=0:
112 raise ValueError,"resolution must be positive"
113 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 raise ValueError,"east south corner must be further down than west north corner"
117 if east_north.lat-west_south.lat < resolution/2:
118 raise ValueError,"latitude length of region must be 2*bigger then resolution"
119 if east_north.long-west_south.long < resolution/2:
120 raise ValueError,"longitude length of region must be 2*bigger then resolution"
121
122 self.west_south=west_south
123 self.east_north=east_north
124 self.resolution=resolution
125
126 def __str__(self):
127 return "RegionOnEarthSurface between %s and %s"%(str(self.west_south),str(self.east_north))
128
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 return self.west_south.long==p.long
134 if face==self.EAST:
135 return p.long==self.east_north.long
136 if face==self.SOUTH:
137 return self.west_south.lat==p.lat
138 if face==self.NORTH:
139 return p.lat==self.east_north.lat
140
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 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
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 if p0.lat<p1.lat:
166 return -1
167 elif p0.lat>p1.lat:
168 return 1
169 else:
170 return 0
171 elif f0 == self.EAST:
172 if p0.lat<p1.lat:
173 return 1
174 elif p0.lat>p1.lat:
175 return -1
176 else:
177 return 0
178 elif f0 == self.NORTH:
179 if p0.long<p1.long:
180 return -1
181 elif p0.long>p1.long:
182 return 1
183 else:
184 return 0
185 else:
186 if p0.long<p1.long:
187 return 1
188 elif p0.long>p1.long:
189 return -1
190 else:
191 return 0
192
193
194 def isInRegion(self,p):
195 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
200 def cutSegment(self,p0,p1):
201 t=None
202 p=None
203 tmp=self.interseptOfSegment(p0,p1,d=0,v=self.west_south.long)
204 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 tmp=self.interseptOfSegment(p0,p1,d=0,v=self.east_north.long)
211 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 tmp=self.interseptOfSegment(p0,p1,d=1,v=self.west_south.lat)
218 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 tmp=self.interseptOfSegment(p0,p1,d=1,v=self.east_north.lat)
225 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 a=p0.long
238 b=p1.long
239 else:
240 a=p0.lat
241 b=p1.lat
242 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 integ+=(p1.lat-p0.lat)*(p1.long+p0.long)-(p1.long-p0.long)*(p1.lat+p0.lat)
326 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 def makeTriangulation(self,west_south_is_water=True,east_south_is_water=True,west_north_is_water=True,east_north_is_water=True):
348 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 if pl[i]-short_pl[-1]>0*EPS+self.region.resolution/10:
358 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 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 mid=v_tmp[0].midPoint(v_tmp[1])
379 n_long=-d_lat/l
380 n_lat=d_long/l
381 s=self.region.resolution
382 for seg in s_tmp:
383 p0=vertices[seg[0]]
384 p1=vertices[seg[1]]
385 a_long=p1.long-p0.long
386 a_lat=p1.lat-p0.lat
387 d=a_lat*n_long-a_long*n_lat
388 if abs(d)>0.:
389 t=((mid.lat-p0.lat)*n_long-(mid.long-p0.long)*n_lat)/d
390 if 0<=t and t<=1:
391 s_tmp=((mid.lat-p0.lat)*a_long-(mid.long-p0.long)*a_lat)/d
392 if s_tmp>EPS: s=min(s,s_tmp)
393 h=PointOnEarthSurface(long=mid.long+s/2*n_long,lat=mid.lat+s/2*n_lat)
394 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 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
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 walking_on_water=True
426 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 for i in range(len(vertices)): f.writelines("%d %e %e\n"%(i,vertices[i].long,vertices[i].lat))
558 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 for i in range(len(holes)): f.writelines("%d %e %e\n"%(i,holes[i].long,holes[i].lat))
562 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 os.remove("%s.msh"%self.fn)
614 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 self.coastline_stream=urllib2.urlopen(c)
646 self.bathymetry_stream=urllib2.urlopen(b)
647
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 while abs(x_grd[nx]-ox)>EPS*diam:
696 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 reg=RegionOnEarthSurface(PointOnEarthSurface(lat=self.south,long=self.west),PointOnEarthSurface(lat=self.north,long=self.east),self.resolution)
741 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 segs.append(PointOnEarthSurface(long=float(x[0]),lat=float(x[1])))
750 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 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 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 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 self.trace("source length = %s"%(dist*2.))
788 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 m=whereNegative(s)
793 f1=(1.-m)*exp(-(s*beta)**2)+m
794 # x_long-direction
795 s=abs(x_long_hat)-dist
796 m=whereNegative(s)
797 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 safety_factor=1.)
822
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 self.__c2=RegionOnEarthSurface.GRAVITY*self.bathymetry/(RegionOnEarthSurface.RADIUS*2*numarray.pi/360.)**2
831 c_max=math.sqrt(Lsup(self.__c2))
832 self.__dt=self.safety_factor*inf(self.domain.getSize()/(sqrt(self.__c2)+EPS*c_max))
833 print c_max,inf(self.domain.getSize())
834 if self.initial_time_step==None: self.initial_time_step=self.__dt
835 self.trace("maximum wave velocity %s m/sec"%c_max)
836 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
859 new_height=self.wave_height+dt*self.wave_velocity+dt*(self.initial_time_step+dt)/2*self.__pde.getSolution()
860
861 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 max_height=1.0,
887 filename="movie.mpg")
888
889 def doInitialization(self):
890 """
891 Initializes the time integration scheme
892 """
893 self.__fn=os.tempnam()+".xml"
894 self.__frame_name=os.tempnam()
895 self.__next_t=self.dt
896 # self.coastline.getVTK()
897 # self.bathymetry.getVTK()
898 # wndow(south,west,north,east)
899
900 # set up the movie parameters
901 self.paramsFileString = "REFERENCE_FRAME DECODED\n"
902 self.paramsFileString += "FRAME_RATE 24\n"
903 self.paramsFileString += "OUTPUT %s\n" % self.filename
904 self.paramsFileString += "PATTERN IBBPBBPBB\n"
905 self.paramsFileString += "FORCE_ENCODE_LAST_FRAME\n"
906 self.paramsFileString += "GOP_SIZE 20\n"
907 self.paramsFileString += "BSEARCH_ALG CROSS2\n"
908 self.paramsFileString += "PSEARCH_ALG TWOLEVEL\n"
909 self.paramsFileString += "IQSCALE 10\n"
910 self.paramsFileString += "PQSCALE 11\n"
911 self.paramsFileString += "BQSCALE 16\n"
912 self.paramsFileString += "RANGE 8\n"
913 self.paramsFileString += "SLICES_PER_FRAME 1\n"
914 self.paramsFileString += "BASE_FILE_FORMAT PNM\n"
915 self.paramsFileString += "INPUT_DIR \n"
916 self.paramsFileString += "INPUT_CONVERT *\n"
917 self.paramsFileString += "INPUT\n"
918
919 self.firstFrame = True
920
921 self.imageFiles = []
922
923 # the bathymmetry colourmap
924 data = []
925 data.append([-8000, 0, 0, 0])
926 data.append([-7000, 0, 5, 25])
927 data.append([-6000, 0, 10, 50])
928 data.append([-5000, 0, 80, 125])
929 data.append([-4000, 0, 150, 200])
930 data.append([-3000, 86, 197, 184])
931 data.append([-2000, 172, 245, 168])
932 data.append([-1000, 211, 250, 211])
933 data.append([0, 16, 48, 16])
934 data.append([300, 70, 150, 50])
935 data.append([500, 146, 126, 60])
936 data.append([1000, 198, 178, 80])
937 data.append([1250, 250, 230, 100])
938 data.append([1500, 250, 234, 126])
939 data.append([1750, 252, 238, 152])
940 data.append([2000, 252, 243, 177])
941 data.append([2250, 253, 249, 216])
942 data.append([2500, 255, 255, 255])
943
944 # the amount to scale the data by
945 scale = 255.0
946 numColours = len(data)
947
948 # convert the colourmap into something vtk is more happy with
949 height = numarray.zeros(numColours, numarray.Float)
950 red = numarray.zeros(numColours, numarray.Float)
951 green = numarray.zeros(numColours, numarray.Float)
952 blue = numarray.zeros(numColours, numarray.Float)
953 for i in range(numColours):
954 (h, r, g, b) = data[i]
955 height[i] = float(h)
956 red[i] = float(r)/scale
957 green[i] = float(g)/scale
958 blue[i] = float(b)/scale
959
960 print "Loading bathymmetry data..."
961 # grab the z data
962 bathZData = self.bathymetry.getData()
963
964 # get the origin
965 bathOrigin = self.bathymetry.getOrigin()
966 # get the delta_x and delta_y
967 bathSpacing = self.bathymetry.getSpacing()
968
969 # now construct the x and y data from the spacing and the origin
970 numXPoints = bathZData.shape[1]
971 numYPoints = bathZData.shape[0]
972 numPoints = numXPoints*numYPoints
973
974 bathXData = numarray.zeros(numXPoints, numarray.Float)
975 bathYData = numarray.zeros(numYPoints, numarray.Float)
976 for i in range(numXPoints):
977 bathXData[i] = bathOrigin[0] + i*bathSpacing[0]
978
979 for i in range(numYPoints):
980 bathYData[i] = bathOrigin[1] + i*bathSpacing[1]
981
982 # calculate the appropriate window size
983 dataWidth = max(bathXData) - min(bathXData)
984 dataHeight = max(bathYData) - min(bathYData)
985 winHeight = 600
986 winWidth = int(dataWidth*float(winHeight)/dataHeight)
987
988 print "Done loading bathymmetry data"
989
990 ### now do the vtk stuff
991
992 # make sure rendering is offscreen
993 factGraphics = vtk.vtkGraphicsFactory()
994 factGraphics.SetUseMesaClasses(1)
995
996 factImage = vtk.vtkImagingFactory()
997 factImage.SetUseMesaClasses(1)
998
999 # make the bathymmetry colourmap
1000 transFunc = vtk.vtkColorTransferFunction()
1001 transFunc.SetColorSpaceToRGB()
1002 for i in range(numColours):
1003 transFunc.AddRGBPoint(height[i], red[i], green[i], blue[i])
1004 h = height[i]
1005 while i < numColours-1 and h < height[i+1]:
1006 h += 1
1007 transFunc.AddRGBPoint(h, red[i], green[i], blue[i])
1008
1009 # set up the lookup table for the wave data
1010 refLut = vtk.vtkLookupTable()
1011 self.lutTrans = vtk.vtkLookupTable()
1012 refLut.Build()
1013 alpha = 0.7 # alpha channel value
1014 for i in range(256):
1015 (r,g,b,a) = refLut.GetTableValue(255-i)
1016 if g == 1.0 and b < 0.5 and r < 0.5:
1017 self.lutTrans.SetTableValue(i, r, g, b, 0.0)
1018 else:
1019 self.lutTrans.SetTableValue(i, r, g-0.2, b, alpha)
1020
1021 print "Generating the bathymmetry vtk object..."
1022
1023 # set up the renderer and the render window
1024 self.ren = vtk.vtkRenderer()
1025 self.renWin = vtk.vtkRenderWindow()
1026 self.renWin.AddRenderer(self.ren)
1027 self.renWin.SetSize(winWidth, winHeight)
1028 self.renWin.OffScreenRenderingOn()
1029
1030 # make an unstructured grid
1031 bathGrid = vtk.vtkUnstructuredGrid()
1032
1033 # make the points for the vtk data
1034 bathPoints = vtk.vtkPoints()
1035 bathPoints.SetNumberOfPoints(numPoints)
1036
1037 # make the vtk bathymmetry data set
1038 bathData = vtk.vtkFloatArray()
1039 bathData.SetNumberOfComponents(1)
1040 bathData.SetNumberOfTuples(numPoints)
1041
1042 # add the points and data
1043 count = 0
1044 for i in range(numXPoints):
1045 for j in range(numYPoints):
1046 bathPoints.InsertPoint(count, bathXData[i], bathYData[j], 0.0)
1047 bathData.InsertTuple1(count, bathZData[j,i])
1048 count += 1
1049
1050 # set the data to the grid
1051 bathGrid.SetPoints(bathPoints)
1052 bathGrid.GetPointData().SetScalars(bathData)
1053
1054 # do a delaunay on the grid
1055 bathDelaunay = vtk.vtkDelaunay2D()
1056 bathDelaunay.SetInput(bathGrid)
1057 bathDelaunay.SetTolerance(0.001)
1058 bathDelaunay.SetAlpha(2.5)
1059
1060 # set up the mapper
1061 bathMapper = vtk.vtkDataSetMapper()
1062 bathMapper.SetInput(bathDelaunay.GetOutput())
1063 bathMapper.SetLookupTable(transFunc)
1064
1065 # set up the actor
1066 bathActor = vtk.vtkActor()
1067 bathActor.SetMapper(bathMapper)
1068
1069 self.ren.AddActor(bathActor)
1070
1071 ### now add the coastline
1072 print "Loading the coastline data..."
1073
1074 # make the coastline points
1075 coastPoints = vtk.vtkPoints()
1076
1077 # make the coastline grid
1078 coastGrid = vtk.vtkUnstructuredGrid()
1079
1080 # now point the points and lines into the grid
1081 totalCoastPoints = 0
1082 for polyline in self.coastline.polylines:
1083 numPoints = len(polyline)
1084 coastLine = vtk.vtkPolyLine()
1085 coastLine.GetPointIds().SetNumberOfIds(numPoints)
1086 j = 0
1087 for point in polyline:
1088 coastLine.GetPointIds().SetId(j, j+totalCoastPoints)
1089 coastPoints.InsertNextPoint(point.long, point.lat, 0.0)
1090 j += 1
1091 coastGrid.InsertNextCell(coastLine.GetCellType(),
1092 coastLine.GetPointIds())
1093 totalCoastPoints += numPoints
1094
1095 coastGrid.SetPoints(coastPoints)
1096
1097 # make the coast's mapper
1098 coastMapper = vtk.vtkDataSetMapper()
1099 coastMapper.SetInput(coastGrid)
1100
1101 # make its actor
1102 coastActor = vtk.vtkActor()
1103 coastActor.SetMapper(coastMapper)
1104 coastActor.GetProperty().SetColor(0,0,0)
1105
1106 # add the actor to the renderer
1107 self.ren.AddActor(coastActor)
1108
1109 # set up the actor for the wave
1110 self.waveActor = vtk.vtkActor()
1111
1112 # add the actor to the renderer
1113 self.ren.AddActor(self.waveActor)
1114
1115 print "Done loading the coastline data"
1116
1117 def doStepPostprocessing(self, dt):
1118 """
1119 Does any necessary postprocessing after each step
1120
1121 @param dt:
1122 """
1123 if self.t>=self.__next_t:
1124 print self.t,"write ",Lsup(self.wave_height)," to ",self.__fn
1125 saveVTK(self.__fn,h=self.wave_height)
1126 # vtkobj=...
1127 # save(self.__frame_name)
1128
1129 # make a reader for the data
1130 waveReader = vtk.vtkXMLUnstructuredGridReader()
1131 waveReader.SetFileName(self.__fn)
1132 waveReader.Update()
1133
1134 # make the grid
1135 waveGrid = waveReader.GetOutput()
1136 waveGrid.Update()
1137
1138 (zMin, zMax) = waveGrid.GetPointData().GetScalars().GetRange()
1139 print "Wave height range %f - %f" % (zMin, zMax)
1140
1141 # make a mapper for the grid
1142 waveMapper = vtk.vtkDataSetMapper()
1143 waveMapper.SetInput(waveGrid)
1144 waveMapper.SetLookupTable(self.lutTrans)
1145 waveMapper.SetScalarRange(-self.max_height*0.3, self.max_height*0.3)
1146
1147 self.waveActor.SetMapper(waveMapper)
1148
1149 # do stuff here that can only be done on the first frame
1150 if self.firstFrame:
1151 # zoom in a bit
1152 self.ren.GetActiveCamera().Zoom(2.0)
1153 self.firstFrame = False
1154
1155 # render the window
1156 self.renWin.Render()
1157
1158 # convert the render window to an image
1159 win2imgFilter = vtk.vtkWindowToImageFilter()
1160 win2imgFilter.SetInput(self.renWin)
1161
1162 # the image name
1163 imgFname = "%s%f.pnm" % (self.__frame_name, self.__next_t)
1164
1165 # save the image to file
1166 outWriter = vtk.vtkPNMWriter()
1167 outWriter.SetInput(win2imgFilter.GetOutput())
1168 outWriter.SetFileName(imgFname)
1169 outWriter.Write()
1170 print "Wrote %s" % imgFname
1171
1172 # helpful for debugging:
1173 #os.system("display %s" % imgFname)
1174
1175 self.paramsFileString += "%s\n" % imgFname
1176 self.imageFiles.append(imgFname)
1177
1178 self.__next_t+=self.dt
1179
1180 def getSafeTimeStepSize(self,dt):
1181 """
1182 returns new step size
1183
1184 @param dt: last time step size used
1185 @type dt: C{float}
1186 @return: time step size that can savely be used
1187 @rtype: C{float}
1188 """
1189 return self.__next_t-self.t
1190
1191 def doFinalization(self):
1192 """
1193 Finalises the visualisation. For instance, makes a movie of the
1194 image files.
1195 """
1196 # make the movie into self.filename
1197 self.paramsFileString += "END_INPUT\n"
1198 self.paramsFileString += "PIXEL HALF\n"
1199 self.paramsFileString += "ASPECT_RATIO 1\n"
1200
1201 # write the params string to file
1202 fp = open("%s.params" % self.filename, "w")
1203 fp.write(self.paramsFileString + '\n')
1204 fp.close()
1205
1206 print "Performing conversion to mpeg"
1207 convertString = "ppmtompeg %s.params" % self.filename
1208 result = os.system(convertString)
1209 if result != 0:
1210 print "An error occurred in mpeg conversion"
1211
1212 # now clean up the image files
1213 if result == 0:
1214 print "Removing temporary image files"
1215 os.unlink("%s.params" % self.filename)
1216 for fname in self.imageFiles:
1217 os.unlink(fname)
1218
1219 def main():
1220 from esys.escript.modelframe import Link,Simulation
1221 from esys.modellib.input import Sequencer
1222
1223 oc=OceanRegionCollector()
1224 oc.north= 26.7
1225 oc.west= 102.9
1226 oc.range360= True
1227 oc.east= 232.6
1228 oc.resolution= 1.
1229 oc.south= -71.3
1230
1231
1232 b=Bathymetry()
1233 b.source=Link(oc,"bathymetry_stream")
1234
1235 oreg=OceanRegion()
1236 oreg.source=Link(oc,"coastline_stream")
1237 oreg.resolution=Link(oc,"resolution")
1238 oreg.south=Link(oc,"south")
1239 oreg.north=Link(oc,"north")
1240 oreg.east=Link(oc,"east")
1241 oreg.west=Link(oc,"west")
1242 oreg.bathymetry_data=Link(b,"bathymetry")
1243
1244 src=TsunamiSource()
1245 src.domain=Link(oreg,"domain")
1246 src.decay_zone= 0.01
1247 src.end_long= 185.
1248 src.end_lat= -37.
1249 src.width=0.5
1250 src.start_long= 174.
1251 src.start_lat= -15.
1252 src.amplitude= 3
1253
1254 ts=TsunamiInDeepWater()
1255 ts.domain=Link(oreg,"domain")
1256 ts.wave_height=Link(src,"wave_height")
1257 ts.wave_velocity=0.
1258 ts.bathymetry=Link(oreg,"bathymetry")
1259
1260 sq=Sequencer()
1261 sq.t_end=15000.
1262
1263 sm=SurfMovie()
1264 sm.bathymetry=Link(b,"bathymetry")
1265 sm.wave_height=Link(ts,"wave_height")
1266 sm.coastline=Link(oreg,"coastline")
1267 sm.t=Link(sq,"t")
1268 sm.filename="mymovie.mpg"
1269 sm.north= 8.7
1270 sm.west= 138.9
1271 sm.dt= 50.
1272 sm.east= 196.6
1273 sm.south= -53.3
1274 sm.max_height=Link(src,"amplitude")
1275
1276 s=Simulation([sq,oc,b,oreg,src,ts,sm])
1277 s.writeXML()
1278 s.run()
1279
1280 if __name__=="__main__":
1281 from esys.modellib import tsunami
1282 tsunami.main()
1283
1284 # vim: expandtab shiftwidth=4:
  ViewVC Help
Powered by ViewVC 1.1.26