DAMASK_EICMD/processing/post/h5_vtkRectilinearGrid.py

136 lines
4.8 KiB
Python
Raw Normal View History

#!/usr/bin/env python2.7
# -*- coding: UTF-8 no BOM -*-
# ------------------------------------------------------------------ #
# NOTE: #
# 1. It might be a good idea to separate IO and calculation. #
# 2. Some of the calculation could be useful in other situations, #
# why not build a math_util, or math_sup module that contains #
# all the useful functions. #
# ------------------------------------------------------------------ #
import os
import vtk
import numpy as np
import damask
from optparse import OptionParser
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName, damask.version])
# ----- HELPER FUNCTION ----- #
def getMeshFromXYZ(xyzArray, mode):
"""calc Vx,Vy,Vz vectors for vtk rectangular mesh"""
# NOTE:
# --> np.unique will automatically sort the list
# --> although not exactly n(1), but since mesh dimension should
# small anyway, so this is still light weight.
dim = xyzArray.shape[1] # 2D:2, 3D:3
coords = [np.unique(xyzArray[:, i]) for i in xrange(dim)]
if mode == 'cell':
# since x, y, z might now have the same number of elements,
# we have to deal with them individually
for ri in xrange(dim):
vctr_pt = coords[ri]
vctr_cell = np.empty(len(vctr_pt)+1)
# calculate first and last end point
vctr_cell[0] = vctr_pt[0] - 0.5*abs(vctr_pt[1] - vctr_pt[0])
vctr_cell[-1] = vctr_pt[-1] + 0.5*abs(vctr_pt[-2] - vctr_pt[-1])
for cj in xrange(1, len(vctr_cell)-1):
vctr_cell[cj] = 0.5*(vctr_pt[cj-1] + vctr_pt[cj])
# update the coords
coords[ri] = vctr_cell
if dim < 3:
coords.append([0]) # expand to a 3D with 0 for z
# auxiliary description
grid = np.array(map(len, coords), 'i')
N = grid.prod() if mode == 'point' else (grid-1).prod()
return coords, grid, N
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
msg = "Create regular voxel grid from points in an ASCIItable."
parser = OptionParser(option_class=damask.extendableOption,
usage='%prog options [file[s]]',
description=msg,
version=scriptID)
parser.add_option('-m',
'--mode',
dest='mode',
metavar='string',
type='choice', choices=['cell', 'point'],
help='cell-centered or point-centered coordinates')
parser.add_option('-p',
'--pos', '--position',
dest='pos',
type='string', metavar='string',
help='label of coordinates [%default]')
parser.set_defaults(mode='cell',
pos='pos')
(options, filenames) = parser.parse_args()
# ----- loop over input files ----- #
for name in filenames:
try:
h5f = damask.H5Table(name, new_file=False)
except:
continue
damask.util.report(scriptName, name)
# ----- read xyzArray from HDF5 file ----- #
xyzArray = h5f.get_data(options.pos)
# ----- figure out size and grid ----- #
coords, grid, N = getMeshFromXYZ(xyzArray, options.mode)
# ----- process data ----- #
rGrid = vtk.vtkRectilinearGrid()
# WARNING: list expansion does not work here as these are
# just pointers for a vtk instance. Simply put,
# DON't USE
# [<VTK_CONSTRUCTOR>] * <NUM_OF_ELEMENTS>
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])
# ----- output result ----- #
dirPath = os.path.split(name)[0]
if name:
writer = vtk.vtkXMLRectilinearGridWriter()
writer.SetCompressorTypeToZLib()
writer.SetDataModeToBinary()
# getting the name is a little bit tricky
vtkFileName = os.path.splitext(os.path.split(name)[1])[0]
vtkFileName += '_{}({})'.format(options.pos, options.mode)
vtkFileName += '.' + writer.GetDefaultFileExtension()
writer.SetFileName(os.path.join(dirPath, vtkFileName))
else:
writer = vtk.vtkDataSetWriter()
writer.SetHeader('# powered by '+scriptID)
writer.WriteToOutputStringOn()
if vtk.VTK_MAJOR_VERSION <= 5:
writer.SetInput(rGrid)
else:
writer.SetInputData(rGrid)
writer.Write()