#!/usr/bin/env python ## Import some modules import vtk, time, os, sys, signal from numpy import * dir1 = "/home/scratch/mstorti/TO-BACKUP/DATA-CD/AHMED-FSI-RESULTS" ## Parameters of run frame1 = 1 if 1: frame2 = 400 nfiles = 174 dire = dir1 + "/STEPS-2009-JUL-03-yaw-cartouche-centerbody" frame_rate = 8 else: frame2 = 700 nfiles = 323 dire = dir1 + "/STEPS-2009-07-06-shiftx" frame_rate = 16 do_cut_mesh = 0 mkvideo = 1 frame_inc = 4 static_cam = 0 xscale = 1.0 # zplane = -0.025 # low zplane = 0.14 # center body approx frame_forced = None if mkvideo: frame_inc = 1 def my_sigint_handler(f,g): global frame_forced print "Press k to kill, anything else to continue" while True: s = raw_input("Enter something: ") if s=='k': print "Exiting" sys.exit() elif s=='c': print "Continuing" frame_forced = frame - 2 if frame_forced<0: frame_forced = 0 break elif s=='n': frame_forced = int(raw_input("Enter frame: ")) print "Continuing from frame %d" % frame break else: print "Unknown answer \"%s\"" % s signal.signal(signal.SIGINT,my_sigint_handler) def draw_axes(xc,L,D): cyl = vtk.vtkCylinderSource() cyl.SetHeight(L) #cyl.SetCenter(xc) cyl.SetRadius(D/2.0) cyl.SetResolution(30) cylMapper = vtk.vtkPolyDataMapper() cylMapper.SetInputConnection(cyl.GetOutputPort() ) ## Create the cylinder actor cylY = vtk.vtkActor() cylY.SetMapper(cylMapper) cylY.SetPosition(xc) cylX = vtk.vtkActor() cylX.SetMapper(cylMapper) cylX.RotateZ(90.0) cylX.SetPosition(xc) cylZ = vtk.vtkActor() cylZ.SetMapper(cylMapper) cylZ.RotateX(90.0) cylZ.SetPosition(xc) cylX.GetProperty().SetColor(1.0,0.0,0.0) cylY.GetProperty().SetColor(0.0,1.0,0.0) cylZ.GetProperty().SetColor(0.0,0.0,1.0) return [cylX,cylY,cylZ] def state_indx(frame): indx = (frame_rate*frame) % (4*(nfiles-1)) xfac = 1.0 if indx>=2*(nfiles-1): xfac = -1.0 indx -= 2*(nfiles-1) if indx>=nfiles: indx = 2*(nfiles-1)-indx return indx,xfac plane = vtk.vtkPlane() plane.SetNormal(0.0, 0.0, -1.0) plane.SetOrigin(0.0, 0.0, zplane) def read_file_cached(filein): assert os.path.exists(filein) t0 = os.path.getmtime(filein) npfile = filein + ".npz" np_ok = os.path.exists(npfile) if np_ok: t1 = os.path.getmtime(npfile) np_ok = t1>t0 if not np_ok: print "saving on file %s" % npfile a = loadtxt(filein) savez(npfile,a) return a else: print "reading file %s" % npfile a = load(npfile) return a['arr_0'] def make_nodes(x,nodes): ids = vtk.vtkIdList() # ids of particles nnod = x.shape[0] for k in range(0,nnod): nodes.InsertNextPoint(x[k]) ids.InsertNextId(k) def set_coords_nodes(x,nodes,xfac): nnod = x.shape[0] for k in range(0,nnod): xx = x[k] xx[0] *= xfac nodes.SetPoint(k,xx) nodes.Modified() def set_mesh_connectivity(mesh,conec,elem): nelem = conec.shape[0] nel = conec.shape[1] # mesh.Allocate(nelem,1000) for k in range(0,nelem): tids = elem.GetPointIds() for l in range(0,nel): tids.SetId(l,conec[k,l]-1) mesh.InsertNextCell(elem.GetCellType(), elem.GetPointIds()) x0 = read_file_cached("./ahmed_125.nod.tmp") print "x0 shape: ", x0.shape nodes = vtk.vtkPoints() # positions of particles make_nodes(x0,nodes) print "after loading nodes ",time.asctime() conec = read_file_cached("./ahmed_125.con") print "conec shape: ", conec.shape print "after reading ",time.asctime() volmesh = vtk.vtkUnstructuredGrid() volmesh.SetPoints(nodes) elem = vtk.vtkTetra() set_mesh_connectivity(volmesh,conec,elem) print "after loading elems ",time.asctime() mapper = vtk.vtkDataSetMapper() mapper.SetInput(volmesh) mapper.ScalarVisibilityOff() mapper.AddClippingPlane(plane) mesh_actor = vtk.vtkActor() mesh_actor.SetMapper(mapper) mesh_actor.GetProperty().SetDiffuseColor(1, 1, 1) cutter = vtk.vtkCutter() cutter.SetInput(volmesh) cutter.SetCutFunction(plane); cutmapper = vtk.vtkDataSetMapper() cutmapper.SetInputConnection(cutter.GetOutputPort()) cutmapper.ScalarVisibilityOff() cutactor = vtk.vtkActor() cutactor.SetMapper(cutmapper) cutactor.GetProperty().SetDiffuseColor(0, 1, 0) prop = cutactor.GetProperty() #prop.SetRepresentationToWireframe() prop.EdgeVisibilityOn() prop.SetLineWidth(2) trimesh = vtk.vtkUnstructuredGrid() tricon = read_file_cached("./ahmed.surf-con2.tmp") print "tri con shape: ", tricon.shape trimesh.SetPoints(nodes) elem = vtk.vtkTriangle() set_mesh_connectivity(trimesh,tricon,elem) print "after loading surf elems ",time.asctime() srfmapper = vtk.vtkDataSetMapper() srfmapper.SetInput(trimesh) srfmapper.ScalarVisibilityOff() srf_actor = vtk.vtkActor() srf_actor.SetMapper(srfmapper) prop = srf_actor.GetProperty() prop.SetDiffuseColor(1, 1, 1) prop.SetOpacity(0.5) if not do_cut_mesh: # prop.EdgeVisibilityOn() # prop.SetLineWidth(1) prop.SetDiffuseColor(0.8, 0.8, 0) prop.SetOpacity(1) ## Draw axes axes = draw_axes([0.,0.5,0.14],10,0.013) # Create the usual rendering stuff. ren = vtk.vtkRenderer() renWin = vtk.vtkRenderWindow() renWin.AddRenderer(ren) renWin.SetSize(640,480) ren.SetBackground(.1, .2, .4) if do_cut_mesh: ren.AddActor(cutactor) ren.AddActor(mesh_actor) ren.AddActor(srf_actor) ren.AddActor(axes[0]) ren.AddActor(axes[1]) ren.AddActor(axes[2]) ## This is used to store the frames ## for creating a movie w2i = vtk.vtkWindowToImageFilter() w2i.SetInput(renWin) w2i.Update() ## The TIFF writer writer = vtk.vtkTIFFWriter() writer.SetInputConnection(w2i.GetOutputPort()) writer.SetCompressionToJPEG() ## Initialize camera cam = ren.GetActiveCamera() to = array([0.,0.5,0.1]) cam.SetFocalPoint(to) cam.SetViewUp(0.,0.,1.); ## The flying camera is as follows. ## It rotates around the focal point but also the distance to the ## object varies sinusoidally between Rcammin and Rcammax Rcammax = 2.5 Rcammin = 1.5 phi0 = 45.0*pi/180.0 ## The speed of rotation of the camera wcam = 0.005 if not do_cut_mesh: wcam = 0.015 ## This is for a static cam if static_cam: Rcammax = 1 Rcammin = Rcammax wcam = 0 ## The speed with which the radius varies wRcam = 1.33*wcam # Initial frame frame = frame1 ## This is added so that it gives time to set ## no border in the OpenGL window and other stuff ## like minimizing other windows. if mkvideo: renWin.Render() raw_input("Enter something to continue: ") # Loop while: rotating the camera and modify # node coordinates while 1: time.sleep(0.3) print "frame ",frame, time.asctime() ## Update camera phi = phi0+ wcam*frame Rcamav =(Rcammin+Rcammax)/2.0 DRcam = (Rcammin-Rcammax)/2.0 Rcam = Rcamav + DRcam*sin(wRcam*frame) Hcam = 1*Rcam view_dir = array([Rcam*cos(phi),Rcam*sin(phi),Hcam]) pfrom = to + view_dir if do_cut_mesh: (indx,xfac) = state_indx(frame) dxfile = "%s/ahmed.mmv-state-%d.tmp.gz" % (dire,indx) if not os.path.exists(dxfile): break dx = loadtxt(dxfile) x = x0 + xscale * dx set_coords_nodes(x,nodes,xfac) cam.SetPosition(pfrom) renWin.Render() ren.ResetCameraClippingRange() ## Save current frame if mkvideo: assert os.path.isdir("./YUV") w2i.Modified() tiff = "./YUV/frame.%d.tiff" % frame yuv = "./YUV/frame.%d.yuv" % frame writer.SetFileName(tiff) writer.Write() os.system("convert %s %s ; gzip -f %s" % (tiff,yuv,yuv)) os.unlink(tiff) if frame == frame2: break ## Update frame counter frame += frame_inc if frame_forced is not None: frame = frame_forced frame_forced = None