/[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 1809 - (show annotations)
Thu Sep 25 06:43:44 2008 UTC (11 years ago) by ksteube
File MIME type: text/x-python
File size: 48777 byte(s)
Copyright updated in all python files

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

  ViewVC Help
Powered by ViewVC 1.1.26