/[escript]/branches/windows_from_1456_trunk_merge_to_1459_trunk/modellib/py_src/tsunami.py
ViewVC logotype

Contents of /branches/windows_from_1456_trunk_merge_to_1459_trunk/modellib/py_src/tsunami.py

Parent Directory Parent Directory | Revision Log Revision Log

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