/[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 264 - (show annotations)
Tue Nov 29 10:06:39 2005 UTC (13 years, 4 months ago) by cochrane
File MIME type: text/x-python
File size: 47890 byte(s)
Removing image files after generation of movie.

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

  ViewVC Help
Powered by ViewVC 1.1.26