/[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 276 by cochrane, Wed Nov 30 04:31:40 2005 UTC
# Line 1  Line 1 
1  # $Id$  #!/usr/bin/env python
2    
3    # $Id$
4    
5  import os  import os, sys
6    import vtk
7  from esys.escript import *  from esys.escript import *
8  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
9  from esys.escript.modelframe import Model  from esys.escript.modelframe import Model
# Line 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             # make a reader for the data
1126             waveReader = vtk.vtkXMLUnstructuredGridReader()
1127             waveReader.SetFileName(self.__fn)
1128             waveReader.Update()
1129    
1130                 # make the grid
1131                 waveGrid = waveReader.GetOutput()
1132             waveGrid.Update()
1133    
1134                 (zMin, zMax) = waveGrid.GetPointData().GetScalars().GetRange()
1135                 print "Wave height range %f - %f" % (zMin, zMax)
1136    
1137                 # make a mapper for the grid
1138                 waveMapper = vtk.vtkDataSetMapper()
1139                 waveMapper.SetInput(waveGrid)
1140             waveMapper.SetLookupTable(self.lutTrans)
1141                 waveMapper.SetScalarRange(zMin, zMax)
1142    
1143                 self.waveActor.SetMapper(waveMapper)
1144    
1145                 # do stuff here that can only be done on the first frame
1146                 if self.firstFrame:
1147                     # zoom in a bit
1148                     self.ren.GetActiveCamera().Zoom(2.0)
1149                     self.firstFrame = False
1150    
1151                 # render the window
1152                 self.renWin.Render()
1153    
1154                 # convert the render window to an image
1155                 win2imgFilter = vtk.vtkWindowToImageFilter()
1156                 win2imgFilter.SetInput(self.renWin)
1157    
1158                 # the image name
1159                 imgFname = "%s%f.pnm" % (self.__frame_name, self.__next_t)
1160      
1161                 # save the image to file
1162                 outWriter = vtk.vtkPNMWriter()
1163                 outWriter.SetInput(win2imgFilter.GetOutput())
1164                 outWriter.SetFileName(imgFname)
1165                 outWriter.Write()
1166                 print "Wrote %s" % imgFname
1167    
1168             # helpful for debugging:
1169             #os.system("display %s" % imgFname)
1170    
1171                 self.paramsFileString += "%s\n" % imgFname
1172             self.imageFiles.append(imgFname)
1173    
1174               self.__next_t+=self.dt               self.__next_t+=self.dt
1175    
1176         def getSafeTimeStepSize(self,dt):         def getSafeTimeStepSize(self,dt):
# Line 922  class SurfMovie(Model): Line 1190  class SurfMovie(Model):
1190            image files.            image files.
1191            """            """
1192            # make the movie into self.filename            # make the movie into self.filename
1193            pass            self.paramsFileString += "END_INPUT\n"
1194              self.paramsFileString += "PIXEL HALF\n"
1195              self.paramsFileString += "ASPECT_RATIO 1\n"
1196    
1197              # write the params string to file
1198              fp = open("%s.params" % self.filename, "w")
1199              fp.write(self.paramsFileString + '\n')
1200              fp.close()
1201    
1202              print "Performing conversion to mpeg"
1203              convertString = "ppmtompeg %s.params" % self.filename
1204              result = os.system(convertString)
1205              if result != 0:
1206                  print "An error occurred in mpeg conversion"
1207    
1208          # now clean up the image files
1209          print "Removing temporary image files"
1210          os.unlink("%s.params" % self.filename)
1211          for fname in self.imageFiles:
1212              os.unlink(fname)
1213    
1214  if __name__=="__main__":  if __name__=="__main__":
1215     from esys.escript.modelframe import Link,Simulation     from esys.escript.modelframe import Link,Simulation
# Line 956  if __name__=="__main__": Line 1242  if __name__=="__main__":
1242     src.end_lat=-12.     src.end_lat=-12.
1243     src.start_long=110.     src.start_long=110.
1244     src.end_long=120.     src.end_long=120.
1245     src.width=0.1     src.width=1.
1246     src.decay_zone=0.01     src.decay_zone=0.01
1247     src.amplitude=1.     src.amplitude=1.
1248    
# Line 980  if __name__=="__main__": Line 1266  if __name__=="__main__":
1266     s=Simulation([sq,oc,b,oreg,src,ts,sm])     s=Simulation([sq,oc,b,oreg,src,ts,sm])
1267     # s.writeXML()     # s.writeXML()
1268     s.run()     s.run()
1269    
1270    # vim: expandtab shiftwidth=4:

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

  ViewVC Help
Powered by ViewVC 1.1.26