[cig-commits] r4450 - in mc/3D/CitcomS/trunk/visual/Mayavi2: .
plugins plugins/filter
maweier at geodynamics.org
maweier at geodynamics.org
Tue Aug 29 16:44:53 PDT 2006
Author: maweier
Date: 2006-08-29 16:44:53 -0700 (Tue, 29 Aug 2006)
New Revision: 4450
Modified:
mc/3D/CitcomS/trunk/visual/Mayavi2/citcoms_hdf2vtk.py
mc/3D/CitcomS/trunk/visual/Mayavi2/plugins/CitcomSHDFUgrid.py
mc/3D/CitcomS/trunk/visual/Mayavi2/plugins/CitcomS_hdf_file_reader.py
mc/3D/CitcomS/trunk/visual/Mayavi2/plugins/filter/CitcomSshowCaps.py
Log:
More stable versions. Some performance enhancements.
Modified: mc/3D/CitcomS/trunk/visual/Mayavi2/citcoms_hdf2vtk.py
===================================================================
--- mc/3D/CitcomS/trunk/visual/Mayavi2/citcoms_hdf2vtk.py 2006-08-29 22:33:43 UTC (rev 4449)
+++ mc/3D/CitcomS/trunk/visual/Mayavi2/citcoms_hdf2vtk.py 2006-08-29 23:44:53 UTC (rev 4450)
@@ -40,11 +40,25 @@
create_bottom = False
create_surface = False
create_ascii = False
-
+nx = None
+ny = None
+nz = None
nx_redu=None
ny_redu=None
nz_redu=None
+el_nx_redu = None
+el_ny_redu = None
+el_nz_redu = None
+radius_inner = None
+radius_outer = None
+nproc_surf = None
+#Filehandler to the HDF file
+f = None
+#####################
+polygons3d = [] # arrays containing connectivity information
+polygons2d = []
+counter=0 #Counts iterations of citcom2vtk
def print_help():
print "Program to convert CitcomS HDF to Vtk files.\n"
@@ -60,114 +74,10 @@
print "-b, --bottom \n\t Set to export Bottom information to Vtk file."
print "-s, --surface \n\t Set to export Surface information to Vtk file."
print "-c, --createtopo \n\t Set to create topography information in bottom and surface Vtk file."
- print "-a, --ascii \n\t Create Vtk ASCII encoded files."
+ print "-a, --ascii \n\t Create Vtk ASCII encoded files instead if binary."
print "-h, --help, -? \n\t Print this help."
-# parse command line parameters
-try:
- opts, args = getopt(sys.argv[1:], "p:o:i:t:x:y:z:bscah?", ['path=','output=','timestep=','x=','y=','z=','bottom','surface','createtopo','ascii', 'help','?'])
-except GetoptError, msg:
- print "Error: %s" % msg
- sys.exit(1)
-
-if len(opts)<=1:
- print_help()
- sys.exit(0)
-for opt,arg in opts:
- if opt in ('-p','--path'):
- path = arg
-
- if opt in ('-o','--output'):
- vtk_path = arg
-
- if opt in ('-i','--initial'):
- try:
- initial = int(arg)
- except ValueError:
- print "Initial is not a number."
- sys.exit(1)
- if opt in ('-t','--timestep'):
- try:
- timesteps = int(arg)
- except ValueError:
- print "Timestep is not a number."
- sys.exit(1)
- if opt in ('-x','--nx_reduce'):
- try:
- nx_redu = int(arg)
- except ValueError:
- print "NX is not a number."
-
- if opt in ('-y','--ny_reduce'):
- try:
- ny_redu = int(arg)
- except ValueError:
- print "NY is not a number."
-
- if opt in ('-z','--nz_reduce'):
- try:
- nz_redu = int(arg)
- except ValueError:
- print "NZ is not a number."
-
- if opt in ('-b','--bottom'):
- create_bottom = True
- if opt in ('-s','--surface'):
- create_surface = True
- if opt in ('-c','--createtopo'):
- create_topo = True
- if opt in ('-a','--ascii'):
- create_ascii = True
- if opt in ('-h','--help'):
- print_help()
- sys.exit(0)
- if opt == '-?':
- print_help()
- sys.exit(0)
-
-#############################globals############################################
-
-f = tables.openFile(path,'r')
-
-
-nx = int(f.root.input._v_attrs.nodex)
-ny = int(f.root.input._v_attrs.nodey)
-nz = int(f.root.input._v_attrs.nodez)
-
-#If not defined as argument read from hdf
-hdf_timesteps = int(f.root.time.nrows)
-
-if timesteps==None or timesteps>hdf_timesteps:
- timesteps = hdf_timesteps
-
-
-if nx_redu:
- nx_redu = nx-1
-if ny_redu:
- ny_redu = ny-1
-if nz_redu:
- nz_redu = nz-1
-
-if nx_redu>=nx:
- nx_redu=nx-1
-if ny_redu>=ny:
- ny_redu=ny-1
-if nz_redu>=nz:
- nz_redu=nz-1
-
-el_nx_redu = nx_redu+1
-el_ny_redu = ny_redu+1
-el_nz_redu = nz_redu+1
-
-radius_inner = float(f.root.input._v_attrs.radius_inner)
-radius_outer = float(f.root.input._v_attrs.radius_outer)
-#rayleigh = float(f.root.input._v_attrs.rayleigh) #Important for creating the right lookup table
-nproc_surf = int(f.root.input._v_attrs.nproc_surf)
-
-
-###############################################################################
-
#Iterator for CitcomDataRepresentation(yxz) to VTK(xyz)
def vtk_iter(nx,ny,nz):
for i in xrange(nx):
@@ -216,15 +126,11 @@
return x, y, z
-polygons3d = [] # arrays containing connectivity information
-polygons2d = []
-counter=0 #Counts iterations of citcom2vtk
-
#Reads Citcom Files and creates a VTK File
def citcom2vtk(t):
print "Timestep:",t
-
+
benchmarkstr = ""
#Assign create_bottom and create_surface to bottom and surface
#to make them valid in methods namespace
@@ -268,7 +174,6 @@
hdf_temperature = cap.temperature[t]
hdf_viscosity = cap.viscosity[t]
-
###Benchmark Point 1 Stop##
#delta = datetime.now() - start
#benchmarkstr += "%.5lf," % (delta.seconds + float(delta.microseconds)/1e6)
@@ -293,12 +198,12 @@
for k in xrange(el_nz_redu):
k_redu = nz_redu_iter.next()
- thet , phi, r = map(float,hdf_coords[i_redu][j_redu][k_redu])
+ thet , phi, r = map(float,hdf_coords[i_redu,j_redu,k_redu])
temp_coords.append((thet,phi,r))
- vel_colat, vel_lon , vel_r = map(float,hdf_velocity[i_redu][j_redu][k_redu])
- temperature = float(hdf_temperature[i_redu][j_redu][k_redu])
- visc = float(hdf_viscosity[i_redu][j_redu][k_redu])
+ vel_colat, vel_lon , vel_r = map(float,hdf_velocity[i_redu,j_redu,k_redu])
+ temperature = float(hdf_temperature[i_redu,j_redu,k_redu])
+ visc = float(hdf_viscosity[i_redu,j_redu,k_redu])
temp_vel.append((vel_colat,vel_lon,vel_r))
temp_temp.append(temperature)
@@ -335,7 +240,7 @@
ordered_temperature.append(temp_temp[iter])
ordered_visc.append(temp_visc[iter])
-
+
###Benchmark Point 3 Stop##
#delta = datetime.now() - start
#benchmarkstr += "%.5lf," % (delta.seconds + float(delta.microseconds)/1e6)
@@ -373,7 +278,7 @@
print "\tCould not find surface information in file.\n \
Set create surface to false"
surface = False
-
+
###Benchmark Point 4 Stop##
#delta = datetime.now() - start
#benchmarkstr += "%.5lf," % (delta.seconds + float(delta.microseconds)/1e6)
@@ -397,6 +302,8 @@
for i in xrange(ny):
botm_mean += numpy.mean(hdf_bottom_topography[i])
botm_mean = botm_mean/ny
+
+
#print "Mean Surface:",surf_mean
###Benchmark Point 5 Stop##
@@ -435,20 +342,20 @@
if surface==True:
#Surface Information
if create_topo==True:
- colat,lon = hdf_surface_coord[i][j]
+ colat,lon = hdf_surface_coord[i,j]
#637100 = Earth radius, 33000 = ?
- x,y,z = RTF2XYZ(colat,lon,radius_outer+float( (hdf_surface_topography[i][j]-surf_mean)*(10**21)/(6371000**2/10**(-6))/(3300*10)/1000 ))
+ x,y,z = RTF2XYZ(colat,lon,radius_outer+float( (hdf_surface_topography[i,j]-surf_mean)*(10**21)/(6371000**2/10**(-6))/(3300*10)/1000 ))
surf_points.append((x,y,z))
else:
- colat, lon = hdf_surface_coord[i][j]
+ colat, lon = hdf_surface_coord[i,j]
x,y,z = RTF2XYZ(colat, lon,radius_outer)
surf_points.append((x,y,z))
#Surface Heatflux
- surf_hflux.append(float(hdf_surface_heatflux[i][j]))
+ surf_hflux.append(float(hdf_surface_heatflux[i,j]))
#Surface Velocity
- vel_colat, vel_lon = map(float,hdf_surface_velocity[i][j])
+ vel_colat, vel_lon = map(float,hdf_surface_velocity[i,j])
x,y,z = velocity2cart(vel_colat,vel_lon, radius_outer, colat, lon, radius_outer)
surf_vec.append((x,y,z))
@@ -528,7 +435,7 @@
#start = datetime.now()
############################
-
+ print 'Writing data to vtk...'
#Surface Points
if surface==True:
struct_coords = pyvtk.UnstructuredGrid(surf_points, pixel=polygons2d)
@@ -601,23 +508,158 @@
#print "\n"
+
+# parse command line parameters
+def initialize():
+ global initial
+ global timesteps
+ global create_topo
+ global create_bottom
+ global create_surface
+ global create_ascii
+ global nx
+ global ny
+ global nz
+ global nx_redu
+ global ny_redu
+ global nz_redu
+ global el_nx_redu
+ global el_ny_redu
+ global el_nz_redu
+ global radius_inner
+ global radius_outer
+ global nproc_surf
+ global f
+
+ try:
+ opts, args = getopt(sys.argv[1:], "p:o:i:t:x:y:z:bscah?", ['path=','output=','timestep=','x=','y=','z=','bottom','surface','createtopo','ascii', 'help','?'])
+ except GetoptError, msg:
+ print "Error: %s" % msg
+ sys.exit(1)
+
+ if len(opts)<=1:
+ print_help()
+ sys.exit(0)
+
+ for opt,arg in opts:
+ if opt in ('-p','--path'):
+ path = arg
+
+ if opt in ('-o','--output'):
+ vtk_path = arg
+
+ if opt in ('-i','--initial'):
+ try:
+ initial = int(arg)
+ except ValueError:
+ print "Initial is not a number."
+ sys.exit(1)
+ if opt in ('-t','--timestep'):
+ try:
+ timesteps = int(arg)
+ except ValueError:
+ print "Timestep is not a number."
+ sys.exit(1)
+ if opt in ('-x','--nx_reduce'):
+ try:
+ nx_redu = int(arg)
+ except ValueError:
+ print "NX is not a number."
+
+ if opt in ('-y','--ny_reduce'):
+ try:
+ ny_redu = int(arg)
+ except ValueError:
+ print "NY is not a number."
+
+ if opt in ('-z','--nz_reduce'):
+ try:
+ nz_redu = int(arg)
+ except ValueError:
+ print "NZ is not a number."
+
+ if opt in ('-b','--bottom'):
+ create_bottom = True
+
+ if opt in ('-s','--surface'):
+ create_surface = True
+
+ if opt in ('-c','--createtopo'):
+ create_topo = True
+
+ if opt in ('-a','--ascii'):
+ create_ascii = True
+
+ if opt in ('-h','--help'):
+ print_help()
+ sys.exit(0)
+ if opt == '-?':
+ print_help()
+ sys.exit(0)
+
+
+ f = tables.openFile(path,'r')
+
+ nx = int(f.root.input._v_attrs.nodex)
+ ny = int(f.root.input._v_attrs.nodey)
+ nz = int(f.root.input._v_attrs.nodez)
+
+ #If not defined as argument read from hdf
+ hdf_timesteps = int(f.root.time.nrows)
+
+ if timesteps==None or timesteps>hdf_timesteps:
+ timesteps = hdf_timesteps
+
+
+ if nx_redu==None:
+ nx_redu = nx-1
+ if ny_redu==None:
+ ny_redu = ny-1
+ if nz_redu==None:
+ nz_redu = nz-1
+
+ if nx_redu>=nx:
+ nx_redu=nx-1
+ if ny_redu>=ny:
+ ny_redu=ny-1
+ if nz_redu>=nz:
+ nz_redu=nz-1
+
+ el_nx_redu = nx_redu+1
+ el_ny_redu = ny_redu+1
+ el_nz_redu = nz_redu+1
+
+ radius_inner = float(f.root.input._v_attrs.radius_inner)
+ radius_outer = float(f.root.input._v_attrs.radius_outer)
+ nproc_surf = int(f.root.input._v_attrs.nproc_surf)
+
+
###############################################################################
-d1 = datetime.now()
+def citcoms_hdf2vtk():
+ global counter
+
+ initialize()
+ d1 = datetime.now()
+ print "Converting Hdf to Vtk"
+ print "Initial:",initial, "Timesteps:",timesteps
+ print "NX:",el_nx_redu, "NY:",el_ny_redu, "NZ:", el_nz_redu
+ print "Create Bottom: ",create_bottom, " Create Surface: ", create_surface
+ print "Create Topography: ", create_topo
-print "Converting Hdf to Vtk"
-print "Initial: ",initial, "Timesteps:",timesteps
-print "NX:",el_nx_redu, "NY:",el_ny_redu, "NZ:", el_nz_redu
-print "Create Bottom: ",create_bottom, " Create Surface: ", create_surface
-print "Create Topography: ", create_topo
-
-for t in xrange(initial,timesteps):
+ for t in xrange(initial,timesteps):
start = datetime.now()
citcom2vtk(t)
counter+=1
delta = datetime.now() - start
print "\t%.3lf sec" % (delta.seconds + float(delta.microseconds)/1e6)
-d2 = datetime.now()
-f.close()
-print "Total: %d seconds" % (d2 - d1).seconds
+ d2 = datetime.now()
+ f.close()
+ print "Total: %d seconds" % (d2 - d1).seconds
###############################################################################
+
+
+
+
+if __name__ == '__main__':
+ citcoms_hdf2vtk()
\ No newline at end of file
Modified: mc/3D/CitcomS/trunk/visual/Mayavi2/plugins/CitcomSHDFUgrid.py
===================================================================
--- mc/3D/CitcomS/trunk/visual/Mayavi2/plugins/CitcomSHDFUgrid.py 2006-08-29 22:33:43 UTC (rev 4449)
+++ mc/3D/CitcomS/trunk/visual/Mayavi2/plugins/CitcomSHDFUgrid.py 2006-08-29 23:44:53 UTC (rev 4450)
@@ -6,6 +6,8 @@
class CitcomSHDFUgrid:
+ """This Class converts CitcomS hdf files to tvtk UnstructuredGrid Dataset Objects """
+
data = None
_nx = None
_ny = None
@@ -18,15 +20,22 @@
timesteps = None
frequency = None
- #Iterator for CitcomDataRepresentation(yxz) to VTK(xyz)
+ #Global because a Unstructured Grid can only hold one scalar value at a time
+ #but our hdf file reader plugin wants to be able to read both
+ __vtkordered_visc = tvtk.FloatArray()
+ __vtkordered_temp = tvtk.FloatArray()
+
+
def vtk_iter(self,nx,ny,nz):
+ """Iterator for CitcomDataRepresentation(yxz) to VTK(xyz)"""
for i in xrange(nx):
for j in xrange(ny):
for k in xrange(nz):
yield k + nz * i + nz * nx * j
- #Reduces the CitcomS grid
+
def reduce_iter(self,n,nredu):
+ """Iterator to reduce the CitcomS grid"""
i=0
n_f=float(n)
nredu_f=float(nredu)
@@ -38,6 +47,7 @@
def velocity2cart(self,vel_colat,vel_long,r, x, y, z):
+ """Converts vectors in spherical to cartesian coordiantes"""
x1 = r*sin(x)*cos(y)+vel_colat*cos(x)*cos(y)-vel_long*sin(y)
y1 = r*sin(x)*sin(y)+vel_colat*cos(x)*sin(y)+vel_long*cos(y)
z1 = r*cos(x)-vel_colat*sin(x)
@@ -46,6 +56,7 @@
#Converts Spherical to Cartesian Coordinates
def RTF2XYZ(self,thet, phi, r):
+ """Converts points from spherical to cartesian coordinates"""
x = r * sin(thet) * cos(phi)
y = r * sin(thet) * sin(phi)
z = r * cos(thet)
@@ -53,24 +64,17 @@
- def citcom2vtk(self,t,f,nproc_surf,nx_redu,ny_redu,nz_redu,bottom,surface):
- #Assign create_bottom and create_surface to bottom and surface
- #to make them valid in methods namespace
+ def __citcom2vtk(self,t,f,nproc_surf,nx_redu,ny_redu,nz_redu,bottom,surface):
+ """Method to convert one timestep from a hdf file to a Vtk file. This Method is used
+ by the method initialize. Initialize reads the necessary meta information from the hdf file"""
+
+ hexagrid = tvtk.UnstructuredGrid()
+ hexagrid.allocate(1,1)
- benchmarkstr = ""
- fd = open('/home/maweier/benchmark.txt','w')
-
- hexagrid = tvtk.UnstructuredGrid()
- surfPixelGrid = tvtk.UnstructuredGrid()
- botmPixelGrid = tvtk.UnstructuredGrid()
-
- vtkordered_temp = tvtk.FloatArray()
vtkordered_velo = tvtk.FloatArray()
- vtkordered_visc = tvtk.FloatArray()
- hexagrid.allocate(1,1)
- surfPixelGrid.allocate(1, 1)
- botmPixelGrid.allocate(1,1)
+ ordered_temp = []
+ ordered_visc = []
nx = self._nx
ny = self._ny
@@ -85,24 +89,9 @@
ordered_velocity = []
ordered_visc = []
- #Surface and Bottom Points
- #Initialize empty sequences
- surf_vec = []
- botm_vec = []
- surf_topo = []
- surf_hflux = []
- botm_topo = []
- botm_hflux = []
-
- surf_points = []
- botm_points = []
-
for capnr in xrange(nproc_surf):
+
- ###Benchmark Point 1 Start##
- d1 = datetime.now()
- ############################
-
cap = f.root._f_getChild("cap%02d" % capnr)
temp_coords = [] # reset Coordinates, Velocity, Temperature Sequence
@@ -111,15 +100,12 @@
temp_visc = []
#Information from hdf
- #This information needs to be read only once
hdf_coords = cap.coord[:]
-
hdf_velocity = cap.velocity[t]
hdf_temperature = cap.temperature[t]
hdf_viscosity = cap.viscosity[t]
-
-
+
#Create Iterator to change data representation
nx_redu_iter = self.reduce_iter(nx,nx_redu)
ny_redu_iter = self.reduce_iter(ny,ny_redu)
@@ -127,14 +113,7 @@
vtk_i = self.vtk_iter(el_nx_redu,el_ny_redu,el_nz_redu)
- ##Benchmark Point 1 Stop##
- delta = datetime.now() - d1
- benchmarkstr += "%.5lf," % (delta.seconds + float(delta.microseconds)/1e6)
-
- ###Benchmark Point 2 Start##
- d1 = datetime.now()
- ############################
-
+
# read citcom data - zxy (z fastest)
for j in xrange(el_ny_redu):
j_redu = ny_redu_iter.next()
@@ -144,12 +123,12 @@
nz_redu_iter = self.reduce_iter(nz,nz_redu)
for k in xrange(el_nz_redu):
k_redu = nz_redu_iter.next()
- thet , phi, r = map(float,hdf_coords[i_redu][j_redu][k_redu])
+ thet , phi, r = map(float,hdf_coords[i_redu,j_redu,k_redu])
temp_coords.append((thet,phi,r))
- vel_colat, vel_lon , vel_r = map(float,hdf_velocity[i][j][k])
- temperature = float(hdf_temperature[i][j][k])
- visc = float(hdf_viscosity[i][j][k])
+ vel_colat, vel_lon , vel_r = map(float,hdf_velocity[i,j,k])
+ temperature = float(hdf_temperature[i,j,k])
+ visc = float(hdf_viscosity[i,j,k])
temp_vel.append((vel_colat,vel_lon,vel_r))
temp_temp.append(temperature)
@@ -161,13 +140,6 @@
del hdf_temperature
del hdf_viscosity
- ##Benchmark Point 2 Stop##
- delta = datetime.now() - d1
- benchmarkstr += "%.5lf," % (delta.seconds + float(delta.microseconds)/1e6)
-
- ###Benchmark Point 3 Start##
- d1 = datetime.now()
- ############################
# rearange vtk data - xyz (x fastest).
for n0 in xrange(el_nz_redu*el_ny_redu*el_nx_redu):
@@ -179,147 +151,25 @@
x_coord, y_coord, z_coord = self.RTF2XYZ(colat,lon,r)
ordered_points.append((x_coord,y_coord,z_coord))
- #Get Vectors in Cartesian Coords from Velocity
+ ordered_temp.append(temp_temp[iter])
+ ordered_visc.append(temp_visc[iter])
+
+
+ #Get Vectors in Cartesian Coords from Velocity field
vel_colat,vel_lon,vel_r = temp_vel[iter]
x_velo, y_velo, z_velo = self.velocity2cart(vel_colat,vel_lon,vel_r, colat,lon , r)
ordered_velocity.append((x_velo,y_velo,z_velo))
################################################
- vtkordered_temp.insert_next_tuple1(temp_temp[iter])
- vtkordered_visc.insert_next_tuple1(temp_visc[iter])
-
- vtkordered_velo.from_array(ordered_velocity)
-
+
+
##Delete Unused Object for GC
del temp_coords
del temp_vel
del temp_temp
del temp_visc
- ##Benchmark Point 3 Stop##
- delta = datetime.now() - d1
- benchmarkstr += "%.5lf," % (delta.seconds + float(delta.microseconds)/1e6)
-
- ###Benchmark Point 4 Start##
- d1 = datetime.now()
- ############################
-
- #Bottom Information from hdf
- if bottom == True:
- try:
- hdf_bottom_coord = cap.botm.coord[:]
- hdf_bottom_heatflux = cap.botm.heatflux[t]
- hdf_bottom_topography = cap.botm.topography[t]
- hdf_bottom_velocity = cap.botm.velocity[t]
- except:
- bottom = False
- #Surface Information from hdf
- if surface==True:
- try:
- hdf_surface_coord = cap.surf.coord[:]
- hdf_surface_heatflux = cap.surf.heatflux[t]
- hdf_surface_topography = cap.surf.topography[t]
- hdf_surface_velocity = cap.surf.velocity[t]
- except:
- surface = False
-
- ##Benchmark Point 4 Stop##
- delta = datetime.now() - d1
- benchmarkstr += "%.5lf," % (delta.seconds + float(delta.microseconds)/1e6)
-
- ###Benchmark Point 5 Start##
- d1 = datetime.now()
- ############################
-
- #Compute surface/bottom topography mean
- if bottom==True or surface==True:
- surf_mean=0.0
- botm_mean=0.0
-
- for i in xrange(ny):
- if surface == True:
- surf_mean += numpy.mean(hdf_surface_topography[i])
- if bottom == True:
- botm_mean += numpy.mean(hdf_bottom_topography[i])
-
- surf_mean = surf_mean/ny
- botm_mean = botm_mean/ny
- #print "Mean Surface:",surf_mean
-
- ##Benchmark Point 5 Stop##
- delta = datetime.now() - d1
- benchmarkstr += "%.5lf," % (delta.seconds + float(delta.microseconds)/1e6)
-
- ###Benchmark Point 6 Start##
- d1 = datetime.now()
- ############################
-
- #Read Surface and Bottom Data
- if bottom==True or surface == True:
- for i in xrange(ny):
- for j in xrange(nx):
-
-
- if bottom==True:
- #Bottom Coordinates
- if create_topo==True:
- colat, lon = hdf_bottom_coord[i][j]
- x,y,z = self.RTF2XYZ(colat,lon,radius_inner+float( (hdf_bottom_topography[i][j]-botm_mean)*(10**21)/(6371000**2/10**(-6))/(3300*10)/1000 ))
- botm_points.append((x,y,z))
- else:
- colat, lon = hdf_bottom_coord[i][j]
- x,y,z = self.RTF2XYZ(colat, lon,radius_inner)
- botm_points.append((x,y,z))
-
- #Bottom Heatflux
- botm_hflux.append(float(hdf_bottom_heatflux[i][j]))
-
- #Bottom Velocity
- vel_colat, vel_lon = map(float,hdf_bottom_velocity[i][j])
- x,y,z = self.velocity2cart(vel_colat,vel_lon, radius_inner, colat, lon, radius_inner)
- botm_vec.append((x,y,z))
-
- if surface==True:
- #Surface Information
- if create_topo==True:
- colat,lon = hdf_surface_coord[i][j]
- #637100 = Earth radius, 33000 = ?
- x,y,z = self.RTF2XYZ(colat,lon,radius_outer+float( (hdf_surface_topography[i][j]-surf_mean)*(10**21)/(6371000**2/10**(-6))/(3300*10)/1000 ))
- surf_points.append((x,y,z))
- else:
- colat, lon = hdf_surface_coord[i][j]
- x,y,z = self.RTF2XYZ(colat, lon,radius_outer)
- surf_points.append((x,y,z))
-
- #Surface Heatflux
- surf_hflux.append(float(hdf_surface_heatflux[i][j]))
-
- #Surface Velocity
- vel_colat, vel_lon = map(float,hdf_surface_velocity[i][j])
- x,y,z = self.velocity2cart(vel_colat,vel_lon, radius_outer, colat, lon, radius_outer)
- surf_vec.append((x,y,z))
-
- vtk_botm_vec.from_array(botm_vec)
- vtk_surf_vec.from_array(surf_vec)
-
- #del variables for GC
- if bottom==True:
- del hdf_bottom_coord
- del hdf_bottom_heatflux
- del hdf_bottom_velocity
- if surface==True:
- del hdf_surface_coord
- del hdf_surface_heatflux
- del hdf_surface_velocity
-
- ##Benchmark Point 6 Stop##
- delta = datetime.now() - d1
- benchmarkstr += "%.5lf," % (delta.seconds + float(delta.microseconds)/1e6)
-
- ###Benchmark Point 7 Start##
- d1 = datetime.now()
- ############################
-##################################################################
+
#Create Connectivity info
if counter==0:
#For 3d Data
@@ -350,44 +200,26 @@
i+=1
- if bottom==True or surface==True:
- #Connectivity for 2d-Data
- i=1
- for n in xrange((nx)*(ny) - nx):
- if i%nx!=0 :
- n0 = n+(capnr*((nx)*(ny)))
- n1 = n0+1
- n2 = n0+ny
- n3 = n2+1
- surfPixelGrid.insert_next_cell(8 , [n0,n1,n2,n3])
- botmPixelGrid.insert_next_cell(8 , [n0,n1,n2,n3])
- i+=1
-
- ##Benchmark Point 7 Stop##
- delta = datetime.now() - d1
- benchmarkstr += "%.5lf \n" % (delta.seconds + float(delta.microseconds)/1e6)
- fd.write(benchmarkstr)
+ #Store Arrays in Vtk conform Datastructures
+ self.__vtkordered_temp.from_array(ordered_temp)
+ self.__vtkordered_visc.from_array(ordered_visc)
+ vtkordered_velo.from_array(ordered_velocity)
- benchmarkstr = '\n\nIO: '
- ###Benchmark Point IO Start##
- d1 = datetime.now()
- ############################
- vtkordered_temp.name = 'Temperature'
- hexagrid.point_data.scalars = vtkordered_temp
+ self.__vtkordered_temp.name = 'Temperature'
+ self.__vtkordered_visc.name = 'Viscosity'
+ hexagrid.point_data.scalars = self.__vtkordered_temp
vtkordered_velo.name = 'Velocity'
hexagrid.point_data.vectors = vtkordered_velo
hexagrid.points = ordered_points
- ##Benchmark Point IO Stop##
- delta = datetime.now() - d1
- benchmarkstr += "%.5lf" % (delta.seconds + float(delta.microseconds)/1e6)
- fd.write(benchmarkstr)
return hexagrid
def initialize(self,filename,timestep,nx_redu,ny_redu,nz_redu,surface,bottom):
-
+ """Call this method to convert a Citcoms Hdf file to a Vtk file"""
+
+ #Read meta-inforamtion
hdf=tables.openFile(filename,'r')
self._nx = int(hdf.root.input._v_attrs.nodex)
self._ny = int(hdf.root.input._v_attrs.nodey)
@@ -401,19 +233,24 @@
if nz_redu==0 or nz_redu>=self._nz:
nz_redu = self._nz-1
-
+ #Make reduction factors global
self._nx_redu = nx_redu
self._ny_redu = ny_redu
self._nz_redu = nz_redu
- #Number of Timesteps in scene \
+
+ #Number of Timesteps in scene
self.timesteps = int(hdf.root.time.nrows)
- #self._radius_inner = float(hdf.root.input._v_attrs.radius_inner) #Only important for displaying data.
- #self._radius_outer = float(hdf.root.input._v_attrs.radius_outer)
+ #Number of caps
nproc_surf = int(hdf.root.input._v_attrs.nproc_surf)
- hexgrid = self.citcom2vtk(timestep,hdf,nproc_surf,nx_redu,ny_redu,nz_redu,surface,bottom)
+ #start computation
+ hexgrid = self.__citcom2vtk(timestep,hdf,nproc_surf,nx_redu,ny_redu,nz_redu,surface,bottom)
hdf.close()
return hexgrid
-
\ No newline at end of file
+ def get_vtk_viscosity(self):
+ return self.__vtkordered_visc
+
+ def get_vtk_temperature(self):
+ return self.__vtkordered_temp
\ No newline at end of file
Modified: mc/3D/CitcomS/trunk/visual/Mayavi2/plugins/CitcomS_hdf_file_reader.py
===================================================================
--- mc/3D/CitcomS/trunk/visual/Mayavi2/plugins/CitcomS_hdf_file_reader.py 2006-08-29 22:33:43 UTC (rev 4449)
+++ mc/3D/CitcomS/trunk/visual/Mayavi2/plugins/CitcomS_hdf_file_reader.py 2006-08-29 23:44:53 UTC (rev 4450)
@@ -17,7 +17,9 @@
from enthought.mayavi.core.common import handle_children_state
from enthought.mayavi.sources.vtk_xml_file_reader import get_all_attributes
import tables
+import re
+
######################################################################
# `CitcomSVTKDataSource` class
######################################################################
@@ -30,20 +32,30 @@
# The VTK dataset to manage.
data = Instance(tvtk.DataSet)
+ # Class to convert Hdf to Vtk Unstructured Grid Objects
citcomshdftougrid = CitcomSHDFUgrid()
current_timestep = Int(0)
+ #To support changing the Scalar values in Mayavi2
+ temperature = Instance(tvtk.FloatArray())
+ viscosity = Instance(tvtk.FloatArray())
+
+ #Resolution comming from Hdf file
nx = Int()
ny = Int()
nz = Int()
+ #Current reduced resolution. User defined
nx_redu = Int()
ny_redu = Int()
nz_redu = Int()
+ #Number of timesteps in Hdf
timesteps = Int()
- frequency = Int()
+
+ #Button to trigger the process of reading a timestep
read_timestep = Button('Read timestep')
filename = Str()
+
########################################
# Dynamic traits: These traits are dummies and are dynamically
# updated depending on the contents of the file.
@@ -83,34 +95,52 @@
######################################################################
# `object` interface
######################################################################
+ #Invoked during process of saving the visualization to a file
def __get_pure_state__(self):
- d = super(VTKDataSource, self).__get_pure_state__()
- data = self.data
- if data:
- w = tvtk.DataSetWriter(write_to_output_string=1)
- warn = w.global_warning_display
- w.set_input(data)
- w.global_warning_display = 0
- w.update()
- w.global_warning_display = warn
- z = gzip_string(w.output_string)
- d['data'] = z
+ d = super(CitcomSHDFFileReader, self).__get_pure_state__()
+ output = "Filename: " +self.filename+ "NX:%d NY:%d NZ:%d" %(self.nx_redu,self.ny_redu,self.nz_redu)
+ z = gzip_string(output)
+ d['data'] = z
return d
-
+
+ #When the Visualisation is opened again this Method is called
def __set_pure_state__(self, state):
z = state.data
if z:
d = gunzip_string(z)
- r = tvtk.DataSetReader(read_from_input_string=1,
- input_string=d)
- r.update()
- self.data = r.output
- # Now set the remaining state without touching the children.
- set_state(self, state, ignore=['children', 'data'])
- # Setup the children.
- handle_children_state(self.children, state.children)
- # Setup the children's state.
- set_state(self, state, first=['children'], ignore=['*'])
+ m = re.search('(?<=Filename:)\w+', header)
+ file_name = m.group(0)
+ m = re.search('(?<=NX:)\w+', d)
+ self.nx_redu = int(m.group(0))
+ m = re.search('(?<=NX:)\w+', d)
+ self.ny_redu = int(m.group(0))
+ m = re.search('(?<=NX:)\w+', d)
+ self.nz_redu = int(m.group(0))
+
+ if not isfile(file_name):
+ msg = 'Could not find file at %s\n'%file_name
+ msg += 'Please move the file there and try again.'
+ raise IOError, msg
+
+ self.filename = file_name
+
+ #self.data = self.citcomshdftougrid.initialize(self.filename,0,0,0,0,False,False)
+ f = tables.openFile(file_name,'r')
+ self.nx = int(f.root.input._v_attrs.nodex)
+ self.ny = int(f.root.input._v_attrs.nodey)
+ self.nz = int(f.root.input._v_attrs.nodez)
+
+ self.timesteps = int(f.root.time.nrows)
+ f.close()
+ self.data = self.citcomshdftougrid.initialize(self.filename,self.current_timestep,self.nx_redu,self.ny_redu,self.nz_redu,False,False)
+ self.data_changed = True
+ self._update_data()
+ # Now set the remaining state without touching the children.
+ set_state(self, state, ignore=['children', '_file_path'])
+ # Setup the children.
+ handle_children_state(self.children, state.children)
+ # Setup the children's state.
+ set_state(self, state, first=['children'], ignore=['*'])
######################################################################
# `Base` interface
@@ -136,6 +166,7 @@
def initialize(self,file_name):
+ """This Methods initializes the reader and reads the meta-information from the Hdf file """
self.filename = file_name
#self.data = self.citcomshdftougrid.initialize(self.filename,0,0,0,0,False,False)
@@ -148,9 +179,10 @@
self.ny_redu = self.ny
self.nz_redu = self.nz
- self.timesteps = int(f.root.input._v_attrs.steps)
- self.frequency = int(f.root.input._v_attrs.monitoringFrequency)
+ self.timesteps = int(f.root.time.nrows)
+ f.close()
+
######################################################################
# `TreeNodeObject` interface
######################################################################
@@ -167,6 +199,7 @@
# Non-public interface
######################################################################
def _data_changed(self, data):
+ """Invoked when the upsteam data sends an data_changed event"""
self._update_data()
self.outputs = [data]
self.data_changed = True
@@ -174,22 +207,31 @@
self.trait_property_changed('name', '',
self.tno_get_label(None))
+ ##Callbacks for our traits
def _current_timestep_changed(self,new_value):
+ """Callback for the current timestep input box"""
if new_value < 0:
current_timestep = 0
if new_value > self.timesteps:
current_timestep = 100
def _read_timestep_fired(self):
+ """Callback for the Button to read one timestep"""
self.data = self.citcomshdftougrid.initialize(self.filename,self.current_timestep,self.nx_redu,self.ny_redu,self.nz_redu,False,False)
+ self.temperature = self.citcomshdftougrid.get_vtk_temperature()
+ self.viscosity = self.citcomshdftougrid.get_vtk_viscosity()
+ self.data_changed = True
+ self._update_data()
def _nx_redu_changed(self, new_value):
+ """callback for the nx_redu input box"""
if new_value < 1:
self.nx_redu = 1
if new_value > self.nx:
self.nx_redu = self.nx
def _ny_redu_changed(self, new_value):
+ """callback for the ny_redu input box"""
if new_value < 1:
self.ny_redu = 1
if new_value > self.ny:
@@ -197,11 +239,28 @@
def _nz_redu_changed(self, new_value):
+ """callback for the nz_redu input box"""
if new_value < 1:
self.nz_redu = 1
if new_value > self.nz:
self.nz_redu = self.nz
+ def _point_scalars_name_changed(self, value):
+ if value == "Temperature":
+ self.data.point_data.scalars = self.temperature
+ if value == "Viscosity":
+ self.data.point_data.scalars = self.viscosity
+ self.data_changed = True
+ self._set_data_name('scalars', 'point', value)
+
+ def _point_vectors_name_changed(self, value):
+ self._set_data_name('vectors', 'point', value)
+
+ def _point_tensors_name_changed(self, value):
+ self._set_data_name('tensors', 'point', value)
+
+
+ ########################Non Public##############################
def _set_data_name(self, data_type, attr_type, value):
if not value:
return
@@ -217,46 +276,13 @@
# Fire an event, so the changes propagate.
self.data_changed = True
- def _point_scalars_name_changed(self, value):
- self._set_data_name('scalars', 'point', value)
-
- def _point_vectors_name_changed(self, value):
- self._set_data_name('vectors', 'point', value)
-
- def _point_tensors_name_changed(self, value):
- self._set_data_name('tensors', 'point', value)
-
- def _cell_scalars_name_changed(self, value):
- self._set_data_name('scalars', 'cell', value)
-
- def _cell_vectors_name_changed(self, value):
- self._set_data_name('vectors', 'cell', value)
-
- def _cell_tensors_name_changed(self, value):
- self._set_data_name('tensors', 'cell', value)
-
+
def _update_data(self):
if not self.data:
return
- pnt_attr, cell_attr = get_all_attributes(self.data)
-
- def _setup_data_traits(obj, attributes, d_type):
- attrs = ['scalars', 'vectors', 'tensors']
- data = getattr(obj.data, '%s_data'%d_type)
- for attr in attrs:
- values = attributes[attr]
- if values:
- default = getattr(obj, '%s_%s_name'%(d_type, attr))
- if default and default in values:
- pass
- else:
- default = values[0]
- trait = Trait(default, TraitPrefixList(values))
- getattr(data, 'set_active_%s'%attr)(default)
- else:
- trait = Trait('', TraitPrefixList(['']))
- obj.add_trait('%s_%s_name'%(d_type, attr), trait)
-
- _setup_data_traits(self, pnt_attr, 'point')
- _setup_data_traits(self, cell_attr, 'cell')
-
+ else:
+ trait = Trait('Temperature', TraitPrefixList('Temperature','Viscosity'))
+ self.add_trait('point_scalars_name', trait)
+ trait = Trait('Velocity', TraitPrefixList('Velocity'))
+ self.add_trait('point_vectors_name', trait)
+
\ No newline at end of file
Modified: mc/3D/CitcomS/trunk/visual/Mayavi2/plugins/filter/CitcomSshowCaps.py
===================================================================
--- mc/3D/CitcomS/trunk/visual/Mayavi2/plugins/filter/CitcomSshowCaps.py 2006-08-29 22:33:43 UTC (rev 4449)
+++ mc/3D/CitcomS/trunk/visual/Mayavi2/plugins/filter/CitcomSshowCaps.py 2006-08-29 23:44:53 UTC (rev 4450)
@@ -1,10 +1,21 @@
-"""A simple filter that thresholds on input data.
+"""A simple filter that thresholds CitcomS Caps from input data."""
-"""
-# Author: Prabhu Ramachandran <prabhu_r at users.sf.net>
-# Copyright (c) 2005, Enthought, Inc.
-# License: BSD Style.
+# Author: Martin Weier
+#
+#Copyright (C) 2006 California Institute of Technology
+#This program is free software; you can redistribute it and/or modify
+#it under the terms of the GNU General Public License as published by
+#the Free Software Foundation; either version 2 of the License, or
+#any later version.
+#This program is distributed in the hope that it will be useful,
+#but WITHOUT ANY WARRANTY; without even the implied warranty of
+#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#GNU General Public License for more details.
+#You should have received a copy of the GNU General Public License
+#along with this program; if not, write to the Free Software
+#Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
# Enthought library imports.
from enthought.traits import Instance, Range, Int
from enthought.traits.ui import View, Group, Item
@@ -40,16 +51,11 @@
Item(name='upper_threshold'))
)
- nx = Int()
- ny = Int()
- nz = Int()
+
n = Int()
caps = Int()
- def setvalues(self,nxin,nyin,nzin):
- self.nx = nxin
- self.ny = nyin
- self.nz = nzin
+
######################################################################
# `Filter` interface.
@@ -81,16 +87,12 @@
# input.
fil = self.ugrid_filter
fil.input = self.inputs[0].outputs[0]
+
+ #Than we calculate how many points belong to one cap
self.caps = 12
self.n = self.inputs[0].outputs[0].number_of_points/12
- #self.inputs[0].outputs[0].number_of_points/(self.nx*self.ny*self.nz)
- # We force the ranges to be reset to the limits of the data.
- # This is because if the data has changed upstream, then the
- # limits of the data must be changed.
- #self._update_ranges(reset=True)
-
- #fil.threshold_between(self.lower_threshold, self.upper_threshold)
- #fil.update()
+
+ #Than we set the output of the filter
self.outputs[0] = fil.output
self.outputs.append(self.inputs[0].outputs[0])
self.pipeline_changed = True
@@ -102,15 +104,7 @@
This method is invoked (automatically) when any of the inputs
sends a `data_changed` event.
"""
-
- #self._update_ranges(reset=True)
-
- # XXX: These are commented since updating the ranges does
- # everything for us.
-
- #fil = self.threshold_filter
- #fil.threshold_between(self.lower_threshold, self.upper_threshold)
- #fil.update()
+ fil.update()
# Propagate the data_changed event.
self.data_changed = True
@@ -118,22 +112,18 @@
# Non-public interface
######################################################################
def _lower_threshold_changed(self,old_value, new_value):
- #if new_value <> self.upper_threshold and new_value<upper_threshold:
+ """Callback interface for the lower threshold slider"""
fil = self.ugrid_filter
fil.point_minimum = (self.lower_threshold)*(self.n)
fil.update()
self.data_changed = True
- #else: #Create a single point to prevent that other filters crash because of missing input
- # self.outputs[0].points = [(100,100,100)]
def _upper_threshold_changed(self, old_value, new_value):
- #if new_value <> self.lower_threshold and new_value>lower_threshold:
+ """Callback interface for the upper threshold slider"""
fil = self.ugrid_filter
fil.point_maximum = self.upper_threshold*(self.n)
fil.update()
self.data_changed = True
- #else:
- #self.outputs[0].points = [(100,100,100)]
More information about the cig-commits
mailing list