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

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