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