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

Legend:
Removed from v.259  
changed lines
  Added in v.264

  ViewVC Help
Powered by ViewVC 1.1.26