for generating rectilinear grid of either cell type or point type from ASCII table.

processing/post/vtk_rectilinearGrid.py

@@ -0,0 +1,125 @@
+#!/usr/bin/env python
+# -*- coding: UTF-8 no BOM -*-
+
+import os,sys,string,re,vtk
+import numpy as np
+from optparse import OptionParser
+import damask
+
+scriptID = '$Id$'
+scriptName = os.path.splitext(scriptID.split()[1])[0]
+
+
+# --------------------------------------------------------------------
+#                                MAIN
+# --------------------------------------------------------------------
+
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+Create regular voxel grid from points in an ASCIItable.
+
+""", version = scriptID)
+
+parser.add_option('-m', '--mode',
+                  dest = 'mode',
+                  type = 'choice', choices = ['cell','point'],
+                  help = 'cell-centered or point-centered coordinates ')
+parser.add_option('-p', '--positions',
+                  dest = 'position',
+                  type = 'string', metavar = 'string',
+                  help = 'coordinate label [%default]')
+parser.set_defaults(position = 'pos',
+                   )
+
+(options, filenames) = parser.parse_args()
+
+if not options.mode: parser.error("No coordinate type specified.")
+# --- loop over input files -------------------------------------------------------------------------
+
+if filenames == []: filenames = [None]
+
+for name in filenames:
+  try:
+    table = damask.ASCIItable(name = name,
+                              buffered = False, readonly = True)
+  except: continue
+  table.croak(damask.util.emph(scriptName)+(': '+name if name else ''))
+
+# --- interpret header ----------------------------------------------------------------------------
+
+  table.head_read()
+  
+# ------------------------------------------ sanity checks ----------------------------------------
+
+  errors  = []
+  remarks = []
+  
+  if table.label_dimension(options.position) != 3:  errors.append('coordinates {} are not a vector.'.format(options.position))
+
+  if remarks != []: table.croak(remarks)
+  if errors  != []:
+    table.croak(errors)
+    table.close(dismiss = True)
+    continue
+
+# --------------- figure out size and grid ---------------------------------------------------------
+
+  table.data_readArray(options.position)
+
+  coords = [np.unique(table.data[:,i]) for i in xrange(3)]
+  if options.mode == 'cell':
+    coords = [0.5 * np.array([3.0 * coords[i][0] - coords[i][0 + len(coords[i]) > 1]] + \
+                             [coords[i][j-1] + coords[i][j] for j in xrange(1,len(coords[i]))] + \
+                             [3.0 * coords[i][-1] - coords[i][-1 - (len(coords[i]) > 1)]]) for i in xrange(3)]
+  grid   = np.array(map(len,coords),'i')
+  N = grid.prod() if options.mode == 'point' else (grid-1).prod()
+  
+  if N != len(table.data): errors.append('data count {} does not match grid {}x{}x{}.'.format(N,*(grid - options.mode == 'cell') ))
+  if errors  != []:
+    table.croak(errors)
+    table.close(dismiss = True)
+    continue
+
+# ------------------------------------------ process data ---------------------------------------  
+
+  rGrid = vtk.vtkRectilinearGrid()
+  coordArray = [vtk.vtkDoubleArray(),
+                vtk.vtkDoubleArray(),
+                vtk.vtkDoubleArray(),
+               ]
+  
+  rGrid.SetDimensions(*grid)
+  for i,points in enumerate(coords):
+    for point in points:
+      coordArray[i].InsertNextValue(point)
+ 
+  rGrid.SetXCoordinates(coordArray[0])
+  rGrid.SetYCoordinates(coordArray[1])
+  rGrid.SetZCoordinates(coordArray[2])
+
+  table.croak('{} points and {} cells...'.format(rGrid.GetNumberOfPoints(),rGrid.GetNumberOfCells(),))
+ 
+# ------------------------------------------ output result ---------------------------------------  
+
+  if name:
+    (dir,filename) = os.path.split(name)
+    writer = vtk.vtkXMLRectilinearGridWriter()
+    writer.SetDataModeToBinary()
+    writer.SetCompressorTypeToZLib()
+    writer.SetFileName(os.path.join(dir,os.path.splitext(filename)[0]+'_{}'.format(options.position) \
+                                        +'.'+writer.GetDefaultFileExtension()))
+    if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(rGrid)
+    else:                          writer.SetInputData(rGrid)
+    writer.Write()
+
+  else:
+    writer = vtk.vtkRectilinearGridWriter()
+    writer.WriteToOutputStringOn()
+    writer.SetFileTypeToASCII()
+    writer.SetHeader('# powered by '+scriptID)
+    if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(rGrid)
+    else:                          writer.SetInputData(rGrid)
+    writer.Write()
+    sys.stdout.write(writer.GetOutputString()[0:writer.GetOutputStringLength()])
+
+  table.close()                                                                                     # close input ASCII table
+