/[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 262 by gross, Tue Nov 29 07:14:41 2005 UTC revision 263 by cochrane, Tue Nov 29 09:55:16 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 881  class SurfMovie(Model): Line 883  class SurfMovie(Model):
883                                   west=15.,                                   west=15.,
884                                   filename="movie.mpg")                                   filename="movie.mpg")
885    
886               self.paramsFileString = "REFERENCE_FRAME DECODED\n"
887               self.paramsFileString += "FRAME_RATE 24\n"
888               self.paramsFileString += "OUTPUT %s\n" % self.filename
889               self.paramsFileString += "PATTERN IBBPBBPBB\n"
890               self.paramsFileString += "FORCE_ENCODE_LAST_FRAME\n"
891               self.paramsFileString += "GOP_SIZE 20\n"
892               self.paramsFileString += "BSEARCH_ALG CROSS2\n"
893               self.paramsFileString += "PSEARCH_ALG TWOLEVEL\n"
894               self.paramsFileString += "IQSCALE 10\n"
895               self.paramsFileString += "PQSCALE 11\n"
896               self.paramsFileString += "BQSCALE 16\n"
897               self.paramsFileString += "RANGE 8\n"
898               self.paramsFileString += "SLICES_PER_FRAME 1\n"
899               self.paramsFileString += "BASE_FILE_FORMAT PNM\n"
900               self.paramsFileString += "INPUT_DIR \n"
901               self.paramsFileString += "INPUT_CONVERT *\n"
902               self.paramsFileString += "INPUT\n"
903    
904               self.firstFrame = True
905    
906         def doInitialization(self):         def doInitialization(self):
907            """            """
908            Initializes the time integartion scheme            Initializes the time integration scheme
909            """            """
910            self.__fn=os.tempnam()+".xml"            self.__fn=os.tempnam()+".xml"
911            self.__frame_name=os.tempnam()            self.__frame_name=os.tempnam()
# Line 892  class SurfMovie(Model): Line 914  class SurfMovie(Model):
914            # self.bathymetry.getVTK()            # self.bathymetry.getVTK()
915            # wndow(south,west,north,east)            # wndow(south,west,north,east)
916    
917              # the bathymmetry colourmap
918          data = []
919          data.append([-8000, 0,   0,   0])
920          data.append([-7000, 0,   5,   25])
921          data.append([-6000, 0,   10,  50])
922          data.append([-5000, 0,   80,  125])
923          data.append([-4000, 0,   150, 200])
924          data.append([-3000, 86,  197, 184])
925          data.append([-2000, 172, 245, 168])
926          data.append([-1000, 211, 250, 211])
927          data.append([0,     16,  48,  16])
928          data.append([300,   70,  150, 50])
929          data.append([500,   146, 126, 60])
930          data.append([1000,  198, 178, 80])
931          data.append([1250,  250, 230, 100])
932          data.append([1500,  250, 234, 126])
933          data.append([1750,  252, 238, 152])
934          data.append([2000,  252, 243, 177])
935          data.append([2250,  253, 249, 216])
936          data.append([2500,  255, 255, 255])
937    
938          # the amount to scale the data by
939              scale = 255.0
940              numColours = len(data)
941    
942          # convert the colourmap into something vtk is more happy with
943              height = numarray.zeros(numColours, numarray.Float)
944              red = numarray.zeros(numColours, numarray.Float)
945              green = numarray.zeros(numColours, numarray.Float)
946              blue = numarray.zeros(numColours, numarray.Float)
947              for i in range(numColours):
948                  (h, r, g, b) = data[i]
949                  height[i] = float(h)
950                  red[i] = float(r)/scale
951                  green[i] = float(g)/scale
952                  blue[i] = float(b)/scale
953    
954              print "Loading bathymmetry data..."
955              # grab the z data
956              bathZData = self.bathymetry.getData()
957    
958              # get the origin
959              bathOrigin = self.bathymetry.getOrigin()
960              # get the delta_x and delta_y
961              bathSpacing = self.bathymetry.getSpacing()
962    
963              # now construct the x and y data from the spacing and the origin
964              numXPoints = bathZData.shape[1]
965              numYPoints = bathZData.shape[0]
966              numPoints = numXPoints*numYPoints
967    
968              bathXData = numarray.zeros(numXPoints, numarray.Float)
969              bathYData = numarray.zeros(numYPoints, numarray.Float)
970              for i in range(numXPoints):
971                  bathXData[i] = bathOrigin[0] + i*bathSpacing[0]
972    
973              for i in range(numYPoints):
974                  bathYData[i] = bathOrigin[1] + i*bathSpacing[1]
975    
976              # calculate the appropriate window size
977              dataWidth = max(bathXData) - min(bathXData)
978              dataHeight = max(bathYData) - min(bathYData)
979              winHeight = 600
980              winWidth = int(dataWidth*float(winHeight)/dataHeight)
981    
982              print "Done loading bathymmetry data"
983    
984              ### now do the vtk stuff
985    
986              # make sure rendering is offscreen
987              factGraphics = vtk.vtkGraphicsFactory()
988              factGraphics.SetUseMesaClasses(1)
989    
990              factImage = vtk.vtkImagingFactory()
991              factImage.SetUseMesaClasses(1)
992    
993              # make the bathymmetry colourmap
994              transFunc = vtk.vtkColorTransferFunction()
995              transFunc.SetColorSpaceToRGB()
996              for i in range(numColours):
997                  transFunc.AddRGBPoint(height[i], red[i], green[i], blue[i])
998                  h = height[i]
999                  while i < numColours-1 and h < height[i+1]:
1000                    h += 1
1001                    transFunc.AddRGBPoint(h, red[i], green[i], blue[i])
1002    
1003          # set up the lookup table for the wave data
1004          refLut = vtk.vtkLookupTable()
1005          self.lutTrans = vtk.vtkLookupTable()
1006          refLut.Build()
1007          alpha = 0.4   # alpha channel value
1008          for i in range(256):
1009              (r,g,b,a) = refLut.GetTableValue(255-i)
1010              if g == 1.0 and b < 0.5 and r < 0.5:
1011              self.lutTrans.SetTableValue(i, r, g, b, 0.0)
1012              else:
1013              self.lutTrans.SetTableValue(i, r, g-0.2, b, alpha)
1014    
1015              print "Generating the bathymmetry vtk object..."
1016    
1017              # set up the renderer and the render window
1018              self.ren = vtk.vtkRenderer()
1019              self.renWin = vtk.vtkRenderWindow()
1020              self.renWin.AddRenderer(self.ren)
1021              self.renWin.SetSize(winWidth, winHeight)
1022              self.renWin.OffScreenRenderingOn()
1023    
1024              # make an unstructured grid
1025              bathGrid = vtk.vtkUnstructuredGrid()
1026    
1027              # make the points for the vtk data
1028              bathPoints = vtk.vtkPoints()
1029              bathPoints.SetNumberOfPoints(numPoints)
1030    
1031              # make the vtk bathymmetry data set
1032              bathData = vtk.vtkFloatArray()
1033              bathData.SetNumberOfComponents(1)
1034              bathData.SetNumberOfTuples(numPoints)
1035    
1036              # add the points and data
1037              count = 0
1038              for i in range(numXPoints):
1039                  for j in range(numYPoints):
1040                      bathPoints.InsertPoint(count, bathXData[i], bathYData[j], 0.0)
1041                      bathData.InsertTuple1(count, bathZData[j,i])
1042                      count += 1
1043            
1044              # set the data to the grid
1045              bathGrid.SetPoints(bathPoints)
1046              bathGrid.GetPointData().SetScalars(bathData)
1047    
1048              # do a delaunay on the grid
1049              bathDelaunay = vtk.vtkDelaunay2D()
1050              bathDelaunay.SetInput(bathGrid)
1051              bathDelaunay.SetTolerance(0.001)
1052              bathDelaunay.SetAlpha(2.5)
1053    
1054              # set up the mapper
1055              bathMapper = vtk.vtkDataSetMapper()
1056              bathMapper.SetInput(bathDelaunay.GetOutput())
1057              bathMapper.SetLookupTable(transFunc)
1058    
1059              # set up the actor
1060              bathActor = vtk.vtkActor()
1061              bathActor.SetMapper(bathMapper)
1062    
1063              self.ren.AddActor(bathActor)
1064    
1065          ### now add the coastline
1066              print "Loading the coastline data..."
1067    
1068          # make the coastline points
1069          coastPoints = vtk.vtkPoints()
1070    
1071          # make the coastline grid
1072          coastGrid = vtk.vtkUnstructuredGrid()
1073    
1074          # now point the points and lines into the grid
1075          totalCoastPoints = 0
1076          for polyline in self.coastline.polylines:
1077              numPoints = len(polyline)
1078              coastLine = vtk.vtkPolyLine()
1079              coastLine.GetPointIds().SetNumberOfIds(numPoints)
1080              j = 0
1081              for point in polyline:
1082              coastLine.GetPointIds().SetId(j, j+totalCoastPoints)
1083              coastPoints.InsertNextPoint(point.long, point.lat, 0.0)
1084              j += 1
1085              coastGrid.InsertNextCell(coastLine.GetCellType(),
1086                  coastLine.GetPointIds())
1087              totalCoastPoints += numPoints
1088    
1089          coastGrid.SetPoints(coastPoints)
1090    
1091          # make the coast's mapper
1092          coastMapper = vtk.vtkDataSetMapper()
1093          coastMapper.SetInput(coastGrid)
1094    
1095          # make its actor
1096          coastActor = vtk.vtkActor()
1097          coastActor.SetMapper(coastMapper)
1098          coastActor.GetProperty().SetColor(0,0,0)
1099    
1100          # add the actor to the renderer
1101          self.ren.AddActor(coastActor)
1102    
1103              # set up the actor for the wave
1104              self.waveActor = vtk.vtkActor()
1105    
1106              # add the actor to the renderer
1107              self.ren.AddActor(self.waveActor)
1108    
1109              print "Done loading the coastline data"
1110              
1111         def doStepPostprocessing(self, dt):         def doStepPostprocessing(self, dt):
1112          """          """
1113          Does any necessary postprocessing after each step          Does any necessary postprocessing after each step
# Line 903  class SurfMovie(Model): Line 1119  class SurfMovie(Model):
1119               saveVTK(self.__fn,h=self.wave_height)               saveVTK(self.__fn,h=self.wave_height)
1120               # vtkobj=...               # vtkobj=...
1121               # save(self.__frame_name)               # save(self.__frame_name)
1122    
1123                 # get the wave data
1124                 waveDomain = self.wave_height.getDomain().getX()
1125                 waveX = waveDomain[0].convertToNumArray()
1126                 waveY = waveDomain[1].convertToNumArray()
1127                 waveZ = self.wave_height.convertToNumArray()
1128    
1129                 numPoints = len(waveZ)
1130    
1131                 # make the points
1132                 wavePoints = vtk.vtkPoints()
1133                 wavePoints.SetNumberOfPoints(numPoints)
1134    
1135                 # make the vtk data array
1136                 waveData = vtk.vtkFloatArray()
1137                 waveData.SetNumberOfComponents(1)
1138                 waveData.SetNumberOfTuples(numPoints)
1139                 waveData.SetName("data")
1140    
1141                 # put the data into the points and array
1142                 for i in range(numPoints):
1143                     wavePoints.InsertPoint(i, waveX[i], waveY[i], 0.0)
1144                     waveData.InsertTuple1(i, waveZ[i])
1145    
1146                 # make the grid
1147                 waveGrid = vtk.vtkUnstructuredGrid()
1148                 waveGrid.SetPoints(wavePoints)
1149                 waveGrid.GetPointData().AddArray(waveData)
1150                 waveGrid.GetPointData().SetActiveScalars("data")
1151    
1152                 (zMin, zMax) = waveGrid.GetPointData().GetScalars().GetRange()
1153                 print "Wave height range %f - %f" % (zMin, zMax)
1154    
1155             # do a delaunay on the data
1156             waveDelaunay = vtk.vtkDelaunay2D()
1157             waveDelaunay.SetInput(waveGrid)
1158             waveDelaunay.SetTolerance(0.001)
1159             waveDelaunay.SetAlpha(2.5)
1160    
1161                 # make a mapper for the grid
1162                 waveMapper = vtk.vtkDataSetMapper()
1163                 waveMapper.SetInput(waveDelaunay.GetOutput())
1164             waveMapper.SetLookupTable(self.lutTrans)
1165                 waveMapper.SetScalarRange(zMin, zMax)
1166    
1167                 self.waveActor.SetMapper(waveMapper)
1168    
1169                 # do stuff here that can only be done on the first frame
1170                 if self.firstFrame:
1171                     # zoom in a bit
1172                     self.ren.GetActiveCamera().Zoom(2.0)
1173                     self.firstFrame = False
1174    
1175                 # render the window
1176                 self.renWin.Render()
1177    
1178                 # convert the render window to an image
1179                 win2imgFilter = vtk.vtkWindowToImageFilter()
1180                 win2imgFilter.SetInput(self.renWin)
1181    
1182                 # the image name
1183                 imgFname = "%s%f.pnm" % (self.__frame_name, self.__next_t)
1184      
1185                 # save the image to file
1186                 outWriter = vtk.vtkPNMWriter()
1187                 outWriter.SetInput(win2imgFilter.GetOutput())
1188                 outWriter.SetFileName(imgFname)
1189                 outWriter.Write()
1190                 print "Wrote %s" % imgFname
1191                 os.system("display %s" % imgFname)
1192                 self.paramsFileString += "%s\n" % imgFname
1193    
1194               self.__next_t+=self.dt               self.__next_t+=self.dt
1195    
1196         def getSafeTimeStepSize(self,dt):         def getSafeTimeStepSize(self,dt):
# Line 922  class SurfMovie(Model): Line 1210  class SurfMovie(Model):
1210            image files.            image files.
1211            """            """
1212            # make the movie into self.filename            # make the movie into self.filename
1213            pass            self.paramsFileString += "END_INPUT\n"
1214              self.paramsFileString += "PIXEL HALF\n"
1215              self.paramsFileString += "ASPECT_RATIO 1\n"
1216    
1217              # write the params string to file
1218              fp = open("%s.params" % self.filename, "w")
1219              fp.write(self.paramsFileString + '\n')
1220              fp.close()
1221    
1222              print "Performing conversion to mpeg"
1223              convertString = "ppmtompeg %s.params" % self.filename
1224              result = os.system(convertString)
1225              if result != 0:
1226                  print "An error occurred in mpeg conversion"
1227    
1228  if __name__=="__main__":  if __name__=="__main__":
1229     from esys.escript.modelframe import Link,Simulation     from esys.escript.modelframe import Link,Simulation
# Line 980  if __name__=="__main__": Line 1280  if __name__=="__main__":
1280     s=Simulation([sq,oc,b,oreg,src,ts,sm])     s=Simulation([sq,oc,b,oreg,src,ts,sm])
1281     # s.writeXML()     # s.writeXML()
1282     s.run()     s.run()
1283    
1284    # vim: expandtab shiftwidth=4:
Legend:
Removed from v.262  
changed lines
  Added in v.263
  ViewVC Help
Powered by ViewVC 1.1.26