/[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 251 by gross, Tue Nov 29 05:54:06 2005 UTC revision 280 by cochrane, Wed Nov 30 06:23:59 2005 UTC
# Line 1  Line 1 
1  # $Id$  #!/usr/bin/env python
2    
3    # $Id$
4    
5  import os  import os, sys
6    import vtk
7  from esys.escript import *  from esys.escript import *
8  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
9  from esys.escript.modelframe import Model  from esys.escript.modelframe import Model
# Line 815  class TsunamiInDeepWater(Model): Line 817  class TsunamiInDeepWater(Model):
817                                   wave_velocity=0.,                                   wave_velocity=0.,
818                                   initial_time_step=None,                                   initial_time_step=None,
819                                   bathymetry=1.,                                   bathymetry=1.,
820                                   safety_factor=0.1)                                   safety_factor=1.)
821    
822         def doInitialization(self):         def doInitialization(self):
823             """             """
# Line 828  class TsunamiInDeepWater(Model): Line 830  class TsunamiInDeepWater(Model):
830             c_max=math.sqrt(Lsup(self.__c2))             c_max=math.sqrt(Lsup(self.__c2))
831             self.__dt=self.safety_factor*inf(self.domain.getSize()/(sqrt(self.__c2)+EPS*c_max))             self.__dt=self.safety_factor*inf(self.domain.getSize()/(sqrt(self.__c2)+EPS*c_max))
832             if self.initial_time_step==None: self.initial_time_step=self.__dt             if self.initial_time_step==None: self.initial_time_step=self.__dt
833             self.trace("maximum wave velocity %s m/sec^2"%c_max)             self.trace("maximum wave velocity %s m/sec"%c_max)
834             self.trace("Time step size is %s sec"%self.__dt)             self.trace("Time step size is %s sec"%self.__dt)
835    
836    
# Line 851  class TsunamiInDeepWater(Model): Line 853  class TsunamiInDeepWater(Model):
853             @type dt: C{float}             @type dt: C{float}
854             """             """
855             self.__pde.setValue(X=-self.__c2*grad(self.wave_height))             self.__pde.setValue(X=-self.__c2*grad(self.wave_height))
856    
857             new_height=self.wave_height+dt*self.wave_velocity+dt*(self.initial_time_step+dt)/2*self.__pde.getSolution()             new_height=self.wave_height+dt*self.wave_velocity+dt*(self.initial_time_step+dt)/2*self.__pde.getSolution()
858    
859             self.wave_velocity=(new_height-self.wave_height)/dt             self.wave_velocity=(new_height-self.wave_height)/dt
860             self.wave_height=new_height             self.wave_height=new_height
861             self.initial_time_step=dt             self.initial_time_step=dt
# Line 877  class SurfMovie(Model): Line 881  class SurfMovie(Model):
881                                   north=5.,                                   north=5.,
882                                   east=3.,                                   east=3.,
883                                   west=15.,                                   west=15.,
884                                     max_height=1.0,
885                                   filename="movie.mpg")                                   filename="movie.mpg")
886    
887               self.paramsFileString = "REFERENCE_FRAME DECODED\n"
888               self.paramsFileString += "FRAME_RATE 24\n"
889               self.paramsFileString += "OUTPUT %s\n" % self.filename
890               self.paramsFileString += "PATTERN IBBPBBPBB\n"
891               self.paramsFileString += "FORCE_ENCODE_LAST_FRAME\n"
892               self.paramsFileString += "GOP_SIZE 20\n"
893               self.paramsFileString += "BSEARCH_ALG CROSS2\n"
894               self.paramsFileString += "PSEARCH_ALG TWOLEVEL\n"
895               self.paramsFileString += "IQSCALE 10\n"
896               self.paramsFileString += "PQSCALE 11\n"
897               self.paramsFileString += "BQSCALE 16\n"
898               self.paramsFileString += "RANGE 8\n"
899               self.paramsFileString += "SLICES_PER_FRAME 1\n"
900               self.paramsFileString += "BASE_FILE_FORMAT PNM\n"
901               self.paramsFileString += "INPUT_DIR \n"
902               self.paramsFileString += "INPUT_CONVERT *\n"
903               self.paramsFileString += "INPUT\n"
904    
905               self.firstFrame = True
906    
907               self.imageFiles = []
908    
909         def doInitialization(self):         def doInitialization(self):
910            """            """
911            Initializes the time integartion scheme            Initializes the time integration scheme
912            """            """
913            self.__fn=os.tempnam()+".xml"            self.__fn=os.tempnam()+".xml"
914            self.__frame_name=os.tempnam()            self.__frame_name=os.tempnam()
# Line 890  class SurfMovie(Model): Line 917  class SurfMovie(Model):
917            # self.bathymetry.getVTK()            # self.bathymetry.getVTK()
918            # wndow(south,west,north,east)            # wndow(south,west,north,east)
919    
920              # the bathymmetry colourmap
921              data = []
922              data.append([-8000, 0,   0,   0])
923              data.append([-7000, 0,   5,   25])
924              data.append([-6000, 0,   10,  50])
925              data.append([-5000, 0,   80,  125])
926              data.append([-4000, 0,   150, 200])
927              data.append([-3000, 86,  197, 184])
928              data.append([-2000, 172, 245, 168])
929              data.append([-1000, 211, 250, 211])
930              data.append([0,     16,  48,  16])
931              data.append([300,   70,  150, 50])
932              data.append([500,   146, 126, 60])
933              data.append([1000,  198, 178, 80])
934              data.append([1250,  250, 230, 100])
935              data.append([1500,  250, 234, 126])
936              data.append([1750,  252, 238, 152])
937              data.append([2000,  252, 243, 177])
938              data.append([2250,  253, 249, 216])
939              data.append([2500,  255, 255, 255])
940    
941              # the amount to scale the data by
942              scale = 255.0
943              numColours = len(data)
944    
945              # convert the colourmap into something vtk is more happy with
946              height = numarray.zeros(numColours, numarray.Float)
947              red = numarray.zeros(numColours, numarray.Float)
948              green = numarray.zeros(numColours, numarray.Float)
949              blue = numarray.zeros(numColours, numarray.Float)
950              for i in range(numColours):
951                  (h, r, g, b) = data[i]
952                  height[i] = float(h)
953                  red[i] = float(r)/scale
954                  green[i] = float(g)/scale
955                  blue[i] = float(b)/scale
956    
957              print "Loading bathymmetry data..."
958              # grab the z data
959              bathZData = self.bathymetry.getData()
960    
961              # get the origin
962              bathOrigin = self.bathymetry.getOrigin()
963              # get the delta_x and delta_y
964              bathSpacing = self.bathymetry.getSpacing()
965    
966              # now construct the x and y data from the spacing and the origin
967              numXPoints = bathZData.shape[1]
968              numYPoints = bathZData.shape[0]
969              numPoints = numXPoints*numYPoints
970    
971              bathXData = numarray.zeros(numXPoints, numarray.Float)
972              bathYData = numarray.zeros(numYPoints, numarray.Float)
973              for i in range(numXPoints):
974                  bathXData[i] = bathOrigin[0] + i*bathSpacing[0]
975    
976              for i in range(numYPoints):
977                  bathYData[i] = bathOrigin[1] + i*bathSpacing[1]
978    
979              # calculate the appropriate window size
980              dataWidth = max(bathXData) - min(bathXData)
981              dataHeight = max(bathYData) - min(bathYData)
982              winHeight = 600
983              winWidth = int(dataWidth*float(winHeight)/dataHeight)
984    
985              print "Done loading bathymmetry data"
986    
987              ### now do the vtk stuff
988    
989              # make sure rendering is offscreen
990              factGraphics = vtk.vtkGraphicsFactory()
991              factGraphics.SetUseMesaClasses(1)
992    
993              factImage = vtk.vtkImagingFactory()
994              factImage.SetUseMesaClasses(1)
995    
996              # make the bathymmetry colourmap
997              transFunc = vtk.vtkColorTransferFunction()
998              transFunc.SetColorSpaceToRGB()
999              for i in range(numColours):
1000                  transFunc.AddRGBPoint(height[i], red[i], green[i], blue[i])
1001                  h = height[i]
1002                  while i < numColours-1 and h < height[i+1]:
1003                    h += 1
1004                    transFunc.AddRGBPoint(h, red[i], green[i], blue[i])
1005    
1006              # set up the lookup table for the wave data
1007              refLut = vtk.vtkLookupTable()
1008              self.lutTrans = vtk.vtkLookupTable()
1009              refLut.Build()
1010              alpha = 0.7   # alpha channel value
1011              for i in range(256):
1012                  (r,g,b,a) = refLut.GetTableValue(255-i)
1013                  if g == 1.0 and b < 0.5 and r < 0.5:
1014                      self.lutTrans.SetTableValue(i, r, g, b, 0.0)
1015                  else:
1016                      self.lutTrans.SetTableValue(i, r, g-0.2, b, alpha)
1017    
1018              print "Generating the bathymmetry vtk object..."
1019    
1020              # set up the renderer and the render window
1021              self.ren = vtk.vtkRenderer()
1022              self.renWin = vtk.vtkRenderWindow()
1023              self.renWin.AddRenderer(self.ren)
1024              self.renWin.SetSize(winWidth, winHeight)
1025              self.renWin.OffScreenRenderingOn()
1026    
1027              # make an unstructured grid
1028              bathGrid = vtk.vtkUnstructuredGrid()
1029    
1030              # make the points for the vtk data
1031              bathPoints = vtk.vtkPoints()
1032              bathPoints.SetNumberOfPoints(numPoints)
1033    
1034              # make the vtk bathymmetry data set
1035              bathData = vtk.vtkFloatArray()
1036              bathData.SetNumberOfComponents(1)
1037              bathData.SetNumberOfTuples(numPoints)
1038    
1039              # add the points and data
1040              count = 0
1041              for i in range(numXPoints):
1042                  for j in range(numYPoints):
1043                      bathPoints.InsertPoint(count, bathXData[i], bathYData[j], 0.0)
1044                      bathData.InsertTuple1(count, bathZData[j,i])
1045                      count += 1
1046            
1047              # set the data to the grid
1048              bathGrid.SetPoints(bathPoints)
1049              bathGrid.GetPointData().SetScalars(bathData)
1050    
1051              # do a delaunay on the grid
1052              bathDelaunay = vtk.vtkDelaunay2D()
1053              bathDelaunay.SetInput(bathGrid)
1054              bathDelaunay.SetTolerance(0.001)
1055              bathDelaunay.SetAlpha(2.5)
1056    
1057              # set up the mapper
1058              bathMapper = vtk.vtkDataSetMapper()
1059              bathMapper.SetInput(bathDelaunay.GetOutput())
1060              bathMapper.SetLookupTable(transFunc)
1061    
1062              # set up the actor
1063              bathActor = vtk.vtkActor()
1064              bathActor.SetMapper(bathMapper)
1065    
1066              self.ren.AddActor(bathActor)
1067    
1068              ### now add the coastline
1069              print "Loading the coastline data..."
1070    
1071              # make the coastline points
1072              coastPoints = vtk.vtkPoints()
1073    
1074              # make the coastline grid
1075              coastGrid = vtk.vtkUnstructuredGrid()
1076    
1077              # now point the points and lines into the grid
1078              totalCoastPoints = 0
1079              for polyline in self.coastline.polylines:
1080                  numPoints = len(polyline)
1081                  coastLine = vtk.vtkPolyLine()
1082                  coastLine.GetPointIds().SetNumberOfIds(numPoints)
1083                  j = 0
1084                  for point in polyline:
1085                      coastLine.GetPointIds().SetId(j, j+totalCoastPoints)
1086                      coastPoints.InsertNextPoint(point.long, point.lat, 0.0)
1087                      j += 1
1088                  coastGrid.InsertNextCell(coastLine.GetCellType(),
1089                          coastLine.GetPointIds())
1090                  totalCoastPoints += numPoints
1091    
1092              coastGrid.SetPoints(coastPoints)
1093    
1094              # make the coast's mapper
1095              coastMapper = vtk.vtkDataSetMapper()
1096              coastMapper.SetInput(coastGrid)
1097    
1098              # make its actor
1099              coastActor = vtk.vtkActor()
1100              coastActor.SetMapper(coastMapper)
1101              coastActor.GetProperty().SetColor(0,0,0)
1102    
1103              # add the actor to the renderer
1104              self.ren.AddActor(coastActor)
1105    
1106              # set up the actor for the wave
1107              self.waveActor = vtk.vtkActor()
1108    
1109              # add the actor to the renderer
1110              self.ren.AddActor(self.waveActor)
1111    
1112              print "Done loading the coastline data"
1113              
1114         def doStepPostprocessing(self, dt):         def doStepPostprocessing(self, dt):
1115          """          """
1116          Does any necessary postprocessing after each step          Does any necessary postprocessing after each step
# Line 901  class SurfMovie(Model): Line 1122  class SurfMovie(Model):
1122               saveVTK(self.__fn,h=self.wave_height)               saveVTK(self.__fn,h=self.wave_height)
1123               # vtkobj=...               # vtkobj=...
1124               # save(self.__frame_name)               # save(self.__frame_name)
1125    
1126                 # make a reader for the data
1127                 waveReader = vtk.vtkXMLUnstructuredGridReader()
1128                 waveReader.SetFileName(self.__fn)
1129                 waveReader.Update()
1130    
1131                 # make the grid
1132                 waveGrid = waveReader.GetOutput()
1133                 waveGrid.Update()
1134    
1135                 (zMin, zMax) = waveGrid.GetPointData().GetScalars().GetRange()
1136                 print "Wave height range %f - %f" % (zMin, zMax)
1137    
1138                 # make a mapper for the grid
1139                 waveMapper = vtk.vtkDataSetMapper()
1140                 waveMapper.SetInput(waveGrid)
1141                 waveMapper.SetLookupTable(self.lutTrans)
1142                 waveMapper.SetScalarRange(-self.max_height, self.max_height)
1143    
1144                 self.waveActor.SetMapper(waveMapper)
1145    
1146                 # do stuff here that can only be done on the first frame
1147                 if self.firstFrame:
1148                     # zoom in a bit
1149                     self.ren.GetActiveCamera().Zoom(2.0)
1150                     self.firstFrame = False
1151    
1152                 # render the window
1153                 self.renWin.Render()
1154    
1155                 # convert the render window to an image
1156                 win2imgFilter = vtk.vtkWindowToImageFilter()
1157                 win2imgFilter.SetInput(self.renWin)
1158    
1159                 # the image name
1160                 imgFname = "%s%f.pnm" % (self.__frame_name, self.__next_t)
1161      
1162                 # save the image to file
1163                 outWriter = vtk.vtkPNMWriter()
1164                 outWriter.SetInput(win2imgFilter.GetOutput())
1165                 outWriter.SetFileName(imgFname)
1166                 outWriter.Write()
1167                 print "Wrote %s" % imgFname
1168    
1169                 # helpful for debugging:
1170                 #os.system("display %s" % imgFname)
1171    
1172                 self.paramsFileString += "%s\n" % imgFname
1173                 self.imageFiles.append(imgFname)
1174    
1175               self.__next_t+=self.dt               self.__next_t+=self.dt
1176    
1177         def getSafeTimeStepSize(self,dt):         def getSafeTimeStepSize(self,dt):
# Line 920  class SurfMovie(Model): Line 1191  class SurfMovie(Model):
1191            image files.            image files.
1192            """            """
1193            # make the movie into self.filename            # make the movie into self.filename
1194            pass            self.paramsFileString += "END_INPUT\n"
1195              self.paramsFileString += "PIXEL HALF\n"
1196              self.paramsFileString += "ASPECT_RATIO 1\n"
1197    
1198              # write the params string to file
1199              fp = open("%s.params" % self.filename, "w")
1200              fp.write(self.paramsFileString + '\n')
1201              fp.close()
1202    
1203              print "Performing conversion to mpeg"
1204              convertString = "ppmtompeg %s.params" % self.filename
1205              result = os.system(convertString)
1206              if result != 0:
1207                  print "An error occurred in mpeg conversion"
1208    
1209              # now clean up the image files
1210              print "Removing temporary image files"
1211              os.unlink("%s.params" % self.filename)
1212              for fname in self.imageFiles:
1213                  os.unlink(fname)
1214    
1215  if __name__=="__main__":  if __name__=="__main__":
1216     from esys.escript.modelframe import Link,Simulation     from esys.escript.modelframe import Link,Simulation
# Line 954  if __name__=="__main__": Line 1243  if __name__=="__main__":
1243     src.end_lat=-12.     src.end_lat=-12.
1244     src.start_long=110.     src.start_long=110.
1245     src.end_long=120.     src.end_long=120.
1246     src.width=0.1     src.width=1.
1247     src.decay_zone=0.01     src.decay_zone=0.01
1248     src.amplitude=1.     src.amplitude=1.
1249    
# Line 974  if __name__=="__main__": Line 1263  if __name__=="__main__":
1263     sm.t=Link(sq,"t")     sm.t=Link(sq,"t")
1264     sm.dt=5000.     sm.dt=5000.
1265     sm.filename="movie.mpg"     sm.filename="movie.mpg"
1266       sm.max_height=Link(src,"amplitude")
1267        
1268     s=Simulation([sq,oc,b,oreg,src,ts,sm])     s=Simulation([sq,oc,b,oreg,src,ts,sm])
1269     s.writeXML()     # s.writeXML()
1270     s.run()     s.run()
1271    
1272    # vim: expandtab shiftwidth=4:

Legend:
Removed from v.251  
changed lines
  Added in v.280

  ViewVC Help
Powered by ViewVC 1.1.26