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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 247 by gross, Tue Nov 29 04:57:38 2005 UTC revision 276 by cochrane, Wed Nov 30 04:31:40 2005 UTC
# Line 1  Line 1 
1  # $Id:$  #!/usr/bin/env python
2    
3    # $Id$
4    
5  import os  import os, sys
6    import vtk
7  from esys.escript import *  from esys.escript import *
8  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
9  from esys.escript.modelframe import Model  from esys.escript.modelframe import Model
# Line 48  class GridData: Line 50  class GridData:
50            i_dx=1            i_dx=1
51            i_dy=1            i_dy=1
52            for i in range(x_array.shape[0]):            for i in range(x_array.shape[0]):
53               x0=x_array[i,0]-ox               x_long=x_array[i,0]-ox
54               x1=x_array[i,1]-oy               x_lat=x_array[i,1]-oy
55               j0=min(max(int(x0/dx),0),data_array.shape[1]-1-i_dy)               j0=min(max(int(x_long/dx),0),data_array.shape[1]-1-i_dy)
56               j1=min(max(int(x1/dy),0),data_array.shape[0]-1-i_dx)               j1=min(max(int(x_lat/dy),0),data_array.shape[0]-1-i_dx)
57               f01=(x0-j0*dx)/dx               f01=(x_long-j0*dx)/dx
58               f00=1.-f01               f00=1.-f01
59               f11=(x1-j1*dy)/dy               f11=(x_lat-j1*dy)/dy
60               f10=1.-f11               f10=1.-f11
61               H00=data_array[j1,j0]               H00=data_array[j1,j0]
62               H01=data_array[j1,j0+i_dx]               H01=data_array[j1,j0+i_dx]
# Line 76  class PointOnEarthSurface: Line 78  class PointOnEarthSurface:
78     """     """
79     coordinates of a  point on the surface of the Earth     coordinates of a  point on the surface of the Earth
80     """     """
81     def __init__(self,phi=0,theta=0):     def __init__(self,long=0,lat=0):
82         self.phi=phi         self.long=long
83         self.theta=theta         self.lat=lat
84     def __str__(self):     def __str__(self):
85         return "(%s,%s)"%(self.phi,self.theta)         return "(%s,%s)"%(self.long,self.lat)
86    
87     def __sub__(self,other):     def __sub__(self,other):
88         return self.dist(other)         return self.dist(other)
89    
90     def split(self,p,t):     def split(self,p,t):
91        phi=self.phi+t*(p.phi-self.phi)        return PointOnEarthSurface(long=self.long+t*(p.long-self.long),lat=self.lat+t*(p.lat-self.lat))
       theta=self.theta+t*(p.theta-self.theta)  
       return PointOnEarthSurface(phi,theta)  
92    
93     def midPoint(self,other):     def midPoint(self,other):
94        return PointOnEarthSurface((self.phi+other.phi)/2,(self.theta+other.theta)/2)        return PointOnEarthSurface(long=(self.long+other.long)/2,lat=(self.lat+other.lat)/2)
95    
96     def dist(self,other):     def dist(self,other):
97        return math.sqrt((self.phi-other.phi)**2+(self.theta-other.theta)**2)        return math.sqrt((self.long-other.long)**2+(self.lat-other.lat)**2)
98    
99  class RegionOnEarthSurface:  class RegionOnEarthSurface:
100       """       """
# Line 106  class RegionOnEarthSurface: Line 106  class RegionOnEarthSurface:
106       NORTH=1       NORTH=1
107       EAST=2       EAST=2
108       SOUTH=3       SOUTH=3
109       def __init__(self,south_east=PointOnEarthSurface(0.,0.),north_west=PointOnEarthSurface(10.,10.),resolution=1.):       def __init__(self,west_south=PointOnEarthSurface(),east_north=PointOnEarthSurface(),resolution=1.):
110            if resolution<=0:            if resolution<=0:
111                raise ValueError,"resolution must be positive"                raise ValueError,"resolution must be positive"
112            if  south_east.phi>=north_west.phi:            if  west_south.long>=east_north.long:
113                raise ValueError,"east south corner must be further east than west north corner"                raise ValueError,"south west corner must be further east than west north corner"
114            if  south_east.theta>=north_west.theta:            if  west_south.lat>=east_north.lat:
115                raise ValueError,"east south corner must be further down than west north corner"                raise ValueError,"east south corner must be further down than west north corner"
116            if  north_west.theta-south_east.theta < resolution/2:            if  east_north.lat-west_south.lat < resolution/2:
117                raise ValueError,"latitude length of region must be 2*bigger then resolution"                raise ValueError,"latitude length of region must be 2*bigger then resolution"
118            if  north_west.phi-south_east.phi < resolution/2:            if  east_north.long-west_south.long < resolution/2:
119                raise ValueError,"longitude length of region must be 2*bigger then resolution"                raise ValueError,"longitude length of region must be 2*bigger then resolution"
120    
121            self.south_east=south_east            self.west_south=west_south
122            self.north_west=north_west            self.east_north=east_north
123            self.resolution=resolution            self.resolution=resolution
124    
125       def __str__(self):       def __str__(self):
126         return "RegionOnEarthSurface between %s and %s"%(str(self.south_east),str(self.north_west))         return "RegionOnEarthSurface between %s and %s"%(str(self.west_south),str(self.east_north))
127    
128       def isOnFace(self,p):       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)           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):       def isOnThisFace(self,p,face):
131          if face==self.WEST:          if face==self.WEST:
132             return self.south_east.phi==p.phi             return self.west_south.long==p.long
133          if face==self.EAST:          if face==self.EAST:
134           return p.phi==self.north_west.phi           return p.long==self.east_north.long
135          if face==self.SOUTH:          if face==self.SOUTH:
136           return self.south_east.theta==p.theta           return self.west_south.lat==p.lat
137          if face==self.NORTH:          if face==self.NORTH:
138             return p.theta==self.north_west.theta             return p.lat==self.east_north.lat
139    
140       def isBoxVertex(self,p):       def isBoxVertex(self,p):
141           return ( self.isOnThisFace(p,self.WEST) and  self.isOnThisFace(p,self.SOUTH) ) or \           return ( self.isOnThisFace(p,self.WEST) and  self.isOnThisFace(p,self.SOUTH) ) or \
# Line 146  class RegionOnEarthSurface: Line 146  class RegionOnEarthSurface:
146    
147       def getFace(self,p):       def getFace(self,p):
148          # order is critical          # order is critical
149          if self.south_east.phi==p.phi: return self.WEST          if self.west_south.long==p.long: return self.WEST
150          if p.theta==self.north_west.theta: return self.NORTH          if p.lat==self.east_north.lat: return self.NORTH
151          if p.phi==self.north_west.phi: return self.EAST          if p.long==self.east_north.long: return self.EAST
152          if self.south_east.theta==p.theta: return self.SOUTH          if self.west_south.lat==p.lat: return self.SOUTH
153    
154       def comparePointsOnFace(self,p0,p1):       def comparePointsOnFace(self,p0,p1):
155           f0=self.getFace(p0)           f0=self.getFace(p0)
# Line 161  class RegionOnEarthSurface: Line 161  class RegionOnEarthSurface:
161              return 1              return 1
162           else:           else:
163              if f0 == self.WEST:              if f0 == self.WEST:
164                  if p0.theta<p1.theta:                  if p0.lat<p1.lat:
165                     return -1                     return -1
166                  elif p0.theta>p1.theta:                  elif p0.lat>p1.lat:
167                     return 1                     return 1
168                  else:                  else:
169                     return 0                     return 0
170              elif f0 == self.EAST:              elif f0 == self.EAST:
171                  if p0.theta<p1.theta:                  if p0.lat<p1.lat:
172                     return 1                     return 1
173                  elif p0.theta>p1.theta:                  elif p0.lat>p1.lat:
174                     return -1                     return -1
175                  else:                  else:
176                     return 0                     return 0
177              elif f0 == self.NORTH:              elif f0 == self.NORTH:
178                  if p0.phi<p1.phi:                  if p0.long<p1.long:
179                     return -1                     return -1
180                  elif p0.phi>p1.phi:                  elif p0.long>p1.long:
181                     return 1                     return 1
182                  else:                  else:
183                     return 0                     return 0
184              else:              else:
185                  if p0.phi<p1.phi:                  if p0.long<p1.long:
186                     return 1                     return 1
187                  elif p0.phi>p1.phi:                  elif p0.long>p1.long:
188                     return -1                     return -1
189                  else:                  else:
190                     return 0                     return 0
191    
192    
193       def isInRegion(self,p):       def isInRegion(self,p):
194           return  self.south_east.phi<=p.phi \           return  self.west_south.long<=p.long \
195               and p.phi<=self.north_west.phi \               and p.long<=self.east_north.long \
196               and self.south_east.theta<=p.theta \               and self.west_south.lat<=p.lat \
197               and p.theta<=self.north_west.theta               and p.lat<=self.east_north.lat
198    
199       def cutSegment(self,p0,p1):       def cutSegment(self,p0,p1):
200              t=None              t=None
201              p=None              p=None
202              tmp=self.interseptOfSegment(p0,p1,d=0,v=self.south_east.phi)              tmp=self.interseptOfSegment(p0,p1,d=0,v=self.west_south.long)
203              if not tmp==None:              if not tmp==None:
204                  p_tmp=p0.split(p1,tmp)                  p_tmp=p0.split(p1,tmp)
205                  if self.isInRegion(p_tmp) and t<tmp:                  if self.isInRegion(p_tmp) and t<tmp:
206                     t=tmp                     t=tmp
207                     p=p_tmp                     p=p_tmp
208    
209              tmp=self.interseptOfSegment(p0,p1,d=0,v=self.north_west.phi)              tmp=self.interseptOfSegment(p0,p1,d=0,v=self.east_north.long)
210              if not tmp==None:              if not tmp==None:
211                  p_tmp=p0.split(p1,tmp)                  p_tmp=p0.split(p1,tmp)
212                  if self.isInRegion(p_tmp) and t<tmp:                  if self.isInRegion(p_tmp) and t<tmp:
213                     t=tmp                     t=tmp
214                     p=p_tmp                     p=p_tmp
215    
216              tmp=self.interseptOfSegment(p0,p1,d=1,v=self.south_east.theta)              tmp=self.interseptOfSegment(p0,p1,d=1,v=self.west_south.lat)
217              if not tmp==None:              if not tmp==None:
218                  p_tmp=p0.split(p1,tmp)                  p_tmp=p0.split(p1,tmp)
219                  if self.isInRegion(p_tmp) and t<tmp:                  if self.isInRegion(p_tmp) and t<tmp:
220                     t=tmp                     t=tmp
221                     p=p_tmp                     p=p_tmp
222    
223              tmp=self.interseptOfSegment(p0,p1,d=1,v=self.north_west.theta)              tmp=self.interseptOfSegment(p0,p1,d=1,v=self.east_north.lat)
224              if not tmp==None:              if not tmp==None:
225                  p_tmp=p0.split(p1,tmp)                  p_tmp=p0.split(p1,tmp)
226                  if self.isInRegion(p_tmp) and t<tmp:                  if self.isInRegion(p_tmp) and t<tmp:
# Line 233  class RegionOnEarthSurface: Line 233  class RegionOnEarthSurface:
233           find t isn [0,1] such that p0+t*(p0-p1) cuts x[d]=v. if it does nit exist None is returned.           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:           if d==0:
236             a=p0.phi             a=p0.long
237             b=p1.phi             b=p1.long
238           else:           else:
239             a=p0.theta             a=p0.lat
240             b=p1.theta             b=p1.lat
241           if b==a:           if b==a:
242             if a==v:             if a==v:
243                t=0                t=0
# Line 321  class Polyline: Line 321  class Polyline:
321           for i in range(len(self.list_of_coordinates)-1):           for i in range(len(self.list_of_coordinates)-1):
322               p0=self.list_of_coordinates[i]               p0=self.list_of_coordinates[i]
323               p1=self.list_of_coordinates[i+1]               p1=self.list_of_coordinates[i+1]
324               integ+=(p1.theta-p0.theta)*(p1.phi+p0.phi)-(p1.phi-p0.phi)*(p1.theta+p0.theta)               integ+=(p1.lat-p0.lat)*(p1.long+p0.long)-(p1.long-p0.long)*(p1.lat+p0.lat)
325           if integ>=0.:           if integ>=0.:
326              return 1.              return 1.
327           else:           else:
# Line 343  class Coastline: Line 343  class Coastline:
343           for pl in self.polylines:           for pl in self.polylines:
344               out+="\n"+str(pl)               out+="\n"+str(pl)
345           return out           return out
346       def makeTriangulation(self,south_east_is_water=True,south_west_is_water=True,north_west_is_water=True,north_east_is_water=True):       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()           self.clean()
348           vertices=[]           vertices=[]
349           segments=[]           segments=[]
# Line 371  class Coastline: Line 371  class Coastline:
371                         vertices+=v_tmp                         vertices+=v_tmp
372                         segments+=s_tmp                         segments+=s_tmp
373                         # find a point in the loop:                         # find a point in the loop:
374                         d_phi=v_tmp[1].phi-v_tmp[0].phi                         d_long=v_tmp[1].long-v_tmp[0].long
375                         d_theta=v_tmp[1].theta-v_tmp[0].theta                         d_lat=v_tmp[1].lat-v_tmp[0].lat
376                         l=math.sqrt(d_phi**2+d_theta**2)                         l=math.sqrt(d_long**2+d_lat**2)
377                         mid=v_tmp[0].midPoint(v_tmp[1])                         mid=v_tmp[0].midPoint(v_tmp[1])
378                         n_phi=-d_theta/l                         n_long=-d_lat/l
379                         n_theta=d_phi/l                         n_lat=d_long/l
380                         s=self.region.resolution                         s=self.region.resolution
381                         for seg in s_tmp:                         for seg in s_tmp:
382                           p0=vertices[seg[0]]                           p0=vertices[seg[0]]
383                           p1=vertices[seg[1]]                           p1=vertices[seg[1]]
384                           a_phi=p1.phi-p0.phi                           a_long=p1.long-p0.long
385                           a_theta=p1.theta-p0.theta                           a_lat=p1.lat-p0.lat
386                           d=a_theta*n_phi-a_phi*n_theta                           d=a_lat*n_long-a_long*n_lat
387                           if abs(d)>0.:                           if abs(d)>0.:
388                              t=((mid.theta-p0.theta)*n_phi-(mid.phi-p0.phi)*n_theta)/d                              t=((mid.lat-p0.lat)*n_long-(mid.long-p0.long)*n_lat)/d
389                              if 0<=t and t<=1:                              if 0<=t and t<=1:
390                                  s_tmp=((mid.theta-p0.theta)*a_phi-(mid.phi-p0.phi)*a_theta)/d                                  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)                                  if s_tmp>EPS: s=min(s,s_tmp)
392                         h=PointOnEarthSurface(mid.phi+s/2*n_phi,mid.theta+s/2*n_theta)                         h=PointOnEarthSurface(long=mid.long+s/2*n_long,lat=mid.lat+s/2*n_lat)
393                         holes.append(h)                         holes.append(h)
394                  else:                  else:
395                     if len(short_pl)>1:                     if len(short_pl)>1:
# Line 401  class Coastline: Line 401  class Coastline:
401                            vertices.append(short_pl[i])                            vertices.append(short_pl[i])
402           # put into the bounding box:           # put into the bounding box:
403           new_vertices=[]           new_vertices=[]
404           if south_east_is_water:           if west_south_is_water:
405              new_vertices.append(PointOnEarthSurface(self.region.south_east.phi,self.region.south_east.theta))              new_vertices.append(PointOnEarthSurface(long=self.region.west_south.long,lat=self.region.west_south.lat))
406           if north_east_is_water:           if east_south_is_water:
407              new_vertices.append(PointOnEarthSurface(self.region.south_east.phi,self.region.north_west.theta))              new_vertices.append(PointOnEarthSurface(long=self.region.east_north.long,lat=self.region.west_south.lat))
408           if north_west_is_water:           if west_north_is_water:
409              new_vertices.append(PointOnEarthSurface(self.region.north_west.phi,self.region.north_west.theta))              new_vertices.append(PointOnEarthSurface(long=self.region.west_south.long,lat=self.region.east_north.lat))
410           if south_west_is_water:           if east_north_is_water:
411              new_vertices.append(PointOnEarthSurface(self.region.north_west.phi,self.region.south_east.theta))              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           # add new vertices if they don't exist yet
414           for q in new_vertices:           for q in new_vertices:
# Line 421  class Coastline: Line 421  class Coastline:
421                  vertices_on_face.append(q)                  vertices_on_face.append(q)
422           vertices_on_face.sort(self.region.comparePointsOnFace)           vertices_on_face.sort(self.region.comparePointsOnFace)
423           index=0           index=0
424           walking_on_water=south_east_is_water           walking_on_water=west_south_is_water
425           l=len(vertices_on_face)           l=len(vertices_on_face)
426           while index<l:           while index<l:
427               p1=vertices_on_face[(index+1)%l]               p1=vertices_on_face[(index+1)%l]
# Line 553  class EarthTriangulation: Line 553  class EarthTriangulation:
553             poly_file=self.fn+".poly"             poly_file=self.fn+".poly"
554             f=open(poly_file,"w")             f=open(poly_file,"w")
555             f.writelines("%d %d %d %d\n"%(len(vertices),2,0,0))             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].phi,vertices[i].theta))             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))             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]))             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)))             f.writelines("%d\n"%(len(holes)))
560             for i in range(len(holes)):  f.writelines("%d %e %e\n"%(i,holes[i].phi,holes[i].theta))             for i in range(len(holes)):  f.writelines("%d %e %e\n"%(i,holes[i].long,holes[i].lat))
561             f.close()             f.close()
562             # start mesh generator:             # start mesh generator:
563             os.system(self.GENERATOR%(resolution**2,poly_file))             os.system(self.GENERATOR%(resolution**2,poly_file))
# Line 736  class OceanRegion(Model): Line 736  class OceanRegion(Model):
736                line=line.strip()                line=line.strip()
737                if line[:7]=="#region":                if line[:7]=="#region":
738                   data=line[line.index("[")+1:line.index("]")].split(",")                   data=line[line.index("[")+1:line.index("]")].split(",")
739                   reg=RegionOnEarthSurface(PointOnEarthSurface(self.south,self.west),PointOnEarthSurface(self.north,self.east),self.resolution)                   reg=RegionOnEarthSurface(PointOnEarthSurface(lat=self.south,long=self.west),PointOnEarthSurface(lat=self.north,long=self.east),self.resolution)
740                   self.trace(str(reg))                   self.trace(str(reg))
741                   self.coastline=Coastline(region=reg,name=data_name)                   self.coastline=Coastline(region=reg,name=data_name)
742                elif line.find("Shore Bin")>-1:                elif line.find("Shore Bin")>-1:
# Line 745  class OceanRegion(Model): Line 745  class OceanRegion(Model):
745                     segs=[]                     segs=[]
746                if not (line=="" or line[0]=="#" or line[0]==">") :                if not (line=="" or line[0]=="#" or line[0]==">") :
747                    x=line.split()                    x=line.split()
748                    segs.append(PointOnEarthSurface(float(x[0]),float(x[1])))                    segs.append(PointOnEarthSurface(long=float(x[0]),lat=float(x[1])))
749                line=f.readline()                line=f.readline()
750             self.coastline.append(Polyline(segs,name))             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]])             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(south_east_is_water=d[0]<0,south_west_is_water=d[1]<0,north_east_is_water=d[2]<0,north_west_is_water=d[3]<0).getFinleyDomain()             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.)             self.bathymetry=maximum(-self.bathymetry_data.interpolate(Function(self.domain).getX()),0.)
754    
755    
# Line 776  class TsunamiSource(Model): Line 776  class TsunamiSource(Model):
776             """             """
777             beta=math.sqrt(-math.log(self.GAMMA))/self.decay_zone             beta=math.sqrt(-math.log(self.GAMMA))/self.decay_zone
778             x=self.domain.getX()             x=self.domain.getX()
779             x0=x[0]             x_long=x[0]
780             x1=x[1]             x_lat=x[1]
781             mid0=(self.start_long+self.end_long)/2             mid_long=(self.start_long+self.end_long)/2
782             mid1=(self.start_lat+self.end_lat)/2             mid_lat=(self.start_lat+self.end_lat)/2
783             dist=math.sqrt((mid0-self.end_long)**2+(mid1-self.end_lat)**2)             dist=math.sqrt((mid_long-self.end_long)**2+(mid_lat-self.end_lat)**2)
784             a=(self.start_long-mid0)/dist             a=(self.start_long-mid_long)/dist
785             b=(self.end_long-mid1)/dist             b=(self.start_lat-mid_lat)/dist
786             self.trace("source length = %s"%(dist*2.))             self.trace("source length = %s"%(dist*2.))
787               x_long_hat= a*(x_long-mid_long)+b*(x_lat-mid_lat)
788             x0_hat= a*(x0-mid0)+b*(x1-mid1)             x_lat_hat=-b*(x_long-mid_long)+a*(x_lat-mid_lat)
789             x1_hat=-b*(x0-mid0)+a*(x1-mid1)             # x_lat-direction
790             # x1-direction             s=abs(x_lat_hat)-self.width
            s=abs(x1_hat)-self.width  
791             m=s.whereNegative()             m=s.whereNegative()
792             f1=(1.-m)*exp(-(s*beta)**2)+m             f1=(1.-m)*exp(-(s*beta)**2)+m
793             # x0-direction             # x_long-direction
794             s=abs(x0_hat)-dist             s=abs(x_long_hat)-dist
795             m=s.whereNegative()             m=s.whereNegative()
796             f0=(1.-m)*exp(-(s*beta)**2)+m             f0=(1.-m)*exp(-(s*beta)**2)+m
797             self.wave_height=f1*f0*self.amplitude             self.wave_height=f1*f0*self.amplitude
# Line 818  class TsunamiInDeepWater(Model): Line 817  class TsunamiInDeepWater(Model):
817                                   wave_velocity=0.,                                   wave_velocity=0.,
818                                   initial_time_step=None,                                   initial_time_step=None,
819                                   bathymetry=1.,                                   bathymetry=1.,
820                                   safety_factor=0.1)                                   safety_factor=1.)
821    
822         def doInitialization(self):         def doInitialization(self):
823             """             """
# Line 831  class TsunamiInDeepWater(Model): Line 830  class TsunamiInDeepWater(Model):
830             c_max=math.sqrt(Lsup(self.__c2))             c_max=math.sqrt(Lsup(self.__c2))
831             self.__dt=self.safety_factor*inf(self.domain.getSize()/(sqrt(self.__c2)+EPS*c_max))             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             if self.initial_time_step==None: self.initial_time_step=self.__dt
833             self.trace("maximum wave velocity %s m/sec^2"%c_max)             self.trace("maximum wave velocity %s m/sec"%c_max)
834             self.trace("Time step size is %s sec"%self.__dt)             self.trace("Time step size is %s sec"%self.__dt)
835    
836    
# Line 854  class TsunamiInDeepWater(Model): Line 853  class TsunamiInDeepWater(Model):
853             @type dt: C{float}             @type dt: C{float}
854             """             """
855             self.__pde.setValue(X=-self.__c2*grad(self.wave_height))             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()             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             self.wave_velocity=(new_height-self.wave_height)/dt
860             self.wave_height=new_height             self.wave_height=new_height
861             self.initial_time_step=dt             self.initial_time_step=dt
# Line 882  class SurfMovie(Model): Line 883  class SurfMovie(Model):
883                                   west=15.,                                   west=15.,
884                                   filename="movie.mpg")                                   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):         def doInitialization(self):
909            """            """
910            Initializes the time integartion scheme            Initializes the time integration scheme
911            """            """
912            self.__fn=os.tempnam()+".xml"            self.__fn=os.tempnam()+".xml"
913            self.__frame_name=os.tempnam()            self.__frame_name=os.tempnam()
# Line 893  class SurfMovie(Model): Line 916  class SurfMovie(Model):
916            # self.bathymetry.getVTK()            # self.bathymetry.getVTK()
917            # wndow(south,west,north,east)            # 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):         def doStepPostprocessing(self, dt):
1114          """          """
1115          Does any necessary postprocessing after each step          Does any necessary postprocessing after each step
# Line 904  class SurfMovie(Model): Line 1121  class SurfMovie(Model):
1121               saveVTK(self.__fn,h=self.wave_height)               saveVTK(self.__fn,h=self.wave_height)
1122               # vtkobj=...               # vtkobj=...
1123               # save(self.__frame_name)               # save(self.__frame_name)
1124    
1125             # make a reader for the data
1126             waveReader = vtk.vtkXMLUnstructuredGridReader()
1127             waveReader.SetFileName(self.__fn)
1128             waveReader.Update()
1129    
1130                 # make the grid
1131                 waveGrid = waveReader.GetOutput()
1132             waveGrid.Update()
1133    
1134                 (zMin, zMax) = waveGrid.GetPointData().GetScalars().GetRange()
1135                 print "Wave height range %f - %f" % (zMin, zMax)
1136    
1137                 # make a mapper for the grid
1138                 waveMapper = vtk.vtkDataSetMapper()
1139                 waveMapper.SetInput(waveGrid)
1140             waveMapper.SetLookupTable(self.lutTrans)
1141                 waveMapper.SetScalarRange(zMin, zMax)
1142    
1143                 self.waveActor.SetMapper(waveMapper)
1144    
1145                 # do stuff here that can only be done on the first frame
1146                 if self.firstFrame:
1147                     # zoom in a bit
1148                     self.ren.GetActiveCamera().Zoom(2.0)
1149                     self.firstFrame = False
1150    
1151                 # render the window
1152                 self.renWin.Render()
1153    
1154                 # convert the render window to an image
1155                 win2imgFilter = vtk.vtkWindowToImageFilter()
1156                 win2imgFilter.SetInput(self.renWin)
1157    
1158                 # the image name
1159                 imgFname = "%s%f.pnm" % (self.__frame_name, self.__next_t)
1160      
1161                 # save the image to file
1162                 outWriter = vtk.vtkPNMWriter()
1163                 outWriter.SetInput(win2imgFilter.GetOutput())
1164                 outWriter.SetFileName(imgFname)
1165                 outWriter.Write()
1166                 print "Wrote %s" % imgFname
1167    
1168             # helpful for debugging:
1169             #os.system("display %s" % imgFname)
1170    
1171                 self.paramsFileString += "%s\n" % imgFname
1172             self.imageFiles.append(imgFname)
1173    
1174               self.__next_t+=self.dt               self.__next_t+=self.dt
1175    
1176         def getSafeTimeStepSize(self,dt):         def getSafeTimeStepSize(self,dt):
# Line 923  class SurfMovie(Model): Line 1190  class SurfMovie(Model):
1190            image files.            image files.
1191            """            """
1192            # make the movie into self.filename            # make the movie into self.filename
1193            pass            self.paramsFileString += "END_INPUT\n"
1194              self.paramsFileString += "PIXEL HALF\n"
1195              self.paramsFileString += "ASPECT_RATIO 1\n"
1196    
1197              # write the params string to file
1198              fp = open("%s.params" % self.filename, "w")
1199              fp.write(self.paramsFileString + '\n')
1200              fp.close()
1201    
1202              print "Performing conversion to mpeg"
1203              convertString = "ppmtompeg %s.params" % self.filename
1204              result = os.system(convertString)
1205              if result != 0:
1206                  print "An error occurred in mpeg conversion"
1207    
1208          # now clean up the image files
1209          print "Removing temporary image files"
1210          os.unlink("%s.params" % self.filename)
1211          for fname in self.imageFiles:
1212              os.unlink(fname)
1213    
1214  if __name__=="__main__":  if __name__=="__main__":
1215     from esys.escript.modelframe import Link,Simulation     from esys.escript.modelframe import Link,Simulation
1216     from esys.modellib.input import Sequencer     from esys.modellib.input import Sequencer
1217    
1218     oc=OceanRegionCollector()     oc=OceanRegionCollector()
1219     oc.resolution=0.5     oc.resolution=2
1220     oc.south=-45.5     oc.south=-45.5
1221     oc.north=-5.4     oc.north=-5.4
1222     oc.east=161.1     oc.east=161.1
# Line 954  if __name__=="__main__": Line 1239  if __name__=="__main__":
1239     src=TsunamiSource()     src=TsunamiSource()
1240     src.domain=Link(oreg,"domain")     src.domain=Link(oreg,"domain")
1241     src.start_lat=-10.     src.start_lat=-10.
    src.start_long=110.  
1242     src.end_lat=-12.     src.end_lat=-12.
1243       src.start_long=110.
1244     src.end_long=120.     src.end_long=120.
1245     src.width=0.1     src.width=1.
1246     src.decay_zone=0.01     src.decay_zone=0.01
1247     src.amplitude=1.     src.amplitude=1.
1248    
# Line 968  if __name__=="__main__": Line 1253  if __name__=="__main__":
1253     ts.bathymetry=Link(oreg,"bathymetry")     ts.bathymetry=Link(oreg,"bathymetry")
1254    
1255     sq=Sequencer()     sq=Sequencer()
1256     sq.t_end=1000.     sq.t_end=100000.
1257    
1258     sm=SurfMovie()     sm=SurfMovie()
1259     sm.bathymetry=Link(b,"bathymetry")     sm.bathymetry=Link(b,"bathymetry")
1260     sm.wave_height=Link(ts,"wave_height")     sm.wave_height=Link(ts,"wave_height")
1261     sm.coastline=Link(oreg,"coastline")     sm.coastline=Link(oreg,"coastline")
1262     sm.t=Link(sq,"t")     sm.t=Link(sq,"t")
1263     sm.dt=100.     sm.dt=5000.
1264     sm.filename="movie.mpg"     sm.filename="movie.mpg"
1265        
1266     s=Simulation([sq,oc,b,oreg,src,ts,sm])     s=Simulation([sq,oc,b,oreg,src,ts,sm])
1267     s.writeXML()     # s.writeXML()
1268     s.run()     s.run()
1269    
1270    # vim: expandtab shiftwidth=4:

Legend:
Removed from v.247  
changed lines
  Added in v.276

  ViewVC Help
Powered by ViewVC 1.1.26